#include "rundeck_opts.h" C23456789012345678901234567890123456789012345678901234567890123456789012 MODULE LAKES 8,5 !@sum LAKES subroutines for Lakes and Rivers !@auth Gavin Schmidt/Gary Russell !@ver 1.0 (based on LB265) USE CONSTANT, only : grav,bygrav,shw,rhow,lhm,shi,teeny USE MODEL_COM, only : im,jm USE DOMAIN_DECOMP, only : HALO_UPDATE, GRID,NORTH,SOUTH USE DOMAIN_DECOMP, only : WRITE_PARALLEL #ifdef TRACERS_WATER USE TRACER_COM, only : trname,ntm #endif IMPLICIT NONE SAVE C**** C**** Changes from Model III: MO -> MWL (kg), G0M -> GML (J), C**** GZM -> TLAKE (deg C) !@var KDIREC directions for river flow C**** (0 no flow, 1-8 anti-clockwise from top RH corner INTEGER, ALLOCATABLE, DIMENSION(:,:) :: KDIREC !@var RATE rate of river flow downslope (fraction) !@var DHORZ horizontal distance to downstream box (m) REAL*8, ALLOCATABLE, DIMENSION(:,:) :: RATE,DHORZ !@var IFLOW,JFLOW grid box indexes for downstream direction INTEGER, ALLOCATABLE, DIMENSION (:,:) :: IFLOW,JFLOW !@param NRVRMX Max No. of specified rivers INTEGER, PARAMETER :: NRVRMX = 42 !@var NRVR actual No. of specified rivers INTEGER :: NRVR !@var IRVRMTH,JRVRMTH indexes for specified river mouths INTEGER, DIMENSION(NRVRMX) :: IRVRMTH,JRVRMTH !@var NAMERVR Names of specified rivers CHARACTER*8, DIMENSION(NRVRMX) :: NAMERVR !@param MINMLD minimum mixed layer depth in lake (m) REAL*8, PARAMETER :: MINMLD = 1. !@param TMAXRHO temperature of maximum density (pure water) (C) REAL*8, PARAMETER :: TMAXRHO = 4. !@param KVLAKE lake diffusion constant at mixed layer depth (m^2/s) REAL*8, PARAMETER :: KVLAKE = 1d-5 !@param TFL freezing temperature for lakes (=0 C) REAL*8, PARAMETER :: TFL = 0. !@param AC1LMIN, AC2LMIN minimum ice thickness for lake ice (kg/m^2) REAL*8, PARAMETER :: AC1LMIN = 0.1, AC2LMIN=0.1 ! (not used) !@param FLEADLK lead fraction for lakes REAL*8, PARAMETER :: FLEADLK = 0. !@param BYZETA reciprocal of solar rad. extinction depth for lake (1/m) REAL*8, PARAMETER :: BYZETA = 1./0.35d0 !@dbparam river_fac Factor to multiply runoff by to balance sea level REAL*8 :: river_fac=1. ! default = 1 !@dbparam init_flake used to make sure FLAKE is properly initialised !@+ when using older restart files INTEGER :: init_flake=1 ! default = 1 !@dbparam variable_lk 1 if lakes are to be variable !@+ (temporary variable for development purposes) INTEGER :: variable_lk=0 ! default = 0 CONTAINS SUBROUTINE LKSOURC (I0,J0,ROICE,MLAKE,ELAKE,RUN0,FODT,FIDT,SROX 1,3 * ,FSR2, #ifdef TRACERS_WATER * TRLAKEL,TRUN0,TREVAP,TRO,TRI, #endif * EVAPO,ENRGFO,ACEFO,ACEFI,ENRGFI) !@sum LKSOURC applies fluxes to lake in ice-covered and ice-free areas !@auth Gary Russell/Gavin Schmidt !@ver 1.0 USE MODEL_COM, only : qcheck IMPLICIT NONE !@var MLAKE,ELAKE mass and energy in lake layers (kg,J /m^2) REAL*8, INTENT(INOUT), DIMENSION(2) :: MLAKE,ELAKE INTEGER, INTENT(IN) :: I0,J0 REAL*8, INTENT(IN) :: ROICE, EVAPO, RUN0 REAL*8, INTENT(IN) :: FODT, FIDT, SROX(2) REAL*8, INTENT(OUT) :: ENRGFO, ACEFO, ENRGFI, ACEFI #ifdef TRACERS_WATER REAL*8, INTENT(INOUT), DIMENSION(NTM,2) :: TRLAKEL REAL*8, INTENT(IN), DIMENSION(NTM) :: TRUN0,TREVAP REAL*8, INTENT(OUT), DIMENSION(NTM) :: TRO,TRI REAL*8, DIMENSION(NTM) :: DTR2,TRUNO,TRUNI,TRF1,TRF2,FRAC #ifdef TRACERS_SPECIAL_O18 REAL*8 fracls INTEGER N #endif #endif !@var emin min energy deficit required before ice forms (J/m^2) REAL*8, PARAMETER :: emin=-1d-10 REAL*8 ENRGF1, ACEF1, ENRGF2, ACEF2, FHO, FHI, FH0, FH1, FH2, FSR2 REAL*8 ENRGI, ENRGI2, ENRGO, ENRGO2, RUNO, RUNI, TLK2, DM2, DH2 REAL*8 FRATO,FRATI,E2O,E2I !@var out_line local variable to hold mixed-type output for parallel I/O character(len=300) :: out_line C**** initiallize output ENRGFO=0. ; ACEFO=0. ; ACEFI=0. ; ENRGFI=0. C**** Calculate heat and mass fluxes to lake ENRGO = FODT-SROX(1)*FSR2 ! in open water ENRGO2= +SROX(1)*FSR2 ! in open water, second layer ENRGI = FIDT-SROX(2)*FSR2 ! under ice ENRGI2= +SROX(2)*FSR2 ! under ice, second layer RUNO =-EVAPO RUNI = RUN0 #ifdef TRACERS_WATER TRUNO(:)=-TREVAP(:) TRUNI(:)= TRUN0(:) FRAC(:)=1. #ifdef TRACERS_SPECIAL_O18 do n=1,ntm FRAC(n)=fracls(n) ! fractionation when freezing end do #endif #endif C**** Bring up mass from second layer if required/allowed IF (MLAKE(1)+RUNO.lt.MINMLD*RHOW.and.MLAKE(2).gt.0) THEN DM2 = MIN(MLAKE(2),MINMLD*RHOW-(MLAKE(1)+RUNO)) DH2 = DM2*(ELAKE(2)+(1.-ROICE)*ENRGO2+ROICE*ENRGI2)/MLAKE(2) #ifdef TRACERS_WATER DTR2(:) = DM2*TRLAKEL(:,2)/MLAKE(2) #endif ELSE DM2 = 0. DH2 = 0. #ifdef TRACERS_WATER DTR2(:) = 0. #endif END IF C**** Apply fluxes to 2nd layer IF (DM2.lt.MLAKE(2)) THEN MLAKE(2)=MLAKE(2) - DM2 ELAKE(2)=ELAKE(2) - DH2 + (1.-ROICE)*ENRGO2 + ROICE*ENRGI2 #ifdef TRACERS_WATER TRLAKEL(:,2)=TRLAKEL(:,2) - DTR2(:) #endif ELSE MLAKE(2)=0. ELAKE(2)=0. #ifdef TRACERS_WATER TRLAKEL(:,2)=0. #endif END IF E2O = 0. ; E2I = 0. C**** Calculate energy in mixed layer (open ocean) IF (ROICE.LT.1d0) THEN FHO=ELAKE(1)+ENRGO+DH2-(MLAKE(1)+DM2+RUNO)*TFL*SHW IF (FHO.LT.emin) THEN ! FLUXES COOL WATER TO FREEZING, FORM ICE ACEFO =FHO/(TFL*(SHI-SHW)-LHM) ACEFO =MIN(ACEFO,MAX(MLAKE(1)+DM2+RUNO-MINMLD*RHOW,0d0)) ENRGFO=ACEFO*(TFL*SHI-LHM) E2O=FHO-ENRGFO END IF END IF IF (ROICE.GT.0) THEN C**** Calculate energy in mixed layer (under ice) FHI=ELAKE(1)+DH2+ENRGI-(MLAKE(1)+DM2+RUNI)*TFL*SHW IF (FHI.LT.emin) THEN ! FLUXES COOL WATER TO FREEZING, FORM ICE ACEFI =FHI/(TFL*(SHI-SHW)-LHM) ACEFI =MIN(ACEFI,MAX(MLAKE(1)+DM2+RUNI-MINMLD*RHOW,0d0)) ENRGFI=ACEFI*(TFL*SHI-LHM) E2I=FHI-ENRGFI END IF END IF #ifdef TRACERS_WATER TRO(:)=ACEFO*FRAC(:)*TRLAKEL(:,1)/MLAKE(1) TRI(:)=ACEFI*FRAC(:)*TRLAKEL(:,1)/MLAKE(1) #endif C**** Update first layer variables MLAKE(1)=MLAKE(1)+DM2+(1.-ROICE)*(RUNO -ACEFO)+ROICE*(RUNI-ACEFI) ELAKE(1)=ELAKE(1)+DH2+(1.-ROICE)*(ENRGO-ENRGFO)+ * ROICE*(ENRGI-ENRGFI) #ifdef TRACERS_WATER TRLAKEL(:,1)=TRLAKEL(:,1)+DTR2(:)+(1.-ROICE)*(TRUNO(:)-TRO(:))+ * ROICE*(TRUNI(:)-TRI(:)) #endif ACEF1=0. ; ACEF2=0. ; ENRGF1=0. ; ENRGF2=0. C**** Take remaining energy and start to freeze second layer FH2= ELAKE(1)-MLAKE(1)*TFL*SHW IF (FH2.LT.emin) THEN IF (MLAKE(2).gt.0) THEN C**** FH2=-ACEF2*(TLK2-TFL)*SHW+ACEF2*LHM TLK2 =ELAKE(2)/(MLAKE(2)*SHW) ACEF2 =-FH2/(TLK2*SHW-TFL*SHI+LHM) ACEF2 =MIN(ACEF2,MLAKE(2)) ENRGF2 =ACEF2*(TFL*SHI-LHM) ELAKE(1)=ELAKE(1)+ACEF2*TLK2*SHW-ENRGF2 ELAKE(2)=ELAKE(2)-ACEF2*TLK2*SHW MLAKE(2)=MLAKE(2)-ACEF2 END IF FH1= ELAKE(1)-MLAKE(1)*TFL*SHW IF (FH1.lt.emin) THEN ! all layer 2 froze, freeze layer 1 ACEF1 =FH1/(TFL*(SHI-SHW)-LHM) C**** limit freezing if lake is between 50 and 20cm depth IF (MLAKE(1).lt.0.5d0*RHOW) THEN ACEF1=MIN(ACEF1,MAX(0.5*(MLAKE(1)-0.2d0*RHOW),0d0)) if (qcheck) print*,"Lake freezing limited",ACEF1/RHOW * ,MLAKE(1)/RHOW END IF ENRGF1 =ACEF1*(TFL*SHI-LHM) ELAKE(1)=ELAKE(1)-ENRGF1 MLAKE(1)=MLAKE(1)-ACEF1 FH0 =ELAKE(1)-MLAKE(1)*TFL*SHW IF (FH0.lt.-1d-8) THEN ! max. amount of lake frozen, cool ice if (qcheck) then WRITE(out_line,*) * "Minimum lake level reached: rsi,mlake,elake",i0,j0 * ,roice,mlake(1)/rhow,elake(1) CALL WRITE_PARALLEL(trim(out_line), UNIT=6) endif ENRGF1 =ENRGF1+FH0 ELAKE(1)=MLAKE(1)*TFL*SHW END IF END IF END IF #ifdef TRACERS_WATER TRF1(:) = ACEF1*FRAC(:)*TRLAKEL(:,1)/(MLAKE(1)+ACEF1) TRLAKEL(:,1)=TRLAKEL(:,1)-TRF1(:) IF (MLAKE(2).gt.0) THEN TRF2(:) = MIN(ACEF2*FRAC(:)/(MLAKE(2)+ACEF2),1d0)*TRLAKEL(:,2) TRLAKEL(:,2)=TRLAKEL(:,2)-TRF2(:) ELSE ! possibility of complete freezing (and so no frac) TRF2(:) = TRLAKEL(:,2) TRLAKEL(:,2) = 0. END IF #endif C**** combine mass and energy fluxes for output C**** Note that output fluxes are over open water/ice covered fractions C**** distribute ice fluxes according to flux amounts FRATO = 1d0 FRATI = 1d0 IF (E2I+E2O.lt.0) THEN FRATO = E2O/(E2I*ROICE+E2O*(1.-ROICE)) FRATI = E2I/(E2I*ROICE+E2O*(1.-ROICE)) END IF ACEFO =ACEFO + (ACEF1 +ACEF2 )*FRATO ACEFI =ACEFI + (ACEF1 +ACEF2 )*FRATI ENRGFO=ENRGFO+ (ENRGF1+ENRGF2)*FRATO ENRGFI=ENRGFI+ (ENRGF1+ENRGF2)*FRATI #ifdef TRACERS_WATER TRO(:)=TRO(:) + (TRF1(:) + TRF2(:))* FRATO TRI(:)=TRI(:) + (TRF1(:) + TRF2(:))* FRATI #endif RETURN END SUBROUTINE LKSOURC SUBROUTINE LKMIX(MLAKE,ELAKE, 1 #ifdef TRACERS_WATER * TRLAKEL, #endif * HLAKE,TKE,ROICE,DTSRC) !@sum LKMIX calculates mixing and entrainment in lakes !@auth Gavin Schmidt !@ver 1.0 IMPLICIT NONE !@var MLAKE,ELAKE mass and energy in lake layers (kg,J /m^2) REAL*8, INTENT(INOUT), DIMENSION(2) :: MLAKE,ELAKE !@var TKE turbulent kinetic energy input at surface of lake (J/m^2) !@var ROICE ice fraction REAL*8, INTENT(IN) :: TKE,ROICE !@var HLAKE sill depth for lake (m) REAL*8, INTENT(IN) :: HLAKE !@var DTSRC source time step (s) REAL*8, INTENT(IN) :: DTSRC #ifdef TRACERS_WATER !@var TRLAKEL tracer mass in lake layers (kg/m^2) REAL*8, INTENT(INOUT), DIMENSION(NTM,2) :: TRLAKEL REAL*8, DIMENSION(NTM) :: DTML,TR1N,TR2N,TRLT #endif !@param MAXRHO,RHO0,BFAC freshwater density function approximation REAL*8, PARAMETER :: MAXRHO=1d3, RHO0=999.842594d0, * BFAC=(MAXRHO-RHO0)/16d0 REAL*8 TLK1, TLK2, HLT, MLT, DTK, E1N, E2N, ATKE, H1, H2, * DRHO, DML, DHML C**** Only mix if there is a second layer! IF (MLAKE(2).gt.0) THEN TLK1=ELAKE(1)/(MLAKE(1)*SHW) TLK2=ELAKE(2)/(MLAKE(2)*SHW) HLT=ELAKE(1)+ELAKE(2) MLT=MLAKE(1)+MLAKE(2) #ifdef TRACERS_WATER TRLT(:)=TRLAKEL(:,1)+TRLAKEL(:,2) #endif C**** Test for static stability C**** DRHO=RHO(TLK2)-RHO(TLK1)~=(TLK2-TLK1)*dRHOdT((TLK1+TLK2)/2) C**** Assumes a parabolic density function going through MAXRHO at C**** TMAXRHO, and RHO0 at T=0. (reasonable up to about 12 C) IF ((TMAXRHO-0.5*(TLK1+TLK2))*(TLK2-TLK1).lt.0) THEN C**** mix uniformly and set MLD to minimum MLAKE(1)=MIN(MLT,MAX(MINMLD*RHOW,MLT-HLAKE*RHOW)) MLAKE(2)=MLT-MLAKE(1) ELAKE(1)=HLT*MLAKE(1)/MLT ELAKE(2)=HLT*MLAKE(2)/MLT #ifdef TRACERS_WATER TRLAKEL(:,1)=TRLT(:)*MLAKE(1)/MLT TRLAKEL(:,2)=TRLT(:)*MLAKE(2)/MLT #endif ELSE ! not unstable, implicitly diffuse heat + entrain C**** reduce mixing if there is ice cover DTK=2.*KVLAKE*(1.-ROICE)*DTSRC*RHOW**2 E1N=(ELAKE(1)+DTK*HLT/(MLT*MLAKE(2)))/ * (1.+DTK/(MLAKE(1)*MLAKE(2))) E2N=(ELAKE(2)+DTK*HLT/(MLT*MLAKE(1)))/ * (1.+DTK/(MLAKE(1)*MLAKE(2))) ELAKE(1)=E1N ELAKE(2)=E2N #ifdef TRACERS_WATER C**** diffuse tracers using same KV as for heat? TR1N(:)=(TRLAKEL(:,1)+DTK*TRLT(:)/(MLT*MLAKE(2)))/ * (1.+DTK/(MLAKE(1)*MLAKE(2))) TR2N(:)=(TRLAKEL(:,2)+DTK*TRLT(:)/(MLT*MLAKE(1)))/ * (1.+DTK/(MLAKE(1)*MLAKE(2))) TRLAKEL(:,1)=TR1N(:) TRLAKEL(:,2)=TR2N(:) #endif C**** entrain deep water if there is available TKE C**** take a factor of TKE and calculate change in PE IF (TKE.gt.0) THEN ATKE=0.2d0*TKE ! 20% of TKE is available for mixing H1=MLAKE(1)/RHOW H2=MLAKE(2)/RHOW DRHO=(TLK2-TLK1)*2d0*BFAC*(TMAXRHO-0.5*(TLK1+TLK2)) DML=ATKE*BYGRAV/(DRHO*0.5*H1) IF (DML*RHOW.lt.MLAKE(2)) THEN DHML=DML*ELAKE(2)/H2 ELAKE(1)=ELAKE(1)+DHML ELAKE(2)=ELAKE(2)-DHML MLAKE(1)=MLAKE(1)+DML*RHOW MLAKE(2)=MLAKE(2)-DML*RHOW #ifdef TRACERS_WATER DTML(:)=DML*TRLAKEL(:,2)/H2 TRLAKEL(:,1)=TRLAKEL(:,1)+DTML(:) TRLAKEL(:,2)=TRLAKEL(:,2)-DTML(:) #endif ELSE ! entire second layer is entrained MLAKE(1)=MLT MLAKE(2)=0 ELAKE(1)=HLT ELAKE(2)=0 #ifdef TRACERS_WATER TRLAKEL(:,1)=TRLT(:) TRLAKEL(:,2)=0. #endif END IF END IF END IF END IF RETURN C**** END SUBROUTINE LKMIX END MODULE LAKES SUBROUTINE ALLOC_LAKES (GRID) 1,3 C23456789012345678901234567890123456789012345678901234567890123456789012 !@SUM To alllocate arrays whose sizes now need to be determined !@+ at run-time !@auth Raul Garza-Robles !@ver 1.0 USE DOMAIN_DECOMP, only: DIST_GRID, GET USE MODEL_COM, only : IM, JM USE LAKES, ONLY: RATE, DHORZ,KDIREC,IFLOW,JFLOW IMPLICIT NONE TYPE (DIST_GRID), INTENT(IN) :: grid ALLOCATE ( KDIREC (IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO), * IFLOW (IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO), * JFLOW (IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO), * RATE (IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO), * DHORZ (IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO) * ) RETURN END SUBROUTINE ALLOC_LAKES SUBROUTINE init_LAKES(inilake,istart) 1,30 !@sum init_LAKES initiallises lake variables !@auth Gary Russell/Gavin Schmidt !@ver 1.0 USE FILEMANAGER USE CONSTANT, only : rhow,shw,tf,pi,grav USE MODEL_COM, only : im,jm,flake0,zatmo,dtsrc,flice,hlake * ,focean,jday,fearth0 USE DOMAIN_DECOMP, only : GRID,WRITE_PARALLEL USE DOMAIN_DECOMP, only : GET,NORTH,SOUTH,HALO_UPDATE c*** USE ESMF_MOD, Only : ESMF_HaloDirection USE GEOM, only : dxyp,dxv,dyv,dxp,dyp,imaxj #ifdef TRACERS_WATER USE TRACER_COM, only : trw0 USE FLUXES, only : gtracer #endif USE FLUXES, only : gtemp,mlhc,gtempr USE SEAICE_COM, only : rsi USE PBLCOM, only : tsavg USE LAKES USE LAKES_COM !USE GHY_COM, only : fearth USE DIAG_COM, only : npts,icon_LKM,icon_LKE,title_con,conpt0 USE PARAM IMPLICIT NONE INTEGER :: FROM,J_0,J_1,J_0H,J_1H,J_0S,J_1S LOGICAL :: HAVE_NORTH_POLE, HAVE_SOUTH_POLE c*** Type (ESMF_HaloDirection) :: direction Integer :: direction ! ESMF_HaloDirection not yet implemented LOGICAL, INTENT(IN) :: inilake INTEGER, INTENT(IN) :: ISTART !@var I,J,I72,IU,JU,ID,JD loop variables INTEGER I,J,I72,IU,JU,ID,JD,INM,KD INTEGER iu_RVR !@var iu_RVR unit number for river direction file CHARACTER TITLEI*80, CDIREC(IM,JM)*1, CONPT(NPTS)*10 REAL*8 SPMIN,SPMAX,SPEED0,SPEED,DZDH,DZDH1,MLK1 LOGICAL :: QCON(NPTS), T=.TRUE. , F=.FALSE. !@var out_line local variable to hold mixed-type output for parallel I/O character(len=300) :: out_line CALL GET(GRID, J_STRT = J_0, J_STOP = J_1, & J_STRT_SKP = J_0S, J_STOP_SKP = J_1S, & J_STRT_HALO= J_0H, J_STOP_HALO= J_1H, & HAVE_SOUTH_POLE = HAVE_SOUTH_POLE, & HAVE_NORTH_POLE = HAVE_NORTH_POLE) C**** C**** LAKECB MWL Mass of water in lake (kg) C**** GML Liquid lake enthalpy (J) C**** TLAKE Temperature of lake surface (C) C**** HLAKE Lake sill depth (m) C**** ALAKE Lake surface area (m^2) C**** TANLK Lake slope (tan(alpha)) (1) C**** C**** FIXDCB FLAKE0 Original Lake fraction (1) C**** C**** Get parameters from rundeck call sync_param("river_fac",river_fac) call sync_param("init_flake",init_flake) call sync_param("variable_lk",variable_lk) C**** initialise FLAKE if requested (i.e. from older restart files) if ((init_flake.eq.1.and.istart.lt.9) .or. INILAKE) THEN print*,"Initialising FLAKE from TOPO file..." FLAKE = FLAKE0 end if C**** Ensure that HLAKE is a minimum of 1m for FLAKE>0 DO J=J_0, J_1 DO I=1,IM IF (FLAKE(I,J).gt.0 .and. HLAKE(I,J).lt.1.) THEN print*,"Warning: Fixing HLAKE",i,j,FLAKE(I,J),FLAKE0(I,J) * ,HLAKE(I,J),"--> 1m" HLAKE(I,J)=1. END IF END DO END DO IF (INILAKE) THEN C**** Set lake variables from surface temperature C**** This is just an estimate for the initiallisation DO J=J_0, J_1 DO I=1,IM IF (FLAKE(I,J).gt.0) THEN TLAKE(I,J) = MAX(0d0,TSAVG(I,J)-TF) MWL(I,J) = RHOW*HLAKE(I,J)*FLAKE(I,J)*DXYP(J) MLK1 = MINMLD*RHOW*FLAKE(I,J)*DXYP(J) GML(I,J) = SHW*(MLK1*TLAKE(I,J) * +(MWL(I,J)-MLK1)*MAX(TLAKE(I,J),4d0)) MLDLK(I,J) = MINMLD #ifdef TRACERS_WATER TRLAKE(:,1,I,J)=MLK1*TRW0(:) TRLAKE(:,2,I,J)=(MWL(I,J)-MLK1)*TRW0(:) #endif ELSE TLAKE(I,J) = 0. MWL(I,J) = 0. GML(I,J) = 0. MLDLK(I,J) = MINMLD #ifdef TRACERS_WATER TRLAKE(:,:,I,J)=0. #endif END IF END DO END DO END IF C**** Set fixed geometric variables C**** TANLK=TAN(ALPHA) = R/H for a conical lake of equivalent volume DO J=J_0, J_1 DO I=1,IM IF (FLAKE0(I,J).gt.0) THEN TANLK(I,J) = SQRT(FLAKE0(I,J)*DXYP(J)/PI)/(3d0*HLAKE(I,J)) ELSE TANLK(I,J) = 2d3 ! reasonable average value END IF END DO END DO CALL PRINTLK("IN") C**** Set GTEMP arrays for lakes IF (ISTART.gt.0) THEN DO J=J_0, J_1 DO I=1,IM IF (FLAKE(I,J).gt.0) THEN GTEMP(1,1,I,J)=TLAKE(I,J) GTEMPR(1,I,J) =TLAKE(I,J)+TF IF (MWL(I,J).gt.(1d-10+MLDLK(I,J))*RHOW*FLAKE(I,J)*DXYP(J)) * THEN GTEMP(2,1,I,J)=(GML(I,J)-TLAKE(I,J)*SHW*MLDLK(I,J)*RHOW * *FLAKE(I,J)*DXYP(J))/(SHW*(MWL(I,J)-MLDLK(I,J) * *RHOW*FLAKE(I,J)*DXYP(J))) C**** If starting from a possibly corrupted rsf file, check Tlk2 IF(GTEMP(2,1,I,J)>TLAKE(I,J)+1.and.GTEMP(2,1,I,J)>10) THEN WRITE(6,*) "Warning: Unphysical Tlk2 fixed",I,J,GTEMP(: * ,1,I,J) GTEMP(2,1,I,J)=GTEMP(1,1,I,J) ! set to Tlk1 GML(I,J)=TLAKE(I,J)*SHW*MWL(I,J) END IF ELSE GTEMP(2,1,I,J)=TLAKE(I,J) END IF #ifdef TRACERS_WATER GTRACER(:,1,I,J)=TRLAKE(:,1,I,J)/(MLDLK(I,J)*RHOW*FLAKE(I,J) * *DXYP(J)) #endif MLHC(I,J)= SHW*MLDLK(I,J)*RHOW END IF END DO END DO END IF C**** C**** Always initiallise River direction and Rate C**** Read in CDIREC: Number = octant direction, Letter = river mouth call openunit("RVR",iu_RVR,.false.,.true.) READ (iu_RVR,910) TITLEI WRITE (out_line,*) 'River Direction file read: ',TITLEI CALL WRITE_PARALLEL(trim(out_line), UNIT=6) READ (iu_RVR,910) DO I72=1,1+(IM-1)/72 DO J=JM, 1, -1 READ (iu_RVR,911) (CDIREC(I,J),I=72*(I72-1)+1,MIN(IM,I72*72)) END DO END DO C**** read in named rivers (if any) READ (iu_RVR,*,END=10) READ (iu_RVR,'(A80)',END=10) TITLEI READ (iu_RVR,*,END=10) IF (TITLEI.eq."Named River Mouths:") THEN DO I=1,NRVRMX,5 READ(iu_RVR,'(5(A8,1X))') NAMERVR(I:MIN(NRVRMX,I+4)) END DO END IF 10 call closeunit (iu_RVR) C**** Create integral direction array KDIREC from CDIREC CALL HALO_UPDATE(GRID, FEARTH0, FROM=NORTH+SOUTH) CALL HALO_UPDATE(GRID, FLICE, FROM=NORTH+SOUTH) CALL HALO_UPDATE(GRID, FLAKE0, FROM=NORTH+SOUTH) ! fixed CALL HALO_UPDATE(GRID, FOCEAN, FROM=NORTH+SOUTH) ! fixed ! Use unusual loop bounds to fill KDIREC in halo DO J=MAX(1,J_0-1),MIN(JM,J_1+1) DO I=1,IM C**** KD: -16 = blank, 0-8 directions >8 named rivers KD= ICHAR(CDIREC(I,J)) - 48 C**** If land but no ocean, and no direction, print warning IF ((FEARTH0(I,J)+FLICE(I,J)+FLAKE0(I,J).gt.0) .and. * FOCEAN(I,J).le.0 .and. (KD.gt.8 .or. KD.lt.0)) THEN WRITE(6,*) "Land box has no river direction I,J: ",I,J * ,FOCEAN(I,J),FLICE(I,J),FLAKE0(I,J),FEARTH0(I,J) END IF C**** Default direction is down (if ocean box), or no outlet (if not) C**** Also ensure that all ocean boxes are done properly IF ((KD.lt.0 .or. KD.gt.8) .or. FOCEAN(I,J).eq.1.) THEN KDIREC(I,J)=0. ELSE KDIREC(I,J) = KD END IF C**** Check for specified river mouths IF (KD.GE.17 .AND. KD.LE.42) THEN IF (FOCEAN(I,J).le.0) THEN WRITE(6,*) * "Warning: Named river outlet must be in ocean",i * ,j,FOCEAN(I,J),FLICE(I,J),FLAKE0(I,J) * ,FEARTH0(I,J) END IF END IF END DO END DO INM=0 DO J=1,JM DO I=1,IM C**** KD: -16 = blank, 0-8 directions >8 named rivers KD= ICHAR(CDIREC(I,J)) - 48 C**** Check for specified river mouths IF (KD.GE.17 .AND. KD.LE.42) THEN INM=INM+1 IRVRMTH(INM)=I JRVRMTH(INM)=J IF (CDIREC(I,J).ne.NAMERVR(INM)(1:1)) THEN WRITE(6,*) * "Warning: Named river in RVR does not correspond" * //" with letter in direction file. Please check" WRITE(6,*) "INM, CDIREC, NAMERVR = ",INM,CDIREC(I,J) * ," ",NAMERVR(INM) NAMERVR(INM)=CDIREC(I,J) ! set default END IF END IF END DO END DO NRVR=INM C**** C**** From each box calculate the downstream river box C**** ! odd bounds to fill IFLOW and JFLOW in halo DO J=MAX(2,J_0H), MIN(JM-1,J_1H) DO I=1,IM SELECT CASE (KDIREC(I,J)) CASE (0) IFLOW(I,J) = I JFLOW(I,J) = J DHORZ(I,J) = 0.5*SQRT(DXP(J)*DXP(J)+DYP(J)*DYP(J)) CASE (1) IFLOW(I,J) = I+1 JFLOW(I,J) = J+1 DHORZ(I,J) = SQRT(DXV(J+1)*DXV(J+1)+DYV(J+1)*DYV(J+1)) ! SQRT(DXV(J)*DXV(J)+DYV(J)*DYV(J)) IF(I.eq.IM) IFLOW(I,J) = 1 CASE (2) IFLOW(I,J) = I JFLOW(I,J) = J+1 DHORZ(I,J) = DYV(J+1) ! DYV(J) CASE (3) IFLOW(I,J) = I-1 JFLOW(I,J) = J+1 DHORZ(I,J) = SQRT(DXV(J+1)*DXV(J+1)+DYV(J+1)*DYV(J+1)) ! SQRT(DXV(J)*DXV(J)+DYV(J)*DYV(J)) IF(I.eq.1) IFLOW(I,J) = IM CASE (4) IFLOW(I,J) = I-1 JFLOW(I,J) = J DHORZ(I,J) = DXP(J) IF(I.eq.1) IFLOW(I,J) = IM CASE (5) IFLOW(I,J) = I-1 JFLOW(I,J) = J-1 DHORZ(I,J) = SQRT(DXV(J)*DXV(J)+DYV(J)*DYV(J)) ! SQRT(DXV(J-1)*DXV(J-1)+DYV(J-1)*DYV(J-1)) IF(I.eq.1) IFLOW(I,J) = IM CASE (6) IFLOW(I,J) = I JFLOW(I,J) = J-1 DHORZ(I,J) = DYV(J) ! DYV(J-1) CASE (7) IFLOW(I,J) = I+1 JFLOW(I,J) = J-1 DHORZ(I,J) = SQRT(DXV(J)*DXV(J)+DYV(J)*DYV(J)) ! SQRT(DXV(J-1)*DXV(J-1)+DYV(J-1)*DYV(J-1)) IF(I.eq.IM) IFLOW(I,J) = 1 CASE (8) IFLOW(I,J) = I+1 JFLOW(I,J) = J DHORZ(I,J) = DXP(J) IF(I.eq.IM) IFLOW(I,J) = 1 END SELECT END DO END DO C**** South Pole is a special case IF (HAVE_SOUTH_POLE) Then DO I=1,IM IF(KDIREC(I,1).eq.2) THEN IFLOW(1,1) = I JFLOW(1,1) = 2 DHORZ(1,1) = DYV(2) ! DYV(1) END IF IF(KDIREC(I,2).eq.6) THEN IFLOW(I,2) = 1 JFLOW(I,2) = 1 END IF END DO END IF C**** C**** Calculate river flow RATE (per source time step) C**** CALL HALO_UPDATE(GRID, zatmo, FROM=NORTH+SOUTH) SPEED0= .35d0 SPMIN = .15d0 SPMAX = 5. DZDH1 = .00005 DO JU = J_0, J_1S DO IU=1,IMAXJ(JU) IF(KDIREC(IU,JU).gt.0) THEN JD=JFLOW(IU,JU) ID=IFLOW(IU,JU) DZDH = (ZATMO(IU,JU)-ZATMO(ID,JD)) / (GRAV*DHORZ(IU,JU)) ELSE DZDH = ZATMO(IU,JU) / (GRAV*DHORZ(IU,JU)) END IF SPEED = SPEED0*DZDH/DZDH1 IF(SPEED.lt.SPMIN) SPEED = SPMIN IF(SPEED.gt.SPMAX) SPEED = SPMAX RATE(IU,JU) = DTsrc*SPEED/DHORZ(IU,JU) END DO END DO C**** Set conservation diagnostics for Lake mass and energy CONPT=CONPT0 CONPT(4)="PREC+LAT M" CONPT(5)="SURFACE" ; CONPT(8)="RIVERS" QCON=(/ F, F, F, T, T, F, F, T, T, F, F/) CALL SET_CON(QCON,CONPT,"LAK MASS","(KG/M^2) ", * "(10**-9 KG/SM^2)",1d0,1d9,icon_LKM) QCON=(/ F, F, F, T, T, F, F, T, T, F, F/) CALL SET_CON(QCON,CONPT,"LAK ENRG","(10**3 J/M^2) ", * "(10**-3 W/M^2) ",1d-3,1d3,icon_LKE) C**** assume that at the start GHY is in balance with LAKES SVFLAKE = FLAKE RETURN C**** 910 FORMAT (A72) 911 FORMAT (72A1) END SUBROUTINE init_LAKES SUBROUTINE RIVERF 1,26 !@sum RIVERF transports lake water from each grid box downstream !@auth Gary Russell/Gavin Schmidt !@ver 1.0 (based on LB265) USE CONSTANT, only : shw,rhow,teeny,bygrav,tf USE MODEL_COM, only : im,jm,focean,zatmo,hlake,itlake,itlkice * ,itocean,itoice,fland,dtsrc USE DOMAIN_DECOMP, only : HALO_UPDATE, GRID,NORTH,SOUTH,GET, * GLOBALSUM, HALO_UPDATE_COLUMN USE GEOM, only : dxyp,bydxyp,imaxj USE DIAG_COM, only : aij=>aij_loc,ij_ervr,ij_mrvr,ij_f0oc, * aj=>aj_loc,aregj=>aregj_loc,jreg,j_rvrd,j_ervr,ij_fwoc USE GHY_COM, only : fearth #ifdef TRACERS_WATER USE TRDIAG_COM, only : taijn =>taijn_loc , tij_rvr USE FLUXES, only : trflowo,gtracer #endif USE FLUXES, only : flowo,eflowo,gtemp,mlhc,gtempr USE LAKES, only : kdirec,rate,iflow,jflow,river_fac USE LAKES_COM, only : tlake,gml,mwl,mldlk,flake #ifdef TRACERS_WATER * ,trlake,ntm #endif USE SEAICE_COM, only : rsi IMPLICIT NONE INTEGER :: FROM,J_0,J_1,J_0H,J_1H,J_0S,J_1S,I_0H,I_1H !@var I,J,IU,JU,ID,JD loop variables INTEGER I,J,IU,JU,ID,JD,JR,ITYPE REAL*8 MWLSILL,DMM,DGM,HLK1,DPE,MWLSILLD,FLFAC REAL*8, DIMENSION(IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: * FLOW,EFLOW !@var URATE upstream fractional rate of river flow per time step !@+ (only for special case) REAL*8 :: URATE = 1d-6 ! roughly 10 day e-folding time #ifdef TRACERS_WATER REAL*8, DIMENSION(NTM) :: DTM REAL*8, DIMENSION(NTM,IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO) * :: TRFLOW #endif LOGICAL :: rvrfl C**** C**** LAKECB MWL Liquid lake mass (kg) C**** GML Liquid lake enthalpy (J) C**** TLAKE Lake surface temperature (C) C**** C**** Calculate net mass and energy changes due to river flow C**** CALL GET(grid, J_STRT=J_0, J_STOP=J_1, & J_STRT_SKP =J_0S, J_STOP_SKP =J_1S, & J_STRT_HALO=J_0H, J_STOP_HALO=J_1H) FLOW = 0. ; EFLOW = 0. FLOWO = 0. ; EFLOWO = 0. #ifdef TRACERS_WATER TRFLOW = 0. TRFLOWO = 0. #endif CALL HALO_UPDATE(GRID, FLAND,FROM=NORTH+SOUTH) ! fixed CALL HALO_UPDATE(GRID,FOCEAN,FROM=NORTH+SOUTH) ! fixed CALL HALO_UPDATE(GRID,FEARTH,FROM=NORTH+SOUTH) CALL HALO_UPDATE(GRID, ZATMO,FROM=NORTH+SOUTH) ! fixed CALL HALO_UPDATE(GRID, HLAKE,FROM=NORTH+SOUTH) CALL HALO_UPDATE(grid, FLAKE,FROM=NORTH+SOUTH) CALL HALO_UPDATE(grid, MWL,FROM=NORTH+SOUTH) CALL HALO_UPDATE(grid, TLAKE,FROM=NORTH+SOUTH) CALL HALO_UPDATE(grid, RATE,FROM=NORTH+SOUTH) ! fixed #ifdef TRACERS_WATER CALL HALO_UPDATE_COLUMN(grid, GTRACER(:,1,:,:),NORTH+SOUTH) CALL HALO_UPDATE_COLUMN(grid, TRLAKE(:,1,:,:),NORTH+SOUTH) CALL HALO_UPDATE(grid, TAIJN(:,:,TIJ_RVR,:),FROM=NORTH+SOUTH) #endif C**** Calculate fluxes downstream if lake mass is above sill height (HLAKE (m)) C**** Also allow flow into ocean fraction of same box if KDIREC=0 C**** SPECIAL CASE: If the downstream box has FLAKE=1 and KDIREC=0 (i.e. C**** no outlet) then the only way to prevent excess water build up is C**** to allow a back-flux. Take account of mean topography change as C**** well. This is mainly an issue for the Caspian and Aral Seas. C**** Loop now includes polar boxes ! note on MPI fixes: since different PEs can influence the downstream ! accumulation of FLOW etc, we loop on the haloed variables to ensure ! that contributions from the halo are included in FLOW/FLOWO etc. ! If downstream box is outside the interior, cycle - this is dealt with on ! a separate PE DO JU=MAX(1,J_0H),MIN(JM,J_1H) DO IU=1,IMAXJ(JU) IF (KDIREC(IU,JU).gt.0 .or. * (FLAND(IU,JU).gt.0 .and. FOCEAN(IU,JU).gt.0)) THEN JD=JFLOW(IU,JU) ID=IFLOW(IU,JU) ! only calculate for downstream interior boxes. IF (JD.gt.J_1H .or. JD.lt.J_0H ) CYCLE C**** MWLSILL/D mass associated with full lake (and downstream) MWLSILL = RHOW*HLAKE(IU,JU)*FLAKE(IU,JU)*DXYP(JU) rvrfl=.false. C**** Check for special case: IF (KDIREC(ID,JD).eq.0 .and. FLAKE(ID,JD).ge. & 0.95d0*(FLAKE(ID,JD)+FEARTH(ID,JD))) THEN MWLSILLD = RHOW*(HLAKE(ID,JD)+BYGRAV*MAX(ZATMO(IU,JU) * -ZATMO(ID,JD),0d0))*DXYP(JD)*FLAKE(ID,JD) IF (MWL(ID,JD)-MWLSILLD.gt.0) THEN ! potential back flux IF (FLAKE(IU,JU).gt.0) THEN IF (MWL(ID,JD)-MWLSILLD-FLAKE(ID,JD)*DXYP(JD)*(MWL(IU,JU * )-MWLSILL)/(FLAKE(IU,JU)*DXYP(JU)).gt.0) THEN DMM=-(MWL(ID,JD)-MWLSILLD-FLAKE(ID,JD)*DXYP(JD) * *(MWL(IU,JU)-MWLSILL)/(FLAKE(IU,JU)*DXYP(JU))) * *URATE*DTsrc ! <0 rvrfl=.true. END IF ELSE DMM=-(MWL(ID,JD)-MWLSILLD)*URATE*DTsrc ! <0 rvrfl=.true. END IF if (rvrfl) then DGM=TLAKE(ID,JD)*DMM*SHW ! TLAKE always defined #ifdef TRACERS_WATER DTM(:) = DMM*GTRACER(:,1,ID,JD) #endif end if END IF END IF C**** Normal downstream flow IF(.not.rvrfl .and. MWL(IU,JU).gt.MWLSILL) THEN rvrfl=.true. DMM = (MWL(IU,JU)-MWLSILL)*RATE(IU,JU) IF (MWL(IU,JU)-DMM.lt.1d-6) DMM=MWL(IU,JU) DMM=MIN(DMM,0.5d0*RHOW*DXYP(JU)) ! minimise 'flood' events! c IF (FLAKE(IU,JU).gt.0) THEN c MLM=RHOW*MLDLK(IU,JU)*FLAKE(IU,JU)*DXYP(JU) c DMM=MIN(DMM,MLM) ! not necessary since MLM>TOTD-HLAKE c END IF DGM=TLAKE(IU,JU)*DMM*SHW ! TLAKE always defined #ifdef TRACERS_WATER if (flake(iu,ju).gt.0) then DTM(:) = DMM*GTRACER(:,1,IU,JU) else DTM(:) = DMM*TRLAKE(:,1,IU,JU)/MWL(IU,JU) end if #endif END IF IF (rvrfl) THEN FLOW(IU,JU) = FLOW(IU,JU) - DMM EFLOW(IU,JU) = EFLOW(IU,JU) - DGM #ifdef TRACERS_WATER TRFLOW(:,IU,JU) = TRFLOW(:,IU,JU) - DTM(:) #endif C**** calculate adjustments for poles IF (JU.eq.1 .or. JU.eq.JM) THEN FLFAC=IM ! pole exception upstream ELSEIF (JD.eq.1 .or. JD.eq.JM) THEN FLFAC=1d0/real(IM) ! pole exception downstream ELSE FLFAC=1. ! default END IF IF(FOCEAN(ID,JD).le.0.) THEN DPE=0. ! DMM*(ZATMO(IU,JU)-ZATMO(ID,JD)) FLOW(ID,JD) = FLOW(ID,JD) + DMM *FLFAC EFLOW(ID,JD) = EFLOW(ID,JD) + (DGM+DPE)*FLFAC #ifdef TRACERS_WATER TRFLOW(:,ID,JD)=TRFLOW(:,ID,JD) +DTM(:)*FLFAC #endif ELSE ! Save river mouth flow to for output to oceans C**** DPE: also add potential energy change to ocean. C**** Normally ocean is at sea level (Duh!), but in some boxes ZATMO C**** may not be zero if there is land as well, while in the Caspian, C**** the ocean level is below zero. C**** Note: this is diasabled until PE of precip is properly calculated C**** in atmosphere as well. Otherwise, there is an energy imbalance. DPE=0. ! DMM*(ZATMO(IU,JU)-MIN(0d0,ZATMO(ID,JD))) C**** possibly adjust mass (not heat) to allow for balancing of sea level DMM=river_fac*DMM FLOWO(ID,JD) = FLOWO(ID,JD)+ DMM *FLFAC EFLOWO(ID,JD)=EFLOWO(ID,JD)+(DGM+DPE)*FLFAC #ifdef TRACERS_WATER DTM(:)=river_fac*DTM(:) TRFLOWO(:,ID,JD)=TRFLOWO(:,ID,JD)+DTM(:)*FLFAC #endif C**** accumulate river runoff diags (moved from ground) AJ(JD,J_RVRD,ITOCEAN)=AJ(JD,J_RVRD,ITOCEAN)+ * (1.-RSI(ID,JD))*DMM*BYDXYP(JD) AJ(JD,J_ERVR,ITOCEAN)=AJ(JD,J_ERVR,ITOCEAN)+ * (1.-RSI(ID,JD))*(DGM+DPE)*BYDXYP(JD) AJ(JD,J_RVRD,ITOICE)=AJ(JD,J_RVRD,ITOICE) + * RSI(ID,JD)*DMM*BYDXYP(JD) AJ(JD,J_ERVR,ITOICE)=AJ(JD,J_ERVR,ITOICE) + * RSI(ID,JD)*(DGM+DPE)*BYDXYP(JD) AIJ(ID,JD,IJ_F0OC)=AIJ(ID,JD,IJ_F0OC)+ * (DGM+DPE)*BYDXYP(JD) AIJ(ID,JD,IJ_FWOC)=AIJ(ID,JD,IJ_FWOC)+DMM*BYDXYP(JD) END IF JR=JREG(ID,JD) AREGJ(JR,JD,J_RVRD)=AREGJ(JR,JD,J_RVRD)+DMM AREGJ(JR,JD,J_ERVR)=AREGJ(JR,JD,J_ERVR)+DGM+DPE AIJ(ID,JD,IJ_MRVR)=AIJ(ID,JD,IJ_MRVR) + DMM AIJ(ID,JD,IJ_ERVR)=AIJ(ID,JD,IJ_ERVR) + DGM+DPE #ifdef TRACERS_WATER TAIJN(ID,JD,TIJ_RVR,:)=TAIJN(ID,JD,TIJ_RVR,:)+DTM(:) * *BYDXYP(JD) #endif END IF END IF END DO END DO C**** C**** Apply net river flow to continental reservoirs C**** DO J=J_0, J_1 DO I=1,IMAXJ(J) IF(FLAND(I,J)+FLAKE(I,J).gt.0.) THEN MWL(I,J) = MWL(I,J) + FLOW(I,J) GML(I,J) = GML(I,J) + EFLOW(I,J) #ifdef TRACERS_WATER TRLAKE(:,1,I,J) = TRLAKE(:,1,I,J) + TRFLOW(:,I,J) #endif C**** remove pathologically small values IF (MWL(I,J).lt.1d-6) THEN MWL(I,J)=0. GML(I,J)=0. #ifdef TRACERS_WATER TRLAKE(:,1:2,I,J) = 0. #endif END IF IF (FLAKE(I,J).gt.0) THEN HLK1=(MLDLK(I,J)*RHOW)*TLAKE(I,J)*SHW MLDLK(I,J)=MLDLK(I,J)+FLOW(I,J)/(RHOW*FLAKE(I,J)*DXYP(J)) TLAKE(I,J)=(HLK1*FLAKE(I,J)*DXYP(J)+EFLOW(I,J)) * /(MLDLK(I,J)*RHOW*FLAKE(I,J)*DXYP(J)*SHW) C**** accumulate some diagnostics AJ(J,J_RVRD,ITLAKE) =AJ(J,J_RVRD,ITLAKE) + FLOW(I,J)* * BYDXYP(J)*(1.-RSI(I,J)) AJ(J,J_ERVR,ITLAKE) =AJ(J,J_ERVR,ITLAKE) +EFLOW(I,J)* * BYDXYP(J)*(1.-RSI(I,J)) AJ(J,J_RVRD,ITLKICE)=AJ(J,J_RVRD,ITLKICE)+ FLOW(I,J)* * BYDXYP(J)*RSI(I,J) AJ(J,J_ERVR,ITLKICE)=AJ(J,J_ERVR,ITLKICE)+EFLOW(I,J)* * BYDXYP(J)*RSI(I,J) ELSE TLAKE(I,J)=GML(I,J)/(SHW*MWL(I,J)+teeny) C**** accounting fix to ensure river flow with no lakes is counted AJ(J,J_RVRD,ITLAKE)=AJ(J,J_RVRD,ITLAKE)+ FLOW(I,J) * *BYDXYP(J) AJ(J,J_ERVR,ITLAKE)=AJ(J,J_ERVR,ITLAKE)+EFLOW(I,J) * *BYDXYP(J) END IF END IF END DO END DO CALL PRINTLK("RV") C**** Set GTEMP array for lakes DO J=J_0, J_1 DO I=1,IM IF (FLAKE(I,J).gt.0) THEN GTEMP(1,1,I,J)=TLAKE(I,J) GTEMPR(1,I,J) =TLAKE(I,J)+TF #ifdef TRACERS_WATER GTRACER(:,1,I,J)=TRLAKE(:,1,I,J)/(MLDLK(I,J)*RHOW*FLAKE(I,J) * *DXYP(J)) #endif MLHC(I,J) = SHW*MLDLK(I,J)*RHOW END IF END DO END DO RETURN C**** END SUBROUTINE RIVERF SUBROUTINE diag_RIVER 1,13 !@sum diag_RIVER prints out the river outflow for various rivers !@auth Gavin Schmidt !@ver 1.0 USE CONSTANT, only : rhow,sday,teeny,undef USE MODEL_COM, only : jyear0,amon0,jdate0,jhour0,jyear,amon * ,jdate,jhour,itime,dtsrc,idacc,itime0,nday,jdpery,jmpery USE DOMAIN_DECOMP, only : HALO_UPDATE, GRID,NORTH,SOUTH, * WRITE_PARALLEL USE GEOM, only : bydxyp USE DIAG_COM, only : aij,ij_mrvr #ifdef TRACERS_WATER USE TRACER_COM, only : ntm,trname,trw0,n_water,itime_tr0 * ,tr_wd_type,nwater USE TRDIAG_COM, only : taijn USE TRDIAG_COM, only : tij_rvr,to_per_mil,units_tij,scale_tij #endif USE LAKES, only : irvrmth,jrvrmth,namervr,nrvr IMPLICIT NONE REAL*8 RVROUT(6), SCALERVR, DAYS INTEGER INM,I,N #ifdef TRACERS_WATER REAL*8 TRRVOUT(6,NTM) #endif !@var out_line local variable to hold mixed-type output for parallel I/O character(len=300) :: out_line DAYS=(Itime-Itime0)/REAL(nday,kind=8) WRITE(out_line,900) JYEAR0,AMON0,JDATE0,JHOUR0,JYEAR,AMON,JDATE, * JHOUR,ITIME,DAYS CALL WRITE_PARALLEL(trim(out_line), UNIT=6) C**** convert kg/(source time step) to km^3/mon SCALERVR = 1d-9*SDAY*JDPERY/(JMPERY*RHOW*DTSRC) DO INM=1,NRVR,6 DO I=1,MIN(6,NRVR+1-INM) RVROUT(I) = SCALERVR*AIJ(IRVRMTH(I-1+INM),JRVRMTH(I-1+INM) * ,IJ_MRVR)/IDACC(1) END DO WRITE(out_line,901) * (NAMERVR(I-1+INM),RVROUT(I),I=1,MIN(6,NRVR+1-INM)) CALL WRITE_PARALLEL(trim(out_line), UNIT=6) END DO #ifdef TRACERS_WATER DO N=1,NTM if (itime.ge.itime_tr0(n) .and. tr_wd_TYPE(n).eq.nWater) then WRITE(out_line,*) "River outflow tracer concentration " * ,trim(units_tij(tij_rvr,n)),":",TRNAME(N) CALL WRITE_PARALLEL(trim(out_line), UNIT=6) DO INM=1,NRVR,6 DO I=1,MIN(6,NRVR+1-INM) IF (AIJ(IRVRMTH(I-1+INM),JRVRMTH(I-1+INM),IJ_MRVR).gt.0) * THEN if (to_per_mil(n).gt.0) then c TRRVOUT(I,N)=1d3*(TAIJN(IRVRMTH(I-1+INM),JRVRMTH(I-1 c * +INM),TIJ_RVR,N)/(trw0(n)*AIJ(IRVRMTH(I-1+INM) c * ,JRVRMTH(I-1+INM),IJ_MRVR)*BYDXYP(JRVRMTH(I-1+INM c * ))) -1.) if (TAIJN(IRVRMTH(I-1+INM),JRVRMTH(I-1+INM),TIJ_RVR * ,N_water).gt.0) then TRRVOUT(I,N)=1d3*(TAIJN(IRVRMTH(I-1+INM),JRVRMTH(I-1 * +INM),TIJ_RVR,N)/(trw0(n)*TAIJN(IRVRMTH(I-1+INM) * ,JRVRMTH(I-1+INM),TIJ_RVR,N_water))-1.) else TRRVOUT(I,N)=undef endif else TRRVOUT(I,N)=scale_tij(TIJ_RVR,n)*TAIJN(IRVRMTH(I-1+INM) * ,JRVRMTH(I-1+INM),TIJ_RVR,N)/(AIJ(IRVRMTH(I-1+INM) * ,JRVRMTH(I-1+INM),IJ_MRVR)*BYDXYP(JRVRMTH(I-1+INM)) * +teeny) end if ELSE TRRVOUT(I,N)=undef END IF END DO WRITE(out_line,901) (NAMERVR(I-1+INM),TRRVOUT(I,N), * I=1,MIN(6,NRVR+1-INM)) CALL WRITE_PARALLEL(trim(out_line), UNIT=6) END DO end if END DO #endif RETURN C**** 900 FORMAT ('1* River Outflow (km^3/mon) ** From:',I6,A6,I2,', Hr' * ,I3,6X,'To:',I6,A6,I2,', Hr',I3,' Model-Time:',I9,5X * ,'Dif:',F7.2,' Days') 901 FORMAT (' ',A8,':',F8.3,5X,A8,':',F8.3,5X,A8,':',F8.3,5X, * A8,':',F8.3,5X,A8,':',F8.3,5X,A8,':',F8.3) END SUBROUTINE diag_RIVER SUBROUTINE CHECKL (SUBR) 1,13 !@sum CHECKL checks whether the lake variables are reasonable. !@auth Gavin Schmidt/Gary Russell !@ver 1.0 (based on LB265) USE CONSTANT, only : rhow USE MODEL_COM, only : im,jm,hlake,qcheck USE DOMAIN_DECOMP, only : HALO_UPDATE, GET, GRID,NORTH,SOUTH USE GEOM, only : dxyp,imaxj #ifdef TRACERS_WATER USE TRACER_COM, only : ntm, trname, t_qlimit #endif USE LAKES USE LAKES_COM USE GHY_COM, only : fearth IMPLICIT NONE INTEGER :: FROM,J_0,J_1,J_0H,J_1H,J_0S,J_1S,I_0H,I_1H INTEGER I,J,N !@var I,J loop variables CHARACTER*6, INTENT(IN) :: SUBR LOGICAL QCHECKL #ifdef TRACERS_WATER integer :: imax,jmax real*8 relerr,errmax #endif CALL GET(grid, J_STRT=J_0, J_STOP=J_1, & J_STRT_SKP=J_0S, J_STOP_SKP=J_1S) C**** Check for NaN/INF in lake data CALL CHECK3(MWL,IM,JM,1,SUBR,'mwl') CALL CHECK3(GML,IM,JM,1,SUBR,'gml') CALL CHECK3(MLDLK,IM,JM,1,SUBR,'mld') CALL CHECK3(TLAKE,IM,JM,1,SUBR,'tlk') QCHECKL = .FALSE. DO J=J_0S, J_1S DO I=1,IM IF(FEARTH(I,J).gt.0.) THEN C**** check for negative mass IF (MWL(I,J).lt.0 .or. MLDLK(I,J).lt.0) THEN WRITE(6,*) 'After ',SUBR,': I,J,TSL,MWL,GML,MLD=', * I,J,TLAKE(I,J),MWL(I,J),GML(I,J),MLDLK(I,J) QCHECKL = .TRUE. END IF C**** check for reasonable lake surface temps IF (TLAKE(I,J).ge.50 .or. TLAKE(I,J).lt.-0.5) THEN WRITE(6,*) 'After ',SUBR,': I,J,TSL=',I,J,TLAKE(I,J) if (TLAKE(I,J).lt.-5.and.FLAKE(I,J).gt.0) QCHECKL = .TRUE. END IF END IF C**** Check total lake mass ( <0.4 m, >20x orig depth) IF(FLAKE(I,J).gt.0.) THEN IF(MWL(I,J).lt.0.4d0*RHOW*DXYP(J)*FLAKE(I,J)) THEN WRITE (6,*) 'After ',SUBR, * ': I,J,FLAKE,HLAKE,lake level low=',I,J,FLAKE(I,J), * HLAKE(I,J),MWL(I,J)/(RHOW*DXYP(J)*FLAKE(I,J)) END IF IF(MWL(I,J).gt.RHOW*MAX(20.*HLAKE(I,J),3d1)*DXYP(J)*FLAKE(I,J) * )THEN WRITE (6,*) 'After ',SUBR, * ': I,J,FLAKE,HLAKE,lake level high=',I,J,FLAKE(I,J), * HLAKE(I,J),MWL(I,J)/(RHOW*DXYP(J)*FLAKE(I,J)) END IF END IF END DO END DO #ifdef TRACERS_WATER do n=1,ntm C**** Check for neg tracers in lake if (t_qlimit(n)) then do j=J_0, J_1 do i=1,imaxj(j) if (fearth(i,j).gt.0) then if (trlake(n,1,i,j).lt.0 .or. trlake(n,2,i,j).lt.0) then print*,"Neg tracer in lake after ",SUBR,i,j,trname(n) * ,trlake(n,:,i,j) QCHECKL=.TRUE. end if end if end do end do end if C**** Check conservation of water tracers in lake if (trname(n).eq.'Water') then errmax = 0. ; imax=1 ; jmax=1 do j=J_0, J_1 do i=1,imaxj(j) if (fearth(i,j).gt.0) then if (flake(i,j).gt.0) then relerr=max( * abs(trlake(n,1,i,j)-mldlk(i,j)*rhow*flake(i,j) * *dxyp(j))/trlake(n,1,i,j),abs(trlake(n,1,i,j) * +trlake(n,2,i,j)-mwl(i,j))/(trlake(n,1,i,j) * +trlake(n,2,i,j))) else if ((mwl(i,j).eq.0 .and. trlake(n,1,i,j)+trlake(n,2,i,j) * .gt.0) .or. (mwl(i,j).gt.0 .and. trlake(n,1,i,j) * +trlake(n,2,i,j).eq.0)) then print*,"CHECKL ",SUBR,i,j,mwl(i,j),trlake(n,1:2,i,j) relerr=0. else if (mwl(i,j).gt.1d-20) then relerr=abs(trlake(n,1,i,j) * +trlake(n,2,i,j)-mwl(i,j))/(trlake(n,1,i,j) * +trlake(n,2,i,j)) else if (mwl(i,j).gt.0) print*,"CHECKL2 ",SUBR,i,j,mwl(i * ,j),trlake(n,1:2,i,j) relerr=0. end if end if end if if (relerr.gt.errmax) then imax=i ; jmax=j ; errmax=relerr end if end if end do end do print*,"Relative error in lake mass after ",trim(subr),":" * ,imax,jmax,errmax,trlake(n,:,imax,jmax),mldlk(imax,jmax) * *rhow*flake(imax,jmax)*dxyp(jmax),mwl(imax,jmax) * -mldlk(imax,jmax)*rhow*flake(imax,jmax)*dxyp(jmax) end if end do #endif IF (QCHECKL) & call stop_model('CHECKL: Lake variables out of bounds',255) RETURN C**** END SUBROUTINE CHECKL SUBROUTINE daily_LAKE 1,14 !@sum daily_LAKE does lake things at the beginning of every day !@auth G. Schmidt !@ver 1.0 USE CONSTANT, only : rhow,by3,pi,lhm,shi,shw,teeny,tf USE MODEL_COM, only : im,fland,flice,focean,itlake,itlkice,hlake USE LAKES, only : kdirec,minmld,variable_lk USE LAKES_COM, only : mwl,flake,tanlk,mldlk,tlake,gml,svflake #ifdef TRACERS_WATER * ,trlake,ntm #endif USE SEAICE_COM, only : rsi,msi,hsi,snowi #ifdef TRACERS_WATER * ,trsi #endif USE SEAICE, only : ace1i,xsi,ac2oim USE GEOM, only : dxyp,imaxj,bydxyp USE GHY_COM, only : fearth USE FLUXES, only : dmwldf,dgml,gtemp,mlhc,gtempr #ifdef TRACERS_WATER * ,dtrl,gtracer #endif USE DIAG_COM, only : aj=>aj_loc,j_run,j_erun,j_imelt,j_hmelt, * aregj=>aregj_loc,jreg USE DOMAIN_DECOMP, only : HALO_UPDATE, GET, GRID,NORTH,SOUTH, * GLOBALSUM IMPLICIT NONE integer i,j,J_0,J_1,jr,itm real*8 new_flake,sumh,msinew,snownew,frac,fmsi2,fmsi3 * ,fmsi4,fhsi2,fhsi3,fhsi4,imlt,hmlt,plake,plkic,hlk * ,frsat,new_MLD,hlkic #ifdef TRACERS_WATER * ,hlk2,ftsi2(ntm),ftsi3(ntm),ftsi4(ntm),sumt,dtr(ntm) & ,tottr(ntm) #endif CALL GET(grid, J_STRT=J_0, J_STOP=J_1) C**** Experimental code: not yet fully functional C**** Update lake fraction as a function of lake mass at end of day C**** Assume lake is conical C**** => A = pi*(h*tanlk)^2, M=(1/3)*pi*rho*h*(h*tanlk)^2 C**** SVFLAKE=FLAKE ! save for ghy purposes if (variable_lk .ne. 0) then DO J=J_0, J_1 DO I=1,IMAXJ(J) JR=JREG(I,J) IF (FLAKE(I,J)+FEARTH(I,J).gt.0 .and.FOCEAN(I,J).eq.0) THEN PLAKE=FLAKE(I,J)*(1.-RSI(I,J)) PLKIC=FLAKE(I,J)* RSI(I,J) new_flake=min(0.95d0*(FLAKE(I,J)+FEARTH(I,J)),(9d0*PI * *(TANLK(I,J)*MWL(I,J)/RHOW)**2)**BY3/DXYP(J)) C**** hack to prevent lakes fooding the snow in GHY C**** (do not flood more than 4.9% of land per day) new_flake=min( new_flake, FLAKE(I,J)+.049d0*FEARTH(I,J) ) hlk=0. hlkic=0. if (new_flake.gt.0) then hlk=MWL(I,J)/(RHOW*new_flake*DXYP(J)) hlkic=(MSI(I,J)+SNOWI(I,J)+ACE1I)*RSI(I,J)/RHOW end if if (new_flake.ne.FLAKE(I,J)) THEN ! something to do IF (FLAKE(I,J).eq.0) HLAKE(I,J)=MAX(1d0,HLAKE(I,J)) IF (new_flake.gt.0 .and. (hlk.gt.0.4 .or. (hlk.gt.0.2 .and * . hlk+hlkic.gt.0.4)) ) THEN ! new or surviving lake HLAKE(I,J)=MAX(HLAKE(I,J),1d0) ! in case it wasn't set C**** adjust for fearth changes FRSAT=0. IF (new_flake.gt.FLAKE(I,J)) THEN ! some water used to saturate if (MWL(I,J).gt.DMWLDF(I,J)*(new_flake * -FLAKE(I,J))*DXYP(J)) THEN FRSAT=DMWLDF(I,J)*(new_flake-FLAKE(I,J))*DXYP(J) * /MWL(I,J) MWL(I,J)=MWL(I,J)*(1.-FRSAT) C**** calculate associated energy/tracer transfer DGML(I,J)=FRSAT*GML(I,J) GML(I,J)=GML(I,J)*(1.-FRSAT) #ifdef TRACERS_WATER DTRL(:,I,J)=FRSAT*(TRLAKE(:,1,I,J)+TRLAKE(:,2,I,J)) TRLAKE(:,:,I,J)=TRLAKE(:,:,I,J)*(1.-FRSAT) #endif MLDLK(I,J)=MLDLK(I,J)*(1.-FRSAT) C**** save some diags AJ(J, J_RUN,ITLAKE) =AJ(J, J_RUN,ITLAKE) +PLAKE* * DMWLDF(I,J)*(new_flake-FLAKE(I,J)) AJ(J, J_RUN,ITLKICE)=AJ(J, J_RUN,ITLKICE)+PLKIC* * DMWLDF(I,J)*(new_flake-FLAKE(I,J)) AJ(J,J_ERUN,ITLAKE) =AJ(J,J_ERUN,ITLAKE) +PLAKE* * DGML(I,J)*BYDXYP(J) AJ(J,J_ERUN,ITLKICE)=AJ(J,J_ERUN,ITLKICE)+PLKIC* * DGML(I,J)*BYDXYP(J) else C**** this is just here to see whether this ever happens. print*,"dont saturate",i,j,(DMWLDF(I,J)*(new_flake * -FLAKE(I,J))*DXYP(J))/MWL(I,J),MWL(I,J) * ,(DMWLDF(I,J)*(new_flake-FLAKE(I,J))*DXYP(J)) * ,(new_flake-FLAKE(I,J)) call stop_model('Not enough water for saturation' * ,255) end if END IF C**** conserve lake ice IF (RSI(I,J)*FLAKE(I,J).gt.new_flake) THEN ! crunch ice up SUMH=PLKIC*SUM(HSI(:,I,J)) FRAC=PLKIC/new_flake SNOWNEW=SNOWI(I,J)*FRAC MSINEW=(MSI(I,J)+ACE1I)*FRAC-ACE1I RSI(I,J)=1. C**** all tracers --> tracer*FRAC, then adjust layering FMSI3=ACE1I*(FRAC-1d0) ! kg/m2 flux over new fraction FMSI2=FMSI3*XSI(1) FMSI4=FMSI3*XSI(4) FHSI2=FMSI2*HSI(1,I,J)/(XSI(1)*(ACE1I+SNOWI(I,J))) IF (FMSI3.LT.FRAC*XSI(2)*(ACE1I+SNOWI(I,J))) THEN FHSI3=FMSI3*HSI(2,I,J)/(XSI(2)*(ACE1I+SNOWI(I,J))) ELSE FHSI3=HSI(2,I,J)*FRAC+(FMSI3-FRAC*XSI(2)*(ACE1I * +SNOWI(I,J)))*HSI(1,I,J)/(XSI(1)*(ACE1I * +SNOWI(I,J))) END IF IF (FMSI4.LT.FRAC*XSI(3)*MSI(I,J)) THEN FHSI4=FMSI4*HSI(3,I,J)/(XSI(3)*MSI(I,J)) ELSE FHSI4=HSI(3,I,J)*FRAC+(FMSI4-FRAC*XSI(3)*MSI(I,J)) * *FHSI3/FMSI3 END IF HSI(1,I,J)=HSI(1,I,J)*(ACE1I+SNOWNEW)/ * (ACE1I+SNOWI(I,J)) HSI(2,I,J)=HSI(2,I,J)*FRAC+FHSI2-FHSI3 HSI(3,I,J)=HSI(3,I,J)*FRAC+FHSI3-FHSI4 HSI(4,I,J)=HSI(4,I,J)*FRAC +FHSI4 #ifdef TRACERS_WATER sumt=rsi(i,j)*flake(I,j)*sum(trsi(1,:,i,j)) FTSI2(:)=FMSI2*TRSI(:,1,I,J)/(XSI(1)*(ACE1I+SNOWI(I,J) * )) IF (FMSI3.LT.FRAC*XSI(2)*(ACE1I+SNOWI(I,J))) THEN FTSI3(:)=FMSI3*TRSI(:,2,I,J)/(XSI(2)*(ACE1I+SNOWI(I * ,J))) ELSE FTSI3(:)=TRSI(:,2,I,J)*FRAC+(FMSI3-FRAC*XSI(2) * *(ACE1I+SNOWI(I,J)))*TRSI(:,1,I,J)/(XSI(1) * *(ACE1I+SNOWI(I,J))) END IF IF (FMSI4.LT.FRAC*XSI(3)*MSI(I,J)) THEN FTSI4(:)=FMSI4*TRSI(:,3,I,J)/(XSI(3)*MSI(I,J)) ELSE FTSI4(:)=TRSI(:,3,I,J)*FRAC+(FMSI4-FRAC*XSI(3)*MSI(I * ,J))*FTSI3(:)/FMSI3 END IF TRSI(:,1,I,J)=TRSI(:,1,I,J)*(ACE1I+SNOWNEW)/ * (ACE1I+SNOWI(I,J)) TRSI(:,2,I,J)=TRSI(:,2,I,J)*FRAC+FTSI2(:)-FTSI3(:) TRSI(:,3,I,J)=TRSI(:,3,I,J)*FRAC+FTSI3(:)-FTSI4(:) TRSI(:,4,I,J)=TRSI(:,4,I,J)*FRAC +FTSI4(:) #endif MSI(I,J)=MSINEW SNOWI(I,J)=SNOWNEW ELSE RSI(I,J)=PLKIC/new_flake END IF C**** adjust layering if necessary HLK=MWL(I,J)/(RHOW*new_flake*DXYP(J)) new_MLD=MIN(MAX(MINMLD,HLK-HLAKE(I,J)),HLK) IF (MLDLK(I,J)*FLAKE(I,J).lt.new_flake*new_MLD) THEN IF (FLAKE(I,J).eq.0 .or. HLK.le.new_MLD) THEN ! new or shallow lake MLDLK(I,J)=new_MLD #ifdef TRACERS_WATER TOTTR(:)=TRLAKE(:,1,I,J)+TRLAKE(:,2,I,J) TRLAKE(:,2,I,J)=TOTTR(:)*(HLK-MLDLK(I,J))/HLK TRLAKE(:,1,I,J)=TOTTR(:)* MLDLK(I,J) /HLK #endif ELSE #ifdef TRACERS_WATER DTR(:)=TRLAKE(:,2,I,J)*(new_flake*new_MLD-MLDLK(I * ,J)*FLAKE(I,J))/(HLK*new_flake-MLDLK(I,J) * *FLAKE(I,J)) TRLAKE(:,1,I,J)=TRLAKE(:,1,I,J)+DTR(:) TRLAKE(:,2,I,J)=TRLAKE(:,2,I,J)-DTR(:) #endif MLDLK(I,J)=new_MLD END IF ELSE MLDLK(I,J)=MLDLK(I,J)*FLAKE(I,J)/new_flake END IF C**** adjust land surface fractions FLAKE(I,J)=new_flake FLAND(I,J)=1.-FLAKE(I,J) FEARTH(I,J)=FLAND(I,J)-FLICE(I,J) ELSE C**** remove/do not create lakes that are too small IF (FLAKE(I,J).gt.0) THEN C**** transfer lake ice mass/energy for accounting purposes IMLT=ACE1I+MSI(I,J)+SNOWI(I,J) HMLT=SUM(HSI(:,I,J)) MWL(I,J)=MWL(I,J)+PLKIC*IMLT*DXYP(J) GML(I,J)=GML(I,J)+PLKIC*HMLT*DXYP(J) #ifdef TRACERS_WATER DO ITM=1,NTM TRLAKE(ITM,1,I,J)=TRLAKE(ITM,1,I,J)+TRLAKE(ITM,2,I,J * )+RSI(I,J)*FLAKE(I,J)*SUM(TRSI(ITM,:,I,J)) * *DXYP(J) TRSI(ITM,:,I,J)=0. END DO #endif C**** save some diags AJ(J,J_IMELT,ITLKICE)=AJ(J,J_IMELT,ITLKICE)+PLKIC*IMLT AJ(J,J_HMELT,ITLKICE)=AJ(J,J_IMELT,ITLKICE)+PLKIC*HMLT C**** Accumulate regional diagnostics AREGJ(JR,J,J_IMELT)=AREGJ(JR,J,J_IMELT)+PLKIC*IMLT * *DXYP(J) AREGJ(JR,J,J_HMELT)=AREGJ(JR,J,J_HMELT)+PLKIC*HMLT * *DXYP(J) C**** RSI(I,J)=0. SNOWI(I,J)=0. HSI(1:2,I,J)=-LHM*XSI(1:2)*ACE1I HSI(3:4,I,J)=-LHM*XSI(3:4)*AC2OIM MSI(I,J)=AC2OIM TLAKE(I,J)=GML(I,J)/(SHW*MWL(I,J)+teeny) MLDLK(I,J)=MINMLD FLAKE(I,J)=0. FLAND(I,J)=1. FEARTH(I,J)=FLAND(I,J)-FLICE(I,J) END IF END IF END IF END IF END DO END DO end if C**** CALL PRINTLK("DY") C**** Set GTEMP array for lakes DO J=J_0, J_1 DO I=1,IMAXJ(J) IF (FLAKE(I,J).gt.0) THEN GTEMP(1,1,I,J)=TLAKE(I,J) GTEMPR(1,I,J) =TLAKE(I,J)+TF #ifdef TRACERS_WATER GTRACER(:,1,I,J)=TRLAKE(:,1,I,J)/(MLDLK(I,J)*RHOW*FLAKE(I,J) * *DXYP(J)) #endif MLHC(I,J) = SHW*MLDLK(I,J)*RHOW END IF END DO END DO C**** RETURN END SUBROUTINE daily_LAKE SUBROUTINE PRECIP_LK 1,10 !@sum PRECIP_LK driver for applying precipitation/melt to lake fraction !@auth Gavin Schmidt !@ver 1.0 USE CONSTANT, only : rhow,shw,teeny,tf USE MODEL_COM, only : im,jm,flice,itlake,itlkice USE DOMAIN_DECOMP, only : HALO_UPDATE, GRID,GET,NORTH,SOUTH USE GEOM, only : imaxj,dxyp,bydxyp USE SEAICE_COM, only : rsi USE LAKES_COM, only : mwl,gml,tlake,mldlk,flake #ifdef TRACERS_WATER * ,trlake,ntm #endif USE FLUXES, only : runpsi,runoli,prec,eprec,gtemp,melti,emelti, * gtempr #ifdef TRACERS_WATER * ,trunpsi,trunoli,trprec,gtracer,trmelti #endif USE DIAG_COM, only : aj=>aj_loc,j_run,aij=>aij_loc,ij_lk IMPLICIT NONE REAL*8 PRCP,ENRGP,PLICE,PLKICE,RUN0,ERUN0,POLAKE,HLK1 INTEGER :: FROM,J_0,J_1,J_0H,J_1H,J_0S,J_1S,I_0H,I_1H INTEGER I,J,ITYPE #ifdef TRACERS_WATER REAL*8, DIMENSION(NTM) :: TRUN0 #endif CALL GET(grid, J_STRT=J_0, J_STOP=J_1, & J_STRT_SKP=J_0S, J_STOP_SKP=J_1S) CALL PRINTLK("PR") DO J=J_0, J_1 DO I=1,IMAXJ(J) IF (FLAKE(I,J)+FLICE(I,J).gt.0) THEN POLAKE=(1.-RSI(I,J))*FLAKE(I,J) PLKICE=RSI(I,J)*FLAKE(I,J) PLICE=FLICE(I,J) PRCP=PREC(I,J) ENRGP=EPREC(I,J) ! energy of precipitation C**** calculate fluxes over whole box RUN0 =POLAKE*PRCP + PLKICE* RUNPSI(I,J) + PLICE* RUNOLI(I,J) ERUN0=POLAKE*ENRGP ! PLKICE*ERUNPSI(I,J) + PLICE*ERUNOLI(I,J) =0 C**** simelt is given as kg, so divide by area IF (FLAKE(I,J).gt.0) THEN RUN0 =RUN0+ MELTI(I,J)*BYDXYP(J) ERUN0=ERUN0+EMELTI(I,J)*BYDXYP(J) END IF MWL(I,J) = MWL(I,J) + RUN0*DXYP(J) GML(I,J) = GML(I,J) + ERUN0*DXYP(J) #ifdef TRACERS_WATER TRUN0(:) = POLAKE*TRPREC(:,I,J)*BYDXYP(J) * + PLKICE*TRUNPSI(:,I,J) + PLICE *TRUNOLI(:,I,J) IF (FLAKE(I,J).gt.0) TRUN0(:)=TRUN0(:)+TRMELTI(:,I,J)*BYDXYP(J) TRLAKE(:,1,I,J)=TRLAKE(:,1,I,J) + TRUN0(:)*DXYP(J) #endif IF (FLAKE(I,J).gt.0) THEN HLK1=TLAKE(I,J)*MLDLK(I,J)*RHOW*SHW MLDLK(I,J)=MLDLK(I,J) + RUN0/(FLAKE(I,J)*RHOW) TLAKE(I,J)=(HLK1*FLAKE(I,J)+ERUN0)/(MLDLK(I,J)*FLAKE(I,J) * *RHOW*SHW) GTEMP(1,1,I,J)=TLAKE(I,J) GTEMPR(1,I,J) =TLAKE(I,J)+TF IF (MWL(I,J).gt.(1d-10+MLDLK(I,J))*RHOW*FLAKE(I,J)*DXYP(J)) * THEN GTEMP(2,1,I,J)=(GML(I,J)-TLAKE(I,J)*SHW*MLDLK(I,J)*RHOW * *FLAKE(I,J)*DXYP(J))/(SHW*(MWL(I,J)-MLDLK(I,J) * *RHOW*FLAKE(I,J)*DXYP(J))) ELSE GTEMP(2,1,I,J)=TLAKE(I,J) END IF #ifdef TRACERS_WATER GTRACER(:,1,I,J)=TRLAKE(:,1,I,J)/(MLDLK(I,J)*RHOW*FLAKE(I,J) * *DXYP(J)) #endif AJ(J,J_RUN,ITLAKE) =AJ(J,J_RUN,ITLAKE) -PLICE*RUNOLI(I,J) * *(1.-RSI(I,J)) AJ(J,J_RUN,ITLKICE)=AJ(J,J_RUN,ITLKICE)-PLICE*RUNOLI(I,J) * *RSI(I,J) ELSE TLAKE(I,J)=GML(I,J)/(MWL(I,J)*SHW+teeny) C**** accounting fix to ensure runoff with no lakes is counted C**** no regional diagnostics required AJ(J,J_RUN,ITLAKE) =AJ(J,J_RUN,ITLAKE)-PLICE*RUNOLI(I,J) END IF C**** save area diag AIJ(I,J,IJ_LK) = AIJ(I,J,IJ_LK) + FLAKE(I,J) END IF END DO END DO RETURN C**** END SUBROUTINE PRECIP_LK SUBROUTINE GROUND_LK 1,16 !@sum GROUND_LK driver for applying surface fluxes to lake fraction !@auth Gavin Schmidt !@ver 1.0 !@calls USE CONSTANT, only : rhow,shw,teeny,tf USE MODEL_COM, only : im,jm,flice,fland,hlake * ,dtsrc,itlake,itlkice USE DOMAIN_DECOMP, only : GRID, GET,GLOBALSUM, HALO_UPDATE, * NORTH,SOUTH USE GEOM, only : imaxj,dxyp USE FLUXES, only : runosi, erunosi, e0, evapor, dmsi, dhsi, dssi, * runoli, runoe, erunoe, solar, dmua, dmva, gtemp, gtempr #ifdef TRACERS_WATER * ,trunoli,trunoe,trevapor,dtrsi,trunosi,gtracer #ifdef TRACERS_DRYDEP * ,trdrydep #endif #endif USE SEAICE_COM, only : rsi USE DIAG_COM, only : aj=>aj_loc,aregj=>aregj_loc,jreg,j_wtr1 * ,j_wtr2,j_run,j_erun USE LAKES_COM, only : mwl,gml,tlake,mldlk,flake #ifdef TRACERS_WATER * ,trlake,ntm USE TRDIAG_COM,only: taijn=>taijn_loc , tij_lk1,tij_lk2 #endif USE LAKES, only : lkmix,lksourc,byzeta,minmld USE GHY_COM, only : fearth IMPLICIT NONE C**** grid box variables REAL*8 ROICE, POLAKE, PLKICE, PEARTH, PLICE !@var MLAKE,ELAKE mass and energy /m^2 for lake model layers REAL*8, DIMENSION(2) :: MLAKE,ELAKE C**** fluxes REAL*8 EVAPO, FIDT, FODT, RUN0, ERUN0, RUNLI, RUNE, ERUNE, * HLK1,TLK1,TLK2,TKE,SROX(2),FSR2, U2RHO C**** output from LKSOURC REAL*8 ENRGFO, ACEFO, ACEFI, ENRGFI #ifdef TRACERS_WATER REAL*8, DIMENSION(NTM) :: TRUN0,TRO,TRI,TREVAP,TOTTRL REAL*8, DIMENSION(NTM,2) :: TRLAKEL #endif INTEGER I,J,JR INTEGER :: J_0,J_1,J_0S,J_1S CALL GET(grid, J_STRT=J_0, J_STOP=J_1, & J_STRT_SKP=J_0S, J_STOP_SKP=J_1S) CALL PRINTLK("GR") DO J=J_0, J_1 DO I=1,IMAXJ(J) JR=JREG(I,J) ROICE=RSI(I,J) PLKICE=FLAKE(I,J)*ROICE POLAKE=FLAKE(I,J)*(1.-ROICE) C**** Add land ice and surface runoff to lake variables IF (FLAND(I,J).gt.0) THEN PLICE =FLICE(I,J) PEARTH=FEARTH(I,J) RUNLI=RUNOLI(I,J) RUNE =RUNOE(I,J) ERUNE=ERUNOE(I,J) C**** calculate flux over whole box RUN0 =RUNLI*PLICE + RUNE*PEARTH ERUN0= ERUNE*PEARTH MWL(I,J) = MWL(I,J) + RUN0*DXYP(J) GML(I,J) = GML(I,J) +ERUN0*DXYP(J) #ifdef TRACERS_WATER TRLAKE(:,1,I,J)=TRLAKE(:,1,I,J)+ * (TRUNOLI(:,I,J)*PLICE+TRUNOE(:,I,J)*PEARTH)*DXYP(J) #endif IF (FLAKE(I,J).gt.0) THEN HLK1=TLAKE(I,J)*MLDLK(I,J)*RHOW*SHW MLDLK(I,J)=MLDLK(I,J) + RUN0/(FLAKE(I,J)*RHOW) TLAKE(I,J)=(HLK1*FLAKE(I,J)+ERUN0)/(MLDLK(I,J)*FLAKE(I,J) * *RHOW*SHW) #ifdef TRACERS_WATER GTRACER(:,1,I,J)=TRLAKE(:,1,I,J)/(MLDLK(I,J)*RHOW*FLAKE(I,J) * *DXYP(J)) #endif AJ(J,J_RUN,ITLAKE)=AJ(J,J_RUN,ITLAKE)- * (RUNE*PEARTH+RUNLI*PLICE)*(1.-RSI(I,J)) AJ(J,J_RUN,ITLKICE)=AJ(J,J_RUN,ITLKICE)- * (RUNE*PEARTH+RUNLI*PLICE)*RSI(I,J) AJ(J,J_ERUN,ITLAKE )=AJ(J,J_ERUN,ITLAKE )-ERUNE*PEARTH* * (1.-RSI(I,J)) AJ(J,J_ERUN,ITLKICE)=AJ(J,J_ERUN,ITLKICE)-ERUNE*PEARTH* * RSI(I,J) ELSE TLAKE(I,J)=GML(I,J)/(MWL(I,J)*SHW+teeny) C**** accounting fix to ensure runoff with no lakes is counted C**** no regional diagnostics required AJ(J,J_RUN,ITLAKE)=AJ(J,J_RUN,ITLAKE)- * (RUNE*PEARTH+RUNLI*PLICE) AJ(J,J_ERUN,ITLAKE )=AJ(J,J_ERUN,ITLAKE )-ERUNE*PEARTH END IF END IF IF (FLAKE(I,J).gt.0) THEN TLK1 =TLAKE(I,J) EVAPO=EVAPOR(I,J,1) ! evap/dew over open lake (kg/m^2) FODT =E0(I,J,1) ! net heat over open lake (J/m^2) SROX(1)=SOLAR(1,I,J) ! solar radiation open lake (J/m^2) SROX(2)=SOLAR(3,I,J) ! solar radiation through ice (J/m^2) FSR2 =EXP(-MLDLK(I,J)*BYZETA) C**** get ice-ocean fluxes from sea ice routine (over ice fraction) RUN0 =RUNOSI(I,J) ! includes ACE2M + basal term FIDT =ERUNOSI(I,J) C**** calculate kg/m^2, J/m^2 from saved variables MLAKE(1)=MLDLK(I,J)*RHOW MLAKE(2)=MAX(MWL(I,J)/(FLAKE(I,J)*DXYP(J))-MLAKE(1),0d0) ELAKE(1)=TLK1*SHW*MLAKE(1) ELAKE(2)=GML(I,J)/(FLAKE(I,J)*DXYP(J))-ELAKE(1) #ifdef TRACERS_WATER TRLAKEL(:,:)=TRLAKE(:,:,I,J)/(FLAKE(I,J)*DXYP(J)) TRUN0(:)=TRUNOSI(:,I,J) TREVAP(:)=TREVAPOR(:,1,I,J) #ifdef TRACERS_DRYDEP * -trdrydep(:,1,i,j) #endif #endif IF (MLAKE(2).lt.1d-10) THEN MLAKE(1)=MLAKE(1)+MLAKE(2) MLAKE(2)=0. ELAKE(1)=ELAKE(1)+ELAKE(2) ELAKE(2)=0. #ifdef TRACERS_WATER TRLAKEL(:,1)=TRLAKEL(:,1)+TRLAKEL(:,2) TRLAKEL(:,2)=0. #endif END IF C**** Limit FSR2 in the case of thin second layer FSR2=MIN(FSR2,MLAKE(2)/(MLAKE(1)+MLAKE(2))) C**** Apply fluxes and calculate the amount of frazil ice formation CALL LKSOURC (I,J,ROICE,MLAKE,ELAKE,RUN0,FODT,FIDT,SROX,FSR2, #ifdef TRACERS_WATER * TRLAKEL,TRUN0,TREVAP,TRO,TRI, #endif * EVAPO,ENRGFO,ACEFO,ACEFI,ENRGFI) C**** Mixing and entrainment C**** Calculate turbulent kinetic energy for lake c U2rho=(1.-ROICE)*SQRT(DMUA(I,J,1)**2+DMVA(I,J,1)**2)/DTSRC c TKE=0.5 * (19.3)^(2/3) * U2rho /rhoair ! (m/s)^2 TKE=0. ! 3.6d0*U2rho/rhoair*MLAKE(1) ! (J/m^2) CALL LKMIX (MLAKE,ELAKE, #ifdef TRACERS_WATER * TRLAKEL, #endif * HLAKE(I,J),TKE,ROICE,DTSRC) C**** Resave prognostic variables MWL(I,J) =(MLAKE(1)+MLAKE(2))*(FLAKE(I,J)*DXYP(J)) GML(I,J) =(ELAKE(1)+ELAKE(2))*(FLAKE(I,J)*DXYP(J)) MLDLK(I,J)= MLAKE(1)/RHOW IF (MLAKE(2).eq.0.) MLDLK(I,J)=MIN(MINMLD,MLDLK(I,J)) TLAKE(I,J)= ELAKE(1)/(SHW*MLAKE(1)) IF (MLAKE(2).gt.0) THEN TLK2 = ELAKE(2)/(SHW*MLAKE(2)) ELSE TLK2 = TLAKE(I,J) END IF #ifdef TRACERS_WATER IF (MLAKE(2).eq.0. .and. MLAKE(1)-MLDLK(I,J)*RHOW.gt.1d-10) THEN TOTTRL(:)=TRLAKEL(:,1) TRLAKEL(:,2)=(MLAKE(1)-MLDLK(I,J)*RHOW)*TRLAKEL(:,1)/MLAKE(1) TRLAKEL(:,1)=TOTTRL(:)-TRLAKEL(:,2) END IF TRLAKE(:,:,I,J)=TRLAKEL(:,:)*(FLAKE(I,J)*DXYP(J)) GTRACER(:,1,I,J)=TRLAKEL(:,1)/(MLDLK(I,J)*RHOW) #endif GTEMP(1,1,I,J)=TLAKE(I,J) GTEMP(2,1,I,J)=TLK2 ! diagnostic only GTEMPR(1,I,J) =TLAKE(I,J)+TF C**** Open lake diagnostics AJ(J,J_WTR1, ITLAKE)=AJ(J,J_WTR1, ITLAKE)+MLAKE(1)*POLAKE AJ(J,J_WTR2, ITLAKE)=AJ(J,J_WTR2, ITLAKE)+MLAKE(2)*POLAKE C**** Ice-covered ocean diagnostics AJ(J,J_WTR1, ITLKICE)=AJ(J,J_WTR1, ITLKICE)+MLAKE(1)*PLKICE AJ(J,J_WTR2, ITLKICE)=AJ(J,J_WTR2, ITLKICE)+MLAKE(2)*PLKICE C**** regional diags AREGJ(JR,J,J_WTR1)=AREGJ(JR,J,J_WTR1)+ * MLAKE(1)*FLAKE(I,J)*DXYP(J) AREGJ(JR,J,J_WTR2)=AREGJ(JR,J,J_WTR2)+ * MLAKE(2)*FLAKE(I,J)*DXYP(J) #ifdef TRACERS_WATER C**** tracer diagnostics TAIJN(I,J,tij_lk1,:)=TAIJN(I,J,tij_lk1,:)+TRLAKEL(:,1) !*PLKICE? TAIJN(I,J,tij_lk2,:)=TAIJN(I,J,tij_lk2,:)+TRLAKEL(:,2) !*PLKICE? #endif C**** Store mass and energy fluxes for formation of sea ice DMSI(1,I,J)=ACEFO DMSI(2,I,J)=ACEFI DHSI(1,I,J)=ENRGFO DHSI(2,I,J)=ENRGFI DSSI(:,I,J)=0. ! always zero salinity #ifdef TRACERS_WATER DTRSI(:,1,I,J)=TRO(:) DTRSI(:,2,I,J)=TRI(:) #endif END IF END DO ! i loop END DO ! j loop CALL PRINTLK("G2") RETURN C**** END SUBROUTINE GROUND_LK SUBROUTINE PRINTLK(STR) 6,8 !@sum PRINTLK print out selected diagnostics from specified lakes !@auth Gavin Schmidt !@ver 1.0 USE CONSTANT, only : lhm,byshi,rhow,shw USE MODEL_COM, only : qcheck USE GEOM, only : dxyp USE LAKES_COM, only : tlake,mwl,mldlk,gml,flake USE SEAICE, only : xsi,ace1i,rhoi USE SEAICE_COM, only : rsi,hsi,msi,snowi USE DOMAIN_DECOMP, only : GRID, GET IMPLICIT NONE CHARACTER*2, INTENT(IN) :: STR INTEGER, PARAMETER :: NDIAG=4 INTEGER I,J,N, J_0, J_1 INTEGER, DIMENSION(NDIAG) :: IDIAG = (/41, 66, 12, 31/), * JDIAG = (/27, 15, 36, 41/) REAL*8 HLK2,TLK2, TSIL(4) IF (.NOT.QCHECK) RETURN CALL GET(grid, J_STRT=J_0, J_STOP=J_1) DO N=1,NDIAG I=IDIAG(N) J=JDIAG(N) if (J.lt. J_0 .or. J.gt. J_1) CYCLE IF (FLAKE(I,J).gt.0) THEN HLK2 = MWL(I,J)/(RHOW*FLAKE(I,J)*DXYP(J)) - MLDLK(I,J) IF (HLK2.gt.0) THEN TLK2 = (GML(I,J)/(SHW*RHOW*FLAKE(I,J)*DXYP(J)) - * TLAKE(I,J)*MLDLK(I,J))/HLK2 ELSE TLK2=0. END IF TSIL(:)=0. IF (RSI(I,J).gt.0) THEN TSIL(1:2) = (HSI(1:2,I,J)/(XSI(1:2)*(ACE1I+SNOWI(I,J)))+LHM) * *BYSHI TSIL(3:4) = (HSI(3:4,I,J)/(XSI(3:4)*MSI(I,J))+LHM)*BYSHI END IF WRITE(99,*) STR,I,J,FLAKE(I,J),TLAKE(I,J),TLK2,MLDLK(I,J),HLK2 * ,RSI(I,J),MSI(I,J)/RHOI,SNOWI(I,J)/RHOW,TSIL(1:4) ELSE WRITE(99,*) STR,I,J,TLAKE(I,J),MWL(I,J) END IF END DO RETURN END SUBROUTINE PRINTLK SUBROUTINE conserv_LKM(LKM),5 !@sum conserv_LKM calculates zonal lake mass !@auth Gary Russell/Gavin Schmidt !@ver 1.0 USE MODEL_COM, only : im,jm,fland,fim USE DOMAIN_DECOMP, only : GRID, GET USE GEOM, only : imaxj,bydxyp USE LAKES_COM, only : mwl,flake IMPLICIT NONE REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: LKM INTEGER :: I,J INTEGER :: J_0,J_1,J_0S,J_1S LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE CALL GET(grid, J_STRT=J_0, J_STOP=J_1, & J_STRT_SKP=J_0S, J_STOP_SKP=J_1S, & HAVE_SOUTH_POLE = HAVE_SOUTH_POLE, & HAVE_NORTH_POLE = HAVE_NORTH_POLE ) C**** C**** LAKE MASS (kg/m^2) C**** DO J=J_0, J_1 LKM(J)=0. DO I=1,IMAXJ(J) IF (FLAND(I,J)+FLAKE(I,J).gt.0) LKM(J)=LKM(J)+MWL(I,J) END DO LKM(J)=LKM(J)*BYDXYP(J) END DO IF (HAVE_SOUTH_POLE) LKM(1) =FIM*LKM(1) IF (HAVE_NORTH_POLE) LKM(JM)=FIM*LKM(JM) RETURN END SUBROUTINE conserv_LKM SUBROUTINE conserv_LKE(LKE),5 !@sum conserv_LKE calculates zonal lake energy !@auth Gary Russell/Gavin Schmidt !@ver 1.0 USE MODEL_COM, only : im,jm,zatmo,fim,fland USE DOMAIN_DECOMP, only : GRID, GET USE GEOM, only : imaxj,bydxyp USE LAKES_COM, only : gml,mwl,flake IMPLICIT NONE REAL*8, DIMENSION(GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: LKE INTEGER :: I,J INTEGER :: FROM,J_0,J_1,J_0S,J_1S LOGICAL :: HAVE_SOUTH_POLE, HAVE_NORTH_POLE CALL GET(grid, J_STRT=J_0, J_STOP=J_1, & J_STRT_SKP=J_0S, J_STOP_SKP=J_1S, & HAVE_SOUTH_POLE=HAVE_SOUTH_POLE, & HAVE_NORTH_POLE=HAVE_NORTH_POLE) C**** C**** LAKE ENERGY (J/m^2) (includes potential energy (DISABLED)) C**** DO J=J_0, J_1 LKE(J)=0. DO I=1,IMAXJ(J) IF (FLAND(I,J)+FLAKE(I,J).gt.0) LKE(J)=LKE(J)+GML(I,J) c * +ZATMO(I,J)*MWL(I,J) END DO LKE(J)=LKE(J)*BYDXYP(J) END DO IF (HAVE_SOUTH_POLE) LKE(1)=FIM*LKE(1) IF (HAVE_NORTH_POLE) LKE(JM)=FIM*LKE(JM) RETURN END SUBROUTINE conserv_LKE