C**********************************************************                  BLSURF.2     
C             SUBROUTINE BLSURF                                              BLSURF.3     
C**********************************************************                  BLSURF.4     
      subroutine BLSURF                                                      BLSURF.5     
C-----------------------------------------------------------------------     BLSURF.6     
C     Perform time-splitting integration of the boundary layer and           BLSURF.7     
C     Surface scheme together for a whole timestep                           BLSURF.8     
C     In the case of LSL or LOC being false, should forget the               BLSURF.9     
C     timesplitting and just do the BL stuff (leave the surface alone)       BLSURF.10    
C                                                                            BLSURF.11    
C-----------------------------------------------------------------------     BLSURF.12    
C                                                                            PARAM1.2     
C     Determines model resolution                                            PARAM1.3     
C                                                                            PARAM1.4     
      PARAMETER(NN=21,MM=21,NHEM=2,NL=5,MOCT=1,MG=64,JG=16,NWJ2=121          NETCDF.35    
     P         ,NCRAY=64,JGL=JG,NTRAC=1,NLEVRF=1)                            NETCDF.36    
                                                                             PARAM1.7     
C                                                                            PARAM1.8     
C                                                                            PARAM2.2     
C     Sets basic constants, especially those needed for array dimensions     PARAM2.3     
C                                                                            PARAM2.4     
      PARAMETER(MH=2,PI=3.14159265359,PI2=2.0*PI                             PARAM2.5     
     +,NNP=NN+1,MGPP=MG+2,JGP=JG+1,JGG=JG*NHEM,JGGP=JGG+1,MJP=NWJ2+NWJ2      PARAM2.6     
     +,NLM=NL-1,NLP=NL+1,NLPP=NL+2,NLA=NL+3,NLB=NL+4,NL2=NL*NL               PARAM2.7     
     +,IDA=(MG+MG+MG)/2+1,IDB=NWJ2*NL,IDC=IDB+IDB,IDD=MGPP*NL                PARAM2.8     
     +,IDE=NL2*NN,IDF=NCRAY*(MG+1),IDG=JG*NL,IDH=JG*MG                       PARAM2.9     
     +,IDI=NNP/2,IDJ=IDI*IDI,IDK=NL*IDI,IDL=MGPP/2,IDM=NNP/2,IDN=IDM*NL      PARAM2.10    
     +,NWW=1+(MM-1)/MOCT)                                                    PARAM2.11    
      PARAMETER(IGA=NWJ2*NHEM,IGB=IDB*NHEM,IGC=MGPP*NHEM,IGD=IDD*NHEM        PARAM2.12    
     +,IGG=IDG*NHEM,IGL=IDL*NHEM,IGM=IDM*NHEM,IGN=IDN*NHEM                   PARAM2.13    
     +,IGO=IGA+IGA,IGP=IGB+IGB,NFTWG=(5+NTRAC)*NL+3                          PARAM2.14    
     +,NFTGW=(6+3*NTRAC)*NL+2,NFTGD=(3+NTRAC)*NL,NLTR=NL*NTRAC)              PARAM2.15    
C     Number of 2D (surface) output fields. This value is                    PARAM2.16    
C     Doubled due to averaged and instantaneous fields.                      PARAM2.17    
      PARAMETER (N2DFLD=21,NGRPAD=N2DFLD*2*IGC)                              PARAM2.18    
C                                                                            PARAM2.19    
C                                                                            BLANK.2     
C     Basic planetary parameters for run plus information about              BLANK.3     
C     vertical grid structure                                                BLANK.4     
C                                                                            BLANK.5     
C     Note that RD and GASCON are identical and CPD is set from RD,AKAP.     BLANK.6     
      COMMON        SQ(NNP),RSQ(NNP),SIGMAH(NLM),SIGMA(NL)                   BLANK.7     
     +              ,T01S2(NLM),T0(NL),ALPHA(NL),DSIGMA(NL),RDSIG(NL)        BLANK.8     
     +              ,TKP(NL),C(NL2),SQH(NNP)                                 BLANK.9     
     +              ,MF,MFP,JZF,NF,NFP                                       BLANK.10    
     +              ,AKAP,GA,GASCON,RADEA,WW,P0,PFAC,EZ,AIOCT                BLANK.11    
     +              ,RD,RV,CPD,CLATNT                                        BLANK.12    
     +              ,LRSTRT,LSHORT,LTVEC,LSTRETCH                            BLANK.13    
     +              ,LFLUX                                                   BLANK.14    
     +              ,LBALAN,LRESTIJ                                          BLANK.15    
     +              ,LCLIM, LPERPET, L22L,LOROG ,LCSFCT                      BLANK.16    
     +              ,LNOISE                                                  BLANK.17    
      COMPLEX EZ,AIOCT                                                       BLANK.18    
      LOGICAL LRSTRT,LSHORT,LTVEC,LSTRETCH,LBALAN,LRESTIJ                    BLANK.19    
     +       ,LFLUX,LNOISE                                                   BLANK.20    
     +       ,LCLIM, LPERPET, L22L,LOROG,LCSFCT                              BLANK.21    
