SUBROUTINE SURFCE 1,30
!@sum SURFCE calculates the surface fluxes which include
!@+   sensible heat, evaporation, thermal radiation, and momentum
!@+   drag.  It also calculates instantaneous surface temperature,
!@+   surface specific humidity, and surface wind components.
!@auth Nobody will claim responsibilty
      USE CONSTANT, only : grav,rgas,kapa,sday,lhm,lhe,lhs,twopi
     *     ,sha,tf,rhow,rhoi,shv,shw,shi,rvap,stbo,bygrav,by6,byshi
     *     ,byrhoi,deltx,byrt3,teeny
      USE DOMAIN_DECOMP, only : GRID, GET, CHECKSUM, HALO_UPDATE, SOUTH
      USE DOMAIN_DECOMP, only : NORTH
      USE MODEL_COM, only : im,jm,lm,fim,dtsrc,nisurf,u,v,t,p,q
     *     ,idacc,dsig,jday,ndasf,jeq,fland,flice,focean
     *     ,fearth,nday,modrd,itime,jhour,sige,byim,itocean
     *     ,itoice,itlake,itlkice,itlandi,qcheck,UOdrag,jdate
      USE GEOM, only : dxyp,imaxj,bydxyp,idjj,idij,rapj,kmaxj,sinip
     *     ,cosip
      USE SOMTQ_COM, only : tmom,qmom,mz,nmom
      USE DYNAMICS, only : pmid,pk,pedn,pek,pdsig,plij,am,byam
      USE RADNCB, only : trhr,fsf,cosz1

C**** Interface to PBL
      USE SOCPBL, only : zgs,ZS1,TGV,TKV,QG_SAT,QG_AVER,HEMI,DTSURF,POLE
     &     ,US,VS,WS,WSM,WSH,TSV,QSRF,PSI,DBL,KHS,KQS !,PPBL ! ,KMS
     &     ,UG,VG,WG,WINT

      USE PBLCOM, only : ipbl,cmgs,chgs,cqgs,tsavg,dclev
      USE PBL_DRV, only : pbl,evap_max,fr_sat,uocean,vocean,psurf,trhr0

      USE DAGCOM, only : oa,aij,tdiurn,aj,areg,adiurn,ndiupt,jreg
     *     ,ij_tsli,ij_shdtli,ij_evhdt,ij_trhdt,ij_shdt,ij_trnfp0
     *     ,ij_srtr,ij_neth,ij_ws,ij_ts,ij_us,ij_vs,ij_taus,ij_tauus
     *     ,ij_tauvs,ij_qs,ij_tg1,ij_evap,ij_evapo,ij_tgo,ij_f0oc
     *     ,ij_f0oi,ij_evapi,ij_f0li,ij_evapli,j_evap,j_evhdt
     *     ,j_tsrf,j_shdt,j_trhdt,j_type,j_tg1,j_tg2,ijdd,idd_spr
     *     ,idd_pt5,idd_pt4,idd_pt3,idd_pt2,idd_pt1,idd_ts,idd_tg1
     *     ,idd_q5,idd_q4,idd_q3,idd_q2,idd_q1,idd_qs,idd_qg,idd_swg
     *     ,idd_lwg,idd_sh,idd_lh,idd_hz0,idd_ug,idd_vg,idd_wg,idd_us
     *     ,idd_vs,idd_ws,idd_cia,idd_cm,idd_ch,idd_cq,idd_eds,idd_dbl
     *     ,idd_ev,idd_ldc,idd_dcf,hdiurn,ij_pblht
      USE LANDICE, only : hc2li,z1e,z2li,hc1li
      USE LANDICE_COM, only : snowli
      USE SEAICE, only : xsi,z1i,ace1i,hc1i,alami,byrli,byrls,
     *     solar_ice_frac,tfrez
      USE SEAICE_COM, only : rsi,msi,snowi,flag_dsws
      USE LAKES_COM, only : mwl,mldlk,gml,flake
      USE LAKES, only : minmld
      USE FLUXES, only : dth1,dq1,e0,e1,evapor,runoe,erunoe,sss
     *     ,solar,dmua,dmva,gtemp,nstype,uflux1,vflux1,tflux1,qflux1
     *     ,uosurf,vosurf,uisurf,visurf

      USE SOIL_DRV, only: earth
      IMPLICIT NONE

      INTEGER I,J,K,KR,JR,NS,NSTEPS,MODDSF,MODDD,ITYPE,IH,IHM,IDTYPE,IM1
      REAL*8 PLAND,PLICE,POICE,POCEAN,PIJ,PS,P1K
     *     ,BETA,ELHX,ACE2,CDTERM,CDENOM,dF1dTG,HCG1,HCG2,EVHDT,F1DT
     *     ,CM,CH,CQ,BETAUP,EVHEAT,F0,F1,DSHDTG,DQGDTG
     *     ,DEVDTG,DTRDTG,DF0DTG,DFDTG,DTG,dSNdTG,dEVdQS,HSDEN,HSCON
     *     ,HSMUL,dHS,dQS,dT2,dTS,DQ1X,EVHDT0,EVAP,F0DT,FTEVAP,PWATER
     *     ,PSK,Q1,THV1,PTYPE,TG1,SRHEAT,SNOW,TG2
     *     ,SHDT,TRHDT,TG,TS,RHOSRF,RCDMWS,RCDHWS,RCDQWS,SHEAT,TRHEAT
     *     ,QSDEN,QSCON,QSMUL,T2DEN,T2CON,T2MUL,TGDEN,FQEVAP
     *     ,Z1BY6L,QZ1,EVAPLIM,F2,FSRI(2),HTLIM

      REAL*8 MA1, MSI1
      REAL*8, DIMENSION(NSTYPE,IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO) ::
     *                                                       TGRND,TGRN2
      REAL*8, PARAMETER :: qmin=1.d-12
      REAL*8, PARAMETER :: S1BYG1 = BYRT3,
     *     Z1IBYL=Z1I/ALAMI, Z2LI3L=Z2LI/(3.*ALAMI), Z1LIBYL=Z1E/ALAMI
      REAL*8 QSAT,DQSATDT
      REAL*8, DIMENSION(7,3,IM,GRID%J_STRT_HALO:GRID%J_STOP_HALO) :: 
     *                                                         AREGIJ
c


      INTEGER :: J_0, J_1, J_0H, J_1H
C****
C**** Extract useful local domain parameters from "grid"
C****
      CALL GET(grid, J_STRT_HALO=J_0H, J_STOP_HALO=J_1H, 
     *               J_STRT=J_0,        J_STOP=J_1)

      NSTEPS=NIsurf*ITime
      DTSURF=DTsrc/NIsurf
      IH=JHOUR+1
      IHM = IH+(JDATE-1)*24

