chem/module_chem_utilities.F

   1 MODULE module_chem_utilities
   2    USE module_domain
   3    USE module_model_constants
   4    USE module_state_description
   5    USE module_configure
   6
   7 CONTAINS
   8    SUBROUTINE chem_prep ( config_flags,                                &  ! input
   9                          u, v, p, pb, alt, ph,                    &  ! input
  10                          phb, t, moist, n_moist,                 &  ! input
  11                          rho, p_phy ,          &  ! output
  12                          u_phy, v_phy, p8w, t_phy, t8w,               &  ! output
  13                          z, z_at_w, dz8w,rh,                          &  ! output
  14                          fzm, fzp,                                    &  ! params
  15                          ids, ide, jds, jde, kds, kde,                &
  16                          ims, ime, jms, jme, kms, kme,                &
  17                          its, ite, jts, jte, kts, kte                )
  18 !----------------------------------------------------------------------
  19    IMPLICIT NONE
  20 !----------------------------------------------------------------------
  21
  22    TYPE(grid_config_rec_type) ,     INTENT(IN   ) :: config_flags
  23    INTEGER ,        INTENT(IN   ) ::   ids, ide, jds, jde, kds, kde, &
  24                                        ims, ime, jms, jme, kms, kme, &
  25                                        its, ite, jts, jte, kts, kte
  26    INTEGER ,          INTENT(IN   ) :: n_moist
  27
  28    REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
  29
  30
  31
  32    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                 &
  33           INTENT(  OUT)                                  ::   u_phy, &
  34                                                               v_phy, &
  35                                                               p_phy, &
  36                                                                 p8w, &
  37                                                               t_phy, &
  38                                                                 t8w, &
  39                                                                 rho, &
  40                                                                   z, &
  41                                                                dz8w, &
  42                                                          rh,  z_at_w
  43
  44    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                 &
  45           INTENT(IN   )                                  ::      pb, &
  46                                                                   p, &
  47                                                                   u, &
  48                                                                   v, &
  49                                                                 alt, &
  50                                                                  ph, &
  51                                                                 phb, &
  52                                                                   t
  53
  54
  55    REAL , DIMENSION( kms:kme ) ,           INTENT(IN   ) ::     fzm,   &
  56                                                                 fzp
  57
  58    INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
  59    INTEGER :: i, j, k
  60    REAL    :: w1, w2, z0, z1, z2
  61
  62 !-----------------------------------------------------------------------
  63
  64 !  set up loop bounds for this grid's boundary conditions
  65
  66     i_start = its
  67     i_end   = min( ite,ide-1 )
  68     j_start = jts
  69     j_end   = min( jte,jde-1 )
  70
  71     k_start = kts
  72     k_end = min( kte, kde-1 )
  73
  74 !  compute thermodynamics and velocities at pressure points
  75     do j = j_start,j_end
  76     do k = k_start, k_end
  77     do i = i_start, i_end
  78
  79       p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
  80       t_phy(i,k,j) = (t(i,k,j)+t0)*(p_phy(i,k,j)/p1000mb)**rcp
  81       rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
  82       u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
  83       v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
  84
  85     enddo
  86     enddo
  87     enddo
  88
  89 !  wig: added to make sure there is no junk in the top level even
  90 !       though it should not be used
  91     do j = j_start,j_end
  92     do i = i_start, i_end
  93       p_phy(i,kte,j) = p_phy(i,k_end,j)
  94       t_phy(i,kte,j) = t_phy(i,k_end,j)
  95       rho(i,kte,j) = rho(i,k_end,j)
  96       u_phy(i,kte,j) = u_phy(i,k_end,j)
  97       v_phy(i,kte,j) = v_phy(i,k_end,j)
  98     enddo
  99     enddo
 100
 101 !  compute z at w points
 102
 103     do j = j_start,j_end
 104     do k = k_start, kte
 105     do i = i_start, i_end
 106       z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
 107     enddo
 108     enddo
 109     enddo
 110
 111     do j = j_start,j_end
 112     do k = k_start, kte-1
 113     do i = i_start, i_end
 114       dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
 115     enddo
 116     enddo
 117     enddo
 118
 119     do j = j_start,j_end
 120     do i = i_start, i_end
 121       dz8w(i,kte,j) = 0.
 122     enddo
 123     enddo
 124
 125 !  compute z at p points (average of z at w points)
 126     do j = j_start,j_end
 127     do k = k_start, k_end
 128     do i = i_start, i_end
 129       z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
 130       rh(i,k,j) = max(.1,MIN( .95, moist(i,k,j,p_qv) / &
 131                (3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
 132                (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j)))))
 133 !      rh(i,k,j)=max(.1,rh(i,k,j))
 134     enddo
 135     enddo
 136     enddo
 137
 138
 139 !  interp t and p at w points
 140
 141     do j = j_start,j_end
 142     do k = 2, k_end
 143     do i = i_start, i_end
 144       p8w(i,k,j) = fzm(k)*p_phy(i,k,j)+fzp(k)*p_phy(i,k-1,j)
 145       t8w(i,k,j) = fzm(k)*t_phy(i,k,j)+fzp(k)*t_phy(i,k-1,j)
 146     enddo
 147     enddo
 148     enddo
 149
 150 !  extrapolate p and t to surface and top.
 151 !  we'll use an extrapolation in z for now
 152
 153     do j = j_start,j_end
 154     do i = i_start, i_end
 155
 156 ! bottom
 157
 158       z0 = z_at_w(i,1,j)
 159       z1 = z(i,1,j)
 160       z2 = z(i,2,j)
 161       w1 = (z0 - z2)/(z1 - z2)
 162       w2 = 1. - w1
 163       p8w(i,1,j) = w1*p_phy(i,1,j)+w2*p_phy(i,2,j)
 164       t8w(i,1,j) = w1*t_phy(i,1,j)+w2*t_phy(i,2,j)
 165
 166 ! top
 167
 168       z0 = z_at_w(i,kte,j)
 169       z1 = z(i,k_end,j)
 170       z2 = z(i,k_end-1,j)
 171       w1 = (z0 - z2)/(z1 - z2)
 172       w2 = 1. - w1
 173
 174 !      p8w(i,kde,j) = w1*p_phy(i,kde-1,j)+w2*p_phy(i,kde-2,j)
 175 !!!  bug fix      extrapolate ln(p) so p is positive definite
 176       p8w(i,kde,j) = exp(w1*log(p_phy(i,kde-1,j))+w2*log(p_phy(i,kde-2,j)))
 177       t8w(i,kde,j) = w1*t_phy(i,kde-1,j)+w2*t_phy(i,kde-2,j)
 178
 179     enddo
 180     enddo
 181 END SUBROUTINE chem_prep
 182
 183    subroutine UPCASE( lstring )
 184 !----------------------------------------------------------------------
 185 !       ... Convert character string lstring to upper case
 186 !----------------------------------------------------------------------
 187    implicit none
 188
 189 !-----------------------------------------------------------------------
 190 !       ... Dummy args
 191 !-----------------------------------------------------------------------
 192    character(len=*), intent(inout) ::  lstring
 193
 194 !-----------------------------------------------------------------------
 195 !       ... Local variables
 196 !-----------------------------------------------------------------------
 197    integer :: i
 198
 199    do i = 1,LEN_TRIM( lstring )
 200      if( ICHAR(lstring(i:i)) >= 97 .and.  ICHAR(lstring(i:i)) <= 122 ) then
 201        lstring(i:i) = CHAR(ICHAR(lstring(i:i)) - 32)
 202      end if
 203    end do
 204
 205    end subroutine UPCASE
 206
 207 END MODULE module_chem_utilities