phys/module_diag_pld.F

   1 #if (NMM_CORE == 1)
   2 MODULE module_diag_pld
   3 CONTAINS
   4    SUBROUTINE diag_pld_stub
   5    END SUBROUTINE diag_pld_stub
   6 END MODULE module_diag_pld
   7 #else
   8 !WRF:MEDIATION_LAYER:PHYSICS
   9 !
  10
  11 MODULE module_diag_pld
  12 CONTAINS
  13
  14    SUBROUTINE pld ( u,v,w,t,qv,zp,zb,pp,pb,p,pw,                    &
  15                     msfux,msfuy,msfvx,msfvy,msftx,msfty,            &
  16                     f,e,                                            &
  17                     use_tot_or_hyd_p,extrap_below_grnd,missing,     &
  18                     num_press_levels,max_press_levels,press_levels, &
  19                     p_pl,u_pl,v_pl,t_pl,rh_pl,ght_pl,s_pl,td_pl,    &
  20                     q_pl,                                           &
  21                     ids,ide, jds,jde, kds,kde,                      &
  22                     ims,ime, jms,jme, kms,kme,                      &
  23                     its,ite, jts,jte, kts,kte                       )
  24
  25       USE module_model_constants
  26
  27       IMPLICIT NONE
  28
  29
  30       !  Input variables
  31
  32       INTEGER, INTENT(IN   )                                          :: ids,ide, jds,jde, kds,kde, &
  33                                                                          ims,ime, jms,jme, kms,kme, &
  34                                                                          its,ite, jts,jte, kts,kte
  35       REAL   , INTENT(IN   ) , DIMENSION(ims:ime , jms:jme)           :: msfux,msfuy,msfvx,msfvy,msftx,msfty, &
  36                                                                          f,e
  37       INTEGER, INTENT(IN   )                                          :: use_tot_or_hyd_p
  38       INTEGER, INTENT(IN   )                                          :: extrap_below_grnd
  39       REAL   , INTENT(IN   )                                          :: missing
  40       REAL   , INTENT(IN   ) , DIMENSION(ims:ime , kms:kme , jms:jme) :: u,v,w,t,qv,zp,zb,pp,pb,p,pw
  41       INTEGER, INTENT(IN   )                                          :: num_press_levels, max_press_levels
  42       REAL   , INTENT(IN   ) , DIMENSION(max_press_levels)            :: press_levels
  43
  44       !  Output variables
  45
  46       REAL   , INTENT(  OUT) ,  DIMENSION(num_press_levels)                     :: p_pl
  47       REAL   , INTENT(  OUT) ,  DIMENSION(ims:ime , num_press_levels , jms:jme) :: u_pl,v_pl,t_pl,rh_pl,ght_pl,s_pl,td_pl,q_pl
  48
  49       !  Local variables
  50
  51       REAL, PARAMETER :: eps = 0.622, t_kelvin = svpt0 , s1 = 243.5, s2 = svp2 , s3 = svp1*10., s4 = 611.0, s5 = 5418.12
  52       REAL, PARAMETER :: zshul=75., tvshul=290.66
  53
  54       INTEGER :: i, j, ke, kp, ke_h, ke_f
  55       REAL    :: pu, pd, pm , &
  56                  tu, td     , &
  57                  su, sd     , &
  58                  uu, ud     , &
  59                  vu, vd     , &
  60                  zu, zd     , &
  61                  qu, qd     , &
  62                  eu, ed, em , &
  63                  du, dd
  64       REAL    :: es, qs
  65       REAL    :: part, gammas, tvu, tvd
  66
  67       !  Silly, but transfer the small namelist.input array into the grid structure for output purposes.
  68
  69       DO kp = 1 , num_press_levels
  70          p_pl(kp) = press_levels(kp)
  71       END DO
  72
  73       !  Initialize pressure level data to un-initialized
  74
  75       DO j = jts , jte
  76          DO kp = 1 , num_press_levels
  77             DO i = its , ite
  78                u_pl  (i,kp,j) = missing
  79                v_pl  (i,kp,j) = missing
  80                t_pl  (i,kp,j) = missing
  81                rh_pl (i,kp,j) = missing
  82                ght_pl(i,kp,j) = missing
  83                s_pl  (i,kp,j) = missing
  84                td_pl (i,kp,j) = missing
  85                q_pl  (i,kp,j) = missing
  86             END DO
  87          END DO
  88       END DO
  89
  90       !  Loop over each i,j location
  91
  92       j_loop : DO j = jts , MIN(jte,jde-1)
  93          i_loop : DO i = its , MIN(ite,ide-1)
  94
  95             !  For each i,j location, loop over the selected
  96             !  pressure levels to find
  97
  98             ke_h = kts
  99             ke_f = kts
 100             kp_loop : DO kp = 1 , num_press_levels
 101
 102                !  For this particular i,j and pressure level, find the
 103                !  eta levels that surround this point on half-levels.
 104
 105                ke_loop_half : DO ke = ke_h , kte-2
 106
 107                   IF      ( use_tot_or_hyd_p .EQ. 1 ) THEN     !  total pressure
 108                      pu = pp(i,ke+1,j)+pb(i,ke+1,j)
 109                      pd = pp(i,ke  ,j)+pb(i,ke  ,j)
 110                   ELSE IF ( use_tot_or_hyd_p .EQ. 2 ) THEN     !  hydrostatic pressure
 111                      pu = p(i,ke+1,j)
 112                      pd = p(i,ke  ,j)
 113                   END IF
 114                   pm = p_pl(kp)
 115
 116                   !  Added option to extrapolate below ground - GAC (AFWA)
 117
 118                   IF ( ( extrap_below_grnd .EQ. 2 ) .AND.  &
 119                      ( ke .EQ. ke_h ) .AND. ( pm .GT. pd )) THEN
 120
 121                      !  Requested pressure level is below ground.
 122                      !  Extrapolate adiabatically if requested in namelist.
 123
 124                      !  Methodology derived from Unified Post Processor (UPP).
 125                      !  Simply conserve first level U, V, and RH below ground.
 126                      !  Assume adiabatic lapse rate of gamma = 6.5 K/km
 127                      !  below ground, using Shuell correction to gamma
 128                      !  ("gammas") to find geopotential height, which is
 129                      !  computed by hydrostatically integrating mean isobaric
 130                      !  virtual temperature downward from the model surface.
 131                      !  Temperature is found by reducing adiabatically
 132                      !  from the first level temperature.
 133                      !  Sources:
 134                      !    Chuang et al, NCEP's WRF Post Processor and
 135                      !      Verification Systems, MM5 Workshop Session 7, 2004.
 136                      !    Unipost source code: MDL2P.f
 137
 138                      !  Z, T, Q, Tv at first half-eta level
 139
 140                      zu = 0.5 * ( zp(i,ke  ,j) + zb(i,ke  ,j) + &
 141                                   zp(i,ke+1,j) + zb(i,ke+1,j) ) / g
 142                      tu = ( t(i,ke,j) + t0 ) * ( pd / p1000mb ) ** rcp
 143                      qu = MAX(qv(i,ke,j),0.)
 144                      tvu = tu * ( 1. + 0.608 * qu )
 145
 146                      !  1. Geopotential height (m)
 147
 148                      IF ( zu .GT. zshul ) THEN
 149                         tvd = tvu + zu * 6.5E-3
 150                         IF ( tvd .GT. tvshul ) THEN
 151                           IF ( tvu .GT. tvshul) THEN
 152                             tvd = tvshul - 5.E-3 * ( tvu - tvshul ) ** 2
 153                           ELSE
 154                             tvd = tvshul
 155                           ENDIF
 156                         ENDIF
 157                         gammas = ( tvu - tvd ) / zu
 158                      ELSE
 159                         gammas = 0.
 160                      ENDIF
 161                      part = ( r_d / g ) * ( ALOG (pm) - ALOG (pd) )
 162                      ght_pl(i,kp,j) = zu - tvu * part / &
 163                                       ( 1. + 0.5 * gammas * part )
 164
 165                      !  2. Temperature (K)
 166
 167                      t_pl(i,kp,j) = tu + ( zu - ght_pl(i,kp,j) ) * 6.5E-3
 168
 169                      !  3. Speed (m s-1)
 170
 171                      s_pl(i,kp,j) = 0.5 * SQRT ( ( u(i,ke  ,j)+ &
 172                                    u(i+1,ke  ,j) )**2 +         &
 173                                    ( v(i,ke  ,j) + v(i,ke  ,j+1) )**2 )
 174
 175                      !  4. U and V (m s-1)
 176
 177                      u_pl(i,kp,j) = 0.5 * ( u(i,ke  ,j) + u(i+1,ke  ,j) )
 178                      v_pl(i,kp,j) = 0.5 * ( v(i,ke  ,j) + v(i,ke  ,j+1) )
 179
 180                      !  5. Relative humidity (%)
 181
 182                      es = s4 * exp(s5 * (1.0 / 273.0 - 1.0 / tu) )
 183                      qs = eps * es / (pd - es)
 184                      rh_pl(i,kp,j)   = MAX(qv(i,ke,j),0.) / qs * 100.
 185
 186                      !  6. Mixing ratio (kg/kg)
 187
 188                      es = s4 * exp(s5 * (1.0 / 273.0 - 1.0 / t_pl(i,kp,j)))
 189                      qs = eps * es / (pm - es)
 190                      q_pl(i,kp,j)   = rh_pl(i,kp,j) * qs / 100.
 191
 192                      !  7. Dewpoint (K) - Use Bolton's approximation
 193
 194                      ed = q_pl(i,kp,j) * pm * 0.01 / ( eps + q_pl(i,kp,j) )
 195                      ed = max(ed, 0.001)   ! water vapor pressure in mb.
 196                      td_pl(i,kp,j) = t_kelvin + (s1 / ((s2 / log(ed/s3)) - 1.0))
 197
 198                      EXIT ke_loop_half
 199                   ELSEIF ( ( pd .GE. pm ) .AND. &
 200                        ( pu .LT. pm ) ) THEN
 201
 202                      !  Found trapping pressure: up, middle, down.
 203                      !  We are doing first order interpolation.
 204                      !  Now we just put in a list of diagnostics for this level.
 205
 206                      !  1. Temperature (K)
 207
 208                      tu = (t(i,ke+1,j)+t0)*(pu/p1000mb)**rcp
 209                      td = (t(i,ke  ,j)+t0)*(pd/p1000mb)**rcp
 210                      t_pl(i,kp,j) = ( tu * (pm-pd) + td * (pu-pm) ) / (pu-pd)
 211
 212                      !  2. Speed (m s-1)
 213
 214                      su = 0.5 * SQRT ( ( u(i,ke+1,j)+u(i+1,ke+1,j) )**2 + &
 215                                        ( v(i,ke+1,j)+v(i,ke+1,j+1) )**2 )
 216                      sd = 0.5 * SQRT ( ( u(i,ke  ,j)+u(i+1,ke  ,j) )**2 + &
 217                                        ( v(i,ke  ,j)+v(i,ke  ,j+1) )**2 )
 218                      s_pl(i,kp,j) = ( su * (pm-pd) + sd * (pu-pm) ) / (pu-pd)
 219
 220                      !  3. U and V (m s-1)
 221
 222                      uu = 0.5 * ( u(i,ke+1,j)+u(i+1,ke+1,j) )
 223                      ud = 0.5 * ( u(i,ke  ,j)+u(i+1,ke  ,j) )
 224                      u_pl(i,kp,j) = ( uu * (pm-pd) + ud * (pu-pm) ) / (pu-pd)
 225
 226                      vu = 0.5 * ( v(i,ke+1,j)+v(i,ke+1,j+1) )
 227                      vd = 0.5 * ( v(i,ke  ,j)+v(i,ke  ,j+1) )
 228                      v_pl(i,kp,j) = ( vu * (pm-pd) + vd * (pu-pm) ) / (pu-pd)
 229
 230                      !  4. Mixing ratio (kg/kg)
 231
 232                      qu = MAX(qv(i,ke+1,j),0.)
 233                      qd = MAX(qv(i,ke  ,j),0.)
 234                      q_pl(i,kp,j) = ( qu * (pm-pd) + qd * (pu-pm) ) / (pu-pd)
 235
 236                      !  5. Dewpoint (K) - Use Bolton's approximation
 237
 238                      eu = qu * pu * 0.01 / ( eps + qu ) ! water vapor press (mb)
 239                      ed = qd * pd * 0.01 / ( eps + qd ) ! water vapor press (mb)
 240                      eu = max(eu, 0.001)
 241                      ed = max(ed, 0.001)
 242
 243                      du = t_kelvin + ( s1 / ((s2 / log(eu/s3)) - 1.0) )
 244                      dd = t_kelvin + ( s1 / ((s2 / log(ed/s3)) - 1.0) )
 245                      td_pl(i,kp,j) = ( du * (pm-pd) + dd * (pu-pm) ) / (pu-pd)
 246
 247
 248                      !  6. Relative humidity (%)
 249
 250                      es = s4 * exp(s5 * (1.0 / 273.0 - 1.0 / t_pl(i,kp,j)))
 251                      qs = eps * es / (pm - es)
 252                      rh_pl(i,kp,j)   = q_pl(i,kp,j) / qs * 100.
 253
 254                      !em = qm * pm * 0.01 / ( eps + qm )                                       ! water vapor pressure at the level.
 255                      !es = s3 * exp( s2 * (t_pl(i,kp,j) - t_kelvin)/(t_pl(i,kp,j) - s4) )      ! sat vapor pressure over liquid water in mb.
 256                      !rh_pl(i,kp,j) = 100. * em * ( pm * 0.01 - es ) / ( es * ( pm * 0.01 - em ) )
 257
 258                      ke_h = ke
 259                      EXIT ke_loop_half
 260                   END IF
 261                END DO ke_loop_half
 262
 263                ke_loop_full : DO ke = ke_f , kte-1
 264
 265                   IF ( ( pw(i,ke  ,j) .GE. p_pl(kp) ) .AND. &
 266                        ( pw(i,ke+1,j) .LT. p_pl(kp) ) ) THEN
 267
 268                      !  Found trapping pressure: up, middle, down.
 269                      !  We are doing first order interpolation.
 270
 271                      pu = LOG(pw(i,ke+1,j))
 272                      pm = LOG(p_pl(kp))
 273                      pd = LOG(pw(i,ke  ,j))
 274
 275                      !  Now we just put in a list of diagnostics for this level.
 276
 277                      !  1. Geopotential height (m)
 278
 279                      zu = ( zp(i,ke+1,j)+zb(i,ke+1,j) ) / g
 280                      zd = ( zp(i,ke  ,j)+zb(i,ke  ,j) ) / g
 281                      ght_pl(i,kp,j) = ( zu * (pm-pd) + zd * (pu-pm) ) / (pu-pd)
 282
 283                      ke_f = ke
 284                      EXIT ke_loop_full
 285                   END IF
 286                END DO ke_loop_full
 287
 288             END DO kp_loop
 289          END DO i_loop
 290       END DO j_loop
 291
 292    END SUBROUTINE pld
 293
 294 END MODULE module_diag_pld
 295 #endif