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 : use time_manager, only: get_curr_date, get_start_date
15 : use time_manager, only: is_leapyear, timemgr_get_calendar_cf, get_calday
16 :
17 : implicit none
18 : private
19 :
20 : ! Public interfaces
21 : public :: aoa_tracers_register ! register constituents
22 : public :: aoa_tracers_implements_cnst ! true if named constituent is implemented by this package
23 : public :: aoa_tracers_init_cnst ! initialize constituent field
24 : public :: aoa_tracers_init ! initialize history fields, datasets
25 : public :: aoa_tracers_timestep_init ! place to perform per timestep initialization
26 : public :: aoa_tracers_timestep_tend ! calculate tendencies
27 : public :: aoa_tracers_readnl ! read namelist options
28 :
29 : ! Private module data
30 :
31 : integer, parameter :: ncnst=3 ! number of constituents implemented by this module
32 :
33 : ! constituent names
34 : character(len=4), parameter :: c_names(ncnst) = (/'AOA1', 'HORZ', 'VERT'/)
35 :
36 : ! constituent source/sink names
37 : character(len=7), parameter :: src_names(ncnst) = (/'AOA1SRC', 'HORZSRC', 'VERTSRC'/)
38 :
39 : integer :: ifirst = -1 ! global index of first constituent
40 : integer :: ixaoa = -1 ! global index for AOA1SRC tracer
41 : integer :: ixht = -1 ! global index for HORZ tracer
42 : integer :: ixvt = -1 ! 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) = -huge(1._r8) ! = -7._r8*log(pref_mid_norm(k)) + vert_offset
70 :
71 : integer :: yr0 = -huge(1)
72 : real(r8) :: calday0 = -huge(1._r8)
73 : real(r8) :: years = -huge(1._r8)
74 :
75 : !===============================================================================
76 : contains
77 : !===============================================================================
78 :
79 : !================================================================================
80 1536 : subroutine aoa_tracers_readnl(nlfile)
81 :
82 : use namelist_utils, only: find_group_name
83 : use cam_abortutils, only: endrun
84 : use spmd_utils, only: mpicom, masterprocid, mpi_logical, mpi_success
85 :
86 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
87 :
88 : ! Local variables
89 : integer :: unitn, ierr
90 : character(len=*), parameter :: subname = 'aoa_tracers_readnl'
91 :
92 : namelist /aoa_tracers_nl/ aoa_tracers_flag, aoa_read_from_ic_file
93 :
94 : !-----------------------------------------------------------------------------
95 :
96 1536 : if (masterproc) then
97 2 : open( newunit=unitn, file=trim(nlfile), status='old' )
98 2 : call find_group_name(unitn, 'aoa_tracers_nl', status=ierr)
99 2 : if (ierr == 0) then
100 2 : read(unitn, aoa_tracers_nl, iostat=ierr)
101 2 : if (ierr /= 0) then
102 0 : call endrun(subname // ':: ERROR reading namelist')
103 : end if
104 : end if
105 2 : close(unitn)
106 : end if
107 :
108 1536 : call mpi_bcast(aoa_tracers_flag, 1, mpi_logical, masterprocid, mpicom, ierr)
109 1536 : if (ierr/=mpi_success) then
110 0 : call endrun(subname//': MPI_BCAST ERROR: aoa_tracers_flag')
111 : end if
112 1536 : call mpi_bcast(aoa_read_from_ic_file, 1, mpi_logical, masterprocid, mpicom, ierr)
113 1536 : if (ierr/=mpi_success) then
114 0 : call endrun(subname//': MPI_BCAST ERROR: aoa_read_from_ic_file')
115 : end if
116 :
117 1536 : endsubroutine aoa_tracers_readnl
118 :
119 : !================================================================================
120 :
121 1536 : subroutine aoa_tracers_register
122 : !-----------------------------------------------------------------------
123 : !
124 : ! Purpose: register advected constituents
125 : !
126 : !-----------------------------------------------------------------------
127 : use physconst, only: cpair, mwdry
128 : !-----------------------------------------------------------------------
129 :
130 : integer :: k
131 :
132 1536 : if (.not. aoa_tracers_flag) return
133 :
134 : call cnst_add(c_names(1), mwdry, cpair, 0._r8, ixaoa, readiv=aoa_read_from_ic_file, &
135 1536 : longname='mixing ratio LB tracer', cam_outfld=.false.)
136 :
137 : call cnst_add(c_names(2), mwdry, cpair, 1._r8, ixht, readiv=aoa_read_from_ic_file, &
138 1536 : longname='horizontal tracer', cam_outfld=.false.)
139 : call cnst_add(c_names(3), mwdry, cpair, 0._r8, ixvt, readiv=aoa_read_from_ic_file, &
140 1536 : longname='vertical tracer', cam_outfld=.false.)
141 :
142 1536 : ifirst = ixaoa
143 :
144 87552 : do k = 1,pver
145 87552 : qrel_vert(k) = -7._r8*log(pref_mid_norm(k)) + vert_offset
146 : enddo
147 :
148 : end subroutine aoa_tracers_register
149 :
150 : !===============================================================================
151 :
152 142848 : function aoa_tracers_implements_cnst(name)
153 : !-----------------------------------------------------------------------
154 : !
155 : ! Purpose: return true if specified constituent is implemented by this package
156 : !
157 : !-----------------------------------------------------------------------
158 :
159 : character(len=*), intent(in) :: name ! constituent name
160 : logical :: aoa_tracers_implements_cnst ! return value
161 :
162 : !---------------------------Local workspace-----------------------------
163 : integer :: m
164 : !-----------------------------------------------------------------------
165 :
166 142848 : aoa_tracers_implements_cnst = .false.
167 :
168 142848 : if (.not. aoa_tracers_flag) return
169 :
170 571392 : do m = 1, ncnst
171 571392 : if (name == c_names(m)) then
172 142848 : aoa_tracers_implements_cnst = .true.
173 : return
174 : end if
175 : end do
176 :
177 : end function aoa_tracers_implements_cnst
178 :
179 : !===============================================================================
180 :
181 0 : subroutine aoa_tracers_init_cnst(name, latvals, lonvals, mask, q)
182 :
183 : !-----------------------------------------------------------------------
184 : !
185 : ! Purpose: initialize test tracers mixing ratio fields
186 : ! This subroutine is called at the beginning of an initial run ONLY
187 : !
188 : !-----------------------------------------------------------------------
189 :
190 : character(len=*), intent(in) :: name
191 : real(r8), intent(in) :: latvals(:) ! lat in degrees (ncol)
192 : real(r8), intent(in) :: lonvals(:) ! lon in degrees (ncol)
193 : logical, intent(in) :: mask(:) ! Only initialize where .true.
194 : real(r8), intent(out) :: q(:,:) ! kg tracer/kg dry air (gcol, plev)
195 :
196 : integer :: m
197 : !-----------------------------------------------------------------------
198 :
199 0 : if (.not. aoa_tracers_flag) return
200 :
201 0 : do m = 1, ncnst
202 0 : if (name == c_names(m)) then
203 : ! pass global constituent index
204 0 : call init_cnst_3d(ifirst+m-1, latvals, lonvals, mask, q)
205 : endif
206 : end do
207 :
208 : end subroutine aoa_tracers_init_cnst
209 :
210 : !===============================================================================
211 :
212 1536 : subroutine aoa_tracers_init
213 :
214 : !-----------------------------------------------------------------------
215 : !
216 : ! Purpose: initialize age of air constituents
217 : ! (declare history variables)
218 : !-----------------------------------------------------------------------
219 :
220 : use cam_history, only: addfld, add_default
221 :
222 : integer :: m, mm
223 : integer :: yr, mon, day, sec, ymd
224 :
225 : !-----------------------------------------------------------------------
226 :
227 1536 : if (.not. aoa_tracers_flag) return
228 :
229 : ! Set names of tendencies and declare them as history variables
230 :
231 6144 : do m = 1, ncnst
232 4608 : mm = ifirst+m-1
233 9216 : call addfld(cnst_name(mm), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(mm))
234 9216 : call addfld(src_names(m), (/ 'lev' /), 'A', 'kg/kg/s', trim(cnst_name(mm))//' source/sink')
235 :
236 4608 : call add_default (cnst_name(mm), 1, ' ')
237 6144 : call add_default (src_names(m), 1, ' ')
238 : end do
239 :
240 1536 : call get_start_date(yr, mon, day, sec)
241 :
242 1536 : ymd = yr*10000 + mon*100 + day
243 :
244 1536 : yr0 = yr
245 1536 : calday0 = get_calday(ymd, sec)
246 :
247 1536 : end subroutine aoa_tracers_init
248 :
249 : !===============================================================================
250 :
251 16128 : subroutine aoa_tracers_timestep_init( phys_state )
252 : !-----------------------------------------------------------------------
253 : ! Provides a place to reinitialize diagnostic constituents HORZ and VERT
254 : !-----------------------------------------------------------------------
255 :
256 1536 : use ppgrid, only: begchunk, endchunk
257 : use physics_types, only: physics_state
258 :
259 : type(physics_state), intent(inout), dimension(begchunk:endchunk), optional :: phys_state
260 :
261 : integer c, i, k, ncol
262 : integer yr, mon, day, tod, ymd
263 : real(r8) :: calday, dpy
264 : !--------------------------------------------------------------------------
265 :
266 16128 : if (.not. aoa_tracers_flag) return
267 :
268 16128 : call get_curr_date (yr,mon,day,tod)
269 :
270 16128 : if ( day == 1 .and. tod == 0) then
271 768 : if (masterproc) then
272 1 : write(iulog,*) 'AGE_OF_AIR_CONSTITUENTS: RE-INITIALIZING HORZ/VERT CONSTITUENTS'
273 : endif
274 :
275 4608 : do c = begchunk, endchunk
276 3840 : ncol = phys_state(c)%ncol
277 219648 : do k = 1, pver
278 3315456 : do i = 1, ncol
279 3096576 : phys_state(c)%q(i,k,ixht) = 2._r8 + sin(phys_state(c)%lat(i))
280 3311616 : phys_state(c)%q(i,k,ixvt) = qrel_vert(k)
281 : end do
282 : end do
283 : end do
284 :
285 : end if
286 :
287 16128 : ymd = yr*10000 + mon*100 + day
288 16128 : calday = get_calday(ymd, tod)
289 :
290 16128 : dpy = 365._r8
291 16128 : if (timemgr_get_calendar_cf() == 'gregorian' .and. is_leapyear(yr)) then
292 16128 : dpy = 366._r8
293 : end if
294 16128 : years = (yr-yr0) + (calday-calday0)/dpy
295 :
296 32256 : end subroutine aoa_tracers_timestep_init
297 :
298 : !===============================================================================
299 :
300 18021120 : subroutine aoa_tracers_timestep_tend(state, ptend, dt)
301 :
302 16128 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
303 : use cam_history, only: outfld
304 :
305 : ! Arguments
306 : type(physics_state), intent(in) :: state ! state variables
307 : type(physics_ptend), intent(out) :: ptend ! package tendencies
308 : real(r8), intent(in) :: dt ! timestep size (sec)
309 :
310 : !----------------- Local workspace-------------------------------
311 :
312 : integer :: i, k
313 : integer :: lchnk ! chunk identifier
314 : integer :: ncol ! no. of column in chunk
315 : real(r8) :: qrel ! value to be relaxed to
316 : real(r8) :: xhorz ! updated value of HORZ
317 : real(r8) :: xvert ! updated value of VERT
318 : logical :: lq(pcnst)
319 : real(r8) :: teul ! relaxation in 1/sec*dt/2 = k*dt/2
320 : real(r8) :: wimp ! 1./(1.+ k*dt/2)
321 : real(r8) :: wsrc ! teul*wimp
322 :
323 : real(r8) :: xmmr
324 : real(r8), parameter :: mmr0 = 1.0e-6_r8 ! initial lower boundary mmr
325 : real(r8), parameter :: per_yr = 0.02_r8 ! fractional increase per year
326 :
327 : real(r8) :: mmr_out(pcols,pver,ncnst)
328 :
329 : !------------------------------------------------------------------
330 :
331 72960 : teul = .5_r8*dt/(86400._r8 * treldays) ! 1/2 for the semi-implicit scheme if dt=time step
332 72960 : wimp = 1._r8/(1._r8 +teul)
333 72960 : wsrc = teul*wimp
334 :
335 72960 : if (.not. aoa_tracers_flag) then
336 0 : call physics_ptend_init(ptend,state%psetcols,'none') !Initialize an empty ptend for use with physics_update
337 : return
338 : end if
339 :
340 72960 : lq(:) = .FALSE.
341 72960 : lq(ixaoa) = .TRUE.
342 72960 : lq(ixht) = .TRUE.
343 72960 : lq(ixvt) = .TRUE.
344 :
345 72960 : call physics_ptend_init(ptend,state%psetcols, 'aoa_tracers', lq=lq)
346 :
347 72960 : lchnk = state%lchnk
348 72960 : ncol = state%ncol
349 :
350 : ! AOA1
351 72960 : xmmr = mmr0*(1._r8 + per_yr*years)
352 1123584 : ptend%q(1:ncol,pver,ixaoa) = (xmmr - state%q(1:ncol,pver,ixaoa)) / dt
353 :
354 4158720 : do k = 1, pver
355 62993664 : do i = 1, ncol
356 :
357 : ! HORZ
358 58834944 : qrel = 2._r8 + sin(state%lat(i)) ! qrel should zonal mean
359 58834944 : xhorz = state%q(i,k,ixht)*wimp + wsrc*qrel ! Xnew = weight*3D-tracer + (1.-weight)*1D-tracer
360 58834944 : ptend%q(i,k,ixht) = (xhorz - state%q(i,k,ixht)) / dt ! Xnew = weight*3D-tracer + (1.-weight)*2D-tracer zonal mean
361 : ! Can be still used .... to diagnose fluxes OT-tracers
362 : ! VERT
363 58834944 : qrel = qrel_vert(k) ! qrel should zonal mean
364 58834944 : xvert = wimp*state%q(i,k,ixvt) + wsrc*qrel
365 62920704 : ptend%q(i,k,ixvt) = (xvert - state%q(i,k,ixvt)) / dt
366 :
367 : end do
368 :
369 : end do
370 :
371 : ! record tendencies on history files
372 72960 : call outfld (src_names(1), ptend%q(:,:,ixaoa), pcols, lchnk)
373 72960 : call outfld (src_names(2), ptend%q(:,:,ixht), pcols, lchnk)
374 72960 : call outfld (src_names(3), ptend%q(:,:,ixvt), pcols, lchnk)
375 :
376 : ! output mixing ratios to history
377 62993664 : mmr_out(:ncol,:,1) = state%q(:ncol,:,ixaoa) + dt*ptend%q(1:ncol,:,ixaoa)
378 62993664 : mmr_out(:ncol,:,2) = state%q(:ncol,:,ixht) + dt*ptend%q(1:ncol,:,ixht)
379 62993664 : mmr_out(:ncol,:,3) = state%q(:ncol,:,ixvt) + dt*ptend%q(1:ncol,:,ixvt)
380 :
381 72960 : call outfld (c_names(1), mmr_out(:,:,1), pcols, lchnk)
382 72960 : call outfld (c_names(2), mmr_out(:,:,2), pcols, lchnk)
383 72960 : call outfld (c_names(3), mmr_out(:,:,3), pcols, lchnk)
384 :
385 72960 : end subroutine aoa_tracers_timestep_tend
386 :
387 : !===========================================================================
388 :
389 0 : subroutine init_cnst_3d(m, latvals, lonvals, mask, q)
390 :
391 : integer, intent(in) :: m ! global constituent index
392 : real(r8), intent(in) :: latvals(:) ! lat in degrees (ncol)
393 : real(r8), intent(in) :: lonvals(:) ! lon in degrees (ncol)
394 : logical, intent(in) :: mask(:) ! Only initialize where .true.
395 : real(r8), intent(out) :: q(:,:) ! kg tracer/kg dry air (gcol,plev)
396 :
397 : integer :: j, k, gsize
398 : !-----------------------------------------------------------------------
399 :
400 0 : if (masterproc) then
401 0 : write(iulog,*) 'AGE-OF-AIR CONSTITUENTS: INITIALIZING ',cnst_name(m),m
402 : end if
403 :
404 0 : if (m == ixaoa) then
405 :
406 : ! AOA1
407 0 : q(:,:) = 0.0_r8
408 :
409 0 : else if (m == ixht) then
410 :
411 : ! HORZ
412 0 : gsize = size(q, 1)
413 0 : do j = 1, gsize
414 0 : q(j,:) = 2._r8 + sin(latvals(j))
415 : end do
416 :
417 0 : else if (m == ixvt) then
418 :
419 : ! VERT
420 0 : do k = 1, pver
421 0 : do j = 1, size(q,1)
422 0 : q(j,k) = qrel_vert(k)
423 : end do
424 : end do
425 :
426 : end if
427 :
428 72960 : end subroutine init_cnst_3d
429 :
430 : !=====================================================================
431 :
432 : end module aoa_tracers
|