C                                                                            BLANK.22    
C                                                                            LEGAU.2     
C     Legendre polynomials and information about gaussian latitudes          LEGAU.3     
C                                                                            LEGAU.4     
      COMMON/LEGAU/ ALPJ(MJP),DALPJ(MJP)                                     LEGAU.5     
     +              ,ALP(NWJ2,2,JGL),DALP(NWJ2,2,JGL)                        LEGAU.6     
     +              ,RLP(NWJ2,2,JGL),RDLP(NWJ2,2,JGL)                        LEGAU.7     
     +              ,SI(JGG),CS(JGG),SISQ(JGG),CSSQ(JGG),SECSQ(JGG)          LEGAU.8     
     +              ,ALAT(JGG),GWT(JGG),AW(JGG),JH,JL,JINC                   LEGAU.9     
C                                                                            LEGAU.10    
C                                                                            GRIDPP.2     
C     Array ordering in GRIDP must correspond to that in SPECTR.             GRIDPP.3     
C     Real arrays: multi-level arrays are 2-dimensional.                     GRIDPP.4     
C     the variables have been renamed to coincide with                       GRIDPP.5     
C     variable names in bgcm5 DGRMLT                                         GRIDPP.6     
C                                                                            GRIDPP.7     
C swapped round VLNG and UNLG to check                                       GRIDPP.8     
      COMMON/GRIDP/ CHIG(IGC,NL),SFG(IGC,NL),UG(IGC,NL),VG(IGC,NL)           GRIDPP.9     
     :              ,TTVD(IGC,NL),QTVD(IGC,NL),TG(IGC,NL)                    GRIDPP.10    
     :              ,TRAG(IGC,NL,NTRAC)                                      GRIDPP.11    
     :              ,PLG(IGC),TYBL(IGC),TXBL(IGC)                            GRIDPP.12    
     :              ,SPG(IGC),VPG(IGC),TTRD(IGC,NL)                          GRIDPP.13    
     :              ,TNLG(IGC,NL),TRANLG(IGC,NL,NTRAC),UNLG(IGC,NL)          GRIDPP.14    
     :              ,VNLG(IGC,NL),TTLR(IGC,NL),UTRAG(IGC,NL,NTRAC)           GRIDPP.15    
     :              ,TTCR(IGC,NL),VTRAG(IGC,NL,NTRAC)                        GRIDPP.16    
     :              ,UTVD(IGC,NL),VTVD(IGC,NL)                               GRIDPP.17    
     :         ,ASSBL(IGC),ASHBL(IGC),ASLBL(IGC),ARRCR(IGC),ARRLR(IGC)       GRIDPP.18    
     :         ,arflux(igc,6),asfld(igc,6),acld(igc,4)                       GRIDPP.19    
     :         ,SSBL(IGC),SHBL(IGC),SLBL(IGC),RRCR(IGC),RRLR(IGC)            GRIDPP.20    
     :         ,rflux(igc,6),sfld(igc,6),cld(igc,4)                          GRIDPP.21    
C                                                                            GRIDPP.22    
C                                                                            BATS.2     
C     Constant arrays and variables associated with time and vertical        BATS.3     
C     differencing. Also counters.                                           BATS.4     
C                                                                            BATS.5     
      COMMON/BATS/  BM1(IDE),AK(NNP),AQ(NL2),G(NL2),TAU(NL2)                 BATS.6     
     +              ,KOUNT,KITS,KSTART,KTOTAL,KRUN,BEGDAY,ITSPD              BATS.7     
     +              ,DELT,DELT2,CV,CG,CT,CQ,PNU,PNU2,PNU21                   BATS.8     
     +              ,NTRACO,CTRA(NTRAC),KOLOUR(NTRAC),RGG(NL2)               BATS.9     
     +              ,BEGDOY,DOY                                              BATS.10    
