Line data Source code
1 : !======================================================================
2 : ! Module implements passive test tracers.
3 : !
4 : ! Two options:
5 : !
6 : ! 1) Specify only the number of tracers desired by setting the -nadv_tt option to configure.
7 : ! This results in setting up the desired number of tracers making use of the tracers_suite
8 : ! module to generate tracer names and initialize mixing ratios.
9 : !
10 : ! 2) Specify both the number of tracers using configure's -nadv_tt option, and specify
11 : ! the same number of tracer names using the test_tracer_names namelist variable. This
12 : ! specifies a set of passive tracers that are either initialized from analytic expressions
13 : ! or by reading values from the IC file. The tracers for which analytic expressions are
14 : ! available are the following:
15 : !
16 : ! test_tracer_names description
17 : ! ----------------- -----------
18 : ! TT_SLOT Non-smooth scalar field (slotted cylinder)
19 : ! TT_SLOT1 1st slotted cylinder for 3 tracer test (Lauritzen and Thuburn, 2012, QJRMS)
20 : ! TT_SLOT2 2nd slotted cylinder for 3 tracer test (Lauritzen and Thuburn, 2012, QJRMS)
21 : ! TT_SLOT3 3rd slotted cylinder for 3 tracer test: TT_SLOT3=1-(TT_SLOT1+TT_SLOT2)
22 : ! TT_GBALL Smooth Gaussian "ball"
23 : ! TT_TANH Zonally constant, tanh function of latitude.
24 : ! TT_EM8 Constant field of size 1.e-8.
25 : ! TT_Y2_2 Approximately Y^2_2 spherical harmonic
26 : ! TT_Y32_16 Approximately Y^32_16 spherical harmonic
27 : ! TT_LATP2 Zonally constant, latitude + 2.
28 : ! TT_LONP2 Meridionally constant, cos(longitude) + 2.
29 : ! TT_COSB Cosine bell (Kent et al., 2012, MWR)
30 : ! TT_CCOSB Non-linearly correlated Cosine bell (Lauritzen and Thuburn,2012, QJRMS)
31 : ! TT_lCCOSB Linearly correlated Cosine bell
32 : !
33 : !======================================================================
34 :
35 : module tracers
36 :
37 : use shr_kind_mod, only: r8 => shr_kind_r8
38 : use shr_sys_mod, only: shr_sys_flush
39 : use spmd_utils, only: masterproc
40 : use ppgrid, only: pver
41 : use physconst, only: mwdry, cpair
42 : use constituents, only: cnst_add, cnst_name, cnst_longname
43 : use tracers_suite, only: get_tracer_name, init_cnst_tr
44 : use cam_history, only: addfld, add_default
45 : use cam_logfile, only: iulog
46 : use cam_abortutils, only: endrun
47 :
48 : implicit none
49 : private
50 : save
51 :
52 : public :: &
53 : tracers_readnl, &! read namelist
54 : tracers_register, &! register constituent
55 : tracers_implements_cnst, &! true if named constituent is implemented by this package
56 : tracers_init_cnst, &! initialize constituent field
57 : tracers_init ! initialize history fields, datasets
58 :
59 : integer, parameter :: num_names_max = 30
60 : integer, parameter :: num_analytic = 14
61 :
62 : ! Data from namelist variables
63 : integer :: test_tracer_num = 0
64 : character(len=16) :: test_tracer_names(num_names_max)
65 :
66 : logical :: tracers_flag = .false. ! true => turn on test tracer code
67 : logical :: tracers_suite_flag = .false. ! true => test tracers provided by tracers_suite module
68 :
69 : integer :: ixtrct=-999 ! index of 1st constituent
70 :
71 : character(len=16), parameter :: analytic_names(num_analytic) = &
72 : (/'TT_SLOT1 ', 'TT_SLOT2 ', 'TT_SLOT3 ', &
73 : 'TT_SLOT ', 'TT_GBALL ', 'TT_TANH ', &
74 : 'TT_EM8 ', 'TT_Y2_2 ', 'TT_Y32_16 ', &
75 : 'TT_LATP2 ', 'TT_LONP2 ', 'TT_COSB ', &
76 : 'TT_CCOSB ', 'TT_lCCOSB '/)
77 :
78 :
79 : logical :: analytic_tracer(num_names_max)
80 :
81 : !======================================================================
82 : contains
83 : !======================================================================
84 :
85 1536 : subroutine tracers_readnl(nlfile)
86 :
87 : use namelist_utils, only: find_group_name
88 : use units, only: getunit, freeunit
89 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_integer, mpi_character
90 :
91 : ! args
92 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
93 :
94 : ! Local variables
95 : integer :: unitn, ierr
96 : integer :: i, j
97 : integer :: num_names
98 : character(len=*), parameter :: subname = 'tracers_readnl'
99 :
100 : namelist /test_tracers_nl/ test_tracer_num, test_tracer_names
101 : !-----------------------------------------------------------------------------
102 :
103 47616 : test_tracer_names = (/ (' ', i=1,num_names_max) /)
104 :
105 1536 : if (masterproc) then
106 2 : unitn = getunit()
107 2 : open( unitn, file=trim(nlfile), status='old' )
108 2 : call find_group_name(unitn, 'test_tracers_nl', status=ierr)
109 2 : if (ierr == 0) then
110 0 : read(unitn, test_tracers_nl, iostat=ierr)
111 0 : if (ierr /= 0) then
112 0 : call endrun(subname // ':: ERROR reading namelist')
113 : end if
114 : end if
115 2 : close(unitn)
116 2 : call freeunit(unitn)
117 : end if
118 :
119 : call mpi_bcast(test_tracer_names, len(test_tracer_names)*num_names_max, mpi_character, &
120 1536 : mstrid, mpicom, ierr)
121 1536 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: test_tracer_names")
122 :
123 1536 : call mpi_bcast(test_tracer_num, 1, mpi_integer, mstrid, mpicom, ierr)
124 1536 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: iradsw")
125 :
126 : ! If any tracers have been specified then turn on the tracers module
127 1536 : if (test_tracer_num > 0) then
128 0 : tracers_flag = .true.
129 : else
130 1536 : return
131 : end if
132 :
133 : ! Determine the number of tracer names supplied:
134 0 : num_names = 0
135 0 : analytic_tracer = .false.
136 0 : do i = 1, num_names_max
137 0 : if (len_trim(test_tracer_names(i)) > 0) then
138 0 : num_names = num_names + 1
139 :
140 : ! Does the tracer have an analytic IC?
141 0 : do j = 1, num_analytic
142 0 : if (trim(test_tracer_names(i)) == trim(analytic_names(j))) then
143 0 : analytic_tracer(i) = .true.
144 0 : exit
145 : end if
146 : end do
147 : else
148 : exit
149 : end if
150 : end do
151 :
152 0 : if (num_names > 0) then
153 : ! If test_tracer_names have been specified, the test_tracer_num should
154 : ! equal the number of names supplied.
155 0 : if (num_names /= test_tracer_num) then
156 0 : write(iulog, *) subname//' number of names, number of tracers: ', num_names, test_tracer_num
157 0 : call endrun(subname // ':: number of names does not match number of tracers')
158 : end if
159 : else
160 : ! If no names have been supplied then
161 : ! the tracers will be provided by the tracers_suite module.
162 0 : tracers_suite_flag = .true.
163 : end if
164 :
165 : ! Print summary to log file
166 0 : if (masterproc) then
167 :
168 0 : write(iulog, *) 'Test Tracers Module'
169 0 : write(iulog, *) ' Number of Test Tracers:', test_tracer_num
170 0 : if (tracers_suite_flag) then
171 0 : write(iulog, *) ' Tracers will be provided by tracers_suite module.'
172 : else
173 0 : do i = 1, num_names
174 0 : if (analytic_tracer(i)) then
175 : write(iulog, *) ' '//trim(test_tracer_names(i))//&
176 0 : ' will be initialized from an analytic expression'
177 : else
178 : write(iulog, *) ' '//trim(test_tracer_names(i))//&
179 0 : ' will be initialized from the IC file'
180 : end if
181 : end do
182 : end if
183 : end if
184 :
185 : end subroutine tracers_readnl
186 :
187 : !======================================================================
188 :
189 1536 : subroutine tracers_register()
190 :
191 : ! Register advected tracers.
192 :
193 : ! Local variables
194 : integer :: m, mm
195 : logical :: read_from_file
196 : real(r8) :: minc
197 : character(len=16) :: name
198 : !-----------------------------------------------------------------------
199 :
200 1536 : if (.not. tracers_flag) return
201 :
202 0 : minc = -1.e36_r8 ! min mixing ratio (disable qneg3)
203 :
204 0 : do m = 1, test_tracer_num
205 :
206 0 : read_from_file = .true.
207 0 : if (tracers_suite_flag) then
208 0 : name = get_tracer_name(m) ! get name from suite file
209 0 : read_from_file = .false.
210 : else
211 0 : name = test_tracer_names(m)
212 0 : if (analytic_tracer(m)) read_from_file = .false.
213 : end if
214 :
215 : ! add constituent name to list of advected, save index number ixtrct
216 : call cnst_add(name, mwdry, cpair, minc, mm, &
217 0 : readiv=read_from_file, mixtype='dry')
218 0 : if (m == 1) ixtrct = mm ! save index number of first tracer
219 : end do
220 :
221 : end subroutine tracers_register
222 :
223 : !======================================================================
224 :
225 0 : function tracers_implements_cnst(name)
226 :
227 : ! return true if specified constituent is implemented by this package
228 :
229 : ! Arguments
230 : character(len=*), intent(in) :: name ! constituent name
231 : logical :: tracers_implements_cnst ! return value
232 :
233 : ! Local variables
234 : integer :: m
235 : character(len=16) :: trc_name
236 : !-----------------------------------------------------------------------
237 :
238 0 : tracers_implements_cnst = .false.
239 0 : if (.not. tracers_flag) return
240 :
241 0 : do m = 1, test_tracer_num
242 :
243 0 : if (tracers_suite_flag) then
244 0 : trc_name = get_tracer_name(m)
245 : else
246 0 : trc_name = test_tracer_names(m)
247 : end if
248 :
249 0 : if (name == trc_name) then
250 0 : tracers_implements_cnst = .true.
251 : return
252 : end if
253 : end do
254 :
255 : end function tracers_implements_cnst
256 :
257 : !===============================================================================
258 :
259 0 : subroutine tracers_init_cnst(name, latvals, lonvals, mask, q, z)
260 :
261 : ! Initialize test tracer mixing ratio
262 :
263 : character(len=*), intent(in) :: name
264 : real(r8), intent(in) :: latvals(:) ! lat in radians (ncol)
265 : real(r8), intent(in) :: lonvals(:) ! lon in radians (ncol)
266 : logical, intent(in) :: mask(:) ! Only initialize where .true.
267 : real(r8), intent(out) :: q(:,:) ! kg tracer/kg dry air (gcol, plev
268 : real(r8),optional,intent(in) :: z(:,:) ! height of full level pressure
269 :
270 : ! Local
271 : integer :: m
272 : logical :: found
273 : character(len=*), parameter :: subname = 'tracers_init_cnst'
274 : !----------------------------------------------------------------------------
275 :
276 0 : if (.not. tracers_flag) return
277 :
278 0 : found = .false.
279 0 : if (tracers_suite_flag) then
280 :
281 0 : do m = 1, test_tracer_num
282 0 : if (name == get_tracer_name(m)) then
283 0 : call init_cnst_tr(m, latvals, lonvals, mask, q)
284 0 : found = .true.
285 0 : exit
286 : endif
287 : end do
288 :
289 : else
290 :
291 0 : do m = 1, test_tracer_num
292 0 : if (name == test_tracer_names(m)) then
293 0 : if (analytic_tracer(m)) then
294 0 : if (present(z)) then
295 0 : call test_func_set(name, latvals, lonvals, mask, q, z=z)
296 : else
297 0 : call test_func_set(name, latvals, lonvals, mask, q)
298 : end if
299 : found = .true.
300 : exit
301 : else
302 : ! The initial values were supposed to be read from the IC file. This
303 : ! call should not have been made in that case, so it appears that a requested
304 : ! tracer is not on the IC file.
305 0 : write(iulog, *) subname//': ERROR: tracer ', trim(name), ' should be on IC file'
306 0 : call endrun(subname//': ERROR: tracer missing from IC file')
307 : end if
308 : end if
309 : end do
310 : end if
311 :
312 : if (.not. found) then
313 : ! unrecognized tracer name
314 0 : write(iulog, *) subname//': ERROR: ', trim(name), ' not recognized'
315 0 : call endrun(subname//': ERROR: tracer name not recognized')
316 : end if
317 :
318 : end subroutine tracers_init_cnst
319 :
320 : !===============================================================================
321 :
322 1536 : subroutine tracers_init()
323 :
324 : ! Add tracers to history output
325 :
326 : ! Local
327 : integer m, mm
328 : character(len=16) :: name ! constituent name
329 :
330 1536 : if (.not. tracers_flag ) return
331 :
332 0 : do m = 1,test_tracer_num
333 0 : mm = ixtrct + m - 1
334 0 : call addfld(cnst_name(mm), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(mm))
335 0 : call add_default(cnst_name(mm), 1, ' ')
336 : end do
337 : end subroutine tracers_init
338 :
339 : !=========================================================================================
340 :
341 0 : subroutine test_func_set(name, latvals, lonvals, mask, q, z)
342 :
343 : ! use test_func code to set array q
344 :
345 : character(len=*), intent(in) :: name
346 : real(r8), intent(in) :: latvals(:) ! lat in radians (ncol)
347 : real(r8), intent(in) :: lonvals(:) ! lon in radians (ncol)
348 : logical, intent(in) :: mask(:) ! Only initialize where .true.
349 : real(r8), intent(out) :: q(:,:) ! kg tracer/kg dry air (gcol, plev
350 : real(r8),optional,intent(in) :: z(:,:) ! height of full level pressure
351 :
352 : ! local variables
353 : integer :: i, k
354 : !----------------------------------------------------------------------------
355 :
356 0 : do i = 1, size(mask)
357 0 : if (mask(i)) then
358 0 : if (present(z)) then
359 0 : do k = 1, size(q, 2)
360 0 : q(i,k) = test_func(name, latvals(i), lonvals(i), k, z=z(i,k))
361 : end do
362 : else
363 0 : do k = 1, size(q, 2)
364 0 : q(i,k) = test_func(name, latvals(i), lonvals(i), k)
365 : end do
366 : end if
367 : end if
368 : end do
369 :
370 0 : end subroutine test_func_set
371 :
372 : !=========================================================================================
373 :
374 0 : function test_func(name, lat, lon, k, z) result(fout)
375 :
376 : ! Analytic test functions.
377 :
378 : use physconst, only: pi
379 : use hycoef, only: hyam, hybm, ps0
380 :
381 : character(len=*), intent(in) :: name
382 : real(r8), intent(in) :: lon ! radians
383 : real(r8), intent(in) :: lat ! radians
384 : integer, intent(in) :: k ! vertical index for layer mid-points
385 : real(r8),optional,intent(in) :: z ! height
386 :
387 : real(r8) :: fout
388 :
389 : real(r8), parameter :: psurf_moist = 100000.0_r8 ! moist surface pressure
390 : real(r8), parameter :: deg2rad = pi/180._r8
391 : real(r8), parameter :: z0 = 4500.0_r8 !Eq 9 in Kent et al. (2012)
392 : real(r8), parameter :: zz = 1000.0_r8 !Eq 9 in Kent et al. (2012)
393 :
394 : real(r8) :: lon1, lat1, R0, Rg1
395 : real(r8) :: eta, eta_c, d1, f1, f2, c
396 : real(r8) :: cos_tmp, sin_tmp
397 :
398 : !----------------------------------------------------------------------------
399 :
400 0 : select case(name)
401 : case('TT_SLOT')
402 : !
403 : ! Non-smooth scalar field (slotted cylinder): min=0.1 and max=1
404 : !
405 0 : R0 = 0.25_r8
406 0 : lon1 = pi/9.0_r8
407 0 : lat1 = 2.0_r8*pi/9.0_r8
408 0 : Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
409 0 : c = 1.0_r8
410 :
411 0 : if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
412 : fout = c
413 : elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
414 0 : .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
415 : fout = c
416 : else
417 0 : fout = 0.1_r8
418 : endif
419 : case('TT_SLOT1')
420 : !
421 : ! Slotted cylinder for 3 tracer test
422 : !
423 0 : R0 = 0.25_r8
424 0 : lon1 = pi/9.0_r8
425 0 : lat1 = 2.0_r8*pi/9.0_r8
426 0 : Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
427 0 : c = 1.0_r8/3.0_r8
428 :
429 0 : if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
430 : fout = c
431 : elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
432 0 : .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
433 : fout = c
434 : else
435 0 : fout = 0.1_r8
436 : endif
437 : case('TT_SLOT2')
438 : !
439 : ! Slotted cylinder for 3 tracer test
440 : !
441 0 : R0 = 0.25_r8
442 0 : lon1 = pi/9.0_r8
443 0 : lat1 = 2.0_r8*pi/9.0_r8+pi/18.0_r8
444 0 : Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
445 0 : c = 2.0_r8/3.0_r8
446 :
447 0 : if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
448 : fout = c
449 : elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
450 0 : .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
451 : fout = c
452 : else
453 0 : fout = 0.1_r8
454 : endif
455 : case('TT_SLOT3')
456 : !
457 : ! Slotted cylinder for 3 tracer test
458 : ! Residual of two slotted cylinders so that TT_SLOT1+TT_SLOT2+TT_SLOT3=1
459 : ! Lauritzen and Thuburn (2012, QJRMS)
460 : !
461 0 : R0 = 0.25_r8
462 0 : lon1 = pi/9.0_r8
463 0 : lat1 = 2.0_r8*pi/9.0_r8
464 0 : Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
465 0 : c = 1.0_r8/3.0_r8
466 :
467 0 : if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
468 0 : f1 = c
469 : elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
470 0 : .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
471 0 : f1 = c
472 : else
473 0 : f1 = 0.1_r8
474 : endif
475 :
476 0 : R0 = 0.25_r8
477 0 : lon1 = pi/9.0_r8
478 0 : lat1 = 2.0_r8*pi/9.0_r8+pi/18.0_r8
479 0 : Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
480 0 : c = 2.0_r8/3.0_r8
481 0 : if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
482 : f2 = c
483 : elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
484 0 : .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
485 : f2 = c
486 : else
487 0 : f2 = 0.1_r8
488 : endif
489 :
490 0 : fout = 1.0_r8-(f1+f2)
491 : case('TT_GBALL')
492 : !
493 : ! Smooth Gaussian "ball"
494 : !
495 0 : R0 = 10.0_r8 ! radius of the perturbation
496 0 : lon1 = 20.0_r8*deg2rad ! longitudinal position, 20E
497 0 : lat1 = 40.0_r8*deg2rad ! latitudinal position, 40N
498 0 : eta_c = 0.6_r8
499 0 : sin_tmp = SIN(lat1)*SIN(lat)
500 0 : cos_tmp = COS(lat1)*COS(lat)
501 :
502 0 : Rg1 = ACOS( sin_tmp + cos_tmp*COS(lon-lon1) ) ! great circle distance
503 0 : eta = (hyam(k)*ps0 + hybm(k)*psurf_moist)/psurf_moist
504 0 : fout = EXP(- ((Rg1*R0)**2 + ((eta-eta_c)/0.1_r8)**2))
505 0 : IF (ABS(fout) < 1.0E-8_r8) fout = 0.0_r8
506 :
507 : case('TT_TANH')
508 : !
509 : !
510 : !
511 0 : fout = 0.5_r8 * ( tanh( 3.0_r8*abs(lat)-pi ) + 1.0_r8)
512 :
513 : case('TT_EM8')
514 0 : fout = 1.0e-8_r8
515 :
516 : case('TT_Y2_2')
517 : !
518 : ! approximately Y^2_2 spherical harmonic
519 : !
520 0 : fout = 0.5_r8 + 0.5_r8*(cos(lat)*cos(lat)*cos(2.0_r8*lon))
521 :
522 : case('TT_Y32_16')
523 : !
524 : ! approximately Y32_16 spherical harmonic
525 : !
526 0 : fout = 0.5_r8 + 0.5_r8*(cos(16*lon)*(sin(2_r8*lat)**16))
527 :
528 : case('TT_LATP2')
529 0 : fout = 2.0_r8 + lat
530 :
531 : case('TT_LONP2')
532 0 : fout = 2.0_r8 + cos(lon)
533 : case('TT_COSB')
534 : !
535 : ! Cosine bell inspired by Kent et al., 2012, MWR; only one bell and location changed
536 : ! https://journals.ametsoc.org/doi/pdf/10.1175/MWR-D-11-00150.1
537 : !
538 0 : R0 = 0.5_r8 ! radius of the perturbation
539 0 : lon1 = pi/9.0_r8
540 0 : lat1 = 2.0_r8*pi/9.0_r8
541 :
542 0 : sin_tmp = SIN(lat1)*SIN(lat)
543 0 : cos_tmp = COS(lat1)*COS(lat)
544 0 : Rg1 = ACOS( sin_tmp + cos_tmp*COS(lon-lon1) ) ! great circle distance
545 :
546 0 : if (present(z)) THEN
547 0 : d1 = MIN(1.0_r8,(Rg1/R0)**2+((z-z0)/zz)**2)
548 : else
549 0 : d1 = MIN(1.0_r8,(Rg1/R0)**2)
550 : end if
551 :
552 0 : if (Rg1 < R0) then
553 0 : fout = 0.1_r8+0.9_r8*0.5_r8*(1.0_r8+COS(pi*d1))
554 : else
555 : fout = 0.1_r8
556 : end if
557 : case('TT_CCOSB')
558 : !
559 : ! Correlated cosine bell
560 : !
561 0 : R0 = 0.5_r8 ! radius of the perturbation
562 0 : lon1 = pi/9.0_r8
563 0 : lat1 = 2.0_r8*pi/9.0_r8
564 :
565 0 : sin_tmp = SIN(lat1)*SIN(lat)
566 0 : cos_tmp = COS(lat1)*COS(lat)
567 0 : Rg1 = ACOS( sin_tmp + cos_tmp*COS(lon-lon1) ) ! great circle distance
568 :
569 0 : if (present(z)) THEN
570 0 : d1 = MIN(1.0_r8,(Rg1/R0)**2+((z-z0)/zz)**2)
571 : else
572 0 : d1 = MIN(1.0_r8,(Rg1/R0)**2)
573 : end if
574 :
575 0 : if (Rg1 < R0) then
576 0 : f1 = 0.1_r8+0.9_r8*0.5_r8*(1.0_r8+COS(pi*d1))
577 : else
578 0 : f1 = 0.1_r8
579 : end if
580 :
581 0 : fout=corr_fct(f1)
582 : case('TT_lCCOSB')
583 : !
584 : ! Correlated cosine bell
585 : !
586 0 : R0 = 0.5_r8 ! radius of the perturbation
587 0 : lon1 = pi/9.0_r8
588 0 : lat1 = 2.0_r8*pi/9.0_r8
589 :
590 0 : sin_tmp = SIN(lat1)*SIN(lat)
591 0 : cos_tmp = COS(lat1)*COS(lat)
592 0 : Rg1 = ACOS( sin_tmp + cos_tmp*COS(lon-lon1) ) ! great circle distance
593 :
594 0 : if (present(z)) THEN
595 0 : d1 = MIN(1.0_r8,(Rg1/R0)**2+((z-z0)/zz)**2)
596 : else
597 0 : d1 = MIN(1.0_r8,(Rg1/R0)**2)
598 : end if
599 :
600 0 : if (Rg1 < R0) then
601 0 : f1 = 0.1_r8+0.9_r8*0.5_r8*(1.0_r8+COS(pi*d1))
602 : else
603 0 : f1 = 0.1_r8
604 : end if
605 :
606 0 : fout=1.0_r8-f1
607 :
608 : case default
609 0 : call endrun("test_func: ERROR: name not recognized")
610 : end select
611 0 : end function test_func
612 :
613 : !
614 : ! correlation function (Lauritzen and Thuburn, 2012, QJRMS)
615 : !
616 0 : function corr_fct(x) result(out)
617 : real(r8), intent(in) :: x
618 : real(r8) :: out
619 0 : out = -0.8_r8*x**2+0.9_r8
620 0 : end function corr_fct
621 : end module tracers
|