C**** ZERO OUT ENERGY AND EVAPORATION FOR GROUND AND INITIALIZE TGRND
      DO J=J_0,J_1
      DO I=1,IM
        TGRND(2,I,J)=GTEMP(1,2,I,J)
        TGRND(3,I,J)=GTEMP(1,3,I,J)
C       TGRND(4,I,J)=GTEMP(1,4,I,J)
        TGRN2(2,I,J)=GTEMP(2,2,I,J)
        TGRN2(3,I,J)=GTEMP(2,3,I,J)
      END DO
      END DO

C**** Zero out fluxes summed over type and surface time step
      E0=0. ; E1=0. ; EVAPOR=0. ; RUNOE=0. ; ERUNOE=0.
      DMUA=0. ; DMVA=0. ; SOLAR=0.
      AREGIJ = 0.


C****
C**** OUTSIDE LOOP OVER TIME STEPS, EXECUTED NIsurf TIMES EVERY HOUR
C****
      DO NS=1,NIsurf
         MODDSF=MOD(NSTEPS+NS-1,NDASF*NIsurf+1)
         IF(MODDSF.EQ.0) IDACC(3)=IDACC(3)+1
         MODDD=MOD(1+ITime/NDAY+NS,NIsurf)   ! 1+ not really needed ??
C**** ZERO OUT FLUXES ACCUMULATED OVER SURFACE TYPES
         DTH1=0. ;  DQ1 =0. ;  uflux1=0. ; vflux1=0.


      call loadbl



      Call CHECKSUM(GRID, uosurf, 193, "SURFACE.f")
      Call CHECKSUM(GRID, vosurf, 194, "SURFACE.f")
      Call CHECKSUM(GRID, uisurf, 195, "SURFACE.f")
      Call CHECKSUM(GRID, visurf, 196, "SURFACE.f")
      Call CHECKSUM(GRID, u     , 197, "SURFACE.f")
      Call CHECKSUM(GRID, v     , 198, "SURFACE.f")

      Call HALO_UPDATE(GRID, uosurf, FROM=SOUTH)
      Call HALO_UPDATE(GRID, vosurf, FROM=SOUTH)
      Call HALO_UPDATE(GRID, uisurf, FROM=SOUTH)
      Call HALO_UPDATE(GRID, visurf, FROM=SOUTH)
      Call HALO_UPDATE(GRID, u     , FROM=NORTH)
      Call HALO_UPDATE(GRID, v     , FROM=NORTH)


C****
C**** OUTSIDE LOOP OVER J AND I, EXECUTED ONCE FOR EACH GRID POINT
C****
!$OMP   PARALLEL DO PRIVATE (ACE2, BETA,BETAUP,   CM,CH,CQ,
!$OMP*  CDTERM,CDENOM,DSHDTG,DQGDTG,DEVDTG,DTRDTG,
!$OMP*  DF0DTG,DFDTG,DTG,DQ1X,DF1DTG,DSNDTG,DEVDQS,
!$OMP*  DHS,DQS,DT2,DTS, EVAP,EVAPLIM,ELHX,EVHDT,EVHEAT,EVHDT0,
!$OMP*  F0DT,F1DT,F0,F1,F2,FSRI, HCG1,HCG2,HSDEN,HSCON,
!$OMP*  HSMUL,HTLIM,I,ITYPE,IDTYPE,IM1, J,K,
!$OMP*  KR, MA1,MSI1, PS,P1K,PLAND,PWATER,
!$OMP*  PLICE,PIJ,POICE,POCEAN,PTYPE,PSK, Q1,QSDEN,
!$OMP*  QSCON,QSMUL, RHOSRF,RCDMWS,RCDHWS,RCDQWS, SHEAT,SRHEAT,
!$OMP*  SNOW,SHDT, T2DEN,T2CON,T2MUL,TGDEN,TS,
!$OMP*  THV1,TG,TG1,TG2,TRHDT,TRHEAT,Z1BY6L

!$OMP*  )
!$OMP*  SCHEDULE(DYNAMIC,2)
C
      DO J=J_0,J_1
      HEMI=1.
      IF(J.LE.JM/2) HEMI=-1.
      POLE= (J.EQ.1 .or. J.EQ.JM)

      IM1=IM
      DO I=1,IMAXJ(J)

      EVAPLIM = 0. ; HTLIM=0.  ! need initialisation

C****
C**** DETERMINE SURFACE CONDITIONS
C****
      PLAND=FLAND(I,J)
      PWATER=1.-PLAND
      PLICE=FLICE(I,J)
      POICE=RSI(I,J)*PWATER
      POCEAN=PWATER-POICE
      PIJ=P(I,J)
      PS=PEDN(1,I,J)
      PSK=PEK(1,I,J)
      P1K=PK(1,I,J)
      Q1=Q(I,J,1)
      THV1=T(I,J,1)*(1.+Q1*deltx)
CCC   JR=JREG(I,J)
      MA1=AM(1,I,J) !@var MA1 mass of lowest atmospheric layer (kg/m^2)
c     MSUM = (PS*100.)/GRAV !@var MSUM total mass of atmosphere (kg/m^2)
c     PGK = (PS*100.)**KAPA
c     PKDN = (GRAV*(MSUM-MA1*0.25))**KAPA