C                                                                            BATS.11    
      COMMON/PHYS/  LBL,LVD,LCR,LLR,LRD,LCUBM,LCBADJ                         PHYSCA.2     
     :              ,TSTAR(IGC,JG),QSTAR(IGC,JG),FRAD(JG,NHEM)               PHYSCA.3     
     :              ,TSTARO(IGC,JG),TDEEPO(IGC,JG),smstar(igc,jg)            PHYSCA.4     
     :              ,tdeep(igc,jg),hsnow(igc,jg),sqstar(igc,jg)              PHYSCA.5     
     :              ,SALB(IGC,JG),SBAL(IGC,JG),BLCD(IGC)                     PHYSCA.6     
     :              ,SVEGE(IGC,JG),CD,DRAG,BLVAD,BLA,BLRH,BLVB(IGC)          PHYSCA.7     
     :              ,AKVV,AKTV,AKQV,ESCONA,ESCONB,EPSIQ,CTQ,CCC              PHYSCA.8     
     : ,ctqi,sdsn,shcs,shcsp,shcsn,skse,sksn,slhf,sd1,sd2,sdw                PHYSCA.9     
     :        ,ssmc,sdsnd,LSL,sasnow,saice,shsstar,shsmax                    PHYSCA.10    
     :     ,LOC,SHCO,SHCI,LNOICE,LOLDBL,LCOND,LNNSK,ITSLL,ITSLO              PHYSCA.11    
     :              ,CCR,RCON,DTBUOY,TSLA,TSLB,TSLC,TSLD,CUT1,CUT2           PHYSCA.12    
     :              ,NLCR,NCUTOP,CURHM,AKTC,AKQC,CUBMT,CBADJT,CBADJP         PHYSCA.13    
     :              ,SKAP(NL),SK(NLM),FWS(NL),CLR(NL),FB(NLM)                PHYSCA.14    
     :              ,TTDC(NL),QTDC(NL),TTMC(NL),QTMC(NL),TC(NL),QC(NL)       PHYSCA.15    
     :              ,CTCR(NL,NHEM),CTLR(NL,NHEM),NAVRD,NAVWT,DELT2C          PHYSCA.16    
      LOGICAL LBL,LVD,LCR,LLR,LRD,LCUBM,LCBADJ,LSL,LOC                       PHYSCA.17    
     :       ,LNOICE,LOLDBL,LCOND,LNNSK                                      PHYSCA.18    
C                                                                            PHYSCA.19    
       COMMON/CPIERS/ICFLAG(IGC,5,2),CFRAC(IGC,5),PNET(IGC,JG)               CPIERS.2     
     :               ,SNET(IGC,JG),RRFLUX(IGC,JG,6)                          CPIERS.3     
     :               ,TTSW(IGC,NL),TTLW(IGC,NL)                              CPIERS.4     
       COMMON/GSG/GSG(IGC,JG)                                                OROG_G.2     
C                                                                            XTRAP.2     
C     Setup moisture variables by equivilencing them to                      XTRAP.3     
C     Tracer No. 1                                                           XTRAP.4     
C                                                                            XTRAP.5     
      REAL QG(IGC,NL),QNLG(IGC,NL),QTLR(IGC,NL),QTCR(IGC,NL)                 XTRAP.6     
      EQUIVALENCE (QG(1,1),TRAG(1,1,1)) , (QNLG(1,1),TRANLG(1,1,1))          XTRAP.7     
     *            ,(QTLR(1,1),UTRAG(1,1,1)),(QTCR(1,1),VTRAG(1,1,1))         XTRAP.8     
C                                                                            XTRAP.9     
       REAL OCFLUX(130,16,12)  ! OCEAN HEAT FLUX MONTHLY-AVE                 BLSURF.23    
C-----------------------------------------------------------------------     BLSURF.24    
C     set up values for reference pressure and various temperatures          BLSURF.25    
C     in non-dimensionalised form                                            BLSURF.26    
C-----------------------------------------------------------------------     BLSURF.27    
      parameter (tzc=273.15)                                                 BLSURF.28    
      logical LSLH,LOCH                                                      BLSURF.29    
      DATA IFIRST/1/                                                         BLSURF.30    
