Line data Source code
1 : module cam3_aero_data
2 : !-----------------------------------------------------------------------
3 : !
4 : ! Purposes:
5 : ! read, store, interpolate, and return fields
6 : ! of aerosols to CAM. The initialization
7 : ! file (mass.nc) is assumed to be a monthly climatology
8 : ! of aerosols from MATCH (on a sigma pressure
9 : ! coordinate system).
10 : ! also provide a "background" aerosol field to correct
11 : ! for any deficiencies in the physical parameterizations
12 : ! This fields is a "tuning" parameter.
13 : ! Public methods:
14 : ! (1) - initialization
15 : ! read aerosol masses from external file
16 : ! also pressure coordinates
17 : ! convert from monthly average values to mid-month values
18 : ! (2) - interpolation (time and vertical)
19 : ! interpolate onto pressure levels of CAM
20 : ! interpolate to time step of CAM
21 : ! return mass of aerosols
22 : !
23 : !-----------------------------------------------------------------------
24 :
25 : use shr_kind_mod, only: r8 => shr_kind_r8
26 : use shr_scam_mod, only: shr_scam_GetCloseLatLon
27 : use spmd_utils, only: masterproc
28 : use ppgrid, only: pcols, pver, pverp, begchunk, endchunk
29 : use phys_grid, only: get_ncols_p, scatter_field_to_chunk
30 : use time_manager, only: get_curr_calday
31 : use infnan, only: nan, assignment(=)
32 : use cam_abortutils, only: endrun
33 : use scamMod, only: scmlon,scmlat,single_column
34 : use error_messages, only: handle_ncerr
35 : use physics_types, only: physics_state
36 : use boundarydata, only: boundarydata_init, boundarydata_type
37 : use perf_mod, only: t_startf, t_stopf
38 : use cam_logfile, only: iulog
39 : use netcdf
40 :
41 : implicit none
42 : private
43 : save
44 :
45 : public :: &
46 : cam3_aero_data_readnl, & ! read namelist
47 : cam3_aero_data_register, & ! register these aerosols with pbuf2d
48 : cam3_aero_data_init, & ! read from file, interpolate onto horiz grid
49 : cam3_aero_data_timestep_init ! update data-aerosols to this timestep
50 :
51 : ! namelist variables
52 : logical, public :: cam3_aero_data_on = .false.
53 : character(len=256) :: bndtvaer = 'bndtvaer' ! full pathname for time-variant aerosol mass climatology dataset
54 :
55 : ! naer is number of species in climatology
56 : integer, parameter :: naer = 11
57 :
58 : real(r8), parameter :: wgt_sscm = 6.0_r8 / 7.0_r8 ! Fraction of total seasalt mass in coarse mode
59 :
60 : ! indices to aerosol array (species portion)
61 : integer, parameter :: &
62 : idxSUL = 1, &
63 : idxSSLTA = 2, & ! accumulation mode
64 : idxSSLTC = 3, & ! coarse mode
65 : idxOCPHO = 8, &
66 : idxBCPHO = 9, &
67 : idxOCPHI = 10, &
68 : idxBCPHI = 11
69 :
70 : ! indices to sections of array that represent
71 : ! groups of aerosols
72 : integer, parameter :: &
73 : idxSSLTfirst = 2, numSSLT = 2, &
74 : idxDUSTfirst = 4, &
75 : numDUST = 4, &
76 : idxCARBONfirst = 8, &
77 : numCARBON = 4
78 :
79 : ! names of aerosols are they are represented in
80 : ! the climatology file.
81 : ! Appended '_V' indicates field has been vertically summed.
82 : character(len=8), parameter :: aerosol_name(naer) = &
83 : (/"MSUL_V "&
84 : ,"MSSLTA_V"&
85 : ,"MSSLTC_V"&
86 : ,"MDUST1_V"&
87 : ,"MDUST2_V"&
88 : ,"MDUST3_V"&
89 : ,"MDUST4_V"&
90 : ,"MOCPHO_V"&
91 : ,"MBCPHO_V"&
92 : ,"MOCPHI_V"&
93 : ,"MBCPHI_V"/)
94 :
95 : ! number of different "groups" of aerosols
96 : integer, parameter :: num_aer_groups=4
97 :
98 : ! which group does each bin belong to?
99 : integer, dimension(naer), parameter :: &
100 : group =(/1,2,2,3,3,3,3,4,4,4,4/)
101 :
102 : ! name of each group
103 : character(len=10), dimension(num_aer_groups), parameter :: &
104 : aerosol_names = (/'sul ','sslt ','dust ','car '/)
105 :
106 : ! this boundarydata_type is used for datasets in the ncols format only.
107 : type(boundarydata_type) :: aerosol_datan
108 :
109 : integer :: aernid = -1 ! netcdf id for aerosol file (init to invalid)
110 : integer :: species_id(naer) = -1 ! netcdf_id of each aerosol species (init to invalid)
111 : integer :: Mpsid ! netcdf id for MATCH PS
112 : integer :: nm = 1 ! index to prv month in array. init to 1 and toggle between 1 and 2
113 : integer :: np = 2 ! index to nxt month in array. init to 2 and toggle between 1 and 2
114 : integer :: mo_nxt = huge(1) ! index to nxt month in file
115 :
116 : real(r8) :: cdaym ! calendar day of prv month
117 : real(r8) :: cdayp ! calendar day of next month
118 :
119 : ! aerosol mass
120 : real(r8), allocatable :: aer_mass(:, :, :, :)
121 :
122 : ! Days into year for mid month date
123 : ! This variable is dumb, the dates are in the dataset to be read in but they are
124 : ! slightly different than this so getting rid of it causes a change which
125 : ! exceeds roundoff.
126 : real(r8) :: Mid(12) = (/16.5_r8, 46.0_r8, 75.5_r8, 106.0_r8, 136.5_r8, 167.0_r8, &
127 : 197.5_r8, 228.5_r8, 259.0_r8, 289.5_r8, 320.0_r8, 350.5_r8 /)
128 :
129 : ! values read from file and temporary values used for interpolation
130 : !
131 : ! aerosolc is:
132 : ! Cumulative Mass at midpoint of each month
133 : ! on CAM's horizontal grid (col)
134 : ! on MATCH's levels (lev)
135 : ! aerosolc
136 : integer, parameter :: paerlev = 28 ! number of levels for aerosol fields (MUST = naerlev)
137 : integer :: naerlev ! size of level dimension in MATCH data
138 : integer :: naerlon
139 : integer :: naerlat
140 : real(r8), pointer :: M_hybi(:) ! MATCH hybi
141 : real(r8), pointer :: M_ps(:,:) ! surface pressure from MATCH file
142 : real(r8), pointer :: aerosolc(:,:,:,:,:) ! Aerosol cumulative mass from MATCH
143 : real(r8), pointer :: M_ps_cam_col(:,:,:) ! PS from MATCH on Cam Columns
144 :
145 : ! indices for fields in the physics buffer
146 : integer :: cam3_sul_idx, cam3_ssam_idx, cam3_sscm_idx, &
147 : cam3_dust1_idx, cam3_dust2_idx, cam3_dust3_idx, cam3_dust4_idx,&
148 : cam3_ocpho_idx, cam3_bcpho_idx, cam3_ocphi_idx, cam3_bcphi_idx
149 :
150 : !================================================================================================
151 : contains
152 : !================================================================================================
153 :
154 1536 : subroutine cam3_aero_data_readnl(nlfile)
155 :
156 : use namelist_utils, only: find_group_name
157 : use units, only: getunit, freeunit
158 : use mpishorthand
159 :
160 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
161 :
162 : ! Local variables
163 : integer :: unitn, ierr
164 : character(len=*), parameter :: subname = 'cam3_aero_data_readnl'
165 :
166 : namelist /cam3_aero_data_nl/ cam3_aero_data_on, bndtvaer
167 : !-----------------------------------------------------------------------------
168 :
169 1536 : if (masterproc) then
170 2 : unitn = getunit()
171 2 : open( unitn, file=trim(nlfile), status='old' )
172 2 : call find_group_name(unitn, 'cam3_aero_data_nl', status=ierr)
173 2 : if (ierr == 0) then
174 0 : read(unitn, cam3_aero_data_nl, iostat=ierr)
175 0 : if (ierr /= 0) then
176 0 : call endrun(subname // ':: ERROR reading namelist')
177 : end if
178 : end if
179 2 : close(unitn)
180 2 : call freeunit(unitn)
181 : end if
182 :
183 : #ifdef SPMD
184 : ! Broadcast namelist variables
185 1536 : call mpibcast(cam3_aero_data_on, 1, mpilog, 0, mpicom)
186 1536 : call mpibcast(bndtvaer, len(bndtvaer), mpichar, 0, mpicom)
187 : #endif
188 :
189 : ! Prevent using these before they are set.
190 1536 : cdaym = nan
191 1536 : cdayp = nan
192 :
193 1536 : end subroutine cam3_aero_data_readnl
194 :
195 : !================================================================================================
196 :
197 0 : subroutine cam3_aero_data_register
198 :
199 : ! register old prescribed aerosols with physics buffer
200 :
201 : use physics_buffer, only: pbuf_add_field, dtype_r8
202 :
203 0 : call pbuf_add_field('cam3_sul', 'physpkg',dtype_r8,(/pcols,pver/),cam3_sul_idx)
204 0 : call pbuf_add_field('cam3_ssam', 'physpkg',dtype_r8,(/pcols,pver/),cam3_ssam_idx)
205 0 : call pbuf_add_field('cam3_sscm', 'physpkg',dtype_r8,(/pcols,pver/),cam3_sscm_idx)
206 0 : call pbuf_add_field('cam3_dust1','physpkg',dtype_r8,(/pcols,pver/),cam3_dust1_idx)
207 0 : call pbuf_add_field('cam3_dust2','physpkg',dtype_r8,(/pcols,pver/),cam3_dust2_idx)
208 0 : call pbuf_add_field('cam3_dust3','physpkg',dtype_r8,(/pcols,pver/),cam3_dust3_idx)
209 0 : call pbuf_add_field('cam3_dust4','physpkg',dtype_r8,(/pcols,pver/),cam3_dust4_idx)
210 0 : call pbuf_add_field('cam3_ocpho','physpkg',dtype_r8,(/pcols,pver/),cam3_ocpho_idx)
211 0 : call pbuf_add_field('cam3_bcpho','physpkg',dtype_r8,(/pcols,pver/),cam3_bcpho_idx)
212 0 : call pbuf_add_field('cam3_ocphi','physpkg',dtype_r8,(/pcols,pver/),cam3_ocphi_idx)
213 0 : call pbuf_add_field('cam3_bcphi','physpkg',dtype_r8,(/pcols,pver/),cam3_bcphi_idx)
214 :
215 0 : end subroutine cam3_aero_data_register
216 :
217 : !================================================================================================
218 :
219 0 : subroutine cam3_aero_data_init(phys_state)
220 : !------------------------------------------------------------------
221 : ! Reads in:
222 : ! file from which to read aerosol Masses on CAM grid. Currently
223 : ! assumed to be MATCH ncep runs, averaged by month.
224 : ! NOTE (Data have been externally interpolated onto CAM grid
225 : ! and backsolved to provide Mid-month values)
226 : !
227 : ! Populates:
228 : ! module variables:
229 : ! aerosolc(pcols,paerlev+1,begchunk:endchunk,naer,2))
230 : ! aerosolc( column_index
231 : ! , level_index (match levels)
232 : ! , chunk_index
233 : ! , species_index
234 : ! , month = 1:2 )
235 : ! M_hybi(level_index = Lev_MATCH) = pressure at mid-level.
236 : ! M_ps_cam_col(column,chunk,month) ! PS from MATCH on Cam Columns
237 : !
238 : ! Method:
239 : ! read data from file
240 : ! allocate memory for storage of aerosol data on CAM horizontal grid
241 : ! distribute data to remote nodes
242 : ! populates the module variables
243 : !
244 : !------------------------------------------------------------------
245 0 : use ioFileMod, only: getfil
246 :
247 : #if ( defined SPMD )
248 : use mpishorthand
249 : #endif
250 : type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
251 :
252 : ! local variables
253 :
254 : integer :: naerlev
255 :
256 : integer dateid ! netcdf id for date variable
257 : integer secid ! netcdf id for seconds variable
258 : integer londimid ! netcdf id for longitude dimension
259 : integer latdimid ! netcdf id for latitude dimension
260 : integer levdimid ! netcdf id for level dimension
261 :
262 : integer timesiz ! number of time samples (=12) in netcdf file
263 : integer latid ! netcdf id for latitude variable
264 : integer Mhybiid ! netcdf id for MATCH hybi
265 : integer timeid ! netcdf id for time variable
266 : integer dimids(nf90_max_var_dims) ! variable shape
267 : integer :: start(4) ! start vector for netcdf calls
268 : integer :: kount(4) ! count vector for netcdf calls
269 : integer mo ! month index
270 : integer m ! constituent index
271 : integer :: n ! loop index
272 : integer :: i,j,k ! spatial indices
273 : integer :: date_aer(12) ! Date on aerosol dataset (YYYYMMDD)
274 : integer :: attnum ! attribute number
275 : integer :: ierr ! netcdf return code
276 : real(r8) :: coldata(paerlev) ! aerosol field read in from dataset
277 : integer :: ret
278 : integer mo_prv ! index to previous month
279 : integer latidx,lonidx
280 :
281 : character(len=8) :: aname ! temporary aerosol name
282 : character(len=8) :: tmp_aero_name(naer) ! name for input to boundary data
283 :
284 : character(len=256) :: locfn ! netcdf local filename to open
285 : !
286 : ! aerosol_data will be read in from the aerosol boundary dataset, then scattered to chunks
287 : ! after filling in the bottom level with zeros
288 : !
289 0 : real(r8), allocatable :: aerosol_data(:,:,:) ! aerosol field read in from dataset
290 0 : real(r8), allocatable :: aerosol_field(:,:,:) ! (plon,paerlev+1,plat) aerosol field to be scattered
291 : real(r8) :: caldayloc ! calendar day of current timestep
292 : real(r8) :: closelat,closelon
293 :
294 : character(len=*), parameter :: subname = 'cam3_aero_data_init'
295 : !------------------------------------------------------------------
296 :
297 0 : call t_startf(subname)
298 :
299 0 : allocate (aer_mass(pcols, pver, naer, begchunk:endchunk) )
300 :
301 : ! set new aerosol names because input file has 1 seasalt bin
302 0 : do m = 1, naer
303 0 : tmp_aero_name(m)=aerosol_name(m)
304 0 : if (aerosol_name(m)=='MSSLTA_V') tmp_aero_name(m) = 'MSSLT_V'
305 0 : if (aerosol_name(m)=='MSSLTC_V') tmp_aero_name(m) = 'MSSLT_V'
306 : end do
307 :
308 0 : allocate (aerosolc(pcols,paerlev+1,begchunk:endchunk,naer,2))
309 0 : aerosolc(:,:,:,:,:) = 0._r8
310 :
311 0 : caldayloc = get_curr_calday ()
312 :
313 0 : if (caldayloc < Mid(1)) then
314 0 : mo_prv = 12
315 0 : mo_nxt = 1
316 0 : else if (caldayloc >= Mid(12)) then
317 0 : mo_prv = 12
318 0 : mo_nxt = 1
319 : else
320 0 : do i = 2 , 12
321 0 : if (caldayloc < Mid(i)) then
322 0 : mo_prv = i-1
323 0 : mo_nxt = i
324 0 : exit
325 : end if
326 : end do
327 : end if
328 :
329 : ! Set initial calendar day values
330 0 : cdaym = Mid(mo_prv)
331 0 : cdayp = Mid(mo_nxt)
332 :
333 0 : if (masterproc) &
334 0 : write(iulog,*) subname//': CAM3 prescribed aerosol dataset is: ', trim(bndtvaer)
335 :
336 0 : call getfil (bndtvaer, locfn, 0)
337 :
338 : call handle_ncerr( nf90_open (locfn, 0, aernid),&
339 0 : subname, __LINE__)
340 :
341 0 : if (single_column) &
342 0 : call shr_scam_GetCloseLatLon(aernid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
343 :
344 : ! Check to see if this dataset is in ncol format.
345 0 : aerosol_datan%isncol=.false.
346 0 : ierr = nf90_inq_dimid( aernid, 'ncol', londimid )
347 0 : if ( ierr==NF90_NOERR ) then
348 :
349 0 : aerosol_datan%isncol=.true.
350 0 : call handle_ncerr(nf90_close(aernid),subname, __LINE__)
351 :
352 : call boundarydata_init(bndtvaer, phys_state, tmp_aero_name, naer, &
353 0 : aerosol_datan, 3)
354 :
355 0 : aerosolc(:,1:paerlev,:,:,:)=aerosol_datan%fields
356 :
357 0 : M_ps_cam_col=>aerosol_datan%ps
358 0 : M_hybi=>aerosol_datan%hybi
359 :
360 : else
361 :
362 : ! Allocate memory for dynamic arrays local to this module
363 0 : allocate (M_ps_cam_col(pcols,begchunk:endchunk,2))
364 0 : allocate (M_hybi(paerlev+1))
365 : ! TBH: HACK to avoid use of uninitialized values when ncols < pcols
366 0 : M_ps_cam_col(:,:,:) = 0._r8
367 :
368 0 : if (masterproc) then
369 :
370 : ! First ensure dataset is CAM-ready
371 :
372 : call handle_ncerr(nf90_inquire_attribute (aernid, nf90_global, 'cam-ready', attnum=attnum),&
373 0 : subname//': interpaerosols needs to be run to create a cam-ready aerosol dataset')
374 :
375 : ! Get and check dimension info
376 :
377 : call handle_ncerr( nf90_inq_dimid( aernid, 'lon', londimid ),&
378 0 : subname, __LINE__)
379 : call handle_ncerr( nf90_inq_dimid( aernid, 'lev', levdimid ),&
380 0 : subname, __LINE__)
381 : call handle_ncerr( nf90_inq_dimid( aernid, 'time', timeid ),&
382 0 : subname, __LINE__)
383 : call handle_ncerr( nf90_inq_dimid( aernid, 'lat', latdimid ),&
384 0 : subname, __LINE__)
385 : call handle_ncerr( nf90_inquire_dimension( aernid, londimid, len=naerlon ),&
386 0 : subname, __LINE__)
387 : call handle_ncerr( nf90_inquire_dimension( aernid, levdimid, len=naerlev ),&
388 0 : subname, __LINE__)
389 : call handle_ncerr( nf90_inquire_dimension( aernid, latdimid, len=naerlat ),&
390 0 : subname, __LINE__)
391 : call handle_ncerr( nf90_inquire_dimension( aernid, timeid, len=timesiz ),&
392 0 : subname, __LINE__)
393 :
394 : call handle_ncerr( nf90_inq_varid( aernid, 'date', dateid ),&
395 0 : subname, __LINE__)
396 : call handle_ncerr( nf90_inq_varid( aernid, 'datesec', secid ),&
397 0 : subname, __LINE__)
398 :
399 0 : do m = 1, naer
400 0 : aname=aerosol_name(m)
401 : ! rename because file has only one seasalt field
402 0 : if (aname=='MSSLTA_V') aname = 'MSSLT_V'
403 0 : if (aname=='MSSLTC_V') aname = 'MSSLT_V'
404 : call handle_ncerr( nf90_inq_varid( aernid, TRIM(aname), species_id(m)), &
405 0 : subname, __LINE__)
406 : end do
407 :
408 : call handle_ncerr( nf90_inq_varid( aernid, 'lat', latid ),&
409 0 : subname, __LINE__)
410 :
411 : ! quick sanity check on one field
412 : call handle_ncerr( nf90_inquire_variable (aernid, species_id(1), dimids=dimids),&
413 0 : subname, __LINE__)
414 :
415 : if ( (dimids(4) /= timeid) .or. &
416 : (dimids(3) /= levdimid) .or. &
417 0 : (dimids(2) /= latdimid) .or. &
418 : (dimids(1) /= londimid) ) then
419 0 : write(iulog,*) subname//': Data must be ordered time, lev, lat, lon'
420 0 : write(iulog,*) 'data are ordered as', dimids(4), dimids(3), dimids(2), dimids(1)
421 0 : write(iulog,*) 'data should be ordered as', timeid, levdimid, latdimid, londimid
422 0 : call endrun ()
423 : end if
424 :
425 : ! use hybi,PS from MATCH
426 : call handle_ncerr( nf90_inq_varid( aernid, 'hybi', Mhybiid ),&
427 0 : subname, __LINE__)
428 : call handle_ncerr( nf90_inq_varid( aernid, 'PS', Mpsid ),&
429 0 : subname, __LINE__)
430 :
431 : ! check dimension order for MATCH's surface pressure
432 : call handle_ncerr( nf90_inquire_variable (aernid, Mpsid, dimids=dimids),&
433 0 : subname, __LINE__)
434 : if ( (dimids(3) /= timeid) .or. &
435 0 : (dimids(2) /= latdimid) .or. &
436 : (dimids(1) /= londimid) ) then
437 0 : write(iulog,*) subname//': Pressure must be ordered time, lat, lon'
438 0 : write(iulog,*) 'data are ordered as', dimids(3), dimids(2), dimids(1)
439 0 : write(iulog,*) 'data should be ordered as', timeid, levdimid, latdimid, londimid
440 0 : call endrun ()
441 : end if
442 :
443 : ! read in hybi from MATCH
444 : call handle_ncerr( nf90_get_var (aernid, Mhybiid, M_hybi),&
445 0 : subname, __LINE__)
446 :
447 : ! Retrieve date and sec variables.
448 : call handle_ncerr( nf90_get_var (aernid, dateid, date_aer),&
449 0 : subname, __LINE__)
450 0 : if (timesiz < 12) then
451 0 : write(iulog,*) subname//': When cycling aerosols, dataset must have 12 consecutive ', &
452 0 : 'months of data starting with Jan'
453 0 : write(iulog,*) 'Current dataset has only ',timesiz,' months'
454 0 : call endrun ()
455 : end if
456 0 : do mo = 1,12
457 0 : if (mod(date_aer(mo),10000)/100 /= mo) then
458 0 : write(iulog,*) subname//': When cycling aerosols, dataset must have 12 consecutive ', &
459 0 : 'months of data starting with Jan'
460 0 : write(iulog,*)'Month ',mo,' of dataset says date=',date_aer(mo)
461 0 : call endrun ()
462 : end if
463 : end do
464 0 : if (single_column) then
465 0 : naerlat=1
466 0 : naerlon=1
467 : endif
468 0 : kount(:) = (/naerlon,naerlat,paerlev,1/)
469 : end if ! masterproc
470 :
471 : ! broadcast hybi to nodes
472 :
473 : #if ( defined SPMD )
474 0 : call mpibcast (M_hybi, paerlev+1, mpir8, 0, mpicom)
475 0 : call mpibcast (kount, 3, mpiint, 0, mpicom)
476 0 : naerlon = kount(1)
477 0 : naerlat = kount(2)
478 : #endif
479 0 : allocate(aerosol_field(kount(1),kount(3)+1,kount(2)))
480 0 : allocate(M_ps(kount(1),kount(2)))
481 0 : if (masterproc) allocate(aerosol_data(kount(1),kount(2),kount(3)))
482 :
483 : ! Retrieve Aerosol Masses (kg/m^2 in each layer), transpose to model order (lon,lev,lat),
484 : ! then scatter to slaves.
485 0 : if (nm /= 1 .or. np /= 2) call endrun (subname//': bad nm or np value')
486 0 : do n=nm,np
487 0 : if (n == 1) then
488 0 : mo = mo_prv
489 : else
490 0 : mo = mo_nxt
491 : end if
492 :
493 0 : do m=1,naer
494 0 : if (masterproc) then
495 0 : if (single_column) then
496 0 : start(:) = (/lonidx,latidx,1,mo/)
497 : else
498 0 : start(:) = (/1,1,1,mo/)
499 : endif
500 0 : kount(:) = (/naerlon,naerlat,paerlev,1/)
501 :
502 0 : call handle_ncerr( nf90_get_var (aernid, species_id(m),aerosol_data, start, kount),&
503 0 : subname, __LINE__)
504 0 : do j=1,naerlat
505 0 : do k=1,paerlev
506 0 : aerosol_field(:,k,j) = aerosol_data(:,j,k)
507 : end do
508 0 : aerosol_field(:,paerlev+1,j) = 0._r8 ! value at bottom
509 : end do
510 :
511 : end if
512 : call scatter_field_to_chunk (1, paerlev+1, 1, naerlon, aerosol_field, &
513 0 : aerosolc(:,:,:,m,n))
514 : end do
515 :
516 : ! Retrieve PS from Match
517 :
518 0 : if (masterproc) then
519 0 : if (single_column) then
520 0 : start(:) = (/lonidx,latidx,mo,-1/)
521 : else
522 0 : start(:) = (/1,1,mo,-1/)
523 : endif
524 0 : kount(:) = (/naerlon,naerlat,1,-1/)
525 : call handle_ncerr( nf90_get_var(aernid, Mpsid, M_ps,start,kount),&
526 0 : subname, __LINE__)
527 : end if
528 0 : call scatter_field_to_chunk (1, 1, 1, naerlon, M_ps(:,:), M_ps_cam_col(:,:,n))
529 : end do ! n=nm,np (=1,2)
530 :
531 0 : if(masterproc) deallocate(aerosol_data)
532 0 : deallocate(aerosol_field)
533 :
534 : end if ! Check to see if this dataset is in ncol format.
535 :
536 0 : call t_stopf(subname)
537 :
538 0 : end subroutine cam3_aero_data_init
539 :
540 : !================================================================================================
541 :
542 0 : subroutine cam3_aero_data_timestep_init(pbuf2d, phys_state)
543 : !------------------------------------------------------------------
544 : !
545 : ! Input:
546 : ! time at which aerosol masses are needed (get_curr_calday())
547 : ! chunk index
548 : ! CAM's vertical grid (pint)
549 : !
550 : ! Output:
551 : ! values for Aerosol Mass at time specified by get_curr_calday
552 : ! on vertical grid specified by pint (aer_mass) :: aerosol at time t
553 : !
554 : ! Method:
555 : ! first determine which indexs of aerosols are the bounding data sets
556 : ! interpolate both onto vertical grid aerm(),aerp().
557 : ! from those two, interpolate in time.
558 : !
559 : !------------------------------------------------------------------
560 :
561 : use interpolate_data, only: get_timeinterp_factors
562 :
563 : use physics_buffer, only: physics_buffer_desc, dtype_r8, pbuf_set_field, pbuf_get_chunk
564 : use cam_logfile, only: iulog
565 : use ppgrid, only: begchunk,endchunk
566 : use physconst, only: gravit
567 :
568 : !
569 : ! aerosol fields interpolated to current time step
570 : ! on pressure levels of this time step.
571 : ! these should be made read-only for other modules
572 : ! Is allocation done correctly here?
573 : !
574 :
575 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
576 : type(physics_state), intent(in), dimension(begchunk:endchunk) :: phys_state
577 :
578 : !
579 : ! Local workspace
580 : !
581 0 : type(physics_buffer_desc), pointer :: phys_buffer_chunk(:)
582 : real(r8) :: pint(pcols,pverp) ! interface pres.
583 : integer :: c ! chunk index
584 : real(r8) caldayloc ! calendar day of current timestep
585 : real(r8) fact1, fact2 ! time interpolation factors
586 :
587 : integer i, k, j ! spatial indices
588 : integer m ! constituent index
589 : integer lats(pcols),lons(pcols) ! latitude and longitudes of column
590 : integer ncol ! number of columns
591 : integer lchnk ! chunk index
592 :
593 : real(r8) speciesmin(naer) ! minimal value for each species
594 : !
595 : ! values before current time step "the minus month"
596 : ! aerosolm(pcols,pver) is value of preceeding month's aerosol masses
597 : ! aerosolp(pcols,pver) is value of next month's aerosol masses
598 : ! (think minus and plus or values to left and right of point to be interpolated)
599 : !
600 0 : real(r8) aerosolm(pcols,pver,naer,begchunk:endchunk) ! aerosol mass from MATCH in column,level at previous (minus) month
601 : !
602 : ! values beyond (or at) current time step "the plus month"
603 : !
604 0 : real(r8) aerosolp(pcols,pver,naer,begchunk:endchunk) ! aerosol mass from MATCH in column,level at next (plus) month
605 : real(r8) :: mass_to_mmr(pcols,pver)
606 :
607 : character(len=*), parameter :: subname = 'cam3_aero_data_timestep_init'
608 :
609 : logical error_found
610 : !------------------------------------------------------------------
611 :
612 0 : call aerint(phys_state)
613 :
614 0 : caldayloc = get_curr_calday ()
615 :
616 : ! Determine time interpolation factors. 1st arg says we are cycling 1 year of data
617 : call get_timeinterp_factors (.true., mo_nxt, cdaym, cdayp, caldayloc, &
618 0 : fact1, fact2, 'GET_AEROSOL:')
619 :
620 : ! interpolate (prv and nxt month) bounding datasets onto cam vertical grid.
621 : ! compute mass mixing ratios on CAMS's pressure coordinate
622 : ! for both the "minus" and "plus" months
623 : !
624 : ! This loop over chunk could probably be removed by working with the whole
625 : ! begchunk:endchunk group at once. It would require a slight generalization
626 : ! in vert_interpolate.
627 0 : do c = begchunk,endchunk
628 :
629 0 : lchnk = phys_state(c)%lchnk
630 0 : pint = phys_state(c)%pint
631 0 : ncol = get_ncols_p(c)
632 :
633 0 : call vert_interpolate (M_ps_cam_col(:,c,nm), pint, nm, aerosolm(:,:,:,c), ncol, c)
634 0 : call vert_interpolate (M_ps_cam_col(:,c,np), pint, np, aerosolp(:,:,:,c), ncol, c)
635 :
636 : ! Time interpolate.
637 0 : do m=1,naer
638 0 : do k=1,pver
639 0 : do i=1,ncol
640 0 : aer_mass(i,k,m,c) = aerosolm(i,k,m,c)*fact1 + aerosolp(i,k,m,c)*fact2
641 : end do
642 : end do
643 : ! Partition seasalt aerosol mass
644 0 : if (m .eq. idxSSLTA) then
645 0 : aer_mass(:ncol,:,m,c) = (1._r8-wgt_sscm)*aer_mass(:ncol,:,m,c) ! fraction of seasalt mass in accumulation mode
646 0 : elseif (m .eq. idxSSLTC) then
647 0 : aer_mass(:ncol,:,m,c) = wgt_sscm*aer_mass(:ncol,:,m,c) ! fraction of seasalt mass in coarse mode
648 : endif
649 : end do
650 :
651 : ! exit if mass is negative (we have previously set
652 : ! cumulative mass to be a decreasing function.)
653 0 : speciesmin(:) = 0._r8 ! speciesmin(m) = 0 is minimum mass for each species
654 :
655 0 : error_found = .false.
656 0 : do m=1,naer
657 0 : do k=1,pver
658 0 : do i=1,ncol
659 0 : if (aer_mass(i, k, m,c) < speciesmin(m)) error_found = .true.
660 : end do
661 : end do
662 : end do
663 0 : if (error_found) then
664 0 : do m=1,naer
665 0 : do k=1,pver
666 0 : do i=1,ncol
667 0 : if (aer_mass(i, k, m,c) < speciesmin(m)) then
668 0 : write(iulog,*) subname//': negative mass mixing ratio, exiting'
669 0 : write(iulog,*) 'm, column, pver',m, i, k ,aer_mass(i, k, m,c)
670 0 : call endrun ()
671 : end if
672 : end do
673 : end do
674 : end do
675 : end if
676 0 : do k = 1, pver
677 0 : mass_to_mmr(1:ncol,k) = gravit/(pint(1:ncol,k+1)-pint(1:ncol,k))
678 : enddo
679 :
680 0 : phys_buffer_chunk => pbuf_get_chunk(pbuf2d, lchnk)
681 :
682 0 : call pbuf_set_field(phys_buffer_chunk, cam3_sul_idx, aer_mass(1:ncol,:, idxSUL,c)*mass_to_mmr(:ncol,:), &
683 0 : start=(/1,1/), kount=(/ncol,pver/))
684 0 : call pbuf_set_field(phys_buffer_chunk, cam3_ssam_idx, aer_mass(1:ncol,:, idxSSLTA,c)*mass_to_mmr(:ncol,:), &
685 0 : start=(/1,1/), kount=(/ncol,pver/))
686 0 : call pbuf_set_field(phys_buffer_chunk, cam3_sscm_idx, aer_mass(1:ncol,:, idxSSLTC,c)*mass_to_mmr(:ncol,:), &
687 0 : start=(/1,1/), kount=(/ncol,pver/))
688 0 : call pbuf_set_field(phys_buffer_chunk, cam3_dust1_idx, aer_mass(1:ncol,:, idxDUSTfirst,c)*mass_to_mmr(:ncol,:), &
689 0 : start=(/1,1/), kount=(/ncol,pver/))
690 0 : call pbuf_set_field(phys_buffer_chunk, cam3_dust2_idx, aer_mass(1:ncol,:,idxDUSTfirst+1,c)*mass_to_mmr(:ncol,:), &
691 0 : start=(/1,1/), kount=(/ncol,pver/))
692 0 : call pbuf_set_field(phys_buffer_chunk, cam3_dust3_idx, aer_mass(1:ncol,:,idxDUSTfirst+2,c)*mass_to_mmr(:ncol,:), &
693 0 : start=(/1,1/), kount=(/ncol,pver/))
694 0 : call pbuf_set_field(phys_buffer_chunk, cam3_dust4_idx, aer_mass(1:ncol,:,idxDUSTfirst+3,c)*mass_to_mmr(:ncol,:), &
695 0 : start=(/1,1/), kount=(/ncol,pver/))
696 0 : call pbuf_set_field(phys_buffer_chunk, cam3_ocpho_idx, aer_mass(1:ncol,:, idxOCPHO,c)*mass_to_mmr(:ncol,:), &
697 0 : start=(/1,1/), kount=(/ncol,pver/))
698 0 : call pbuf_set_field(phys_buffer_chunk, cam3_bcpho_idx, aer_mass(1:ncol,:, idxBCPHO,c)*mass_to_mmr(:ncol,:), &
699 0 : start=(/1,1/), kount=(/ncol,pver/))
700 0 : call pbuf_set_field(phys_buffer_chunk, cam3_ocphi_idx, aer_mass(1:ncol,:, idxOCPHI,c)*mass_to_mmr(:ncol,:), &
701 0 : start=(/1,1/), kount=(/ncol,pver/))
702 0 : call pbuf_set_field(phys_buffer_chunk, cam3_bcphi_idx, aer_mass(1:ncol,:, idxBCPHI,c)*mass_to_mmr(:ncol,:), &
703 0 : start=(/1,1/), kount=(/ncol,pver/))
704 :
705 : enddo ! c = begchunk:endchunk
706 :
707 0 : end subroutine cam3_aero_data_timestep_init
708 :
709 : !================================================================================================
710 :
711 0 : subroutine vert_interpolate (Match_ps, pint, n, aerosol_mass, ncol, c)
712 : !--------------------------------------------------------------------
713 : ! Input: match surface pressure, cam interface pressure,
714 : ! month index, number of columns, chunk index
715 : !
716 : ! Output: Aerosol mass mixing ratio (aerosol_mass)
717 : !
718 : ! Method:
719 : ! interpolate column mass (cumulative) from match onto
720 : ! cam's vertical grid (pressure coordinate)
721 : ! convert back to mass mixing ratio
722 : !
723 : !--------------------------------------------------------------------
724 :
725 : real(r8), intent(out) :: aerosol_mass(pcols,pver,naer) ! aerosol mass from MATCH
726 : real(r8), intent(in) :: Match_ps(pcols) ! surface pressure at a particular month
727 : real(r8), intent(in) :: pint(pcols,pverp) ! interface pressure from CAM
728 :
729 : integer, intent(in) :: ncol,c ! chunk index and number of columns
730 : integer, intent(in) :: n ! prv or nxt month index
731 : !
732 : ! Local workspace
733 : !
734 : integer m ! index to aerosol species
735 : integer kupper(pcols) ! last upper bound for interpolation
736 : integer i, k, kk, kkstart, kount ! loop vars for interpolation
737 : integer isv, ksv, msv ! loop indices to save
738 :
739 : logical bad ! indicates a bad point found
740 : logical lev_interp_comp ! interpolation completed for a level
741 : logical error_found
742 :
743 : real(r8) aerosol(pcols,pverp,naer) ! cumulative mass of aerosol in column beneath upper
744 : ! interface of level in column at particular month
745 : real(r8) dpl, dpu ! lower and upper intepolation factors
746 : real(r8) v_coord ! vertical coordinate
747 : real(r8) AER_diff ! temp var for difference between aerosol masses
748 :
749 : character(len=*), parameter :: subname = 'cam3_aero_data.vert_interpolate'
750 : !-----------------------------------------------------------------------
751 :
752 0 : call t_startf ('vert_interpolate')
753 : !
754 : ! Initialize index array
755 : !
756 0 : do i=1,ncol
757 0 : kupper(i) = 1
758 : end do
759 : !
760 : ! assign total mass to topmost level
761 : !
762 0 : aerosol(:,1,:) = aerosolc(:,1,c,:,n)
763 : !
764 : ! At every pressure level, interpolate onto that pressure level
765 : !
766 0 : do k=2,pver
767 : !
768 : ! Top level we need to start looking is the top level for the previous k
769 : ! for all longitude points
770 : !
771 0 : kkstart = paerlev+1
772 0 : do i=1,ncol
773 0 : kkstart = min0(kkstart,kupper(i))
774 : end do
775 : kount = 0
776 : !
777 : ! Store level indices for interpolation
778 : !
779 : ! for the pressure interpolation should be comparing
780 : ! pint(column,lev) with M_hybi(lev)*M_ps_cam_col(month,column,chunk)
781 : !
782 : lev_interp_comp = .false.
783 0 : do kk=kkstart,paerlev
784 0 : if(.not.lev_interp_comp) then
785 0 : do i=1,ncol
786 0 : v_coord = pint(i,k)
787 0 : if (M_hybi(kk)*Match_ps(i) .lt. v_coord .and. v_coord .le. M_hybi(kk+1)*Match_ps(i)) then
788 0 : kupper(i) = kk
789 0 : kount = kount + 1
790 : end if
791 : end do
792 : !
793 : ! If all indices for this level have been found, do the interpolation and
794 : ! go to the next level
795 : !
796 : ! Interpolate in pressure.
797 : !
798 0 : if (kount.eq.ncol) then
799 0 : do m=1,naer
800 0 : do i=1,ncol
801 0 : dpu = pint(i,k) - M_hybi(kupper(i))*Match_ps(i)
802 0 : dpl = M_hybi(kupper(i)+1)*Match_ps(i) - pint(i,k)
803 0 : aerosol(i,k,m) = &
804 0 : (aerosolc(i,kupper(i) ,c,m,n)*dpl + &
805 0 : aerosolc(i,kupper(i)+1,c,m,n)*dpu)/(dpl + dpu)
806 : enddo !i
807 : end do
808 : lev_interp_comp = .true.
809 : end if
810 : end if
811 : end do
812 : !
813 : ! If we've fallen through the kk=1,levsiz-1 loop, we cannot interpolate and
814 : ! must extrapolate from the bottom or top pressure level for at least some
815 : ! of the longitude points.
816 : !
817 :
818 0 : if(.not.lev_interp_comp) then
819 0 : do m=1,naer
820 0 : do i=1,ncol
821 0 : if (pint(i,k) .lt. M_hybi(1)*Match_ps(i)) then
822 0 : aerosol(i,k,m) = aerosolc(i,1,c,m,n)
823 0 : else if (pint(i,k) .gt. M_hybi(paerlev+1)*Match_ps(i)) then
824 0 : aerosol(i,k,m) = 0.0_r8
825 : else
826 0 : dpu = pint(i,k) - M_hybi(kupper(i))*Match_ps(i)
827 0 : dpl = M_hybi(kupper(i)+1)*Match_ps(i) - pint(i,k)
828 0 : aerosol(i,k,m) = &
829 0 : (aerosolc(i,kupper(i) ,c,m,n)*dpl + &
830 0 : aerosolc(i,kupper(i)+1,c,m,n)*dpu)/(dpl + dpu)
831 : end if
832 : end do
833 : end do
834 :
835 0 : if (kount.gt.ncol) then
836 0 : call endrun (subname//': Bad data: non-monotonicity suspected in dependent variable')
837 : end if
838 : end if
839 : end do
840 :
841 : ! call t_startf ('vi_checks')
842 : !
843 : ! aerosol mass beneath lowest interface (pverp) must be 0
844 : !
845 0 : aerosol(1:ncol,pverp,:) = 0._r8
846 : !
847 : ! Set mass in layer to zero whenever it is less than
848 : ! 1.e-40 kg/m^2 in the layer
849 : !
850 0 : do m = 1, naer
851 0 : do k = 1, pver
852 0 : do i = 1, ncol
853 0 : if (aerosol(i,k,m) < 1.e-40_r8) aerosol(i,k,m) = 0._r8
854 : end do
855 : end do
856 : end do
857 : !
858 : ! Set mass in layer to zero whenever it is less than
859 : ! 10^-15 relative to column total mass
860 : !
861 0 : error_found = .false.
862 0 : do m = 1, naer
863 0 : do k = 1, pver
864 0 : do i = 1, ncol
865 0 : AER_diff = aerosol(i,k,m) - aerosol(i,k+1,m)
866 0 : if( abs(AER_diff) < 1e-15_r8*aerosol(i,1,m)) then
867 0 : AER_diff = 0._r8
868 : end if
869 0 : aerosol_mass(i,k,m)= AER_diff
870 0 : if (aerosol_mass(i,k,m) < 0) error_found = .true.
871 : end do
872 : end do
873 : end do
874 0 : if (error_found) then
875 0 : do m = 1, naer
876 0 : do k = 1, pver
877 0 : do i = 1, ncol
878 0 : if (aerosol_mass(i,k,m) < 0) then
879 0 : write(iulog,*) subname//': mass < 0, m, col, lev, mass',m, i, k, aerosol_mass(i,k,m)
880 0 : write(iulog,*) subname//': aerosol(k),(k+1)',aerosol(i,k,m),aerosol(i,k+1,m)
881 0 : write(iulog,*) subname//': pint(k+1),(k)',pint(i,k+1),pint(i,k)
882 0 : write(iulog,*)'n,c',n,c
883 0 : call endrun()
884 : end if
885 : end do
886 : end do
887 : end do
888 : end if
889 :
890 0 : call t_stopf ('vert_interpolate')
891 :
892 0 : return
893 0 : end subroutine vert_interpolate
894 :
895 : !================================================================================================
896 :
897 0 : subroutine aerint (phys_state)
898 :
899 : type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
900 :
901 : integer :: ntmp ! used in index swapping
902 : integer :: start(4) ! start vector for netcdf calls
903 : integer :: kount(4) ! count vector for netcdf calls
904 : integer :: i,j,k ! spatial indices
905 : integer :: m ! constituent index
906 : integer :: cols, cole
907 : integer :: lchnk, ncol
908 : real(r8) :: caldayloc ! calendar day of current timestep
909 0 : real(r8) :: aerosol_data(naerlon,naerlat,paerlev) ! aerosol field read in from dataset
910 0 : real(r8) :: aerosol_field(naerlon,paerlev+1,naerlat) ! aerosol field to be scattered
911 : integer latidx,lonidx
912 : real(r8) closelat,closelon
913 :
914 : character(len=*), parameter :: subname = 'cam3_aero_data.aerint'
915 : !-----------------------------------------------------------------------
916 :
917 0 : if (single_column) &
918 0 : call shr_scam_GetCloseLatLon(aernid,scmlat,scmlon,closelat,closelon,latidx,lonidx)
919 :
920 : !
921 : ! determine if need to read in next month data
922 : ! also determine time interpolation factors
923 : !
924 0 : caldayloc = get_curr_calday ()
925 : !
926 : ! If model time is past current forward timeslice, then
927 : ! masterproc reads in the next timeslice for time interpolation. Messy logic is
928 : ! for interpolation between December and January (mo_nxt == 1). Just like
929 : ! ozone_data_timestep_init, sstint.
930 : !
931 0 : if (caldayloc > cdayp .and. .not. (mo_nxt == 1 .and. caldayloc >= cdaym)) then
932 0 : mo_nxt = mod(mo_nxt,12) + 1
933 0 : cdaym = cdayp
934 0 : cdayp = Mid(mo_nxt)
935 : !
936 : ! Check for valid date info
937 : !
938 0 : if (.not. (mo_nxt == 1 .or. caldayloc <= cdayp)) then
939 0 : call endrun (subname//': Non-monotonicity suspected in input aerosol data')
940 : end if
941 :
942 0 : ntmp = nm
943 0 : nm = np
944 0 : np = ntmp
945 :
946 0 : if(aerosol_datan%isncol) then
947 0 : do lchnk=begchunk,endchunk
948 0 : ncol=phys_state(lchnk)%ncol
949 0 : cols=1
950 0 : cole=cols+aerosol_datan%count(cols,lchnk)-1
951 0 : do while(cole<=ncol)
952 0 : start=(/aerosol_datan%start(cols,lchnk),mo_nxt,1,-1/)
953 0 : kount=(/aerosol_datan%count(cols,lchnk),1,-1,-1/)
954 : call handle_ncerr( nf90_get_var(aerosol_datan%ncid, aerosol_datan%psid , &
955 0 : aerosol_datan%ps(cols:cole,lchnk,np), start(1:2), &
956 : kount(1:2)),&
957 0 : subname, __LINE__)
958 0 : start(2)=1
959 0 : start(3)=mo_nxt
960 0 : kount(2)=paerlev
961 0 : kount(3)=1
962 0 : do m=1,naer
963 0 : call handle_ncerr( nf90_get_var(aerosol_datan%ncid, aerosol_datan%dataid(m) , &
964 0 : aerosol_datan%fields(cols:cole,:,lchnk,m,np), &
965 : start(1:3), kount(1:3)),&
966 0 : subname, __LINE__)
967 :
968 : end do
969 0 : if(cols==ncol) exit
970 0 : cols=cols+aerosol_datan%count(cols,lchnk)
971 0 : cole=cols+aerosol_datan%count(cols,lchnk)-1
972 : end do
973 : end do
974 0 : aerosolc(:,1:paerlev,:,:,np)=aerosol_datan%fields(:,:,:,:,np)
975 : else
976 0 : do m=1,naer
977 0 : if (masterproc) then
978 0 : if (single_column) then
979 0 : naerlon=1
980 0 : naerlat=1
981 0 : start(:) = (/lonidx,latidx,1,mo_nxt/)
982 : else
983 0 : start(:) = (/1,1,1,mo_nxt/)
984 : endif
985 0 : kount(:) = (/naerlon,naerlat,paerlev,1/)
986 0 : call handle_ncerr( nf90_get_var (aernid, species_id(m), aerosol_data, start, kount),&
987 0 : subname, __LINE__)
988 :
989 0 : do j=1,naerlat
990 0 : do k=1,paerlev
991 0 : aerosol_field(:,k,j) = aerosol_data(:,j,k)
992 : end do
993 0 : aerosol_field(:,paerlev+1,j) = 0._r8 ! value at bottom
994 : end do
995 : end if
996 : call scatter_field_to_chunk (1, paerlev+1, 1, naerlon, aerosol_field, &
997 0 : aerosolc(:,:,:,m,np))
998 : end do
999 : !
1000 : ! Retrieve PS from Match
1001 : !
1002 0 : if (masterproc) then
1003 0 : if (single_column) then
1004 0 : naerlon=1
1005 0 : naerlat=1
1006 0 : start(:) = (/lonidx,latidx,mo_nxt,-1/)
1007 : else
1008 0 : start(:) = (/1,1,mo_nxt,-1/)
1009 : endif
1010 0 : kount(:) = (/naerlon,naerlat,1,-1/)
1011 : call handle_ncerr( nf90_get_var (aernid, Mpsid, M_ps, start, kount),&
1012 0 : subname, __LINE__)
1013 0 : write(iulog,*) subname//': Read aerosols data for julian day', Mid(mo_nxt)
1014 : end if
1015 0 : call scatter_field_to_chunk (1, 1, 1, naerlon, M_ps(:,:), M_ps_cam_col(:,:,np))
1016 : end if
1017 : end if
1018 :
1019 0 : end subroutine aerint
1020 :
1021 : end module cam3_aero_data
|
