MODULE LAKES 6,2
!@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

      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, DIMENSION(IM,JM) :: KDIREC
!@var RATE rate of river flow downslope (fraction)
!@var DHORZ horizontal distance to downstream box (m)
      REAL*8, DIMENSION(IM,JM) :: RATE,DHORZ
!@var IFLOW,JFLOW grid box indexes for downstream direction
      INTEGER IFLOW(IM,JM),JFLOW(IM,JM)
!@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

      CONTAINS


      SUBROUTINE LKSOURC (I0,J0,ROICE,MLAKE,ELAKE,RUN0,FODT,FIDT,SROX 1,1
     *     ,FSR2,

     *     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

!@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
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


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)

      ELSE
        DM2 = 0.
        DH2 = 0.

      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

      ELSE
        MLAKE(2)=0.
        ELAKE(2)=0.

      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


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)


      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) WRITE(6,*)
     *           "Minimum lake level reached: rsi,mlake,elake",i0,j0
     *           ,roice,mlake(1)/rhow,elake(1)
            ENRGF1  =ENRGF1+FH0
            ELAKE(1)=MLAKE(1)*TFL*SHW
          END IF
        END IF
      END IF


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


      RETURN
      END SUBROUTINE LKSOURC


      SUBROUTINE LKMIX(MLAKE,ELAKE, 1

     *     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

!@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)