C                                                                            BLSURF.31    
C Lookup table for vegetation to Z0: from Unified Model Doc. no. 70          BLSURF.32    
C Snow free lookup table                            MMJ April 2000           BLSURF.33    
C                                                                            BLSURF.34    
       DIMENSION SZLOOK(24)                                                  BLSURF.35    
      save apfrac,mpth1,mpth2                                                NETCDF.52    
       DATA SZLOOK /0.001,1.0E-4,3.0E-4,1.0,1.2,1.0,                         BLSURF.36    
     $ 1.0,1.2,1.0,0.4,0.4,0.4,                                              BLSURF.37    
     $ 0.12,0.12,0.12,0.12,0.12,0.12,                                        BLSURF.38    
     $ 0.12,0.12,1.5,0.12,0.12,3.0E-3/                                       BLSURF.39    
C                                                                            BLSURF.40    
C   THIS BIT READS IN RUN 1.0 (10 year ocean flux)                           BLSURF.41    
C 30-3-98  PMF                                                               BLSURF.42    
      IF (IFIRST.EQ.1.AND.LOC) THEN                                          BLSURF.43    
        IFIRST=0                                                             BLSURF.44    
C OCEAN FLUX is Wm-2 +ve Downwards                                           BLSURF.45    
C FILE PRODUCED BY AVERAGING 10 YEAR UTF SURFACE BALANCE TERM                BLSURF.46    
C OF RUN 1.0                                                                 BLSURF.47    
C                                                                            BLSURF.48    
        READ(41,*)                                                           BLSURF.49    
        READ(41,*) OCFLUX                                                    BLSURF.50    
        CLOSE(41)                                                            BLSURF.51    
        PRINT *,' READ IN OCEAN FLUXES'                                      BLSURF.52    
      ENDIF                                                                  BLSURF.53    
c                                                                            BLSURF.54    
C  SORT OUT TIMINGS                                                          BLSURF.55    
c                                                                            BLSURF.56    
      IF (JH.EQ.1.AND.LOC) THEN                                              BLSURF.57    
         CALL CALNDR(DOY,Mpth1,ApFRAC)                                       BLSURF.58    
         MPth2=mpth1+1                                                       BLSURF.59    
         if (mpth2.eq.13) mpth2=1                                            BLSURF.60    
      ENDIF                                                                  BLSURF.61    
      tp1=(tzc+1.)/ct                                                        BLSURF.62    
      tz=tzc/ct                                                              BLSURF.63    
      tm2=(tzc-2.)/ct                                                        BLSURF.64    
      tm3=(tzc-3.)/ct                                                        BLSURF.65    
                                                                             BLSURF.66    
C     Set up the number of timesplitting levels for land and ocean,          BLSURF.67    
C     also the timesteps.                                                    BLSURF.68    
                                                                             BLSURF.69    
      DELTL=DELT2C/ITSLL                                                     BLSURF.70    
      DELTO=DELT2C/ITSLO                                                     BLSURF.71    
                                                                             BLSURF.72    
                                                                             BLSURF.73    
C     Set up some more constants for the soil scheme                         BLSURF.74    
                                                                             BLSURF.75    
      ssmc23=ssmc*(2./3.)                                                    BLSURF.76    
      ssmc13=ssmc*(1./3.)                                                    BLSURF.77    
      srsum=sd1/skse+sd2/skse                                                BLSURF.78    
                                                                             BLSURF.79    
C     And some for the BL scheme                                             BLSURF.80    
                                                                             BLSURF.81    
      RCSJ=SECSQ(JH)                                                         BLSURF.82    
      SQRC=SQRT(RCSJ)                                                        BLSURF.83    
      RSIGF=2.*(1-SIGMA(NL))/(1+SIGMA(NL))                                   BLSURF.84    
C     5cm roughness length over land - note includes factor of g             BLSURF.85    
      RGZ0L=2.2726E-6                                                        BLSURF.86    
C     1mm over ocean                                                         BLSURF.87    
      RGZ0O=4.5453E-8                                                        BLSURF.88    
                                                                             BLSURF.89    
                                                                             BLSURF.90    
C     Now begin looping over gridpoints                                      BLSURF.91    
                                                                             BLSURF.92    
      do ihem=1,nhem                                                         BLSURF.93    
         do i=1,mg                                                           BLSURF.94    
            j=(ihem-1)*mgpp+i                                                BLSURF.95    
                                                                             BLSURF.96    
