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 calculates vertrical advection rates using
6 : !! Piecewise Polynomial Method [Colela and Woodard, J. Comp. Phys.,
7 : !! 54, 174-201, 1984]
8 : !!
9 : !! @author Eric Jensen
10 : !! @version Dec-1996
11 145973900 : subroutine vertadv(carma, cstate, vtrans, cvert, itbnd, ibbnd, cvert_tbnd, cvert_bbnd, vertadvu, vertadvd, rc)
12 :
13 : ! types
14 : use carma_precision_mod
15 : use carma_enums_mod
16 : use carma_constants_mod
17 : use carma_types_mod
18 : use carmastate_mod
19 : use carma_mod
20 :
21 : implicit none
22 :
23 : type(carma_type), intent(in) :: carma !! the carma object
24 : type(carmastate_type), intent(inout) :: cstate !! the carma state object
25 : real(kind=f), intent(inout) :: vtrans(NZP1) !! vertical velocity
26 : real(kind=f), intent(in) :: cvert(NZ) !! quantity being transported
27 : integer, intent(in) :: itbnd !! top boundary condition
28 : integer, intent(in) :: ibbnd !! bottom boundary condition
29 : real(kind=f), intent(in) :: cvert_tbnd !! quantity at top boundary
30 : real(kind=f), intent(in) :: cvert_bbnd !! quantity at bottom boundary
31 : real(kind=f), intent(out) :: vertadvu(NZP1) !! upward vertical transport rate into level k from level k-1 [cm/s]
32 : real(kind=f), intent(out) :: vertadvd(NZP1) !! downward vertical transport rate into level k from level k-1 [cm/s]
33 : integer, intent(inout) :: rc !! return code, negative indicates failure
34 :
35 : ! Local declarations
36 : integer :: k
37 : integer :: nzm1
38 : integer :: nzm2
39 : integer :: itwo
40 291947800 : real(kind=f) :: dela(NZ)
41 291947800 : real(kind=f) :: delma(NZ)
42 291947800 : real(kind=f) :: aju(NZ)
43 291947800 : real(kind=f) :: ar(NZ)
44 291947800 : real(kind=f) :: al(NZ)
45 145973900 : real(kind=f) :: a6(NZ)
46 : real(kind=f) :: dpc, dpc1, dpcm1
47 : real(kind=f) :: ratt1, ratt2, ratt3, rat1, rat2, rat3, rat4, den1
48 : real(kind=f) :: com2, x, xpos
49 : real(kind=f) :: cvert0, cvertnzp1
50 :
51 :
52 : ! Initialize fluxes to zero
53 10510120800 : vertadvu(:) = 0._f
54 10510120800 : vertadvd(:) = 0._f
55 :
56 : ! If doing explicit sedimentation then do a simple sorting of positive and negative
57 : ! velocities into up and down components.
58 145973900 : if (do_explised) then
59 0 : where (vtrans < 0._f)
60 : vertadvd = -vtrans
61 : elsewhere
62 : vertadvu = vtrans
63 : end where
64 : else
65 :
66 :
67 145973900 : if( ibbnd .eq. I_FLUX_SPEC ) vtrans(1) = 0._f
68 145973900 : if( itbnd .eq. I_FLUX_SPEC ) vtrans(NZP1) = 0._f
69 :
70 : ! Set some constants
71 145973900 : nzm1 = max( 1, NZ-1 )
72 145973900 : nzm2 = max( 1, NZ-2 )
73 145973900 : itwo = min( 2, NZ )
74 :
75 : ! First, use cubic fits to estimate concentration values at layer
76 : ! boundaries
77 10072199100 : do k = 2,NZ-1
78 9926225200 : dpc = cvert(k) / dz(k)
79 9926225200 : dpc1 = cvert(k+1) / dz(k+1)
80 9926225200 : dpcm1 = cvert(k-1) / dz(k-1)
81 9926225200 : ratt1 = dz(k) / ( dz(k-1) + dz(k) + dz(k+1) )
82 9926225200 : ratt2 = ( 2._f*dz(k-1) + dz(k) ) / ( dz(k+1) + dz(k) )
83 9926225200 : ratt3 = ( 2._f*dz(k+1) + dz(k) ) / ( dz(k-1) + dz(k) )
84 9926225200 : dela(k) = ratt1 * ( ratt2*(dpc1-dpc) + ratt3*(dpc-dpcm1) )
85 :
86 10072199100 : if( (dpc1-dpc)*(dpc-dpcm1) .gt. 0._f .and. dela(k) .ne. 0._f ) then
87 8423047515 : delma(k) = min(abs(dela(k)), 2._f*abs(dpc-dpc1), 2._f*abs(dpc-dpcm1)) * abs(dela(k))/dela(k)
88 : else
89 1503177685 : delma(k) = 0._f
90 : endif
91 : enddo ! k = 2,NZ-2
92 :
93 9926225200 : do k = 2,NZ-2
94 9780251300 : dpc = cvert(k) / dz(k)
95 9780251300 : dpc1 = cvert(k+1) / dz(k+1)
96 9780251300 : dpcm1 = cvert(k-1) / dz(k-1)
97 9780251300 : rat1 = dz(k) / ( dz(k) + dz(k+1) )
98 9780251300 : rat2 = 2._f * dz(k+1) * dz(k) / ( dz(k) + dz(k+1) )
99 9780251300 : rat3 = ( dz(k-1) + dz(k) ) / ( 2._f*dz(k) + dz(k+1) )
100 9780251300 : rat4 = ( dz(k+2) + dz(k+1) ) / ( 2._f*dz(k+1) + dz(k) )
101 9780251300 : den1 = dz(k-1) + dz(k) + dz(k+1) + dz(k+2)
102 :
103 : ! <aju(k)> is the estimate for concentration (dn/dz) at layer
104 : ! boundary <k>+1/2.
105 : aju(k) = dpc + rat1*(dpc1-dpc) + 1._f/den1 * ( rat2*(rat3-rat4)*(dpc1-dpc) - &
106 9926225200 : dz(k)*rat3*delma(k+1) + dz(k+1)*rat4*delma(k) )
107 :
108 : enddo ! k = 2,NZ-2
109 :
110 : ! Now construct polynomial functions in each layer
111 9780251300 : do k = 3,NZ-2
112 9634277400 : al(k) = aju(k-1)
113 9780251300 : ar(k) = aju(k)
114 : enddo
115 :
116 : ! Use linear functions in first two and last two layers
117 145973900 : ar(itwo) = aju(itwo)
118 145973900 : al(itwo) = cvert(1)/dz(1) + (zl(itwo)-zc(1)) / &
119 145973900 : (zc(itwo)-zc(1)) * (cvert(itwo)/dz(itwo)-cvert(1)/dz(1))
120 145973900 : ar(1) = al(itwo)
121 0 : al(1) = cvert(1)/dz(1) - (zc(1)-zl(1)) / &
122 145973900 : (zc(itwo)-zc(1)) * (cvert(itwo)/dz(itwo)-cvert(1)/dz(1))
123 :
124 145973900 : al(nzm1) = aju(nzm2)
125 0 : ar(nzm1) = cvert(nzm1)/dz(nzm1) + (zl(NZ)-zc(nzm1)) &
126 145973900 : / (zc(NZ)-zc(nzm1)) * (cvert(NZ)/dz(NZ)-cvert(nzm1)/dz(nzm1))
127 145973900 : al(NZ) = ar(nzm1)
128 145973900 : ar(NZ) = cvert(nzm1)/dz(nzm1) + (zl(NZ+1)-zc(nzm1)) &
129 291947800 : / (zc(NZ)-zc(nzm1)) * (cvert(NZ)/dz(NZ)-cvert(nzm1)/dz(nzm1))
130 :
131 : ! Ensure that boundary values are not negative
132 145973900 : al(1) = max( al(1), 0._f )
133 145973900 : ar(NZ) = max( ar(NZ), 0._f )
134 :
135 : ! Next, ensure that polynomial functions do not deviate beyond the
136 : ! range [<al(k)>,<ar(k)>]
137 10364146900 : do k = 1,NZ
138 10218173000 : dpc = cvert(k) / dz(k)
139 10218173000 : if( (ar(k)-dpc)*(dpc-al(k)) .le. 0._f ) then
140 1503177685 : al(k) = dpc
141 1503177685 : ar(k) = dpc
142 : endif
143 :
144 10218173000 : if( (ar(k)-al(k))*( dpc - 0.5_f*(al(k)+ar(k)) ) .gt. 1._f/6._f*(ar(k)-al(k))**2 ) &
145 778173010 : al(k) = 3._f*dpc - 2._f*ar(k)
146 :
147 10218173000 : if( (ar(k)-al(k))*( dpc - 0.5_f*(al(k)+ar(k)) ) .lt. -1._f/6._f*(ar(k)-al(k))**2 ) &
148 3234918887 : ar(k) = 3._f*dpc - 2._f*al(k)
149 : enddo
150 :
151 : ! Calculate fluxes across each layer boundary
152 10364146900 : do k = 1,NZ
153 10218173000 : dpc = cvert(k) / dz(k)
154 10218173000 : dela(k) = ar(k) - al(k)
155 10364146900 : a6(k) = 6._f * ( dpc - 0.5_f*(ar(k)+al(k)) )
156 : enddo
157 :
158 10218173000 : do k = 1,NZ-1
159 10072199100 : com2 = ( dz(k) + dz(k+1) ) / 2._f
160 10072199100 : x = vtrans(k+1)*dtime/dz(k)
161 10072199100 : xpos = abs(x)
162 :
163 : ! Upward transport rate
164 10218173000 : if( vtrans(k+1) .gt. 0._f )then
165 :
166 10072199100 : if( x .lt. 1._f .and. cvert(k) .ne. 0._f )then
167 : vertadvu(k+1) = ( vtrans(k+1) * com2 ) * ( ( ar(k) - 0.5_f*dela(k)*x + &
168 8029066293 : (x/2._f - (x**2)/3._f)*a6(k) ) / cvert(k) )
169 :
170 : ! If Courant # > 1, use upwind advection
171 : else
172 2043132807 : vertadvu(k+1) = vtrans(k+1)
173 : endif
174 :
175 : ! Downward transport rate
176 0 : elseif( vtrans(k+1) .lt. 0._f )then
177 :
178 0 : if( x .gt. -1._f .and. cvert(k+1) .ne. 0._f )then
179 : vertadvd(k+1) = ( -vtrans(k+1) * com2 ) * &
180 : ( ( al(k+1) + 0.5_f*dela(k+1)*xpos + &
181 0 : ( xpos/2._f - (xpos**2)/3._f)*a6(k+1) ) / cvert(k+1) )
182 : else
183 0 : vertadvd(k+1) = -vtrans(k+1)
184 : endif
185 : endif
186 :
187 : enddo ! k = 1,NZ-1
188 :
189 : ! Lower boundary transport rates: If I_FIXED_CONC boundary
190 : ! condtion is selected, then use concentration assumed just beyond
191 : ! the lowest layer edge to calculate the transport rate across
192 : ! the bottom boundary of the model.
193 145973900 : if( ibbnd .eq. I_FIXED_CONC ) then
194 :
195 145973900 : com2 = ( dz(1) + dz(itwo) ) / 2._f
196 145973900 : x = vtrans(1)*dtime/dz(1)
197 145973900 : xpos = abs(x)
198 145973900 : cvert0 = cvert_bbnd
199 145973900 : if( vtrans(1) .gt. 0._f )then
200 :
201 145973900 : if( x .lt. 1._f .and. cvert0 .ne. 0._f )then
202 : vertadvu(1) = vtrans(1)/cvert0*com2 &
203 : * ( ar(1) - 0.5_f*dela(1)*x + &
204 0 : (x/2._f - (x**2)/3._f)*a6(1) )
205 : else
206 145973900 : vertadvu(1) = vtrans(1)
207 : endif
208 :
209 0 : elseif( vtrans(1) .lt. 0._f )then
210 :
211 0 : if( x .gt. -1._f .and. cvert(1) .ne. 0._f )then
212 : vertadvd(1) = -vtrans(1)/ &
213 : cvert(1)*com2 &
214 : * ( al(1) + 0.5_f*dela(1)*xpos + &
215 0 : (xpos/2._f - (xpos**2)/3._f)*a6(1) )
216 : else
217 0 : vertadvd(1) = -vtrans(1)
218 : endif
219 : endif
220 : endif
221 :
222 : ! Upper boundary transport rates
223 145973900 : if( itbnd .eq. I_FIXED_CONC ) then
224 :
225 145973900 : com2 = ( dz(NZ) + dz(nzm1) ) / 2._f
226 145973900 : x = vtrans(NZ+1)*dtime/dz(NZ)
227 145973900 : xpos = abs(x)
228 145973900 : cvertnzp1 = cvert_tbnd
229 :
230 145973900 : if( vtrans(NZ+1) .gt. 0._f )then
231 :
232 145973900 : if( x .lt. 1._f .and. cvert(NZ) .ne. 0._f )then
233 : vertadvu(NZ+1) = vtrans(NZ+1)/cvert(NZ)*com2 &
234 : * ( ar(NZ) - 0.5_f*dela(NZ)*x + &
235 142727855 : (x/2._f - (x**2)/3._f)*a6(NZ) )
236 : else
237 3246045 : vertadvu(NZ+1) = vtrans(NZ+1)
238 : endif
239 :
240 0 : elseif( vtrans(NZ+1) .lt. 0._f )then
241 :
242 0 : if( x .gt. -1._f .and. cvertnzp1 .ne. 0._f )then
243 : vertadvd(NZ+1) = -vtrans(NZ+1)/ &
244 : cvertnzp1*com2 &
245 0 : * ( al(NZ) + 0.5_f*dela(NZ)*xpos + &
246 0 : (xpos/2._f - (xpos**2)/3._f)*a6(NZ) )
247 : else
248 0 : vertadvd(NZ+1) = -vtrans(NZ+1)
249 : endif
250 : endif
251 : endif
252 : endif
253 :
254 : ! Return to caller with vertical transport rates.
255 145973900 : return
256 145973900 : end
|