C**** QUANTITIES ACCUMULATED HOURLY FOR DIAGDD
         IF(MODDD.EQ.0) THEN
         DO KR=1,NDIUPT
           IF(I.EQ.IJDD(1,KR).AND.J.EQ.IJDD(2,KR)) THEN
             ADIURN(IH,IDD_SPR,KR)=ADIURN(IH,IDD_SPR,KR)+PS
             ADIURN(IH,IDD_PT5,KR)=ADIURN(IH,IDD_PT5,KR)+PSK*T(I,J,5)
             ADIURN(IH,IDD_PT4,KR)=ADIURN(IH,IDD_PT4,KR)+PSK*T(I,J,4)
             ADIURN(IH,IDD_PT3,KR)=ADIURN(IH,IDD_PT3,KR)+PSK*T(I,J,3)
             ADIURN(IH,IDD_PT2,KR)=ADIURN(IH,IDD_PT2,KR)+PSK*T(I,J,2)
             ADIURN(IH,IDD_PT1,KR)=ADIURN(IH,IDD_PT1,KR)+PSK*T(I,J,1)
             ADIURN(IH,IDD_Q5,KR)=ADIURN(IH,IDD_Q5,KR)+Q(I,J,5)
             ADIURN(IH,IDD_Q4,KR)=ADIURN(IH,IDD_Q4,KR)+Q(I,J,4)
             ADIURN(IH,IDD_Q3,KR)=ADIURN(IH,IDD_Q3,KR)+Q(I,J,3)
             ADIURN(IH,IDD_Q2,KR)=ADIURN(IH,IDD_Q2,KR)+Q(I,J,2)
             ADIURN(IH,IDD_Q1,KR)=ADIURN(IH,IDD_Q1,KR)+Q1
             HDIURN(IHM,IDD_SPR,KR)=HDIURN(IHM,IDD_SPR,KR)+PS
             HDIURN(IHM,IDD_PT5,KR)=HDIURN(IHM,IDD_PT5,KR)+PSK*T(I,J,5)
             HDIURN(IHM,IDD_PT4,KR)=HDIURN(IHM,IDD_PT4,KR)+PSK*T(I,J,4)
             HDIURN(IHM,IDD_PT3,KR)=HDIURN(IHM,IDD_PT3,KR)+PSK*T(I,J,3)
             HDIURN(IHM,IDD_PT2,KR)=HDIURN(IHM,IDD_PT2,KR)+PSK*T(I,J,2)
             HDIURN(IHM,IDD_PT1,KR)=HDIURN(IHM,IDD_PT1,KR)+PSK*T(I,J,1)
             HDIURN(IHM,IDD_Q5,KR)=HDIURN(IHM,IDD_Q5,KR)+Q(I,J,5)
             HDIURN(IHM,IDD_Q4,KR)=HDIURN(IHM,IDD_Q4,KR)+Q(I,J,4)
             HDIURN(IHM,IDD_Q3,KR)=HDIURN(IHM,IDD_Q3,KR)+Q(I,J,3)
             HDIURN(IHM,IDD_Q2,KR)=HDIURN(IHM,IDD_Q2,KR)+Q(I,J,2)
             HDIURN(IHM,IDD_Q1,KR)=HDIURN(IHM,IDD_Q1,KR)+Q1
           END IF
         END DO
         END IF
C****
      DO ITYPE=1,3       ! no earth type
      ipbl(i,j,itype)=0
!!!      SELECT CASE (ITYPE)
C****
C**** OPEN OCEAN/LAKE
C****
!!!      CASE (1)
      if ( ITYPE == 1 ) then

      PTYPE=POCEAN
      IF (PTYPE.gt.0) THEN
      TG1=GTEMP(1,1,I,J)
      TG2=GTEMP(2,1,I,J)   ! diagnostic only
      IF (FLAKE(I,J).gt.0) THEN
C**** limit evap/cooling if between MINMLD and 40cm, no evap below 40cm
        IF (MWL(I,J).lt.MINMLD*RHOW*FLAKE(I,J)*DXYP(J)) THEN
          EVAPLIM=MAX(0.5*(MWL(I,J)/(FLAKE(I,J)*DXYP(J))-0.4d0*RHOW),
     *         0d0)
        ELSE
          EVAPLIM=MWL(I,J)/(FLAKE(I,J)*DXYP(J))-(0.5*MINMLD+0.2d0)*RHOW
        END IF

        HTLIM = GML(I,J)/(FLAKE(I,J)*DXYP(J)) + 0.5*LHM*EVAPLIM
        IDTYPE=ITLAKE
        uocean = 0. ; vocean = 0. ! no velocities for lakes
      ELSE
        IDTYPE=ITOCEAN
        if (UOdrag.eq.1) then ! use uocean for drag calculation
C**** Convert UOSURF,VOSURF from C grid to A grid
C**** Note that uosurf,vosurf start with j=1, (not j=2 as in atm winds)
          if (pole) then
            uocean = 0. ; vocean = 0.
            do k=1,kmaxj(j)
              uocean = uocean + rapj(k,j)*(uosurf(idij(k,i,j),idjj(k,j)
     *             -1)*cosip(k)-hemi*vosurf(idij(k,i,j),idjj(k,j)-1)
     *             *sinip(k))
              vocean = vocean + rapj(k,j)*(vosurf(idij(k,i,j),idjj(k,j)
     *             -1)*cosip(k)+hemi*uosurf(idij(k,i,j),idjj(k,j)-1)
     *             *sinip(k))
            end do
          else
            uocean = 0.5*(uosurf(i,j)+uosurf(im1,j))
            vocean = 0.5*(vosurf(i,j)+vosurf(i,j-1))
          end if
        else
          uocean=0. ; vocean=0.
        end if
      END IF
      SRHEAT=FSF(ITYPE,I,J)*COSZ1(I,J)
      SOLAR(1,I,J)=SOLAR(1,I,J)+DTSURF*SRHEAT
            OA(I,J,5)=OA(I,J,5)+SRHEAT*DTSURF
      BETA=1.
      ELHX=LHE
      END IF
C****
C**** OCEAN/LAKE ICE
C****
!!!      CASE (2)
      else if ( ITYPE == 2 ) then

      PTYPE=POICE
      IF (PTYPE.gt.0) THEN
      IF (FLAKE(I,J).gt.0) THEN
        IDTYPE=ITLKICE
        uocean = 0. ; vocean = 0. ! no dynamic ice for lakes
      ELSE
        IDTYPE=ITOICE
        if (UOdrag.eq.1) then ! use ice velcoities in drag calculation
C**** Convert UISURF,VISURF from C grid to A grid
C**** Note that uisurf,visurf start with j=1, (not j=2 as in atm winds)
          if (pole) then
            uocean = 0. ; vocean = 0.
            do k=1,kmaxj(j)
              uocean = uocean + rapj(k,j)*(uisurf(idij(k,i,j),idjj(k,j)
     *             -1)*cosip(k)-hemi*visurf(idij(k,i,j),idjj(k,j)-1)
     *             *sinip(k))
              vocean = vocean + rapj(k,j)*(visurf(idij(k,i,j),idjj(k,j)
     *             -1)*cosip(k)+hemi*uisurf(idij(k,i,j),idjj(k,j)-1)
     *             *sinip(k))
            end do
          else
            uocean = 0.5*(uisurf(i,j)+uisurf(im1,j))
            vocean = 0.5*(visurf(i,j)+visurf(i,j-1))
          end if
        else
          uocean = 0. ; vocean =0.
        end if
      END IF
      SNOW=SNOWI(I,J)
      TG1=TGRND(2,I,J)
      TG2=TGRN2(2,I,J)
      MSI1 = SNOW+ACE1I ! snow and first layer ice mass (kg/m^2)
      ACE2=MSI(I,J) ! second (physical) layer ice mass (kg/m^2)
      dF1dTG = 2./(ACE1I*BYRLI+SNOW*BYRLS)
      HCG1 = SHI*XSI(1)*MSI1 ! heat capacity of top ice layer (J/C*m^2)
      HCG2 = SHI*XSI(2)*MSI1 ! heat capacity of second layer ice
      SRHEAT=FSF(ITYPE,I,J)*COSZ1(I,J)
      SOLAR(2,I,J)=SOLAR(2,I,J)+DTSURF*SRHEAT
