Line data Source code
1 : !===============================================================================
2 : ! Age of air test tracers
3 : ! provides dissipation rate and surface fluxes for diagnostic constituents
4 : !===============================================================================
5 :
6 : module aoa_tracers
7 :
8 : use shr_kind_mod, only: r8 => shr_kind_r8
9 : use spmd_utils, only: masterproc
10 : use ppgrid, only: pcols, pver
11 : use constituents, only: pcnst, cnst_add, cnst_name, cnst_longname
12 : use cam_logfile, only: iulog
13 : use ref_pres, only: pref_mid_norm
14 :
15 : implicit none
16 : private
17 : save
18 :
19 : ! Public interfaces
20 : public :: aoa_tracers_register ! register constituents
21 : public :: aoa_tracers_implements_cnst ! true if named constituent is implemented by this package
22 : public :: aoa_tracers_init_cnst ! initialize constituent field
23 : public :: aoa_tracers_init ! initialize history fields, datasets
24 : public :: aoa_tracers_timestep_init ! place to perform per timestep initialization
25 : public :: aoa_tracers_timestep_tend ! calculate tendencies
26 : public :: aoa_tracers_readnl ! read namelist options
27 :
28 : ! Private module data
29 :
30 : integer, parameter :: ncnst=4 ! number of constituents implemented by this module
31 :
32 : ! constituent names
33 : character(len=8), parameter :: c_names(ncnst) = (/'AOA1', 'AOA2', 'HORZ', 'VERT'/)
34 :
35 : ! constituent source/sink names
36 : character(len=8), parameter :: src_names(ncnst) = (/'AOA1SRC', 'AOA2SRC', 'HORZSRC', 'VERTSRC'/)
37 :
38 : integer :: ifirst ! global index of first constituent
39 : integer :: ixaoa1 ! global index for AOA1 tracer
40 : integer :: ixaoa2 ! global index for AOA2 tracer
41 : integer :: ixht ! global index for HORZ tracer
42 : integer :: ixvt ! global index for VERT tracer
43 :
44 : ! Data from namelist variables
45 : logical :: aoa_tracers_flag = .false. ! true => turn on test tracer code, namelist variable
46 : logical :: aoa_read_from_ic_file = .true. ! true => tracers initialized from IC file
47 :
48 : real(r8), parameter :: treldays = 15._r8
49 : real(r8), parameter :: vert_offset = 10._r8
50 :
51 : ! 15-days used for diagnostic of transport circulation and K-tensors
52 : ! relaxation (in the original papers PM-1987 and YSGD-2000) => Zonal Mean
53 : ! to evaluate eddy-fluxes for 2D-diagnostics, here relaxation to the GLOBAL MEAN IC
54 : ! it may help to keep gradients but will rule-out 2D-transport diagnostics
55 : ! in km to avoid negative values of vertical tracers
56 : ! VERT(k) = -7._r8*alog(hyam(k)+hybm(k)) + vert_offset
57 :
58 : ! PM-1987:
59 : ! Plumb, R. A., and J. D. Mahlman (1987), The zonally averaged transport
60 : ! characteristics of the GFDL general circulation/transport model,
61 : ! J. Atmos.Sci.,44, 298-327
62 :
63 : ! YSGD-2000:
64 : ! Yudin, Valery A., Sergey P. Smyshlyaev, Marvin A. Geller, Victor L. Dvortsov, 2000:
65 : ! Transport Diagnostics of GCMs and Implications for 2D Chemistry-Transport Model of
66 : ! Troposphere and Stratosphere. J. Atmos. Sci., 57, 673-699.
67 : ! doi: http://dx.doi.org/10.1175/1520-0469(2000)057<0673:TDOGAI>2.0.CO;2
68 :
69 : real(r8) :: qrel_vert(pver) ! = -7._r8*log(pref_mid_norm(k)) + vert_offset
70 :
71 : !===============================================================================
72 : contains
73 : !===============================================================================
74 :
75 : !================================================================================
76 1536 : subroutine aoa_tracers_readnl(nlfile)
77 :
78 : use namelist_utils, only: find_group_name
79 : use units, only: getunit, freeunit
80 : use mpishorthand
81 : use cam_abortutils, only: endrun
82 :
83 : implicit none
84 :
85 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
86 :
87 : ! Local variables
88 : integer :: unitn, ierr
89 : character(len=*), parameter :: subname = 'aoa_tracers_readnl'
90 :
91 :
92 : namelist /aoa_tracers_nl/ aoa_tracers_flag, aoa_read_from_ic_file
93 :
94 : !-----------------------------------------------------------------------------
95 :
96 1536 : if (masterproc) then
97 2 : unitn = getunit()
98 2 : open( unitn, file=trim(nlfile), status='old' )
99 2 : call find_group_name(unitn, 'aoa_tracers_nl', status=ierr)
100 2 : if (ierr == 0) then
101 0 : read(unitn, aoa_tracers_nl, iostat=ierr)
102 0 : if (ierr /= 0) then
103 0 : call endrun(subname // ':: ERROR reading namelist')
104 : end if
105 : end if
106 2 : close(unitn)
107 2 : call freeunit(unitn)
108 : end if
109 :
110 : #ifdef SPMD
111 1536 : call mpibcast(aoa_tracers_flag, 1, mpilog, 0, mpicom)
112 1536 : call mpibcast(aoa_read_from_ic_file, 1, mpilog, 0, mpicom)
113 : #endif
114 :
115 1536 : endsubroutine aoa_tracers_readnl
116 :
117 : !================================================================================
118 :
119 1536 : subroutine aoa_tracers_register
120 : !-----------------------------------------------------------------------
121 : !
122 : ! Purpose: register advected constituents
123 : !
124 : !-----------------------------------------------------------------------
125 : use physconst, only: cpair, mwdry
126 : !-----------------------------------------------------------------------
127 :
128 1536 : if (.not. aoa_tracers_flag) return
129 :
130 : call cnst_add(c_names(1), mwdry, cpair, 0._r8, ixaoa1, readiv=aoa_read_from_ic_file, &
131 0 : longname='Age-of_air tracer 1')
132 0 : ifirst = ixaoa1
133 : call cnst_add(c_names(2), mwdry, cpair, 0._r8, ixaoa2, readiv=aoa_read_from_ic_file, &
134 0 : longname='Age-of_air tracer 2')
135 : call cnst_add(c_names(3), mwdry, cpair, 1._r8, ixht, readiv=aoa_read_from_ic_file, &
136 0 : longname='horizontal tracer')
137 : call cnst_add(c_names(4), mwdry, cpair, 0._r8, ixvt, readiv=aoa_read_from_ic_file, &
138 0 : longname='vertical tracer')
139 :
140 : end subroutine aoa_tracers_register
141 :
142 : !===============================================================================
143 :
144 0 : function aoa_tracers_implements_cnst(name)
145 : !-----------------------------------------------------------------------
146 : !
147 : ! Purpose: return true if specified constituent is implemented by this package
148 : !
149 : !-----------------------------------------------------------------------
150 :
151 : character(len=*), intent(in) :: name ! constituent name
152 : logical :: aoa_tracers_implements_cnst ! return value
153 :
154 : !---------------------------Local workspace-----------------------------
155 : integer :: m
156 : !-----------------------------------------------------------------------
157 :
158 0 : aoa_tracers_implements_cnst = .false.
159 :
160 0 : if (.not. aoa_tracers_flag) return
161 :
162 0 : do m = 1, ncnst
163 0 : if (name == c_names(m)) then
164 0 : aoa_tracers_implements_cnst = .true.
165 : return
166 : end if
167 : end do
168 :
169 : end function aoa_tracers_implements_cnst
170 :
171 : !===============================================================================
172 :
173 0 : subroutine aoa_tracers_init_cnst(name, latvals, lonvals, mask, q)
174 :
175 : !-----------------------------------------------------------------------
176 : !
177 : ! Purpose: initialize test tracers mixing ratio fields
178 : ! This subroutine is called at the beginning of an initial run ONLY
179 : !
180 : !-----------------------------------------------------------------------
181 :
182 : character(len=*), intent(in) :: name
183 : real(r8), intent(in) :: latvals(:) ! lat in degrees (ncol)
184 : real(r8), intent(in) :: lonvals(:) ! lon in degrees (ncol)
185 : logical, intent(in) :: mask(:) ! Only initialize where .true.
186 : real(r8), intent(out) :: q(:,:) ! kg tracer/kg dry air (gcol, plev)
187 :
188 : integer :: m
189 : !-----------------------------------------------------------------------
190 :
191 0 : if (.not. aoa_tracers_flag) return
192 :
193 0 : do m = 1, ncnst
194 0 : if (name == c_names(m)) then
195 : ! pass global constituent index
196 0 : call init_cnst_3d(ifirst+m-1, latvals, lonvals, mask, q)
197 : endif
198 : end do
199 :
200 : end subroutine aoa_tracers_init_cnst
201 :
202 : !===============================================================================
203 :
204 1536 : subroutine aoa_tracers_init
205 :
206 : !-----------------------------------------------------------------------
207 : !
208 : ! Purpose: initialize age of air constituents
209 : ! (declare history variables)
210 : !-----------------------------------------------------------------------
211 :
212 : use cam_history, only: addfld, add_default
213 :
214 : integer :: m, mm, k
215 : !-----------------------------------------------------------------------
216 :
217 1536 : if (.not. aoa_tracers_flag) return
218 :
219 : ! Set names of tendencies and declare them as history variables
220 :
221 0 : do m = 1, ncnst
222 0 : mm = ifirst+m-1
223 0 : call addfld(cnst_name(mm), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(mm))
224 0 : call addfld(src_names(m), (/ 'lev' /), 'A', 'kg/kg/s', trim(cnst_name(mm))//' source/sink')
225 :
226 0 : call add_default (cnst_name(mm), 1, ' ')
227 0 : call add_default (src_names(m), 1, ' ')
228 : end do
229 :
230 0 : do k = 1,pver
231 0 : qrel_vert(k) = -7._r8*log(pref_mid_norm(k)) + vert_offset
232 : enddo
233 :
234 1536 : end subroutine aoa_tracers_init
235 :
236 : !===============================================================================
237 :
238 16128 : subroutine aoa_tracers_timestep_init( phys_state )
239 : !-----------------------------------------------------------------------
240 : ! Provides a place to reinitialize diagnostic constituents HORZ and VERT
241 : !-----------------------------------------------------------------------
242 :
243 1536 : use time_manager, only: get_curr_date
244 : use ppgrid, only: begchunk, endchunk
245 : use physics_types, only: physics_state
246 :
247 : type(physics_state), intent(inout), dimension(begchunk:endchunk), optional :: phys_state
248 :
249 :
250 : integer c, i, k, ncol
251 : integer yr, mon, day, tod
252 : !--------------------------------------------------------------------------
253 :
254 16128 : if (.not. aoa_tracers_flag) return
255 :
256 0 : call get_curr_date (yr,mon,day,tod)
257 :
258 0 : if ( day == 1 .and. tod == 0) then
259 0 : if (masterproc) then
260 0 : write(iulog,*) 'AGE_OF_AIR_CONSTITUENTS: RE-INITIALIZING HORZ/VERT CONSTITUENTS'
261 : endif
262 :
263 0 : do c = begchunk, endchunk
264 0 : ncol = phys_state(c)%ncol
265 0 : do k = 1, pver
266 0 : do i = 1, ncol
267 0 : phys_state(c)%q(i,k,ixht) = 2._r8 + sin(phys_state(c)%lat(i))
268 0 : phys_state(c)%q(i,k,ixvt) = qrel_vert(k)
269 : end do
270 : end do
271 : end do
272 :
273 : end if
274 :
275 32256 : end subroutine aoa_tracers_timestep_init
276 :
277 : !===============================================================================
278 :
279 2470608 : subroutine aoa_tracers_timestep_tend(state, ptend, cflx, landfrac, dt)
280 :
281 16128 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
282 : use cam_history, only: outfld
283 : use time_manager, only: get_nstep
284 :
285 : ! Arguments
286 : type(physics_state), intent(in) :: state ! state variables
287 : type(physics_ptend), intent(out) :: ptend ! package tendencies
288 : real(r8), intent(inout) :: cflx(pcols,pcnst) ! Surface constituent flux (kg/m^2/s)
289 : real(r8), intent(in) :: landfrac(pcols) ! Land fraction
290 : real(r8), intent(in) :: dt ! timestep
291 :
292 : !----------------- Local workspace-------------------------------
293 :
294 : integer :: i, k
295 : integer :: lchnk ! chunk identifier
296 : integer :: ncol ! no. of column in chunk
297 : integer :: nstep ! current timestep number
298 : real(r8) :: qrel ! value to be relaxed to
299 : real(r8) :: xhorz ! updated value of HORZ
300 : real(r8) :: xvert ! updated value of VERT
301 : logical :: lq(pcnst)
302 : real(r8) :: teul ! relaxation in 1/sec*dt/2 = k*dt/2
303 : real(r8) :: wimp ! 1./(1.+ k*dt/2)
304 : real(r8) :: wsrc ! teul*wimp
305 : !------------------------------------------------------------------
306 :
307 58824 : teul = .5_r8*dt/(86400._r8 * treldays) ! 1/2 for the semi-implicit scheme if dt=time step
308 58824 : wimp = 1._r8/(1._r8 +teul)
309 58824 : wsrc = teul*wimp
310 :
311 58824 : if (.not. aoa_tracers_flag) then
312 58824 : call physics_ptend_init(ptend,state%psetcols,'none') !Initialize an empty ptend for use with physics_update
313 : return
314 : end if
315 :
316 0 : lq(:) = .FALSE.
317 0 : lq(ixaoa1) = .TRUE.
318 0 : lq(ixaoa2) = .TRUE.
319 0 : lq(ixht) = .TRUE.
320 0 : lq(ixvt) = .TRUE.
321 0 : call physics_ptend_init(ptend,state%psetcols, 'aoa_tracers', lq=lq)
322 :
323 0 : nstep = get_nstep()
324 0 : lchnk = state%lchnk
325 0 : ncol = state%ncol
326 :
327 0 : do k = 1, pver
328 0 : do i = 1, ncol
329 :
330 : ! AOA1
331 0 : ptend%q(i,k,ixaoa1) = 0.0_r8
332 :
333 : ! AOA2
334 0 : ptend%q(i,k,ixaoa2) = 0.0_r8
335 :
336 : ! HORZ
337 0 : qrel = 2._r8 + sin(state%lat(i)) ! qrel should zonal mean
338 0 : xhorz = state%q(i,k,ixht)*wimp + wsrc*qrel ! Xnew = weight*3D-tracer + (1.-weight)*1D-tracer
339 0 : ptend%q(i,k,ixht) = (xhorz - state%q(i,k,ixht)) / dt ! Xnew = weight*3D-tracer + (1.-weight)*2D-tracer zonal mean
340 : ! Can be still used .... to diagnose fluxes OT-tracers
341 : ! VERT
342 0 : qrel = qrel_vert(k) ! qrel should zonal mean
343 0 : xvert = wimp*state%q(i,k,ixvt) + wsrc*qrel
344 0 : ptend%q(i,k,ixvt) = (xvert - state%q(i,k,ixvt)) / dt
345 :
346 : end do
347 : end do
348 :
349 : ! record tendencies on history files
350 0 : call outfld (src_names(1), ptend%q(:,:,ixaoa1), pcols, lchnk)
351 0 : call outfld (src_names(2), ptend%q(:,:,ixaoa2), pcols, lchnk)
352 0 : call outfld (src_names(3), ptend%q(:,:,ixht), pcols, lchnk)
353 0 : call outfld (src_names(4), ptend%q(:,:,ixvt), pcols, lchnk)
354 :
355 : ! Set tracer fluxes
356 0 : do i = 1, ncol
357 :
358 : ! AOA1
359 0 : cflx(i,ixaoa1) = 1.e-6_r8
360 :
361 : ! AOA2
362 0 : if (landfrac(i) .eq. 1._r8 .and. state%lat(i) .gt. 0.35_r8) then
363 0 : cflx(i,ixaoa2) = 1.e-6_r8 + 1e-6_r8*0.0434_r8*real(nstep,r8)*dt/(86400._r8*365._r8)
364 : else
365 0 : cflx(i,ixaoa2) = 0._r8
366 : endif
367 :
368 : ! HORZ
369 0 : cflx(i,ixht) = 0._r8
370 :
371 : ! VERT
372 0 : cflx(i,ixvt) = 0._r8
373 :
374 : end do
375 :
376 58824 : end subroutine aoa_tracers_timestep_tend
377 :
378 : !===========================================================================
379 :
380 0 : subroutine init_cnst_3d(m, latvals, lonvals, mask, q)
381 :
382 : integer, intent(in) :: m ! global constituent index
383 : real(r8), intent(in) :: latvals(:) ! lat in degrees (ncol)
384 : real(r8), intent(in) :: lonvals(:) ! lon in degrees (ncol)
385 : logical, intent(in) :: mask(:) ! Only initialize where .true.
386 : real(r8), intent(out) :: q(:,:) ! kg tracer/kg dry air (gcol,plev)
387 :
388 : integer :: j, k, gsize
389 : !-----------------------------------------------------------------------
390 :
391 0 : if (masterproc) then
392 0 : write(iulog,*) 'AGE-OF-AIR CONSTITUENTS: INITIALIZING ',cnst_name(m),m
393 : end if
394 :
395 0 : if (m == ixaoa1) then
396 :
397 0 : q(:,:) = 0.0_r8
398 :
399 0 : else if (m == ixaoa2) then
400 :
401 0 : q(:,:) = 0.0_r8
402 :
403 0 : else if (m == ixht) then
404 :
405 0 : gsize = size(q, 1)
406 0 : do j = 1, gsize
407 0 : q(j,:) = 2._r8 + sin(latvals(j))
408 : end do
409 :
410 0 : else if (m == ixvt) then
411 :
412 0 : do k = 1, pver
413 0 : do j = 1, size(q,1)
414 0 : q(j,k) = qrel_vert(k)
415 : end do
416 : end do
417 :
418 : end if
419 :
420 58824 : end subroutine init_cnst_3d
421 :
422 : !=====================================================================
423 :
424 :
425 : end module aoa_tracers
|
