phys/module_madwrf.F

   1   module module_madwrf
   2
   3   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   4   !                                          !
   5   ! Purpose: MAD-WRF model components        !
   6   !                                          !
   7   ! Author: Pedro A. Jimenez                 !
   8   !                                          !
   9   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  10
  11     use module_model_constants, only : G, T0, RCP
  12 #if ( ! NMM_CORE == 1 )
  13     use module_soil_pre, only : Skip_middle_points_t
  14 #endif
  15
  16     implicit none
  17
  18     private
  19     public :: Init_madwrf_tracers, Init_madwrf_clouds
  20
  21   contains
  22
  23     function Calc_cldtopz_from_brtemp (cldmask, ts, dzs, brtemp, tropoz, ht, kts, kte) result (return_value)
  24
  25     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  26     !                                                                                         !
  27     ! Purpose: Calculate cloud top height from brigtness temperature                          !
  28     !                                                                                         !
  29     ! Authors: Pedro A. Jimenez version of Greg Thompson code                                 !
  30     !                                                                                         !
  31     ! Comments: Staring at model top, go downwards first to the level of the tropopause,      !
  32     !           then keep going until the model temperature is greater (warmer) than the      !
  33     !           incoming GOES satellite longwave IR tempeature (channel13 on GOES-R series).  !
  34     !           As version1.0, this will be the highest altitude of potential cloud based on  !
  35     !           IR temperature whereas an improvement is possible for inversions in which the !
  36     !           cloud top could be lower altitude that might be detected using the RH field.  !
  37     !                                                                                         !
  38     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  39
  40       implicit none
  41
  42       real :: return_value
  43
  44       real, intent (in) :: cldmask, brtemp, tropoz, ht
  45       integer, intent (in) :: kts, kte
  46       real, dimension (kts:kte), intent (in) :: ts, dzs
  47
  48       real, parameter :: MISSING_CLDTOPZ = -999.9, MISSING_BRTEMP = -999.9
  49       real :: this_height
  50       integer :: k
  51
  52
  53       return_value = MISSING_CLDTOPZ
  54       if (brtemp == MISSING_BRTEMP) return
  55       if (cldmask > 0.0) then
  56         this_height = sum (dzs(kts:kte - 1)) + ht
  57         do k = kte - 2, kts, -1
  58           this_height = this_height - dzs(k + 1)
  59           if (this_height > tropoz) cycle
  60           if (ts(k) > brtemp) then
  61             return_value = this_height + dzs(k) * (brtemp - ts(k)) / (ts(k + 1) - ts(k))
  62             if ((return_value < this_height) .or. (return_value > this_height + dzs(k))) &
  63                 return_value = tropoz
  64             exit
  65           end if
  66         end do
  67       end if
  68
  69     end function Calc_cldtopz_from_brtemp
  70
  71     subroutine Init_madwrf_clouds (moist, p_qv, p_qc, p_qi, p_qs, p00, t, p, ph_2, phb, alt, xland, cldmask, cldtopz, cldbasez, &
  72          brtemp, ht,  dx, dy, flag_cldmask, flag_cldtopz, flag_cldbasez, flag_brtemp, em_width, hold_ups, ids, ide, jds, jde,   &
  73          its, ims, ime, jms, jme, kms, kme, ite, jts, jte, kts, kte, cldfra)
  74
  75     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  76     !                                          !
  77     ! Purpose: Cloud initialization            !
  78     !                                          !
  79     ! Author: Pedro A. Jimenez                 !
  80     !                                          !
  81     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  82
  83       implicit none
  84
  85         ! t, p, ph_2, phb, alt, xland
  86       real, dimension (ims:ime, kms:kme, jms:jme, *), intent (inout) :: moist
  87       real, dimension (ims:ime, kms:kme, jms:jme), intent (in) :: t, p, ph_2, phb, alt
  88       real, dimension (ims:ime, jms:jme), intent (in) :: xland, cldbasez, brtemp, ht
  89       real, dimension (ims:ime, jms:jme), intent (inout) :: cldmask, cldtopz
  90       logical, intent (in) :: hold_ups
  91       integer, intent (in) :: p_qv, p_qc, p_qi, p_qs, flag_cldmask, flag_cldtopz, flag_cldbasez, flag_brtemp, em_width, &
  92           ids, ide, jds, jde, ims, ime, jms, jme, kms, kme, its, ite, jts, jte, kts, kte
  93       real, intent (in) :: p00, dx, dy
  94       real, dimension (ims:ime, kms:kme, jms:jme), intent (out) :: cldfra
  95
  96       real, parameter :: CONVERT_M_TO_KM = 0.001, MISSING_CLDTOPZ = -999.9, MISSING_CLDMASK = -999.9
  97       real, dimension (kts:kte - 1) :: dzs, ts
  98       real :: gridkm, tropoz, cldtopz_agl
  99       integer :: i, j, k, insert_clouds, k_tropo
 100       logical, parameter :: LOCAL_DEBUG = .false.
 101
 102
 103         ! Identify the cloud initialization method
 104       insert_clouds = 0
 105       if (flag_cldmask == 1 .and. flag_cldtopz == 0 .and. flag_brtemp == 0) insert_clouds = 1
 106       if ((flag_cldmask == 1 .and. flag_brtemp == 1) .or. flag_cldtopz == 1) insert_clouds = 2
 107       if (insert_clouds == 2 .and. flag_cldbasez == 1) insert_clouds = 3
 108       if (LOCAL_DEBUG) then
 109         print *, 'flag_cldmask = ', flag_cldmask
 110         print *, 'flag_cldtopz = ', flag_cldtopz
 111         print *, 'flag_brtemp = ', flag_brtemp
 112         print *, 'insert_clouds = ', insert_clouds
 113       end if
 114
 115       gridkm = sqrt (dx * dy) * CONVERT_M_TO_KM
 116
 117       Loop_j: do j = jts, min (jte, jde - 1)
 118         Loop_i: do i = its, min (ite, ide - 1)
 119 #if ( ! NMM_CORE == 1 )
 120           if (Skip_middle_points_t (ids, ide, jds, jde, i, j, em_width, hold_ups)) cycle
 121 #endif
 122
 123           if (LOCAL_DEBUG) print *, 'Calculations for i, j = ', i, j
 124
 125           Loop_k: do k = kts, kte - 1
 126               ! Convert theta pert to T
 127             ts(k) = (t(i, k, j) + T0) / ((p00 / p(i, k, j)) ** RCP)
 128               ! Calc layer thickness
 129             dzs(k) = (ph_2(i, k + 1, j) + phb(i, k + 1, j) - (ph_2(i, k, j) + phb(i, k, j))) / G
 130           end do Loop_k
 131
 132             ! Calc troposphere height
 133           call Calc_tropo_height (ts, p(i, :, j), dzs, kts, kte, LOCAL_DEBUG, k_tropo, tropoz)
 134           tropoz = tropoz + ht(i, j)
 135
 136             ! Calc cloud top height agl if necessary
 137           if (insert_clouds > 0) then
 138             cldtopz_agl = MISSING_CLDTOPZ
 139             if (flag_cldtopz == 1) then
 140                 ! cloud mask
 141               if (cldtopz(i, j) > 0.0) then
 142                 cldmask(i, j) = 1.0
 143               elseif (cldtopz(i, j) == 0.0) then
 144                 cldmask(i, j) = 0.0
 145               else
 146                 cldmask(i, j) = MISSING_CLDMASK
 147               end if
 148                 ! cloud top height
 149               if (cldtopz(i, j) > 0.0) cldtopz_agl = max (0.0, cldtopz(i, j) -  ht(i, j))
 150             else if (flag_brtemp == 1) then
 151                 ! Only cloud top height needed
 152               cldtopz(i, j) = Calc_cldtopz_from_brtemp (cldmask(i, j), ts, dzs, brtemp(i, j), tropoz, ht(i, j), kts, kte)
 153               if (cldtopz(i, j) > 0.0) cldtopz_agl = max (0.0, cldtopz(i, j) -  ht(i, j))
 154             end if
 155           end if
 156
 157           select_cld_impro: select case (insert_clouds)
 158             case (0)
 159                 ! Use the analysis field
 160               call cal_cldfra3_madwrf (cldfra(i, :, j), moist(i, :, j, p_qv), moist(i, :, j, p_qc), &
 161                   moist(i, :, j, p_qi), moist(i, :, j, p_qs), dzs, p(i, :, j), ts(:),   &
 162                   xland(i, j), gridkm, .true., 1.5, tropoz, kts, kte, .false.)
 163
 164             case (1)
 165                 ! Remove clouds if clear (cldmask = 0)
 166               call cal_cldfra3_madwrf (cldfra(i, :, j), moist(i, :, j, p_qv), moist(i, :, j, p_qc), &
 167                   moist(i, :, j, p_qi), moist(i, :, j, p_qs), dzs, p(i, :, j), ts(:),   &
 168                   xland(i, j), gridkm, .true., 1.5, tropoz, kts, kte, .false., &
 169                   cldmask = cldmask(i, j))
 170
 171             case (2)
 172                 ! Remove clouds if clear (cldmask = 0)
 173                 ! Reduce / extend cloud top heights to match observations
 174                 ! Add clouds over clear sky regions (cldmask = 1)
 175               call cal_cldfra3_madwrf (cldfra(i, :, j), moist(i, :, j, p_qv), moist(i, :, j, p_qc), &
 176                   moist(i, :, j, p_qi), moist(i, :, j, p_qs), dzs, p(i, :, j), ts(:),   &
 177                   xland(i, j), gridkm, .true., 1.5, tropoz, kts, kte, .false., &
 178                   cldmask = cldmask(i, j), cldtopz = cldtopz_agl)
 179
 180             case (3)
 181                 ! Remove clouds if clear (cldmask = 0)
 182                 ! Reduce / extend cloud top / base heights to match observations
 183               call cal_cldfra3_madwrf (cldfra(i, :, j), moist(i, :, j, p_qv), moist(i, :, j, p_qc), &
 184                   moist(i, :, j, p_qi), moist(i, :, j, p_qs), dzs, p(i, :, j), ts(:),   &
 185                   xland(i, j), gridkm, .true., 1.5, tropoz, kts, kte, .false.,  &
 186                   cldmask = cldmask(i, j), cldtopz = cldtopz_agl, cldbasez = cldbasez(i, j))
 187
 188           end select select_cld_impro
 189
 190 !               do k = kts, kte - 1
 191 !                  moist(i,k,j,P_QV) = MAX(temp_Qv(k), moist(i,k,j,P_QV))
 192 !                  moist(i,k,j,P_QC) = temp_Qc(k) / temp_R(k)
 193 !                  moist(i,k,j,P_QI) = temp_Qi(k) / temp_R(k)
 194 !                  if (P_QNI .gt. 1) then
 195 !                     scalar(i,k,j,P_QNI) = temp_Ni(k) / temp_R(k)
 196 !                  endif
 197 !                  if (P_QNC .gt. 1) then
 198 !                     scalar(i,k,j,P_QNC) = temp_Nc(k) / temp_R(k)
 199 !                  endif
 200 !               enddo
 201         end do Loop_i
 202       end do Loop_j
 203
 204     end subroutine Init_madwrf_clouds
 205
 206     subroutine Init_madwrf_tracers (tracer, moist, p_qc, p_qi, p_qs, p_tr_qc, p_tr_qi, p_tr_qs, &
 207         ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, its , ite , jts , jte , kts , kte)
 208
 209     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 210     !                                          !
 211     ! Purpose: Initializes MAD-WRF tracers     !
 212     !                                          !
 213     ! Author: Pedro A. Jimenez                 !
 214     !                                          !
 215     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 216
 217       implicit none
 218
 219       real, dimension (ims:ime, kms:kme, jms:jme, *), intent (inout) :: tracer
 220       real, dimension (ims:ime, kms:kme, jms:jme, *), intent (in) :: moist
 221       integer, intent (in) :: p_qc, p_qi, p_qs, p_tr_qc, p_tr_qi, p_tr_qs, ids, ide, jds, jde, &
 222           kds, kde, ims, ime, jms, jme, kms, kme, its , ite , jts , jte , kts , kte
 223
 224       integer :: i, j, k
 225
 226
 227       do j = jts, min (jte, jde - 1)
 228         do k = kts, kte - 1
 229           do i = its, min (ite, ide - 1)
 230             tracer(i, k, j, p_tr_qc) = moist(i, k, j, p_qc)
 231             tracer(i, k, j, p_tr_qi) = moist(i, k, j, p_qi)
 232             tracer(i, k, j, p_tr_qs) = moist(i, k, j, p_qs)
 233           end do
 234         end do
 235       end do
 236
 237     end subroutine Init_madwrf_tracers
 238
 239 !+---+-----------------------------------------------------------------+
 240
 241       SUBROUTINE cal_cldfra3_madwrf(CLDFRA, qv, qc, qi, qs, dz,                &
 242      &                 p, t, XLAND, gridkm,                             &
 243      &                 modify_qvapor, max_relh,                         &
 244      &                 tropo_z, kts,kte, debug_flag, k_tropo, cldmask,  &
 245      &                 cldtopz, cldbasez)
 246 !
 247       USE module_mp_thompson   , ONLY : rsif, rslf
 248       IMPLICIT NONE
 249 !
 250       INTEGER, INTENT(IN):: kts, kte
 251       LOGICAL, INTENT(IN):: modify_qvapor
 252       REAL, DIMENSION(kts:kte), INTENT(INOUT):: qv, qc, qi, qs, cldfra
 253       REAL, DIMENSION(kts:kte), INTENT(IN):: p, t, dz
 254       REAL, INTENT(IN):: gridkm, XLAND, max_relh
 255       REAL, INTENT(INOUT):: tropo_z
 256       LOGICAL, INTENT(IN):: debug_flag
 257       integer, intent(in), optional :: k_tropo
 258       real, intent(in), optional :: cldmask, cldtopz, cldbasez
 259
 260 !..Local vars.
 261       REAL:: RH_00L, RH_00O, RH_00
 262       REAL:: entrmnt=0.5
 263       INTEGER:: k
 264       REAL:: TC, qvsi, qvsw, RHUM, delz, cldbasek_tmp
 265       REAL, DIMENSION(kts:kte):: qvs, rh, rhoa
 266
 267       character*512 dbg_msg
 268       logical :: is_tropo_init, impose_cldmask, impose_cldtopz, &
 269           impose_cldbasetopz, keep_clouds_below_lcl
 270       INTEGER :: cldtopk, cldbasek, cldthick_bot_k, cldfra_thresh_k
 271       REAL, PARAMETER :: CLDTHICK_DEF = -999.9  ! Default thickness [m] of new clouds
 272       REAL, PARAMETER :: CLDFRA_DEF = 0.5  ! Default cloud fraction to insert in new clouds
 273       REAL, PARAMETER :: CLDFRA_THRESH = 0.0  ! Default threshold cloud fraction to determine presence of clouds
 274       REAL, PARAMETER :: RH_NOCLOUD = 0.3
 275       logical, parameter :: LOCAL_DEBUG = .false.
 276
 277
 278         ! Define logical vars to handle optional
 279         ! arguments
 280       if (present(k_tropo)) then
 281         is_tropo_init = .true.
 282       else
 283         is_tropo_init = .false.
 284       end if
 285
 286       keep_clouds_below_lcl = .false.
 287       impose_cldmask = .false.
 288       impose_cldtopz = .false.
 289       impose_cldbasetopz = .false.
 290
 291       if (present(cldmask) .and. .not. present(cldtopz) .and. .not. present(cldbasez)) then
 292         impose_cldmask = .true.
 293       elseif (present(cldmask) .and. present(cldtopz) .and. .not. present(cldbasez)) then
 294         impose_cldtopz = .true.
 295       elseif (present(cldmask) .and. present(cldtopz) .and. present(cldbasez)) then
 296         impose_cldbasetopz = .true.
 297       end if
 298
 299       !! Initialize variables
 300       cldtopk = -999
 301       cldbasek = -999
 302       cldthick_bot_k = -999
 303       cldfra_thresh_k = -999
 304
 305         ! Remove hydrometeors if necessary and
 306         ! calculate height related variables
 307       if (impose_cldmask) then
 308         if (cldmask == 0.0) then
 309           do k = kts, kte - 1
 310             cldfra(k) = 0.0
 311             qc(k) = 0.0
 312             qi(k) = 0.0
 313             qs(k) = 0.0
 314           end do
 315           return
 316         end if
 317       end if
 318
 319       if_impose_cldtopz: if (impose_cldtopz .or. impose_cldbasetopz) then
 320          if_cldmask: if (cldmask == 0.0) then
 321             do k = kts, kte - 1
 322                cldfra(k) = 0.0
 323                qc(k) = 0.0
 324                qi(k) = 0.0
 325                qs(k) = 0.0
 326             end do
 327             return
 328          else if (cldmask > 0.0 .and. cldmask <= 1.0 .and. cldtopz > 0.0) then
 329               ! Zero out hydrometeors above cldtopz
 330             call find_k_from_z_agl(cldtopz, cldtopk, dz, kts, kte)
 331             if (cldtopk > 0 .and. cldtopk < kte) then
 332               do k = cldtopk + 1, kte - 1
 333                  cldfra(k) = 0.0
 334                  qc(k) = 0.0
 335                  qi(k) = 0.0
 336                  qs(k) = 0.0
 337               end do
 338             end if
 339
 340               ! Calculate the cloud bottom k if the specified thickness were to be applied
 341               ! If that would be below ground, then set the k level to 1
 342             if (CLDTHICK_DEF > 0.0) then
 343               if (cldtopz - CLDTHICK_DEF > 0.0) then
 344                  call find_k_from_z_agl(cldtopz - CLDTHICK_DEF, cldthick_bot_k, dz, kts, kte)
 345               else
 346                  cldthick_bot_k = 1
 347               end if
 348             end if
 349          end if if_cldmask
 350       end if if_impose_cldtopz
 351
 352
 353       if (impose_cldbasetopz) then
 354          !! If cloud base information (cldbzin) is also provided in met_em files, then...
 355          !! Zero out hydrometeors below observed cloud base in columns where cldmask > 0.0
 356          !! Note that hydrometeors have already been zeroed out in columns where cldmask = 0.0
 357          if (cldmask > 0.0 .and. cldmask <= 1.0 .and. cldbasez > 0.0) then
 358             ! First, find out level of observed cloud base
 359             call find_k_from_z_agl(cldbasez, cldbasek, dz, kts, kte)
 360             !! Trust satellite estimations of cloud top height more than METAR estimates/extrapolations of cloud base height
 361             !! Thus, if cldtopk < cldbasek, then pretend we don't have cldbasek
 362             if (cldtopk < cldbasek) cldbasek = -999
 363             if (cldbasek > 0 .and. cldbasek > kts) then
 364                do k = kts, cldbasek-1
 365                   cldfra(k) = 0.0
 366                   qc(k) = 0.0
 367                   qi(k) = 0.0
 368                   qs(k) = 0.0
 369                end do
 370             end if
 371          end if
 372       end if
 373
 374 !+---+
 375
 376 !..Initialize cloud fraction, compute RH, and rho-air.
 377
 378          DO k = kts, kte - 1
 379             CLDFRA(K) = 0.0
 380             tc = t(k) - 273.15
 381
 382             qvsw = rslf(P(k), t(k))
 383             if (debug_flag) print *, 'k, p, t, qvsw, qv =', k, p(k), t(k), qvsw, qv(k)
 384             if (tc .lt. 0.0) then
 385                qvsi = rsif(P(k), t(k)+0.025)    !..To ensure a tiny bit ice supersaturated
 386             else
 387                qvsi = qvsw
 388             endif
 389
 390             if (tc .ge. -12.0) then
 391                qvs(k) = qvsw
 392             elseif (tc .lt. -35.0) then
 393                qvs(k) = qvsi
 394             else
 395                qvs(k) = qvsw - (qvsw-qvsi)*(-12.0-tc)/(-12.0+35.)
 396             endif
 397
 398             if (modify_qvapor) then
 399                if (qc(k).gt.1.E-8) then
 400                   qv(k) = MAX(qv(k), qvsw)
 401                   qvs(k) = qvsw
 402                endif
 403                if (qc(k).le.1.E-8 .and. qi(k).ge.1.E-9) then
 404                   qv(k) = MAX(qv(k), qvsi)
 405                   qvs(k) = qvsi
 406                endif
 407             endif
 408
 409
 410             rh(k) = MAX(0.01, qv(k)/qvs(k))
 411             if (debug_flag) print *, '  qv, qvs, qvsi, qc, qi, rh =', qv(k), qvs(k), qvsi, qc(k), qi(k), rh(k)
 412             rhoa(k) = p(k)/(287.0*t(k))
 413          ENDDO
 414
 415           ! Over ocean, find the level at which RH is lower than RH_NOCLOUD starting from cloud top
 416           ! to avoid putting clouds in very dry layers
 417         if ((impose_cldtopz .or. impose_cldbasetopz) .and. (XLAND - 1.5) > 0.0) then
 418           if (cldtopk > 0 .and. cldtopk < kte) then
 419             if (cldbasek < 0 .or. cldbasek > cldtopk) then
 420               if (rh(cldtopk) < RH_NOCLOUD) then
 421                 cldbasek = cldtopk + 1
 422               else
 423                 cldbasek = cldtopk
 424               end if
 425             else
 426               cldbasek_tmp = cldtopk + 1
 427               do k = cldtopk, kts, -1
 428                 if (rh(k) < RH_NOCLOUD) exit
 429                   cldbasek_tmp = k
 430               end do
 431               if (cldbasek_tmp > cldbasek) cldbasek = cldbasek_tmp
 432             end if
 433           end if
 434         end if
 435
 436 !..First cut scale-aware. Higher resolution should require closer to
 437 !.. saturated grid box for higher cloud fraction.  Simple functions
 438 !.. chosen based on Mocko and Cotton (1995) starting point and desire
 439 !.. to get near 100% RH as grid spacing moves toward 1.0km, but higher
 440 !.. RH over ocean required as compared to over land.
 441
 442       DO k = kts, kte - 1
 443
 444          delz = MAX(100., dz(k))
 445          RH_00L = 0.65 + SQRT(1./(25.0+gridkm*gridkm*delz*0.01))
 446          RH_00O = 0.81 + SQRT(1./(50.0+gridkm*gridkm*delz*0.01))
 447          RHUM = rh(k)
 448
 449          if (qc(k).gt.1.E-7 .or. qi(k).ge.1.E-7                         &
 450      &                    .or. (qs(k).gt.1.E-6 .and. t(k).lt.273.)) then
 451             CLDFRA(K) = 1.0
 452          else
 453
 454             IF ((XLAND-1.5).GT.0.) THEN                                  !--- Ocean
 455                RH_00 = RH_00O
 456             ELSE                                                         !--- Land
 457                RH_00 = RH_00L
 458             ENDIF
 459
 460             tc = t(k) - 273.15
 461             if (tc .lt. -12.0) RH_00 = RH_00L
 462
 463             if (tc .ge. 25.0) then
 464                CLDFRA(K) = 0.0
 465             elseif (tc .ge. -12.0) then
 466                RHUM = MIN(rh(k), 1.0)
 467                CLDFRA(K) = MAX(0., 1.0-SQRT((1.005-RHUM)/(1.005-RH_00)))
 468             else
 469                if (max_relh.gt.1.12 .or. (.NOT.(modify_qvapor)) ) then
 470 !..For HRRR model, the following look OK.
 471                   RHUM = MIN(rh(k), 1.45)
 472                   RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+100.)
 473                   CLDFRA(K) = MAX(0., 1.0-SQRT((1.5-RHUM)/(1.5-RH_00)))
 474                else
 475 !..but for the GFS model, RH is way lower.
 476                   RHUM = MIN(rh(k), 1.05)
 477                   RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+100.)
 478                   CLDFRA(K) = MAX(0., 1.0-SQRT((1.05-RHUM)/(1.05-RH_00)))
 479                endif
 480             endif
 481             if (CLDFRA(K).gt.0.) CLDFRA(K) = MAX(0.01, MIN(CLDFRA(K),0.9))
 482
 483             if (debug_flag) then
 484               WRITE (dbg_msg,*) 'DEBUG-GT: cloud fraction: ', RH_00, RHUM, CLDFRA(K)
 485               CALL wrf_debug (150, dbg_msg)
 486             endif
 487
 488          endif
 489       ENDDO
 490
 491       if ((impose_cldtopz .or. impose_cldbasetopz) .and. cldtopk > 0 .and. cldtopk < kte) then
 492           ! Remove clouds above observed cloud top
 493         CLDFRA(cldtopk + 1:kte) = 0.0
 494
 495         keep_clouds_below_lcl = .true.
 496
 497           ! Find the first cloudy level
 498         call find_thresh_k_downward(cldfra, CLDFRA_THRESH, cldfra_thresh_k, cldtopk, kts, kts, kte)
 499         if (cldfra_thresh_k > 0) then
 500              ! Extend cloud until cldtopz
 501            if (cldfra_thresh_k < cldtopk) cldfra(cldfra_thresh_k + 1:cldtopk) = CLDFRA_DEF
 502         else
 503             ! No clouds present
 504           if (CLDTHICK_DEF > 0.0) then
 505             cldfra(cldthick_bot_k:cldtopk) = CLDFRA_DEF
 506 !          else
 507 !            if (rh(cldtopk) > RH_NOCLOUD) cldfra(cldtopk) = CLDFRA_DEF
 508           end if
 509         end if
 510       end if
 511
 512       if (impose_cldbasetopz .and. cldbasek > kts .and. cldbasek <= kte) then
 513          !! Remove clouds below observed cloud base
 514          CLDFRA(kts:cldbasek-1) = 0.0
 515
 516          !! If necessary, add clouds with default CLDFRA from observed cloud base up to base of already-existing cloud
 517          cldfra_thresh_k = -999  ! reset this variable
 518          call find_thresh_k_upward(cldfra, CLDFRA_THRESH, cldfra_thresh_k, cldbasek, kte, kts, kte)
 519          if (cldfra_thresh_k > cldbasek) then
 520            do k = cldbasek, cldfra_thresh_k - 1
 521              if (rh(k) > RH_NOCLOUD) cldfra(k) = CLDFRA_DEF
 522            end do
 523          end if
 524       end if
 525
 526
 527       if (is_tropo_init) then
 528         call find_cloudLayers(qvs, cldfra, T, P, Dz, entrmnt,             &
 529        &                      debug_flag, qc, qi, qs, tropo_z, kts,kte, keep_clouds_below_lcl, k_tropo)
 530       else
 531         call find_cloudLayers(qvs, cldfra, T, P, Dz, entrmnt,             &
 532        &                      debug_flag, qc, qi, qs, tropo_z, kts,kte, keep_clouds_below_lcl)
 533       end if
 534
 535 !..Do a final total column adjustment since we may have added more than 1mm
 536 !.. LWP/IWP for multiple cloud decks.
 537
 538       call adjust_cloudFinal(cldfra, qc, qi, rhoa, dz, kts,kte)
 539
 540       if (debug_flag) then
 541         WRITE (dbg_msg,*) 'DEBUG-GT:  Made-up fake profile of clouds'
 542         CALL wrf_debug (150, dbg_msg)
 543         do k = kte - 1, kts, -1
 544            write(dbg_msg,'(f7.2, 2x, f7.2, 2x, f6.4, 2x, f7.3, 2x,  f15.7, 2x, f15.7)') &
 545      &          T(k)-273.15, P(k)*0.01, rh(k), cldfra(k)*100., qc(k)*1000.,qi(k)*1000.
 546            CALL wrf_debug (150, dbg_msg)
 547         enddo
 548       endif
 549
 550       if (modify_qvapor) then
 551          DO k = kts, kte - 1
 552             if (cldfra(k).gt.0.20 .and. cldfra(k).lt.1.0) then
 553                qv(k) = MAX(qv(k),qvs(k))
 554             endif
 555          ENDDO
 556       endif
 557
 558       END SUBROUTINE cal_cldfra3_madwrf
 559
 560       SUBROUTINE find_k_from_z_agl(z_agl, k_lev, dz, kts, kte)
 561
 562       IMPLICIT NONE
 563
 564       INTEGER, INTENT(IN)  :: kts, kte
 565       REAL, INTENT(IN)     :: z_agl
 566       INTEGER, INTENT(OUT) :: k_lev
 567       REAL, DIMENSION(kts:kte), INTENT(IN) :: dz
 568       INTEGER :: k
 569       REAL    :: z_this, z_next, z_full_lev
 570
 571       !! Identify the k-level of a given height AGL, starting from ground level
 572       z_full_lev = 0.0
 573       do k = kts, kte
 574          z_this = z_full_lev
 575          z_next = z_this + dz(k)
 576          if (z_agl > 0.0) then
 577             if (z_agl > z_this .and. z_agl <= z_next) then
 578                k_lev = k
 579             end if
 580          end if
 581          z_full_lev = z_next
 582       end do
 583
 584       END SUBROUTINE find_k_from_z_agl
 585
 586       SUBROUTINE find_thresh_k_downward(var, var_thresh, k_lev, k_top, k_bot, kts, kte)
 587
 588       IMPLICIT NONE
 589
 590       INTEGER, INTENT(IN)  :: kts, kte
 591       INTEGER, INTENT(IN)  :: k_top, k_bot
 592       REAL, DIMENSION(kts:kte), INTENT(IN) :: var
 593       REAL, INTENT(IN)     :: var_thresh
 594       INTEGER, INTENT(OUT) :: k_lev
 595       INTEGER :: k
 596
 597       !! Identify the k-level where the input variable first exceeds a threshold, searching downward over a given range
 598       do k = k_top, k_bot, -1
 599          if (var(k) > var_thresh) then
 600             k_lev = k
 601             exit
 602          end if
 603       end do
 604
 605       END SUBROUTINE find_thresh_k_downward
 606
 607       SUBROUTINE find_thresh_k_upward(var, var_thresh, k_lev, k_bot, k_top, kts, kte)
 608
 609       IMPLICIT NONE
 610
 611       INTEGER, INTENT(IN)  :: kts, kte
 612       INTEGER, INTENT(IN)  :: k_top, k_bot
 613       REAL, DIMENSION(kts:kte), INTENT(IN) :: var
 614       REAL, INTENT(IN)     :: var_thresh
 615       INTEGER, INTENT(OUT) :: k_lev
 616       INTEGER :: k
 617
 618       !! Identify the k-level where the input variable first exceeds a threshold, searching upward over a given range
 619       do k = k_bot, k_top
 620          if (var(k) > var_thresh) then
 621             k_lev = k
 622             exit
 623          end if
 624       end do
 625
 626       END SUBROUTINE
 627
 628 !+---+-----------------------------------------------------------------+
 629 !..From cloud fraction array, find clouds of multi-level depth and compute
 630 !.. a reasonable value of LWP or IWP that might be contained in that depth,
 631 !.. unless existing LWC/IWC is already there.
 632
 633       SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,&
 634      &                            debugfl, qc1d, qi1d, qs1d,            &
 635      &                            tropo_z, kts,kte, keep_clouds_below_lcl, ktropo)
 636 !
 637       IMPLICIT NONE
 638 !
 639       INTEGER, INTENT(IN):: kts, kte
 640       LOGICAL, INTENT(IN):: debugfl, keep_clouds_below_lcl
 641       REAL, INTENT(IN):: entrmnt
 642       REAL, INTENT(INOUT):: tropo_z
 643       REAL, DIMENSION(kts:kte), INTENT(IN):: qs1d,qvs1d,T1d,P1d,Dz1d
 644       REAL, DIMENSION(kts:kte), INTENT(INOUT):: cfr1d, qc1d, qi1d
 645       integer, intent(in), optional :: ktropo
 646
 647 !..Local vars.
 648       REAL, DIMENSION(kts:kte):: theta
 649       REAL:: theta1, theta2, delz
 650       INTEGER:: k, k2, k_tropo, k_m12C, k_m40C, k_cldb, k_cldt, kbot
 651       LOGICAL:: in_cloud
 652       character*512 dbg_msg
 653       logical :: is_tropo_init
 654
 655 !+---+
 656
 657       if (present(ktropo)) then
 658         is_tropo_init = .true.
 659       else
 660         is_tropo_init = .false.
 661       end if
 662
 663       k_m12C = 0
 664       k_m40C = 0
 665       DO k = kte - 1, kts, -1
 666          theta(k) = T1d(k)*((100000.0/P1d(k))**(287.05/1004.))
 667          if (T1d(k)-273.16 .gt. -40.0 .and. P1d(k).gt.7000.0) k_m40C = MAX(k_m40C, k)
 668          if (T1d(k)-273.16 .gt. -12.0 .and. P1d(k).gt.10100.0) k_m12C = MAX(k_m12C, k)
 669       ENDDO
 670       if (k_m40C .le. kts) k_m40C = kts
 671       if (k_m12C .le. kts) k_m12C = kts
 672
 673       if (k_m40C.gt.kte-2 .OR. k_m12C.gt.kte-3) then
 674         WRITE (dbg_msg,*) 'DEBUG-GT: WARNING, no possible way neg40 or neg12C can occur this high up: ', k_m40C, k_m12C
 675         CALL wrf_debug (0, dbg_msg)
 676         do k = kte - 1, kts, -1
 677            WRITE (dbg_msg,*) 'DEBUG-GT,   P, T : ', P1d(k)*0.01,T1d(k)-273.16
 678            CALL wrf_debug (0, dbg_msg)
 679         enddo
 680         call wrf_error_fatal ('FATAL ERROR, problem in temperature profile.')
 681       endif
 682
 683       if (is_tropo_init) then
 684        k_tropo = ktropo
 685       else
 686         call Calc_tropo_height (T1d, P1d, dz1d, kts, &
 687           kte, debugfl, k_tropo, tropo_z)
 688       end if
 689
 690 !..Eliminate possible fractional clouds above supposed tropopause.
 691       DO k = k_tropo+1, kte - 1
 692          if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.1.0) then
 693 !            print *, 'Removing stratospheric cloud', k
 694             cfr1d(k) = 0.
 695          endif
 696       ENDDO
 697
 698 !..We would like to prevent fractional clouds below LCL in idealized
 699 !.. situation with deep well-mixed convective PBL, that otherwise is
 700 !.. likely to get clouds in more realistic capping inversion layer.
 701       kbot = kts+2
 702       DO k = kbot, k_m12C
 703          if ( (theta(k)-theta(k-1)) .gt. 0.010E-3*Dz1d(k)) EXIT
 704       ENDDO
 705       kbot = MAX(kts+1, k-2)
 706       if (.not. keep_clouds_below_lcl) then
 707         DO k = kts, kbot
 708 !           print *, 'Reducing cloud below LCL', k
 709            if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.1.0) cfr1d(k) = 0.5*cfr1d(k)
 710         ENDDO
 711       else
 712         kbot = kts + 1
 713       end if
 714
 715
 716 !..Starting below tropo height, if cloud fraction greater than 1 percent,
 717 !.. compute an approximate total layer depth of cloud, determine a total
 718 !.. liquid water/ice path (LWP/IWP), then reduce that amount with tuning
 719 !.. parameter to represent entrainment factor, then divide up LWP/IWP
 720 !.. into delta-Z weighted amounts for individual levels per cloud layer.
 721
 722       k_cldb = k_tropo
 723       in_cloud = .false.
 724       k = k_tropo
 725       DO WHILE (.not. in_cloud .AND. k.gt.k_m12C+1)
 726          k_cldt = 0
 727          if (cfr1d(k).ge.0.01) then
 728             in_cloud = .true.
 729             k_cldt = MAX(k_cldt, k)
 730          endif
 731          if (in_cloud) then
 732             DO k2 = k_cldt-1, k_m12C, -1
 733                if (cfr1d(k2).lt.0.01 .or. k2.eq.k_m12C) then
 734                   k_cldb = k2+1
 735                   goto 87
 736                endif
 737             ENDDO
 738  87         continue
 739             in_cloud = .false.
 740          endif
 741          if ((k_cldt - k_cldb + 1) .ge. 2) then
 742       if (debugfl) then
 743         WRITE (dbg_msg,*) 'DEBUG-GT: An ice cloud layer is found between ', k_cldt, k_cldb, P1d(k_cldt)*0.01, P1d(k_cldb)*0.01
 744         CALL wrf_debug (150, dbg_msg)
 745       endif
 746             call adjust_cloudIce(cfr1d, qi1d, qs1d, qvs1d, T1d, Dz1d,   &
 747      &                           entrmnt, k_cldb,k_cldt,kts,kte)
 748             k = k_cldb
 749          elseif ((k_cldt - k_cldb + 1) .eq. 1) then
 750       if (debugfl) then
 751         WRITE (dbg_msg,*) 'DEBUG-GT: A single-layer ice cloud layer is found on ', k_cldb, P1d(k_cldb)*0.01
 752         CALL wrf_debug (150, dbg_msg)
 753       endif
 754             if (cfr1d(k_cldb).gt.0.and.cfr1d(k_cldb).lt.1.)             &
 755      &               qi1d(k_cldb)=0.05*qvs1d(k_cldb)*cfr1d(k_cldb)*cfr1d(k_cldb)
 756             k = k_cldb
 757          endif
 758          k = k - 1
 759       ENDDO
 760
 761
 762       k_cldb = k_m12C + 3
 763       in_cloud = .false.
 764       k = k_m12C + 2
 765       DO WHILE (.not. in_cloud .AND. k.gt.kbot)
 766          k_cldt = 0
 767          if (cfr1d(k).ge.0.01) then
 768             in_cloud = .true.
 769             k_cldt = MAX(k_cldt, k)
 770          endif
 771          if (in_cloud) then
 772             DO k2 = k_cldt-1, kbot, -1
 773                if (cfr1d(k2).lt.0.01 .or. k2.eq.kbot) then
 774                   k_cldb = k2+1
 775                   goto 88
 776                endif
 777             ENDDO
 778  88         continue
 779             in_cloud = .false.
 780          endif
 781          if ((k_cldt - k_cldb + 1) .ge. 2) then
 782       if (debugfl) then
 783         WRITE (dbg_msg,*) 'DEBUG-GT: A water cloud layer is found between ', k_cldt, k_cldb, P1d(k_cldt)*0.01, P1d(k_cldb)*0.01
 784         CALL wrf_debug (150, dbg_msg)
 785       endif
 786             call adjust_cloudH2O(cfr1d, qc1d, qvs1d, T1d, Dz1d,         &
 787      &                           entrmnt, k_cldb,k_cldt,kts,kte)
 788             k = k_cldb
 789          elseif ((k_cldt - k_cldb + 1) .eq. 1) then
 790             if (cfr1d(k_cldb).gt.0.and.cfr1d(k_cldb).lt.1.)             &
 791      &               qc1d(k_cldb)=0.05*qvs1d(k_cldb)*cfr1d(k_cldb)*cfr1d(k_cldb)
 792             k = k_cldb
 793          endif
 794          k = k - 1
 795       ENDDO
 796
 797       END SUBROUTINE find_cloudLayers
 798
 799       subroutine Calc_tropo_height (T1d, P1d, dz1d, kts, kte, debugfl, k_tropo, tropo_z)
 800
 801         !..Find tropopause height, best surrogate, because we would not really
 802         !.. wish to put fake clouds into the stratosphere.  The 10/1500 ratio
 803         !.. d(Theta)/d(Z) approximates a vertical line on typical SkewT chart
 804         !.. near typical (mid-latitude) tropopause height.  Since messy data
 805         !.. could give us a false signal of such a transition, do the check over
 806         !.. three K-level change, not just a level-to-level check.  This method
 807         !.. has potential failure in arctic-like conditions with extremely low
 808         !.. tropopause height, as would any other diagnostic, so ensure resulting
 809         !.. k_tropo level is above 700hPa.
 810
 811         implicit none
 812
 813         REAL, DIMENSION(kts:kte), INTENT(IN) :: T1d, P1d, Dz1d
 814         integer, intent(in) :: kts, kte
 815         LOGICAL, INTENT(IN):: debugfl
 816         integer, intent(out) :: k_tropo
 817         real, intent (out) :: tropo_z
 818
 819           ! Local vars
 820         integer :: k, k_p200
 821         real :: theta1, theta2, delz
 822         character*512 dbg_msg
 823         REAL, DIMENSION(kts:kte):: theta
 824
 825
 826         k_p200 = 0
 827         DO k = kte - 1, kts, -1
 828           theta(k) = T1d(k)*((100000.0/P1d(k))**(287.05/1004.))
 829           if (P1d(k).gt.19999.0 .and. k_p200.eq.0) k_p200 = k
 830         END DO
 831
 832         if ( (kte-k_p200) .lt. 3) k_p200 = kte-3
 833         DO k = k_p200-2, kts, -1
 834            theta1 = theta(k)
 835            theta2 = theta(k+2)
 836            delz = 0.5*dz1d(k) + dz1d(k+1) + 0.5*dz1d(k+2)
 837            if ( (((theta2-theta1)/delz).lt.10./1500.) .OR. P1d(k).gt.70000.) EXIT
 838         ENDDO
 839         k_tropo = MAX(kts+2, MIN(k+2, kte-1))
 840
 841         if (k_tropo .gt. k_p200) then
 842            DO k = kte-3, k_p200-2, -1
 843               theta1 = theta(k)
 844               theta2 = theta(k+2)
 845               delz = 0.5*dz1d(k) + dz1d(k+1) + 0.5*dz1d(k+2)
 846               if ( (((theta2-theta1)/delz).lt.10./1500.) .AND. P1d(k).gt.9500.) EXIT
 847            ENDDO
 848            k_tropo = MAX(k_p200-1, MIN(k+2, kte-1))
 849         endif
 850         tropo_z = SUM(dz1d(kts:k_tropo))
 851
 852         if (k_tropo.gt.kte-2) then
 853           WRITE (dbg_msg,*) 'DEBUG-GT: CAUTION, tropopause appears to be very high up: ', k_tropo
 854           CALL wrf_debug (150, dbg_msg)
 855           do k = kte - 1, kts, -1
 856              WRITE (dbg_msg,*) 'DEBUG-GT,   P, T : ', k,P1d(k)*0.01,T1d(k)-273.16
 857              CALL wrf_debug (150, dbg_msg)
 858           enddo
 859         elseif (debugfl) then
 860           WRITE (dbg_msg,*) 'DEBUG-GT: FOUND TROPOPAUSE k=', k_tropo
 861           CALL wrf_debug (150, dbg_msg)
 862         endif
 863
 864         if (debugfl) then
 865           print *, 'FOUND TROPOPAUSE k, height=', k_tropo, tropo_z
 866         end if
 867
 868       end subroutine Calc_tropo_height
 869
 870       SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte)
 871 !
 872       IMPLICIT NONE
 873 !
 874       INTEGER, INTENT(IN):: k1,k2, kts,kte
 875       REAL, INTENT(IN):: entr
 876       REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, qs, qvs, T, dz
 877       REAL, DIMENSION(kts:kte), INTENT(INOUT):: qi
 878       REAL:: iwc, max_iwc, tdz, this_iwc, this_dz
 879       INTEGER:: k
 880
 881       tdz = 0.
 882       do k = k1, k2
 883          tdz = tdz + dz(k)
 884       enddo
 885       max_iwc = ABS(qvs(k2)-qvs(k1))
 886 !     print*, ' max_iwc = ', max_iwc, ' over DZ=',tdz
 887
 888       do k = k1, k2
 889          max_iwc = MAX(1.E-6, max_iwc - (qi(k)+qs(k)))
 890       enddo
 891       max_iwc = MIN(1.E-3, max_iwc)
 892
 893       this_dz = 0.0
 894       do k = k1, k2
 895          if (k.eq.k1) then
 896             this_dz = this_dz + 0.5*dz(k)
 897          else
 898             this_dz = this_dz + dz(k)
 899          endif
 900          this_iwc = max_iwc*this_dz/tdz
 901          iwc = MAX(1.E-6, this_iwc*(1.-entr))
 902          if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.203.16) then
 903             qi(k) = qi(k) + cfr(k)*cfr(k)*iwc
 904          endif
 905       enddo
 906
 907       END SUBROUTINE adjust_cloudIce
 908
 909 !+---+-----------------------------------------------------------------+
 910
 911       SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte)
 912 !
 913       IMPLICIT NONE
 914 !
 915       INTEGER, INTENT(IN):: k1,k2, kts,kte
 916       REAL, INTENT(IN):: entr
 917       REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, qvs, T, dz
 918       REAL, DIMENSION(kts:kte), INTENT(INOUT):: qc
 919       REAL:: lwc, max_lwc, tdz, this_lwc, this_dz
 920       INTEGER:: k
 921
 922       tdz = 0.
 923       do k = k1, k2
 924          tdz = tdz + dz(k)
 925       enddo
 926       max_lwc = ABS(qvs(k2)-qvs(k1))
 927 !     print*, ' max_lwc = ', max_lwc, ' over DZ=',tdz
 928
 929       do k = k1, k2
 930          max_lwc = MAX(1.E-6, max_lwc - qc(k))
 931       enddo
 932       max_lwc = MIN(1.E-3, max_lwc)
 933
 934       this_dz = 0.0
 935       do k = k1, k2
 936          if (k.eq.k1) then
 937             this_dz = this_dz + 0.5*dz(k)
 938          else
 939             this_dz = this_dz + dz(k)
 940          endif
 941          this_lwc = max_lwc*this_dz/tdz
 942          lwc = MAX(1.E-6, this_lwc*(1.-entr))
 943          if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.253.16) then
 944             qc(k) = qc(k) + cfr(k)*cfr(k)*lwc
 945          endif
 946       enddo
 947
 948       END SUBROUTINE adjust_cloudH2O
 949
 950 !+---+-----------------------------------------------------------------+
 951
 952 !..Do not alter any grid-explicitly resolved hydrometeors, rather only
 953 !.. the supposed amounts due to the cloud fraction scheme.
 954
 955       SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte)
 956 !
 957       IMPLICIT NONE
 958 !
 959       INTEGER, INTENT(IN):: kts,kte
 960       REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, Rho, dz
 961       REAL, DIMENSION(kts:kte), INTENT(INOUT):: qc, qi
 962       REAL:: lwp, iwp, xfac
 963       INTEGER:: k
 964
 965       lwp = 0.
 966       iwp = 0.
 967       do k = kts, kte - 1
 968          if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then
 969             lwp = lwp + qc(k)*Rho(k)*dz(k)
 970             iwp = iwp + qi(k)*Rho(k)*dz(k)
 971          endif
 972       enddo
 973
 974       if (lwp .gt. 1.0) then
 975          xfac = 1.0/lwp
 976          do k = kts, kte - 1
 977             if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then
 978                qc(k) = qc(k)*xfac
 979             endif
 980          enddo
 981       endif
 982
 983       if (iwp .gt. 1.0) then
 984          xfac = 1.0/iwp
 985          do k = kts, kte - 1
 986             if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then
 987                qi(k) = qi(k)*xfac
 988             endif
 989          enddo
 990       endif
 991
 992       END SUBROUTINE adjust_cloudFinal
 993
 994   end module module_madwrf