C**** fraction of solar radiation leaving layer 1 and 2
      IF (SRHEAT.gt.0) THEN ! only bother if there is sun
        call solar_ice_frac(SNOW,ACE2,FLAG_DSWS(I,J),FSRI,2)
      ELSE
        FSRI(1:2) = 0
      END IF
            OA(I,J,12)=OA(I,J,12)+SRHEAT*DTSURF
      BETA=1.
      ELHX=LHS

      Z1BY6L=(Z1LIBYL+SNOW*BYRLS)*BY6
      CDENOM=1./(2.*Z1BY6L+Z2LI3L)

      END IF
C****
C**** LAND ICE
C****
!!!      CASE (3)
      else if ( ITYPE == 3 ) then

      PTYPE=PLICE
      IF (PTYPE.gt.0) THEN
      IDTYPE=ITLANDI
      SNOW=SNOWLI(I,J)
      TG1=TGRND(3,I,J)
      TG2=GTEMP(2,3,I,J)
      SRHEAT=FSF(ITYPE,I,J)*COSZ1(I,J)
      Z1BY6L=(Z1LIBYL+SNOW*BYRLS)*BY6
      CDTERM=TG2
      CDENOM=1./(2.*Z1BY6L+Z2LI3L)
      HCG1=HC1LI+SNOW*SHI
      BETA=1.
      ELHX=LHS
      uocean = 0. ; vocean = 0. ! no land ice velocity
      END IF
!!!      END SELECT
      endif
C****
      IF (PTYPE.gt.0) THEN
C****
C**** BOUNDARY LAYER INTERACTION
C****
      TKV=THV1*PSK  ! TKV is referenced to the surface pressure
      ZS1=.5*DSIG(1)*RGAS*BYGRAV*TKV*PIJ/PMID(1,I,J)
      SHDT=0.
      EVHDT=0.
      TRHDT=0.
      F1DT=0.

      TG=TG1+TF
      QG_SAT=QSAT(TG,ELHX,PS)
      IF (ITYPE.eq.1 .and. focean(i,j).gt.0) QG_SAT=0.98d0*QG_SAT
      QG_AVER=QG_SAT
      TGV=TG*(1.+QG_SAT*deltx)
      psurf=PS   ! extra values to pass to PBL, possibly temporary
      trhr0 = TRHR(0,I,J)

C =====================================================================
      fr_sat = 1. ! entire surface is saturated
      evap_max = 1.
      CALL PBL(I,J,ITYPE,PTYPE)
      CM = cmgs(i,j,itype)
      CH = chgs(i,j,itype)
      CQ = cqgs(i,j,itype)
C =====================================================================
      TS=TSV/(1.+QSRF*deltx)
C**** CALCULATE RHOSRF*CM*WS AND RHOSRF*CH*WS
      RHOSRF=100.*PS/(RGAS*TSV)
      RCDMWS=CM*WSM*RHOSRF
      RCDHWS=CH*WSH*RHOSRF
      RCDQWS=CQ*WSH*RHOSRF
C**** CALCULATE FLUXES OF SENSIBLE HEAT, LATENT HEAT, THERMAL
C****   RADIATION, AND CONDUCTION HEAT (WATTS/M**2)
      SHEAT=SHA*RCDHWS*(TS-TG)
      BETAUP = BETA
      IF (QSRF .GT. QG_SAT) BETAUP = 1.
      EVHEAT=(LHE+TG1*SHV)*BETAUP*RCDQWS*(QSRF-QG_SAT)
      TRHEAT=TRHR(0,I,J)-STBO*(TG*TG)*(TG*TG)
C****
!!!      SELECT CASE (ITYPE)

!!!      CASE (1) ! FLUXES USING IMPLICIT TIME STEP FOR OCEAN POINTS
      if ( ITYPE == 1) then
c       DSHDTG=-RCDHWS*SHA
c       dEVdQS = LHE*RCDQWS
c       dHS = -(1.+2.*S1BYG1)*DTSURF*PGK*SHEAT/
c    A       ((1.+2.*S1BYG1)*DTSURF*PGK*RCDHWS+MA1*PKDN)
c       dTS = -(1.+2.*S1BYG1)*DTSURF*PGK*SHEAT/
c    A       (MA1*PKDN*SHA-(1.+2.*S1BYG1)*DTSURF*PGK*DSHDTG)
c       dQS = -(1.+2.*S1BYG1)*DTSURF*EVHEAT/
c    A       ((1.+2.*S1BYG1)*DTSURF*dEVdQS+MA1*LHE)
c       SHDT = DTSURF*(SHEAT-dTS*DSHDTG)
c       EVHDT=DTSURF*(EVHEAT+dQS*dEVdQS) ! latent heat flux
        SHDT = DTSURF*SHEAT
        EVHDT=DTSURF*EVHEAT              ! latent heat flux
        TRHDT=DTSURF*TRHEAT

C****
!!!      CASE (2) ! FLUXES USING IMPLICIT TIME STEP FOR ICE POINTS
      else if ( ITYPE == 2 ) then

! heat flux on first/second/third layers (W/m^2)
        F1 = (TG1-TG2)*dF1dTG + SRHEAT*FSRI(1)
        F2 = SRHEAT*FSRI(2)
        EVHEAT = LHE*RCDQWS*(QSRF-QG_SAT) ! latent heat flux (W/m^2)
        F0=SRHEAT+TRHEAT+SHEAT+EVHEAT
        dSNdTG=-RCDHWS*SHA
        dQGdTG=QG_SAT*DQSATDT(TG,ELHX) ! d(QG)/dTG
        dEVdTG = -dQGdTG*LHE*RCDQWS ! d(EVHEAT)/dTG
        dTRdTG = -4*STBO*TG*TG*TG ! d(TRHEAT)/dTG
        dF0dTG = dSNdTG+dEVdTG+dTRdTG ! d(F0)/dTG
