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