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

1     ! Colarco, May 24, 2007
2     ! sed command to replace F77 comments with F90: sed 's/^c/\!/g'
3     ! F90-version of original CARMA initgas.f routine (see comments below from
4     ! original routine header).
5     
6           subroutine initgas ( 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 :: igas, j, iz, iy, ix, &
17                      iztop, izbot, kb, ke, idk
18           real(kind=f) :: rh_init, rhi_init, rvap, xyzmet
19     
20     #include "carma_globaer.h"
21     
22           rc = 0
23     
24     #ifdef DEBUG
25            write(*,*) '+ initgas'
26     #endif
27     
28     !       subroutine initgas
29     !
30     !
31     !  @(#) initgas.f  Ackerman  Dec-1995
32     !  This routine initializes the atmospheric profiles of all gases.
33     !
34     !    gc       Gas concentration at layer mid-point [g/cm^3]
35     !
36     !  Presently the only vertical coordinate is altitude, zl [cm].
37     !
38     !
39     !  Argument list input:
40     !    None.
41     !
42     !  Argument list output:
43     !    None.
44     !
45     !
46     !  Include global constants and variables
47     !
48     !      include 'globaer.h'
49     !
50     !
51     !   Define formats
52     !
53         1 format(/,'Gas concentrations for ',a,'(initgas)',//, &
54             a3, 1x, 4(a11,4x), /) 
55         2 format(i3,1x,1p,3(e11.3,4x),0p,f11.3)
56     
57     !
58     !-------------------------------------------------------------------------------
59     !
60     !  Announce entry to this routine
61     !
62     !      if( DEBUG ) write(LUNOPRT,'(/,a)') 'Enter initgas'
63     !
64     !-------------------------------------------------------------------------------
65     !
66     !
67     !   Calculate vapor pressures
68     !
69           do iz = 1,NZ
70            do iy = 1,NY
71             do ix = 1,NX
72              call vaporp ( ix, iy, iz, carma, rc )
73             enddo
74            enddo
75           enddo
76     !
77     !
78     !   Parameters for water vapor profile
79     !
80           igas = 1
81     !
82     !
83     !   <rh_init> is initial relative humidity [%]
84     !
85           rh_init = 40._f
86           rhi_init = 60._f
87           rh_init = 100._f
88           rhi_init = 100._f
89           rh_init = 80._f
90     !
91     !
92     !   Gas constant for water vapor
93     !
94           rvap = RGAS/gwtmol(igas)
95     !
96     !
97     !   Loop over all spatial dimensions and gases
98     !
99           do igas = 1,NGAS
100            do ix = 1,NX
101             do iy = 1,NY
102              do iz = 1,NZ
103     
104                if( igas .eq. 1 )then
105     !
106     !
107     !   Water vapor concentration from relative humidity and vapor pressure
108     !
109     
110                  if( zc(ix,iy,iz) .ge. 5.e3 .and. &
111                      zc(ix,iy,iz) .le. 6.e3 ) then
112                   gc(ix,iy,iz,igas) = rh_init/100._f*pvapl(ix,iy,iz,igas) &
113                                   / ( rvap*t(ix,iy,iz) ) 
114                  else
115                   gc(ix,iy,iz,igas) = 0.25_f*pvapi(ix,iy,iz,igas) &
116                                   / ( rvap*t(ix,iy,iz) )
117                  endif
118     
119     !             if( t(ix,iy,iz) .ge. T0 ) then
120     !               gc(ix,iy,iz,igas) = rh_init/100._f*pvapl(ix,iy,iz,igas) &
121     !                                   / ( rvap*t(ix,iy,iz) )
122     !             else
123     !               gc(ix,iy,iz,igas) = rhi_init/100._f*pvapi(ix,iy,iz,igas) &
124     !                                   / ( rvap*t(ix,iy,iz) )
125     !             endif
126     
127                else
128                  write(LUNOPRT,'(/,a)') 'invalid <igas> in initgas.f'
129                  stop 1
130                endif
131     
132              enddo
133             enddo
134            enddo
135           enddo
136     !
137     !  PRC August 3, 2007
138     !  The convention for top and bottom needs to be agreed upon (concur with JAS)
139     !  I choose that "top" refers to level NZ and "bottom" refers to level 1 in 
140     !  all vertical coordinate 
141     !
142     !  Specify fluxes at top and bottom of model [kg/m^2/s]
143     !
144           ftopgas(:,:,:) = 0._f
145           fbotgas(:,:,:) = 0._f
146     !
147     !
148     !  Scale particle concentrations and boundary fluxes from 
149     !  cartesian coordinates to coordinate system specified by <igrid>
150     !
151     !
152     !  Pick indices for top and bottom layers
153     !
154           iztop = NZ
155           izbot = 1
156     
157           do igas = 1, NGAS
158     
159             gc(:,:,:,igas ) = gc(:,:,:,igas) * xmet*ymet*zmet
160             ftopgas(:,:,igas) = ftopgas(:,:,igas) * xmet(:,:,iztop)*ymet(:,:,iztop)
161             fbotgas(:,:,igas) = fbotgas(:,:,igas) * xmet(:,:,izbot)*ymet(:,:,izbot)
162     
163           enddo
164     !
165     !
166     !  Initialize <supsati> and <supsatl>
167     !
168           do iz = 1,NZ
169            do iy = 1,NY
170             do ix = 1,NX
171              call supersat ( ix, iy, iz, carma, rc )
172             enddo
173            enddo
174           enddo
175     !
176     !
177     !  Print gas concentrations at 1 horizontal grid point
178     !
179           ix = 1
180           iy = 1
181     !
182     !
183     !  Set vertical loop index to increment downwards
184     !
185           if( igridv .eq. I_CART )then
186             kb  = NZ
187             ke  = 1
188             idk = -1
189           else 
190             kb  = 1
191             ke  = NZ
192             idk = 1
193           endif
194     
195           do igas = 1,NGAS
196     
197             write(LUNOPRT,1) gasname(igas), &
198                              'iz','zc','gc [kg/m^3]','supsat','T [K]'
199     
200             do iz = kb,ke,idk
201               xyzmet = xmet(ix,iy,iz)*ymet(ix,iy,iz)*zmet(ix,iy,iz)
202               write(LUNOPRT,2) iz, zc(ix,iy,iz), gc(ix,iy,iz,igas)/xyzmet, &
203                  supsatl(ix,iy,iz,igas), t(ix,iy,iz)
204             enddo
205     
206           enddo
207     !
208     !
209     !  Specify the values of <gc> assumed just above(below) the top(bottom)
210     !  of the model domain.  
211     !
212           do igas = 1,NGAS
213            do iy = 1, NY
214             do ix = 1, NX
215               gc_topbnd(ix,iy,igas) = gc(ix,iy,NZ,igas)
216               gc_botbnd(ix,iy,igas) = gc(ix,iy,1,igas)
217             enddo
218            enddo
219           enddo
220     !
221     !  Return to caller with gas concentrations initialized.
222     !
223           return
224           end
225