phys/module_bl_gfs.F

   1 !LWRF:MODEL_LAYER:PHYSICS
   2 !
   3 MODULE module_bl_gfs
   4
   5 CONTAINS
   6
   7 !-------------------------------------------------------------------
   8    SUBROUTINE BL_GFS(U3D,V3D,TH3D,T3D,QV3D,QC3D,QI3D,P3D,PI3D,     &
   9                   RUBLTEN,RVBLTEN,RTHBLTEN,                        &
  10                   RQVBLTEN,RQCBLTEN,RQIBLTEN,                      &
  11                   CP,G,ROVCP,R,ROVG,P_QI,P_FIRST_SCALAR,           &
  12                   dz8w,z,PSFC,                                     &
  13                   UST,PBL,PSIM,PSIH,                               &
  14                   HFX,QFX,TSK,GZ1OZ0,WSPD,BR,                      &
  15                   DT,KPBL2D,EP1,KARMAN,                            &
  16 #if (NMM_CORE==1)
  17                   DISHEAT,                                         &
  18 #endif
  19 #if (HWRF==1)
  20                   ALPHA,                                           &
  21                   HPBL2D, EVAP2D, HEAT2D,                          &    !Kwon add FOR SHAL. CON.
  22                   VAR_RIC,                                         &    !Kwon for variable Ric
  23                   U10,V10,ZNT,MZNT,rc2d,                           &    !Kwon for variable Ric
  24                   DKU3D,DKT3D,coef_ric_l,coef_ric_s,xland,         &    !Kwon for variable Ric
  25                   msang,scurx,scury,iwavecpl,lcurr_sf,                      &
  26                   pert_pbl, ens_random_seed, ens_pblamp,           &
  27 #endif
  28                   ids,ide, jds,jde, kds,kde,                       &
  29                   ims,ime, jms,jme, kms,kme,                       &
  30                   its,ite, jts,jte, kts,kte                        )
  31 !--------------------------------------------------------------------
  32       USE MODULE_GFS_MACHINE , ONLY : kind_phys
  33 !-------------------------------------------------------------------
  34       IMPLICIT NONE
  35 !-------------------------------------------------------------------
  36 !-- U3D         3D u-velocity interpolated to theta points (m/s)
  37 !-- V3D         3D v-velocity interpolated to theta points (m/s)
  38 !-- TH3D        3D potential temperature (K)
  39 !-- T3D         temperature (K)
  40 !-- QV3D        3D water vapor mixing ratio (Kg/Kg)
  41 !-- QC3D        3D cloud mixing ratio (Kg/Kg)
  42 !-- QI3D        3D ice mixing ratio (Kg/Kg)
  43 !-- P3D         3D pressure (Pa)
  44 !-- PI3D        3D exner function (dimensionless)
  45 !-- rr3D        3D dry air density (kg/m^3)
  46 !-- RUBLTEN     U tendency due to
  47 !               PBL parameterization (m/s^2)
  48 !-- RVBLTEN     V tendency due to
  49 !               PBL parameterization (m/s^2)
  50 !-- RTHBLTEN    Theta tendency due to
  51 !               PBL parameterization (K/s)
  52 !-- RQVBLTEN    Qv tendency due to
  53 !               PBL parameterization (kg/kg/s)
  54 !-- RQCBLTEN    Qc tendency due to
  55 !               PBL parameterization (kg/kg/s)
  56 !-- RQIBLTEN    Qi tendency due to
  57 !               PBL parameterization (kg/kg/s)
  58 !-- CP          heat capacity at constant pressure for dry air (J/kg/K)
  59 !-- G           acceleration due to gravity (m/s^2)
  60 !-- ROVCP       R/CP
  61 !-- R           gas constant for dry air (J/kg/K)
  62 !-- ROVG        R/G
  63 !-- P_QI        species index for cloud ice
  64 !-- dz8w        dz between full levels (m)
  65 !-- z           height above sea level (m)
  66 !-- PSFC        pressure at the surface (Pa)
  67 !-- UST         u* in similarity theory (m/s)
  68 !-- PBL         PBL height (m)
  69 !-- PSIM        similarity stability function for momentum
  70 !-- PSIH        similarity stability function for heat
  71 !-- HFX         upward heat flux at the surface (W/m^2)
  72 !-- QFX         upward moisture flux at the surface (kg/m^2/s)
  73 !-- TSK         surface temperature (K)
  74 !-- GZ1OZ0      log(z/z0) where z0 is roughness length
  75 !-- WSPD        wind speed at lowest model level (m/s)
  76 !-- BR          bulk Richardson number in surface layer
  77 !-- DT          time step (s)
  78 !-- rvovrd      R_v divided by R_d (dimensionless)
  79 !-- EP1         constant for virtual temperature (R_v/R_d - 1) (dimensionless)
  80 !-- KARMAN      Von Karman constant
  81 !-- ALPHA       boundary depth scaling factor
  82 !-- VAR_RIC     Flag for using variable Ric or not (=1: variable Ric, =0: constant Ric)
  83 !-- RO          Surface Rossby number
  84 !-- ids         start index for i in domain
  85 !-- ide         end index for i in domain
  86 !-- jds         start index for j in domain
  87 !-- jde         end index for j in domain
  88 !-- kds         start index for k in domain
  89 !-- kde         end index for k in domain
  90 !-- ims         start index for i in memory
  91 !-- ime         end index for i in memory
  92 !-- jms         start index for j in memory
  93 !-- jme         end index for j in memory
  94 !-- kms         start index for k in memory
  95 !-- kme         end index for k in memory
  96 !-- its         start index for i in tile
  97 !-- ite         end index for i in tile
  98 !-- jts         start index for j in tile
  99 !-- jte         end index for j in tile
 100 !-- kts         start index for k in tile
 101 !-- kte         end index for k in tile
 102 !-------------------------------------------------------------------
 103
 104 #if (NMM_CORE==1)
 105       LOGICAL , INTENT(IN)::            DISHEAT                         !  gopal's doing
 106 #endif
 107
 108 #if (HWRF==1)
 109       INTEGER , INTENT(IN)          ::  iwavecpl
 110       LOGICAL , INTENT(IN)          ::  lcurr_sf
 111       REAL,  DIMENSION(ims:ime, jms:jme), INTENT(OUT) ::              &
 112                                         HPBL2D,                         &    !ADDED BY KWON FOR SHALLOW CONV.
 113                                         EVAP2D,                         &    !ADDED BY KWON FOR SHALLOW CONV.
 114                                         HEAT2D,RC2D,MZNT                     !ADDED BY KWON FOR SHALLOW CONV.
 115       REAL,  DIMENSION(ims:ime, jms:jme), INTENT(IN) ::              &
 116                                         U10,                         &    !ADDED BY KWON FOR VARIABLE Ric
 117                                         V10,XLAND,                   &    !ADDED BY KWON FOR VARIABLE Ric
 118                                         ZNT                               !ADDED BY KWON FOR VARIABLE Ric
 119       REAL,  DIMENSION(ims:ime, jms:jme, kms:kme), INTENT(OUT) :: DKU3D,DKT3D
 120       REAL,    INTENT(IN) :: VAR_RIC,coef_ric_l,coef_ric_s                   !ADDED BY KWON
 121       REAL,  DIMENSION(ims:ime, jms:jme), INTENT(IN) ::                 &
 122                                         SCURX,                          &
 123                                         SCURY,                          &
 124                                         MSANG
 125       integer,intent(in) :: ens_random_seed
 126       real,intent(in) :: ens_pblamp
 127       logical,intent(in) :: pert_pbl
 128 #endif
 129
 130
 131
 132       INTEGER, INTENT(IN) ::            ids,ide, jds,jde, kds,kde,      &
 133                                         ims,ime, jms,jme, kms,kme,      &
 134                                         its,ite, jts,jte, kts,kte,      &
 135                                         P_QI,P_FIRST_SCALAR
 136
 137       REAL,    INTENT(IN) ::                                            &
 138                                         CP,                             &
 139                                         DT,                             &
 140                                         EP1,                            &
 141                                         G,                              &
 142                                         KARMAN,                         &
 143                                         R,                              &
 144                                         ROVCP,                          &
 145                                         ROVG
 146
 147       REAL,    DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) ::      &
 148                                         DZ8W,                           &
 149                                         P3D,                            &
 150                                         PI3D,                           &
 151                                         QC3D,                           &
 152                                         QI3D,                           &
 153                                         QV3D,                           &
 154                                         T3D,                            &
 155                                         TH3D,                           &
 156                                         U3D,                            &
 157                                         V3D,                            &
 158                                         Z
 159
 160
 161       REAL,    DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::   &
 162                                         RTHBLTEN,                       &
 163                                         RQCBLTEN,                       &
 164                                         RQIBLTEN,                       &
 165                                         RQVBLTEN,                       &
 166                                         RUBLTEN,                        &
 167                                         RVBLTEN
 168
 169       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(IN) ::               &
 170                                         BR,                             &
 171                                         GZ1OZ0,                         &
 172                                         HFX,                            &
 173                                         PSFC,                           &
 174                                         PSIM,                           &
 175                                         PSIH,                           &
 176                                         QFX,                            &
 177                                         TSK
 178
 179       REAL,    DIMENSION(ims:ime, jms:jme), INTENT(INOUT) ::            &
 180                                         PBL,                            &
 181                                         UST,                            &
 182                                         WSPD
 183
 184       INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(OUT) ::              &
 185                                         KPBL2D
 186
 187
 188 !--------------------------- LOCAL VARS ------------------------------
 189
 190
 191       REAL     (kind=kind_phys), DIMENSION(its:ite, kts:kte) ::         &
 192                                         DEL,                            &
 193                                         DU,                             &
 194                                         DV,                             &
 195                                         PHIL,                           &
 196                                         PRSL,                           &
 197                                         PRSLK,                          &
 198                                         T1,                             &
 199                                         TAU,                            &
 200                                         dishx,                          &
 201 #if (HWRF==1)
 202                                         dku,dkt,   &       !Kwon for diffusivity
 203 #endif
 204                                         U1,                             &
 205                                         V1
 206
 207       REAL     (kind=kind_phys), DIMENSION(its:ite, kts:kte+1) ::       &
 208                                         PHII,                           &
 209                                         PRSI
 210
 211       REAL     (kind=kind_phys), DIMENSION(its:ite, kts:kte, 3) ::      &
 212                                         Q1,                             &
 213                                         RTG
 214
 215       REAL     (kind=kind_phys), DIMENSION(its:ite) ::                  &
 216                                         DQSFC,                          &
 217                                         DTSFC,                          &
 218                                         DUSFC,                          &
 219                                         DVSFC,                          &
 220                                         EVAP,                           &
 221                                         FH,                             &
 222                                         FM,                             &
 223                                         HEAT,                           &
 224                                         HGAMQ,                          &
 225                                         HGAMT,                          &
 226                                         HPBL,                           &
 227                                         PSK,                            &
 228                                         QSS,                            &
 229                                         RBSOIL,                         &
 230                                         RCL,                            &
 231                                         XLAND1,                         &
 232                                         SPD1,                           &
 233                                         STRESS,                         &
 234                                         RO,rbcr,                        &  !Kwon for variablr Ric(surface Rossby number)
 235                                         TSEA
 236
 237       REAL     (kind=kind_phys) ::                                      &
 238                                         CPM,                            &
 239                                         cpmikj,                         &
 240                                         DELTIM,                         &
 241                                         FMTMP,                          &
 242                                         RRHOX
 243
 244 #if (HWRF == 1)
 245       REAL      ::      ALPHA
 246 #else
 247       REAL, PARAMETER:: ALPHA=1.0
 248 #endif
 249
 250 #if (HWRF == 1)
 251       REAL :: UBOT, VBOT, UBOT1, VBOT1
 252 #endif
 253       INTEGER, DIMENSION( its:ite ) ::                                  &
 254                                         KPBL
 255
 256       INTEGER ::                                                        &
 257                                         I,                              &
 258                                         IM,                             &
 259                                         J,                              &
 260                                         K,                              &
 261                                         KM,                             &
 262                                         KTEM,                           &
 263                                         KTEP,                           &
 264                                         KX,                             &
 265                                         L,                              &
 266                                         NTRAC
 267
 268    IM=ITE-ITS+1
 269    KX=KTE-KTS+1
 270    KTEM=KTE-1
 271    KTEP=KTE+1
 272    NTRAC=2
 273    DELTIM=DT
 274    IF (P_QI.ge.P_FIRST_SCALAR) NTRAC=3
 275
 276
 277    DO J=jts,jte
 278
 279       DO i=its,ite
 280         RRHOX=(R*T3D(I,KTS,J)*(1.+EP1*QV3D(I,KTS,J)))/PSFC(I,J)
 281         CPM=CP*(1.+0.8*QV3D(i,kts,j))
 282         FMTMP=GZ1OZ0(i,j)-PSIM(i,j)
 283         PSK(i)=(PSFC(i,j)*.00001)**ROVCP
 284         FM(i)=FMTMP
 285         FH(i)=GZ1OZ0(i,j)-PSIH(i,j)
 286         TSEA(i)=TSK(i,j)
 287         QSS(i)=QV3D(i,kts,j)               ! not used in moninp so set to qv3d for now
 288         HEAT(i)=HFX(i,j)/CPM*RRHOX
 289         EVAP(i)=QFX(i,j)*RRHOX
 290         XLAND1(i) = 0.0
 291 #if (HWRF==1)
 292 ! Kwon FOR NEW SHALLOW CONVECTION
 293         HEAT2D(i,j)=HFX(i,j)/CPM*RRHOX
 294         EVAP2D(i,j)=QFX(i,j)*RRHOX
 295         XLAND1(i) = XLAND(I,J)
 296 #endif
 297 !
 298 #if (HWRF==1)
 299         IF ( XLAND(I,J)  > 1.99 ) THEN
 300           IF ( LCURR_SF ) THEN
 301             UBOT = U3D(I,KTS,J)-SCURX(I,J)
 302             VBOT = V3D(I,KTS,J)-SCURY(I,J)
 303           ELSE
 304             UBOT = U3D(I,KTS,J)
 305             VBOT = V3D(I,KTS,J)
 306           ENDIF
 307           IF ( IWAVECPL .eq. 1 ) THEN
 308             UBOT1 = ( UBOT * COS(MSANG(I,J))  -               &
 309                      VBOT * SIN(MSANG(I,J)) )                 &
 310                                   * COS(MSANG(I,J))
 311             VBOT1 = ( VBOT * COS(MSANG(I,J))  -               &
 312                      UBOT * SIN(MSANG(I,J)) )                &
 313                                   * COS(MSANG(I,J))
 314
 315             WSPD(i,j) = SQRT(UBOT1*UBOT1+VBOT1*VBOT1) + 1.E-9
 316           ELSE
 317             WSPD(i,j) = SQRT(UBOT*UBOT+VBOT*VBOT) + 1.E-9
 318           ENDIF
 319         ENDIF
 320 #endif
 321         STRESS(i)=KARMAN*KARMAN*WSPD(i,j)*WSPD(i,j)/(FMTMP*FMTMP)
 322         SPD1(i)=WSPD(i,j)
 323         PRSI(i,kts)=PSFC(i,j)*.001
 324         PHII(I,kts)=0.
 325         RCL(i)=1.
 326         RBSOIL(I)=BR(i,j)
 327 #if (HWRF==1)
 328 ! Kwon for variable Ric   : Ro=W10/(f*zo): surface Rossby number
 329        Ro(I)=SQRT(U10(I,J)**2 + V10(I,J)**2) / (1.E-4 * MZNT(I,J))
 330 #endif
 331       ENDDO
 332
 333       DO k=kts,kte
 334         DO i=its,ite
 335           DV(I,K) = 0.
 336           DU(I,K) = 0.
 337           TAU(I,K) = 0.
 338 #if (HWRF==1)
 339           IF ( XLAND(I,J) > 1.99  .AND. k == KTS ) THEN
 340             IF ( LCURR_SF ) THEN
 341               UBOT = U3D(i,k,j) - SCURX(I,J)
 342               VBOT = V3D(i,k,j) - SCURY(I,J)
 343             ELSE
 344               UBOT = U3D(i,k,j)
 345               VBOT = V3D(i,k,j)
 346             ENDIF
 347             IF ( IWAVECPL .eq. 1 ) THEN
 348               U1(I,K) = ( UBOT * COS(MSANG(I,J))  -            &
 349                           VBOT * SIN(MSANG(I,J)) )             &
 350                                   * COS(MSANG(I,J))
 351               V1(I,K) = ( VBOT * COS(MSANG(I,J))  -            &
 352                           UBOT * SIN(MSANG(I,J)) )             &
 353                                   * COS(MSANG(I,J))
 354             ELSE
 355               U1(I,K) =  UBOT
 356               V1(I,K) =  VBOT
 357             ENDIF
 358           ELSE
 359             U1(I,K) = U3D(i,k,j)
 360             V1(I,K) = V3D(i,k,j)
 361           ENDIF
 362 #else
 363           U1(I,K) = U3D(i,k,j)
 364           V1(I,K) = V3D(i,k,j)
 365 #endif
 366           T1(I,K) = T3D(i,k,j)
 367           Q1(I,K,1) = QV3D(i,k,j)/(1.+QV3D(i,k,j))
 368           Q1(I,K,2) = QC3D(i,k,j)/(1.+QC3D(i,k,j))
 369           PRSL(I,K)=P3D(i,k,j)*.001
 370         ENDDO
 371       ENDDO
 372
 373       DO k=kts,kte
 374         DO i=its,ite
 375           PRSLK(I,K)=(PRSL(i,k)*.01)**ROVCP
 376         ENDDO
 377       ENDDO
 378
 379
 380       DO k=kts+1,kte
 381         km=k-1
 382         DO i=its,ite
 383           DEL(i,km)=PRSL(i,km)/ROVG*dz8w(i,km,j)/T3D(i,km,j)
 384           PRSI(i,k)=PRSI(i,km)-DEL(i,km)
 385           PHII(I,K)=(Z(i,k,j)-Z(i,kts,j))*G
 386           PHIL(I,KM)=0.5*(Z(i,k,j)+Z(i,km,j)-2.*Z(i,kts,j))*G
 387         ENDDO
 388       ENDDO
 389
 390       DO i=its,ite
 391         DEL(i,kte)=DEL(i,ktem)
 392         PRSI(i,ktep)=PRSI(i,kte)-DEL(i,ktem)
 393         PHII(I,KTEP)=PHII(I,KTE)+dz8w(i,kte,j)*G
 394         PHIL(I,KTE)=PHII(I,KTE)-PHIL(I,KTEM)+PHII(I,KTE)
 395       ENDDO
 396
 397       IF (P_QI.ge.P_FIRST_SCALAR) THEN
 398         DO k=kts,kte
 399           DO i=its,ite
 400             Q1(I,K,3) = QI3D(i,k,j)/(1.+QI3D(i,k,j))
 401           ENDDO
 402         ENDDO
 403       ENDIF
 404
 405       DO l=1,ntrac
 406         DO k=kts,kte
 407           DO i=its,ite
 408             RTG(I,K,L) = 0.
 409           ENDDO
 410         ENDDO
 411       ENDDO
 412 !
 413       CALL MONINP(IM,IM,KX,NTRAC,DV,DU,TAU,RTG,U1,V1,T1,Q1,             &
 414                   PSK,RBSOIL,FM,FH,TSEA,QSS,HEAT,EVAP,STRESS,           &
 415                   SPD1,KPBL,PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL,          &
 416                   DELTIM,DUSFC,DVSFC,DTSFC,DQSFC,HPBL,HGAMT,            &
 417 #if (HWRF==1)
 418                VAR_RIC,Ro,DKU,DKT,coef_ric_l,coef_ric_s,xland1,  &
 419               pert_pbl, ens_random_seed, ens_pblamp,             &
 420 #endif
 421                 RBCR,HGAMQ,ALPHA)
 422
 423 !============================================================================
 424 !    ADD  IN  DISSIPATIVE HEATING .... v*dv. This is Bob's doing
 425 !============================================================================
 426
 427 #if (NMM_CORE==1)
 428
 429       IF(DISHEAT)THEN
 430        DO k=kts,kte
 431          DO i=its,ite
 432           dishx(i,k)=u1(i,k)*du(i,k) + v1(i,k)*dv(i,k)
 433           cpmikj=CP*(1.+0.8*QV3D(i,k,j))
 434           dishx(i,k)=-dishx(i,k)/cpmikj
 435 !         IF(k==1)WRITE(0,*)'ADDITIONAL DISSIPATIVE HEATING',tau(i,k),dishx(i,k)
 436           tau(i,k)=tau(i,k)+dishx(i,k)
 437          ENDDO
 438        ENDDO
 439       ENDIF
 440 #endif
 441
 442 !=============================================================================
 443
 444
 445       DO k=kts,kte
 446         DO i=its,ite
 447           RVBLTEN(I,K,J)=DV(I,K)
 448           RUBLTEN(I,K,J)=DU(I,K)
 449           RTHBLTEN(I,K,J)=TAU(I,K)/PI3D(I,K,J)
 450           RQVBLTEN(I,K,J)=RTG(I,K,1)/(1.-Q1(I,K,1))**2
 451           RQCBLTEN(I,K,J)=RTG(I,K,2)/(1.-Q1(I,K,2))**2
 452         ENDDO
 453       ENDDO
 454
 455       IF (P_QI.ge.P_FIRST_SCALAR) THEN
 456         DO k=kts,kte
 457           DO i=its,ite
 458             RQIBLTEN(I,K,J)=RTG(I,K,3)/(1.-Q1(I,K,3))**2
 459           ENDDO
 460         ENDDO
 461       ENDIF
 462
 463       DO i=its,ite
 464         UST(i,j)=SQRT(STRESS(i))
 465         WSPD(i,j)=SQRT(U3D(I,KTS,J)*U3D(I,KTS,J)+                       &
 466                        V3D(I,KTS,J)*V3D(I,KTS,J))+1.E-9
 467         PBL(i,j)=HPBL(i)
 468 #if (HWRF==1)
 469 !Kwon For new shallow convection
 470         HPBL2D(i,j)=HPBL(i)
 471         rc2D(i,j)=rbcr(i)
 472 #endif
 473 !
 474         KPBL2D(i,j)=kpbl(i)
 475       ENDDO
 476 ! INITIALIZE DKU3D and DKT3D  (3D momentum and thermal diffusivity for
 477 ! diagnostics)
 478 !
 479 #if (HWRF==1)
 480      DO i=its,ite
 481      DO k=kts,kte
 482       DKU3D(I,J,K) = 0.
 483       DKT3D(I,J,K) = 0.
 484      ENDDO
 485      ENDDO
 486
 487      DO i=its,ite
 488      DO k=kts,kte-1
 489       DKU3D(I,J,K) = DKU(I,K)
 490       DKT3D(I,J,K) = DKT(I,K)
 491      ENDDO
 492      ENDDO
 493 #endif
 494
 495     ENDDO
 496
 497
 498    END SUBROUTINE BL_GFS
 499
 500 !===================================================================
 501    SUBROUTINE gfsinit(RUBLTEN,RVBLTEN,RTHBLTEN,RQVBLTEN,           &
 502                       RQCBLTEN,RQIBLTEN,P_QI,P_FIRST_SCALAR,       &
 503                       restart,                                     &
 504                       allowed_to_read,                             &
 505                       ids, ide, jds, jde, kds, kde,                &
 506                       ims, ime, jms, jme, kms, kme,                &
 507                       its, ite, jts, jte, kts, kte                 )
 508 !-------------------------------------------------------------------
 509    IMPLICIT NONE
 510 !-------------------------------------------------------------------
 511    LOGICAL , INTENT(IN)          ::  allowed_to_read,restart
 512    INTEGER , INTENT(IN)          ::  ids, ide, jds, jde, kds, kde, &
 513                                      ims, ime, jms, jme, kms, kme, &
 514                                      its, ite, jts, jte, kts, kte
 515    INTEGER , INTENT(IN)          ::  P_QI,P_FIRST_SCALAR
 516
 517    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) ::         &
 518                                                          RUBLTEN, &
 519                                                          RVBLTEN, &
 520                                                          RTHBLTEN, &
 521                                                          RQVBLTEN, &
 522                                                          RQCBLTEN, &
 523                                                          RQIBLTEN
 524    INTEGER :: i, j, k, itf, jtf, ktf
 525
 526    jtf=min0(jte,jde-1)
 527    ktf=min0(kte,kde-1)
 528    itf=min0(ite,ide-1)
 529
 530    IF(.not.restart)THEN
 531      DO j=jts,jtf
 532      DO k=kts,ktf
 533      DO i=its,itf
 534         RUBLTEN(i,k,j)=0.
 535         RVBLTEN(i,k,j)=0.
 536         RTHBLTEN(i,k,j)=0.
 537         RQVBLTEN(i,k,j)=0.
 538         RQCBLTEN(i,k,j)=0.
 539      ENDDO
 540      ENDDO
 541      ENDDO
 542    ENDIF
 543
 544    IF (P_QI .ge. P_FIRST_SCALAR .and. .not.restart) THEN
 545       DO j=jts,jtf
 546       DO k=kts,ktf
 547       DO i=its,itf
 548          RQIBLTEN(i,k,j)=0.
 549       ENDDO
 550       ENDDO
 551       ENDDO
 552    ENDIF
 553
 554    IF (P_QI .ge. P_FIRST_SCALAR) THEN
 555       DO j=jts,jtf
 556       DO k=kts,ktf
 557       DO i=its,itf
 558          RQIBLTEN(i,k,j)=0.
 559       ENDDO
 560       ENDDO
 561       ENDDO
 562    ENDIF
 563
 564
 565    END SUBROUTINE gfsinit
 566
 567 ! --------------------------------------------------------------
 568 !FPP$ NOCONCUR R
 569       SUBROUTINE MONINP(IX,IM,KM,ntrac,DV,DU,TAU,RTG,                   &
 570      &     U1,V1,T1,Q1,                                                 &
 571      &     PSK,RBSOIL,FM,FH,TSEA,QSS,HEAT,EVAP,STRESS,SPD1,KPBL,        &
 572      &     PRSI,DEL,PRSL,PRSLK,PHII,PHIL,RCL,DELTIM,                    &
 573      &     DUSFC,DVSFC,DTSFC,DQSFC,HPBL,HGAMT,                     &
 574 #if (HWRF==1)
 575            VAR_RIC,Ro,DKU,DKT,coef_ric_l,coef_ric_s,xland1,   &
 576            pert_pbl, ens_random_seed, ens_pblamp,             &
 577 #endif
 578            RBCR,HGAMQ,ALPHA)
 579 !
 580       USE MODULE_GFS_MACHINE, ONLY : kind_phys
 581       USE MODULE_GFS_PHYSCONS, grav => con_g, RD => con_RD, CP => con_CP &
 582      &,             HVAP => con_HVAP, ROG => con_ROG, FV => con_FVirt
 583
 584       implicit none
 585 !
 586 !     include 'constant.h'
 587 !
 588 !
 589 !     Arguments
 590 !
 591       integer IX, IM, KM, ntrac, KPBL(IM)
 592 !
 593       real(kind=kind_phys) DELTIM
 594       real :: ALPHA
 595
 596 #if (HWRF==1)
 597       real :: VAR_RIC,coef_ric_l,coef_ric_s
 598       integer,intent(in) :: ens_random_seed
 599       real,intent(in) :: ens_pblamp
 600       logical,intent(in) :: pert_pbl
 601 #endif
 602       real(kind=kind_phys) DV(IM,KM),     DU(IM,KM),                    &
 603      &                     TAU(IM,KM),    RTG(IM,KM,ntrac),             &
 604      &                     U1(IX,KM),     V1(IX,KM),                    &
 605      &                     T1(IX,KM),     Q1(IX,KM,ntrac),              &
 606      &                     PSK(IM),       RBSOIL(IM),                   &
 607 !    &                     CD(IM),        CH(IM),                       &
 608      &                     FM(IM),        FH(IM),                       &
 609      &                     TSEA(IM),      QSS(IM),                      &
 610      &                                    SPD1(IM),                     &
 611 !    &                     DPHI(IM),      SPD1(IM),                     &
 612      &                     PRSI(IX,KM+1), DEL(IX,KM),                   &
 613      &                     PRSL(IX,KM),   PRSLK(IX,KM),                 &
 614      &                     PHII(IX,KM+1), PHIL(IX,KM),                  &
 615      &                     RCL(IM),       DUSFC(IM),                    &
 616      &                     dvsfc(IM),     dtsfc(IM),                    &
 617      &                     DQSFC(IM),     HPBL(IM),                     &
 618      &                     HGAMT(IM),     hgamq(IM), RBCR(IM)
 619 #if (HWRF==1)
 620 real(kind=kind_phys) RO(IM),xland1(IM)
 621 #endif
 622 !
 623 !    Locals
 624 !
 625       integer i,iprt,is,iun,k,kk,kmpbl,lond
 626 !     real(kind=kind_phys) betaq(IM), betat(IM),   betaw(IM),           &
 627       real(kind=kind_phys) evap(IM),  heat(IM),    phih(IM),            &
 628      &                     phim(IM),  rbdn(IM),    rbup(IM),            &
 629      &                     the1(IM),  stress(im),  beta(im),            &
 630      &                     the1v(IM), thekv(IM),   thermal(IM),         &
 631      &                     thesv(IM), ustar(IM),   wscale(IM)
 632 !    &                     thesv(IM), ustar(IM),   wscale(IM),  zl1(IM)
 633 !
 634       real(kind=kind_phys) RDZT(IM,KM-1),                               &
 635      &                     ZI(IM,KM+1),     ZL(IM,KM),                  &
 636      &                     DKO(IM,KM-1),                                &
 637      &                     AL(IM,KM-1),     AD(IM,KM),                  &
 638      &                     AU(IM,KM-1),     A1(IM,KM),                  &
 639      &                     A2(IM,KM),       THETA(IM,KM),               &
 640      &                     AT(IM,KM*(ntrac-1)),DKU(IM,KM-1),DKT(IM,KM-1),WSPM(IM,KM-1) ! RGF added WSPM
 641       logical              pblflg(IM),   sfcflg(IM), stable(IM)
 642 !
 643       real(kind=kind_phys) aphi16,  aphi5,  bet1,   bvf2,               &
 644      &                     cfac,    conq,   cont,   conw,               &
 645      &                     conwrc,  dk,     dkmax,  dkmin,              &
 646      &                     dq1,     dsdz2,  dsdzq,  dsdzt,              &
 647      &                     dsig,    dt,     dthe1,  dtodsd,             &
 648      &                     dtodsu,  dw2,    dw2min, g,                  &
 649      &                     gamcrq,  gamcrt, gocp,   gor, gravi,         &
 650      &                     hol,     pfac,   prmax,  prmin, prinv,       &
 651      &                     prnum,   qmin,   qtend,                      &
 652      &                     rbint,   rdt,    rdz,    rdzt1,              &
 653      &                     ri,      rimin,  rl2,    rlam,               &
 654      &                     rone,   rzero,  sfcfrac,                     &
 655      &                     sflux,   shr2,   spdk2,  sri,                &
 656      &                     tem,     ti,     ttend,  tvd,                &
 657      &                     tvu,     utend,  vk,     vk2,                &
 658      &                     vpert,   vtend,  xkzo,   zfac,               &
 659      &                     zfmin,   zk,     tem1
 660       integer kLOC ! RGF
 661       real xDKU    ! RGF
 662 !
 663       PARAMETER(g=grav)
 664       PARAMETER(GOR=G/RD,GOCP=G/CP)
 665       PARAMETER(CONT=1000.*CP/G,CONQ=1000.*HVAP/G,CONW=1000./G)
 666       PARAMETER(RLAM=150.,VK=0.4,VK2=VK*VK,PRMIN=1.0,PRMAX=4.)
 667       PARAMETER(DW2MIN=0.0001,DKMIN=1.0,DKMAX=1000.,RIMIN=-100.)
 668       PARAMETER(CFAC=7.8,PFAC=2.0,SFCFRAC=0.1)
 669       PARAMETER(QMIN=1.E-8,XKZO=1.0,ZFMIN=1.E-8,APHI5=5.,APHI16=16.)
 670 !     PARAMETER(GAMCRT=3.,GAMCRQ=2.E-3)
 671       PARAMETER(GAMCRT=3.,GAMCRQ=0.)
 672       PARAMETER(RZERO=0.,RONE=1.)
 673       PARAMETER(IUN=84)
 674 #if HWRF==1
 675       real*8 :: ran1           !zhang
 676       real :: rr
 677       logical,save :: pert_pbl_local !zhang
 678       integer,save :: ens_random_seed_local !zhang
 679       real,save :: ens_pblamp_local   !zhang
 680       data ens_random_seed_local/0/
 681 !zz      print*, 'zhang in pbl==========='
 682       if ( ens_random_seed_local .eq. 0 ) then
 683          pert_pbl_local=pert_pbl
 684          ens_random_seed_local=ens_random_seed
 685          ens_pblamp_local=ens_pblamp
 686 !zz      print*, "zhang in pbl= one time ", pert_pbl_local, ens_random_seed_local, ens_pblamp_local
 687       endif
 688 !zz      print*, "zhang in pbl=",pert_pbl_local, ens_random_seed_local, ens_pblamp_local
 689 #endif
 690 !
 691 !
 692 !-----------------------------------------------------------------------
 693 !
 694  601  FORMAT(1X,' MONINP LAT LON STEP HOUR ',3I6,F6.1)
 695  602      FORMAT(1X,'    K','        Z','        T','       TH',        &
 696      &     '      TVH','        Q','        U','        V',             &
 697      &     '       SP')
 698  603      FORMAT(1X,I5,8F9.1)
 699  604      FORMAT(1X,'  SFC',9X,F9.1,18X,F9.1)
 700  605      FORMAT(1X,'    K      ZL    SPD2   THEKV   THE1V'             &
 701      &         ,' THERMAL    RBUP')
 702  606      FORMAT(1X,I5,6F8.2)
 703  607      FORMAT(1X,' KPBL    HPBL      FM      FH   HGAMT',            &
 704      &         '   HGAMQ      WS   USTAR      CD      CH')
 705  608      FORMAT(1X,I5,9F8.2)
 706  609      FORMAT(1X,' K PR DKT DKU ',I5,3F8.2)
 707  610      FORMAT(1X,' K PR DKT DKU ',I5,3F8.2,' L2 RI T2',              &
 708      &         ' SR2  ',2F8.2,2E10.2)
 709 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 710 !     COMPUTE PRELIMINARY VARIABLES
 711 !
 712
 713       if (IX .lt. im) stop
 714 !
 715       IPRT = 0
 716       IF(IPRT.EQ.1) THEN
 717 !!!   LATD = 0
 718       LOND = 0
 719       ELSE
 720 !!!   LATD = 0
 721       LOND = 0
 722       ENDIF
 723 !
 724 ! define critical Richardson number    by KWON: Vickers and Mahart(2004) J. Appl. Meteo.
 725 !  coef_ric=0.16 originally but it may too small: controled by namelist=0.25
 726 !  Land and Ocean points are treated differently
 727 !  by Kwon
 728 !
 729      do i=1,im
 730      RBCR(I) = 0.25
 731 #if (HWRF==1)
 732      IF(var_ric.eq.1.) THEN
 733       IF(xland1(i).eq.1)  RBCR(I) = coef_ric_l*(1.E-7*Ro(I))**(-0.18)
 734       IF(xland1(i).eq.2)  RBCR(I) = coef_ric_s*(1.E-7*Ro(I))**(-0.18)
 735 !      write(0,*) 'xland1 coef_ric_l coef_ric_s ',xland1(i),coef_ric_l,coef_ric_s
 736      ENDIF
 737      IF(RBCR(I).GT.0.5) RBCR(I)=0.5    !set upper limit Suggsted by Han
 738 #endif
 739      enddo
 740 !
 741       gravi = 1.0 / grav
 742       DT    = 2. * DELTIM
 743       RDT   = 1. / DT
 744       KMPBL = KM / 2
 745 !
 746       do k=1,km
 747         do i=1,im
 748           zi(i,k) = phii(i,k) * gravi
 749           zl(i,k) = phil(i,k) * gravi
 750         enddo
 751       enddo
 752 !
 753       do k=1,kmpbl
 754         do i=1,im
 755           theta(i,k) = t1(i,k) * psk(i) / prslk(i,k)
 756         enddo
 757       enddo
 758 !
 759       DO K = 1,KM-1
 760         DO I=1,IM
 761           RDZT(I,K) = GOR * PRSI(I,K+1) / (PRSL(I,K) - PRSL(I,K+1))
 762         ENDDO
 763       ENDDO
 764 !
 765       DO I = 1,IM
 766          DUSFC(I) = 0.
 767          DVSFC(I) = 0.
 768          DTSFC(I) = 0.
 769          DQSFC(I) = 0.
 770          HGAMT(I) = 0.
 771          HGAMQ(I) = 0.
 772          WSCALE(I) = 0.
 773          KPBL(I) = 1
 774          HPBL(I) = ZI(I,2)
 775          PBLFLG(I) = .TRUE.
 776          SFCFLG(I) = .TRUE.
 777          IF(RBSOIL(I).GT.0.0) SFCFLG(I) = .FALSE.
 778       ENDDO
 779 !!
 780       DO I=1,IM
 781          RDZT1    = GOR * prSL(i,1) / DEL(i,1)
 782 !        BET1     = DT*RDZT1*SPD1(I)/T1(I,1)
 783          BETA(I)  = DT*RDZT1/T1(I,1)
 784 !        BETAW(I) = BET1*CD(I)
 785 !        BETAT(I) = BET1*CH(I)
 786 !        BETAQ(I) = DPHI(I)*BETAT(I)
 787       ENDDO
 788 !
 789       DO I=1,IM
 790 !        ZL1(i) = 0.-(T1(I,1)+TSEA(I))/2.*LOG(PRSL(I,1)/PRSI(I,1))*ROG
 791 !        USTAR(I) = SQRT(CD(I)*SPD1(I)**2)
 792          USTAR(I) = SQRT(STRESS(I))
 793       ENDDO
 794 !
 795       DO I=1,IM
 796          THESV(I)   = TSEA(I)*(1.+FV*MAX(QSS(I),QMIN))
 797          THE1(I)    = THETA(I,1)
 798          THE1V(I)   = THE1(I)*(1.+FV*MAX(Q1(I,1,1),QMIN))
 799          THERMAL(I) = THE1V(I)
 800 !        DTHE1      = (THE1(I)-TSEA(I))
 801 !        DQ1        = (MAX(Q1(I,1,1),QMIN) - MAX(QSS(I),QMIN))
 802 !        HEAT(I)    = -CH(I)*SPD1(I)*DTHE1
 803 !        EVAP(I)    = -CH(I)*SPD1(I)*DQ1
 804       ENDDO
 805 !
 806 !
 807 !     COMPUTE THE FIRST GUESS OF PBL HEIGHT
 808 !
 809       DO I=1,IM
 810          STABLE(I) = .FALSE.
 811 !        ZL(i,1) = ZL1(i)
 812          RBUP(I) = RBSOIL(I)
 813       ENDDO
 814       DO K = 2, KMPBL
 815         DO I = 1, IM
 816           IF(.NOT.STABLE(I)) THEN
 817              RBDN(I)   = RBUP(I)
 818 !            ZL(I,k)   = ZL(I,K-1) - (T1(i,k)+T1(i,K-1))/2 *
 819 !    &                   LOG(PRSL(I,K)/PRSL(I,K-1)) * ROG
 820              THEKV(I)  = THETA(i,k)*(1.+FV*MAX(Q1(i,k,1),QMIN))
 821              SPDK2 = MAX(RCL(i)*(U1(i,k)*U1(i,k)+V1(i,k)*V1(i,k)),RONE)
 822              RBUP(I)   = (THEKV(I)-THE1V(I))*(G*ZL(I,k)/THE1V(I))/SPDK2
 823              KPBL(I)   = K
 824              STABLE(I) = RBUP(I).GT.RBCR(I)
 825           ENDIF
 826         ENDDO
 827       ENDDO
 828 !
 829       DO I = 1,IM
 830          K = KPBL(I)
 831          IF(RBDN(I).GE.RBCR(I)) THEN
 832             RBINT = 0.
 833          ELSEIF(RBUP(I).LE.RBCR(I)) THEN
 834             RBINT = 1.
 835          ELSE
 836             RBINT = (RBCR(I)-RBDN(I))/(RBUP(I)-RBDN(I))
 837          ENDIF
 838          HPBL(I) = ZL(I,K-1) + RBINT*(ZL(I,K)-ZL(I,K-1))
 839          IF(HPBL(I).LT.ZI(I,KPBL(I))) KPBL(I) = KPBL(I) - 1
 840       ENDDO
 841 !!
 842       DO I=1,IM
 843            HOL = MAX(RBSOIL(I)*FM(I)*FM(I)/FH(I),RIMIN)
 844            IF(SFCFLG(I)) THEN
 845               HOL = MIN(HOL,-ZFMIN)
 846            ELSE
 847               HOL = MAX(HOL,ZFMIN)
 848            ENDIF
 849 !
 850 !          HOL = HOL*HPBL(I)/ZL1(I)*SFCFRAC
 851            HOL = HOL*HPBL(I)/ZL(I,1)*SFCFRAC
 852            IF(SFCFLG(I)) THEN
 853 !             PHIM = (1.-APHI16*HOL)**(-1./4.)
 854 !             PHIH = (1.-APHI16*HOL)**(-1./2.)
 855               TEM  = 1.0 / (1. - APHI16*HOL)
 856               PHIH(I) = SQRT(TEM)
 857               PHIM(I) = SQRT(PHIH(I))
 858            ELSE
 859               PHIM(I) = (1.+APHI5*HOL)
 860               PHIH(I) = PHIM(I)
 861            ENDIF
 862            WSCALE(I) = USTAR(I)/PHIM(I)
 863            WSCALE(I) = MIN(WSCALE(I),USTAR(I)*APHI16)
 864            WSCALE(I) = MAX(WSCALE(I),USTAR(I)/APHI5)
 865       ENDDO
 866 !
 867 !     COMPUTE THE SURFACE VARIABLES FOR PBL HEIGHT ESTIMATION
 868 !     UNDER UNSTABLE CONDITIONS
 869 !
 870       DO I = 1,IM
 871          SFLUX  = HEAT(I) + EVAP(I)*FV*THE1(I)
 872          IF(SFCFLG(I).AND.SFLUX.GT.0.0) THEN
 873            HGAMT(I)   = MIN(CFAC*HEAT(I)/WSCALE(I),GAMCRT)
 874            HGAMQ(I)   = MIN(CFAC*EVAP(I)/WSCALE(I),GAMCRQ)
 875            VPERT      = HGAMT(I) + FV*THE1(I)*HGAMQ(I)
 876            VPERT      = MIN(VPERT,GAMCRT)
 877            THERMAL(I) = THERMAL(I) + MAX(VPERT,RZERO)
 878            HGAMT(I)   = MAX(HGAMT(I),RZERO)
 879            HGAMQ(I)   = MAX(HGAMQ(I),RZERO)
 880          ELSE
 881            PBLFLG(I) = .FALSE.
 882          ENDIF
 883       ENDDO
 884 !
 885       DO I = 1,IM
 886          IF(PBLFLG(I)) THEN
 887             KPBL(I) = 1
 888             HPBL(I) = ZI(I,2)
 889          ENDIF
 890       ENDDO
 891 !
 892 !     ENHANCE THE PBL HEIGHT BY CONSIDERING THE THERMAL
 893 !
 894       DO I = 1, IM
 895          IF(PBLFLG(I)) THEN
 896             STABLE(I) = .FALSE.
 897             RBUP(I) = RBSOIL(I)
 898          ENDIF
 899       ENDDO
 900       DO K = 2, KMPBL
 901         DO I = 1, IM
 902           IF(.NOT.STABLE(I).AND.PBLFLG(I)) THEN
 903             RBDN(I)   = RBUP(I)
 904 !           ZL(I,k)   = ZL(I,K-1) - (T1(i,k)+T1(i,K-1))/2 *
 905 !    &                  LOG(PRSL(I,K)/PRSL(I,K-1)) * ROG
 906             THEKV(I)  = THETA(i,k)*(1.+FV*MAX(Q1(i,k,1),QMIN))
 907             SPDK2   = MAX(RCL(i)*(U1(i,k)*U1(i,k)+V1(i,k)*V1(i,k)),RONE)
 908             RBUP(I)   = (THEKV(I)-THERMAL(I))*(G*ZL(I,k)/THE1V(I))/SPDK2
 909             KPBL(I)   = K
 910             STABLE(I) = RBUP(I).GT.RBCR(I)
 911           ENDIF
 912         ENDDO
 913       ENDDO
 914 !
 915       DO I = 1,IM
 916          IF(PBLFLG(I)) THEN
 917             K = KPBL(I)
 918             IF(RBDN(I).GE.RBCR(I)) THEN
 919                RBINT = 0.
 920             ELSEIF(RBUP(I).LE.RBCR(I)) THEN
 921                RBINT = 1.
 922             ELSE
 923                RBINT = (RBCR(I)-RBDN(I))/(RBUP(I)-RBDN(I))
 924             ENDIF
 925             HPBL(I) = ZL(I,K-1) + RBINT*(ZL(I,k)-ZL(I,K-1))
 926 #if (HWRF==1)
 927 !zhang adding PBL perturtion
 928             if ( pert_pbl_local ) then
 929             rr=(2.0*ens_pblamp_local*ran1(ens_random_seed_local)-ens_pblamp_local)
 930             print*, "zhang inside the loop", rr, ens_pblamp_local,ens_random_seed_local
 931             HPBL(I) = HPBL(I)*(1.0+rr)
 932             endif
 933 #endif
 934             IF(HPBL(I).LT.ZI(I,KPBL(I))) KPBL(I) = KPBL(I) - 1
 935             IF(KPBL(I).LE.1) PBLFLG(I) = .FALSE.
 936          ENDIF
 937       ENDDO
 938 !!
 939 !
 940 !     COMPUTE DIFFUSION COEFFICIENTS BELOW PBL
 941 !
 942
 943 #if (HWRF==1)
 944 ! -------------------------------------------------------------------------------------
 945 ! begin RGF modifications
 946 ! this is version MOD05
 947
 948
 949 ! RGF determine wspd at roughly 500 m above surface, or as close as possible, reuse SPDK2
 950 !  zi(i,k) is AGL, right?  May not matter if applied only to water grid points
 951       if(ALPHA.lt.0)then
 952
 953        DO I=1,IM
 954          SPDK2 = 0.
 955          WSPM(i,1) = 0.
 956          DO K = 1, KMPBL ! kmpbl is like a max possible pbl height
 957           if(zi(i,k).le.500.and.zi(i,k+1).gt.500.)then ! find level bracketing 500 m
 958            SPDK2 = SQRT(U1(i,k)*U1(i,k)+V1(i,k)*V1(i,k)) ! wspd near 500 m
 959            WSPM(i,1) = SPDK2/0.6  ! now the Km limit for 500 m.  just store in K=1
 960            WSPM(i,2) = float(k)  ! height of level at gridpoint i. store in K=2
 961 !           if(i.eq.25) print *,' IK ',i,k,' ZI ',zi(i,k), ' WSPM1 ',wspm(i,1),' KMPBL ',kmpbl,' KPBL ',kpbl(i)
 962           endif
 963          ENDDO
 964        ENDDO ! i
 965
 966       endif ! ALPHA < 0
 967 #endif
 968
 969
 970       DO I=1,IM                 ! RGF SWAPPED ORDER.  TESTED - NO IMPACT
 971          DO K = 1, KMPBL
 972
 973
 974
 975 ! First guess at DKU.  If alpha >= 0, this is the only loop executed
 976
 977             IF(KPBL(I).GT.K) THEN ! first guess DKU, this is original loop
 978                PRINV = 1.0 / (PHIH(I)/PHIM(I)+CFAC*VK*.1)
 979                PRINV = MIN(PRINV,PRMAX)
 980                PRINV = MAX(PRINV,PRMIN)
 981 !              ZFAC = MAX((1.-(ZI(I,K+1)-ZL1(I))/                       &
 982 !    &                (HPBL(I)-ZL1(I))), ZFMIN)
 983                ZFAC = MAX((1.-(ZI(I,K+1)-ZL(I,1))/                      &
 984      &                (HPBL(I)-ZL(I,1))), ZFMIN)
 985 !
 986 !              alpha factor (0-1.0) is multiplied to profile function to reduce the effect
 987 !              height of the Hurricane Boundary Layer. This is gopal's doing
 988 !
 989
 990
 991                DKU(i,k) = XKZO + WSCALE(I)*VK*ZI(I,K+1)                 &
 992      &                         *ABS(ALPHA)* ZFAC**PFAC
 993 ! if alpha = -1, the above provides the first guess for DKU, based on assumption alpha = +1
 994 !               (other values of alpha < 0 can also be applied)
 995 ! if alpha > 0, the above applies the alpha suppression factor and we are finished
 996 !               if(i.eq.25) print *,' I25 K ',k,' ORIG DKU ',dku(i,k)
 997                DKT(i,k) = DKU(i,k)*PRINV
 998                DKO(i,k) = (DKU(i,k)-XKZO)*PRINV
 999                DKU(i,k) = MIN(DKU(i,k),DKMAX)
