phys/module_sf_sfclay.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 !
   3 MODULE module_sf_sfclay
   4
   5  REAL    , PARAMETER ::  VCONVC=1.
   6  REAL    , PARAMETER ::  CZO=0.0185
   7  REAL    , PARAMETER ::  OZO=1.59E-5
   8
   9  REAL,   DIMENSION(0:1000 ),SAVE          :: PSIMTB,PSIHTB
  10
  11 CONTAINS
  12
  13 !-------------------------------------------------------------------
  14    SUBROUTINE SFCLAY(U3D,V3D,T3D,QV3D,P3D,dz8w,                    &
  15                      CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM,      &
  16                      ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
  17                      FM,FH,                                        &
  18                      XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
  19                      U10,V10,TH2,T2,Q2,                            &
  20                      GZ1OZ0,WSPD,BR,ISFFLX,DX,                     &
  21                      SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
  22                      KARMAN,EOMEG,STBOLT,                          &
  23                      P1000mb,LAKEMASK,                             &
  24                      ids,ide, jds,jde, kds,kde,                    &
  25                      ims,ime, jms,jme, kms,kme,                    &
  26                      its,ite, jts,jte, kts,kte,                    &
  27                      ustm,ck,cka,cd,cda,                           &
  28                      isftcflx,iz0tlnd,scm_force_flux)
  29 !-------------------------------------------------------------------
  30       IMPLICIT NONE
  31 !-------------------------------------------------------------------
  32 !   Changes in V3.7 over water surfaces:
  33 !          1. for ZNT/Cd, replacing constant OZO with 0.11*1.5E-5/UST(I)
  34 !             the COARE 3.5 (Edson et al. 2013) formulation is also available
  35 !          2. for VCONV, reducing magnitude by half
  36 !          3. for Ck, replacing Carlson-Boland with COARE 3
  37 !-------------------------------------------------------------------
  38 !-- U3D         3D u-velocity interpolated to theta points (m/s)
  39 !-- V3D         3D v-velocity interpolated to theta points (m/s)
  40 !-- T3D         temperature (K)
  41 !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
  42 !-- P3D         3D pressure (Pa)
  43 !-- dz8w        dz between full levels (m)
  44 !-- CP          heat capacity at constant pressure for dry air (J/kg/K)
  45 !-- G           acceleration due to gravity (m/s^2)
  46 !-- ROVCP       R/CP
  47 !-- R           gas constant for dry air (J/kg/K)
  48 !-- XLV         latent heat of vaporization for water (J/kg)
  49 !-- PSFC        surface pressure (Pa)
  50 !-- ZNT         roughness length (m)
  51 !-- UST         u* in similarity theory (m/s)
  52 !-- USTM        u* in similarity theory (m/s) without vconv correction
  53 !               used to couple with TKE scheme
  54 !-- PBLH        PBL height from previous time (m)
  55 !-- MAVAIL      surface moisture availability (between 0 and 1)
  56 !-- ZOL         z/L height over Monin-Obukhov length
  57 !-- MOL         T* (similarity theory) (K)
  58 !-- REGIME      flag indicating PBL regime (stable, unstable, etc.)
  59 !-- PSIM        similarity stability function for momentum
  60 !-- PSIH        similarity stability function for heat
  61 !-- FM          integrated stability function for momentum
  62 !-- FH          integrated stability function for heat
  63 !-- XLAND       land mask (1 for land, 2 for water)
  64 !-- HFX         upward heat flux at the surface (W/m^2)
  65 !-- QFX         upward moisture flux at the surface (kg/m^2/s)
  66 !-- LH          net upward latent heat flux at surface (W/m^2)
  67 !-- TSK         surface temperature (K)
  68 !-- FLHC        exchange coefficient for heat (W/m^2/K)
  69 !-- FLQC        exchange coefficient for moisture (kg/m^2/s)
  70 !-- CHS         heat/moisture exchange coefficient for LSM (m/s)
  71 !-- QGH         lowest-level saturated mixing ratio
  72 !-- QSFC        ground saturated mixing ratio
  73 !-- U10         diagnostic 10m u wind
  74 !-- V10         diagnostic 10m v wind
  75 !-- TH2         diagnostic 2m theta (K)
  76 !-- T2          diagnostic 2m temperature (K)
  77 !-- Q2          diagnostic 2m mixing ratio (kg/kg)
  78 !-- GZ1OZ0      log(z/z0) where z0 is roughness length
  79 !-- WSPD        wind speed at lowest model level (m/s)
  80 !-- BR          bulk Richardson number in surface layer
  81 !-- ISFFLX      isfflx=1 for surface heat and moisture fluxes
  82 !-- DX          horizontal grid size (m)
  83 !-- SVP1        constant for saturation vapor pressure (kPa)
  84 !-- SVP2        constant for saturation vapor pressure (dimensionless)
  85 !-- SVP3        constant for saturation vapor pressure (K)
  86 !-- SVPT0       constant for saturation vapor pressure (K)
  87 !-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
  88 !-- EP2         constant for specific humidity calculation
  89 !               (R_d/R_v) (dimensionless)
  90 !-- KARMAN      Von Karman constant
  91 !-- EOMEG       angular velocity of earth's rotation (rad/s)
  92 !-- STBOLT      Stefan-Boltzmann constant (W/m^2/K^4)
  93 !-- ck          enthalpy exchange coeff at 10 meters
  94 !-- cd          momentum exchange coeff at 10 meters
  95 !-- cka         enthalpy exchange coeff at the lowest model level
  96 !-- cda         momentum exchange coeff at the lowest model level
  97 !-- isftcflx    =0, (Charnock and Carlson-Boland); =1, AHW Ck, Cd, =2 Garratt
  98 !-- iz0tlnd     =0 Carlson-Boland, =1 Czil_new
  99 !-- ids         start index for i in domain
 100 !-- ide         end index for i in domain
 101 !-- jds         start index for j in domain
 102 !-- jde         end index for j in domain
 103 !-- kds         start index for k in domain
 104 !-- kde         end index for k in domain
 105 !-- ims         start index for i in memory
 106 !-- ime         end index for i in memory
 107 !-- jms         start index for j in memory
 108 !-- jme         end index for j in memory
 109 !-- kms         start index for k in memory
 110 !-- kme         end index for k in memory
 111 !-- its         start index for i in tile
 112 !-- ite         end index for i in tile
 113 !-- jts         start index for j in tile
 114 !-- jte         end index for j in tile
 115 !-- kts         start index for k in tile
 116 !-- kte         end index for k in tile
 117 !-------------------------------------------------------------------
 118       INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
 119                                         ims,ime, jms,jme, kms,kme, &
 120                                         its,ite, jts,jte, kts,kte
 121 !
 122       INTEGER,  INTENT(IN )   ::        ISFFLX
 123       REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
 124       REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
 125       REAL,     INTENT(IN )   ::        P1000mb
 126 !
 127       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 128                 INTENT(IN   )   ::                           dz8w
 129
 130       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 131                 INTENT(IN   )   ::                           QV3D, &
 132                                                               P3D, &
 133                                                               T3D
 134
 135       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 136                 INTENT(IN   )               ::             MAVAIL, &
 137                                                              PBLH, &
 138                                                             XLAND, &
 139                                                          LAKEMASK, &
 140                                                               TSK
 141
 142       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 143                 INTENT(INOUT)               ::             REGIME, &
 144                                                               HFX, &
 145                                                               QFX, &
 146                                                                LH, &
 147                                                           MOL,RMOL
 148 !
 149       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 150                 INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
 151                                                   PSIM,PSIH,FM,FH
 152
 153       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 154                 INTENT(IN   )   ::                            U3D, &
 155                                                               V3D
 156
 157       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 158                 INTENT(IN   )               ::               PSFC
 159
 160       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 161                 INTENT(INOUT)   ::                            ZNT, &
 162                                                               ZOL, &
 163                                                               UST, &
 164                                                               CPM, &
 165                                                              CHS2, &
 166                                                              CQS2, &
 167                                                               CHS
 168
 169       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 170                 INTENT(INOUT)   ::                      FLHC,FLQC
 171
 172       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 173                 INTENT(INOUT)   ::                                 &
 174                                                               QGH
 175       REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV
 176
 177       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 178                 INTENT(IN   )               ::                 DX
 179
 180       REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme )              , &
 181                 INTENT(OUT)     ::                  ck,cka,cd,cda
 182
 183       REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme )              , &
 184                 INTENT(INOUT)   ::                           USTM
 185
 186       INTEGER,  OPTIONAL,  INTENT(IN )   ::     ISFTCFLX, IZ0TLND
 187       INTEGER,  OPTIONAL,  INTENT(IN )   ::     SCM_FORCE_FLUX
 188
 189       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 190                 INTENT(INOUT  )             ::               QSFC
 191
 192       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 193                 INTENT(OUT  )               ::                U10, &
 194                                                               V10, &
 195                                                               TH2, &
 196                                                                T2, &
 197                                                                Q2
 198
 199 ! LOCAL VARS
 200
 201       REAL,     DIMENSION( its:ite ) ::                       U1D, &
 202                                                               V1D, &
 203                                                              QV1D, &
 204                                                               P1D, &
 205                                                               T1D
 206
 207       REAL,     DIMENSION( its:ite ) ::                    dz8w1d
 208
 209       REAL,     DIMENSION( its:ite ) ::                      DX2D
 210
 211       INTEGER ::  I,J
 212
 213       DO J=jts,jte
 214
 215         DO i=its,ite
 216            DX2D(i)=DX(i,j)
 217         ENDDO
 218
 219         DO i=its,ite
 220           dz8w1d(I) = dz8w(i,1,j)
 221         ENDDO
 222
 223         DO i=its,ite
 224            U1D(i) =U3D(i,1,j)
 225            V1D(i) =V3D(i,1,j)
 226            QV1D(i)=QV3D(i,1,j)
 227            P1D(i) =P3D(i,1,j)
 228            T1D(i) =T3D(i,1,j)
 229         ENDDO
 230
 231         !  Sending array starting locations of optional variables may cause
 232         !  troubles, so we explicitly change the call.
 233
 234         CALL SFCLAY1D(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,               &
 235                 CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),&
 236                 CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j),   &
 237                 ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j),    &
 238                 MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j),  &
 239                 FM(ims,j),FH(ims,j),                               &
 240                 XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j),     &
 241                 U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j),        &
 242                 Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j),      &
 243                 QSFC(ims,j),LH(ims,j),                             &
 244                 GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX2D,   &
 245                 SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,  &
 246                 P1000mb,LAKEMASK(ims,j),                           &
 247                 ids,ide, jds,jde, kds,kde,                         &
 248                 ims,ime, jms,jme, kms,kme,                         &
 249                 its,ite, jts,jte, kts,kte                          &
 250 #if ( ( EM_CORE == 1 ) || ( defined(mpas) ) )
 251                 ,isftcflx,iz0tlnd,scm_force_flux,                  &
 252                 USTM(ims,j),CK(ims,j),CKA(ims,j),                  &
 253                 CD(ims,j),CDA(ims,j)                               &
 254 #endif
 255                                                                    )
 256       ENDDO
 257
 258
 259    END SUBROUTINE SFCLAY
 260
 261
 262 !-------------------------------------------------------------------
 263    SUBROUTINE SFCLAY1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d,                &
 264                      CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
 265                      ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH,FM,FH,&
 266                      XLAND,HFX,QFX,TSK,                            &
 267                      U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH,              &
 268                      QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX,             &
 269                      SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
 270                      KARMAN,EOMEG,STBOLT,                          &
 271                      P1000mb,LAKEMASK,                             &
 272                      ids,ide, jds,jde, kds,kde,                    &
 273                      ims,ime, jms,jme, kms,kme,                    &
 274                      its,ite, jts,jte, kts,kte,                    &
 275                      isftcflx, iz0tlnd, scm_force_flux,            &
 276                      ustm,ck,cka,cd,cda                            )
 277 !-------------------------------------------------------------------
 278       IMPLICIT NONE
 279 !-------------------------------------------------------------------
 280       REAL,     PARAMETER     ::        XKA=2.4E-5
 281       REAL,     PARAMETER     ::        PRT=1.
 282       REAL,     PARAMETER     ::        SALINITY_FACTOR=0.98
 283
 284       INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
 285                                         ims,ime, jms,jme, kms,kme, &
 286                                         its,ite, jts,jte, kts,kte, &
 287                                         J
 288 !
 289       INTEGER,  INTENT(IN )   ::        ISFFLX
 290       REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
 291       REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
 292       REAL,     INTENT(IN )   ::        P1000mb
 293
 294 !
 295       REAL,     DIMENSION( ims:ime )                             , &
 296                 INTENT(IN   )               ::             MAVAIL, &
 297                                                              PBLH, &
 298                                                             XLAND, &
 299                                                          LAKEMASK, &
 300                                                               TSK
 301 !
 302       REAL,     DIMENSION( ims:ime )                             , &
 303                 INTENT(IN   )               ::             PSFCPA
 304
 305       REAL,     DIMENSION( ims:ime )                             , &
 306                 INTENT(INOUT)               ::             REGIME, &
 307                                                               HFX, &
 308                                                               QFX, &
 309                                                          MOL,RMOL
 310 !m the following 5 are changed to memory size---
 311 !
 312       REAL,     DIMENSION( ims:ime )                             , &
 313                 INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
 314                                                   PSIM,PSIH,FM,FH
 315
 316       REAL,     DIMENSION( ims:ime )                             , &
 317                 INTENT(INOUT)   ::                            ZNT, &
 318                                                               ZOL, &
 319                                                               UST, &
 320                                                               CPM, &
 321                                                              CHS2, &
 322                                                              CQS2, &
 323                                                               CHS
 324
 325       REAL,     DIMENSION( ims:ime )                             , &
 326                 INTENT(INOUT)   ::                      FLHC,FLQC
 327
 328       REAL,     DIMENSION( ims:ime )                             , &
 329                 INTENT(INOUT)   ::                                 &
 330                                                          QSFC,QGH
 331
 332       REAL,     DIMENSION( ims:ime )                             , &
 333                 INTENT(OUT)     ::                        U10,V10, &
 334                                                      TH2,T2,Q2,LH
 335
 336       REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV
 337
 338       REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::      DX
 339
 340 ! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE SFCLAY
 341       REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::  dz8w1d
 342
 343       REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::      UX, &
 344                                                                VX, &
 345                                                              QV1D, &
 346                                                               P1D, &
 347                                                               T1D
 348
 349       REAL, OPTIONAL, DIMENSION( ims:ime )                       , &
 350                 INTENT(OUT)     ::                  ck,cka,cd,cda
 351       REAL, OPTIONAL, DIMENSION( ims:ime )                       , &
 352                 INTENT(INOUT)   ::                           USTM
 353
 354       INTEGER,  OPTIONAL,  INTENT(IN )   ::     ISFTCFLX, IZ0TLND
 355       INTEGER,  OPTIONAL,  INTENT(IN )   ::     SCM_FORCE_FLUX
 356
 357 ! LOCAL VARS
 358
 359       REAL,     DIMENSION( its:ite )        ::                 ZA, &
 360                                                         THVX,ZQKL, &
 361                                                            ZQKLP1, &
 362                                                            THX,QX, &
 363                                                             PSIH2, &
 364                                                             PSIM2, &
 365                                                            PSIH10, &
 366                                                            PSIM10, &
 367                                                            DENOMQ, &
 368                                                           DENOMQ2, &
 369                                                           DENOMT2, &
 370                                                             WSPDI, &
 371                                                            GZ2OZ0, &
 372                                                            GZ10OZ0
 373 !
 374       REAL,     DIMENSION( its:ite )        ::                     &
 375                                                       RHOX,GOVRTH, &
 376                                                             TGDSA
 377 !
 378       REAL,     DIMENSION( its:ite)         ::          SCR3,SCR4
 379       REAL,     DIMENSION( its:ite )        ::         THGB, PSFC
 380 !
 381       INTEGER                               ::                 KL
 382
 383       INTEGER ::  N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10
 384
 385       REAL    ::  PL,THCON,TVCON,E1
 386       REAL    ::  ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
 387       REAL    ::  DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2,PSIQ10
 388       REAL    ::  FLUXC,VSGD,Z0Q,VISC,RESTAR,CZIL,GZ0OZQ,GZ0OZT
 389       REAL    ::  ZW, ZN1, ZN2
 390       REAL    ::  Z0T, CZC
 391 !-------------------------------------------------------------------
 392       KL=kte
 393
 394       DO i=its,ite
 395 ! PSFC cb
 396          PSFC(I)=PSFCPA(I)/1000.
 397       ENDDO
 398 !
 399 !----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:
 400 !
 401       DO 5 I=its,ite
 402         TGDSA(I)=TSK(I)
 403 ! PSFC cb
 404 !        THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP
 405         THGB(I)=TSK(I)*(P1000mb/PSFCPA(I))**ROVCP
 406     5 CONTINUE
 407 !
 408 !-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
 409 !     T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.
 410 !
 411 !     *** NOTE ***
 412 !         THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,
 413 !         SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
 414 !         TENDENCIES.
 415 !
 416    10 CONTINUE
 417
 418 !     DO 24 I=its,ite
 419 !        UX(I)=U1D(I)
 420 !        VX(I)=V1D(I)
 421 !  24 CONTINUE
 422
 423    26 CONTINUE
 424
 425 !.....SCR3(I,K) STORE TEMPERATURE,
 426 !     SCR4(I,K) STORE VIRTUAL TEMPERATURE.
 427
 428       DO 30 I=its,ite
 429 ! PL cb
 430          PL=P1D(I)/1000.
 431          SCR3(I)=T1D(I)
 432 !         THCON=(100./PL)**ROVCP
 433          THCON=(P1000mb*0.001/PL)**ROVCP
 434          THX(I)=SCR3(I)*THCON
 435          SCR4(I)=SCR3(I)
 436          THVX(I)=THX(I)
 437          QX(I)=0.
 438    30 CONTINUE
 439 !
 440       DO I=its,ite
 441          QGH(I)=0.
 442          FLHC(I)=0.
 443          FLQC(I)=0.
 444          CPM(I)=CP
 445       ENDDO
 446 !
 447 !     IF(IDRY.EQ.1)GOTO 80
 448       DO 50 I=its,ite
 449          QX(I)=QV1D(I)
 450          TVCON=(1.+EP1*QX(I))
 451          THVX(I)=THX(I)*TVCON
 452          SCR4(I)=SCR3(I)*TVCON
 453    50 CONTINUE
 454 !
 455       DO 60 I=its,ite
 456         E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
 457 !  for land points QSFC can come from previous time step
 458 !  the saturation vapor pressure for salty water is on average 2% lower
 459         if(xland(i).gt.1.5 .and. lakemask(i).eq.0.) E1=E1*SALINITY_FACTOR
 460         if(xland(i).gt.1.5.or.qsfc(i).le.0.0)QSFC(I)=EP2*E1/(PSFC(I)-E1)
 461 ! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE
 462 ! Q2SAT = QGH IN LSM
 463         E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3))
 464         PL=P1D(I)/1000.
 465         QGH(I)=EP2*E1/(PL-E1)
 466         CPM(I)=CP*(1.+0.8*QX(I))
 467    60 CONTINUE
 468    80 CONTINUE
 469
 470 !-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND
 471 !     LEVEL, AND THE LAYER THICKNESSES.
 472
 473       DO 90 I=its,ite
 474         ZQKLP1(I)=0.
 475         RHOX(I)=PSFC(I)*1000./(R*SCR4(I))
 476    90 CONTINUE
 477 !
 478       DO 110 I=its,ite
 479            ZQKL(I)=dz8w1d(I)+ZQKLP1(I)
 480   110 CONTINUE
 481 !
 482       DO 120 I=its,ite
 483          ZA(I)=0.5*(ZQKL(I)+ZQKLP1(I))
 484   120 CONTINUE
 485 !
 486       DO 160 I=its,ite
 487         GOVRTH(I)=G/THX(I)
 488   160 CONTINUE
 489
 490 !-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO
 491 !     AKB(1976), EQ(12).
 492
 493       DO 260 I=its,ite
 494         GZ1OZ0(I)=ALOG(ZA(I)/ZNT(I))
 495         GZ2OZ0(I)=ALOG(2./ZNT(I))
 496         GZ10OZ0(I)=ALOG(10./ZNT(I))
 497         IF((XLAND(I)-1.5).GE.0)THEN
 498           ZL=ZNT(I)
 499         ELSE
 500           ZL=0.01
 501         ENDIF
 502         WSPD(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))
 503
 504         TSKV=THGB(I)*(1.+EP1*QSFC(I))
 505         DTHVDZ=(THVX(I)-TSKV)
 506 !  Convective velocity scale Vc and subgrid-scale velocity Vsg
 507 !  following Beljaars (1994, QJRMS) and Mahrt and Sun (1995, MWR)
 508 !                                ... HONG Aug. 2001
 509 !
 510 !       VCONV = 0.25*sqrt(g/tskv*pblh(i)*dthvm)
 511 !      Use Beljaars over land, old MM5 (Wyngaard) formula over water
 512         if (xland(i).lt.1.5) then
 513         fluxc = max(hfx(i)/rhox(i)/cp                    &
 514               + ep1*tskv*qfx(i)/rhox(i),0.)
 515         VCONV = vconvc*(g/tgdsa(i)*pblh(i)*fluxc)**.33
 516         else
 517         IF(-DTHVDZ.GE.0)THEN
 518           DTHVM=-DTHVDZ
 519         ELSE
 520           DTHVM=0.
 521         ENDIF
 522 !       VCONV = 2.*SQRT(DTHVM)
 523 ! V3.7: reducing contribution in calm conditions
 524         VCONV = SQRT(DTHVM)
 525         endif
 526 ! Mahrt and Sun low-res correction
 527         VSGD = 0.32 * (max(dx(i)/5000.-1.,0.))**.33
 528         WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV+vsgd*vsgd)
 529         WSPD(I)=AMAX1(WSPD(I),0.1)
 530         BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I))
 531 !  IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2
 532         IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0)
 533 !jdf
 534         RMOL(I)=-GOVRTH(I)*DTHVDZ*ZA(I)*KARMAN
 535 !jdf
 536
 537   260 CONTINUE
 538
 539 !
 540 !-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
 541 !
 542 !
 543 !     THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)
 544 !     AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).
 545 !
 546 !     CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
 547 !
 548 !        1. BR .GE. 0.2;
 549 !               REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
 550 !
 551 !        2. BR .LT. 0.2 .AND. BR .GT. 0.0;
 552 !               REPRESENTS DAMPED MECHANICAL TURBULENT CONDITIONS
 553 !               (REGIME=2),
 554 !
 555 !        3. BR .EQ. 0.0
 556 !               REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
 557 !
 558 !        4. BR .LT. 0.0
 559 !               REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
 560 !
 561 !CCCCC
 562
 563       DO 320 I=its,ite
 564 !CCCCC
 565 !CC     REMOVE REGIME 3 DEPENDENCE ON PBL HEIGHT
 566 !CC          IF(BR(I).LT.0..AND.HOL(I,J).GT.1.5)GOTO 310
 567         IF(BR(I).LT.0.)GOTO 310
 568 !
 569 !-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
 570 !
 571         IF(BR(I).LT.0.2)GOTO 270
 572         REGIME(I)=1.
 573         PSIM(I)=-10.*GZ1OZ0(I)
 574 !    LOWER LIMIT ON PSI IN STABLE CONDITIONS
 575         PSIM(I)=AMAX1(PSIM(I),-10.)
 576         PSIH(I)=PSIM(I)
 577         PSIM10(I)=10./ZA(I)*PSIM(I)
 578         PSIM10(I)=AMAX1(PSIM10(I),-10.)
 579         PSIH10(I)=PSIM10(I)
 580         PSIM2(I)=2./ZA(I)*PSIM(I)
 581         PSIM2(I)=AMAX1(PSIM2(I),-10.)
 582         PSIH2(I)=PSIM2(I)
 583
 584 !       1.0 over Monin-Obukhov length
 585         IF(UST(I).LT.0.01)THEN
 586            RMOL(I)=BR(I)*GZ1OZ0(I) !ZA/L
 587         ELSE
 588            RMOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I)) !ZA/L
 589         ENDIF
 590         RMOL(I)=AMIN1(RMOL(I),9.999) ! ZA/L
 591         RMOL(I) = RMOL(I)/ZA(I) !1.0/L
 592
 593         GOTO 320
 594 !
 595 !-----CLASS 2; DAMPED MECHANICAL TURBULENCE:
 596 !
 597   270   IF(BR(I).EQ.0.0)GOTO 280
 598         REGIME(I)=2.
 599         PSIM(I)=-5.0*BR(I)*GZ1OZ0(I)/(1.1-5.0*BR(I))
 600 !    LOWER LIMIT ON PSI IN STABLE CONDITIONS
 601         PSIM(I)=AMAX1(PSIM(I),-10.)
 602 !.....AKB(1976), EQ(16).
 603         PSIH(I)=PSIM(I)
 604         PSIM10(I)=10./ZA(I)*PSIM(I)
 605         PSIM10(I)=AMAX1(PSIM10(I),-10.)
 606         PSIH10(I)=PSIM10(I)
 607         PSIM2(I)=2./ZA(I)*PSIM(I)
 608         PSIM2(I)=AMAX1(PSIM2(I),-10.)
 609         PSIH2(I)=PSIM2(I)
 610
 611         ! Linear form: PSIM = -0.5*ZA/L; e.g, see eqn 16 of
 612         ! Blackadar, Modeling the nocturnal boundary layer, Preprints,
 613         ! Third Symposium on Atmospheric Turbulence Diffusion and Air Quality,
 614         ! Raleigh, NC, 1976
 615         ZOL(I) = BR(I)*GZ1OZ0(I)/(1.00001-5.0*BR(I))
 616
 617         if ( ZOL(I) .GT. 0.5 ) then ! linear form ok
 618            ! Holtslag and de Bruin, J. App. Meteor 27, 689-704, 1988;
 619            ! see also, Launiainen, Boundary-Layer Meteor 76,165-179, 1995
 620            ! Eqn (8) of Launiainen, 1995
 621            ZOL(I) = ( 1.89*GZ1OZ0(I) + 44.2 ) * BR(I)*BR(I)    &
 622                 + ( 1.18*GZ1OZ0(I) - 1.37 ) * BR(I)
 623            ZOL(I)=AMIN1(ZOL(I),9.999)
 624         end if
 625
 626         ! 1.0 over Monin-Obukhov length
 627         RMOL(I)= ZOL(I)/ZA(I)
 628
 629         GOTO 320
 630 !
 631 !-----CLASS 3; FORCED CONVECTION:
 632 !
 633   280   REGIME(I)=3.
 634         PSIM(I)=0.0
 635         PSIH(I)=PSIM(I)
 636         PSIM10(I)=0.
 637         PSIH10(I)=PSIM10(I)
 638         PSIM2(I)=0.
 639         PSIH2(I)=PSIM2(I)
 640
 641
 642         IF(UST(I).LT.0.01)THEN
 643           ZOL(I)=BR(I)*GZ1OZ0(I)
 644         ELSE
 645           ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I))
 646         ENDIF
 647
 648         RMOL(I) = ZOL(I)/ZA(I)
 649
 650         GOTO 320
 651 !
 652 !-----CLASS 4; FREE CONVECTION:
 653 !
 654   310   CONTINUE
 655         REGIME(I)=4.
 656         IF(UST(I).LT.0.01)THEN
 657           ZOL(I)=BR(I)*GZ1OZ0(I)
 658         ELSE
 659           ZOL(I)=KARMAN*GOVRTH(I)*ZA(I)*MOL(I)/(UST(I)*UST(I))
 660         ENDIF
 661         ZOL10=10./ZA(I)*ZOL(I)
 662         ZOL2=2./ZA(I)*ZOL(I)
 663         ZOL(I)=AMIN1(ZOL(I),0.)
 664         ZOL(I)=AMAX1(ZOL(I),-9.9999)
 665         ZOL10=AMIN1(ZOL10,0.)
 666         ZOL10=AMAX1(ZOL10,-9.9999)
 667         ZOL2=AMIN1(ZOL2,0.)
 668         ZOL2=AMAX1(ZOL2,-9.9999)
 669         NZOL=INT(-ZOL(I)*100.)
 670         RZOL=-ZOL(I)*100.-NZOL
 671         NZOL10=INT(-ZOL10*100.)
 672         RZOL10=-ZOL10*100.-NZOL10
 673         NZOL2=INT(-ZOL2*100.)
 674         RZOL2=-ZOL2*100.-NZOL2
 675         PSIM(I)=PSIMTB(NZOL)+RZOL*(PSIMTB(NZOL+1)-PSIMTB(NZOL))
 676         PSIH(I)=PSIHTB(NZOL)+RZOL*(PSIHTB(NZOL+1)-PSIHTB(NZOL))
 677         PSIM10(I)=PSIMTB(NZOL10)+RZOL10*(PSIMTB(NZOL10+1)-PSIMTB(NZOL10))
 678         PSIH10(I)=PSIHTB(NZOL10)+RZOL10*(PSIHTB(NZOL10+1)-PSIHTB(NZOL10))
 679         PSIM2(I)=PSIMTB(NZOL2)+RZOL2*(PSIMTB(NZOL2+1)-PSIMTB(NZOL2))
 680         PSIH2(I)=PSIHTB(NZOL2)+RZOL2*(PSIHTB(NZOL2+1)-PSIHTB(NZOL2))
 681
 682 !---LIMIT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS
 683 !---  THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL
 684 !       PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
 685 !       PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
 686         PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
 687         PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
 688         PSIH2(I)=AMIN1(PSIH2(I),0.9*GZ2OZ0(I))
 689         PSIM10(I)=AMIN1(PSIM10(I),0.9*GZ10OZ0(I))
 690 ! AHW: mods to compute ck, cd
 691         PSIH10(I)=AMIN1(PSIH10(I),0.9*GZ10OZ0(I))
 692
 693         RMOL(I) = ZOL(I)/ZA(I)
 694
 695   320 CONTINUE
 696 !
 697 !-----COMPUTE THE FRICTIONAL VELOCITY:
 698 !     ZA(1982) EQS(2.60),(2.61).
 699 !
 700       DO 330 I=its,ite
 701         DTG=THX(I)-THGB(I)
 702         PSIX=GZ1OZ0(I)-PSIM(I)
 703         PSIX10=GZ10OZ0(I)-PSIM10(I)
 704 !     LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL
 705 !     ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0
 706         PSIT=AMAX1(GZ1OZ0(I)-PSIH(I),2.)
 707
 708         IF((XLAND(I)-1.5).GE.0)THEN
 709           ZL=ZNT(I)
 710         ELSE
 711           ZL=0.01
 712         ENDIF
 713         PSIQ=ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I)
 714         PSIT2=GZ2OZ0(I)-PSIH2(I)
 715         PSIQ2=ALOG(KARMAN*UST(I)*2./XKA+2./ZL)-PSIH2(I)
 716 ! AHW: mods to compute ck, cd
 717         PSIQ10=ALOG(KARMAN*UST(I)*10./XKA+10./ZL)-PSIH10(I)
 718
 719 ! V3.7: using Fairall 2003 to compute z0q and z0t over water:
 720 !       adapted from module_sf_mynn.F
 721         IF ( (XLAND(I)-1.5).GE.0. ) THEN
 722               VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
 723 !             VISC=1.326e-5*(1. + 6.542e-3*SCR3(I) + 8.301e-6*SCR3(I)*SCR3(I) &
 724 !                 - 4.84e-9*SCR3(I)*SCR3(I)*SCR3(I))
 725               RESTAR=UST(I)*ZNT(I)/VISC
 726               Z0T = (5.5e-5)*(RESTAR**(-0.60))
 727               Z0T = MIN(Z0T,1.0e-4)
 728               Z0T = MAX(Z0T,2.0e-9)
 729               Z0Q = Z0T
 730
 731               PSIQ=max(ALOG((ZA(I)+Z0Q)/Z0Q)-PSIH(I), 2.)
 732               PSIT=max(ALOG((ZA(I)+Z0T)/Z0T)-PSIH(I), 2.)
 733               PSIQ2=max(ALOG((2.+Z0Q)/Z0Q)-PSIH2(I), 2.)
 734               PSIT2=max(ALOG((2.+Z0T)/Z0T)-PSIH2(I), 2.)
 735               PSIQ10=max(ALOG((10.+Z0Q)/Z0Q)-PSIH10(I), 2.)
 736         ENDIF
 737
 738         IF ( PRESENT(ISFTCFLX) ) THEN
 739            IF ( ISFTCFLX.EQ.1 .AND. (XLAND(I)-1.5).GE.0. ) THEN
 740 ! v3.1
 741 !             Z0Q = 1.e-4 + 1.e-3*(MAX(0.,UST(I)-1.))**2
 742 ! hfip1
 743 !             Z0Q = 0.62*2.0E-5/UST(I) + 1.E-3*(MAX(0.,UST(I)-1.5))**2
 744 ! v3.2
 745               Z0Q = 1.e-4
 746               PSIQ=ALOG(ZA(I)/Z0Q)-PSIH(I)
 747               PSIT=PSIQ
 748               PSIQ2=ALOG(2./Z0Q)-PSIH2(I)
 749               PSIQ10=ALOG(10./Z0Q)-PSIH10(I)
 750               PSIT2=PSIQ2
 751            ENDIF
 752            IF ( ISFTCFLX.EQ.2 .AND. (XLAND(I)-1.5).GE.0. ) THEN
 753 ! AHW: Garratt formula: Calculate roughness Reynolds number
 754 !        Kinematic viscosity of air (linear approc to
 755 !                 temp dependence at sea level)
 756 ! GZ0OZT and GZ0OZQ are based off formulas from Brutsaert (1975), which
 757 ! Garratt (1992) used with values of k = 0.40, Pr = 0.71, and Sc = 0.60
 758               VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
 759 !!            VISC=1.5E-5
 760               RESTAR=UST(I)*ZNT(I)/VISC
 761               GZ0OZT=0.40*(7.3*SQRT(SQRT(RESTAR))*SQRT(0.71)-5.)
 762               GZ0OZQ=0.40*(7.3*SQRT(SQRT(RESTAR))*SQRT(0.60)-5.)
 763               PSIT=GZ1OZ0(I)-PSIH(I)+GZ0OZT
 764               PSIQ=GZ1OZ0(I)-PSIH(I)+GZ0OZQ
 765               PSIT2=GZ2OZ0(I)-PSIH2(I)+GZ0OZT
 766               PSIQ2=GZ2OZ0(I)-PSIH2(I)+GZ0OZQ
 767               PSIQ10=GZ10OZ0(I)-PSIH(I)+GZ0OZQ
 768            ENDIF
 769         ENDIF
 770         IF(PRESENT(ck) .and. PRESENT(cd) .and. PRESENT(cka) .and. PRESENT(cda)) THEN
 771            Ck(I)=(karman/psix10)*(karman/psiq10)
 772            Cd(I)=(karman/psix10)*(karman/psix10)
 773            Cka(I)=(karman/psix)*(karman/psiq)
 774            Cda(I)=(karman/psix)*(karman/psix)
 775         ENDIF
 776         IF ( PRESENT(IZ0TLND) ) THEN
 777            IF ( IZ0TLND.GE.1 .AND. (XLAND(I)-1.5).LE.0. ) THEN
 778               ZL=ZNT(I)
 779 !             CZIL RELATED CHANGES FOR LAND
 780               VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
 781               RESTAR=UST(I)*ZL/VISC
 782 !             Modify CZIL according to Chen & Zhang, 2009 if iz0tlnd = 1
 783 !             If iz0tlnd = 2, use traditional value
 784
 785               IF ( IZ0TLND.EQ.1 ) THEN
 786                  CZIL = 10.0 ** ( -0.40 * ( ZL / 0.07 ) )
 787               ELSE IF ( IZ0TLND.EQ.2 ) THEN
 788                  CZIL = 0.1
 789               END IF
 790
 791               PSIT=GZ1OZ0(I)-PSIH(I)+CZIL*KARMAN*SQRT(RESTAR)
 792               PSIQ=GZ1OZ0(I)-PSIH(I)+CZIL*KARMAN*SQRT(RESTAR)
 793               PSIT2=GZ2OZ0(I)-PSIH2(I)+CZIL*KARMAN*SQRT(RESTAR)
 794               PSIQ2=GZ2OZ0(I)-PSIH2(I)+CZIL*KARMAN*SQRT(RESTAR)
 795
 796            ENDIF
 797         ENDIF
 798 ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
 799         UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX
 800 ! TKE coupling: compute ust without vconv for use in tke scheme
 801         WSPDI(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))
 802         IF ( PRESENT(USTM) ) THEN
 803         USTM(I)=0.5*USTM(I)+0.5*KARMAN*WSPDI(I)/PSIX
 804         ENDIF
 805         U10(I)=UX(I)*PSIX10/PSIX
 806         V10(I)=VX(I)*PSIX10/PSIX
 807         TH2(I)=THGB(I)+DTG*PSIT2/PSIT
 808         Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ
 809 !        T2(I) = TH2(I)*(PSFC(I)/100.)**ROVCP
 810         T2(I) = TH2(I)*(PSFCPA(I)/P1000mb)**ROVCP
 811 !       LATER Q2 WILL BE OVERWRITTEN FOR LAND POINTS IN SURFCE
 812 !       QA2(I,J) = Q2(I)
 813 !       UA10(I,J) = U10(I)
 814 !       VA10(I,J) = V10(I)
 815 !       write(*,1002)UST(I),KARMAN*WSPD(I),PSIX,KARMAN*WSPD(I)/PSIX
 816 !
 817         IF((XLAND(I)-1.5).LT.0.)THEN
 818           UST(I)=AMAX1(UST(I),0.1)
 819         ENDIF
 820         MOL(I)=KARMAN*DTG/PSIT/PRT
 821         DENOMQ(I)=PSIQ
 822         DENOMQ2(I)=PSIQ2
 823         DENOMT2(I)=PSIT2
 824         FM(I)=PSIX
 825         FH(I)=PSIT
 826   330 CONTINUE
 827 !
 828   335 CONTINUE
 829
 830 !-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES:
 831       IF ( PRESENT(SCM_FORCE_FLUX) ) THEN
 832          IF (SCM_FORCE_FLUX.EQ.1) GOTO 350
 833       ENDIF
 834       DO i=its,ite
 835         QFX(i)=0.
 836         HFX(i)=0.
 837       ENDDO
 838   350 CONTINUE
 839
 840       IF (ISFFLX.EQ.0) GOTO 410
 841
 842 !-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST).
 843
 844       DO 360 I=its,ite
 845         IF((XLAND(I)-1.5).GE.0)THEN
 846 !         ZNT(I)=CZO*UST(I)*UST(I)/G+OZO
 847 ! Since V3.7 (ref: EC Physics document for Cy36r1)
 848           ZNT(I)=CZO*UST(I)*UST(I)/G+0.11*1.5E-5/UST(I)
 849 ! V3.9: Add limit as in isftcflx = 1,2
 850           ZNT(I)=MIN(ZNT(I),2.85e-3)
 851 ! COARE 3.5 (Edson et al. 2013)
 852 !         CZC = 0.0017*WSPD(I)-0.005
 853 !         CZC = min(CZC,0.028)
 854 !         ZNT(I)=CZC*UST(I)*UST(I)/G+0.11*1.5E-5/UST(I)
 855 ! AHW: change roughness length, and hence the drag coefficients Ck and Cd
 856           IF ( PRESENT(ISFTCFLX) ) THEN
 857              IF ( ISFTCFLX.NE.0 ) THEN
 858 !               ZNT(I)=10.*exp(-9.*UST(I)**(-.3333))
 859 !               ZNT(I)=10.*exp(-9.5*UST(I)**(-.3333))
 860 !               ZNT(I)=ZNT(I) + 0.11*1.5E-5/AMAX1(UST(I),0.01)
 861 !               ZNT(I)=0.011*UST(I)*UST(I)/G+OZO
 862 !               ZNT(I)=MAX(ZNT(I),3.50e-5)
 863 ! AHW 2012:
 864                 ZW  = MIN((UST(I)/1.06)**(0.3),1.0)
 865                 ZN1 = 0.011*UST(I)*UST(I)/G + OZO
 866                 ZN2 = 10.*exp(-9.5*UST(I)**(-.3333)) + &
 867                        0.11*1.5E-5/AMAX1(UST(I),0.01)
 868                 ZNT(I)=(1.0-ZW) * ZN1 + ZW * ZN2
 869                 ZNT(I)=MIN(ZNT(I),2.85e-3)
 870                 ZNT(I)=MAX(ZNT(I),1.27e-7)
 871              ENDIF
 872           ENDIF
 873           ZL = ZNT(I)
 874         ELSE
 875           ZL = 0.01
 876         ENDIF
 877         FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/DENOMQ(I)
 878 !       FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/(   &
 879 !               ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I))
 880         DTTHX=ABS(THX(I)-THGB(I))
 881         IF(DTTHX.GT.1.E-5)THEN
 882           FLHC(I)=CPM(I)*RHOX(I)*UST(I)*MOL(I)/(THX(I)-THGB(I))
 883 !         write(*,1001)FLHC(I),CPM(I),RHOX(I),UST(I),MOL(I),THX(I),THGB(I),I
 884  1001   format(f8.5,2x,f12.7,2x,f12.10,2x,f12.10,2x,f13.10,2x,f12.8,f12.8,2x,i3)
 885         ELSE
 886           FLHC(I)=0.
 887         ENDIF
 888   360 CONTINUE
 889
 890 !
 891 !-----COMPUTE SURFACE MOIST FLUX:
 892 !
 893 !     IF(IDRY.EQ.1)GOTO 390
 894      IF ( PRESENT(SCM_FORCE_FLUX) ) THEN
 895         IF (SCM_FORCE_FLUX.EQ.1) GOTO 405
 896      ENDIF
 897 !
 898       DO 370 I=its,ite
 899         QFX(I)=FLQC(I)*(QSFC(I)-QX(I))
 900 !       QFX(I)=AMAX1(QFX(I),0.)
 901         LH(I)=XLV*QFX(I)
 902   370 CONTINUE
 903
 904 !-----COMPUTE SURFACE HEAT FLUX:
 905 !
 906   390 CONTINUE
 907       DO 400 I=its,ite
 908         IF(XLAND(I)-1.5.GT.0.)THEN
 909           HFX(I)=FLHC(I)*(THGB(I)-THX(I))
 910 !         IF ( PRESENT(ISFTCFLX) ) THEN
 911 !            IF ( ISFTCFLX.NE.0 ) THEN
 912 ! AHW: add dissipative heating term (commented out in 3.6.1)
 913 !               HFX(I)=HFX(I)+RHOX(I)*USTM(I)*USTM(I)*WSPDI(I)
 914 !            ENDIF
 915 !         ENDIF
 916         ELSEIF(XLAND(I)-1.5.LT.0.)THEN
 917           HFX(I)=FLHC(I)*(THGB(I)-THX(I))
 918 !         HFX(I)=AMAX1(HFX(I),-250.)
 919         ENDIF
 920   400 CONTINUE
 921
 922   405 CONTINUE
 923
 924       DO I=its,ite
 925          IF((XLAND(I)-1.5).GE.0)THEN
 926            ZL=ZNT(I)
 927          ELSE
 928            ZL=0.01
 929          ENDIF
 930          CHS(I)=UST(I)*KARMAN/DENOMQ(I)
 931 !        GZ2OZ0(I)=ALOG(2./ZNT(I))
 932 !        PSIM2(I)=-10.*GZ2OZ0(I)
 933 !        PSIM2(I)=AMAX1(PSIM2(I),-10.)
 934 !        PSIH2(I)=PSIM2(I)
 935          CQS2(I)=UST(I)*KARMAN/DENOMQ2(I)
 936          CHS2(I)=UST(I)*KARMAN/DENOMT2(I)
 937       ENDDO
 938
 939   410 CONTINUE
 940 !jdf
 941 !     DO I=its,ite
 942 !       IF(UST(I).GE.0.1) THEN
 943 !         RMOL(I)=RMOL(I)*(-FLHC(I))/(UST(I)*UST(I)*UST(I))
 944 !       ELSE
 945 !         RMOL(I)=RMOL(I)*(-FLHC(I))/(0.1*0.1*0.1)
 946 !       ENDIF
 947 !     ENDDO
 948 !jdf
 949
 950 !
 951    END SUBROUTINE SFCLAY1D
 952
 953 !====================================================================
 954    SUBROUTINE sfclayinit( allowed_to_read )
 955
 956    LOGICAL , INTENT(IN)      ::      allowed_to_read
 957    INTEGER                   ::      N
 958    REAL                      ::      ZOLN,X,Y
 959
 960    DO N=0,1000
 961       ZOLN=-FLOAT(N)*0.01
 962       X=(1-16.*ZOLN)**0.25
 963       PSIMTB(N)=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))- &
 964                 2.*ATAN(X)+2.*ATAN(1.)
 965       Y=(1-16*ZOLN)**0.5
 966       PSIHTB(N)=2*ALOG(0.5*(1+Y))
 967    ENDDO
 968
 969    END SUBROUTINE sfclayinit
 970
 971 !-------------------------------------------------------------------
 972
 973 END MODULE module_sf_sfclay