phys/module_sf_sfclayrev.F

   1 !WRF:MODEL_LAYER:PHYSICS
   2 !
   3 MODULE module_sf_sfclayrev
   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 :: psim_stab,psim_unstab,psih_stab,psih_unstab
  10
  11 CONTAINS
  12
  13 !-------------------------------------------------------------------
  14    SUBROUTINE SFCLAYREV(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,                                      &
  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,isftcflx,iz0tlnd,          &
  28                      shalwater_z0,water_depth,shalwater_depth,     &
  29                      scm_force_flux                                )
  30 !-------------------------------------------------------------------
  31       IMPLICIT NONE
  32 !-------------------------------------------------------------------
  33 !   Changes in V3.7 over water surfaces:
  34 !          1. for ZNT/Cd, replacing constant OZO with 0.11*1.5E-5/UST(I)
  35 !             the COARE 3.5 (Edson et al. 2013) formulation is also available
  36 !          2. for VCONV, reducing magnitude by half
  37 !          3. for Ck, replacing Carlson-Boland with COARE 3
  38 !-------------------------------------------------------------------
  39 !-- U3D         3D u-velocity interpolated to theta points (m/s)
  40 !-- V3D         3D v-velocity interpolated to theta points (m/s)
  41 !-- T3D         temperature (K)
  42 !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
  43 !-- P3D         3D pressure (Pa)
  44 !-- dz8w        dz between full levels (m)
  45 !-- CP          heat capacity at constant pressure for dry air (J/kg/K)
  46 !-- G           acceleration due to gravity (m/s^2)
  47 !-- ROVCP       R/CP
  48 !-- R           gas constant for dry air (J/kg/K)
  49 !-- XLV         latent heat of vaporization for water (J/kg)
  50 !-- PSFC        surface pressure (Pa)
  51 !-- ZNT         roughness length (m)
  52 !-- UST         u* in similarity theory (m/s)
  53 !-- USTM        u* in similarity theory (m/s) without vconv correction
  54 !               used to couple with TKE scheme
  55 !-- PBLH        PBL height from previous time (m)
  56 !-- MAVAIL      surface moisture availability (between 0 and 1)
  57 !-- ZOL         z/L height over Monin-Obukhov length
  58 !-- MOL         T* (similarity theory) (K)
  59 !-- REGIME      flag indicating PBL regime (stable, unstable, etc.)
  60 !-- PSIM        similarity stability function for momentum
  61 !-- PSIH        similarity stability function for heat
  62 !-- FM          integrated stability function for momentum
  63 !-- FH          integrated stability function for heat
  64 !-- XLAND       land mask (1 for land, 2 for water)
  65 !-- HFX         upward heat flux at the surface (W/m^2)
  66 !-- QFX         upward moisture flux at the surface (kg/m^2/s)
  67 !-- LH          net upward latent heat flux at surface (W/m^2)
  68 !-- TSK         surface temperature (K)
  69 !-- FLHC        exchange coefficient for heat (W/m^2/K)
  70 !-- FLQC        exchange coefficient for moisture (kg/m^2/s)
  71 !-- CHS         heat/moisture exchange coefficient for LSM (m/s)
  72 !-- QGH         lowest-level saturated mixing ratio
  73 !-- QSFC        ground saturated mixing ratio
  74 !-- U10         diagnostic 10m u wind
  75 !-- V10         diagnostic 10m v wind
  76 !-- TH2         diagnostic 2m theta (K)
  77 !-- T2          diagnostic 2m temperature (K)
  78 !-- Q2          diagnostic 2m mixing ratio (kg/kg)
  79 !-- GZ1OZ0      log(z/z0) where z0 is roughness length
  80 !-- WSPD        wind speed at lowest model level (m/s)
  81 !-- BR          bulk Richardson number in surface layer
  82 !-- ISFFLX      isfflx=1 for surface heat and moisture fluxes
  83 !-- DX          horizontal grid size (m)
  84 !-- SVP1        constant for saturation vapor pressure (kPa)
  85 !-- SVP2        constant for saturation vapor pressure (dimensionless)
  86 !-- SVP3        constant for saturation vapor pressure (K)
  87 !-- SVPT0       constant for saturation vapor pressure (K)
  88 !-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
  89 !-- EP2         constant for specific humidity calculation
  90 !               (R_d/R_v) (dimensionless)
  91 !-- KARMAN      Von Karman constant
  92 !-- EOMEG       angular velocity of earth's rotation (rad/s)
  93 !-- STBOLT      Stefan-Boltzmann constant (W/m^2/K^4)
  94 !-- ck          enthalpy exchange coeff at 10 meters
  95 !-- cd          momentum exchange coeff at 10 meters
  96 !-- cka         enthalpy exchange coeff at the lowest model level
  97 !-- cda         momentum exchange coeff at the lowest model level
  98 !-- isftcflx    =0, (Charnock and Carlson-Boland); =1, AHW Ck, Cd, =2 Garratt
  99 !-- iz0tlnd     =0 Carlson-Boland, =1 Czil_new
 100 !-- ids         start index for i in domain
 101 !-- ide         end index for i in domain
 102 !-- jds         start index for j in domain
 103 !-- jde         end index for j in domain
 104 !-- kds         start index for k in domain
 105 !-- kde         end index for k in domain
 106 !-- ims         start index for i in memory
 107 !-- ime         end index for i in memory
 108 !-- jms         start index for j in memory
 109 !-- jme         end index for j in memory
 110 !-- kms         start index for k in memory
 111 !-- kme         end index for k in memory
 112 !-- its         start index for i in tile
 113 !-- ite         end index for i in tile
 114 !-- jts         start index for j in tile
 115 !-- jte         end index for j in tile
 116 !-- kts         start index for k in tile
 117 !-- kte         end index for k in tile
 118 !-------------------------------------------------------------------
 119       INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
 120                                         ims,ime, jms,jme, kms,kme, &
 121                                         its,ite, jts,jte, kts,kte
 122 !
 123       INTEGER,  INTENT(IN )   ::        ISFFLX
 124       REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
 125       REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
 126       REAL,     INTENT(IN )   ::        P1000mb
 127 !
 128       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 129                 INTENT(IN   )   ::                           dz8w
 130
 131       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 132                 INTENT(IN   )   ::                           QV3D, &
 133                                                               P3D, &
 134                                                               T3D
 135
 136       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 137                 INTENT(IN   )               ::             MAVAIL, &
 138                                                              PBLH, &
 139                                                             XLAND, &
 140                                                               TSK
 141       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 142                 INTENT(OUT  )               ::                U10, &
 143                                                               V10, &
 144                                                               TH2, &
 145                                                                T2, &
 146                                                                Q2
 147
 148       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 149                 INTENT(INOUT)               ::             REGIME, &
 150                                                               HFX, &
 151                                                               QFX, &
 152                                                                LH, &
 153                                                              QSFC, &
 154                                                           MOL,RMOL
 155 !m the following 5 are change to memory size
 156 !
 157       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 158                 INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
 159                                                   PSIM,PSIH,FM,FH
 160
 161       REAL,     DIMENSION( ims:ime, kms:kme, jms:jme )           , &
 162                 INTENT(IN   )   ::                            U3D, &
 163                                                               V3D
 164
 165       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 166                 INTENT(IN   )               ::               PSFC
 167
 168       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 169                 INTENT(INOUT)   ::                            ZNT, &
 170                                                               ZOL, &
 171                                                               UST, &
 172                                                               CPM, &
 173                                                              CHS2, &
 174                                                              CQS2, &
 175                                                               CHS
 176
 177       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 178                 INTENT(INOUT)   ::                      FLHC,FLQC
 179
 180       REAL,     DIMENSION( ims:ime, jms:jme )                    , &
 181                 INTENT(INOUT)   ::                                 &
 182                                                               QGH
 183
 184       REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV,DX
 185
 186       REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme )              , &
 187                 INTENT(OUT)     ::                  ck,cka,cd,cda
 188
 189       REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme )              , &
 190                 INTENT(INOUT)   ::                           USTM
 191
 192       INTEGER,  OPTIONAL,  INTENT(IN )   ::     ISFTCFLX, IZ0TLND
 193       INTEGER,  OPTIONAL,  INTENT(IN )   ::     SCM_FORCE_FLUX
 194
 195       INTEGER,  INTENT(IN )   ::     shalwater_z0
 196       REAL,     INTENT(IN )   ::     shalwater_depth
 197       REAL,     DIMENSION( ims:ime, jms:jme ), INTENT(IN)  :: water_depth
 198 ! LOCAL VARS
 199
 200       REAL,     DIMENSION( its:ite ) ::                       U1D, &
 201                                                               V1D, &
 202                                                              QV1D, &
 203                                                               P1D, &
 204                                                               T1D
 205
 206       REAL,     DIMENSION( its:ite ) ::                    dz8w1d
 207
 208       INTEGER ::  I,J
 209
 210       DO J=jts,jte
 211         DO i=its,ite
 212           dz8w1d(I) = dz8w(i,1,j)
 213         ENDDO
 214
 215         DO i=its,ite
 216            U1D(i) =U3D(i,1,j)
 217            V1D(i) =V3D(i,1,j)
 218            QV1D(i)=QV3D(i,1,j)
 219            P1D(i) =P3D(i,1,j)
 220            T1D(i) =T3D(i,1,j)
 221         ENDDO
 222
 223         !  Sending array starting locations of optional variables may cause
 224         !  troubles, so we explicitly change the call.
 225
 226         CALL SFCLAYREV1D(J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,               &
 227                 CP,G,ROVCP,R,XLV,PSFC(ims,j),CHS(ims,j),CHS2(ims,j),&
 228                 CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j),   &
 229                 ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j),    &
 230                 MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j),  &
 231                 FM(ims,j),FH(ims,j),                               &
 232                 XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j),     &
 233                 U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j),        &
 234                 Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j),      &
 235                 QSFC(ims,j),LH(ims,j),                             &
 236                 GZ1OZ0(ims,j),WSPD(ims,j),BR(ims,j),ISFFLX,DX,     &
 237                 SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN,EOMEG,STBOLT,  &
 238                 P1000mb,                                           &
 239                 shalwater_z0,water_depth(ims,j),shalwater_depth,   &
 240                 ids,ide, jds,jde, kds,kde,                         &
 241                 ims,ime, jms,jme, kms,kme,                         &
 242                 its,ite, jts,jte, kts,kte                          &
 243 #if ( EM_CORE == 1 )
 244                 ,isftcflx,iz0tlnd,scm_force_flux,                               &
 245                 USTM(ims,j),CK(ims,j),CKA(ims,j),                  &
 246                 CD(ims,j),CDA(ims,j)                               &
 247 #endif
 248                                                                    )
 249       ENDDO
 250
 251
 252    END SUBROUTINE SFCLAYREV
 253
 254
 255 !-------------------------------------------------------------------
 256    SUBROUTINE SFCLAYREV1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d,                &
 257                      CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
 258                      ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH,FM,FH,&
 259                      XLAND,HFX,QFX,TSK,                            &
 260                      U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH,              &
 261                      QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX,             &
 262                      SVP1,SVP2,SVP3,SVPT0,EP1,EP2,                 &
 263                      KARMAN,EOMEG,STBOLT,                          &
 264                      P1000mb,                                      &
 265                      shalwater_z0,water_depth,shalwater_depth,     &
 266                      ids,ide, jds,jde, kds,kde,                    &
 267                      ims,ime, jms,jme, kms,kme,                    &
 268                      its,ite, jts,jte, kts,kte,                    &
 269                      isftcflx, iz0tlnd,scm_force_flux,                            &
 270                      ustm,ck,cka,cd,cda                            )
 271 !-------------------------------------------------------------------
 272       IMPLICIT NONE
 273 !-------------------------------------------------------------------
 274       REAL,     PARAMETER     ::        XKA=2.4E-5
 275       REAL,     PARAMETER     ::        PRT=1.
 276
 277       INTEGER,  INTENT(IN )   ::        ids,ide, jds,jde, kds,kde, &
 278                                         ims,ime, jms,jme, kms,kme, &
 279                                         its,ite, jts,jte, kts,kte, &
 280                                         J
 281 !
 282       INTEGER,  INTENT(IN )   ::        ISFFLX
 283       REAL,     INTENT(IN )   ::        SVP1,SVP2,SVP3,SVPT0
 284       REAL,     INTENT(IN )   ::        EP1,EP2,KARMAN,EOMEG,STBOLT
 285       REAL,     INTENT(IN )   ::        P1000mb
 286
 287 !
 288       REAL,     DIMENSION( ims:ime )                             , &
 289                 INTENT(IN   )               ::             MAVAIL, &
 290                                                              PBLH, &
 291                                                             XLAND, &
 292                                                               TSK
 293 !
 294       REAL,     DIMENSION( ims:ime )                             , &
 295                 INTENT(IN   )               ::             PSFCPA
 296
 297       REAL,     DIMENSION( ims:ime )                             , &
 298                 INTENT(INOUT)               ::             REGIME, &
 299                                                               HFX, &
 300                                                               QFX, &
 301                                                          MOL,RMOL
 302 !m the following 5 are changed to memory size---
 303 !
 304       REAL,     DIMENSION( ims:ime )                             , &
 305                 INTENT(INOUT)   ::                 GZ1OZ0,WSPD,BR, &
 306                                                   PSIM,PSIH,FM,FH
 307
 308       REAL,     DIMENSION( ims:ime )                             , &
 309                 INTENT(INOUT)   ::                            ZNT, &
 310                                                               ZOL, &
 311                                                               UST, &
 312                                                               CPM, &
 313                                                              CHS2, &
 314                                                              CQS2, &
 315                                                               CHS
 316
 317       REAL,     DIMENSION( ims:ime )                             , &
 318                 INTENT(INOUT)   ::                      FLHC,FLQC
 319
 320       REAL,     DIMENSION( ims:ime )                             , &
 321                 INTENT(INOUT)   ::                                 &
 322                                                          QSFC,QGH
 323
 324       REAL,     DIMENSION( ims:ime )                             , &
 325                 INTENT(OUT)     ::                        U10,V10, &
 326                                                      TH2,T2,Q2,LH
 327
 328
 329       REAL,     INTENT(IN   )               ::   CP,G,ROVCP,R,XLV,DX
 330
 331       INTEGER,  INTENT(IN )   ::     shalwater_z0
 332       REAL,     INTENT(IN )   ::     shalwater_depth
 333       REAL,     DIMENSION( ims:ime ), INTENT(IN)  :: water_depth
 334 ! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE SFCLAY
 335       REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::  dz8w1d
 336
 337       REAL,     DIMENSION( its:ite ),  INTENT(IN   )   ::      UX, &
 338                                                                VX, &
 339                                                              QV1D, &
 340                                                               P1D, &
 341                                                               T1D
 342
 343       REAL, OPTIONAL, DIMENSION( ims:ime )                       , &
 344                 INTENT(OUT)     ::                  ck,cka,cd,cda
 345       REAL, OPTIONAL, DIMENSION( ims:ime )                       , &
 346                 INTENT(INOUT)   ::                           USTM
 347
 348       INTEGER,  OPTIONAL,  INTENT(IN )   ::     ISFTCFLX, IZ0TLND
 349       INTEGER,  OPTIONAL,  INTENT(IN )   ::     SCM_FORCE_FLUX
 350
 351 ! LOCAL VARS
 352
 353       REAL,     DIMENSION( its:ite )        ::                 ZA, &
 354                                                         THVX,ZQKL, &
 355                                                            ZQKLP1, &
 356                                                            THX,QX, &
 357                                                             PSIH2, &
 358                                                             PSIM2, &
 359                                                            PSIH10, &
 360                                                            PSIM10, &
 361                                                            DENOMQ, &
 362                                                           DENOMQ2, &
 363                                                           DENOMT2, &
 364                                                             WSPDI, &
 365                                                            GZ2OZ0, &
 366                                                            GZ10OZ0
 367 !
 368       REAL,     DIMENSION( its:ite )        ::                     &
 369                                                       RHOX,GOVRTH, &
 370                                                             TGDSA
 371 !
 372       REAL,     DIMENSION( its:ite)         ::          SCR3,SCR4
 373       REAL,     DIMENSION( its:ite )        ::         THGB, PSFC
 374 !
 375       INTEGER                               ::                 KL
 376
 377       INTEGER ::  N,I,K,KK,L,NZOL,NK,NZOL2,NZOL10
 378
 379       REAL    ::  PL,THCON,TVCON,E1
 380       REAL    ::  ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
 381       REAL    ::  DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2,PSIQ10
 382       REAL    ::  FLUXC,VSGD,Z0Q,VISC,RESTAR,CZIL,GZ0OZQ,GZ0OZT
 383       REAL    ::  ZW, ZN1, ZN2
 384 !
 385 ! .... paj ...
 386 !
 387       REAL    :: zolzz,zol0
 388 !     REAL    :: zolri,zolri2
 389 !     REAL    :: psih_stable,psim_stable,psih_unstable,psim_unstable
 390 !     REAL    :: psih_stable_full,psim_stable_full,psih_unstable_full,psim_unstable_full
 391       REAL    :: zl2,zl10,z0t
 392       REAL,     DIMENSION( its:ite )        ::         pq,pq2,pq10
 393
 394
 395 !-------------------------------------------------------------------
 396       KL=kte
 397
 398       DO i=its,ite
 399 ! PSFC cb
 400          PSFC(I)=PSFCPA(I)/1000.
 401       ENDDO
 402 !
 403 !----CONVERT GROUND TEMPERATURE TO POTENTIAL TEMPERATURE:
 404 !
 405       DO 5 I=its,ite
 406         TGDSA(I)=TSK(I)
 407 ! PSFC cb
 408 !        THGB(I)=TSK(I)*(100./PSFC(I))**ROVCP
 409         THGB(I)=TSK(I)*(P1000mb/PSFCPA(I))**ROVCP
 410     5 CONTINUE
 411 !
 412 !-----DECOUPLE FLUX-FORM VARIABLES TO GIVE U,V,T,THETA,THETA-VIR.,
 413 !     T-VIR., QV, AND QC AT CROSS POINTS AND AT KTAU-1.
 414 !
 415 !     *** NOTE ***
 416 !         THE BOUNDARY WINDS MAY NOT BE ADEQUATELY AFFECTED BY FRICTION,
 417 !         SO USE ONLY INTERIOR VALUES OF UX AND VX TO CALCULATE
 418 !         TENDENCIES.
 419 !
 420    10 CONTINUE
 421
 422 !     DO 24 I=its,ite
 423 !        UX(I)=U1D(I)
 424 !        VX(I)=V1D(I)
 425 !  24 CONTINUE
 426
 427    26 CONTINUE
 428
 429 !.....SCR3(I,K) STORE TEMPERATURE,
 430 !     SCR4(I,K) STORE VIRTUAL TEMPERATURE.
 431
 432       DO 30 I=its,ite
 433 ! PL cb
 434          PL=P1D(I)/1000.
 435          SCR3(I)=T1D(I)
 436 !         THCON=(100./PL)**ROVCP
 437          THCON=(P1000mb*0.001/PL)**ROVCP
 438          THX(I)=SCR3(I)*THCON
 439          SCR4(I)=SCR3(I)
 440          THVX(I)=THX(I)
 441          QX(I)=0.
 442    30 CONTINUE
 443 !
 444       DO I=its,ite
 445          QGH(I)=0.
 446          FLHC(I)=0.
 447          FLQC(I)=0.
 448          CPM(I)=CP
 449       ENDDO
 450 !
 451 !     IF(IDRY.EQ.1)GOTO 80
 452       DO 50 I=its,ite
 453          QX(I)=QV1D(I)
 454          TVCON=(1.+EP1*QX(I))
 455          THVX(I)=THX(I)*TVCON
 456          SCR4(I)=SCR3(I)*TVCON
 457    50 CONTINUE
 458 !
 459       DO 60 I=its,ite
 460         E1=SVP1*EXP(SVP2*(TGDSA(I)-SVPT0)/(TGDSA(I)-SVP3))
 461 !  for land points QSFC can come from previous time step
 462         if(xland(i).gt.1.5.or.qsfc(i).le.0.0)QSFC(I)=EP2*E1/(PSFC(I)-E1)
 463 ! QGH CHANGED TO USE LOWEST-LEVEL AIR TEMP CONSISTENT WITH MYJSFC CHANGE
 464 ! Q2SAT = QGH IN LSM
 465         E1=SVP1*EXP(SVP2*(T1D(I)-SVPT0)/(T1D(I)-SVP3))
 466         PL=P1D(I)/1000.
 467         QGH(I)=EP2*E1/(PL-E1)
 468         CPM(I)=CP*(1.+0.8*QX(I))
 469    60 CONTINUE
 470    80 CONTINUE
 471
 472 !-----COMPUTE THE HEIGHT OF FULL- AND HALF-SIGMA LEVELS ABOVE GROUND
 473 !     LEVEL, AND THE LAYER THICKNESSES.
 474
 475       DO 90 I=its,ite
 476         ZQKLP1(I)=0.
 477         RHOX(I)=PSFC(I)*1000./(R*SCR4(I))
 478    90 CONTINUE
 479 !
 480       DO 110 I=its,ite
 481            ZQKL(I)=dz8w1d(I)+ZQKLP1(I)
 482   110 CONTINUE
 483 !
 484       DO 120 I=its,ite
 485          ZA(I)=0.5*(ZQKL(I)+ZQKLP1(I))
 486   120 CONTINUE
 487 !
 488       DO 160 I=its,ite
 489         GOVRTH(I)=G/THX(I)
 490   160 CONTINUE
 491
 492 !-----CALCULATE BULK RICHARDSON NO. OF SURFACE LAYER, ACCORDING TO
 493 !     AKB(1976), EQ(12).
 494
 495       DO 260 I=its,ite
 496         GZ1OZ0(I)=ALOG((ZA(I)+ZNT(I))/ZNT(I))   ! log((z+z0)/z0)
 497         GZ2OZ0(I)=ALOG((2.+ZNT(I))/ZNT(I))      ! log((2+z0)/z0)
 498         GZ10OZ0(I)=ALOG((10.+ZNT(I))/ZNT(I))    ! log((10+z0)z0)
 499         IF((XLAND(I)-1.5).GE.0)THEN
 500           ZL=ZNT(I)
 501         ELSE
 502           ZL=0.01
 503         ENDIF
 504         WSPD(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))
 505
 506         TSKV=THGB(I)*(1.+EP1*QSFC(I))
 507         DTHVDZ=(THVX(I)-TSKV)
 508 !  Convective velocity scale Vc and subgrid-scale velocity Vsg
 509 !  following Beljaars (1994, QJRMS) and Mahrt and Sun (1995, MWR)
 510 !                                ... HONG Aug. 2001
 511 !
 512 !       VCONV = 0.25*sqrt(g/tskv*pblh(i)*dthvm)
 513 !      Use Beljaars over land, old MM5 (Wyngaard) formula over water
 514         if (xland(i).lt.1.5) then
 515         fluxc = max(hfx(i)/rhox(i)/cp                    &
 516               + ep1*tskv*qfx(i)/rhox(i),0.)
 517         VCONV = vconvc*(g/tgdsa(i)*pblh(i)*fluxc)**.33
 518         else
 519         IF(-DTHVDZ.GE.0)THEN
 520           DTHVM=-DTHVDZ
 521         ELSE
 522           DTHVM=0.
 523         ENDIF
 524 !       VCONV = 2.*SQRT(DTHVM)
 525 ! V3.7: reducing contribution in calm conditions
 526         VCONV = SQRT(DTHVM)
 527         endif
 528 ! Mahrt and Sun low-res correction
 529         VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
 530         WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV+vsgd*vsgd)
 531         WSPD(I)=AMAX1(WSPD(I),0.1)
 532         BR(I)=GOVRTH(I)*ZA(I)*DTHVDZ/(WSPD(I)*WSPD(I))
 533 !  IF PREVIOUSLY UNSTABLE, DO NOT LET INTO REGIMES 1 AND 2
 534         IF(MOL(I).LT.0.)BR(I)=AMIN1(BR(I),0.0)
 535 !jdf
 536         RMOL(I)=-GOVRTH(I)*DTHVDZ*ZA(I)*KARMAN
 537 !jdf
 538
 539   260 CONTINUE
 540
 541 !
 542 !-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
 543 !
 544 !
 545 !     THE STABILITY CLASSES ARE DETERMINED BY BR (BULK RICHARDSON NO.)
 546 !     AND HOL (HEIGHT OF PBL/MONIN-OBUKHOV LENGTH).
 547 !
 548 !     CRITERIA FOR THE CLASSES ARE AS FOLLOWS:
 549 !
 550 !        1. BR .GE. 0.0;
 551 !               REPRESENTS NIGHTTIME STABLE CONDITIONS (REGIME=1),
 552 !
 553 !        3. BR .EQ. 0.0
 554 !               REPRESENTS FORCED CONVECTION CONDITIONS (REGIME=3),
 555 !
 556 !        4. BR .LT. 0.0
 557 !               REPRESENTS FREE CONVECTION CONDITIONS (REGIME=4).
 558 !
 559 !CCCCC
 560
 561       DO 320 I=its,ite
 562 !
 563       if (br(I).gt.0) then
 564         if (br(I).gt.250.0) then
 565         zol(I)=zolri(250.0,ZA(I),ZNT(I))
 566         else
 567         zol(I)=zolri(br(I),ZA(I),ZNT(I))
 568         endif
 569       endif
 570 !
 571       if (br(I).lt.0) then
 572        IF(UST(I).LT.0.001)THEN
 573           ZOL(I)=BR(I)*GZ1OZ0(I)
 574         ELSE
 575         if (br(I).lt.-250.0) then
 576         zol(I)=zolri(-250.0,ZA(I),ZNT(I))
 577         else
 578         zol(I)=zolri(br(I),ZA(I),ZNT(I))
 579         endif
 580        ENDIF
 581       endif
 582 !
 583 ! ... paj: compute integrated similarity functions.
 584 !
 585         zolzz=zol(I)*(za(I)+znt(I))/za(I) ! (z+z0/L
 586         zol10=zol(I)*(10.+znt(I))/za(I)   ! (10+z0)/L
 587         zol2=zol(I)*(2.+znt(I))/za(I)     ! (2+z0)/L
 588         zol0=zol(I)*znt(I)/za(I)          ! z0/L
 589         ZL2=(2.)/ZA(I)*ZOL(I)             ! 2/L
 590         ZL10=(10.)/ZA(I)*ZOL(I)           ! 10/L
 591
 592         IF((XLAND(I)-1.5).LT.0.)THEN
 593         ZL=(0.01)/ZA(I)*ZOL(I)   ! (0.01)/L
 594         ELSE
 595         ZL=ZOL0                     ! z0/L
 596         ENDIF
 597
 598         IF(BR(I).LT.0.)GOTO 310  ! go to unstable regime (class 4)
 599         IF(BR(I).EQ.0.)GOTO 280  ! go to neutral regime (class 3)
 600 !
 601 !-----CLASS 1; STABLE (NIGHTTIME) CONDITIONS:
 602 !
 603         REGIME(I)=1.
 604 !
 605 ! ... paj: psim and psih. Follows Cheng and Brutsaert 2005 (CB05).
 606 !
 607         psim(I)=psim_stable(zolzz)-psim_stable(zol0)
 608         psih(I)=psih_stable(zolzz)-psih_stable(zol0)
 609 !
 610         psim10(I)=psim_stable(zol10)-psim_stable(zol0)
 611         psih10(I)=psih_stable(zol10)-psih_stable(zol0)
 612 !
 613         psim2(I)=psim_stable(zol2)-psim_stable(zol0)
 614         psih2(I)=psih_stable(zol2)-psih_stable(zol0)
 615 !
 616 ! ... paj: preparations to compute PSIQ. Follows CB05+Carlson Boland JAM 1978.
 617 !
 618         pq(I)=psih_stable(zol(I))-psih_stable(zl)
 619         pq2(I)=psih_stable(zl2)-psih_stable(zl)
 620         pq10(I)=psih_stable(zl10)-psih_stable(zl)
 621 !
 622 !       1.0 over Monin-Obukhov length
 623         RMOL(I)=ZOL(I)/ZA(I)
 624 !
 625
 626         GOTO 320
 627 !
 628 !-----CLASS 3; FORCED CONVECTION:
 629 !
 630   280   REGIME(I)=3.
 631         PSIM(I)=0.0
 632         PSIH(I)=PSIM(I)
 633         PSIM10(I)=0.
 634         PSIH10(I)=PSIM10(I)
 635         PSIM2(I)=0.
 636         PSIH2(I)=PSIM2(I)
 637 !
 638 ! paj: preparations to compute PSIQ.
 639 !
 640         pq(I)=PSIH(I)
 641         pq2(I)=PSIH2(I)
 642         pq10(I)=0.
 643 !
 644         ZOL(I)=0.
 645         RMOL(I) = ZOL(I)/ZA(I)
 646
 647         GOTO 320
 648 !
 649 !-----CLASS 4; FREE CONVECTION:
 650 !
 651   310   CONTINUE
 652         REGIME(I)=4.
 653 !
 654 ! ... paj: PSIM and PSIH ...
 655 !
 656         psim(I)=psim_unstable(zolzz)-psim_unstable(zol0)
 657         psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
 658 !
 659         psim10(I)=psim_unstable(zol10)-psim_unstable(zol0)
 660         psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
 661 !
 662         psim2(I)=psim_unstable(zol2)-psim_unstable(zol0)
 663         psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
 664 !
 665 ! ... paj: preparations to compute PSIQ
 666 !
 667         pq(I)=psih_unstable(zol(I))-psih_unstable(zl)
 668         pq2(I)=psih_unstable(zl2)-psih_unstable(zl)
 669         pq10(I)=psih_unstable(zl10)-psih_unstable(zl)
 670 !
 671 !---LIMIOT PSIH AND PSIM IN THE CASE OF THIN LAYERS AND HIGH ROUGHNESS
 672 !---  THIS PREVENTS DENOMINATOR IN FLUXES FROM GETTING TOO SMALL
 673         PSIH(I)=AMIN1(PSIH(I),0.9*GZ1OZ0(I))
 674         PSIM(I)=AMIN1(PSIM(I),0.9*GZ1OZ0(I))
 675         PSIH2(I)=AMIN1(PSIH2(I),0.9*GZ2OZ0(I))
 676         PSIM10(I)=AMIN1(PSIM10(I),0.9*GZ10OZ0(I))
 677 !
 678 ! AHW: mods to compute ck, cd
 679         PSIH10(I)=AMIN1(PSIH10(I),0.9*GZ10OZ0(I))
 680
 681         RMOL(I) = ZOL(I)/ZA(I)
 682
 683   320 CONTINUE
 684 !
 685 !-----COMPUTE THE FRICTIONAL VELOCITY:
 686 !     ZA(1982) EQS(2.60),(2.61).
 687 !
 688       DO 330 I=its,ite
 689         DTG=THX(I)-THGB(I)
 690         PSIX=GZ1OZ0(I)-PSIM(I)
 691         PSIX10=GZ10OZ0(I)-PSIM10(I)
 692
 693 !     LOWER LIMIT ADDED TO PREVENT LARGE FLHC IN SOIL MODEL
 694 !     ACTIVATES IN UNSTABLE CONDITIONS WITH THIN LAYERS OR HIGH Z0
 695 !       PSIT=AMAX1(GZ1OZ0(I)-PSIH(I),2.)
 696        PSIT=GZ1OZ0(I)-PSIH(I)
 697        PSIT2=GZ2OZ0(I)-PSIH2(I)
 698 !
 699         IF((XLAND(I)-1.5).GE.0)THEN
 700           ZL=ZNT(I)
 701         ELSE
 702           ZL=0.01
 703         ENDIF
 704 !
 705         PSIQ=ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-pq(I)
 706         PSIQ2=ALOG(KARMAN*UST(I)*2./XKA+2./ZL)-pq2(I)
 707
 708 ! AHW: mods to compute ck, cd
 709         PSIQ10=ALOG(KARMAN*UST(I)*10./XKA+10./ZL)-pq10(I)
 710
 711 ! V3.7: using Fairall 2003 to compute z0q and z0t over water:
 712 !       adapted from module_sf_mynn.F
 713         IF ( (XLAND(I)-1.5).GE.0. ) THEN
 714               VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
 715               RESTAR=UST(I)*ZNT(I)/VISC
 716               Z0T = (5.5e-5)*(RESTAR**(-0.60))
 717               Z0T = MIN(Z0T,1.0e-4)
 718               Z0T = MAX(Z0T,2.0e-9)
 719               Z0Q = Z0T
 720
 721 ! following paj:
 722            zolzz=zol(I)*(za(I)+z0t)/za(I)    ! (z+z0t)/L
 723            zol10=zol(I)*(10.+z0t)/za(I)   ! (10+z0t)/L
 724            zol2=zol(I)*(2.+z0t)/za(I)     ! (2+z0t)/L
 725            zol0=zol(I)*z0t/za(I)          ! z0t/L
 726 !
 727               if (zol(I).gt.0.) then
 728               psih(I)=psih_stable(zolzz)-psih_stable(zol0)
 729               psih10(I)=psih_stable(zol10)-psih_stable(zol0)
 730               psih2(I)=psih_stable(zol2)-psih_stable(zol0)
 731               else
 732                 if (zol(I).eq.0) then
 733                 psih(I)=0.
 734                 psih10(I)=0.
 735                 psih2(I)=0.
 736                 else
 737                 psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
 738                 psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
 739                 psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
 740                 endif
 741               endif
 742               PSIT=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I)
 743               PSIT2=ALOG((2.+z0t)/Z0t)-PSIH2(I)
 744
 745            zolzz=zol(I)*(za(I)+z0q)/za(I)    ! (z+z0q)/L
 746            zol10=zol(I)*(10.+z0q)/za(I)   ! (10+z0q)/L
 747            zol2=zol(I)*(2.+z0q)/za(I)     ! (2+z0q)/L
 748            zol0=zol(I)*z0q/za(I)          ! z0q/L
 749 !
 750               if (zol(I).gt.0.) then
 751               psih(I)=psih_stable(zolzz)-psih_stable(zol0)
 752               psih10(I)=psih_stable(zol10)-psih_stable(zol0)
 753               psih2(I)=psih_stable(zol2)-psih_stable(zol0)
 754               else
 755                 if (zol(I).eq.0) then
 756                 psih(I)=0.
 757                 psih10(I)=0.
 758                 psih2(I)=0.
 759                 else
 760                 psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
 761                 psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
 762                 psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
 763                 endif
 764               endif
 765 !
 766               PSIQ=ALOG((ZA(I)+z0q)/Z0q)-PSIH(I)
 767               PSIQ2=ALOG((2.+z0q)/Z0q)-PSIH2(I)
 768               PSIQ10=ALOG((10.+z0q)/Z0q)-PSIH10(I)
 769         ENDIF
 770
 771         IF ( PRESENT(ISFTCFLX) ) THEN
 772            IF ( ISFTCFLX.EQ.1 .AND. (XLAND(I)-1.5).GE.0. ) THEN
 773 ! v3.1
 774 !             Z0Q = 1.e-4 + 1.e-3*(MAX(0.,UST(I)-1.))**2
 775 ! hfip1
 776 !             Z0Q = 0.62*2.0E-5/UST(I) + 1.E-3*(MAX(0.,UST(I)-1.5))**2
 777 ! v3.2
 778               Z0Q = 1.e-4
 779 !
 780 ! ... paj: recompute psih for z0q
 781 !
 782            zolzz=zol(I)*(za(I)+z0q)/za(I)    ! (z+z0q)/L
 783            zol10=zol(I)*(10.+z0q)/za(I)   ! (10+z0q)/L
 784            zol2=zol(I)*(2.+z0q)/za(I)     ! (2+z0q)/L
 785            zol0=zol(I)*z0q/za(I)          ! z0q/L
 786 !
 787               if (zol(I).gt.0.) then
 788               psih(I)=psih_stable(zolzz)-psih_stable(zol0)
 789               psih10(I)=psih_stable(zol10)-psih_stable(zol0)
 790               psih2(I)=psih_stable(zol2)-psih_stable(zol0)
 791               else
 792                 if (zol(I).eq.0) then
 793                 psih(I)=0.
 794                 psih10(I)=0.
 795                 psih2(I)=0.
 796                 else
 797                 psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
 798                 psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
 799                 psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
 800                 endif
 801               endif
 802 !
 803               PSIQ=ALOG((ZA(I)+z0q)/Z0Q)-PSIH(I)
 804               PSIT=PSIQ
 805               PSIQ2=ALOG((2.+z0q)/Z0Q)-PSIH2(I)
 806               PSIQ10=ALOG((10.+z0q)/Z0Q)-PSIH10(I)
 807               PSIT2=PSIQ2
 808            ENDIF
 809           IF ( ISFTCFLX.EQ.2 .AND. (XLAND(I)-1.5).GE.0. ) THEN
 810 ! AHW: Garratt formula: Calculate roughness Reynolds number
 811 !        Kinematic viscosity of air (linear approc to
 812 !                 temp dependence at sea level)
 813 ! GZ0OZT and GZ0OZQ are based off formulas from Brutsaert (1975), which
 814 ! Garratt (1992) used with values of k = 0.40, Pr = 0.71, and Sc = 0.60
 815               VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
 816 !!            VISC=1.5E-5
 817               RESTAR=UST(I)*ZNT(I)/VISC
 818               GZ0OZT=0.40*(7.3*SQRT(SQRT(RESTAR))*SQRT(0.71)-5.)
 819 !
 820 ! ... paj: compute psih for z0t for temperature ...
 821 !
 822               z0t=znt(I)/exp(GZ0OZT)
 823 !
 824            zolzz=zol(I)*(za(I)+z0t)/za(I)    ! (z+z0t)/L
 825            zol10=zol(I)*(10.+z0t)/za(I)   ! (10+z0t)/L
 826            zol2=zol(I)*(2.+z0t)/za(I)     ! (2+z0t)/L
 827            zol0=zol(I)*z0t/za(I)          ! z0t/L
 828 !
 829               if (zol(I).gt.0.) then
 830               psih(I)=psih_stable(zolzz)-psih_stable(zol0)
 831               psih10(I)=psih_stable(zol10)-psih_stable(zol0)
 832               psih2(I)=psih_stable(zol2)-psih_stable(zol0)
 833               else
 834                 if (zol(I).eq.0) then
 835                 psih(I)=0.
 836                 psih10(I)=0.
 837                 psih2(I)=0.
 838                 else
 839                 psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
 840                 psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
 841                 psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
 842                 endif
 843               endif
 844 !
 845 !              PSIT=GZ1OZ0(I)-PSIH(I)+RESTAR2
 846 !              PSIT2=GZ2OZ0(I)-PSIH2(I)+RESTAR2
 847               PSIT=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I)
 848               PSIT2=ALOG((2.+z0t)/Z0t)-PSIH2(I)
 849 !
 850               GZ0OZQ=0.40*(7.3*SQRT(SQRT(RESTAR))*SQRT(0.60)-5.)
 851               z0q=znt(I)/exp(GZ0OZQ)
 852 !
 853            zolzz=zol(I)*(za(I)+z0q)/za(I)    ! (z+z0q)/L
 854            zol10=zol(I)*(10.+z0q)/za(I)   ! (10+z0q)/L
 855            zol2=zol(I)*(2.+z0q)/za(I)     ! (2+z0q)/L
 856            zol0=zol(I)*z0q/za(I)          ! z0q/L
 857 !
 858               if (zol(I).gt.0.) then
 859               psih(I)=psih_stable(zolzz)-psih_stable(zol0)
 860               psih10(I)=psih_stable(zol10)-psih_stable(zol0)
 861               psih2(I)=psih_stable(zol2)-psih_stable(zol0)
 862               else
 863                 if (zol(I).eq.0) then
 864                 psih(I)=0.
 865                 psih10(I)=0.
 866                 psih2(I)=0.
 867                 else
 868                 psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
 869                 psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
 870                 psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
 871                 endif
 872               endif
 873 !
 874               PSIQ=ALOG((ZA(I)+z0q)/Z0q)-PSIH(I)
 875               PSIQ2=ALOG((2.+z0q)/Z0q)-PSIH2(I)
 876               PSIQ10=ALOG((10.+z0q)/Z0q)-PSIH10(I)
 877 !              PSIQ=GZ1OZ0(I)-PSIH(I)+2.28*SQRT(SQRT(RESTAR))-2.
 878 !              PSIQ2=GZ2OZ0(I)-PSIH2(I)+2.28*SQRT(SQRT(RESTAR))-2.
 879 !              PSIQ10=GZ10OZ0(I)-PSIH(I)+2.28*SQRT(SQRT(RESTAR))-2.
 880            ENDIF
 881         ENDIF
 882         IF(PRESENT(ck) .and. PRESENT(cd) .and. PRESENT(cka) .and. PRESENT(cda)) THEN
 883            Ck(I)=(karman/psix10)*(karman/psiq10)
 884            Cd(I)=(karman/psix10)*(karman/psix10)
 885            Cka(I)=(karman/psix)*(karman/psiq)
 886            Cda(I)=(karman/psix)*(karman/psix)
 887         ENDIF
 888         IF ( PRESENT(IZ0TLND) ) THEN
 889            IF ( IZ0TLND.GE.1 .AND. (XLAND(I)-1.5).LE.0. ) THEN
 890               ZL=ZNT(I)
 891 !             CZIL RELATED CHANGES FOR LAND
 892               VISC=(1.32+0.009*(SCR3(I)-273.15))*1.E-5
 893               RESTAR=UST(I)*ZL/VISC
 894 !             Modify CZIL according to Chen & Zhang, 2009 if iz0tlnd = 1
 895 !             If iz0tlnd = 2, use traditional value
 896
 897               IF ( IZ0TLND.EQ.1 ) THEN
 898                  CZIL = 10.0 ** ( -0.40 * ( ZL / 0.07 ) )
 899               ELSE IF ( IZ0TLND.EQ.2 ) THEN
 900                  CZIL = 0.1
 901               END IF
 902 !
 903 ! ... paj: compute phish for z0t over land
 904 !
 905               z0t=znt(I)/exp(CZIL*KARMAN*SQRT(RESTAR))
 906 !
 907            zolzz=zol(I)*(za(I)+z0t)/za(I)    ! (z+z0t)/L
 908            zol10=zol(I)*(10.+z0t)/za(I)   ! (10+z0t)/L
 909            zol2=zol(I)*(2.+z0t)/za(I)     ! (2+z0t)/L
 910            zol0=zol(I)*z0t/za(I)          ! z0t/L
 911 !
 912               if (zol(I).gt.0.) then
 913               psih(I)=psih_stable(zolzz)-psih_stable(zol0)
 914               psih10(I)=psih_stable(zol10)-psih_stable(zol0)
 915               psih2(I)=psih_stable(zol2)-psih_stable(zol0)
 916               else
 917                 if (zol(I).eq.0) then
 918                 psih(I)=0.
 919                 psih10(I)=0.
 920                 psih2(I)=0.
 921                 else
 922                 psih(I)=psih_unstable(zolzz)-psih_unstable(zol0)
 923                 psih10(I)=psih_unstable(zol10)-psih_unstable(zol0)
 924                 psih2(I)=psih_unstable(zol2)-psih_unstable(zol0)
 925                 endif
 926               endif
 927 !
 928               PSIQ=ALOG((ZA(I)+z0t)/Z0t)-PSIH(I)
 929               PSIQ2=ALOG((2.+z0t)/Z0t)-PSIH2(I)
 930               PSIT=PSIQ
 931               PSIT2=PSIQ2
 932 !
 933 !              PSIT=GZ1OZ0(I)-PSIH(I)+CZIL*KARMAN*SQRT(RESTAR)
 934 !              PSIQ=GZ1OZ0(I)-PSIH(I)+CZIL*KARMAN*SQRT(RESTAR)
 935 !              PSIT2=GZ2OZ0(I)-PSIH2(I)+CZIL*KARMAN*SQRT(RESTAR)
 936 !              PSIQ2=GZ2OZ0(I)-PSIH2(I)+CZIL*KARMAN*SQRT(RESTAR)
 937
 938            ENDIF
 939         ENDIF
 940 ! TO PREVENT OSCILLATIONS AVERAGE WITH OLD VALUE
 941         UST(I)=0.5*UST(I)+0.5*KARMAN*WSPD(I)/PSIX
 942 ! TKE coupling: compute ust without vconv for use in tke scheme
 943         WSPDI(I)=SQRT(UX(I)*UX(I)+VX(I)*VX(I))
 944         IF ( PRESENT(USTM) ) THEN
 945         USTM(I)=0.5*USTM(I)+0.5*KARMAN*WSPDI(I)/PSIX
 946         ENDIF
 947
 948         U10(I)=UX(I)*PSIX10/PSIX
 949         V10(I)=VX(I)*PSIX10/PSIX
 950         TH2(I)=THGB(I)+DTG*PSIT2/PSIT
 951         Q2(I)=QSFC(I)+(QX(I)-QSFC(I))*PSIQ2/PSIQ
 952         T2(I) = TH2(I)*(PSFCPA(I)/P1000mb)**ROVCP
 953 !
 954         IF((XLAND(I)-1.5).LT.0.)THEN
 955           UST(I)=AMAX1(UST(I),0.001)
 956         ENDIF
 957         MOL(I)=KARMAN*DTG/PSIT/PRT
 958         DENOMQ(I)=PSIQ
 959         DENOMQ2(I)=PSIQ2
 960         DENOMT2(I)=PSIT2
 961         FM(I)=PSIX
 962         FH(I)=PSIT
 963   330 CONTINUE
 964 !
 965   335 CONTINUE
 966
 967 !-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES:
 968       IF ( PRESENT(SCM_FORCE_FLUX) ) THEN
 969          IF (SCM_FORCE_FLUX.EQ.1) GOTO 350
 970       ENDIF
 971       DO i=its,ite
 972         QFX(i)=0.
 973         HFX(i)=0.
 974       ENDDO
 975   350 CONTINUE
 976
 977       IF (ISFFLX.EQ.0) GOTO 410
 978
 979 !-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST).
 980
 981       DO 360 I=its,ite
 982         IF((XLAND(I)-1.5).GE.0)THEN
 983 !         ZNT(I)=CZO*UST(I)*UST(I)/G+OZO
 984           ! PSH - formulation for depth-dependent roughness from
 985           ! ... Jimenez and Dudhia, 2018
 986           IF ( shalwater_z0 .eq. 1 ) THEN
 987              ZNT(I) = depth_dependent_z0(water_depth(I),ZNT(I),UST(I))
 988           ELSE
 989              ! Since V3.7 (ref: EC Physics document for Cy36r1)
 990              ZNT(I)=CZO*UST(I)*UST(I)/G+0.11*1.5E-5/UST(I)
 991              ! V3.9: Add limit as in isftcflx = 1,2
 992              ZNT(I)=MIN(ZNT(I),2.85e-3)
 993           ENDIF
 994 ! COARE 3.5 (Edson et al. 2013)
 995 !         CZC = 0.0017*WSPD(I)-0.005
 996 !         CZC = min(CZC,0.028)
 997 !         ZNT(I)=CZC*UST(I)*UST(I)/G+0.11*1.5E-5/UST(I)
 998 ! AHW: change roughness length, and hence the drag coefficients Ck and Cd
 999           IF ( PRESENT(ISFTCFLX) ) THEN
