File: C:\NOAA\NEMS_11731\src\chem\gocart\src\Components\GOCART_GridComp\CARMA_GridComp\hydrostat.F90

1     ! Colarco, May 28, 2007
2     ! sed command to replace F77 comments with F90: sed 's/^c/\!/g'
3     ! F90-version of original CARMA hydrostat.f routine (see comments below from
4     ! original routine header).
5     
6            subroutine hydrostat ( carma, rc )
7     
8     !     types
9           use carma_types_mod
10     
11           implicit none
12     
13           integer, intent(out) :: rc
14     
15     !     Local declarations
16           integer :: iz, ios
17           real(kind=f), allocatable, dimension(:,:)  :: dp, ptop, pstar, &
18                                                         xymet, rhoa_mks
19     
20     #include "carma_globaer.h"
21     
22           rc = 0
23           allocate(dp(carma%NX,carma%NY), ptop(carma%NX,carma%NY), &
24                    pstar(carma%NX,carma%NY), xymet(carma%NX,carma%NY), &
25                    rhoa_mks(carma%NX,carma%NY), stat=ios)
26           if(ios /= 0) then
27            rc = 1
28            return
29           endif
30     
31     #ifdef DEBUG
32           write(*,*) '+ hydrostat'
33     #endif
34     
35     !
36     !
37     !  @(#) hydrostat.f  Ackerman  Jul-1997
38     !
39     !    This routine updates pressure by hydrostatic integration.
40     !    In sigma coordinates, it also updates vertical metric scale
41     !    factor <zmet>, temperature <t>, and scaled air density <rhoa>.
42     !
43     !  Argument list input:
44     !
45     !    None.
46     !
47     !  Argument list output:
48     !    None.
49     !
50     !
51     !  Include global constants and variables
52     !
53     !      include 'globaer.h'
54     !
55     !-------------------------------------------------------------------------------
56     !
57     !  Announce entry to this routine
58     !
59     !      if( DEBUG ) write(LUNOPRT,'(/,a)') 'Enter hydrostat'
60     !
61     !
62     !-------------------------------------------------------------------------------
63     !
64     !
65     !  cartesian coordinates
66     !
67             if( igridv .eq. I_CART )then
68     !
69     !
70     !  <ptop> is pressure at top of layer
71     !
72               iz = NZ
73               ptop = p_top
74               dp = dz(:,:,iz)*GRAV*rhoa(:,:,iz)
75               p(:,:,iz) = ptop * sqrt( 1. + dp/ptop )
76        
77               do iz = NZ-1,1,-1
78                 ptop = ptop + dp
79                 dp = dz(:,:,iz)*GRAV*rhoa(:,:,iz)
80                 p(:,:,iz) = ptop * sqrt( 1. + dp/ptop )
81               enddo
82               p_surf = ptop + dp
83     !
84     !
85     !  Sigma coordinates: first integrate for total column mass, then
86     !  update pressures and scaled air density, then temperatures, and then
87     !  air density (in cgs units) and vertical metric scale factor.
88     !
89             else if( igridv .eq. I_SIG )then
90     
91               pstar = 0._f
92               do iz = 1,NZ
93                 xymet = xmet(:,:,iz) * ymet(:,:,iz)
94                 pstar = pstar + dz(:,:,iz) * GRAV * rhoa(:,:,iz) / xymet
95               enddo
96               p_surf = p_top + pstar
97     !PRC
98           pstar = p_surf
99     
100               do iz = 1,NZ
101                 xymet = xmet(:,:,iz) * ymet(:,:,iz)
102     !            p(:,:,iz) = p_top(:,:) + zc(:,:,iz) * pstar
103     !PRC - change definition of p, but this balances exactly with input profile
104                 p(:,:,iz) = zc(:,:,iz) * p_surf
105     !PRC - it's here the troubles begin.  WHat is PTC anyway?
106                 rhoa(:,:,iz) = pstar / GRAV * xymet
107                 t(:,:,iz) = ptc(:,:,iz) / rhoa(:,:,iz) * &
108                             ( p(:,:,iz) / PREF )**(R_AIR/CP)
109                 rhoa_mks = p(:,:,iz) / ( R_AIR * t(:,:,iz) )
110                 zmet(:,:,iz) = pstar / ( GRAV * rhoa_mks )
111               enddo
112     
113             endif
114     !
115     !
116     !  Return to caller with pressures hydrostatically balanced.
117     !
118           deallocate(dp, ptop, pstar, xymet, rhoa_mks, stat=ios)
119           if(ios /= 0) rc = 1
120     
121           return
122           end
123