1000                DKU(i,k) = MAX(DKU(i,k),DKMIN)
1001                DKT(i,k) = MIN(DKT(i,k),DKMAX)
1002                DKT(i,k) = MAX(DKT(i,k),DKMIN)
1003                DKO(i,k) = MAX(RZERO, MIN(DKMAX, DKO(i,k)))
1004               endif ! KPBL
1005         ENDDO ! K
1006
1007
1008 #if (HWRF==1)
1009 ! possible modification of first guess DKU, under certain conditions
1010 ! (1) this applies only to columns over water
1011
1012         IF(xland1(i).eq.2)then ! sea only
1013
1014 ! (2) alpha test
1015 ! if alpha < 0, find alpha for each column and do the loop again
1016 ! if alpha > 0, we are finished
1017
1018
1019         if(alpha.lt.0)then      ! variable alpha test
1020
1021 ! k-level of layer around 500 m
1022             kLOC = INT(WSPM(i,2))
1023 !            print *,' kLOC ',kLOC,' KPBL ',KPBL(I)
1024
1025 ! (3) only do  this IF KPBL(I) >= kLOC.  Otherwise, we are finished, with DKU as if alpha = +1
1026
1027           if(KPBL(I).gt.kLOC)then
1028
1029             xDKU = DKU(i,kLOC)     ! Km at k-level
1030
1031 ! (4) DKU check.
1032 ! WSPM(i,1) is the KM cap for the 500-m level.
1033 !  if DKU at 500-m level < WSPM(i,1), do not limit Km ANYWHERE.  Alpha = abs(alpha).  No need to recalc.
1034 !  if DKU at 500-m level > WSPM(i,1), then alpha = WSPM(i,1)/xDKU for entire column
1035             if(xDKU.ge.WSPM(i,1)) then ! ONLY if DKU at 500-m exceeds cap, otherwise already done
1036
1037
1038             WSPM(i,3) = WSPM(i,1)/xDKU  ! ratio of cap to Km at k-level, store in WSPM(i,3)
1039             WSPM(i,4) = amin1(WSPM(I,3),1.0) ! this is new column alpha. cap at 1. ! should never be needed
1040
1041 !    if(i.eq.25) print *,' I25 kLOC ',kLOC,' xDKU ',xDKU,' WSPM1 ',WSPM(i,1),' WSPM3 ',WSPM(i,3),' WSPM4 ',WSPM(i,4)
1042
1043             DO K = 1, KMPBL ! now go through K loop again
1044              IF(KPBL(I).GT.K) THEN
1045                PRINV = 1.0 / (PHIH(I)/PHIM(I)+CFAC*VK*.1)
1046                PRINV = MIN(PRINV,PRMAX)
1047                PRINV = MAX(PRINV,PRMIN)
1048 !              ZFAC = MAX((1.-(ZI(I,K+1)-ZL1(I))/                       &
1049 !    &                (HPBL(I)-ZL1(I))), ZFMIN)
1050                ZFAC = MAX((1.-(ZI(I,K+1)-ZL(I,1))/                      &
1051      &                (HPBL(I)-ZL(I,1))), ZFMIN)
1052
1053                DKU(i,k) = XKZO + WSCALE(I)*VK*ZI(I,K+1)                 &
1054      &                         *WSPM(i,4)* ZFAC**PFAC ! recalculated DKU using column alpha
1055 !               if(i.eq.25) print *,' I25 K ',k,' DKU AFTER ',dku(i,k)
1056
1057
1058                DKT(i,k) = DKU(i,k)*PRINV
1059                DKO(i,k) = (DKU(i,k)-XKZO)*PRINV
1060                DKU(i,k) = MIN(DKU(i,k),DKMAX)
1061                DKU(i,k) = MAX(DKU(i,k),DKMIN)
1062                DKT(i,k) = MIN(DKT(i,k),DKMAX)
1063                DKT(i,k) = MAX(DKT(i,k),DKMIN)
1064                DKO(i,k) = MAX(RZERO, MIN(DKMAX, DKO(i,k)))
1065             ENDIF ! KPBL
1066           ENDDO ! DO K (RGF ALTERED)
1067          endif ! xDKU.ge.WSPM(i,1)
1068          endif ! KPBL(I).ge.kLOC
1069         endif ! alpha < 0
1070         endif ! xland1 = 2
1071 #endif
1072       ENDDO ! DO I
1073
1074 ! end RGF modifications
1075 ! -------------------------------------------------------------------------------------
1076
1077 !
1078 !     COMPUTE DIFFUSION COEFFICIENTS OVER PBL (FREE ATMOSPHERE)
1079 !
1080       DO K = 1, KM-1
1081          DO I=1,IM
1082             IF(K.GE.KPBL(I)) THEN
1083 !              TI   = 0.5*(T1(i,k)+T1(i,K+1))
1084                TI   = 2.0 / (T1(i,k)+T1(i,K+1))
1085 !              RDZ  = RDZT(I,K)/TI
1086                RDZ  = RDZT(I,K) * TI
1087 !              RDZ  = RDZT(I,K)
1088                DW2  = RCL(i)*((U1(i,k)-U1(i,K+1))**2                    &
1089      &                      + (V1(i,k)-V1(i,K+1))**2)
1090                SHR2 = MAX(DW2,DW2MIN)*RDZ**2
1091                TVD  = T1(i,k)*(1.+FV*MAX(Q1(i,k,1),QMIN))
1092                TVU  = T1(i,K+1)*(1.+FV*MAX(Q1(i,K+1,1),QMIN))
1093 !              BVF2 = G*(GOCP+RDZ*(TVU-TVD))/TI
1094                BVF2 = G*(GOCP+RDZ*(TVU-TVD)) * TI
1095                RI   = MAX(BVF2/SHR2,RIMIN)
1096                ZK   = VK*ZI(I,K+1)
1097 !              RL2  = (ZK*RLAM/(RLAM+ZK))**2
1098 !              DK   = RL2*SQRT(SHR2)
1099                RL2  = ZK*RLAM/(RLAM+ZK)
1100                DK   = RL2*RL2*SQRT(SHR2)
1101                IF(RI.LT.0.) THEN ! UNSTABLE REGIME
1102                   SRI = SQRT(-RI)
1103                   DKU(i,k) = XKZO + DK*(1+8.*(-RI)/(1+1.746*SRI))
1104 !                 DKT(i,k) = XKZO + DK*(1+8.*(-RI)/(1+1.286*SRI))
1105                   tem      =        DK*(1+8.*(-RI)/(1+1.286*SRI))
1106                   DKT(i,k) = XKZO + tem
1107                   DKO(i,k) =        tem
1108                ELSE             ! STABLE REGIME
1109 !                 DKT(i,k)  = XKZO + DK/(1+5.*RI)**2
1110                   tem       =        DK/(1+5.*RI)**2
1111                   DKT(i,k)  = XKZO + tem
1112                   DKO(i,k)  =        tem
1113                   PRNUM     = 1.0 + 2.1*RI
1114                   PRNUM     = MIN(PRNUM,PRMAX)
1115                   DKU(i,k)  = (DKT(i,k)-XKZO)*PRNUM + XKZO
1116                ENDIF
1117 !
1118                DKU(i,k) = MIN(DKU(i,k),DKMAX)
1119                DKU(i,k) = MAX(DKU(i,k),DKMIN)
1120                DKT(i,k) = MIN(DKT(i,k),DKMAX)
1121                DKT(i,k) = MAX(DKT(i,k),DKMIN)
1122                DKO(i,k) = MAX(RZERO, MIN(DKMAX, DKO(i,k)))
1123 !
1124 !!!   IF(I.EQ.LOND.AND.LAT.EQ.LATD) THEN
1125 !!!   PRNUM = DKU(k)/DKT(k)
1126 !!!   WRITE(IUN,610) K,PRNUM,DKT(k),DKU(k),RL2,RI,
1127 !!!   1              BVF2,SHR2
1128 !!!   ENDIF
1129 !
1130             ENDIF
1131          ENDDO
1132       ENDDO
1133
1134 !
1135 !     COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR HEAT AND MOISTURE
1136 !
1137       DO I=1,IM
1138          AD(I,1) = 1.
1139          A1(I,1) = T1(i,1)   + BETA(i) * HEAT(I)
1140          A2(I,1) = Q1(i,1,1) + BETA(i) * EVAP(I)
1141 !        A1(I,1) = T1(i,1)-BETAT(I)*(THETA(i,1)-TSEA(I))
1142 !        A2(I,1) = Q1(i,1,1)-BETAQ(I)*
1143 !    &           (MAX(Q1(i,1,1),QMIN)-MAX(QSS(I),QMIN))
1144       ENDDO
1145 !
1146       DO K = 1,KM-1
1147         DO I = 1,IM
1148           DTODSD = DT/DEL(I,K)
1149           DTODSU = DT/DEL(I,K+1)
1150           DSIG   = PRSL(I,K)-PRSL(I,K+1)
1151           RDZ    = RDZT(I,K)*2./(T1(i,k)+T1(i,K+1))
1152 !         RDZ    = RDZT(I,K)
1153           tem1   = DSIG * DKT(i,k) * RDZ
1154           IF(PBLFLG(I).AND.K.LT.KPBL(I)) THEN
1155 !            DSDZT = DSIG*DKT(i,k)*RDZ*(GOCP-HGAMT(I)/HPBL(I))
1156 !            DSDZQ = DSIG*DKT(i,k)*RDZ*(-HGAMQ(I)/HPBL(I))
1157              tem   = 1.0 / HPBL(I)
1158              DSDZT = tem1 * (GOCP-HGAMT(I)*tem)
1159              DSDZQ = tem1 * (-HGAMQ(I)*tem)
1160              A2(I,k)   = A2(I,k)+DTODSD*DSDZQ
1161              A2(I,k+1) = Q1(i,k+1,1)-DTODSU*DSDZQ
1162           ELSE
1163 !            DSDZT = DSIG*DKT(i,k)*RDZ*(GOCP)
1164              DSDZT = tem1 * GOCP
1165              A2(I,k+1) = Q1(i,k+1,1)
1166           ENDIF
1167 !         DSDZ2 = DSIG*DKT(i,k)*RDZ*RDZ
1168           DSDZ2     = tem1 * RDZ
1169           AU(I,k)   = -DTODSD*DSDZ2
1170           AL(I,k)   = -DTODSU*DSDZ2
1171           AD(I,k)   = AD(I,k)-AU(I,k)
1172           AD(I,k+1) = 1.-AL(I,k)
1173           A1(I,k)   = A1(I,k)+DTODSD*DSDZT
1174           A1(I,k+1) = T1(i,k+1)-DTODSU*DSDZT
1175         ENDDO
1176       ENDDO
1177 !
1178 !     SOLVE TRIDIAGONAL PROBLEM FOR HEAT AND MOISTURE
1179 !
1180       CALL TRIDIN(IM,KM,1,AL,AD,AU,A1,A2,AU,A1,A2)
1181 !
1182 !     RECOVER TENDENCIES OF HEAT AND MOISTURE
1183 !
1184       DO  K = 1,KM
1185          DO I = 1,IM
1186             TTEND      = (A1(I,k)-T1(i,k))*RDT
1187             QTEND      = (A2(I,k)-Q1(i,k,1))*RDT
1188             TAU(i,k)   = TAU(i,k)+TTEND
1189             RTG(I,k,1) = RTG(i,k,1)+QTEND
1190             DTSFC(I)   = DTSFC(I)+CONT*DEL(I,K)*TTEND
1191             DQSFC(I)   = DQSFC(I)+CONQ*DEL(I,K)*QTEND
1192          ENDDO
1193       ENDDO
1194 !
1195 !     COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR MOMENTUM
1196 !
1197       DO I=1,IM
1198 !        AD(I,1) = 1.+BETAW(I)
1199          AD(I,1) = 1.0 + BETA(i) * STRESS(I) / SPD1(I)
1200          A1(I,1) = U1(i,1)
1201          A2(I,1) = V1(i,1)
1202 !        AD(I,1) = 1.0
1203 !        tem     = 1.0 + BETA(I) * STRESS(I) / SPD1(I)
1204 !        A1(I,1) = U1(i,1) * tem
1205 !        A2(I,1) = V1(i,1) * tem
1206       ENDDO
1207 !
1208       DO K = 1,KM-1
1209         DO I=1,IM
1210           DTODSD    = DT/DEL(I,K)
1211           DTODSU    = DT/DEL(I,K+1)
1212           DSIG      = PRSL(I,K)-PRSL(I,K+1)
1213           RDZ       = RDZT(I,K)*2./(T1(i,k)+T1(i,k+1))
1214 !         RDZ       = RDZT(I,K)
1215           DSDZ2     = DSIG*DKU(i,k)*RDZ*RDZ
1216           AU(I,k)   = -DTODSD*DSDZ2
1217           AL(I,k)   = -DTODSU*DSDZ2
1218           AD(I,k)   = AD(I,k)-AU(I,k)
1219           AD(I,k+1) = 1.-AL(I,k)
1220           A1(I,k+1) = U1(i,k+1)
1221           A2(I,k+1) = V1(i,k+1)
1222         ENDDO
1223       ENDDO
1224 !
1225 !     SOLVE TRIDIAGONAL PROBLEM FOR MOMENTUM
1226 !
1227       CALL TRIDI2(IM,KM,AL,AD,AU,A1,A2,AU,A1,A2)
1228 !
1229 !     RECOVER TENDENCIES OF MOMENTUM
1230 !
1231       DO K = 1,KM
1232          DO I = 1,IM
1233             CONWRC = CONW*SQRT(RCL(i))
1234             UTEND = (A1(I,k)-U1(i,k))*RDT
1235             VTEND = (A2(I,k)-V1(i,k))*RDT
1236             DU(i,k)  = DU(i,k)+UTEND
1237             DV(i,k)  = DV(i,k)+VTEND
1238             DUSFC(I) = DUSFC(I)+CONWRC*DEL(I,K)*UTEND
1239             DVSFC(I) = DVSFC(I)+CONWRC*DEL(I,K)*VTEND
1240          ENDDO
1241       ENDDO
1242 !!
1243 !
1244 !     COMPUTE TRIDIAGONAL MATRIX ELEMENTS FOR TRACERS
1245 !
1246       if (ntrac .ge. 2) then
1247         DO I=1,IM
1248          AD(I,1) = 1.
1249         ENDDO
1250         do k = 2, ntrac
1251           is = (k-2) * km
1252           do i = 1, im
1253             AT(I,1+is) = Q1(i,1,k)
1254           enddo
1255         enddo
1256 !
1257         DO K = 1,KM-1
1258           DO I = 1,IM
1259             DTODSD = DT/DEL(I,K)
1260             DTODSU = DT/DEL(I,K+1)
1261             DSIG   = PRSL(I,K)-PRSL(I,K+1)
1262             RDZ    = RDZT(I,K)*2./(T1(i,k)+T1(i,K+1))
1263             tem1   = DSIG * DKT(i,k) * RDZ
1264             DSDZ2     = tem1 * RDZ
1265             AU(I,k)   = -DTODSD*DSDZ2
1266             AL(I,k)   = -DTODSU*DSDZ2
1267             AD(I,k)   = AD(I,k)-AU(I,k)
1268             AD(I,k+1) = 1.-AL(I,k)
1269           ENDDO
1270         ENDDO
1271         do kk = 2, ntrac
1272           is = (kk-2) * km
1273           do k = 1, km - 1
1274             do i = 1, im
1275               AT(I,k+1+is) = Q1(i,k+1,kk)
1276             enddo
1277           enddo
1278         enddo
1279 !
1280 !     SOLVE TRIDIAGONAL PROBLEM FOR TRACERS
1281 !
1282         CALL TRIDIT(IM,KM,ntrac-1,AL,AD,AU,AT,AU,AT)
1283 !
1284 !     RECOVER TENDENCIES OF TRACERS
1285 !
1286         do kk = 2, ntrac
1287           is = (kk-2) * km
1288           do k = 1, km
1289             do i = 1, im
1290               QTEND = (AT(I,K+is)-Q1(i,K,kk))*RDT
1291               RTG(i,K,kk) = RTG(i,K,kk) + QTEND
1292             enddo
1293           enddo
1294         enddo
1295       endif
1296 !!
1297       RETURN
1298       END SUBROUTINE MONINP
1299 !FPP$ NOCONCUR R
1300 !-----------------------------------------------------------------------
1301       SUBROUTINE TRIDI2(L,N,CL,CM,CU,R1,R2,AU,A1,A2)
1302 !sela %INCLUDE DBTRIDI2;
1303 !
1304       USE MODULE_GFS_MACHINE, ONLY : kind_phys
1305       implicit none
1306       integer             k,n,l,i
1307       real(kind=kind_phys) fk
1308 !
1309       real(kind=kind_phys) CL(L,2:N),CM(L,N),CU(L,N-1),R1(L,N),R2(L,N), &
1310      &          AU(L,N-1),A1(L,N),A2(L,N)
1311 !-----------------------------------------------------------------------
1312       DO I=1,L
1313         FK      = 1./CM(I,1)
1314         AU(I,1) = FK*CU(I,1)
1315         A1(I,1) = FK*R1(I,1)
1316         A2(I,1) = FK*R2(I,1)
1317       ENDDO
1318       DO K=2,N-1
1319         DO I=1,L
1320           FK      = 1./(CM(I,K)-CL(I,K)*AU(I,K-1))
1321           AU(I,K) = FK*CU(I,K)
1322           A1(I,K) = FK*(R1(I,K)-CL(I,K)*A1(I,K-1))
1323           A2(I,K) = FK*(R2(I,K)-CL(I,K)*A2(I,K-1))
1324         ENDDO
1325       ENDDO
1326       DO I=1,L
1327         FK      = 1./(CM(I,N)-CL(I,N)*AU(I,N-1))
1328         A1(I,N) = FK*(R1(I,N)-CL(I,N)*A1(I,N-1))
1329         A2(I,N) = FK*(R2(I,N)-CL(I,N)*A2(I,N-1))
1330       ENDDO
1331       DO K=N-1,1,-1
1332         DO I=1,L
1333           A1(I,K) = A1(I,K)-AU(I,K)*A1(I,K+1)
1334           A2(I,K) = A2(I,K)-AU(I,K)*A2(I,K+1)
1335         ENDDO
1336       ENDDO
1337 !-----------------------------------------------------------------------
1338       RETURN
1339       END SUBROUTINE TRIDI2
1340 !FPP$ NOCONCUR R
1341 !-----------------------------------------------------------------------
1342       SUBROUTINE TRIDIN(L,N,nt,CL,CM,CU,R1,R2,AU,A1,A2)
1343 !sela %INCLUDE DBTRIDI2;
1344 !
1345       USE MODULE_GFS_MACHINE, ONLY : kind_phys
1346       implicit none
1347       integer             is,k,kk,n,nt,l,i
1348       real(kind=kind_phys) fk(L)
1349 !
1350       real(kind=kind_phys) CL(L,2:N), CM(L,N), CU(L,N-1),               &
1351      &                     R1(L,N),   R2(L,N*nt),                       &
1352      &                     AU(L,N-1), A1(L,N), A2(L,N*nt),              &
1353      &                     FKK(L,2:N-1)
1354 !-----------------------------------------------------------------------
1355       DO I=1,L
1356         FK(I)   = 1./CM(I,1)
1357         AU(I,1) = FK(I)*CU(I,1)
1358         A1(I,1) = FK(I)*R1(I,1)
1359       ENDDO
1360       do k = 1, nt
1361         is = (k-1) * n
1362         do i = 1, l
1363           a2(i,1+is) = fk(I) * r2(i,1+is)
1364         enddo
1365       enddo
1366       DO K=2,N-1
1367         DO I=1,L
1368           FKK(I,K) = 1./(CM(I,K)-CL(I,K)*AU(I,K-1))
1369           AU(I,K)  = FKK(I,K)*CU(I,K)
1370           A1(I,K)  = FKK(I,K)*(R1(I,K)-CL(I,K)*A1(I,K-1))
1371         ENDDO
1372       ENDDO
1373       do kk = 1, nt
1374         is = (kk-1) * n
1375         DO K=2,N-1
1376           DO I=1,L
1377             A2(I,K+is) = FKK(I,K)*(R2(I,K+is)-CL(I,K)*A2(I,K+is-1))
1378           ENDDO
1379         ENDDO
1380       ENDDO
1381       DO I=1,L
1382         FK(I)   = 1./(CM(I,N)-CL(I,N)*AU(I,N-1))
1383         A1(I,N) = FK(I)*(R1(I,N)-CL(I,N)*A1(I,N-1))
1384       ENDDO
1385       do k = 1, nt
1386         is = (k-1) * n
1387         do i = 1, l
1388           A2(I,N+is) = FK(I)*(R2(I,N+is)-CL(I,N)*A2(I,N+is-1))
1389         enddo
1390       enddo
1391       DO K=N-1,1,-1
1392         DO I=1,L
1393           A1(I,K) = A1(I,K) - AU(I,K)*A1(I,K+1)
1394         ENDDO
1395       ENDDO
1396       do kk = 1, nt
1397         is = (kk-1) * n
1398         DO K=n-1,1,-1
1399           DO I=1,L
1400             A2(I,K+is) = A2(I,K+is) - AU(I,K)*A2(I,K+is+1)
1401           ENDDO
1402         ENDDO
1403       ENDDO
1404 !-----------------------------------------------------------------------
1405       RETURN
1406       END SUBROUTINE TRIDIN
1407       SUBROUTINE TRIDIT(L,N,nt,CL,CM,CU,RT,AU,AT)
1408 !sela %INCLUDE DBTRIDI2;
1409 !
1410       USE MODULE_GFS_MACHINE, ONLY : kind_phys
1411       implicit none
1412       integer             is,k,kk,n,nt,l,i
1413       real(kind=kind_phys) fk(L)
1414 !
1415       real(kind=kind_phys) CL(L,2:N), CM(L,N), CU(L,N-1),               &
1416      &                     RT(L,N*nt),                                  &
1417      &                     AU(L,N-1), AT(L,N*nt),                       &
1418      &                     FKK(L,2:N-1)
1419 !-----------------------------------------------------------------------
1420       DO I=1,L
1421         FK(I)   = 1./CM(I,1)
1422         AU(I,1) = FK(I)*CU(I,1)
1423       ENDDO
1424       do k = 1, nt
1425         is = (k-1) * n
1426         do i = 1, l
1427           at(i,1+is) = fk(I) * rt(i,1+is)
1428         enddo
1429       enddo
1430       DO K=2,N-1
1431         DO I=1,L
1432           FKK(I,K) = 1./(CM(I,K)-CL(I,K)*AU(I,K-1))
1433           AU(I,K)  = FKK(I,K)*CU(I,K)
1434         ENDDO
1435       ENDDO
1436       do kk = 1, nt
1437         is = (kk-1) * n
1438         DO K=2,N-1
1439           DO I=1,L
1440             AT(I,K+is) = FKK(I,K)*(RT(I,K+is)-CL(I,K)*AT(I,K+is-1))
1441           ENDDO
1442         ENDDO
1443       ENDDO
1444       DO I=1,L
1445         FK(I)   = 1./(CM(I,N)-CL(I,N)*AU(I,N-1))
1446       ENDDO
1447       do k = 1, nt
1448         is = (k-1) * n
1449         do i = 1, l
1450           AT(I,N+is) = FK(I)*(RT(I,N+is)-CL(I,N)*AT(I,N+is-1))
1451         enddo
1452       enddo
1453       do kk = 1, nt
1454         is = (kk-1) * n
1455         DO K=n-1,1,-1
1456           DO I=1,L
1457             AT(I,K+is) = AT(I,K+is) - AU(I,K)*AT(I,K+is+1)
1458           ENDDO
1459         ENDDO
1460       ENDDO
1461 !-----------------------------------------------------------------------
1462       RETURN
1463       END SUBROUTINE TRIDIT
1464
1465 !-----------------------------------------------------------------------
1466
1467       END MODULE module_bl_gfs