phys/module_diag_zld.F

   1 #if (NMM_CORE == 1)
   2 MODULE module_diag_zld
   3 CONTAINS
   4    SUBROUTINE diag_zld_stub
   5    END SUBROUTINE diag_zld_stub
   6 END MODULE module_diag_zld
   7 #else
   8 !WRF:MEDIATION_LAYER:PHYSICS
   9 !
  10
  11 MODULE module_diag_zld
  12 CONTAINS
  13
  14    SUBROUTINE zld ( u,v,w,t,qv,zp,zb,pp,pb,p,pw,                    &
  15                     msfux,msfuy,msfvx,msfvy,msftx,msfty,            &
  16                     f,e,ht,                                         &
  17                     use_tot_or_hyd_p,extrap_below_grnd,missing,     &
  18                     num_z_levels,max_z_levels,z_levels,             &
  19                     z_zl,u_zl,v_zl,t_zl,rh_zl,ght_zl,s_zl,td_zl,    &
  20                     q_zl,p_zl,                                      &
  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,ht
  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_z_levels, max_z_levels
  42       REAL   , INTENT(IN   ) , DIMENSION(max_z_levels)                :: z_levels
  43
  44       !  Output variables
  45
  46       REAL   , INTENT(  OUT) ,  DIMENSION(num_z_levels)                     :: z_zl
  47       REAL   , INTENT(  OUT) ,  DIMENSION(ims:ime , num_z_levels , jms:jme) :: u_zl,v_zl,t_zl,rh_zl,ght_zl,s_zl,td_zl,q_zl,p_zl
  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, kz, ke_h, ke_f
  55       REAL    :: zu, zd, zm , &
  56                  tu, td     , &
  57                  su, sd     , &
  58                  uu, ud     , &
  59                  vu, vd     , &
  60                  qu, qd     , &
  61                  eu, ed, em , &
  62                  pu, pd, pm , &
  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 kz = 1 , num_z_levels
  70          z_zl(kz) = z_levels(kz)
  71       END DO
  72
  73       !  Initialize height level data to un-initialized
  74
  75       DO j = jts , jte
  76          DO kz = 1 , num_z_levels
  77             DO i = its , ite
  78                u_zl  (i,kz,j) = missing
  79                v_zl  (i,kz,j) = missing
  80                t_zl  (i,kz,j) = missing
  81                rh_zl (i,kz,j) = missing
  82                ght_zl(i,kz,j) = missing
  83                s_zl  (i,kz,j) = missing
  84                td_zl (i,kz,j) = missing
  85                q_zl  (i,kz,j) = missing
  86                p_zl  (i,kz,j) = missing
  87             END DO
  88          END DO
  89       END DO
  90
  91       !  Loop over each i,j location
  92
  93       j_loop : DO j = jts , MIN(jte,jde-1)
  94          i_loop : DO i = its , MIN(ite,ide-1)
  95
  96             !  For each i,j location, loop over the selected
  97             !  pressure levels to find
  98
  99             ke_h = kts
 100             ke_f = kts
 101             kz_loop : DO kz = 1 , num_z_levels
 102
 103                !  For this particular i,j and height level, find the
 104                !  eta levels that surround this point on half-levels.
 105                !  Negative heights are a flag to do AGL.
 106
 107                ke_loop_half : DO ke = ke_h , kte-2
 108
 109                   zu = ( zp(i,ke+1,j)+zb(i,ke+1,j) + zp(i,ke+2,j)+zb(i,ke+2,j) ) / 2.0 / g
 110                   zd = ( zp(i,ke  ,j)+zb(i,ke  ,j) + zp(i,ke+1,j)+zb(i,ke+1,j) ) / 2.0 / g
 111                   IF ( z_zl(kz) .LT. 1 ) THEN
 112                      zm = ABS(z_zl(kz)) + ht(i,j)
 113                   ELSE
 114                      zm = z_zl(kz)
 115                   END IF
 116
 117                   IF ( ( zd .LE. zm ) .AND. ( zu .GT. zm ) ) THEN
 118
 119                      pu = pp(i,ke+1,j)+pb(i,ke+1,j)
 120                      pd = pp(i,ke  ,j)+pb(i,ke  ,j)
 121
 122                      !  Found trapping height: up, middle, down.
 123                      !  We are doing first order interpolation.
 124                      !  Now we just put in a list of diagnostics for this level.
 125
 126                      !  0. Pressure (Pa)
 127                      !  Note that it is ln(p) that varies linearly with height, not p itself
 128
 129                      pm = exp( ( log(pu) * (zm-zd) + log(pd) * (zu-zm) ) / (zu-zd) )
 130                      p_zl(i,kz,j) = pm
 131
 132                      !  1. Temperature (K)
 133
 134                      tu = (t(i,ke+1,j)+t0)*(pu/p1000mb)**rcp
 135                      td = (t(i,ke  ,j)+t0)*(pd/p1000mb)**rcp
 136                      t_zl(i,kz,j) = ( tu * (zm-zd) + td * (zu-zm) ) / (zu-zd)
 137
 138                      !  2. Speed (m s-1)
 139
 140                      su = 0.5 * SQRT ( ( u(i,ke+1,j)+u(i+1,ke+1,j) )**2 + &
 141                                        ( v(i,ke+1,j)+v(i,ke+1,j+1) )**2 )
 142                      sd = 0.5 * SQRT ( ( u(i,ke  ,j)+u(i+1,ke  ,j) )**2 + &
 143                                        ( v(i,ke  ,j)+v(i,ke  ,j+1) )**2 )
 144                      s_zl(i,kz,j) = ( su * (zm-zd) + sd * (zu-zm) ) / (zu-zd)
 145
 146                      !  3. U and V (m s-1)
 147
 148                      uu = 0.5 * ( u(i,ke+1,j)+u(i+1,ke+1,j) )
 149                      ud = 0.5 * ( u(i,ke  ,j)+u(i+1,ke  ,j) )
 150                      u_zl(i,kz,j) = ( uu * (zm-zd) + ud * (zu-zm) ) / (zu-zd)
 151
 152                      vu = 0.5 * ( v(i,ke+1,j)+v(i,ke+1,j+1) )
 153                      vd = 0.5 * ( v(i,ke  ,j)+v(i,ke  ,j+1) )
 154                      v_zl(i,kz,j) = ( vu * (zm-zd) + vd * (zu-zm) ) / (zu-zd)
 155
 156                      !  4. Mixing ratio (kg/kg)
 157
 158                      qu = MAX(qv(i,ke+1,j),0.)
 159                      qd = MAX(qv(i,ke  ,j),0.)
 160                      q_zl(i,kz,j) = ( qu * (zm-zd) + qd * (zu-zm) ) / (zu-zd)
 161
 162                      !  5. Dewpoint (K) - Use Bolton's approximation
 163
 164                      eu = qu * pu * 0.01 / ( eps + qu ) ! water vapor press (mb)
 165                      ed = qd * pd * 0.01 / ( eps + qd ) ! water vapor press (mb)
 166                      eu = max(eu, 0.001)
 167                      ed = max(ed, 0.001)
 168
 169                      du = t_kelvin + ( s1 / ((s2 / log(eu/s3)) - 1.0) )
 170                      dd = t_kelvin + ( s1 / ((s2 / log(ed/s3)) - 1.0) )
 171                      td_zl(i,kz,j) = ( du * (zm-zd) + dd * (zu-zm) ) / (zu-zd)
 172
 173
 174                      !  6. Relative humidity (%)
 175
 176                      es = s4 * exp(s5 * (1.0 / 273.0 - 1.0 / t_zl(i,kz,j)))
 177                      qs = eps * es / (pm - es)
 178                      rh_zl(i,kz,j)   = q_zl(i,kz,j) / qs * 100.
 179
 180                      ke_h = ke
 181                      EXIT ke_loop_half
 182                   END IF
 183                END DO ke_loop_half
 184
 185                ke_loop_full : DO ke = ke_f , kte-1
 186
 187                   zu = ( zp(i,ke+1,j)+zb(i,ke+1,j) ) / g
 188                   zd = ( zp(i,ke  ,j)+zb(i,ke  ,j) ) / g
 189                   IF ( z_zl(kz) .LT. 1 ) THEN
 190                      zm = ABS(z_zl(kz)) + ht(i,j)
 191                   ELSE
 192                      zm = z_zl(kz)
 193                   END IF
 194
 195                   IF ( ( zd .LE. zm ) .AND. &
 196                        ( zu .GT. zm) ) THEN
 197
 198                      !  Found trapping height: up, middle, down.
 199                      !  We are doing first order interpolation.
 200
 201                      !  Now we just put in a list of diagnostics for this level.
 202
 203                      !  1. Geopotential height (m)
 204
 205                      ght_zl(i,kz,j) = ( zu * (zm-zd) + zd * (zu-zm) ) / (zu-zd)
 206
 207                      ke_f = ke
 208                      EXIT ke_loop_full
 209                   END IF
 210                END DO ke_loop_full
 211
 212             END DO kz_loop
 213          END DO i_loop
 214       END DO j_loop
 215
 216    END SUBROUTINE zld
 217
 218 END MODULE module_diag_zld
 219 #endif