cc      dSNdHS = RCDHWS ! d(SHEAT)/dHS - kg/(sec*m^2)
c       dEVdQS = LHE*RCDQWS     ! d(EVHEAT)/dQS
c       HSDEN = -(1.+2.*S1BYG1)*DTSURF*PGK*BETA*dSNdTG+MA1*PKDN*SHA
c       HSCON = -(1.+2.*S1BYG1)*DTSURF*PGK*SHEAT/HSDEN ! (J*sec)/kg
c       HSMUL=-(1.+2.*S1BYG1)*DTSURF*PGK*BETA*dSNdTG/HSDEN ! J/(kg*degC)
c       QSDEN = (1.+2.*S1BYG1)*BETA*DTSURF*dEVdQS+MA1*LHE
c       QSCON = -(1.+2.*S1BYG1)*DTSURF*EVHEAT/QSDEN
c       QSMUL = -(1.+2.*S1BYG1)*DTSURF*BETA*dEVdTG/QSDEN
        T2DEN = HCG2+BETA*DTSURF*dF1dTG
        T2CON = DTSURF*(F1-F2)/T2DEN
        T2MUL = BETA*DTSURF*dF1dTG/T2DEN
c       TGDEN = HCG1-BETA*DTSURF*(dF0dTG-dF1dTG-
c    A       HSMUL*dSNdTG+QSMUL*dEVdQS+T2MUL*dF1dTG) ! W/(m^2*degC)
c       dTG = DTSURF*(F0-F1+BETA*
c    A       (QSCON*dEVdQS-HSCON*dSNdTG+T2CON*dF1dTG))/TGDEN ! degC

        DFDTG=DF0DTG-(1.-DF0DTG*Z1BY6L)*CDENOM
        DTG=(F0-F1)*DTSURF/(HCG1-DTSURF*DFDTG)

        IF (TG1+dTG .GT. 0.) dTG = -TG1
c       dHS = HSCON+HSMUL*dTG   ! (J*sec)/kg
c       dQS = QSCON+QSMUL*dTG
        dT2 = T2CON+T2MUL*dTG
c       SHDT = DTSURF*(SHEAT+BETA*((dTG-dHS)*dSNdTG)) ! sensible
c       EVHDT = DTSURF*(EVHEAT+BETA*(dTG*dEVdTG+dQS*dEVdQS)) ! latent
        SHDT = DTSURF*(SHEAT+BETA*dTG*dSNdTG) ! sensible
        EVHDT = DTSURF*(EVHEAT+BETA*dTG*dEVdTG) ! latent

        TRHDT = DTSURF*(TRHEAT+BETA*dTG*dTRdTG) ! thermal flux (J/m^2)
        F1DT = DTSURF*(F1+BETA*(dTG*dF1dTG-dT2*dF1dTG))
        TG1 = TG1+dTG          ! first layer sea ice temperature (degC)
        TG2 = TG2+dT2          ! second layer sea ice temperature (degC)
        TGRN2(ITYPE,I,J) = TG2
C****
!!!      CASE (3) ! IMPLICIT TIME STEP OVER LANDICE
      else if ( ITYPE == 3 ) then

        F0=SRHEAT+TRHEAT+SHEAT+EVHEAT
        F1=(TG1-CDTERM-F0*Z1BY6L)*CDENOM
        DSHDTG=-RCDHWS*SHA
        DQGDTG=QG_SAT*DQSATDT(TG,ELHX)
        DEVDTG=-RCDQWS*LHE*BETAUP*DQGDTG
        DTRDTG=-4.*STBO*TG*TG*TG
        DF0DTG=DSHDTG+DEVDTG+DTRDTG
        DFDTG=DF0DTG-(1.-DF0DTG*Z1BY6L)*CDENOM
        DTG=(F0-F1)*DTSURF/(HCG1-DTSURF*DFDTG)
        SHDT=DTSURF*(SHEAT+DTG*DSHDTG)
        EVHDT=DTSURF*(EVHEAT+DTG*DEVDTG)
        TRHDT=DTSURF*(TRHEAT+DTG*DTRDTG)
        F1DT=DTSURF*(TG1-CDTERM-(F0+DTG*DFDTG)*Z1BY6L)*CDENOM
        TG1=TG1+DTG

!!!      END SELECT
      endif

C**** CALCULATE EVAPORATION
      DQ1X =EVHDT/((LHE+TG1*SHV)*MA1)
      EVHDT0=EVHDT
C**** Limit evaporation if lake mass is at minimum
      IF (ITYPE.EQ.1 .and. FLAKE(I,J).GT.0 .and.
     *     (EVAPOR(I,J,1)-DQ1X*MA1).gt.EVAPLIM) THEN
        if (QCHECK) WRITE(99,*) "Lake EVAP limited: I,J,EVAP,MWL",I,J
     *     ,EVAPOR(I,J,1)-DQ1X*MA1,MWL(I,J)/(RHOW*FLAKE(I,J)*DXYP(J))
        DQ1X=(EVAPOR(I,J,1)-EVAPLIM)*BYAM(1,I,J)
      ELSEIF (DQ1X.GT.Q1+DQ1(I,J)) THEN
        DQ1X=(Q1+DQ1(I,J))
      ELSE
        GO TO 3720
      END IF
      EVHDT=DQ1X*(LHE+TG1*SHV)*MA1
      IF (ITYPE.NE.1) TG1=TG1+(EVHDT-EVHDT0)/HCG1
 3720 EVAP=-DQ1X*MA1




C**** ACCUMULATE SURFACE FLUXES AND PROGNOSTIC AND DIAGNOSTIC QUANTITIES
      F0DT=DTSURF*SRHEAT+TRHDT+SHDT+EVHDT
C**** Limit heat fluxes out of lakes if near minimum depth
      IF (ITYPE.eq.1 .and. FLAKE(I,J).gt.0 .and.
     *     E0(I,J,1)+F0DT+HTLIM.lt.0) THEN
        if (QCHECK) write(6,*) "Limiting heat flux from lake",i,j,SHDT
     *       ,F0DT,E0(I,J,1),DTSURF*SRHEAT,TRHDT,EVHDT,HTLIM
        SHDT = -(HTLIM+E0(I,J,1)+DTSURF*SRHEAT+TRHDT+EVHDT)
        F0DT = -E0(I,J,1)-HTLIM
        if (QCHECK) write(6,*) "New SHDT,F0DT",i,j,SHDT,F0DT
      END IF
      E0(I,J,ITYPE)=E0(I,J,ITYPE)+F0DT
      E1(I,J,ITYPE)=E1(I,J,ITYPE)+F1DT
      EVAPOR(I,J,ITYPE)=EVAPOR(I,J,ITYPE)+EVAP
      TGRND(ITYPE,I,J)=TG1
      DTH1(I,J)=DTH1(I,J)-SHDT*PTYPE/(SHA*MA1*P1K)
      DQ1(I,J) =DQ1(I,J) -DQ1X*PTYPE
      DMUA(I,J,ITYPE)=DMUA(I,J,ITYPE)+PTYPE*DTSURF*RCDMWS*(US-UOCEAN)
      DMVA(I,J,ITYPE)=DMVA(I,J,ITYPE)+PTYPE*DTSURF*RCDMWS*(VS-VOCEAN)
      uflux1(i,j)=uflux1(i,j)+PTYPE*RCDMWS*(US-UOCEAN)
      vflux1(i,j)=vflux1(i,j)+PTYPE*RCDMWS*(VS-UOCEAN)
