Line data Source code
1 : !-------------------------------------------------------------------
2 : ! Manages reading and interpolation of prescribed aerosol concentrations.
3 : ! This places the concentration fields in the physics buffer so that
4 : ! radiation package can access them.
5 : !
6 : ! This has been generalized so that the field names in the data files
7 : ! and the field names in the physics buffer can be arbitrary.
8 : !
9 : ! The prescribed_aero_specifier namelist variable specifies a list of
10 : ! variable names of the concentration fields in the netCDF dataset (ncdf_fld_name)
11 : ! and the corresponding names used in the physics buffer:
12 : !
13 : ! prescribed_aero_specifier = 'pbuf_name1:ncdf_fld_name1','pbuf_name2:ncdf_fld_name2', ...
14 : !
15 : ! If there is no ":" then the specified name is used as both the
16 : ! pbuf_name and ncdf_fld_name
17 : !
18 : ! Created by: Francis Vitt
19 : !-------------------------------------------------------------------
20 : module prescribed_aero
21 :
22 : use shr_kind_mod, only : r8 => shr_kind_r8
23 : use cam_abortutils, only : endrun
24 : use spmd_utils, only : masterproc
25 : use tracer_data, only : trfld, trfile
26 : use cam_logfile, only : iulog
27 : use pio, only : var_desc_t
28 :
29 : implicit none
30 : private
31 : save
32 :
33 : type(trfld), pointer :: fields(:)
34 : type(trfile) :: file
35 :
36 : public :: prescribed_aero_init
37 : public :: prescribed_aero_adv
38 : public :: write_prescribed_aero_restart
39 : public :: read_prescribed_aero_restart
40 : public :: has_prescribed_aero
41 : public :: prescribed_aero_register
42 : public :: init_prescribed_aero_restart
43 : public :: prescribed_aero_readnl
44 :
45 : logical :: has_prescribed_aero = .false.
46 :
47 : ! Decides if its a modal aerosol simulation or not
48 : logical :: clim_modal_aero = .false.
49 :
50 : integer, parameter, public :: N_AERO = 50
51 :
52 : integer :: number_flds
53 :
54 : character(len=256) :: filename = 'NONE'
55 : character(len=256) :: filelist = ' '
56 : character(len=256) :: datapath = ' '
57 : character(len=32) :: datatype = 'SERIAL'
58 : logical :: rmv_file = .false.
59 : integer :: cycle_yr = 0
60 : integer :: fixed_ymd = 0
61 : integer :: fixed_tod = 0
62 : character(len=32) :: specifier(N_AERO) = ''
63 :
64 : ! prescribed aerosol names
65 : character(len=16), allocatable :: pbuf_names(:)
66 :
67 : integer :: aero_cnt
68 : integer :: aero_cnt_c = 0 ! clound borne species count for modal aerosols (zero otherwise)
69 :
70 : ! Normal random number which persists from one tiemstep to the next
71 : real(r8) :: randn_persists
72 :
73 : ! Following definitions are added to
74 : ! allow randn_persists to persist during restart runs
75 : type(var_desc_t) :: randn_persists_desc
76 : character(len=16), parameter :: randn_persists_name = 'prescraero_randn'
77 :
78 : contains
79 :
80 : !-------------------------------------------------------------------
81 : ! registers aerosol fields to the phys buffer
82 : !-------------------------------------------------------------------
83 2304 : subroutine prescribed_aero_register()
84 :
85 : use ppgrid, only: pver,pcols
86 :
87 : use physics_buffer, only : pbuf_add_field, dtype_r8
88 : integer :: i,idx
89 :
90 1536 : if (has_prescribed_aero) then
91 0 : do i = 1,aero_cnt
92 0 : call pbuf_add_field(pbuf_names(i),'physpkg',dtype_r8,(/pcols,pver/),idx)
93 : enddo
94 : endif
95 :
96 1536 : endsubroutine prescribed_aero_register
97 :
98 : !-------------------------------------------------------------------
99 : ! parses prescribed_aero_specifier namelist option
100 : !-------------------------------------------------------------------
101 1536 : subroutine prescribed_aero_init()
102 :
103 1536 : use tracer_data, only : trcdata_init
104 : use cam_history, only : addfld
105 : use ppgrid, only : pver
106 : use error_messages, only: handle_err
107 : use ppgrid, only: pcols, pver, begchunk, endchunk
108 : use physics_buffer, only : physics_buffer_desc
109 :
110 : implicit none
111 :
112 : ! local vars
113 : character(len=32) :: spec_a
114 : integer :: ndx, istat, i, i_c, j
115 :
116 1536 : if ( has_prescribed_aero ) then
117 0 : if ( masterproc ) then
118 0 : write(iulog,*) 'aero is prescribed in :'//trim(filename)
119 : endif
120 : else
121 : return
122 : endif
123 :
124 :
125 0 : allocate (file%in_pbuf(size(specifier)))
126 0 : if (clim_modal_aero) then
127 0 : file%in_pbuf(:) = .false.
128 0 : do i = 1,N_AERO
129 0 : do j=1,aero_cnt
130 0 : if(specifier(i) .eq. pbuf_names(j)) then
131 0 : file%in_pbuf(i) = .true.
132 0 : exit
133 : endif
134 : enddo
135 : enddo
136 : else
137 0 : file%in_pbuf(:) = .true.
138 : endif
139 : call trcdata_init( specifier, filename, filelist, datapath, fields, file, &
140 0 : rmv_file, cycle_yr, fixed_ymd, fixed_tod, datatype)
141 :
142 0 : number_flds = 0
143 0 : if (associated(fields)) number_flds = size( fields )
144 :
145 0 : if( number_flds < 1 ) then
146 0 : if ( masterproc ) then
147 0 : write(iulog,*) 'There are no prescribed aerosols'
148 0 : write(iulog,*) ' '
149 : endif
150 0 : return
151 : end if
152 :
153 : ! Following loop add fields for output. For modal aerosols, first the cld borne aersols
154 : ! are added and then their interstitial counterparts are added. The loop exits once all the cld borne
155 : ! aerosols (and their interstitial counterparts) are added. Note that the units(fields(i)%units),
156 : ! will be same for both interstitial and cloud borne species.
157 : ! All other aerosol treatments(bulk) are left unchanged
158 : i_c = 0
159 0 : fldloop:do i = 1,number_flds
160 0 : if(clim_modal_aero .and. index(trim(fields(i)%fldnam),'_c') > 1) then ! Only cloud borne species
161 0 : call addfld(trim(fields(i)%fldnam)//'_D', (/ 'lev' /), 'A',trim(fields(i)%units), 'prescribed aero' )
162 0 : call spec_c_to_a(trim(fields(i)%fldnam),spec_a)
163 0 : call addfld(trim(spec_a)//'_D', (/ 'lev' /), 'A',trim(fields(i)%units), 'prescribed aero' )
164 0 : i_c = i_c + 1
165 0 : if(i_c >= aero_cnt_c) exit fldloop
166 : else
167 0 : call addfld(trim(fields(i)%fldnam)//'_D', (/ 'lev' /), 'A',trim(fields(i)%units), 'prescribed aero' )
168 : endif
169 : enddo fldloop
170 :
171 1536 : end subroutine prescribed_aero_init
172 :
173 : !-------------------------------------------------------------------
174 : ! reads namelist options
175 : !-------------------------------------------------------------------
176 1536 : subroutine prescribed_aero_readnl(nlfile)
177 :
178 1536 : use namelist_utils, only: find_group_name
179 : use units, only: getunit, freeunit
180 : use mpishorthand
181 : use rad_constituents, only: rad_cnst_get_info ! Added to query if it is a modal aero sim or not
182 :
183 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
184 :
185 : ! Local variables
186 : integer :: unitn, ierr, nmodes, aero_loop_end
187 : logical :: skip_spec
188 : character(len=*), parameter :: subname = 'prescribed_aero_readnl'
189 :
190 : character(len=32) :: prescribed_aero_specifier(N_AERO)
191 : character(len=256) :: prescribed_aero_file
192 : character(len=256) :: prescribed_aero_filelist
193 : character(len=256) :: prescribed_aero_datapath
194 : character(len=32) :: prescribed_aero_type
195 : logical :: prescribed_aero_rmfile
196 : integer :: prescribed_aero_cycle_yr
197 : integer :: prescribed_aero_fixed_ymd
198 : integer :: prescribed_aero_fixed_tod
199 : integer :: i,k
200 :
201 : namelist /prescribed_aero_nl/ &
202 : prescribed_aero_specifier, &
203 : prescribed_aero_file, &
204 : prescribed_aero_filelist, &
205 : prescribed_aero_datapath, &
206 : prescribed_aero_type, &
207 : prescribed_aero_rmfile, &
208 : prescribed_aero_cycle_yr, &
209 : prescribed_aero_fixed_ymd, &
210 : prescribed_aero_fixed_tod
211 : !-----------------------------------------------------------------------------
212 :
213 : ! Initialize namelist variables from local module variables.
214 78336 : prescribed_aero_specifier= specifier
215 1536 : prescribed_aero_file = filename
216 1536 : prescribed_aero_filelist = filelist
217 1536 : prescribed_aero_datapath = datapath
218 1536 : prescribed_aero_type = datatype
219 1536 : prescribed_aero_rmfile = rmv_file
220 1536 : prescribed_aero_cycle_yr = cycle_yr
221 1536 : prescribed_aero_fixed_ymd= fixed_ymd
222 1536 : prescribed_aero_fixed_tod= fixed_tod
223 :
224 : ! Read namelist
225 1536 : if (masterproc) then
226 2 : unitn = getunit()
227 2 : open( unitn, file=trim(nlfile), status='old' )
228 2 : call find_group_name(unitn, 'prescribed_aero_nl', status=ierr)
229 2 : if (ierr == 0) then
230 2 : read(unitn, prescribed_aero_nl, iostat=ierr)
231 2 : if (ierr /= 0) then
232 0 : call endrun(subname // ':: ERROR reading namelist')
233 : end if
234 : end if
235 2 : close(unitn)
236 2 : call freeunit(unitn)
237 : end if
238 :
239 : #ifdef SPMD
240 : ! Broadcast namelist variables
241 1536 : call mpibcast(prescribed_aero_specifier,len(prescribed_aero_specifier(1))*N_AERO, mpichar, 0, mpicom)
242 1536 : call mpibcast(prescribed_aero_file, len(prescribed_aero_file), mpichar, 0, mpicom)
243 1536 : call mpibcast(prescribed_aero_filelist, len(prescribed_aero_filelist), mpichar, 0, mpicom)
244 1536 : call mpibcast(prescribed_aero_datapath, len(prescribed_aero_datapath), mpichar, 0, mpicom)
245 1536 : call mpibcast(prescribed_aero_type, len(prescribed_aero_type), mpichar, 0, mpicom)
246 1536 : call mpibcast(prescribed_aero_rmfile, 1, mpilog, 0, mpicom)
247 1536 : call mpibcast(prescribed_aero_cycle_yr, 1, mpiint, 0, mpicom)
248 1536 : call mpibcast(prescribed_aero_fixed_ymd,1, mpiint, 0, mpicom)
249 1536 : call mpibcast(prescribed_aero_fixed_tod,1, mpiint, 0, mpicom)
250 : #endif
251 :
252 : ! Update module variables with user settings.
253 78336 : specifier = prescribed_aero_specifier
254 1536 : filename = prescribed_aero_file
255 1536 : filelist = prescribed_aero_filelist
256 1536 : datapath = prescribed_aero_datapath
257 1536 : datatype = prescribed_aero_type
258 1536 : rmv_file = prescribed_aero_rmfile
259 1536 : cycle_yr = prescribed_aero_cycle_yr
260 1536 : fixed_ymd = prescribed_aero_fixed_ymd
261 1536 : fixed_tod = prescribed_aero_fixed_tod
262 :
263 : ! Turn on prescribed aerosols if user has specified an input dataset.
264 1536 : has_prescribed_aero = len_trim(filename) > 0 .and. filename.ne.'NONE'
265 :
266 1536 : if ( .not. has_prescribed_aero) return
267 :
268 : ! Determine whether its a 'modal' aerosol simulation or not
269 0 : call rad_cnst_get_info(0, nmodes=nmodes)
270 0 : clim_modal_aero = (nmodes > 0)
271 :
272 : ! For modal aerosols, interstitial species(*_a) are diagnosed from
273 : ! their *_logm and *_logv counterparts (e.g. soa_a1 is diagnosed from
274 : ! soa_a1_logm and soa_a1_logv). Therefore, only *_logm and *_logv and cloud
275 : ! borne (*_c) species are specified in the build-namelist. In the following
276 : ! cnt_loop, we will count the cloud borne species and *_logm species(in lieu
277 : ! of *_a species). We will skip *_logv species. This will ensure that aero_cnt
278 : ! variable is the sum of cloud borne and interstitial species (which we will
279 : ! manually add in pbuf_names array later). We are also counting cloud borne
280 : ! (*_c) species which will help adding the same number of interstitial species
281 : ! in pbuf_names array
282 :
283 : ! determine which prescibed aerosols are specified
284 0 : aero_cnt = 0
285 0 : aero_cnt_c = 0 ! cloud borne species count
286 0 : cnt_loop: do i = 1,N_AERO
287 0 : if ( len_trim(specifier(i))==0 ) exit cnt_loop
288 0 : skip_spec = .FALSE.
289 0 : if(clim_modal_aero) then
290 : ! For modal aerosols, skip counting species ending with *_logv
291 0 : if(index(specifier(i),'_c') >= 1) aero_cnt_c = aero_cnt_c + 1
292 0 : if(index(specifier(i),'_logv') >= 1) skip_spec = .TRUE.
293 : endif
294 0 : if(.NOT.skip_spec)aero_cnt = aero_cnt+1
295 : end do cnt_loop
296 :
297 0 : has_prescribed_aero = aero_cnt>0
298 0 : if ( .not. has_prescribed_aero) return
299 :
300 0 : allocate(pbuf_names(aero_cnt))
301 :
302 : ! For modal aerosols, the following loop will add the cloud borne
303 : ! species directly and interstitial species through the "add_interstitial_spec"
304 : ! call. Interstitial species are added at the end of the cloud borne species in
305 : ! pbuf_names array.
306 : ! Note that aero_cnt_c should be zero for all other aerosol trearments
307 : ! except the modal aerosols (e.g. bulk)
308 :
309 0 : if(.NOT. clim_modal_aero) aero_cnt_c = 0
310 0 : aero_loop_end = aero_cnt - aero_cnt_c
311 :
312 0 : do i = 1,aero_loop_end
313 0 : k = scan( specifier(i),':')
314 0 : if (k>1) then
315 0 : pbuf_names(i) = trim(specifier(i)(:k-1))
316 0 : if(clim_modal_aero)call add_interstitial_spec(aero_loop_end,i)
317 : else
318 0 : pbuf_names(i) = trim(specifier(i))
319 0 : if(clim_modal_aero)call add_interstitial_spec(aero_loop_end,i)
320 : endif
321 : enddo
322 :
323 1536 : end subroutine prescribed_aero_readnl
324 :
325 : !-------------------------------------------------------------------
326 : ! Add interstitial aerosols in pbuf_names array for modal aerosols
327 : !-------------------------------------------------------------------
328 0 : subroutine add_interstitial_spec(aero_loop_end,i_in)
329 : implicit none
330 :
331 : !Arguments
332 : integer, intent(in) :: i_in, aero_loop_end
333 :
334 : !Local
335 : character(len=32) :: spec_a
336 :
337 : ! Replace 'c' with 'a' in species name
338 0 : call spec_c_to_a(pbuf_names(i_in), spec_a)
339 0 : pbuf_names(aero_loop_end+i_in) = spec_a
340 1536 : end subroutine add_interstitial_spec
341 :
342 : !-------------------------------------------------------------------
343 : ! A generic subroutine which replaces 'c' in the cloud borne
344 : ! species name with 'a' to make it interstital species
345 : !-------------------------------------------------------------------
346 0 : subroutine spec_c_to_a (spec_c_in,spec_a_out)
347 : implicit none
348 :
349 : !Arguments
350 : character(len=*), intent(in) :: spec_c_in
351 : character(len=*), intent(out) :: spec_a_out
352 :
353 : !Local
354 : character(len=32) :: name
355 : character(len=1000) :: errMsg
356 : integer :: k_c, k_cp1
357 :
358 0 : k_c = index(trim(adjustl(spec_c_in)),'_c')
359 0 : if(k_c >= 1) then
360 0 : k_cp1 = k_c + 1
361 0 : name = trim(adjustl(spec_c_in))
362 0 : name(k_cp1:k_cp1) = 'a'
363 0 : spec_a_out = trim(adjustl(name))
364 : else
365 0 : write(errMsg,*) "Input string (",trim(spec_c_in)," is not a cld borne aerosol,", &
366 0 : "cannot form interstitial species name"
367 0 : call endrun(trim(errMsg))
368 : endif
369 0 : end subroutine spec_c_to_a
370 :
371 : !-------------------------------------------------------------------
372 : ! advances the aerosol fields to the current time step
373 : !-------------------------------------------------------------------
374 741888 : subroutine prescribed_aero_adv( state, pbuf2d )
375 :
376 : use tracer_data, only : advance_trcdata
377 : use physics_types,only : physics_state
378 : use ppgrid, only : begchunk, endchunk
379 : use ppgrid, only : pcols, pver
380 : use string_utils, only : to_lower, GLC
381 : use cam_history, only : outfld
382 :
383 : use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk, &
384 : pbuf_get_index
385 :
386 : implicit none
387 :
388 : type(physics_state), intent(in) :: state(begchunk:endchunk)
389 :
390 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
391 :
392 370944 : type(physics_buffer_desc), pointer :: pbuf_chnk(:)
393 :
394 : character(len=32) :: spec_a
395 : integer :: c,ncol, i, i_c, spec_a_ndx, errcode
396 370944 : real(r8),pointer :: outdata(:,:)
397 : logical :: cld_borne_aero = .FALSE.
398 :
399 370944 : if( .not. has_prescribed_aero ) return
400 :
401 0 : call advance_trcdata( fields, file, state, pbuf2d )
402 :
403 : ! Diagnose interstital species using random sampling
404 0 : if ( clim_modal_aero ) then
405 0 : call rand_sample_prescribed_aero(state, pbuf2d )
406 : endif
407 :
408 : ! set the tracer fields with the correct units
409 0 : i_c = 0
410 0 : fldloop : do i = 1,number_flds
411 :
412 : !$OMP PARALLEL DO PRIVATE (C, NCOL, OUTDATA,PBUF_CHNK)
413 0 : do c = begchunk,endchunk
414 0 : ncol = state(c)%ncol
415 0 : pbuf_chnk => pbuf_get_chunk(pbuf2d, c)
416 0 : if(clim_modal_aero .and. index(trim(fields(i)%fldnam),'_c') > 1) then ! Only cloud borne species
417 0 : call pbuf_get_field(pbuf_chnk, fields(i)%pbuf_ndx, outdata)
418 0 : call outfld( trim(fields(i)%fldnam)//'_D', outdata(:ncol,:), ncol, state(c)%lchnk )
419 :
420 0 : call spec_c_to_a(trim(fields(i)%fldnam),spec_a)
421 0 : spec_a_ndx = pbuf_get_index(trim(fields(i)%fldnam),errcode)
422 0 : call pbuf_get_field(pbuf_chnk, spec_a_ndx, outdata)
423 0 : call outfld( trim(spec_a)//'_D', outdata(:ncol,:), ncol, state(c)%lchnk )
424 0 : cld_borne_aero = .TRUE.
425 : else
426 0 : call pbuf_get_field(pbuf_chnk, fields(i)%pbuf_ndx, outdata)
427 0 : call outfld( trim(fields(i)%fldnam)//'_D', outdata(:ncol,:), ncol, state(c)%lchnk )
428 : endif
429 : enddo
430 0 : if(cld_borne_aero)then
431 0 : i_c = i_c + 1
432 0 : cld_borne_aero = .FALSE.
433 0 : if(i_c >= aero_cnt_c) exit fldloop
434 : endif
435 : enddo fldloop
436 :
437 370944 : end subroutine prescribed_aero_adv
438 :
439 : !-------------------------------------------------------------------
440 0 : subroutine rand_sample_prescribed_aero(state, pbuf2d)
441 :
442 : !Purpose: Generates log normal distribution for the interstitial species.
443 : !Note that only the interstitial aerosols are diagnosed here
444 : !
445 : !Written by Balwinder Singh (12/14/2012)
446 : !Adapted from Jin-Ho Yoon
447 : !
448 : !Update log:
449 :
450 370944 : use physics_types, only : physics_state
451 : use ppgrid, only : begchunk, endchunk, pver
452 : use physics_buffer, only : physics_buffer_desc, pbuf_get_field, pbuf_get_chunk, &
453 : pbuf_get_index
454 : use ppgrid, only : pcols, pver
455 :
456 : !Arguments
457 : type(physics_state), intent(in) :: state(begchunk:endchunk)
458 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
459 :
460 : !Local
461 : real(r8), parameter :: mean_max_val = 5.0_r8
462 : real(r8), parameter :: std_max_val = 3.0_r8
463 :
464 : integer :: c, ncol, i, kc, klog, logv_ndx, logm_ndx, a_ndx, icol, kpver
465 : real(r8) :: logm2, variance, std, mean_max, std_max, randn
466 0 : type(physics_buffer_desc), pointer :: pbuf_chnk(:)
467 0 : real(r8), pointer :: data(:,:)
468 : real(r8) :: logm_data(pcols,pver),logv_data(pcols,pver)
469 :
470 0 : randn = randn_prescribed_aero()
471 0 : do i = 1, aero_cnt
472 : !Species with '_a' are updated using random sampling.
473 : !Cloud borne species (ends with _c1,_c2, etc.) have to be skipped
474 0 : kc = index(pbuf_names(i),'_c')
475 : ! if kc>1, species is cloud borne species
476 0 : if( kc > 1) cycle
477 :
478 : !If species is ending with _a1, a2 etc., then find indicies of their
479 : !logv and lom conterparts
480 0 : logv_ndx = logvm_get_index(pbuf_names(i),'v')
481 0 : logm_ndx = logvm_get_index(pbuf_names(i),'m')
482 0 : a_ndx = pbuf_get_index(trim(adjustl(pbuf_names(i))))
483 :
484 0 : do c = begchunk, endchunk
485 0 : ncol = state(c)%ncol
486 0 : pbuf_chnk => pbuf_get_chunk(pbuf2d, c)
487 :
488 : !Get the species mixing ratio nad its logv and logm values
489 0 : call pbuf_get_field(pbuf_chnk, a_ndx, data)
490 0 : logv_data = fields(logv_ndx)%data(:,:,c)
491 0 : logm_data = fields(logm_ndx)%data(:,:,c)
492 :
493 0 : do icol = 1, ncol
494 0 : do kpver = 1, pver
495 0 : logm2 = logm_data(icol,kpver) * logm_data(icol,kpver)
496 :
497 : !Compute variance
498 0 : variance = max( 0.0_r8, (logv_data(icol,kpver) - logm2) )
499 :
500 : !Standard deviation
501 0 : std = sqrt(variance)
502 :
503 : !Bounds to keep mixing ratios from going unphysical
504 0 : mean_max = exp(logm_data(icol,kpver)) * mean_max_val
505 0 : std_max = exp(logm_data(icol,kpver) + std_max_val * std )
506 :
507 0 : data(icol,kpver) = min(exp(logm_data(icol,kpver)+randn*std),mean_max,std_max)
508 :
509 : enddo !pver
510 : enddo !col
511 : enddo !chunk
512 : enddo !flds
513 :
514 0 : end subroutine rand_sample_prescribed_aero
515 : !-------------------------------------------------------------------
516 0 : function randn_prescribed_aero()
517 : !Pupose: This function generates a new normally distributed random
518 : ! number at end end of each day. This random number then stays the same
519 : ! for the whole day
520 : !
521 : !Written by Balwinder Singh (12/14/2012)
522 : !Adapted from Jin-Ho Yoon
523 : !
524 : !Update log:
525 :
526 0 : use time_manager, only : is_end_curr_day, is_first_step, get_nstep
527 :
528 : integer, parameter :: rconst1_1 = 5000000
529 : integer, parameter :: rconst1_2 = 50
530 : integer, parameter :: rconst2_1 = 10000000
531 : integer, parameter :: rconst2_2 = 10
532 :
533 : integer :: i, seed_size, nstep
534 0 : integer, allocatable :: seed(:)
535 :
536 : real(r8) :: randn_prescribed_aero
537 : real(r8) :: randu1, randu2
538 :
539 : !Use same random number for the entire day and generate a new normally
540 : !distributed random number at the start of the new day
541 0 : if(is_first_step() .or. is_end_curr_day()) then
542 : !Generate two uniformly distributed random numbers (between 0 and 1)
543 0 : call random_seed(size=seed_size)
544 0 : allocate(seed(seed_size))
545 :
546 : ! Using nstep as a seed to generate same sequence
547 0 : nstep = get_nstep()
548 0 : do i = 1, seed_size
549 0 : seed(i) = rconst1_1*nstep + rconst1_2*(i-1)
550 : end do
551 0 : call random_seed(put=seed)
552 0 : call random_number (randu1)
553 :
554 0 : do i = 1, seed_size
555 0 : seed(i) = rconst2_1*nstep + rconst2_2*(i-1)
556 : end do
557 0 : call random_seed(put=seed)
558 0 : call random_number (randu2)
559 0 : deallocate(seed)
560 :
561 : !Normal distribution (Mean = 0, standard dev = 1)
562 0 : randn_prescribed_aero = boxMuller(randu1,randu2)
563 0 : randn_persists = randn_prescribed_aero
564 : else
565 : !Use the previously generated random number
566 0 : randn_prescribed_aero = randn_persists
567 : endif
568 0 : end function randn_prescribed_aero
569 : !-------------------------------------------------------------------
570 0 : function logvm_get_index(name,type) result(index)
571 : implicit none
572 :
573 : !Args
574 : character(len=*), intent(in) :: name, type
575 :
576 : !Loc
577 : character(len=64) :: tmp_name
578 : character(len=1000) :: msgStr
579 : integer :: index, i
580 :
581 :
582 0 : index = -1
583 0 : tmp_name = trim(adjustl(name))//'_log'//trim(adjustl(type))
584 :
585 0 : fldloop: do i = 1, number_flds
586 0 : if(fields(i)%fldnam == tmp_name) then
587 : index = i
588 : exit fldloop
589 : endif
590 : enddo fldloop
591 :
592 0 : if(index == -1) then
593 0 : write(msgStr,*) "Prescribed_aero.F90: ",tmp_name," doesn't exist in the fields%fldnam"
594 0 : call endrun(msgStr)
595 : endif
596 :
597 0 : end function logvm_get_index
598 : !-------------------------------------------------------------------
599 0 : function boxMuller(ru1,ru2) result(rn)
600 : use physconst, only : pi
601 : implicit none
602 :
603 : !Args
604 : real(r8), intent(in) :: ru1, ru2
605 :
606 : !Loc
607 : real(r8), parameter :: pi2 = 2._r8 * pi
608 : real(r8) :: ur, theta, rn
609 :
610 : !Based on Box Muller Method, generate normally distributed random numbers
611 0 : ur = sqrt(-2._r8 * log(ru1))
612 0 : theta = pi2 * ru2
613 :
614 : !Normal distribution (Mean = 0, standard dev = 1)
615 0 : rn = ur * cos(theta)
616 :
617 0 : end function boxMuller
618 :
619 : !-------------------------------------------------------------------
620 : !-------------------------------------------------------------------
621 1536 : subroutine init_prescribed_aero_restart( piofile )
622 : use pio, only : file_desc_t, pio_def_var, pio_double
623 : use tracer_data, only : init_trc_restart
624 : implicit none
625 : type(file_desc_t),intent(inout) :: pioFile ! pio File pointer
626 : integer :: ierr
627 :
628 : ! For allowing randn_persists to persist during reststarts
629 1536 : ierr = pio_def_var(piofile, randn_persists_name, pio_double, randn_persists_desc)
630 :
631 1536 : call init_trc_restart( 'prescribed_aero', piofile, file )
632 :
633 1536 : end subroutine init_prescribed_aero_restart
634 : !-------------------------------------------------------------------
635 1536 : subroutine write_prescribed_aero_restart( piofile )
636 1536 : use tracer_data, only : write_trc_restart
637 : use pio, only : file_desc_t, pio_put_var
638 : implicit none
639 :
640 : type(file_desc_t) :: piofile
641 : integer :: ierr
642 :
643 : ! For allowing randn_persists to persist during reststarts
644 3072 : ierr = pio_put_var(piofile, randn_persists_desc, (/randn_persists/))
645 :
646 1536 : call write_trc_restart( piofile, file )
647 :
648 1536 : end subroutine write_prescribed_aero_restart
649 :
650 : !-------------------------------------------------------------------
651 : !-------------------------------------------------------------------
652 768 : subroutine read_prescribed_aero_restart( pioFile )
653 1536 : use tracer_data, only : read_trc_restart
654 : use pio, only : file_desc_t, pio_inq_varid, pio_get_var
655 : implicit none
656 :
657 : type(file_desc_t) :: piofile
658 : integer :: ierr
659 :
660 : ! For allowing randn_persists to persist during reststarts
661 768 : ierr = pio_inq_varid(pioFile, randn_persists_name, randn_persists_desc)
662 768 : ierr = pio_get_var(pioFile, randn_persists_desc, randn_persists)
663 :
664 768 : call read_trc_restart( 'prescribed_aero', piofile, file )
665 :
666 768 : end subroutine read_prescribed_aero_restart
667 :
668 : end module prescribed_aero
|
