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=6), parameter :: c_names(ncnst) = (/'AOAMF ', 'HORZ ', 'VERT '/)
35 :
36 : ! constituent source/sink names
37 : character(len=8), parameter :: src_names(ncnst) = (/'AOAMFSRC', 'HORZSRC ', 'VERTSRC '/)
38 :
39 : integer :: ifirst = -1 ! global index of first constituent
40 : integer :: ixaoa = -1 ! global index for AOAMFSRC 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 0 : read(unitn, aoa_tracers_nl, iostat=ierr)
101 0 : 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 0 : longname='mixing ratio LB tracer')
136 :
137 : call cnst_add(c_names(2), mwdry, cpair, 1._r8, ixht, readiv=aoa_read_from_ic_file, &
138 0 : longname='horizontal tracer')
139 : call cnst_add(c_names(3), mwdry, cpair, 0._r8, ixvt, readiv=aoa_read_from_ic_file, &
140 0 : longname='vertical tracer')
141 :
142 0 : ifirst = ixaoa
143 :
144 0 : do k = 1,pver
145 0 : 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 0 : 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 0 : aoa_tracers_implements_cnst = .false.
167 :
168 0 : if (.not. aoa_tracers_flag) return
169 :
170 0 : do m = 1, ncnst
171 0 : if (name == c_names(m)) then
172 0 : 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 0 : do m = 1, ncnst
232 0 : mm = ifirst+m-1
233 0 : call addfld(cnst_name(mm), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(mm))
234 0 : call addfld(src_names(m), (/ 'lev' /), 'A', 'kg/kg/s', trim(cnst_name(mm))//' source/sink')
235 :
236 0 : call add_default (cnst_name(mm), 1, ' ')
237 0 : call add_default (src_names(m), 1, ' ')
238 : end do
239 :
240 0 : call get_start_date(yr, mon, day, sec)
241 :
242 0 : ymd = yr*10000 + mon*100 + day
243 :
244 0 : yr0 = yr
245 0 : calday0 = get_calday(ymd, sec)
246 :
247 1536 : end subroutine aoa_tracers_init
248 :
249 : !===============================================================================
250 :
251 370944 : 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 370944 : if (.not. aoa_tracers_flag) return
267 :
268 0 : call get_curr_date (yr,mon,day,tod)
269 :
270 0 : if ( day == 1 .and. tod == 0) then
271 0 : if (masterproc) then
272 0 : write(iulog,*) 'AGE_OF_AIR_CONSTITUENTS: RE-INITIALIZING HORZ/VERT CONSTITUENTS'
273 : endif
274 :
275 0 : do c = begchunk, endchunk
276 0 : ncol = phys_state(c)%ncol
277 0 : do k = 1, pver
278 0 : do i = 1, ncol
279 0 : phys_state(c)%q(i,k,ixht) = 2._r8 + sin(phys_state(c)%lat(i))
280 0 : 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 0 : ymd = yr*10000 + mon*100 + day
288 0 : calday = get_calday(ymd, tod)
289 :
290 0 : dpy = 365._r8
291 0 : if (timemgr_get_calendar_cf() == 'gregorian' .and. is_leapyear(yr)) then
292 0 : dpy = 366._r8
293 : end if
294 0 : years = (yr-yr0) + (calday-calday0)/dpy
295 :
296 741888 : end subroutine aoa_tracers_timestep_init
297 :
298 : !===============================================================================
299 :
300 62545392 : subroutine aoa_tracers_timestep_tend(state, ptend, dt)
301 :
302 370944 : 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 : !------------------------------------------------------------------
328 :
329 1489176 : teul = .5_r8*dt/(86400._r8 * treldays) ! 1/2 for the semi-implicit scheme if dt=time step
330 1489176 : wimp = 1._r8/(1._r8 +teul)
331 1489176 : wsrc = teul*wimp
332 :
333 1489176 : if (.not. aoa_tracers_flag) then
334 1489176 : call physics_ptend_init(ptend,state%psetcols,'none') !Initialize an empty ptend for use with physics_update
335 : return
336 : end if
337 :
338 0 : lq(:) = .FALSE.
339 0 : lq(ixaoa) = .TRUE.
340 0 : lq(ixht) = .TRUE.
341 0 : lq(ixvt) = .TRUE.
342 :
343 0 : call physics_ptend_init(ptend,state%psetcols, 'aoa_tracers', lq=lq)
344 :
345 0 : lchnk = state%lchnk
346 0 : ncol = state%ncol
347 :
348 : ! AOAMF
349 0 : xmmr = mmr0*(1._r8 + per_yr*years)
350 0 : ptend%q(1:ncol,pver,ixaoa) = (xmmr - state%q(1:ncol,pver,ixaoa)) / dt
351 :
352 0 : do k = 1, pver
353 0 : do i = 1, ncol
354 :
355 : ! HORZ
356 0 : qrel = 2._r8 + sin(state%lat(i)) ! qrel should zonal mean
357 0 : xhorz = state%q(i,k,ixht)*wimp + wsrc*qrel ! Xnew = weight*3D-tracer + (1.-weight)*1D-tracer
358 0 : ptend%q(i,k,ixht) = (xhorz - state%q(i,k,ixht)) / dt ! Xnew = weight*3D-tracer + (1.-weight)*2D-tracer zonal mean
359 : ! Can be still used .... to diagnose fluxes OT-tracers
360 : ! VERT
361 0 : qrel = qrel_vert(k) ! qrel should zonal mean
362 0 : xvert = wimp*state%q(i,k,ixvt) + wsrc*qrel
363 0 : ptend%q(i,k,ixvt) = (xvert - state%q(i,k,ixvt)) / dt
364 :
365 : end do
366 :
367 : end do
368 :
369 : ! record tendencies on history files
370 0 : call outfld (src_names(1), ptend%q(:,:,ixaoa), pcols, lchnk)
371 0 : call outfld (src_names(2), ptend%q(:,:,ixht), pcols, lchnk)
372 0 : call outfld (src_names(3), ptend%q(:,:,ixvt), pcols, lchnk)
373 :
374 1489176 : end subroutine aoa_tracers_timestep_tend
375 :
376 : !===========================================================================
377 :
378 0 : subroutine init_cnst_3d(m, latvals, lonvals, mask, q)
379 :
380 : integer, intent(in) :: m ! global constituent index
381 : real(r8), intent(in) :: latvals(:) ! lat in degrees (ncol)
382 : real(r8), intent(in) :: lonvals(:) ! lon in degrees (ncol)
383 : logical, intent(in) :: mask(:) ! Only initialize where .true.
384 : real(r8), intent(out) :: q(:,:) ! kg tracer/kg dry air (gcol,plev)
385 :
386 : integer :: j, k, gsize
387 : !-----------------------------------------------------------------------
388 :
389 0 : if (masterproc) then
390 0 : write(iulog,*) 'AGE-OF-AIR CONSTITUENTS: INITIALIZING ',cnst_name(m),m
391 : end if
392 :
393 0 : if (m == ixaoa) then
394 :
395 : ! AOAMF
396 0 : q(:,:) = 0.0_r8
397 :
398 0 : else if (m == ixht) then
399 :
400 : ! HORZ
401 0 : gsize = size(q, 1)
402 0 : do j = 1, gsize
403 0 : q(j,:) = 2._r8 + sin(latvals(j))
404 : end do
405 :
406 0 : else if (m == ixvt) then
407 :
408 : ! VERT
409 0 : do k = 1, pver
410 0 : do j = 1, size(q,1)
411 0 : q(j,k) = qrel_vert(k)
412 : end do
413 : end do
414 :
415 : end if
416 :
417 1489176 : end subroutine init_cnst_3d
418 :
419 : !=====================================================================
420 :
421 : end module aoa_tracers
|