C****
C**** ACCUMULATE DIAGNOSTICS FOR EACH SURFACE TIME STEP AND ITYPE
C****
        AJ(J,J_EVAP ,IDTYPE)=AJ(J,J_EVAP ,IDTYPE)+ EVAP*PTYPE
        AJ(J,J_EVHDT,IDTYPE)=AJ(J,J_EVHDT,IDTYPE)+EVHDT*PTYPE
        AJ(J,J_SHDT ,IDTYPE)=AJ(J,J_SHDT ,IDTYPE)+ SHDT*PTYPE
        AJ(J,J_TRHDT,IDTYPE)=AJ(J,J_TRHDT,IDTYPE)+TRHDT*PTYPE
        IF(MODDSF.EQ.0) THEN
          AJ(J,J_TSRF,IDTYPE)=AJ(J,J_TSRF,IDTYPE)+(TS-TF)*PTYPE
          AJ(J,J_TYPE,IDTYPE)=AJ(J,J_TYPE,IDTYPE)+        PTYPE
          AJ(J,J_TG1 ,IDTYPE)=AJ(J,J_TG1 ,IDTYPE)+    TG1*PTYPE
          AJ(J,J_TG2 ,IDTYPE)=AJ(J,J_TG2 ,IDTYPE)+    TG2*PTYPE
        END IF
C**** QUANTITIES ACCUMULATED FOR REGIONS IN DIAGJ
CCC     AREG(JR,J_TRHDT)=AREG(JR,J_TRHDT)+TRHDT*PTYPE*DXYP(J)
CCC     AREG(JR,J_SHDT )=AREG(JR,J_SHDT )+SHDT *PTYPE*DXYP(J)
CCC     AREG(JR,J_EVHDT)=AREG(JR,J_EVHDT)+EVHDT*PTYPE*DXYP(J)
CCC     AREG(JR,J_EVAP )=AREG(JR,J_EVAP )+EVAP *PTYPE*DXYP(J)
CCC     IF(MODDSF.EQ.0)
CCC  *       AREG(JR,J_TSRF)=AREG(JR,J_TSRF)+(TS-TF)*PTYPE*DXYP(J)
        AREGIJ(1,ITYPE,I,J)=TRHDT*PTYPE*DXYP(J)
        AREGIJ(2,ITYPE,I,J)=SHDT *PTYPE*DXYP(J)
        AREGIJ(3,ITYPE,I,J)=EVHDT*PTYPE*DXYP(J)
        AREGIJ(4,ITYPE,I,J)=EVAP *PTYPE*DXYP(J)
        IF(MODDSF.EQ.0) THEN
          AREGIJ(5,ITYPE,I,J)=(TS-TF)*PTYPE*DXYP(J)
          AREGIJ(6,ITYPE,I,J)=    TG1*PTYPE*DXYP(J)
          AREGIJ(7,ITYPE,I,J)=    TG2*PTYPE*DXYP(J)
        END IF
C
C**** QUANTITIES ACCUMULATED FOR LATITUDE-LONGITUDE MAPS IN DIAGIJ
        AIJ(I,J,IJ_SHDT)=AIJ(I,J,IJ_SHDT)+SHDT*PTYPE
        IF(MODRD.EQ.0) AIJ(I,J,IJ_TRNFP0)=AIJ(I,J,IJ_TRNFP0)+TRHDT
     *       *PTYPE/DTSRC
        AIJ(I,J,IJ_SRTR)=AIJ(I,J,IJ_SRTR)+(SRHEAT*DTSURF+TRHDT)*PTYPE
        AIJ(I,J,IJ_NETH)=AIJ(I,J,IJ_NETH)+(SRHEAT*DTSURF+TRHDT+SHDT
     *       +EVHDT)*PTYPE
        AIJ(I,J,IJ_EVAP)=AIJ(I,J,IJ_EVAP)+EVAP*PTYPE
        IF(MODDSF.EQ.0) THEN
          AIJ(I,J,IJ_WS)=AIJ(I,J,IJ_WS)+WS*PTYPE
          AIJ(I,J,IJ_TS)=AIJ(I,J,IJ_TS)+(TS-TF)*PTYPE
          AIJ(I,J,IJ_US)=AIJ(I,J,IJ_US)+US*PTYPE
          AIJ(I,J,IJ_VS)=AIJ(I,J,IJ_VS)+VS*PTYPE
          AIJ(I,J,IJ_TAUS)=AIJ(I,J,IJ_TAUS)+RCDMWS*WSM*PTYPE
          AIJ(I,J,IJ_TAUUS)=AIJ(I,J,IJ_TAUUS)+RCDMWS*(US-UOCEAN)*PTYPE
          AIJ(I,J,IJ_TAUVS)=AIJ(I,J,IJ_TAUVS)+RCDMWS*(VS-VOCEAN)*PTYPE
          AIJ(I,J,IJ_QS)=AIJ(I,J,IJ_QS)+QSRF*PTYPE
          AIJ(I,J,IJ_TG1)=AIJ(I,J,IJ_TG1)+TG1*PTYPE
          AIJ(I,J,IJ_PBLHT)=AIJ(I,J,IJ_PBLHT)+dbl*PTYPE
        END IF