C**** Test for static stability
        IF ((TMAXRHO-TLK1)*(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

        ELSE ! not unstable, implicitly diffuse heat + entrain
C**** reduce mixing if there is ice cover, no mixing if temperature
C**** gradient is negative and deep water is at max density.
          IF (TLK2.lt.TLK1 .and. TLK2.gt.TMAXRHO) THEN
            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

          END IF
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
C**** DRHO=RHO(TLK2)-RHO(TLK1)~=(TLK2-TLK1)*dRHOdT(TLK1)
C**** Assumes a parabolic density function going through MAXRHO at
C**** TMAXRHO, and RHO0 at T=0. (reasonable up to about 12 C)
            DRHO=(TLK2-TLK1)*2d0*BFAC*(TMAXRHO-TLK1)
            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

            ELSE                ! entire second layer is entrained
              MLAKE(1)=MLT
              MLAKE(2)=0
              ELAKE(1)=HLT
              ELAKE(2)=0

            END IF
          END IF
        END IF
      END IF
      RETURN
C****
      END SUBROUTINE LKMIX

      END MODULE LAKES


      SUBROUTINE init_LAKES(inilake,istart) 1,15
!@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,fearth,jday
      USE GEOM, only : dxyp,dxv,dyv,dxp,dyp,imaxj

      USE FLUXES, only : gtemp,mlhc
      USE SEAICE_COM, only : rsi
      USE PBLCOM, only : tsavg
      USE LAKES
      USE LAKES_COM
      USE DAGCOM, only : npts,icon_LKM,icon_LKE,title_con,conpt0
      IMPLICIT NONE

      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.
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****
      FLAKE = FLAKE0   ! this is only here until FLAKE is variable

C**** Ensure that HLAKE is a minimum of 1m for FLAKE>0
      DO J=1,JM
        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),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=1,JM
          DO I=1,IM
c            FLAKE(I,J) = FLAKE0(I,J)
            IF (FLAKE(I,J).gt.0) THEN
              IF (HLAKE(I,J).lt.1.) print*,
     *             "HLAKE too small for lake fraction",i,j,HLAKE(I,J)
     *             ,FLAKE(I,J)
              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

            ELSE
              TLAKE(I,J) = 0.
              MWL(I,J)   = 0.
              GML(I,J)   = 0.
              MLDLK(I,J) = MINMLD

            END IF
          END DO
        END DO
      END IF

C**** Set geometric variables
C**** TANLK=TAN(ALPHA) = R/H for a conical lake of equivalent volume
      DO J=1,JM
        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=1,JM
        DO I=1,IM
          IF (FLAKE(I,J).gt.0) THEN
            GTEMP(1,1,I,J)=TLAKE(I,J)
            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

            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 (6,*) 'River Direction file read: ',TITLEI
      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,'(        END DO
      END IF
 10   call closeunit (iu_RVR)
C**** Create integral direction array KDIREC from CDIREC
      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**** If land but no ocean, and no direction, print warning
        IF ((FEARTH(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),FEARTH(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
          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
          IF (FOCEAN(I,J).le.0) THEN
            WRITE(6,*) "Warning: Named river outlet must be in ocean",i
     *           ,j,NAMERVR(INM),FOCEAN(I,J),FLICE(I,J),FLAKE0(I,J)
     *           ,FEARTH(I,J) 
          END IF
        END IF
      END DO
      END DO
      NRVR=INM
C****
C**** From each box calculate the downstream river box
C****
      DO J=2,JM-1
        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
      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
C****
C**** Calculate river flow RATE (per source time step)
C****
      SPEED0= .35d0
      SPMIN = .15d0
      SPMAX = 5.
      DZDH1 = .00005
      DO JU=1,JM-1
        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)

      RETURN
C****
 910  FORMAT (A72)
 911  FORMAT (72A1)
      END SUBROUTINE init_LAKES


      SUBROUTINE RIVERF 1,9
!@sum  RIVERF transports lake water from each grid box downstream
!@auth Gary Russell/Gavin Schmidt
!@ver  1.0 (based on LB265)
      USE CONSTANT, only : grav,shw,rhow,teeny
      USE MODEL_COM, only : im,jm,focean,zatmo,hlake,itlake,itlkice
     *     ,itocean,itoice,fland
      USE GEOM, only : dxyp,bydxyp
      USE DAGCOM, only : aij,ij_ervr,ij_mrvr,ij_f0oc,aj,areg,jreg,j_rvrd
     *     ,j_ervr

      USE FLUXES, only : flowo,eflowo,gtemp,mlhc
      USE LAKES, only : kdirec,rate,iflow,jflow
      USE LAKES_COM, only : tlake,gml,mwl,mldlk,flake

      USE SEAICE_COM, only : rsi
      IMPLICIT NONE

!@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
      REAL*8, DIMENSION(IM,JM) :: FLOW,EFLOW

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****
      FLOW = 0. ; EFLOW = 0.
      FLOWO = 0. ; EFLOWO = 0.

      DO JU=2,JM-1
        DO IU=1,IM
C**** Also allow flow into ocean fraction of same box if KDIREC=0
          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)
C**** Only overflow if lake mass is above sill height (HLAKE (m))
            MWLSILL = RHOW*HLAKE(IU,JU)*FLAKE(IU,JU)*DXYP(JU)
            IF(MWL(IU,JU).gt.MWLSILL) THEN
              DMM = (MWL(IU,JU)-MWLSILL)*RATE(IU,JU)
              IF (MWL(IU,JU)-DMM.lt.1d-20) DMM=MWL(IU,JU)
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
              FLOW(IU,JU) =  FLOW(IU,JU) - DMM
              EFLOW(IU,JU) = EFLOW(IU,JU) - DGM

              IF(FOCEAN(ID,JD).le.0.) THEN
                DPE=0.  ! DMM*(ZATMO(IU,JU)-ZATMO(ID,JD))
                FLOW(ID,JD) =  FLOW(ID,JD) + DMM
                EFLOW(ID,JD) = EFLOW(ID,JD) + DGM+DPE

              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)))
                FLOWO(ID,JD)=FLOWO(ID,JD)+DMM
                EFLOWO(ID,JD)=EFLOWO(ID,JD)+DGM+DPE

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)+
     *               (1.-RSI(ID,JD))*(DGM+DPE)*BYDXYP(JD)
              END IF
              JR=JREG(ID,JD)
              AREG(JR,J_RVRD)=AREG(JR,J_RVRD)+DMM
              AREG(JR,J_ERVR)=AREG(JR,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

            END IF
          END IF
        END DO
      END DO
C****
C**** Calculate river flow at the South Pole
C****
       FLOW(1,1) =  FLOW(1,1)/IM
      EFLOW(1,1) = EFLOW(1,1)/IM

C**** Only overflow if lake mass is above sill height (HLAKE (m))
      MWLSILL = RHOW*HLAKE(1,1)*FLAKE(1,1)*DXYP(1)
      JD=JFLOW(1,1)
      ID=IFLOW(1,1)
      IF(MWL(1,1).gt.MWLSILL) THEN
        DMM = (MWL(1,1)-MWLSILL)*RATE(1,1)
c        IF (FLAKE(1,1).gt.0) THEN
c         MLM=RHOW*MLDLK(1,1)*FLAKE(1,1)*DXYP(1)
c         DMM=MIN(DMM,MLM)   ! not necessary since MLM>TOTD-HLAKE
c        END IF
        DGM=TLAKE(1,1)*DMM*SHW ! TLAKE always defined
        DPE=0.  ! DMM*(ZATMO(1,1)-ZATMO(ID,JD))
        FLOW(1,1) =  FLOW(1,1) - DMM
        EFLOW(1,1) = EFLOW(1,1) - DGM

        AIJ(ID,JD,IJ_MRVR)=AIJ(ID,JD,IJ_MRVR) + DMM
        AIJ(ID,JD,IJ_ERVR)=AIJ(ID,JD,IJ_ERVR) + DGM+DPE
         FLOW(ID,JD) =  FLOW(ID,JD) + IM*DMM
        EFLOW(ID,JD) = EFLOW(ID,JD) + IM*(DGM+DPE)

      END IF
C****
C**** Apply net river flow to continental reservoirs
C****
      DO J=2,JM-1
        DO I=1,IM
          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)

