phys/module_diag_solar.F

   1 #if (NMM_CORE == 1)
   2 MODULE module_diag_solar
   3 CONTAINS
   4    SUBROUTINE solar_diag
   5    END SUBROUTINE solar_diag
   6 END MODULE module_diag_solar
   7 #else
   8 !WRF:MEDIATION_LAYER:PHYSICS
   9 !
  10
  11 MODULE module_diag_solar
  12
  13 CONTAINS
  14
  15 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  16 !                                                                                        !
  17 ! Diagnostics code to output variables relevant for WRF-Solar applications               !
  18 !                                                                                        !
  19 ! Written by Timothy W. Juliano and Pedro A. Jimenez (NCAR/RAL)                          !
  20 ! Some of the code has been modified from FARMS (module_ra_farms)                        !
  21 !                                                                                        !
  22 ! Activated using [solar_diagnostics = 1] option &diags section of namelist.input file   !
  23 !                                                                                        !
  24 ! Note that output variables with 'tot' include unresolved hydrometeors                  !
  25 ! Only subgrid qcloud and qice are accounted for                                         !
  26 !                                                                                        !
  27 ! Corresponding variables are located in registry.solar_fields                           !
  28 !                                                                                        !
  29 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  30
  31    SUBROUTINE solar_diag (dt, rho, dz8w, ph, phb, cldfrac3d, coszen, swdnb, swdnt, &
  32              param_first_scalar, p_qc, p_qi, p_qs, qv, qc, qi, qs, &
  33              qc_tot, qi_tot, has_reqc, has_reqi, has_reqs, f_qv, f_qc, f_qi, f_qs, &
  34              re_cloud, re_ice, re_snow, clrnidx, sza, ghi_accum, cldfrac2d, wvp2d, lwp2d, iwp2d, swp2d, &
  35              wp2d_sum, lwp2d_tot, iwp2d_tot, wp2d_tot_sum, re_cloud_path, re_ice_path, re_snow_path, &
  36              re_cloud_path_tot, re_ice_path_tot, tau_qc, tau_qi, tau_qs, tau_qc_tot, tau_qi_tot, &
  37              cbase, ctop, cbase_tot, ctop_tot, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, &
  38              kms, kme, ips, ipe, jps, jpe, kps, kpe, i_start, i_end, j_start, j_end, kts, kte, num_tiles)
  39
  40 !----------------------------------------------------------------------
  41
  42      USE module_model_constants, ONLY: G
  43
  44      IMPLICIT NONE
  45
  46      REAL, INTENT(IN) :: dt
  47      REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: coszen, swdnb, swdnt
  48      REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: ph, phb, cldfrac3d, qv, qc, qi, qs, qc_tot, qi_tot, &
  49           re_cloud, re_ice, re_snow, rho, dz8w
  50      INTEGER, INTENT(IN) :: param_first_scalar, p_qc, p_qi, p_qs, has_reqc, has_reqi, has_reqs
  51      LOGICAL, INTENT(IN) :: f_qv, f_qc, f_qi, f_qs
  52      REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: clrnidx, sza, cldfrac2d, wvp2d, lwp2d, iwp2d, swp2d, wp2d_sum, &
  53           lwp2d_tot, iwp2d_tot, wp2d_tot_sum, re_cloud_path, re_ice_path, re_snow_path, re_cloud_path_tot, re_ice_path_tot, &
  54           tau_qc, tau_qi, tau_qs, tau_qc_tot, tau_qi_tot, cbase, ctop, cbase_tot, ctop_tot
  55      REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: ghi_accum
  56      INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, kms, kme, &
  57              ips, ipe, jps, jpe, kps, kpe, kts, kte, num_tiles
  58      INTEGER, DIMENSION(num_tiles), INTENT(IN) :: i_start, i_end, j_start, j_end
  59
  60      INTEGER :: i, j, k, its, ite, jts, jte, ij
  61      REAL    :: pmw, swp, iwp, lwp, q_aux, wc
  62      REAL, PARAMETER :: THREE_OVER_TWO = 3.0 / 2.0
  63      REAL, PARAMETER :: PI = 4.0 * ATAN(1.0)
  64      REAL, PARAMETER :: DEGRAD = 180.0 / PI
  65      REAL, PARAMETER :: Q_MIN = 0.0
  66      REAL, PARAMETER :: WC_MIN = 1E-5
  67      REAL, PARAMETER :: RE_MIN = 0.0
  68      REAL, PARAMETER :: MISSING = -999.0
  69      REAL, DIMENSION(kms:kme) :: rhodz
  70      INTEGER, PARAMETER :: TAU_ICE_METHOD = 1
  71      INTEGER, PARAMETER :: I_TO_PLOT = 299, J_TO_PLOT = 84
  72
  73 !      !$OMP PARALLEL DO   &
  74 !      !$OMP PRIVATE ( ij )
  75        DO ij = 1, num_tiles
  76          its = i_start(ij)
  77          ite = i_end(ij)
  78          jts = j_start(ij)
  79          jte = j_end(ij)
  80          DO j = jts, jte
  81            DO i = its, ite
  82
  83              !!! SOLAR ZENITH ANGLE
  84              sza(i, j) = DEGRAD * MAX( ACOS( coszen(i, j) ), 0.0 )
  85
  86              !!! CLEARNESS INDEX
  87              if ( (sza(i, j) < 90.0) .and. (swdnt(i, j) > 0.0) ) then
  88                clrnidx(i, j) = swdnb(i, j) / swdnt(i, j)
  89              else
  90                clrnidx(i, j) = MISSING
  91              end if
  92
  93              !!! ACCUMULATED GHI
  94              ghi_accum(i, j) = ghi_accum(i, j) + swdnb(i, j) * dt
  95
  96              !!! 2-D CLOUD FRACTION
  97              cldfrac2d(i, j) = MAXVAL( cldfrac3d(i, kts:kte-1, j) )
  98
  99              !!! PREPARE FOR VARIABLE INTEGRATIONS
 100              rhodz(:) = rho(i, :, j) * dz8w(i, :, j) / (1. + qv(i, :, j))
 101
 102              !!! WATER VAPOR
 103              if (f_qv) then
 104                  ! Calc water vapor path
 105                pmw = integrate_1var (rhodz, qv(i, :, j), kms, kme, kts, kte)
 106                wvp2d(i, j) = pmw
 107              end if
 108
 109              !!! LIQUID WATER
 110              if (f_qc .and. p_qc > param_first_scalar) then
 111                  !!! RESOLVED !!!
 112                  ! Calc liquid water path
 113                q_aux = integrate_1var (rhodz, qc(i, :, j), kms, kme, kts, kte)
 114                lwp = q_aux
 115                lwp2d(i, j) = SIGN( MAX( lwp, 0.0 ), 1.0 )
 116
 117                if (has_reqc == 1) then
 118                    ! Calc effective radius water
 119                  if (q_aux > Q_MIN) then
 120                    re_cloud_path(i, j) = integrate_2var (rhodz, qc(i, :, j), &
 121                        re_cloud(i, :, j), kms, kme, kts, kte)
 122                    re_cloud_path(i, j) = re_cloud_path(i, j) / q_aux
 123                  else
 124                    re_cloud_path(i, j) = 0.0
 125                  end if
 126
 127                    ! Calc optical thickness water
 128                  if (re_cloud_path(i, j) > RE_MIN) then
 129                    tau_qc(i, j) = THREE_OVER_TWO * lwp / re_cloud_path(i, j) / 1000.0
 130                  else
 131                    tau_qc(i, j) = 0.0
 132                  end if
 133                end if
 134
 135                  !!! TOTAL (RESOLVED + UNRESOLVED) !!!
 136                  ! Calc liquid water path
 137                q_aux = integrate_1var (rhodz, qc_tot(i, :, j), kms, kme, kts, kte)
 138                lwp = q_aux
 139                lwp2d_tot(i, j) = SIGN( MAX( lwp, 0.0 ), 1.0 )
 140
 141                if (has_reqc == 1) then
 142                    ! Calc effective radius water
 143                  if (q_aux > Q_MIN) then
 144                    re_cloud_path_tot(i, j) = integrate_2var (rhodz, qc_tot(i, :, j), &
 145                        re_cloud(i, :, j), kms, kme, kts, kte)
 146                    re_cloud_path_tot(i, j) = re_cloud_path_tot(i, j) / q_aux
 147                  else
 148                    re_cloud_path_tot(i, j) = 0.0
 149                  end if
 150
 151                    ! Calc optical thickness water
 152                  if (re_cloud_path_tot(i, j) > RE_MIN) then
 153                    tau_qc_tot(i, j) = THREE_OVER_TWO * lwp / re_cloud_path_tot(i, j) / 1000.0
 154                  else
 155                    tau_qc_tot(i, j) = 0.0
 156                  end if
 157                end if
 158              else
 159                lwp2d(i, j) = MISSING
 160                re_cloud_path(i, j) = MISSING
 161                tau_qc(i, j) = MISSING
 162                lwp2d_tot(i, j) = MISSING
 163                re_cloud_path_tot(i, j) = MISSING
 164                tau_qc_tot(i, j) = MISSING
 165              end if
 166
 167              !!! ICE
 168              if (f_qi .and. p_qi > param_first_scalar) then
 169                  !!! RESOLVED !!!
 170                  ! Calc ice water path
 171                q_aux = integrate_1var (rhodz, qi(i, :, j), kms, kme, kts, kte)
 172                iwp = q_aux
 173                iwp2d(i, j) = SIGN( MAX( iwp, 0.0 ), 1.0 )
 174
 175                if (has_reqi == 1) then
 176                    ! Calc effective radius ice
 177                  if (q_aux > Q_MIN) then
 178                    re_ice_path(i, j) = integrate_2var (rhodz, qi(i, :, j), &
 179                        re_ice(i, :, j), kms, kme, kts, kte)
 180                    re_ice_path(i, j) = re_ice_path(i, j) / q_aux
 181                  else
 182                    re_ice_path(i, j) = 0.0
 183                  end if
 184
 185                    ! Calc optical thickness ice
 186                  if (re_ice_path(i, j) > RE_MIN) then
 187                    if (TAU_ICE_METHOD == 1) then
 188                        ! Eq 10 in Matrosov et al. (2002)
 189                      tau_qi(i, j) = iwp * 1000.0 * (0.02 + 4.2 / (2.0 * re_ice_path(i, j) * 1.0e+6))
 190                    else
 191                      tau_qi(i, j) = iwp * 1000.0 * (-0.006656 + 3.686 / (2.0 * re_ice_path(i, j) * 1.0e+6))
 192                    end if
 193                  else
 194                    tau_qi(i, j) = 0.0
 195                  end if
 196                end if
 197
 198                  !!! TOTAL (RESOLVED + UNRESOLVED) !!!
 199                  ! Calc ice water path
 200                q_aux = integrate_1var (rhodz, qi_tot(i, :, j), kms, kme, kts, kte)
 201                iwp = q_aux
 202                iwp2d_tot(i, j) = SIGN( MAX( iwp, 0.0 ), 1.0 )
 203
 204                if (has_reqi == 1) then
 205                    ! Calc effective radius ice
 206                  if (q_aux > Q_MIN) then
 207                    re_ice_path_tot(i, j) = integrate_2var (rhodz, qi_tot(i, :, j), &
 208                        re_ice(i, :, j), kms, kme, kts, kte)
 209                    re_ice_path_tot(i, j) = re_ice_path_tot(i, j) / q_aux
 210                  else
 211                    re_ice_path_tot(i, j) = 0.0
 212                  end if
 213
 214                    ! Calc optical thickness ice
 215                  if (re_ice_path_tot(i, j) > RE_MIN) then
 216                    if (TAU_ICE_METHOD == 1) then
 217                        ! Eq 10 in Matrosov et al. (2002)
 218                      tau_qi_tot(i, j) = iwp * 1000.0 * (0.02 + 4.2 / (2.0 * re_ice_path_tot(i, j) * 1.0e+6))
 219                    else
 220                      tau_qi_tot(i, j) = iwp * 1000.0 * (-0.006656 + 3.686 / (2.0 * re_ice_path_tot(i, j) * 1.0e+6))
 221                    end if
 222                  else
 223                    tau_qi_tot(i, j) = 0.0
 224                  end if
 225                end if
 226              else
 227                iwp2d(i, j) = MISSING
 228                re_ice_path(i, j) = MISSING
 229                tau_qi(i, j) = MISSING
 230                iwp2d_tot(i, j) = MISSING
 231                re_ice_path_tot(i, j) = MISSING
 232                tau_qi_tot(i, j) = MISSING
 233              end if
 234
 235              !!! SNOW
 236              if (f_qs .and. p_qs > param_first_scalar) then
 237                  !!! RESOLVED !!!
 238                  ! Calc effective radius snow
 239                q_aux = integrate_1var (rhodz, qs(i, :, j), kms, kme, kts, kte)
 240                swp = q_aux
 241                swp2d(i, j) = SIGN( MAX( swp, 0.0 ), 1.0 )
 242
 243                if (has_reqs == 1) then
 244                  if (q_aux > Q_MIN) then
 245                    re_snow_path(i, j) = integrate_2var (rhodz, qs(i, :, j), &
 246                        re_snow(i, :, j), kms, kme, kts, kte)
 247                    re_snow_path(i, j) = re_snow_path(i, j) / q_aux
 248                  else
 249                    re_snow_path(i, j) = 0.0
 250                  end if
 251
 252                    ! Optical thickness snow
 253                  if (re_snow_path(i, j) > RE_MIN) then
 254                    if (TAU_ICE_METHOD == 1) then
 255                        ! Eq 10 in Matrosov et al. (2002)
 256                      tau_qs(i, j) = swp * 1000.0 * (0.02 + 4.2 / (2.0 * re_snow_path(i, j) * 1.0e+6))
 257                    else
 258                      tau_qs(i, j) = swp * 1000.0 * (-0.006656 + 3.686 / (2.0 * re_snow_path(i, j) * 1.0e+6))
 259                    end if
 260                  else
 261                    tau_qs(i, j) = 0.0
 262                  end if
 263                end if
 264              else
 265                swp2d(i, j) = MISSING
 266                re_snow_path(i, j) = MISSING
 267                tau_qs(i, j) = MISSING
 268              end if
 269
 270              if ( (f_qc .or. f_qi .or. f_qs) .and. &
 271                   (p_qc > param_first_scalar .or. &
 272                    p_qi > param_first_scalar .or. &
 273                    p_qs > param_first_scalar) ) then
 274                  !!! RESOLVED !!!
 275                  ! Sum water paths (cloud liquid + ice + snow)
 276                wp2d_sum(i, j) = MAX( lwp2d(i, j), 0.0 ) + MAX( iwp2d(i, j), 0.0 ) + MAX( swp2d(i, j), 0.0 )
 277
 278                  !!! CLOUD BASE & TOP HEIGHTS
 279                  ! Cloud base first
 280                if (wp2d_sum(i, j) > Q_MIN) then
 281                  k = kts
 282                  wc = 0.0
 283                  if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc(i, k, j)
 284                  if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi(i, k, j)
 285                  if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j)
 286                  do while ( (wc < WC_MIN) .and. (k < kte) )
 287                    k = k + 1
 288                    wc = 0.0
 289                    if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc(i, k, j)
 290                    if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi(i, k, j)
 291                    if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j)
 292                  end do
 293
 294                  if (k == kte) then
 295                    cbase(i, j) = MISSING
 296                  else
 297                    cbase(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G)
 298                  end if
 299
 300                    ! Cloud top second
 301                  k = kte
 302                  wc = 0.0
 303                  if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc(i, k, j)
 304                  if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi(i, k, j)
 305                  if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j)
 306                  do while ( (wc < WC_MIN) .and. (k > kts) )
 307                    k = k - 1
 308                    wc = 0.0
 309                    if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc(i, k, j)
 310                    if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi(i, k, j)
 311                    if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j)
 312                  end do
 313
 314                  if (k == kts) then
 315                    ctop(i, j) = MISSING
 316                  else
 317                    ctop(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G)
 318                  end if
 319                else
 320                  cbase(i, j) = MISSING
 321                  ctop(i, j) = MISSING
 322                end if
 323
 324                  !!! TOTAL (RESOLVED + UNRESOLVED) !!!
 325                  ! Sum water paths (cloud liquid + ice + snow)
 326                  ! Note that snow is from resolved only
 327                wp2d_tot_sum(i, j) = MAX( lwp2d_tot(i, j), 0.0 ) + MAX( iwp2d_tot(i, j), 0.0 ) + MAX( swp2d(i, j), 0.0 )
 328
 329                  ! Cloud base first
 330                if (wp2d_tot_sum(i, j) > Q_MIN) then
 331                  k = kts
 332                  wc = 0.0
 333                  if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc_tot(i, k, j)
 334                  if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi_tot(i, k, j)
 335                  if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j)
 336                  do while ( (wc < WC_MIN) .and. (k < kte) )
 337                    k = k + 1
 338                    wc = 0.0
 339                    if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc_tot(i, k, j)
 340                    if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi_tot(i, k, j)
 341                    if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j)
 342                  end do
 343
 344                  if (k == kte) then
 345                    cbase_tot(i, j) = MISSING
 346                  else
 347                    cbase_tot(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G)
 348                  end if
 349
 350                    ! Cloud top second
 351                  k = kte
 352                  wc = 0.0
 353                  if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc_tot(i, k, j)
 354                  if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi_tot(i, k, j)
 355                  if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j)
 356                  do while ( (wc < WC_MIN) .and. (k > kts) )
 357                    k = k - 1
 358                    wc = 0.0
 359                    if (f_qc .and. p_qc > param_first_scalar) wc = wc + qc_tot(i, k, j)
 360                    if (f_qi .and. p_qi > param_first_scalar) wc = wc + qi_tot(i, k, j)
 361                    if (f_qs .and. p_qs > param_first_scalar) wc = wc + qs(i, k, j)
 362                  end do
 363
 364                  if (k == kts) then
 365                    ctop_tot(i, j) = MISSING
 366                  else
 367                    ctop_tot(i, j) = ( (ph(i, k, j) + phb(i, k, j)) + (ph(i, k + 1, j) + phb(i, k + 1, j)) ) / (2.0 * G)
 368                  end if
 369                else
 370                  cbase_tot(i, j) = MISSING
 371                  ctop_tot(i, j) = MISSING
 372                end if
 373              else
 374                wp2d_sum(i, j) = MISSING
 375                cbase(i, j) = MISSING
 376                ctop(i, j) = MISSING
 377                wp2d_tot_sum(i, j) = MISSING
 378                cbase_tot(i, j) = MISSING
 379                ctop_tot(i, j) = MISSING
 380              end if
 381
 382            ENDDO
 383          ENDDO
 384        ENDDO
 385 !      !$OMP END PARALLEL DO
 386
 387    END SUBROUTINE solar_diag
 388
 389
 390    FUNCTION Integrate_1var (rhodz, var1_1d, kms, kme, kts, kte) &
 391         RESULT (return_value)
 392
 393       IMPLICIT NONE
 394
 395       INTEGER, INTENT(IN) :: kts, kte, kms, kme
 396       REAL, DIMENSION(kms:kme), INTENT(IN) :: var1_1d, rhodz
 397
 398         ! Local
 399       REAL :: return_value
 400       INTEGER :: k
 401
 402       return_value = 0.0
 403       do k = kts, kte - 1
 404         return_value = return_value + var1_1d(k) * rhodz(k)
 405       end do
 406
 407     END FUNCTION Integrate_1var
 408
 409
 410     FUNCTION Integrate_2var (rhodz, var1_1d, var2_1d, kms, kme, kts, kte) &
 411         RESULT (return_value)
 412
 413       IMPLICIT NONE
 414
 415       INTEGER, INTENT(IN) :: kts, kte, kms, kme
 416       REAL, DIMENSION(kms:kme), INTENT(IN) :: var1_1d, var2_1d, rhodz
 417
 418         ! Local
 419       REAL :: return_value
 420       INTEGER :: k
 421
 422       return_value = 0.0
 423       do k = kts, kte - 1
 424         return_value = return_value + var1_1d(k) * var2_1d(k) * rhodz(k)
 425       end do
 426
 427     END FUNCTION Integrate_2var
 428
 429 END MODULE module_diag_solar
 430 #endif