1000              IF ( ISFTCFLX.NE.0 ) THEN
1001 !               ZNT(I)=10.*exp(-9.*UST(I)**(-.3333))
1002 !               ZNT(I)=10.*exp(-9.5*UST(I)**(-.3333))
1003 !               ZNT(I)=ZNT(I) + 0.11*1.5E-5/AMAX1(UST(I),0.01)
1004 !               ZNT(I)=0.011*UST(I)*UST(I)/G+OZO
1005 !               ZNT(I)=MAX(ZNT(I),3.50e-5)
1006 ! AHW 2012:
1007                 ZW  = MIN((UST(I)/1.06)**(0.3),1.0)
1008                 ZN1 = 0.011*UST(I)*UST(I)/G + OZO
1009                 ZN2 = 10.*exp(-9.5*UST(I)**(-.3333)) + &
1010                        0.11*1.5E-5/AMAX1(UST(I),0.01)
1011                 ZNT(I)=(1.0-ZW) * ZN1 + ZW * ZN2
1012                 ZNT(I)=MIN(ZNT(I),2.85e-3)
1013                 ZNT(I)=MAX(ZNT(I),1.27e-7)
1014              ENDIF
1015           ENDIF
1016           ZL = ZNT(I)
1017         ELSE
1018           ZL = 0.01
1019         ENDIF
1020         FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/DENOMQ(I)
1021 !       FLQC(I)=RHOX(I)*MAVAIL(I)*UST(I)*KARMAN/(   &
1022 !               ALOG(KARMAN*UST(I)*ZA(I)/XKA+ZA(I)/ZL)-PSIH(I))
1023         DTTHX=ABS(THX(I)-THGB(I))
1024         IF(DTTHX.GT.1.E-5)THEN
1025           FLHC(I)=CPM(I)*RHOX(I)*UST(I)*MOL(I)/(THX(I)-THGB(I))
1026 !         write(*,1001)FLHC(I),CPM(I),RHOX(I),UST(I),MOL(I),THX(I),THGB(I),I
1027  1001   format(f8.5,2x,f12.7,2x,f12.10,2x,f12.10,2x,f13.10,2x,f12.8,f12.8,2x,i3)
1028         ELSE
1029           FLHC(I)=0.
1030         ENDIF
1031   360 CONTINUE
1032
1033 !
1034 !-----COMPUTE SURFACE MOIST FLUX:
1035 !
1036 !     IF(IDRY.EQ.1)GOTO 390
1037 !
1038      IF ( PRESENT(SCM_FORCE_FLUX) ) THEN
1039         IF (SCM_FORCE_FLUX.EQ.1) GOTO 405
1040      ENDIF
1041
1042       DO 370 I=its,ite
1043         QFX(I)=FLQC(I)*(QSFC(I)-QX(I))
1044         QFX(I)=AMAX1(QFX(I),0.)
1045         LH(I)=XLV*QFX(I)
1046   370 CONTINUE
1047
1048 !-----COMPUTE SURFACE HEAT FLUX:
1049 !
1050   390 CONTINUE
1051       DO 400 I=its,ite
1052         IF(XLAND(I)-1.5.GT.0.)THEN
1053           HFX(I)=FLHC(I)*(THGB(I)-THX(I))
1054 !         IF ( PRESENT(ISFTCFLX) ) THEN
1055 !            IF ( ISFTCFLX.NE.0 ) THEN
1056 ! AHW: add dissipative heating term (commented out in 3.6.1)
1057 !               HFX(I)=HFX(I)+RHOX(I)*USTM(I)*USTM(I)*WSPDI(I)
1058 !            ENDIF
1059 !         ENDIF
1060         ELSEIF(XLAND(I)-1.5.LT.0.)THEN
1061           HFX(I)=FLHC(I)*(THGB(I)-THX(I))
1062           HFX(I)=AMAX1(HFX(I),-250.)
1063         ENDIF
1064   400 CONTINUE
1065
1066   405 CONTINUE
1067
1068       DO I=its,ite
1069          IF((XLAND(I)-1.5).GE.0)THEN
1070            ZL=ZNT(I)
1071          ELSE
1072            ZL=0.01
1073          ENDIF
1074 !v3.1.1
1075 !         CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) &
1076 !                /XKA+ZA(I)/ZL)-PSIH(I))
1077          CHS(I)=UST(I)*KARMAN/DENOMQ(I)
1078 !        GZ2OZ0(I)=ALOG(2./ZNT(I))
1079 !        PSIM2(I)=-10.*GZ2OZ0(I)
1080 !        PSIM2(I)=AMAX1(PSIM2(I),-10.)
1081 !        PSIH2(I)=PSIM2(I)
1082 ! v3.1.1
1083 !         CQS2(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*2.0  &
1084 !               /XKA+2.0/ZL)-PSIH2(I))
1085 !         CHS2(I)=UST(I)*KARMAN/(GZ2OZ0(I)-PSIH2(I))
1086          CQS2(I)=UST(I)*KARMAN/DENOMQ2(I)
1087          CHS2(I)=UST(I)*KARMAN/DENOMT2(I)
1088       ENDDO
1089
1090   410 CONTINUE
1091 !jdf
1092 !     DO I=its,ite
1093 !       IF(UST(I).GE.0.1) THEN
1094 !         RMOL(I)=RMOL(I)*(-FLHC(I))/(UST(I)*UST(I)*UST(I))
1095 !       ELSE
1096 !         RMOL(I)=RMOL(I)*(-FLHC(I))/(0.1*0.1*0.1)
1097 !       ENDIF
1098 !     ENDDO
1099 !jdf
1100
1101 !
1102    END SUBROUTINE SFCLAYREV1D
1103
1104 !====================================================================
1105    SUBROUTINE sfclayrevinit(ims,ime,jms,jme,                    &
1106                             its,ite,jts,jte,                    &
1107                             bathymetry_flag, shalwater_z0,      &
1108                             shalwater_depth, water_depth,       &
1109                             xland,LakeModel,lake_depth,lakemask )
1110
1111     INTEGER                   ::      N
1112     REAL                      ::      zolf
1113
1114    INTEGER, INTENT(IN)      ::   ims,ime,jms,jme,its,ite,jts,jte
1115    INTEGER, INTENT(IN)      ::   shalwater_z0
1116    REAL,    INTENT(IN)      ::   shalwater_depth
1117    INTEGER, INTENT(IN)      ::   bathymetry_flag
1118    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)    ::  water_depth
1119    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::  xland
1120    INTEGER         ::     LakeModel
1121    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::  lake_depth
1122    REAL,    DIMENSION( ims:ime, jms:jme )                   ::  lakemask
1123
1124     DO N=0,1000
1125 ! stable function tables
1126        zolf = float(n)*0.01
1127        psim_stab(n)=psim_stable_full(zolf)
1128        psih_stab(n)=psih_stable_full(zolf)
1129
1130 ! unstable function tables
1131        zolf = -float(n)*0.01
1132        psim_unstab(n)=psim_unstable_full(zolf)
1133        psih_unstab(n)=psih_unstable_full(zolf)
1134
1135     ENDDO
1136     IF ( shalwater_z0 .EQ. 1 ) THEN
1137        CALL shalwater_init(ims,ime,jms,jme,            &
1138                  its,ite,jts,jte,                      &
1139                  bathymetry_flag, shalwater_z0,        &
1140                  shalwater_depth, water_depth,         &
1141                  xland,LakeModel,lake_depth,lakemask   )
1142     END IF
1143
1144    END SUBROUTINE sfclayrevinit
1145
1146    SUBROUTINE shalwater_init(ims,ime,jms,jme,                    &
1147                              its,ite,jts,jte,                    &
1148                              bathymetry_flag, shalwater_z0,      &
1149                              shalwater_depth, water_depth,       &
1150                              xland,LakeModel,lake_depth,lakemask )
1151
1152    INTEGER, INTENT(IN)      ::   ims,ime,jms,jme,its,ite,jts,jte
1153    INTEGER, INTENT(IN)      ::   shalwater_z0
1154    REAL,    INTENT(IN)      ::   shalwater_depth
1155    INTEGER, INTENT(IN)      ::   bathymetry_flag
1156    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(INOUT)    ::  water_depth
1157    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::  xland
1158    INTEGER         ::     LakeModel
1159    REAL,    DIMENSION( ims:ime, jms:jme ), INTENT(IN   )    ::  lake_depth
1160    REAL,    DIMENSION( ims:ime, jms:jme )                   ::  lakemask
1161
1162    ! Local
1163    LOGICAL :: overwrite_water_depth
1164
1165    overwrite_water_depth = .False.
1166
1167    IF ( bathymetry_flag .eq. 1 ) THEN
1168       IF ( shalwater_depth .LE. 0.0 ) THEN
1169          IF ( LakeModel .ge. 1 ) THEN
1170             DO j = jts,jte
1171                DO i = its,ite
1172                   IF ( lakemask(i,j) .EQ. 1 ) THEN
1173                      water_depth(i,j) = lake_depth(i,j)
1174                   END IF
1175                END DO
1176             END DO
1177          END IF
1178       ELSE
1179          overwrite_water_depth = .True.
1180       END IF
1181    ELSE
1182       IF ( shalwater_depth .GT. 0.0 ) THEN
1183          overwrite_water_depth = .True.
1184       ELSE
1185          CALL wrf_error_fatal('No bathymetry data detected and shalwater_depth not greater than 0.0. Re-run WPS to get bathymetry data or set shalwater_depth > 0.0')
1186       END IF
1187    END IF
1188
1189    IF (overwrite_water_depth) THEN
1190       DO j = jts,jte
1191          DO i = its,ite
1192             IF((XLAND(i,j)-1.5).GE.0)THEN
1193                water_depth(i,j) = shalwater_depth
1194             ELSE
1195                water_depth(i,j) = -2.0
1196             END IF
1197          END DO
1198       END DO
1199    END IF
1200
1201    END SUBROUTINE shalwater_init
1202
1203       function zolri(ri,z,z0)
1204 !
1205       if (ri.lt.0.)then
1206         x1=-5.
1207         x2=0.
1208       else
1209         x1=0.
1210         x2=5.
1211       endif
1212 !
1213       fx1=zolri2(x1,ri,z,z0)
1214       fx2=zolri2(x2,ri,z,z0)
1215       Do While (abs(x1 - x2) > 0.01)
1216 ! check added for potential divide by zero (2019/11)
1217       if(fx1.eq.fx2)return
1218       if(abs(fx2).lt.abs(fx1))then
1219         x1=x1-fx1/(fx2-fx1)*(x2-x1)
1220         fx1=zolri2(x1,ri,z,z0)
1221         zolri=x1
1222       else
1223         x2=x2-fx2/(fx2-fx1)*(x2-x1)
1224         fx2=zolri2(x2,ri,z,z0)
1225         zolri=x2
1226       endif
1227 !
1228       enddo
1229 !
1230
1231       return
1232       end function
1233
1234 !
1235 ! -----------------------------------------------------------------------
1236 !
1237       function zolri2(zol2,ri2,z,z0)
1238 !
1239       if(zol2*ri2 .lt. 0.)zol2=0.  ! limit zol2 - must be same sign as ri2
1240 !
1241       zol20=zol2*z0/z ! z0/L
1242       zol3=zol2+zol20 ! (z+z0)/L
1243 !
1244       if (ri2.lt.0) then
1245       psix2=log((z+z0)/z0)-(psim_unstable(zol3)-psim_unstable(zol20))
1246       psih2=log((z+z0)/z0)-(psih_unstable(zol3)-psih_unstable(zol20))
1247       else
1248       psix2=log((z+z0)/z0)-(psim_stable(zol3)-psim_stable(zol20))
1249       psih2=log((z+z0)/z0)-(psih_stable(zol3)-psih_stable(zol20))
1250       endif
1251 !
1252       zolri2=zol2*psih2/psix2**2-ri2
1253 !
1254       return
1255       end function
1256 !
1257 ! ... integrated similarity functions ...
1258 !
1259       function psim_stable_full(zolf)
1260         psim_stable_full=-6.1*log(zolf+(1+zolf**2.5)**(1./2.5))
1261       return
1262       end function
1263
1264       function psih_stable_full(zolf)
1265         psih_stable_full=-5.3*log(zolf+(1+zolf**1.1)**(1./1.1))
1266       return
1267       end function
1268
1269       function psim_unstable_full(zolf)
1270         x=(1.-16.*zolf)**.25
1271         psimk=2*ALOG(0.5*(1+X))+ALOG(0.5*(1+X*X))-2.*ATAN(X)+2.*ATAN(1.)
1272 !
1273         ym=(1.-10.*zolf)**0.33
1274         psimc=(3./2.)*log((ym**2.+ym+1.)/3.)-sqrt(3.)*ATAN((2.*ym+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.)
1275 !
1276         psim_unstable_full=(psimk+zolf**2*(psimc))/(1+zolf**2.)
1277
1278       return
1279       end function
1280
1281       function psih_unstable_full(zolf)
1282         y=(1.-16.*zolf)**.5
1283         psihk=2.*log((1+y)/2.)
1284 !
1285         yh=(1.-34.*zolf)**0.33
1286         psihc=(3./2.)*log((yh**2.+yh+1.)/3.)-sqrt(3.)*ATAN((2.*yh+1)/sqrt(3.))+4.*ATAN(1.)/sqrt(3.)
1287 !
1288         psih_unstable_full=(psihk+zolf**2*(psihc))/(1+zolf**2.)
1289
1290       return
1291       end function
1292
1293 ! look-up table functions
1294       function psim_stable(zolf)
1295       integer :: nzol
1296       real    :: rzol
1297         nzol = int(zolf*100.)
1298         rzol = zolf*100. - nzol
1299         if(nzol+1 .lt. 1000)then
1300            psim_stable = psim_stab(nzol) + rzol*(psim_stab(nzol+1)-psim_stab(nzol))
1301         else
1302            psim_stable = psim_stable_full(zolf)
1303         endif
1304       return
1305       end function
1306
1307       function psih_stable(zolf)
1308       integer :: nzol
1309       real    :: rzol
1310         nzol = int(zolf*100.)
1311         rzol = zolf*100. - nzol
1312         if(nzol+1 .lt. 1000)then
1313            psih_stable = psih_stab(nzol) + rzol*(psih_stab(nzol+1)-psih_stab(nzol))
1314         else
1315            psih_stable = psih_stable_full(zolf)
1316         endif
1317       return
1318       end function
1319
1320       function psim_unstable(zolf)
1321       integer :: nzol
1322       real    :: rzol
1323         nzol = int(-zolf*100.)
1324         rzol = -zolf*100. - nzol
1325         if(nzol+1 .lt. 1000)then
1326            psim_unstable = psim_unstab(nzol) + rzol*(psim_unstab(nzol+1)-psim_unstab(nzol))
1327         else
1328            psim_unstable = psim_unstable_full(zolf)
1329         endif
1330       return
1331       end function
1332
1333       function psih_unstable(zolf)
1334       integer :: nzol
1335       real    :: rzol
1336         nzol = int(-zolf*100.)
1337         rzol = -zolf*100. - nzol
1338         if(nzol+1 .lt. 1000)then
1339            psih_unstable = psih_unstab(nzol) + rzol*(psih_unstab(nzol+1)-psih_unstab(nzol))
1340         else
1341            psih_unstable = psih_unstable_full(zolf)
1342         endif
1343       return
1344       end function
1345
1346       function depth_dependent_z0(water_depth,z0,UST)
1347       real :: depth_b
1348       real :: effective_depth
1349          IF ( water_depth .lt. 10.0 ) THEN
1350            effective_depth = 10.0
1351          ELSEIF ( water_depth .gt. 100.0 ) THEN
1352            effective_depth = 100.0
1353          ELSE
1354            effective_depth = water_depth
1355          ENDIF
1356
1357          depth_b = 1 / 30.0 * log (1260.0 / effective_depth)
1358          depth_dependent_z0 = exp((2.7 * ust - 1.8 / depth_b) / (ust + 0.17 / depth_b) )
1359          depth_dependent_z0 = MIN(depth_dependent_z0,0.1)
1360      return
1361      end function
1362 !-------------------------------------------------------------------
1363
1364 END MODULE module_sf_sfclayrev
1365
1366 !
1367 ! ----------------------------------------------------------
1368 !
1369
1370