Line data    Source code

       1             : ! Include shortname defintions, so that the F77 code does not have to be modified to
       2             : ! reference the CARMA structure.
       3             : #include "carma_globaer.h"
       4             : 
       5             : !! This routine setups up parameters related to the atmospheric state. It assumes that the
       6             : !! pressure, temperature, and dimensional fields (xc, dx, yc, dy, zc, zl) have already been
       7             : !! specified and all state arrays allocated via CARMASTATE_Create().
       8             : !!
       9             : !! @author Chuck Bardeen
      10             : !! @ version Feb-1995
      11             : !! @see CARMASTATE_Create
      12     1050624 : subroutine setupatm(carma, cstate, rescale, rc)
      13             : 
      14             :   ! types
      15             :   use carma_precision_mod
      16             :   use carma_enums_mod
      17             :   use carma_constants_mod
      18             :   use carma_types_mod
      19             :   use carmastate_mod
      20             :   use carma_mod
      21             : 
      22             :   implicit none
      23             : 
      24             :   type(carma_type), intent(in)         :: carma       !! the carma object
      25             :   type(carmastate_type), intent(inout) :: cstate      !! the carma state object
      26             :   logical, intent(in)                  :: rescale     !! rescale the fall velocity for zmet change, this is instead of realculating
      27             :   integer, intent(inout)               :: rc          !! return code, negative indicates failure
      28             : 
      29             :   ! Local declarations
      30             :   !--
      31             :   ! For air viscosity calculations
      32             :   ! Air viscosity <rmu> is from Sutherland's equation (using Smithsonian
      33             :   !   Meteorological Tables, in which there is a misprint -- T is deg_K, not
      34             :   !   deg_C.
      35             :   real(kind=f), parameter :: rmu_0 = 1.8325e-4_f
      36             :   real(kind=f), parameter :: rmu_t0 = 296.16_f
      37             :   real(kind=f), parameter :: rmu_c = 120._f
      38             :   real(kind=f), parameter :: rmu_const = rmu_0 * (rmu_t0 + rmu_c)
      39             : 
      40             :   integer :: ielem, ibin, i, j, ix, iy, iz, ie, ig, ip, igrp, jgrp, igroup
      41             : 
      42             : 
      43             :   ! Calculate the dry air density at each level, using the ideal gas
      44             :   ! law. This will be used to calculate zmet.
      45    34670592 :   rhoa(:) = p(:) / (R_AIR * t(:))
      46             : 
      47             :   ! Calculate the dimensions and the dimensional metrics.
      48    34670592 :   dz(:) = abs(zl(2:NZP1) - zl(1:NZ))
      49             : 
      50             :   ! Put the fall velocity back into cgs units, so that we can determine
      51             :   ! new metrics and then scale it back. This is optional and is done instead
      52             :   ! of recalculating everything from scratch to improve performance.
      53     1050624 :   if (rescale .and. (igridv /= I_CART)) then
      54           0 :     do ibin = 1, NBIN
      55           0 :       do igroup = 1, NGROUP
      56           0 :         vf(:, ibin, igroup)   = vf(:, ibin, igroup)  * zmetl(:)
      57           0 :         dkz(:, ibin, igroup)  = dkz(:, ibin, igroup) * (zmetl(:)**2)
      58             :       end do
      59             :     end do
      60             :   end if
      61             : 
      62             : 
      63             :   ! Vertical Metrics
      64     1050624 :   select case(igridv)
      65             :     ! Cartesian
      66             :     case (I_CART)
      67           0 :       zmet = 1._f
      68             : 
      69             :     ! Sigma
      70             :     case (I_SIG)
      71           0 :       zmet(:) = abs(((pl(1:NZ) - pl(2:NZP1)) / (zl(1:NZ) - zl(2:NZP1))) / &
      72           0 :         (GRAV * rhoa(:)))
      73             : 
      74             :     ! Hybrid
      75             :     case (I_HYBRID)
      76     4202496 :       zmet(:) = abs(((pl(1:NZ) - pl(2:NZP1)) / (zl(1:NZ) - zl(2:NZP1))) / &
      77    38873088 :         (GRAV * rhoa(:)))
      78             : 
      79             :     case default
      80           0 :       if (do_print) write(LUNOPRT,*) "setupatm:: ERROR - The specified vertical grid type (", igridv, &
      81           0 :         ") is not supported."
      82     1050624 :       rc = -1
      83             :   end select
      84             : 
      85             :   ! Interpolate the z metric to the grid box edges.
      86     1050624 :   if (NZ == 1) then
      87           0 :     zmetl(:) = zmet(1)
      88             :   else
      89             : 
      90             :     ! Extrpolate the top and bottom.
      91     1050624 :     zmetl(1)    = zmet(1)  + (zmet(2) - zmet(1))     / (zc(2) - zc(1))     * (zl(1) - zc(1))
      92     1050624 :     zmetl(NZP1) = zmet(NZ) + (zmet(NZ) - zmet(NZ-1)) / (zc(NZ) - zc(NZ-1)) * (zl(NZP1) - zc(NZ))
      93             : 
      94             :     ! Interpolate the middles.
      95     1050624 :     if (NZ > 2) then
      96    33619968 :       do iz = 2, NZ
      97    33619968 :         zmetl(iz) = zmet(iz-1) + (zmet(iz) - zmet(iz-1)) / (zc(iz) - zc(iz-1)) * (zl(iz) - zc(iz-1))
      98             :       end do
      99             :     end if
     100             :   end if
     101             : 
     102             : 
     103             :   ! Determine the z metrics at the grid box edges and then use this to put the
     104             :   ! fall velocity back into /x/y/z units.
     105     1050624 :   if (rescale .and. (igridv /= I_CART)) then
     106           0 :     do ibin = 1, NBIN
     107           0 :       do igroup = 1, NGROUP
     108           0 :         vf(:, ibin, igroup)   = vf(:, ibin, igroup)  / zmetl(:)
     109           0 :         dkz(:, ibin, igroup)  = dkz(:, ibin, igroup) / (zmetl(:)**2)
     110             :       end do
     111             :     end do
     112             :   end if
     113             : 
     114             : 
     115             :   ! Scale the density into the units carma wants (i.e. /x/y/z)
     116    34670592 :   rhoa(:) = rhoa(:) * zmet(:)
     117             : 
     118             :   ! Use the pressure difference across the cell and the fact that the
     119             :   ! atmosphere is hydrostatic to caclulate an average density in the
     120             :   ! grid box.
     121    34670592 :   rhoa_wet(:) = abs((pl(2:NZP1) - pl(1:NZ))) / (GRAV)
     122    34670592 :   rhoa_wet(:) = rhoa_wet(:) / dz(:)
     123             : 
     124             :   ! Calculate the thermal properties of the atmosphere.
     125    34670592 :   rmu(:)     = rmu_const / ( t(:) + rmu_c ) * (t(:) / rmu_t0 )**1.5_f
     126    34670592 :   thcond(:)  = (5.69_f + .017_f*(t(:) - T0)) * 4.186e2_f
     127     1050624 : end subroutine