C**** QUANTITIES ACCUMULATED HOURLY FOR DIAGDD
        IF(MODDD.EQ.0) THEN
          DO KR=1,NDIUPT
            IF(I.EQ.IJDD(1,KR).AND.J.EQ.IJDD(2,KR)) THEN
              ADIURN(IH,IDD_TS,KR)=ADIURN(IH,IDD_TS,KR)+TS*PTYPE
              ADIURN(IH,IDD_TG1,KR)=ADIURN(IH,IDD_TG1,KR)+(TG1+TF)*PTYPE
              ADIURN(IH,IDD_QS,KR)=ADIURN(IH,IDD_QS,KR)+QSRF*PTYPE
              ADIURN(IH,IDD_QG,KR)=ADIURN(IH,IDD_QG,KR)+QG_SAT*PTYPE
              ADIURN(IH,IDD_SWG,KR)=ADIURN(IH,IDD_SWG,KR)+SRHEAT*DTSURF
     *             *PTYPE
              ADIURN(IH,IDD_LWG,KR)=ADIURN(IH,IDD_LWG,KR)+TRHDT*PTYPE
              ADIURN(IH,IDD_SH,KR)=ADIURN(IH,IDD_SH,KR)+SHDT*PTYPE
              ADIURN(IH,IDD_LH,KR)=ADIURN(IH,IDD_LH,KR)+EVHDT*PTYPE
              ADIURN(IH,IDD_HZ0,KR)=ADIURN(IH,IDD_HZ0,KR)
     *             +(SRHEAT*DTSURF+TRHDT+SHDT+EVHDT)*PTYPE
              ADIURN(IH,IDD_UG,KR)=ADIURN(IH,IDD_UG,KR)+UG*PTYPE
              ADIURN(IH,IDD_VG,KR)=ADIURN(IH,IDD_VG,KR)+VG*PTYPE
              ADIURN(IH,IDD_WG,KR)=ADIURN(IH,IDD_WG,KR)+WG*PTYPE
              ADIURN(IH,IDD_US,KR)=ADIURN(IH,IDD_US,KR)+US*PTYPE
              ADIURN(IH,IDD_VS,KR)=ADIURN(IH,IDD_VS,KR)+VS*PTYPE
              ADIURN(IH,IDD_WS,KR)=ADIURN(IH,IDD_WS,KR)+WS*PTYPE
              ADIURN(IH,IDD_CIA,KR)=ADIURN(IH,IDD_CIA,KR)+PSI*PTYPE
              ADIURN(IH,IDD_CM,KR)=ADIURN(IH,IDD_CM,KR)+CM*PTYPE
              ADIURN(IH,IDD_CH,KR)=ADIURN(IH,IDD_CH,KR)+CH*PTYPE
              ADIURN(IH,IDD_CQ,KR)=ADIURN(IH,IDD_CQ,KR)+CQ*PTYPE
              ADIURN(IH,IDD_EDS,KR)=ADIURN(IH,IDD_EDS,KR)+KHS*PTYPE
              ADIURN(IH,IDD_DBL,KR)=ADIURN(IH,IDD_DBL,KR)+DBL*PTYPE
              ADIURN(IH,IDD_EV,KR)=ADIURN(IH,IDD_EV,KR)+EVAP*PTYPE
              HDIURN(IHM,IDD_TS,KR)=HDIURN(IHM,IDD_TS,KR)+TS*PTYPE
              HDIURN(IHM,IDD_TG1,KR)=HDIURN(IHM,IDD_TG1,KR)+(TG1+TF)
     *             *PTYPE
              HDIURN(IHM,IDD_QS,KR)=HDIURN(IHM,IDD_QS,KR)+QSRF*PTYPE
              HDIURN(IHM,IDD_QG,KR)=HDIURN(IHM,IDD_QG,KR)+QG_SAT*PTYPE
              HDIURN(IHM,IDD_SWG,KR)=HDIURN(IHM,IDD_SWG,KR)+SRHEAT
     *             *DTSURF*PTYPE
              HDIURN(IHM,IDD_LWG,KR)=HDIURN(IHM,IDD_LWG,KR)+TRHDT*PTYPE
              HDIURN(IHM,IDD_SH,KR)=HDIURN(IHM,IDD_SH,KR)+SHDT*PTYPE
              HDIURN(IHM,IDD_LH,KR)=HDIURN(IHM,IDD_LH,KR)+EVHDT*PTYPE
              HDIURN(IHM,IDD_HZ0,KR)=HDIURN(IHM,IDD_HZ0,KR)
     *             +(SRHEAT*DTSURF+TRHDT+SHDT+EVHDT)*PTYPE
              HDIURN(IHM,IDD_UG,KR)=HDIURN(IHM,IDD_UG,KR)+UG*PTYPE
              HDIURN(IHM,IDD_VG,KR)=HDIURN(IHM,IDD_VG,KR)+VG*PTYPE
              HDIURN(IHM,IDD_WG,KR)=HDIURN(IHM,IDD_WG,KR)+WG*PTYPE
              HDIURN(IHM,IDD_US,KR)=HDIURN(IHM,IDD_US,KR)+US*PTYPE
              HDIURN(IHM,IDD_VS,KR)=HDIURN(IHM,IDD_VS,KR)+VS*PTYPE
              HDIURN(IHM,IDD_WS,KR)=HDIURN(IHM,IDD_WS,KR)+WS*PTYPE
              HDIURN(IHM,IDD_CIA,KR)=HDIURN(IHM,IDD_CIA,KR)+PSI*PTYPE
              HDIURN(IHM,IDD_CM,KR)=HDIURN(IHM,IDD_CM,KR)+CM*PTYPE
              HDIURN(IHM,IDD_CH,KR)=HDIURN(IHM,IDD_CH,KR)+CH*PTYPE
              HDIURN(IHM,IDD_CQ,KR)=HDIURN(IHM,IDD_CQ,KR)+CQ*PTYPE
              HDIURN(IHM,IDD_EDS,KR)=HDIURN(IHM,IDD_EDS,KR)+KHS*PTYPE
              HDIURN(IHM,IDD_DBL,KR)=HDIURN(IHM,IDD_DBL,KR)+DBL*PTYPE
              HDIURN(IHM,IDD_EV,KR)=HDIURN(IHM,IDD_EV,KR)+EVAP*PTYPE
            END IF
          END DO
        END IF
C****

C****
C**** SAVE SOME TYPE DEPENDENT FLUXES/DIAGNOSTICS
C****
!!!      SELECT CASE (ITYPE)
!!!      CASE (1)  ! ocean
      if ( ITYPE == 1 ) then
        OA(I,J,6)=OA(I,J,6)+TRHDT
        OA(I,J,7)=OA(I,J,7)+SHDT
        OA(I,J,8)=OA(I,J,8)+EVHDT
        IF (MODDSF.eq.0) AIJ(I,J,IJ_TGO)=AIJ(I,J,IJ_TGO)+TG1
        AIJ(I,J,IJ_EVAPO)=AIJ(I,J,IJ_EVAPO)+EVAP*PTYPE
        IF (FOCEAN(I,J).gt.0) AIJ(I,J,IJ_F0OC)=AIJ(I,J,IJ_F0OC) +F0DT
     *       *PTYPE