C                                                                            BLSURF.97    
C Get land roughness lengths from vegetation index                           BLSURF.98    
C                                                                            BLSURF.99    
            RGZ0L=SZLOOK(NINT(SVEGE(J,JH)))/2.2E4                            BLSURF.100   
C                                                                            BLSURF.101   
C     Now do a single gridpoint                                              BLSURF.102   
                                                                             BLSURF.103   
            LSLH=(LSL.and.(gsg(j,jh).gt.0.))                                 BLSURF.104   
            LOCH=(LOC.and.(gsg(j,jh).le.0.))                                 BLSURF.105   
                                                                             BLSURF.106   
            if (gsg(j,jh).gt.0.) then                                        BLSURF.107   
               MAXITER=ITSLL                                                 BLSURF.108   
               DELTU=DELTL                                                   BLSURF.109   
               RGZ0=RGZ0L                                                    BLSURF.110   
               if (LSL) then                                                 BLSURF.111   
                  TSCUR=TSTARO(J,JH)                                         BLSURF.112   
                  TDCUR=TDEEPO(J,JH)                                         BLSURF.113   
               else                                                          BLSURF.114   
                  TSCUR=TSTAR(J,JH)                                          BLSURF.115   
                  QCUR=QSTAR(J,JH)                                           BLSURF.116   
                  SQCUR=SQSTAR(J,JH)                                         BLSURF.117   
               endif                                                         BLSURF.118   
            else                                                             BLSURF.119   
               MAXITER=ITSLO                                                 BLSURF.120   
               DELTU=DELTO                                                   BLSURF.121   
               RGZ0=RGZ0O                                                    BLSURF.122   
               if (LOC) then                                                 BLSURF.123   
                  TSCUR=TSTARO(J,JH)                                         BLSURF.124   
               else                                                          BLSURF.125   
                  TSCUR=TSTAR(J,JH)                                          BLSURF.126   
                  QCUR=QSTAR(J,JH)                                           BLSURF.127   
                  SQCUR=SQSTAR(J,JH)                                         BLSURF.128   
               endif                                                         BLSURF.129   
            endif                                                            BLSURF.130   
                                                                             BLSURF.131   
            QNLGS=0.                                                         BLSURF.132   
            TNLGS=0.                                                         BLSURF.133   
            UNLGS=0.                                                         BLSURF.134   
            VNLGS=0.                                                         BLSURF.135   
            SHBLS=0.                                                         BLSURF.136   
            SLBLS=0.                                                         BLSURF.137   
            BLCDS=0.                                                         BLSURF.138   
            BLVBS=0.                                                         BLSURF.139   
            TXBLS=0.                                                         BLSURF.140   
            TYBLS=0.                                                         BLSURF.141   
                                                                             BLSURF.142   
            do iter=1,MAXITER                                                BLSURF.143   
                                                                             BLSURF.144   
C     if we have interactive soil, calculate surface humidities              BLSURF.145   
                                                                             BLSURF.146   
               if (LSLH) then                                                BLSURF.147   
                  if (smstar(j,jh).gt.ssmc) then                             BLSURF.148   
                     srh=1.                                                  BLSURF.149   
                     smstar(j,jh)=ssmc                                       BLSURF.150   
                  else if ((smstar(j,jh).gt.ssmc23)                          BLSURF.151   
     $                    .or.(hsnow(j,jh).gt.0.)) then                      BLSURF.152   
                     srh=1.                                                  BLSURF.153   
                  else if (smstar(j,jh).gt.ssmc13) then                      BLSURF.154   
                     srh=(smstar(j,jh)/ssmc13)-1.                            BLSURF.155   
                  else                                                       BLSURF.156   
                     srh=0.                                                  BLSURF.157   
                  endif                                                      BLSURF.158   
                                                                             BLSURF.159   
                  sqcur=pqsat(tscur)/plg(j)                                  BLSURF.160   
                                                                             BLSURF.161   
                  qcur=srh*sqcur                                             BLSURF.162   
                                                                             BLSURF.163   
               endif                                                         BLSURF.164   
                                                                             BLSURF.165   
C     Ditto for interactive ocean                                            BLSURF.166   
                                                                             BLSURF.167   
               if (LOCH) then                                                BLSURF.168   
                  sqcur=pqsat(tscur)/plg(j)                                  BLSURF.169   
                  qcur=sqcur                                                 BLSURF.170   
               endif                                                         BLSURF.171   
                                                                             BLSURF.172   
               THNL=TG(J,NL)/SKAP(NL)                                        BLSURF.173   
               THBAR=0.5*(THNL+TSCUR)                                        BLSURF.174   
               RZZ0=RSIGF*THBAR/RGZ0                                         BLSURF.175   
               RLZZ02=(LOG(RZZ0))**2                                         BLSURF.176   