C**** remove pathologically small values
            IF (MWL(I,J).lt.1d-20) THEN
              MWL(I,J)=0.
              GML(I,J)=0.

            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
      MWL(1,1) = MWL(1,1) +  FLOW(1,1)
      GML(1,1) = GML(1,1) + EFLOW(1,1)

      IF (MWL(1,1).lt.1d-20) THEN
        MWL(1,1)=0.
        GML(1,1)=0.

      END IF
      IF (FLAKE(1,1).gt.0) THEN
        HLK1=(MLDLK(1,1)*RHOW)*TLAKE(1,1)*SHW
        MLDLK(1,1)=MLDLK(1,1)+FLOW(1,1)/(RHOW*FLAKE(1,1)*DXYP(1))
        TLAKE(1,1)=(HLK1*FLAKE(1,1)*DXYP(1)+EFLOW(1,1))
     *       /(MLDLK(1,1)*RHOW*FLAKE(1,1)*DXYP(1)*SHW)
C**** accumulate some diagnostics
        AJ(1,J_RVRD,ITLAKE) =AJ(1,J_RVRD,ITLAKE) + FLOW(1,1)*
     *             BYDXYP(1)*(1.-RSI(1,1))
        AJ(1,J_ERVR,ITLAKE) =AJ(1,J_ERVR,ITLAKE) +EFLOW(1,1)*
     *       BYDXYP(1)*(1.-RSI(1,1))
        AJ(1,J_RVRD,ITLKICE)=AJ(1,J_RVRD,ITLKICE)+ FLOW(1,1)*
     *       BYDXYP(1)*RSI(1,1)
        AJ(1,J_ERVR,ITLKICE)=AJ(1,J_ERVR,ITLKICE)+EFLOW(1,1)*
     *       BYDXYP(1)*RSI(1,1)
      ELSE
        TLAKE(1,1)=GML(1,1)/(SHW*MWL(1,1)+teeny)
        AJ(1,J_RVRD,ITLAKE)=AJ(1,J_RVRD,ITLAKE)+ FLOW(1,1)*BYDXYP(1)
        AJ(1,J_ERVR,ITLAKE)=AJ(1,J_ERVR,ITLAKE)+EFLOW(1,1)*BYDXYP(1)
      END IF

      CALL PRINTLK("RV")