C****
!!!      CASE (2)  ! seaice
      else if ( ITYPE == 2 ) then
        IF (TG1.GT.TDIURN(I,J,7)) TDIURN(I,J,7) = TG1
        OA(I,J,9)=OA(I,J,9)+TRHDT
        OA(I,J,10)=OA(I,J,10)+SHDT
        OA(I,J,11)=OA(I,J,11)+EVHDT
        AIJ(I,J,IJ_F0OI) =AIJ(I,J,IJ_F0OI) +F0DT*PTYPE
        AIJ(I,J,IJ_EVAPI)=AIJ(I,J,IJ_EVAPI)+EVAP*PTYPE
C****
!!!      CASE (3) ! land ice
      else if ( ITYPE == 3 ) then
        IF (TG1.GT.TDIURN(I,J,8)) TDIURN(I,J,8) = TG1
        IF (MODDSF.eq.0) AIJ(I,J,IJ_TSLI)=AIJ(I,J,IJ_TSLI)+(TS-TF)
        AIJ(I,J,IJ_SHDTLI)=AIJ(I,J,IJ_SHDTLI)+SHDT
        AIJ(I,J,IJ_EVHDT)=AIJ(I,J,IJ_EVHDT)+EVHDT
        AIJ(I,J,IJ_TRHDT)=AIJ(I,J,IJ_TRHDT)+TRHDT
        AIJ(I,J,IJ_F0LI)=AIJ(I,J,IJ_F0LI)+F0DT
        AIJ(I,J,IJ_EVAPLI)=AIJ(I,J,IJ_EVAPLI)+EVAP
C****
!!!      END SELECT
      endif
C****
      END IF
      END DO   ! end of itype loop
      IM1=I
      END DO   ! end of I loop

      END DO   ! end of J loop
!$OMP  END PARALLEL DO

      Call CHECKSUM(GRID, uosurf, 952, "SURFACE.f")
      Call CHECKSUM(GRID, vosurf, 953, "SURFACE.f")
      Call CHECKSUM(GRID, uisurf, 954, "SURFACE.f")
      Call CHECKSUM(GRID, visurf, 955, "SURFACE.f")

      DO J=J_0,J_1
      DO I=1,IMAXJ(J)
        JR=JREG(I,J)
        DO K=1,3
          AREG(JR,J_TRHDT)=AREG(JR,J_TRHDT)+AREGIJ(1,K,I,J)
          AREG(JR,J_SHDT )=AREG(JR,J_SHDT )+AREGIJ(2,K,I,J)
          AREG(JR,J_EVHDT)=AREG(JR,J_EVHDT)+AREGIJ(3,K,I,J)
          AREG(JR,J_EVAP )=AREG(JR,J_EVAP )+AREGIJ(4,K,I,J)
          IF(MODDSF.EQ.0) THEN
            AREG(JR,J_TSRF )=AREG(JR,J_TSRF )+AREGIJ(5,K,I,J)
            AREG(JR,J_TG1)=AREG(JR,J_TG1)+AREGIJ(6,K,I,J)
            AREG(JR,J_TG2)=AREG(JR,J_TG2)+AREGIJ(7,K,I,J)
          END IF
        END DO
      END DO
      END DO
C****
C**** EARTH
C****
      CALL EARTH(NS,MODDSF,MODDD)
C****
C**** UPDATE FIRST LAYER QUANTITIES
C****
!$OMP  PARALLEL DO PRIVATE (I,J,FTEVAP,FQEVAP,P1K)
!$OMP*          SCHEDULE(DYNAMIC,2)
      DO J=J_0,J_1
      DO I=1,IMAXJ(J)
        FTEVAP=0
        IF (DTH1(I,J)*T(I,J,1).lt.0) FTEVAP=-DTH1(I,J)/T(I,J,1)
        FQEVAP=0
        IF (DQ1(I,J).lt.0.and.Q(I,J,1).gt.0) FQEVAP=-DQ1(I,J)/Q(I,J,1)
! Z-moments should be set from PBL
        TMOM(:,I,J,1) = TMOM(:,I,J,1)*(1.-FTEVAP)
        QMOM(:,I,J,1) = QMOM(:,I,J,1)*(1.-FQEVAP)
        IF ( Q(I,J,1)+DQ1(I,J) .LT. qmin ) THEN
          WRITE(99,*)
     &         ITime,'I,J:',I,J,' Q1:',Q(I,J,1)+DQ1(I,J),'->',qmin
          dq1(i,j)=qmin-q(i,j,1)
          QMOM(:,I,J,1)=0.
        ENDIF
c****   retrieve fluxes
        P1K=PK(1,I,J)
        tflux1(i,j)=dth1(i,j)*(-AM(1,I,J)*P1K)/(dtsurf)
        qflux1(i,j)=dq1(i,j)*(-AM(1,I,J))/(dtsurf)
C**** Diurnal cycle of temperature diagnostics
        tdiurn(i,j,5)=tdiurn(i,j,5)+(tsavg(i,j)-tf)
        if(tsavg(i,j).gt.tdiurn(i,j,6)) tdiurn(i,j,6)=tsavg(i,j)
        if(tsavg(i,j).lt.tdiurn(i,j,9)) tdiurn(i,j,9)=tsavg(i,j)
      END DO
      END DO
!$OMP  END PARALLEL DO

c****
c**** apply surface fluxes to the first layer of the atmosphere
c****  (replaced with dummy subroutine when ATURB is used)
c****
      call apply_fluxes_to_atm(dtsurf)

C**** Call dry convection or aturb depending on rundeck
      CALL ATM_DIFFUS(1,1,dtsurf)
C****
C**** ACCUMULATE SOME ADDITIONAL BOUNDARY LAYER DIAGNOSTICS
C****
      IF(MODDD.EQ.0) THEN
        DO KR=1,NDIUPT
C**** CHECK IF DRY CONV HAS HAPPENED FOR THIS DIAGNOSTIC
          IF(DCLEV(IJDD(1,KR),IJDD(2,KR)).GT.1.) THEN
            ADIURN(IH,IDD_DCF,KR)=ADIURN(IH,IDD_DCF,KR)+1.
            ADIURN(IH,IDD_LDC,KR)=ADIURN(IH,IDD_LDC,KR)
     *           +DCLEV(IJDD(1,KR),IJDD(2,KR))
            HDIURN(IHM,IDD_DCF,KR)=HDIURN(IHM,IDD_DCF,KR)+1.
            HDIURN(IHM,IDD_LDC,KR)=HDIURN(IHM,IDD_LDC,KR)
     *           +DCLEV(IJDD(1,KR),IJDD(2,KR))
          END IF
        END DO
      END IF
C****
      END DO   ! end of surface time step



      RETURN
C****
      END SUBROUTINE SURFCE