C     CD=(K/LN(Z/Z0))**2, K=0.41                                             BLSURF.177   
               CD=0.1681/RLZZ02                                              BLSURF.178   
               blcds=blcds+cd/maxiter                                        BLSURF.179   
               CTAU=CD/TSCUR                                                 BLSURF.180   
               CSH=CTAU/AKAP                                                 BLSURF.181   
               DRAGJ=CD*DRAG/TSCUR                                           BLSURF.182   
C     CALCULATE WIND SPEED, INCLUDING MINIMUM                                BLSURF.183   
               VM=SQRT(RCSJ*(UG(J,NL)*UG(J,NL)+VG(J,NL)*VG(J,NL)))           BLSURF.184   
               VMP=VM+BLVAD                                                  BLSURF.185   
C     CALCULATE SURFACE STRESS AND LOWEST LEVEL FRICTION                     BLSURF.186   
               CVM=DRAGJ*VMP                                                 BLSURF.187   
               UNLG(J,NL)=-CVM*UG(J,NL)                                      BLSURF.188   
               UNLGS=UNLGS+UNLG(J,NL)/MAXITER                                BLSURF.189   
               VNLG(J,NL)=-CVM*VG(J,NL)                                      BLSURF.190   
               VNLGS=VNLGS+VNLG(J,NL)/MAXITER                                BLSURF.191   
               CVM=CTAU*PLG(J)*VMP*SQRC                                      BLSURF.192   
               TXBL(J)=CVM*UG(J,NL)                                          BLSURF.193   
               TXBLS=TXBLS+TXBL(J)/MAXITER                                   BLSURF.194   
               TYBL(J)=CVM*VG(J,NL)                                          BLSURF.195   
               TYBLS=TYBLS+TYBL(J)/MAXITER                                   BLSURF.196   
C     CALCULATE SENSIBLE AND LATENT HEAT FLUXES                              BLSURF.197   
               QCUR=MIN(SQCUR,MAX(QCUR,QG(J,NL)))                            BLSURF.198   
               DTH=TSCUR-THNL                                                BLSURF.199   
               DQ=BLRH*(QCUR-QG(J,NL))                                       BLSURF.200   
               BLVB(J)=0.                                                    BLSURF.201   
               IF (DTH.GT.0.) THEN                                           BLSURF.202   
                  BLVB(J)=5.95*RLZZ02*SQRT(RGZ0*DTH/THBAR)                   BLSURF.203   
C                  BLVB(J)=5.95*RLZZ02*SQRT(RGZ0L*DTH/THBAR)                 BLSURF.204   
C                  if (gsg(j,jh).le.0.) then                                 BLSURF.205   
C                     BLVB(J)=500./464.*SQRT(DTH/THBAR)                      BLSURF.206   
C                  endif                                                     BLSURF.207   
C     RESTRICT ENHANCEMENT TO A FACTOR OF 4 OVER NEUTRAL VALUES              BLSURF.208   
C                  IF (BLVB(J)/VM.GT.4.) BLVB(J)=VM*4.                       BLSURF.209   
               ENDIF                                                         BLSURF.210   
C                                                                            BLSURF.211   
C Cd_h is multiplied by 0.2: I think this underestimates                     BLSURF.212   
C Cd_h in unstable conditions SO                                             BLSURF.213   
C MAKE Cd_h = Cd, not 0.2*Cd in unstable conditions                          BLSURF.214   
C                                                                            BLSURF.215   
C               VMP=VM+BLVB(J)                                               BLSURF.216   
               VMP=VM+5.0*BLVB(J)                                            BLSURF.217   
               CVM=0.2*DRAGJ*VMP                                             BLSURF.218   
               TNLG(J,NL)=CVM*DTH                                            BLSURF.219   
               TNLGS=TNLGS+TNLG(J,NL)/MAXITER                                BLSURF.220   
               QNLG(J,NL)=CVM*DQ                                             BLSURF.221   
               QNLGS=QNLGS+QNLG(J,NL)/MAXITER                                BLSURF.222   
               CVM=0.2*CSH*PLG(J)*VMP                                        BLSURF.223   
               SHBL(J)=CVM*DTH                                               BLSURF.224   
               SHBLS=SHBLS+SHBL(J)/MAXITER                                   BLSURF.225   