C**** Set GTEMP array for lakes
      DO J=1,JM
        DO I=1,IM
          IF (FLAKE(I,J).gt.0) THEN
            GTEMP(1,1,I,J)=TLAKE(I,J)

            MLHC(I,J) = SHW*MLDLK(I,J)*RHOW
          END IF
        END DO
      END DO

      RETURN
C****
      END SUBROUTINE RIVERF


      SUBROUTINE diag_RIVER 1,5
!@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 GEOM, only : bydxyp
      USE DAGCOM, only : aij,ij_mrvr

      USE LAKES, only : irvrmth,jrvrmth,namervr,nrvr
      IMPLICIT NONE
      REAL*8 RVROUT(6), SCALERVR, DAYS
      INTEGER INM,I,N


      DAYS=(Itime-Itime0)/REAL(nday,kind=8)
      WRITE(6,900) JYEAR0,AMON0,JDATE0,JHOUR0,JYEAR,AMON,JDATE,JHOUR
     *     ,ITIME,DAYS
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(6,901) (NAMERVR(I-1+INM),RVROUT(I),I=1,MIN(6,NRVR+1-INM))
      END DO



      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,10
!@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,fearth,qcheck
      USE GEOM, only : dxyp,imaxj

      USE LAKES
      USE LAKES_COM
      IMPLICIT NONE

      INTEGER I,J,N !@var I,J loop variables
      CHARACTER*6, INTENT(IN) :: SUBR
      LOGICAL QCHECKL


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=2,JM-1
      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)
c            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



      IF (QCHECKL)
     &     call stop_model('CHECKL: Lake variables out of bounds',255)
      RETURN
C****
      END SUBROUTINE CHECKL


      SUBROUTINE daily_LAKE 1,1
!@sum  daily_LAKE does lake things at the beginning of every day
!@auth G. Schmidt
!@ver  1.0
      IMPLICIT NONE

C**** Experimental code: not yet 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****
c        PRINT*,"TEST FLAKE CHANGE"
c        DO J=1,JM
c        DO I=1,IMAXJ(J)
c          IF (FLAKE(I,J).gt.0) THEN
c            PRINT*,"tan(a)",I,J,TANLK(I,J)
c            PRINT*,FLAKE(I,J),(9d0*PI*(TANLK(I,J)*MWL(I,J)/RHOW)**2)
c     *           **BY3/DXYP(J)
C**** remove lakes that are too small ??
c            IF (FLAKE(I,J).lt.0.005) THEN
c          PRINT*,"Lake->0 (too small)",I,J,MWL(I,J),GML(I,J),MLDLK(I,J)
c          MLDLK(I,J)=MIN(MLDLK(I,J),MWL(I,J)/(RHOW*FLAKE(I,J)*DXYP(J)))
c              FLAKE(I,J) = 0.
c          END IF
c        END DO
c        END DO
C****
      CALL PRINTLK("DY")
C****
      RETURN
      END SUBROUTINE daily_LAKE


      SUBROUTINE PRECIP_LK 1,8
!@sum  PRECIP_LK driver for applying precipitation/melt to lake fraction
!@auth Gavin Schmidt
!@ver  1.0
      USE CONSTANT, only : rhow,shw,teeny
      USE MODEL_COM, only : im,jm,flice,itlake,itlkice
      USE GEOM, only : imaxj,dxyp,bydxyp
      USE SEAICE_COM, only : rsi
      USE LAKES_COM, only : mwl,gml,tlake,mldlk,flake

      USE FLUXES, only : runpsi,runoli,prec,eprec,gtemp,melti,emelti

      USE DAGCOM, only : aj,j_run
      IMPLICIT NONE

      REAL*8 PRCP,ENRGP,PLICE,PLKICE,RUN0,ERUN0,POLAKE,HLK1
      INTEGER I,J,ITYPE


      CALL PRINTLK("PR")

      DO J=1,JM
      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)


        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)
          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

          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

      END IF
      END DO
      END DO
      RETURN
