Line data Source code
1 : module majorsp_diffusion
2 :
3 : !--------------------------------------------------------------------------
4 : ! This module computes the diffusion of major species (O2 and O) mass mixing
5 : ! ratio. This routine computes both the molecular and eddy diffusivity. This
6 : ! is adapted from the major species diffusion calculation of TIME-GCM.
7 : !
8 : ! Calling sequence:
9 : ! initialization:
10 : ! init
11 : ! call mspd_init
12 : !
13 : ! interfacing:
14 : ! tphysac
15 : ! (after vertical_diffusion_tend)
16 : ! call mspd_intr
17 : ! call mspdiff
18 : !
19 : !---------------------------Code history--------------------------------
20 : ! Adapted from TIME-GCM (comp.F): H.-L. Liu, Nov 2003
21 : !--------------------------------------------------------------------------
22 :
23 :
24 : use shr_kind_mod, only: r8 => shr_kind_r8
25 : use ppgrid, only: pcols, pver, pverp
26 : use constituents, only: pcnst, cnst_name, cnst_get_ind, cnst_mw
27 : use cam_history, only: outfld
28 :
29 : implicit none
30 :
31 : private ! Make default type private to the module
32 : save
33 : !-----------------------
34 : ! Public interfaces
35 : !-----------------------
36 : public mspd_init ! Initialization
37 : public mspd_intr ! Full routine
38 : !-----------------------
39 : ! Private data
40 : !-----------------------
41 :
42 : real(r8) :: rmass_o2, rmass_o1, rmass_n2 ! molecular weight kg/kmol
43 : real(r8) :: rmassinv_o2, rmassinv_o1, rmassinv_n2 ! 1/rmass_o2...
44 : real(r8) :: phi(2,3) ! mutual diffusion constants of
45 : ! major constituents
46 : real(r8) :: delta(2,2) ! unit matrix
47 :
48 : real(r8), parameter :: t00=273._r8 ! reference temperature
49 : real(r8), parameter :: ptref=5.e-5_r8 ! thermosphere reference pressure (Pa)
50 : real(r8), parameter :: tau=1.86e3_r8 ! diffusive time constant (sec).
51 : real(r8), parameter :: protonmass=1.6726e-27_r8 ! Proton mass (kg)
52 : real(r8), parameter :: mmrMin=1.e-20_r8 ! lower limit of o2 and o mixing ratio
53 : real(r8), parameter :: N2mmrMin=1.e-6_r8 ! lower limit of o2, o, and h mixing ratios
54 :
55 : integer :: indx_O2 ! cnst index for o2
56 : integer :: indx_O ! cnst index for o
57 : integer :: indx_H ! cnst index for h
58 : integer, parameter :: io2=1, io1=2 ! local indices to o2 , o respectively
59 : logical :: fixed_ubc(2) ! flag for fixed upper boundary condition
60 :
61 : real(r8) :: o2mmr_ubc(pcols) ! MMR of O2 at top boundary (specified)
62 : real(r8) :: ommr_ubc(pcols) ! MMR of O at top boundary
63 :
64 : character(len=8), private :: mjdiffnam(2) ! names of v-diff tendencies
65 :
66 : contains
67 :
68 : !===============================================================================
69 768 : subroutine mspd_init()
70 :
71 : !-------------------------------------------------------------------------------
72 : ! Define constants and coeficient matrices, phi and delta, in the initialization.
73 : !-------------------------------------------------------------------------------
74 : use constituents, only: cnst_mw, cnst_fixed_ubc
75 : use cam_history, only: addfld, add_default
76 : use phys_control, only: phys_getopts
77 :
78 : !------------------------------Arguments--------------------------------
79 :
80 : !---------------------------Local storage-------------------------------
81 : logical :: history_waccmx
82 :
83 768 : call phys_getopts(history_waccmx_out=history_waccmx)
84 :
85 : !-----------------------------------------------------------
86 : ! Get required molecular weights
87 : !-----------------------------------------------------------
88 768 : call cnst_get_ind('O2', indx_O2, abort=.true.)
89 768 : call cnst_get_ind('O', indx_O, abort=.true.)
90 768 : call cnst_get_ind('H', indx_H, abort=.true.)
91 :
92 768 : rmass_o2 = cnst_mw(indx_O2)
93 768 : rmass_o1 = cnst_mw(indx_O)
94 768 : rmass_n2 = 28._r8
95 :
96 768 : rmassinv_o2 = 1._r8/rmass_o2
97 768 : rmassinv_o1 = 1._r8/rmass_o1
98 768 : rmassinv_n2 = 1._r8/rmass_n2
99 :
100 : !--------------------------------------------------------------------
101 : ! Get fixed upper boundary flags and set vertical range for diffusion
102 : !--------------------------------------------------------------------
103 768 : fixed_ubc(io2) = cnst_fixed_ubc(indx_O2)
104 768 : fixed_ubc(io1) = cnst_fixed_ubc(indx_O)
105 :
106 : !------------------------------------------------
107 : ! Set diffusion constants and setup matrix
108 : !------------------------------------------------
109 2304 : phi(:,1)=(/0._r8 ,0.673_r8/)
110 2304 : phi(:,2)=(/1.35_r8,0._r8 /)
111 2304 : phi(:,3)=(/1.11_r8,0.769_r8/)
112 2304 : delta(:,1)=(/1._r8,0._r8/)
113 2304 : delta(:,2)=(/0._r8,1._r8/)
114 :
115 : ! Set names of major diffusion tendencies and declare them as history variables
116 768 : mjdiffnam(1) = 'MD'//cnst_name(indx_O2)
117 1536 : call addfld (mjdiffnam(1),(/ 'lev' /), 'A','kg/kg/s','Major diffusion of '//cnst_name(indx_O2))
118 768 : mjdiffnam(2) = 'MD'//cnst_name(indx_O)
119 1536 : call addfld (mjdiffnam(2),(/ 'lev' /), 'A','kg/kg/s','Major diffusion of '//cnst_name(indx_O))
120 :
121 1536 : call addfld ('MBARV' , (/ 'lev' /),'I','g/mole','Variable Mean Mass')
122 :
123 768 : if (history_waccmx) then
124 768 : call add_default (mjdiffnam(1), 1, ' ')
125 768 : call add_default (mjdiffnam(2), 1, ' ')
126 768 : call add_default ('MBARV', 1, ' ')
127 : end if
128 :
129 768 : end subroutine mspd_init
130 :
131 : !===============================================================================
132 21888 : subroutine mspd_intr(ztodt ,state ,ptend)
133 :
134 : !-------------------------------------------------------------------------------
135 : ! interface routine. output tendency.
136 : !-------------------------------------------------------------------------------
137 768 : use physics_types, only: physics_state, physics_ptend
138 : use upper_bc, only: ubc_get_vals
139 : use air_composition, only: rairv, mbarv
140 :
141 : !------------------------------Arguments--------------------------------
142 : real(r8), intent(in) :: ztodt ! 2 delta-t
143 : type(physics_state), intent(in) :: state ! Physics state variables
144 : type(physics_ptend), intent(inout) :: ptend ! indivdual parameterization tendencies
145 : !---------------------------Local storage-------------------------------
146 : real(r8) :: rztodt ! 1/ztodt
147 : real(r8) :: tendo2o(pcols,pver,2) ! temporary array for o2 and o tendency
148 : real(r8) :: ubc_mmr(pcols,pcnst) ! upper bndy mixing ratios (kg/kg)
149 : real(r8) :: ubc_t(pcols) ! upper bndy temperature (K)
150 : integer :: lchnk ! chunk identifier
151 : integer :: ncol ! number of atmospheric columns
152 : integer :: i, k ! indexing integers
153 :
154 : !--------------------------------------------------------------------------------------------
155 : ! local constants
156 : !--------------------------------------------------------------------------------------------
157 21888 : rztodt = 1._r8/ztodt
158 21888 : lchnk = state%lchnk
159 21888 : ncol = state%ncol
160 :
161 : !----------------------------------------------------------------------------------------------
162 : ! Store the o2 and o tendencies calculated from vertical_diffusion (due to eddy diffusion only)
163 : !----------------------------------------------------------------------------------------------
164 37012608 : tendo2o(:ncol,:,io2) = ptend%q(:ncol,:,indx_O2)
165 37012608 : tendo2o(:ncol,:,io1) = ptend%q(:ncol,:,indx_O)
166 :
167 : !----------------------------------------------------------------------
168 : ! Operate on copies of the input states, convert to tendencies at end.
169 : !----------------------------------------------------------------------
170 37012608 : ptend%q(:ncol,:,indx_O2) = state%q(:ncol,:,indx_O2)
171 37012608 : ptend%q(:ncol,:,indx_O) = state%q(:ncol,:,indx_O)
172 :
173 21888 : if (fixed_ubc(io2) .or. fixed_ubc(io1)) then
174 : !-------------------------------------------
175 : ! set upper boundary values of O2 and O MMR.
176 : !-------------------------------------------
177 0 : call ubc_get_vals( lchnk, ncol, state%pint, state%zi, ubc_t, ubc_mmr )
178 0 : o2mmr_ubc(:ncol) = ubc_mmr(:ncol,indx_O2)
179 0 : ommr_ubc(:ncol) = ubc_mmr(:ncol,indx_O)
180 : endif
181 :
182 : ! Since this is a combined tendency, retain the old name for output
183 : ! and debugging purposes.
184 21888 : ptend%name = trim(ptend%name)//"+mspd"
185 21888 : ptend%lq(indx_O2) = .TRUE.
186 21888 : ptend%lq(indx_O) = .TRUE.
187 : !---------------------------------------------
188 : ! Call the major species diffusion subroutine.
189 : !---------------------------------------------
190 : call mspdiff (lchnk ,ncol , &
191 : state%t ,ptend%q ,state%pmid ,state%pint , &
192 21888 : state%pdel ,ztodt ,rairv(:,:,lchnk), mbarv(:,:,lchnk))
193 :
194 : !---------------------------------------------
195 : ! Update O2 and O tendencies and output
196 : !---------------------------------------------
197 2867328 : do k=1,pver
198 37012608 : do i=1,ncol
199 34145280 : ptend%q(i,k,indx_O2) = (ptend%q(i,k,indx_O2)-state%q(i,k,indx_O2))*rztodt &
200 68290560 : +tendo2o(i,k,io2)
201 34145280 : ptend%q(i,k,indx_O) = (ptend%q(i,k,indx_O)-state%q(i,k,indx_O))*rztodt &
202 71136000 : +tendo2o(i,k,io1)
203 : enddo
204 : enddo
205 :
206 21888 : call outfld(mjdiffnam(1),ptend%q(1,1,indx_O2),pcols,lchnk)
207 21888 : call outfld(mjdiffnam(2),ptend%q(1,1,indx_O),pcols,lchnk)
208 :
209 21888 : end subroutine mspd_intr
210 :
211 :
212 : !===============================================================================
213 21888 : subroutine mspdiff (lchnk ,ncol , &
214 : t ,q ,pmid ,pint , &
215 : pdel ,ztodt ,rairv ,mbarv)
216 : !-----------------------------------------------------------------------
217 : ! Driver routine to compute major species diffusion (O2 and O).
218 :
219 : ! Turbulent diffusivities and boundary layer nonlocal transport terms are
220 : ! obtained from the turbulence module.
221 : !---------------------------Arguments------------------------------------
222 21888 : use ref_pres, only: lev0 => nbot_molec
223 : use physconst, only: gravit
224 :
225 : integer, intent(in) :: lchnk ! chunk identifier
226 : integer, intent(in) :: ncol ! number of atmospheric columns
227 : real(r8), intent(in) :: t(pcols,pver) ! temperature input
228 : real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressures
229 : real(r8), intent(in) :: pint(pcols,pverp) ! interface pressures
230 : real(r8), intent(in) :: pdel(pcols,pver) ! thickness between interfaces
231 : real(r8), intent(in) :: ztodt ! 2 delta-t
232 : real(r8), intent(in) :: rairv(pcols,pver) ! composition dependent gas "constant"
233 : real(r8), intent(in) :: mbarv(pcols,pver) ! composition dependent mean mass
234 :
235 : real(r8), intent(inout) :: q(pcols,pver,pcnst) ! constituents
236 :
237 : !---------------------------Local storage-------------------------------
238 : real(r8) :: o2(pcols,pver), o1(pcols,pver) ! o2, o1 mixing ratio (kg/kg moist air)
239 : real(r8) :: h_atom(pcols,pver) ! H mixing ratio
240 : real(r8) :: rztodt ! 1/ztodt
241 : real(r8) :: dz(pcols,pver) ! log-pressure interval between interfaces
242 : real(r8) :: rdz(pcols,pver) ! 1./dz, defined on midpoints
243 : real(r8) :: dzmid(pcols,pverp) ! log-pressure interval between midpoints
244 : real(r8) :: rdzmid(pcols,pverp) ! 1./dzmid, defined on interfaces
245 : real(r8) :: ak(pcols,2,2,2) ! coefficient matrix "Alfa"
246 : real(r8) :: ep(pcols,2,2) ! coefficient matrix
247 : real(r8) :: difk(pcols,pverp) ! eddy diffusion normalized by scale height (1/sec)
248 : real(r8) :: expzm(pcols,pver) ! exp(-z)=pmid/ptref
249 : real(r8) :: expzi(pcols,pverp) ! exp(-z)=pint/ptref
250 : real(r8) :: wks1(pcols)
251 : real(r8) :: wks3(pcols),wks4(pcols) ! temporary working arrays
252 : real(r8) :: psclht(pcols,pverp) ! pressure scale height
253 : real(r8) :: p_ubc(pcols) ! extrapolated pressure at upper boundary level, pmid(1)^2=pmid(2)*p_ubc
254 : real(r8) :: rair_ubc(pcols) ! extrapolated rair at upper boundary level
255 : real(r8) :: mbar_ubc(pcols) ! extrapolated mbar at upper boundary level
256 : real(r8) :: flb(pcols,2) ! lower boundary condition for o2 and o, now
257 : ! calculated locally.
258 : real(r8) :: fub(pcols,2) ! upper boundary condition for o2 and o, now
259 : ! calculated locally.
260 : real(r8) :: fk(pcols,2) ! temporary working array for rhs
261 : real(r8) :: pk(pcols,2,2) ! temp array for coefficients on lower diagonal
262 : real(r8) :: rk(pcols,2,2) ! temp array for coefficients on upper diagonal
263 : real(r8) :: qk(pcols,2,2) ! temp array for coefficients on diagonal
264 : real(r8) :: apk(2,2,pcols,pver) ! coefficients on lower diagonal
265 : real(r8) :: ark(2,2,pcols,pver) ! coefficients on upper diagonal
266 : real(r8) :: aqk(2,2,pcols,pver) ! coefficients on diagonal
267 : real(r8) :: rfk(2,pcols,pver) ! rhs of the array equation
268 : real(r8) :: betawk(2,2,pcols,pver) ! working arrays for blktri solver.
269 : real(r8) :: gammawk(2,2,pcols,pver) ! working arrays for blktri solver.
270 : real(r8) :: ywk(2,pcols,pver) ! working arrays for blktri solver.
271 : real(r8) :: xwk(2,pcols,pver) ! working arrays for blktri solver.
272 : integer :: nlevs ! number of levels
273 : real(r8) :: t_ubc(pcols) ! Temperature at top boundary
274 : integer :: i, k, km, kp, m, ktmp, isp, kk, kr
275 :
276 : !---------------------------------------------------
277 : ! Set vertical grid and get time step for diffusion
278 : !---------------------------------------------------
279 21888 : nlevs = lev0
280 21888 : rztodt = 1._r8/ztodt
281 :
282 : !------------------------------------------------------
283 : ! Get species to diffuse and set upper/lower boundaries
284 : !------------------------------------------------------
285 37012608 : o2(:ncol,:) = q(:ncol,:,indx_O2)
286 37012608 : o1(:ncol,:) = q(:ncol,:,indx_O)
287 37012608 : h_atom(:ncol,:) = q(:ncol,:,indx_H)
288 :
289 284544 : flb(:ncol,1) = o2(:ncol,lev0+1) ! fixed lower boundary condition
290 284544 : flb(:ncol,2) = o1(:ncol,lev0+1)
291 21888 : if(fixed_ubc(io2).or.fixed_ubc(io1)) then
292 0 : fub(:ncol,1) = o2mmr_ubc(:ncol) ! fixed upper boundary condition
293 0 : fub(:ncol,2) = ommr_ubc(:ncol)
294 : endif
295 :
296 : !------------------------------------------------------------------
297 : ! Get log-pressure intervals between midpoints and interface points
298 : !------------------------------------------------------------------
299 37012608 : dz(:ncol,:) = pdel(:ncol,:)/pmid(:ncol,:)
300 37012608 : rdz(:ncol,:) = 1._r8/dz(:ncol,:)
301 2845440 : do k=2,pver
302 36728064 : do i=1,ncol
303 33882624 : dzmid(i,k) = (pmid(i,k)-pmid(i,k-1))/pint(i,k)
304 36706176 : rdzmid(i,k) = 1._r8/dzmid(i,k)
305 : enddo
306 : enddo
307 284544 : do i=1,ncol
308 262656 : p_ubc(i) = pmid(i,1)*pmid(i,1)/pmid(i,2)
309 262656 : dzmid(i,1) = (pmid(i,1)-p_ubc(i))/pint(i,1)
310 284544 : rdzmid(i,1) = 1._r8/dzmid(i,1)
311 : enddo
312 284544 : do i=1,ncol
313 262656 : dzmid(i,pverp) = dzmid(i,pver)
314 284544 : rdzmid(i,pverp) = 1._r8/dzmid(i,pverp)
315 : enddo
316 :
317 : !------------------------------------------------------------------
318 : ! Get log-pressure intervals between midpoints and interface points
319 : !------------------------------------------------------------------
320 37297152 : expzi(:ncol,:) = pint(:ncol,:)/ptref
321 37012608 : expzm(:ncol,:) = pmid(:ncol,:)/ptref
322 :
323 : !------------------------------------------------------------------
324 : ! Get pressure scale height
325 : !------------------------------------------------------------------
326 2845440 : do k=2,pver
327 36728064 : do i=1,ncol
328 36706176 : psclht(i,k) = .5_r8*(rairv(i,k)*t(i,k)+rairv(i,k-1)*t(i,k-1))/gravit
329 : enddo
330 : enddo
331 284544 : do i=1,ncol
332 262656 : rair_ubc(i) = 1.5_r8*rairv(i,1)-.5_r8*rairv(i,2)
333 262656 : t_ubc(i) = 1.5_r8*t(i,1)-.5_r8*t(i,2)
334 262656 : psclht(i,1) = .5_r8*(rairv(i,1)*t(i,1)+rair_ubc(i)*t_ubc(i))/gravit
335 284544 : psclht(i,pverp) = psclht(i,pver)
336 : enddo
337 :
338 : !------------------------------------------------------------------
339 : ! Initialize scale height normalized eddy diffusion
340 : !------------------------------------------------------------------
341 2889216 : do k=1,pverp
342 37297152 : do i=1,ncol
343 37275264 : difk(i,k) = 0._r8 ! eddy diffusion already calculated in vertical_diffusion
344 : enddo
345 : enddo
346 :
347 21888 : call outfld ('MBARV', mbarv(:,:), pcols, lchnk)
348 :
349 : !------------------------------------------------------------------
350 : ! Set up mean mass working array
351 : !------------------------------------------------------------------
352 : ! ep, ak at the interface level immediately below midpoint level nbot_molec
353 :
354 : ! WKS4 = .5*(DMBAR/DZ)/MBAR
355 284544 : do i=1,ncol
356 787968 : wks4(i) = (mbarv(i,lev0)-mbarv(i,lev0+1))/ &
357 1072512 : (dzmid(i,lev0+1)*(mbarv(i,lev0+1)+mbarv(i,lev0)))
358 : enddo
359 :
360 : !-----------------------------------
361 : ! Calculate coefficient matrices
362 : !-----------------------------------
363 284544 : km = 1
364 284544 : kp = 2
365 284544 : do i=1, ncol
366 787968 : ep(i,io2,kp) = 1._r8-(2._r8/(mbarv(i,lev0+1)+mbarv(i,lev0)))* &
367 1050624 : (rmass_o2+(mbarv(i,lev0)-mbarv(i,lev0+1))*rdzmid(i,lev0+1))
368 : ep(i,io1,kp) = 1._r8-(2._r8/(mbarv(i,lev0+1)+mbarv(i,lev0)))* &
369 284544 : (rmass_o1+(mbarv(i,lev0)-mbarv(i,lev0+1))*rdzmid(i,lev0+1))
370 : enddo
371 :
372 :
373 65664 : do m=1,2
374 590976 : do i=1,ncol
375 1050624 : ak(i,io2,m,kp) = &
376 : -delta(io2,m)*(phi(io1,3)+(phi(io1,io2)-phi(io1,3))* &
377 1050624 : .5_r8*(o2(i,lev0+1)+o2(i,lev0)))-(1._r8-delta(io2,m))* &
378 2626560 : (phi(io2,m)-phi(io2,3))*.5_r8*(o2(i,lev0+1)+o2(i,lev0))
379 : ak(i,io1,m,kp) = &
380 : -delta(io1,m)*(phi(io2,3)+(phi(io2,io1)-phi(io2,3))* &
381 : .5_r8*(o1(i,lev0+1)+o1(i,lev0)))-(1._r8-delta(io1,m))* &
382 569088 : (phi(io1,m)-phi(io1,3))*.5_r8*(o1(i,lev0+1)+o1(i,lev0))
383 : enddo
384 : enddo
385 : !
386 : ! WKS1=MBAR/M3*(T00/(T0+T))*0.25/(TAU*DET(ak)) ak at the interface level
387 : ! immediately below midpoint level nbot_molec.
388 284544 : do i=1,ncol
389 787968 : wks1(i) = 0.5_r8*(mbarv(i,lev0+1)+mbarv(i,lev0))*rmassinv_n2* &
390 : (2._r8*t00/(t(i,lev0+1)+t(i,lev0)))**0.25_r8/ &
391 1072512 : (tau*(ak(i,1,1,kp)*ak(i,2,2,kp)-ak(i,1,2,kp)*ak(i,2,1,kp)))
392 : enddo
393 : !
394 : ! Complete calculation of ak at the interface level immediately below midpoint
395 : ! level nbot_molec.
396 65664 : do m=1,2
397 590976 : do i=1,ncol
398 525312 : ak(i,io2,m,kp) = ak(i,io2,m,kp)*wks1(i)
399 569088 : ak(i,io1,m,kp) = ak(i,io1,m,kp)*wks1(i)
400 : enddo
401 : enddo
402 :
403 21888 : km = 1
404 21888 : kp = 2
405 1838592 : do k=lev0,2,-1
406 23617152 : ktmp = km
407 23617152 : km = kp
408 23617152 : kp = ktmp
409 23617152 : do i=1,ncol
410 87201792 : ep(i,io2,kp) = 1._r8-(2._r8/(mbarv(i,k)+mbarv(i,k-1)))*(rmass_o2+ &
411 109002240 : (mbarv(i,k-1)-mbarv(i,k))*rdzmid(i,k))
412 : ep(i,io1,kp) = 1._r8-(2._r8/(mbarv(i,k)+mbarv(i,k-1)))*(rmass_o1+ &
413 23617152 : (mbarv(i,k-1)-mbarv(i,k))*rdzmid(i,k))
414 : enddo
415 :
416 5450112 : do m=1,2
417 49051008 : do i=1,ncol
418 130802688 : ak(i,io2,m,kp) = &
419 : -delta(io2,m)*(phi(io1,3)+(phi(io1,io2)-phi(io1,3))* &
420 87201792 : .5_r8*(o2(i,k)+o2(i,k-1)))- &
421 43600896 : (1._r8-delta(io2,m))*(phi(io2,m)-phi(io2,3))* &
422 218004480 : .5_r8*(o2(i,k)+o2(i,k-1))
423 :
424 : ak(i,io1,m,kp) = &
425 : -delta(io1,m)*(phi(io2,3)+(phi(io2,io1)-phi(io2,3))* &
426 : .5_r8*(o1(i,k)+o1(i,k-1)))- &
427 : (1._r8-delta(io1,m))*(phi(io1,m)-phi(io1,3))* &
428 47234304 : .5_r8*(o1(i,k)+o1(i,k-1))
429 :
430 : enddo
431 : enddo
432 :
433 : !---------------------------------------------
434 : ! Calculate coefficients for diagonals and rhs
435 : !---------------------------------------------
436 : !
437 : ! WKS1=MBAR/M3*(T00/(T0+T))**0.25/(TAU*DET(ALFA))
438 23617152 : do i=1,ncol
439 65401344 : wks1(i) = 0.5_r8*(mbarv(i,k)+mbarv(i,k-1))*rmassinv_n2* &
440 : (2._r8*t00/(t(i,k)+t(i,k-1)))**0.25_r8/ &
441 87201792 : (tau*(ak(i,1,1,kp)*ak(i,2,2,kp)-ak(i,1,2,kp)*ak(i,2,1,kp)))
442 21800448 : wks3(i) = wks4(i)
443 : wks4(i) = (mbarv(i,k-1)-mbarv(i,k))/ &
444 23617152 : (dzmid(i,k)*(mbarv(i,k)+mbarv(i,k-1)))
445 : enddo
446 :
447 : !
448 : ! FINISH CALCULATING AK(K+1/2) AND GENERATE PK, QK, RK
449 5450112 : do m=1,2
450 12716928 : do isp=io2,io1
451 98102016 : do i=1,ncol
452 87201792 : ak(i,isp,m,kp) = ak(i,isp,m,kp)*wks1(i)
453 :
454 174403584 : pk(i,isp,m) = (ak(i,isp,m,km)*(rdzmid(i,k+1)+ep(i,m,km)/2._r8)- &
455 : expzi(i,k+1)*difk(i,k+1)*(rdzmid(i,k+1)- &
456 261605376 : wks3(i))*delta(isp,m))*rdz(i,k)
457 :
458 87201792 : rk(i,isp,m) = (ak(i,isp,m,kp)*(rdzmid(i,k)-ep(i,m,kp)/2._r8)- &
459 : expzi(i,k)*difk(i,k)*(rdzmid(i,k)+ &
460 87201792 : wks4(i))*delta(isp,m))*rdz(i,k)
461 :
462 : qk(i,isp,m) = -(ak(i,isp,m,km)*(rdzmid(i,k+1)-ep(i,m,km)/2._r8)+ &
463 : ak(i,isp,m,kp)*(rdzmid(i,k)+ep(i,m,kp)/2._r8))*rdz(i,k)+ &
464 : ((expzi(i,k)*difk(i,k)*(rdzmid(i,k)-wks4(i))+ &
465 : expzi(i,k+1)*difk(i,k+1)*(rdzmid(i,k+1)+wks3(i)))* &
466 94468608 : rdz(i,k)+expzm(i,k)*rztodt)*delta(isp,m)
467 :
468 : enddo
469 : enddo
470 : enddo
471 :
472 23617152 : do i=1,ncol
473 21800448 : fk(i,io2) = expzm(i,k)*o2(i,k)*rztodt
474 23617152 : fk(i,io1) = expzm(i,k)*o1(i,k)*rztodt
475 : enddo
476 :
477 : !----------------------------
478 : ! Lower boundary
479 : !----------------------------
480 1816704 : if (k==lev0) then
481 65664 : do m=1,2
482 590976 : do i=1,ncol
483 525312 : fk(i,io2) = fk(i,io2)-pk(i,io2,m)*flb(i,m)
484 525312 : fk(i,io1) = fk(i,io1)-pk(i,io1,m)*flb(i,m)
485 1619712 : pk(i,:,m) = 0._r8
486 : enddo
487 : enddo
488 : endif
489 :
490 1816704 : kr = lev0-k+1
491 :
492 23617152 : do i=1,ncol
493 67218048 : do m=1,2
494 152603136 : do kk=1,2
495 87201792 : apk(kk,m,i,kr) = pk(i,kk,m)
496 87201792 : aqk(kk,m,i,kr) = qk(i,kk,m)
497 130802688 : ark(kk,m,i,kr) = rk(i,kk,m)
498 : enddo
499 : enddo
500 : enddo
501 :
502 23639040 : do i=1,ncol
503 21800448 : rfk(io2,i,kr) = fk(i,io2)
504 23617152 : rfk(io1,i,kr) = fk(i,io1)
505 : enddo
506 :
507 : enddo
508 :
509 : !----------------------------
510 : ! Upper boundary
511 : !----------------------------
512 21888 : k=1
513 21888 : ktmp = km
514 21888 : km = kp
515 21888 : kp = ktmp
516 21888 : if(fixed_ubc(io2).or.fixed_ubc(io1)) then
517 0 : do i=1,ncol
518 0 : mbar_ubc(i) = 1._r8/(o2mmr_ubc(i)*rmassinv_o2+ommr_ubc(i)*rmassinv_o1+ &
519 0 : (1._r8-o2mmr_ubc(i)-ommr_ubc(i))*rmassinv_n2)
520 0 : ep(i,io2,kp) = 1._r8-(2._r8/(mbarv(i,k)+mbar_ubc(i)))*(rmass_o2+ &
521 0 : (mbar_ubc(i)-mbarv(i,k))*rdzmid(i,k))
522 : ep(i,io1,kp) = 1._r8-(2._r8/(mbarv(i,k)+mbar_ubc(i)))*(rmass_o1+ &
523 0 : (mbar_ubc(i)-mbarv(i,k))*rdzmid(i,k))
524 : enddo
525 :
526 0 : do m=1,2
527 0 : do i=1,ncol
528 0 : ak(i,io2,m,kp) = &
529 : -delta(io2,m)*(phi(io1,3)+(phi(io1,io2)-phi(io1,3))* &
530 : .5_r8*(o2(i,k)+o2mmr_ubc(i)))- &
531 0 : (1._r8-delta(io2,m))*(phi(io2,m)-phi(io2,3))* &
532 0 : .5_r8*(o2(i,k)+o2mmr_ubc(i))
533 :
534 : ak(i,io1,m,kp) = &
535 : -delta(io1,m)*(phi(io2,3)+(phi(io2,io1)-phi(io2,3))* &
536 : .5_r8*(o1(i,k)+ommr_ubc(i)))- &
537 : (1._r8-delta(io1,m))*(phi(io1,m)-phi(io1,3))* &
538 0 : .5_r8*(o1(i,k)+ommr_ubc(i))
539 :
540 : enddo
541 : enddo
542 :
543 : !
544 : ! WKS1=MBAR/M3*(T00/(T0+T))**0.25/(TAU*DET(ALFA))
545 0 : do i=1,ncol
546 0 : wks1(i) = 0.5_r8*(mbarv(i,k)+mbar_ubc(i))*rmassinv_n2* &
547 : (2._r8*t00/(t(i,k)+t_ubc(i)))**0.25_r8/ &
548 0 : (tau*(ak(i,1,1,kp)*ak(i,2,2,kp)-ak(i,1,2,kp)*ak(i,2,1,kp)))
549 0 : wks3(i) = wks4(i)
550 : wks4(i) = (mbar_ubc(i)-mbarv(i,k))/ &
551 0 : (dzmid(i,k)*(mbarv(i,k)+mbar_ubc(i)))
552 : enddo
553 :
554 : !
555 : ! FINISH CALCULATING AK(K+1/2) AND GENERATE PK, QK, RK
556 0 : do m=1,2
557 0 : do isp=io2,io1
558 0 : do i=1,ncol
559 0 : ak(i,isp,m,kp) = ak(i,isp,m,kp)*wks1(i)
560 :
561 0 : pk(i,isp,m) = (ak(i,isp,m,km)*(rdzmid(i,k+1)+ep(i,m,km)/2._r8)- &
562 : expzi(i,k+1)*difk(i,k+1)*(rdzmid(i,k+1)- &
563 0 : wks3(i))*delta(isp,m))*rdz(i,k)
564 :
565 : rk(i,isp,m) = (ak(i,isp,m,kp)*(rdzmid(i,k)-ep(i,m,kp)/2._r8)- &
566 : expzi(i,k)*difk(i,k)*(rdzmid(i,k)+ &
567 0 : wks4(i))*delta(isp,m))*rdz(i,k)
568 :
569 : qk(i,isp,m) = -(ak(i,isp,m,km)*(rdzmid(i,k+1)-ep(i,m,km)/2._r8)+ &
570 : ak(i,isp,m,kp)*(rdzmid(i,k)+ep(i,m,kp)/2._r8))*rdz(i,k)+ &
571 : ((expzi(i,k)*difk(i,k)*(rdzmid(i,k)-wks4(i))+ &
572 : expzi(i,k+1)*difk(i,k+1)*(rdzmid(i,k+1)+wks3(i)))* &
573 0 : rdz(i,k)+expzm(i,k)*rztodt)*delta(isp,m)
574 :
575 : enddo
576 : enddo
577 : enddo
578 :
579 0 : do i=1,ncol
580 0 : fk(i,io2) = expzm(i,k)*o2(i,k)*rztodt
581 0 : fk(i,io1) = expzm(i,k)*o1(i,k)*rztodt
582 : enddo
583 0 : do m=1,2
584 0 : do i=1,ncol
585 0 : fk(i,io2) = fk(i,io2)-rk(i,io2,m)*fub(i,m)
586 0 : fk(i,io1) = fk(i,io1)-rk(i,io1,m)*fub(i,m)
587 0 : rk(i,:,m) = 0._r8
588 : enddo
589 : enddo
590 :
591 : else
592 :
593 284544 : do i=1,ncol
594 284544 : wks3(i) = wks4(i)
595 : enddo
596 65664 : do m=1,2
597 153216 : do isp=io2,io1
598 1181952 : do i=1,ncol
599 4202496 : pk(i,isp,m) = (ak(i,isp,m,km)*(rdzmid(i,k+1)+ep(i,m,km)/2._r8)- &
600 : expzi(i,k+1)*difk(i,k+1)*(rdzmid(i,k+1)- &
601 4202496 : wks3(i))*delta(isp,m))*rdz(i,k)
602 :
603 : qk(i,isp,m) = -(ak(i,isp,m,km)*(rdzmid(i,k+1)-ep(i,m,km)/2._r8)) &
604 : *rdz(i,k)+ &
605 : (expzi(i,k+1)*difk(i,k+1)*(rdzmid(i,k+1)+wks3(i))* &
606 1138176 : rdz(i,k)+expzm(i,k)*rztodt)*delta(isp,m)
607 :
608 : enddo
609 : enddo
610 : enddo
611 :
612 284544 : do i=1,ncol
613 262656 : fk(i,io2) = expzm(i,k)*o2(i,k)*rztodt
614 284544 : fk(i,io1) = expzm(i,k)*o1(i,k)*rztodt
615 : enddo
616 :
617 : endif
618 :
619 21888 : kr = lev0-k+1
620 :
621 284544 : do i=1,ncol
622 809856 : do m=1,2
623 1838592 : do kk=1,2
624 1050624 : apk(kk,m,i,kr) = pk(i,kk,m)
625 1050624 : aqk(kk,m,i,kr) = qk(i,kk,m)
626 1575936 : ark(kk,m,i,kr) = rk(i,kk,m)
627 : enddo
628 : enddo
629 : enddo
630 :
631 284544 : do i=1,ncol
632 262656 : rfk(io2,i,kr) = fk(i,io2)
633 284544 : rfk(io1,i,kr) = fk(i,io1)
634 : enddo
635 :
636 : !------------------------------------
637 : ! Call solver to get diffused species
638 : !------------------------------------
639 : call blktri(apk,aqk,ark,rfk,pcols,1,ncol,pver,1,nlevs, &
640 21888 : betawk, gammawk, ywk, xwk)
641 :
642 1860480 : do k=lev0,1,-1
643 1838592 : kr = lev0-k+1
644 23923584 : do i=1,ncol
645 22063104 : o2(i,k) = xwk(1,i,kr)
646 23901696 : o1(i,k) = xwk(2,i,kr)
647 : enddo
648 : enddo
649 :
650 : !---------------------------------------------------------------
651 : ! Ensure non-negative O2 and O and check for N2 greater than one
652 : !---------------------------------------------------------------
653 284544 : do i=1,ncol
654 22347648 : do k=1,lev0
655 :
656 22063104 : if (o2(i,k) < mmrMin) o2(i,k) = mmrMin
657 22063104 : if (o1(i,k) < mmrMin) o1(i,k) = mmrMin
658 :
659 22325760 : if(1._r8-mmrMin-o2(i,k)-o1(i,k)-h_atom(i,k) < 0._r8) then
660 1728 : o2(i,k) = o2(i,k)*((1._r8-N2mmrMin-h_atom(i,k))/(o2(i,k)+o1(i,k)))
661 1728 : o1(i,k) = o1(i,k)*((1._r8-N2mmrMin-h_atom(i,k))/(o2(i,k)+o1(i,k)))
662 : endif
663 : enddo
664 : enddo
665 :
666 37012608 : q(:ncol,:,indx_O2) = o2(:ncol,:)
667 37012608 : q(:ncol,:,indx_O) = o1(:ncol,:)
668 :
669 21888 : end subroutine mspdiff
670 :
671 :
672 : !===============================================================================
673 21888 : SUBROUTINE BLKTRI(A,B,C,F,IF,I1,I2,KF,K1,K2,BETA,GAMMA,Y,X)
674 : implicit none
675 : ! ****
676 : ! **** This procedure solves (I2-I1+1) tridiagonal block matrix
677 : ! **** systems in which all blocks are 2 x 2 matrices.
678 : ! ****
679 : ! **** Each system may be written:
680 : ! ****
681 : ! **** A(K) * X(K-1) + B(K) * X(K) + Z(K) * X(K+1) = F(K)
682 : ! ****
683 : ! **** where:
684 : ! ****
685 : ! **** K = K1,K2,1
686 : ! ****
687 : ! **** A(K), B(K), C(K) are given (2 x 2) matrices.
688 : ! ****
689 : ! **** The F(k) are given two componente vectors.
690 : ! ****
691 : ! **** The system is to be solved for the two component
692 : ! **** vectors, X(K).
693 : ! ****
694 : ! **** A(K1) = C(K2) = 0.
695 : ! ****
696 : ! **** BETA(K), GAMMA(K), (K = K1,K2,1), are work space for
697 : ! **** (2 x 2) matrices.
698 : ! ****
699 : ! **** Y(K), (K = K1,K2,1), is work space for two component
700 : ! **** vectors.
701 : ! ****
702 : ! **** Algorithm: (See Isaacson and Keller p55)
703 : ! ****
704 : ! **** Forward sweep from K = K1 to K = K2:
705 : ! ****
706 : ! **** BETA(K1) = B(K1)**(-1)
707 : ! ****
708 : ! **** Y(K1) = BETA(K1)*F(K1)
709 : ! ****
710 : ! **** GAMMA(K) = BETA(K)*C(K), K = K1,(K2-1),1
711 : ! ****
712 : ! **** BETA(K) = (B(K) - A(K)*GAMMA(K-1))**(-1),
713 : ! **** K = K1+1,K2,1
714 : ! ****
715 : ! **** Y(K) = BETA(K)*(F(K) - A(K)*Y(K-1)),
716 : ! **** K = K1+1,K2,1
717 : ! ****
718 : ! **** Backward sweep, K = K2,K1,-1
719 : ! ****
720 : ! **** X(K2) = Y(K2)
721 : ! ****
722 : ! **** X(K) = Y(K) - GAMMA(K)*X(K+1), K = K2-1,K1,-1
723 : ! ****
724 : ! **** Dimension statements:
725 : ! ****
726 : ! **** Block matrices are dimensioned thus:
727 : ! ****
728 : ! **** MATRIX(2,2,IF,KF)
729 : ! ****
730 : ! **** Two component vectors are similarly treated:
731 : ! ****
732 : ! **** VECTOR(2,IF,KF)
733 : ! ****
734 : ! **** Our block matrix scheme spans the range, (K = K1,K2,1),
735 : ! **** where (1 .LE. K1 .LT. K2 .LE. KF)
736 : ! ****
737 : ! **** Similarly, we are solving (I1-I2+1) systems
738 : ! **** simultaneously as the index, I, spans the range,
739 : ! **** (I = I1,I2,1), where (1 .LE. I1 .LT. I2 .LE. IF)
740 : ! ****
741 : ! ****
742 : ! **** Dimension statements:
743 : ! SUBROUTINE BLKTRI(A,B,C,F,IF,I1,I2,KF,K1,K2,BETA,GAMMA,Y,X)
744 : ! ****
745 : ! Args:
746 : integer,intent(in) :: if,i1,i2,kf,k1,k2
747 : real(r8),intent(in) :: a(2,2,if,kf), b(2,2,if,kf), c(2,2,if,kf), &
748 : f(2,if,kf)
749 : real(r8),intent(out) :: beta(2,2,if,kf), gamma(2,2,if,kf), &
750 : y(2,if,kf), x(2,if,kf)
751 : !
752 : ! DIMENSION A(2,2,IF,KF), B(2,2,IF,KF), C(2,2,IF,KF), F(2,IF,KF),
753 : ! 1 BETA(2,2,IF,KF), GAMMA(2,2,IF,KF), Y(2,IF,KF), X(2,IF,KF)
754 : !
755 : ! Local:
756 : integer :: i,k
757 : ! ****
758 : ! **** Lower boundary at K = K1
759 : ! ****
760 284544 : DO I = I1,I2
761 : ! ****
762 : ! **** Y(1,I,K1) = determinant(B(K1))
763 : ! ****
764 262656 : Y(1,I,K1) = B(1,1,I,K1)*B(2,2,I,K1) - B(1,2,I,K1)*B(2,1,I,K1)
765 : ! ****
766 : ! **** BETA(K1) = B(K1)**(-1)
767 : ! ****
768 262656 : BETA(1,1,I,K1) = B(2,2,I,K1)/Y(1,I,K1)
769 262656 : BETA(1,2,I,K1) = -B(1,2,I,K1)/Y(1,I,K1)
770 262656 : BETA(2,1,I,K1) = -B(2,1,I,K1)/Y(1,I,K1)
771 262656 : BETA(2,2,I,K1) = B(1,1,I,K1)/Y(1,I,K1)
772 : ! ****
773 : ! **** Y(K1) = BETA(K1)*F(K1)
774 : ! ****
775 262656 : Y(1,I,K1) = BETA(1,1,I,K1)*F(1,I,K1) + BETA(1,2,I,K1)*F(2,I,K1)
776 284544 : Y(2,I,K1) = BETA(2,1,I,K1)*F(1,I,K1) + BETA(2,2,I,K1)*F(2,I,K1)
777 : ENDDO
778 : ! ****
779 : ! **** Now deal with levels (K1+1),K2,1
780 : ! ****
781 1838592 : DO K = K1+1,K2
782 23639040 : DO I = I1,I2
783 : ! ****
784 : ! **** GAMMA(K-1) = BETA(K-1)*C(K-1)
785 : ! ****
786 43600896 : GAMMA(1,1,I,K-1) = BETA(1,1,I,K-1)*C(1,1,I,K-1) + &
787 43600896 : BETA(1,2,I,K-1)*C(2,1,I,K-1)
788 : GAMMA(1,2,I,K-1) = BETA(1,1,I,K-1)*C(1,2,I,K-1) + &
789 21800448 : BETA(1,2,I,K-1)*C(2,2,I,K-1)
790 : GAMMA(2,1,I,K-1) = BETA(2,1,I,K-1)*C(1,1,I,K-1) + &
791 21800448 : BETA(2,2,I,K-1)*C(2,1,I,K-1)
792 : GAMMA(2,2,I,K-1) = BETA(2,1,I,K-1)*C(1,2,I,K-1) + &
793 21800448 : BETA(2,2,I,K-1)*C(2,2,I,K-1)
794 : ! ****
795 : ! **** GAMMA(K) = B(K) - A(K)*GAMMA(K-1)
796 : ! ****
797 21800448 : GAMMA(1,1,I,K) = B(1,1,I,K) - A(1,1,I,K)*GAMMA(1,1,I,K-1) - &
798 21800448 : A(1,2,I,K)*GAMMA(2,1,I,K-1)
799 : GAMMA(1,2,I,K) = B(1,2,I,K) - A(1,1,I,K)*GAMMA(1,2,I,K-1) - &
800 21800448 : A(1,2,I,K)*GAMMA(2,2,I,K-1)
801 : GAMMA(2,1,I,K) = B(2,1,I,K) - A(2,1,I,K)*GAMMA(1,1,I,K-1) - &
802 21800448 : A(2,2,I,K)*GAMMA(2,1,I,K-1)
803 : GAMMA(2,2,I,K) = B(2,2,I,K) - A(2,1,I,K)*GAMMA(1,2,I,K-1) - &
804 21800448 : A(2,2,I,K)*GAMMA(2,2,I,K-1)
805 : ! ****
806 : ! **** Y(1,I,K) = determinant(GAMMA(K))
807 : ! ****
808 : Y(1,I,K) = GAMMA(1,1,I,K)*GAMMA(2,2,I,K) - &
809 21800448 : GAMMA(1,2,I,K)*GAMMA(2,1,I,K)
810 : ! ****
811 : ! **** BETA(K) = GAMMA(K)**(-1)
812 : ! ****
813 21800448 : BETA(1,1,I,K) = GAMMA(2,2,I,K)/Y(1,I,K)
814 21800448 : BETA(1,2,I,K) = -GAMMA(1,2,I,K)/Y(1,I,K)
815 21800448 : BETA(2,1,I,K) = -GAMMA(2,1,I,K)/Y(1,I,K)
816 21800448 : BETA(2,2,I,K) = GAMMA(1,1,I,K)/Y(1,I,K)
817 : ! ****
818 : ! **** X(K) = F(K) - A(K)*Y(K-1)
819 : ! ****
820 : X(1,I,K) = F(1,I,K) - A(1,1,I,K)*Y(1,I,K-1) - &
821 21800448 : A(1,2,I,K)*Y(2,I,K-1)
822 : X(2,I,K) = F(2,I,K) - A(2,1,I,K)*Y(1,I,K-1) - &
823 21800448 : A(2,2,I,K)*Y(2,I,K-1)
824 : ! ****
825 : ! **** Y(K) = BETA(K)*X(K)
826 : ! ****
827 21800448 : Y(1,I,K) = BETA(1,1,I,K)*X(1,I,K) + BETA(1,2,I,K)*X(2,I,K)
828 23617152 : Y(2,I,K) = BETA(2,1,I,K)*X(1,I,K) + BETA(2,2,I,K)*X(2,I,K)
829 : ENDDO
830 : ENDDO
831 : ! ****
832 : ! **** Backward sweep to determine final solution, X(K) for
833 : ! **** K = K2,K1,-1
834 : ! ****
835 : ! **** X(K2) = Y(K2)
836 : ! ****
837 284544 : DO I = I1,I2
838 262656 : X(1,I,K2) = Y(1,I,K2)
839 284544 : X(2,I,K2) = Y(2,I,K2)
840 : ENDDO
841 : ! ****
842 : ! **** X(K) = Y(K) - GAMMA(K)*X(K+1)
843 : ! ****
844 1838592 : DO K = K2-1,K1,-1
845 23639040 : DO I = I1,I2
846 65401344 : X(1,I,K) = Y(1,I,K) - GAMMA(1,1,I,K)*X(1,I,K+1) - &
847 65401344 : GAMMA(1,2,I,K)*X(2,I,K+1)
848 : X(2,I,K) = Y(2,I,K) - GAMMA(2,1,I,K)*X(1,I,K+1) - &
849 23617152 : GAMMA(2,2,I,K)*X(2,I,K+1)
850 :
851 : ENDDO
852 : ENDDO
853 21888 : RETURN
854 21888 : END SUBROUTINE BLKTRI
855 :
856 : end module majorsp_diffusion
|