C     CHOOSE APPROPRIATE LATENT HEAT                                         BLSURF.226   
               IF (TSCUR.GT.TZ) THEN                                         BLSURF.227   
                  SLBL(J)=CVM*CTQ*DQ                                         BLSURF.228   
               ELSE                                                          BLSURF.229   
                  SLBL(J)=CVM*CTQI*DQ                                        BLSURF.230   
               ENDIF                                                         BLSURF.231   
               SLBLS=SLBLS+SLBL(J)/MAXITER                                   BLSURF.232   
               UG(J,NL)=UG(J,NL)+DELTU*UNLG(J,NL)                            BLSURF.233   
               VG(J,NL)=VG(J,NL)+DELTU*VNLG(J,NL)                            BLSURF.234   
               TG(J,NL)=TG(J,NL)+DELTU*TNLG(J,NL)                            BLSURF.235   
               QG(J,NL)=QG(J,NL)+DELTU*QNLG(J,NL)                            BLSURF.236   
                                                                             BLSURF.237   
               if (LSLH) then                                                BLSURF.238   
                                                                             BLSURF.239   
C                  if ((j.eq.18).and.(jh.eq.10))then                         BLSURF.240   
C                     write (*,*) 'Initial TS=',tscur*                       BLSURF.241   
C     $                    752.-273.15                                       BLSURF.242   
C                  endif                                                     BLSURF.243   
                                                                             BLSURF.244   
C     CHECK SIGNS of shbl,slbl,snet                                          BLSURF.245   
                                                                             BLSURF.246   
                  SFC=-shbl(j)-slbl(j)+snet(j,jh)/(CV*P0)                    BLSURF.247   
                                                                             BLSURF.248   
                  rkappa=2./(srsum+hsnow(j,jh)/sksn)                         BLSURF.249   
                  pc1=shcs                                                   BLSURF.250   
                  pc2=shcs                                                   BLSURF.251   
                  if ((tscur.gt.tm3).and.(tscur.lt.tp1))                     BLSURF.252   
     $                 pc1=pc1+shcsp                                         BLSURF.253   
                  if ((tdcur.gt.tm3).and.(tdcur.lt.tp1))                     BLSURF.254   
     $                 pc1=pc1+shcsp                                         BLSURF.255   
                  hc1=pc1*sd1+min(hsnow(j,jh),shsmax)*shcsn                  BLSURF.256   
                  hc2=pc2*sd2                                                BLSURF.257   
                                                                             BLSURF.258   
                  trans=rkappa*(tscur-tdcur)                                 BLSURF.259   
                  rct1=(sfc-trans)/hc1                                       BLSURF.260   
                  rct2=trans/hc2                                             BLSURF.261   
                                                                             BLSURF.262   
                  tscur=tscur+rct1*deltu                                     BLSURF.263   
                  tdcur=tdcur+rct2*deltu                                     BLSURF.264   
                                                                             BLSURF.265   
                  if (tscur.gt.tz) then                                      BLSURF.266   
                     if (hsnow(j,jh).gt.0.) then                             BLSURF.267   
C     The factor of 2 is *REALLY* meant to be here, since we                 BLSURF.268   
C     are dealing with a forward time single level variable (hsnow)          BLSURF.269   
C     within an (effectively) centre time timestep - thus each               BLSURF.270   
C     segment of time is covered twice in this routine...                    BLSURF.271   
                        hsnow(j,jh)=hsnow(j,jh)                              BLSURF.272   
     $                       -(tscur-tz)*hc1/(slhf*sdsn)/2.                  BLSURF.273   
                        tscur=tz                                             BLSURF.274   
                     endif                                                   BLSURF.275   
                  endif                                                      BLSURF.276   
                                                                             BLSURF.277   
               endif                                                         BLSURF.278   
                                                                             BLSURF.279   
               if (LOCH) then                                                BLSURF.280   
c                                                                            BLSURF.281   
C INTEROP between months                                                     BLSURF.282   
          OCF=(1.0-ApFRAC)*OCFLUX(J,JH,MpTH1)+                               BLSURF.283   
     $        ApFRAC*OCFLUX(J,JH,MpTH2)                                      BLSURF.284   