C****
      END SUBROUTINE PRECIP_LK


      SUBROUTINE GROUND_LK 1,12
!@sum  GROUND_LK driver for applying surface fluxes to lake fraction
!@auth Gavin Schmidt
!@ver  1.0
!@calls
      USE CONSTANT, only : rhow,shw,teeny
      USE MODEL_COM, only : im,jm,flice,fland,hlake
     *     ,fearth,dtsrc,itlake,itlkice
      USE GEOM, only : imaxj,dxyp
      USE FLUXES, only : runosi, erunosi, e0, evapor, dmsi, dhsi, dssi,
     *     runoli, runoe, erunoe, solar, dmua, dmva, gtemp

      USE SEAICE_COM, only : rsi
      USE DAGCOM, only : aj,areg,jreg,j_wtr1,j_wtr2,j_run,j_erun
      USE LAKES_COM, only : mwl,gml,tlake,mldlk,flake

      USE LAKES, only : lkmix,lksourc,byzeta,minmld
      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

      INTEGER I,J,JR

      CALL PRINTLK("GR")

      DO J=1,JM
      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)

        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)

          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)

        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.

        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,

     *       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,

     *       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

        GTEMP(1,1,I,J)=TLAKE(I,J)
        GTEMP(2,1,I,J)=TLK2       ! diagnostic only

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
        AREG(JR,J_WTR1)=AREG(JR,J_WTR1)+MLAKE(1)*FLAKE(I,J)*DXYP(J)
        AREG(JR,J_WTR2)=AREG(JR,J_WTR2)+MLAKE(2)*FLAKE(I,J)*DXYP(J)


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

      END IF
      END DO
      END DO

      CALL PRINTLK("G2")

      RETURN
C****
      END SUBROUTINE GROUND_LK



      SUBROUTINE PRINTLK(STR) 6,6
!@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
      IMPLICIT NONE
      CHARACTER*2, INTENT(IN) :: STR
      INTEGER, PARAMETER :: NDIAG=1
      INTEGER I,J,N
      INTEGER, DIMENSION(NDIAG) :: IDIAG = (/46/),
     *                             JDIAG = (/32/)
      REAL*8 HLK2,TLK2, TSIL(4)

      IF (.NOT.QCHECK) RETURN

      DO N=1,NDIAG
        I=IDIAG(N)
        J=JDIAG(N)
        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)
      END DO

      RETURN
      END  SUBROUTINE PRINTLK


      SUBROUTINE conserv_LKM(LKM),3
!@sum  conserv_LKM calculates zonal lake mass
!@auth Gary Russell/Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : im,jm,fland,fearth,fim
      USE GEOM, only : imaxj,bydxyp
      USE LAKES_COM, only : mwl,flake
      IMPLICIT NONE
      REAL*8, DIMENSION(JM) :: LKM
      INTEGER :: I,J
C****
C**** LAKE MASS (kg/m^2)
C****
      DO J=1,JM
        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
      LKM(1) =FIM*LKM(1)
      LKM(JM)=FIM*LKM(JM)
      RETURN
      END SUBROUTINE conserv_LKM


      SUBROUTINE conserv_LKE(LKE),3
!@sum  conserv_LKE calculates zonal lake energy
!@auth Gary Russell/Gavin Schmidt
!@ver  1.0
      USE MODEL_COM, only : im,jm,zatmo,fim,fland,fearth
      USE GEOM, only : imaxj,bydxyp
      USE LAKES_COM, only : gml,mwl,flake
      IMPLICIT NONE
      REAL*8, DIMENSION(JM) :: LKE
      INTEGER :: I,J
      REAL*8 AREA
C****
C**** LAKE ENERGY (J/m^2) (includes potential energy (DISABLED))
C****
      DO J=1,JM
        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
      LKE(1) =FIM*LKE(1)
      LKE(JM)=FIM*LKE(JM)
      RETURN
      END SUBROUTINE conserv_LKE