c                                                                            BLSURF.285   
                  if (tscur.lt.tm2) then !ice                                BLSURF.286   
              if (ocf.gt.0) Then  ! flux corr out of ice                     BLSURF.287   
                    ocf=0   !ice can't cool to ocean                         BLSURF.288   
              endif                                                          BLSURF.289   
                     hc=shci                                                 BLSURF.290   
                  else                                                       BLSURF.291   
                     hc=shco                                                 BLSURF.292   
                  endif                                                      BLSURF.293   
c                                                                            BLSURF.294   
C add OCF term to ocean heat flux                                            BLSURF.295   
                SFC=-shbl(j)-slbl(j)+(snet(j,jh)-OCF)/(CV*P0)                BLSURF.296   
c               SFC=-shbl(j)-slbl(j)+(snet(j,jh))/(CV*P0)                    BLSURF.297   
c                                                                            BLSURF.298   
                  dts=sfc/hc*DELTU                                           BLSURF.299   
c                                                                            BLSURF.300   
                  if ((tscur.lt.tm2).and.(tscur+dts                          BLSURF.301   
     $                 .ge.tm2)) then                                        BLSURF.302   
                     tscur=tm2+(tscur+dts-tm2)*shci/shco                     BLSURF.303   
                  else if ((tscur.ge.tm2).and.(                              BLSURF.304   
     $                    tscur+dts.lt.tm2)) then                            BLSURF.305   
                     tscur=tm2+(tscur+dts-tm2)*shco/shci                     BLSURF.306   
                  else                                                       BLSURF.307   
                     tscur=tscur+dts                                         BLSURF.308   
                  endif                                                      BLSURF.309   
C                                                                            BLSURF.310   
                  if (tscur.ge.tm2) then                                     BLSURF.311   
                     salb(j,jh)=sbal(j,jh)                                   BLSURF.312   
                  else                                                       BLSURF.313   
                     salb(j,jh)=saice                                        BLSURF.314   
                  endif                                                      BLSURF.315   
                                                                             BLSURF.316   
               endif                                                         BLSURF.317   
            enddo                                                            BLSURF.318   
                                                                             BLSURF.319   
C     End of iteration process, done full double timestep.                   BLSURF.320   
                                                                             BLSURF.321   
            if (LSLH) then                                                   BLSURF.322   
               tdeepo(j,jh)=tdeep(j,jh)*(1.-2.*pnu)                          BLSURF.323   
     $              +pnu*(tdeepo(j,jh)+tdcur)                                BLSURF.324   
               tdeep(j,jh)=tdcur                                             BLSURF.325   
               tstaro(j,jh)=tstar(j,jh)*(1.-2.*pnu)                          BLSURF.326   
     $              +pnu*(tstaro(j,jh)+tscur)                                BLSURF.327   
               tstar(j,jh)=tscur                                             BLSURF.328   
               sqstar(j,jh)=sqcur                                            BLSURF.329   
               qstar(j,jh)=qcur                                              BLSURF.330   
            endif                                                            BLSURF.331   
                                                                             BLSURF.332   
            if (LOCH) then                                                   BLSURF.333   
               tstaro(j,jh)=tstar(j,jh)*(1.-2.*pnu)                          BLSURF.334   
     $              +pnu*(tstaro(j,jh)+tscur)                                BLSURF.335   
               tstar(j,jh)=tscur                                             BLSURF.336   
               sqstar(j,jh)=sqcur                                            BLSURF.337   
               qstar(j,jh)=qcur                                              BLSURF.338   
            endif                                                            BLSURF.339   
                                                                             BLSURF.340   
            QNLG(J,NL)=QNLGS                                                 BLSURF.341   
            TNLG(J,NL)=TNLGS                                                 BLSURF.342   
            UNLG(J,NL)=UNLGS                                                 BLSURF.343   
            VNLG(J,NL)=VNLGS                                                 BLSURF.344   
            SHBL(J)=SHBLS                                                    BLSURF.345   
            SLBL(J)=SLBLS                                                    BLSURF.346   
            BLCD(J)=BLCDS                                                    BLSURF.347   
            BLVB(J)=BLVBS                                                    BLSURF.348   
            TXBL(J)=TXBLS                                                    BLSURF.349   
            TYBL(J)=TYBLS                                                    BLSURF.350   
                                                                             BLSURF.351   
         enddo                                                               BLSURF.352   
      enddo                                                                  BLSURF.353   
                                                                             BLSURF.354   
      end                                                                    BLSURF.355   
