Line data Source code
1 : module mo_fstrat
2 : !---------------------------------------------------------------
3 : ! ... variables for the upper boundary values
4 : !---------------------------------------------------------------
5 :
6 : use shr_kind_mod, only : r8 => shr_kind_r8
7 : use ppgrid, only : pcols, pver,pverp, begchunk, endchunk
8 : use chem_mods, only : gas_pcnst
9 : use cam_abortutils, only : endrun
10 : use cam_pio_utils, only : cam_pio_openfile
11 : use pio
12 : use cam_logfile, only : iulog
13 :
14 : implicit none
15 :
16 : private
17 : public :: fstrat_inti
18 : public :: set_fstrat_vals
19 : public :: set_fstrat_h2o
20 : public :: has_fstrat
21 :
22 : save
23 :
24 : real(r8), parameter :: taurelax = 864000._r8 ! 10 days
25 : integer :: table_nox_ndx = -1
26 : integer :: table_h2o_ndx = -1
27 : integer :: table_ox_ndx = -1
28 : integer :: no_ndx
29 : integer :: no2_ndx
30 : integer :: h2o_ndx
31 : integer :: ox_ndx
32 : integer :: o3s_ndx
33 : integer :: o3inert_ndx
34 : integer :: o3a_ndx
35 : integer :: xno_ndx
36 : integer :: xno2_ndx
37 : integer :: o3rad_ndx
38 : real(r8) :: facrelax
39 : real(r8) :: days(12)
40 : real(r8), allocatable :: ub_plevs(:) ! table midpoint pressure levels (Pa)
41 : real(r8), allocatable :: ub_plevse(:) ! table edge pressure levels (Pa)
42 : integer :: ub_nlevs ! # of levs in ubc file
43 : integer :: ub_nlat ! # of lats in ubc file
44 : integer :: ub_nspecies ! # of species in ubc file
45 : integer :: ub_nmonth ! # of months in ubc file
46 : real(r8), allocatable :: mr_ub(:,:,:,:,:) ! vmr
47 : integer, allocatable :: map(:) ! species indices for ubc species
48 : logical :: sim_has_nox
49 : integer :: dtime ! model time step (s)
50 : logical :: has_fstrat(gas_pcnst)
51 :
52 : contains
53 :
54 0 : subroutine fstrat_inti( fstrat_file, fstrat_list )
55 : !------------------------------------------------------------------
56 : ! ... initialize upper boundary values
57 : !------------------------------------------------------------------
58 :
59 : use mo_constants, only : d2r
60 : use phys_grid, only : get_ncols_p, get_rlat_all_p
61 :
62 : use time_manager, only : get_step_size
63 : use time_manager, only : get_calday
64 : use ioFileMod, only : getfil
65 : use spmd_utils, only : masterproc
66 : use mo_tracname, only : solsym
67 : use chem_mods, only : gas_pcnst
68 : use mo_chem_utls, only : get_spc_ndx, get_inv_ndx
69 : use constituents, only : pcnst
70 : use interpolate_data, only : interp_type, lininterp_init, lininterp
71 :
72 : character(len=*), intent(in) :: fstrat_file
73 : character(len=*), intent(in) :: fstrat_list(:)
74 :
75 : !------------------------------------------------------------------
76 : ! ... local variables
77 : !------------------------------------------------------------------
78 : real(r8), parameter :: mb2pa = 100._r8
79 :
80 : integer :: i, j, nchar
81 : integer :: spcno, lev, month, ierr
82 : integer :: vid, ndims
83 : type(file_desc_t) :: ncid
84 : integer :: dimid_lat, dimid_lev, dimid_species, dimid_month
85 : integer :: dimid(4)
86 : integer :: start(4)
87 : integer :: count(4)
88 : integer, parameter :: dates(12) = (/ 116, 214, 316, 415, 516, 615, &
89 : 716, 816, 915, 1016, 1115, 1216 /)
90 : integer :: ncols, c
91 0 : real(r8), allocatable :: mr_ub_in(:,:,:,:)
92 0 : real(r8), allocatable :: lat(:)
93 : real(r8) :: to_lats(pcols)
94 : character(len=80) :: attribute
95 : character(len=8) :: wrk_name
96 0 : character(len=25), allocatable :: ub_species_names(:)
97 : character(len=256) :: locfn
98 : type(interp_type) :: lat_wgts
99 :
100 :
101 : !-----------------------------------------------------------------------
102 : ! ... get species indicies
103 : !-----------------------------------------------------------------------
104 0 : no_ndx = get_spc_ndx( 'NO' )
105 0 : no2_ndx = get_spc_ndx( 'NO2' )
106 0 : sim_has_nox = no_ndx > 0 .or. no2_ndx > 0
107 0 : ox_ndx = get_spc_ndx( 'OX' )
108 0 : if( ox_ndx < 1 ) then
109 0 : ox_ndx = get_spc_ndx( 'O3' )
110 : end if
111 0 : o3s_ndx = get_spc_ndx( 'O3S' )
112 0 : o3inert_ndx = get_spc_ndx( 'O3INERT' )
113 0 : o3rad_ndx = get_spc_ndx( 'O3RAD' )
114 0 : o3a_ndx = get_spc_ndx( 'O3A' )
115 0 : xno_ndx = get_spc_ndx( 'XNO' )
116 0 : xno2_ndx = get_spc_ndx( 'XNO2' )
117 :
118 0 : if(masterproc) then
119 0 : if (ox_ndx > 0) then
120 0 : if ( .not. any(fstrat_list == 'O3') ) then
121 0 : write(iulog,*) 'fstrat_inti: ***WARNING*** O3 is not include in fstrat_list namelist variable'
122 : endif
123 : endif
124 : end if
125 0 : h2o_ndx = get_spc_ndx( 'H2O' )
126 0 : if( h2o_ndx < 0 ) then
127 0 : h2o_ndx = get_inv_ndx( 'H2O' )
128 : end if
129 :
130 0 : has_fstrat(:) = .false.
131 :
132 0 : do i = 1,pcnst
133 :
134 0 : if ( len_trim(fstrat_list(i))==0 ) exit
135 :
136 0 : j = get_spc_ndx(fstrat_list(i))
137 :
138 0 : if ( j > 0 ) then
139 0 : has_fstrat(j) = .true.
140 : else
141 0 : write(iulog,*) 'fstrat_inti: '//trim(fstrat_list(i))//' is not included in species set'
142 0 : call endrun('fstrat_inti: invalid fixed stratosphere species')
143 : endif
144 :
145 : enddo
146 :
147 : if (.not. any(has_fstrat(:))) return
148 :
149 : !-----------------------------------------------------------------------
150 : ! ... open netcdf file
151 : !-----------------------------------------------------------------------
152 : call getfil (fstrat_file, locfn, 0)
153 : call cam_pio_openfile (ncid, trim(locfn), PIO_NOWRITE)
154 : !-----------------------------------------------------------------------
155 : ! ... get latitude
156 : !-----------------------------------------------------------------------
157 : ierr = pio_inq_dimid( ncid, 'lat', dimid_lat )
158 : ierr = pio_inq_dimlen( ncid, dimid_lat, ub_nlat )
159 : allocate( lat(ub_nlat), stat=ierr )
160 : if( ierr /= 0 ) then
161 : write(iulog,*) 'fstrat_inti: lat allocation error = ',ierr
162 : call endrun
163 : end if
164 : ierr = pio_inq_varid( ncid, 'lat', vid )
165 : ierr = pio_get_var( ncid, vid, lat )
166 : lat(:ub_nlat) = lat(:ub_nlat) * d2r
167 :
168 : if( ierr /= 0 ) then
169 : write(iulog,*) 'fstrat_inti: failed to deallocate lat; ierr = ',ierr
170 : call endrun
171 : end if
172 :
173 : !-----------------------------------------------------------------------
174 : ! ... get vertical coordinate (if necessary, convert units to pa)
175 : !-----------------------------------------------------------------------
176 : ierr = pio_inq_dimid( ncid, 'lev', dimid_lev )
177 : ierr = pio_inq_dimlen( ncid, dimid_lev, ub_nlevs )
178 : allocate( ub_plevs(ub_nlevs), stat=ierr )
179 : if( ierr /= 0 ) then
180 : write(iulog,*) 'fstrat_inti: ub_plevs allocation error = ',ierr
181 : call endrun
182 : end if
183 : ierr = pio_inq_varid( ncid, 'lev', vid )
184 : ierr = pio_get_var( ncid, vid, ub_plevs )
185 : attribute(:) = ' '
186 : call pio_seterrorhandling(ncid, pio_BCAST_error)
187 : ierr = pio_get_att( ncid, vid, 'units', attribute )
188 : call pio_seterrorhandling(ncid, pio_INTERNAL_error)
189 : if( ierr == PIO_noerr )then
190 : if( trim(attribute) == 'mb' .or. trim(attribute) == 'hpa' )then
191 : if (masterproc) write(iulog,*) 'fstrat_inti: units for lev = ',trim(attribute),'... converting to pa'
192 : ub_plevs(:) = mb2pa * ub_plevs(:)
193 : else if( trim(attribute) /= 'pa' .and. trim(attribute) /= 'pa' )then
194 : write(iulog,*) 'fstrat_inti: unknown units for lev, units=*',trim(attribute),'*'
195 : write(iulog,*) 'fstrat_inti: ',attribute=='mb',trim(attribute)=='mb',attribute(1:2)=='mb'
196 : call endrun
197 : end if
198 : else
199 : write(iulog,*) 'fstrat_inti: warning! units attribute for lev missing, assuming mb'
200 : ub_plevs(:) = mb2pa * ub_plevs(:)
201 : end if
202 : !-----------------------------------------------------------------------
203 : ! ... get time and species dimensions
204 : !-----------------------------------------------------------------------
205 : ierr = pio_inq_dimid( ncid, 'month', dimid_month )
206 : ierr = pio_inq_dimlen( ncid, dimid_month, ub_nmonth )
207 : if( ub_nmonth /= 12 )then
208 : write(iulog,*) 'fstrat_inti: error! number of months = ',ub_nmonth,', expecting 12'
209 : call endrun
210 : end if
211 : ierr = pio_inq_dimid( ncid, 'species', dimid_species )
212 : ierr = pio_inq_dimlen( ncid, dimid_species, ub_nspecies )
213 :
214 : !------------------------------------------------------------------
215 : ! ... allocate arrays
216 : !------------------------------------------------------------------
217 : allocate( mr_ub_in(ub_nlat,ub_nspecies,ub_nmonth,ub_nlevs), stat=ierr )
218 : if( ierr /= 0 ) then
219 : write(iulog,*) 'fstrat_inti: mr_ub_in allocation error = ',ierr
220 : call endrun
221 : end if
222 : allocate( mr_ub(pcols,ub_nspecies,ub_nmonth,ub_nlevs,begchunk:endchunk), stat=ierr )
223 : if( ierr /= 0 ) then
224 : write(iulog,*) 'fstrat_inti: mr_ub allocation error = ',ierr
225 : call endrun
226 : end if
227 :
228 : !------------------------------------------------------------------
229 : ! ... read in the species names
230 : !------------------------------------------------------------------
231 :
232 : ierr = pio_inq_varid( ncid, 'specname', vid )
233 : ierr = pio_inq_vardimid( ncid, vid, dimid )
234 : ierr = pio_inq_dimlen( ncid, dimid(1), nchar )
235 : allocate( ub_species_names(ub_nspecies), stat=ierr )
236 : if( ierr /= 0 ) then
237 : write(iulog,*) 'fstrat_inti: ub_species_names allocation error = ',ierr
238 : call endrun
239 : end if
240 : allocate( map(ub_nspecies), stat=ierr )
241 : if( ierr /= 0 ) then
242 : write(iulog,*) 'fstrat_inti: map allocation error = ',ierr
243 : call endrun
244 : end if
245 :
246 : table_loop : do i = 1,ub_nspecies
247 : start(:2) = (/ 1, i /)
248 : count(:2) = (/ nchar, 1 /)
249 : ub_species_names(i)(:) = ' '
250 : ierr = pio_get_var( ncid, vid, start(:2), count(:2), ub_species_names(i:i))
251 : if( trim(ub_species_names(i)) == 'NOX' ) then
252 : table_nox_ndx = i
253 : else if( trim(ub_species_names(i)) == 'H2O' ) then
254 : table_h2o_ndx = i
255 : else if( trim(ub_species_names(i)) == 'OX' ) then
256 : table_ox_ndx = i
257 : end if
258 : map(i) = 0
259 : do j = 1,gas_pcnst
260 : if( trim(ub_species_names(i)) == trim(solsym(j)) ) then
261 : if( has_fstrat(j) ) then
262 : map(i) = j
263 : if( masterproc ) write(iulog,*) 'fstrat_inti: '//trim(solsym(j))//' is fixed in stratosphere'
264 : exit
265 : end if
266 : endif
267 : end do
268 : if( map(i) == 0 ) then
269 : if( trim(ub_species_names(i)) == 'OX' ) then
270 : if( o3rad_ndx > 0 ) then
271 : wrk_name = 'O3RAD'
272 : else
273 : wrk_name = 'O3'
274 : end if
275 : do j = 1,gas_pcnst
276 : if( trim(wrk_name) == trim(solsym(j)) ) then
277 : if( has_fstrat(j) ) then
278 : if( masterproc ) write(iulog,*) 'fstrat_inti: '//trim(solsym(j))//' is fixed in stratosphere'
279 : map(i) = j
280 : exit
281 : end if
282 : end if
283 : end do
284 : if( map(i) == 0 ) then
285 : write(iulog,*) 'fstrat_inti: ubc table species ',trim(ub_species_names(i)), ' not used'
286 : end if
287 : else if( (trim(ub_species_names(i)) /= 'NOX') .and. (trim(ub_species_names(i)) /= 'H2O') ) then
288 : write(iulog,*) 'fstrat_inti: ubc table species ',trim(ub_species_names(i)), ' not used'
289 : end if
290 : end if
291 : end do table_loop
292 :
293 : if( table_nox_ndx > 0 ) then
294 : if ( any(fstrat_list(:) == 'NO') .or. any(fstrat_list(:) == 'NO2') ) then
295 : map(table_nox_ndx) = gas_pcnst + 1
296 : else
297 : write(iulog,*) 'fstrat_inti: ubc table species ',trim(ub_species_names(table_nox_ndx)), ' not used'
298 : endif
299 : end if
300 : if( table_h2o_ndx > 0 ) then
301 : if ( h2o_ndx > 0 ) then
302 : map(table_h2o_ndx) = gas_pcnst + 2
303 : else
304 : write(iulog,*) 'fstrat_inti: ubc table species ',trim(ub_species_names(table_h2o_ndx)), ' not used'
305 : endif
306 : end if
307 :
308 : if (masterproc) write(iulog,*) 'fstrat_inti: h2o_ndx, table_h2o_ndx = ', h2o_ndx, table_h2o_ndx
309 :
310 : !------------------------------------------------------------------
311 : ! ... check dimensions for vmr variable
312 : !------------------------------------------------------------------
313 : ierr = pio_inq_varid( ncid, 'vmr', vid )
314 : ierr = pio_inq_varndims( ncid, vid, ndims )
315 : if( ndims /= 4 ) then
316 : write(iulog,*) 'fstrat_inti: error! variable vmr has ndims = ',ndims,', expecting 4'
317 : call endrun
318 : end if
319 : ierr = pio_inq_vardimid( ncid, vid, dimid )
320 : if( dimid(1) /= dimid_lat .or. dimid(2) /= dimid_species .or. &
321 : dimid(3) /= dimid_month .or. dimid(4) /= dimid_lev )then
322 : write(iulog,*) 'fstrat_inti: error! dimensions in wrong order for variable vmr,'// &
323 : 'expecting (lat,species,month,lev)'
324 : call endrun
325 : end if
326 :
327 : !------------------------------------------------------------------
328 : ! ... read in the ub mixing ratio values
329 : !------------------------------------------------------------------
330 : start = (/ 1, 1, 1, 1 /)
331 : count = (/ ub_nlat, ub_nspecies, ub_nmonth, ub_nlevs /)
332 :
333 : ierr = pio_get_var( ncid, vid, start, count, mr_ub_in )
334 : call pio_closefile (ncid)
335 : !--------------------------------------------------------------------
336 : ! ... regrid
337 : !--------------------------------------------------------------------
338 : do c=begchunk,endchunk
339 : ncols = get_ncols_p(c)
340 : call get_rlat_all_p(c, pcols, to_lats)
341 : call lininterp_init(lat,ub_nlat, to_lats, ncols, 1, lat_wgts)
342 :
343 : do lev = 1,ub_nlevs
344 : do month = 1,ub_nmonth
345 : do spcno = 1,ub_nspecies
346 : if( map(spcno) > 0 ) then
347 :
348 : call lininterp(mr_ub_in(:,spcno, month, lev), ub_nlat, mr_ub(:,spcno,month,lev,c), &
349 : ncols, lat_wgts)
350 : #ifdef DEBUG
351 : if( lev == 25 .and. month == 1 .and. spcno == 1 ) then
352 : write(iulog,*) 'mr_ub_in='
353 : write(iulog,'(10f7.1)') mr_ub_in(:,spcno,month,lev)*1.e9_r8
354 : write(iulog,*) 'mr_ub='
355 : write(iulog,'(10f7.1)') mr_ub(:,spcno,month,lev,c)*1.e9_r8
356 : end if
357 : #endif
358 : end if
359 :
360 : end do
361 : end do
362 : end do
363 : end do
364 : !--------------------------------------------------------
365 : ! ... initialize the monthly day of year times
366 : !--------------------------------------------------------
367 : do month = 1,12
368 : days(month) = get_calday( dates(month), 0 )
369 : end do
370 :
371 : !--------------------------------------------------------
372 : ! ... set up the relaxation for lower stratosphere
373 : !--------------------------------------------------------
374 : ! ... taurelax = relaxation timescale (in sec)
375 : ! facrelax = fractional relaxation towards ubc
376 : ! 1 => use ubc
377 : ! 0 => ignore ubc, use model concentrations
378 : !--------------------------------------------------------
379 : dtime = get_step_size()
380 : facrelax = 1._r8 - exp( -real(dtime)/taurelax )
381 :
382 : !--------------------------------------------------------
383 : ! ... setup conserving interp for OX
384 : !--------------------------------------------------------
385 : if( table_ox_ndx > 0 ) then
386 : allocate( ub_plevse(ub_nlevs-1), stat=ierr )
387 : if( ierr /= 0 ) then
388 : write(iulog,*) 'fstrat_inti: ub_plevse allocation error = ',ierr
389 : call endrun
390 : end if
391 : ub_plevse(1:ub_nlevs-1) = .5_r8*(ub_plevs(1:ub_nlevs-1) + ub_plevs(2:ub_nlevs))
392 : end if
393 :
394 0 : end subroutine fstrat_inti
395 :
396 0 : subroutine set_fstrat_vals( vmr, pmid, pint, ltrop, calday, ncol,lchnk )
397 :
398 : !--------------------------------------------------------------------
399 : ! ... set the upper boundary values for :
400 : ! ox, nox, hno3, ch4, co, n2o, n2o5 & stratospheric o3
401 : !--------------------------------------------------------------------
402 :
403 : !--------------------------------------------------------------------
404 : ! ... dummy args
405 : !--------------------------------------------------------------------
406 : integer, intent(in) :: lchnk ! chunk number
407 : integer, intent(in) :: ncol ! columns in chunk
408 : integer, intent(in) :: ltrop(pcols) ! tropopause vertical index
409 : real(r8), intent(in) :: calday ! day of year including fraction
410 : real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressure (Pa)
411 : real(r8), intent(in) :: pint(pcols,pverp) ! interface pressure (Pa)
412 : real(r8), intent(inout) :: vmr(ncol,pver,gas_pcnst) ! species concentrations (mol/mol)
413 :
414 : !--------------------------------------------------------------------
415 : ! ... local variables
416 : !--------------------------------------------------------------------
417 : integer, parameter :: zlower = pver
418 :
419 : integer :: m, last, next, i, k, k1, km
420 : integer :: astat
421 0 : integer :: kmax(ncol)
422 : integer :: levrelax
423 0 : integer :: kl(ncol,zlower)
424 0 : integer :: ku(ncol,zlower)
425 : real(r8) :: vmrrelax
426 : real(r8) :: pinterp
427 : real(r8) :: nox_ubc, xno, xno2, rno
428 : real(r8) :: dels
429 0 : real(r8) :: delp(ncol,zlower)
430 : real(r8) :: pint_vals(2)
431 0 : real(r8), allocatable :: table_ox(:)
432 :
433 : if (.not. any(has_fstrat(:))) return
434 :
435 : !--------------------------------------------------------
436 : ! ... setup the time interpolation
437 : !--------------------------------------------------------
438 : if( calday < days(1) ) then
439 : next = 1
440 : last = 12
441 : dels = (365._r8 + calday - days(12)) / (365._r8 + days(1) - days(12))
442 : else if( calday >= days(12) ) then
443 : next = 1
444 : last = 12
445 : dels = (calday - days(12)) / (365._r8 + days(1) - days(12))
446 : else
447 : do m = 11,1,-1
448 : if( calday >= days(m) ) then
449 : exit
450 : end if
451 : end do
452 : last = m
453 : next = m + 1
454 : dels = (calday - days(m)) / (days(m+1) - days(m))
455 : end if
456 : dels = max( min( 1._r8,dels ),0._r8 )
457 :
458 : !--------------------------------------------------------
459 : ! ... setup the pressure interpolation
460 : !--------------------------------------------------------
461 : do k = 1,zlower
462 : do i = 1,ncol
463 : if( pmid(i,k) <= ub_plevs(1) ) then
464 : kl(i,k) = 1
465 : ku(i,k) = 1
466 : delp(i,k) = 0._r8
467 : else if( pmid(i,k) >= ub_plevs(ub_nlevs) ) then
468 : kl(i,k) = ub_nlevs
469 : ku(i,k) = ub_nlevs
470 : delp(i,k) = 0._r8
471 : else
472 : pinterp = pmid(i,k)
473 : do k1 = 2,ub_nlevs
474 : if( pinterp <= ub_plevs(k1) ) then
475 : ku(i,k) = k1
476 : kl(i,k) = k1 - 1
477 : delp(i,k) = log( pinterp/ub_plevs(kl(i,k)) ) &
478 : / log( ub_plevs(ku(i,k))/ub_plevs(kl(i,k)) )
479 : exit
480 : end if
481 : end do
482 : end if
483 : end do
484 : end do
485 :
486 : !--------------------------------------------------------
487 : ! ... find max level less than 50 mb
488 : ! fix UB vals from top of model to this level
489 : !--------------------------------------------------------
490 : do i = 1,ncol
491 : do k = 2,pver
492 : if( pmid(i,k) > 50.e2_r8 ) then
493 : kmax(i) = k
494 : exit
495 : end if
496 : end do
497 : end do
498 :
499 : !--------------------------------------------------------
500 : ! ... setup for ox conserving interp
501 : !--------------------------------------------------------
502 : if( table_ox_ndx > 0 ) then
503 : if( map(table_ox_ndx) > 0 ) then
504 : allocate( table_ox(ub_nlevs-2),stat=astat )
505 : if( astat /= 0 ) then
506 : write(iulog,*) 'set_fstrat_vals: table_ox allocation error = ',astat
507 : call endrun
508 : end if
509 : #ifdef UB_DEBUG
510 : write(iulog,*) ' '
511 : write(iulog,*) 'set_fstrat_vals: ub_nlevs = ',ub_nlevs
512 : write(iulog,*) 'set_fstrat_vals: ub_plevse'
513 : write(iulog,'(1p5g15.7)') ub_plevse(:)
514 : write(iulog,*) ' '
515 : #endif
516 : end if
517 : end if
518 :
519 : !--------------------------------------------------------
520 : ! ... set the mixing ratios at upper boundary
521 : !--------------------------------------------------------
522 : species_loop : do m = 1,ub_nspecies
523 : ub_overwrite : if( map(m) > 0 ) then
524 : if( m == table_ox_ndx ) then
525 : do i = 1,ncol
526 :
527 : table_ox(1:ub_nlevs-2) = mr_ub(i,m,last,2:ub_nlevs-1,lchnk) &
528 : + dels*(mr_ub(i,m,next,2:ub_nlevs-1,lchnk) &
529 : - mr_ub(i,m,last,2:ub_nlevs-1,lchnk))
530 : #ifdef UB_DEBUG
531 : write(iulog,*) 'set_fstrat_vals: table_ox @ i = ',i
532 : write(iulog,'(1p5g15.7)') table_ox(:)
533 : write(iulog,*) ' '
534 : #endif
535 :
536 : km = kmax(i)
537 : #ifdef UB_DEBUG
538 : write(iulog,*) 'set_fstrat_vals: pint with km = ',km
539 : write(iulog,'(1p5g15.7)') pint(i,:km+1)
540 : write(iulog,*) ' '
541 : write(iulog,*) 'set_fstrat_vals: pmid with km = ',km
542 : write(iulog,'(1p5g15.7)') pmid(i,:km)
543 : write(iulog,*) ' '
544 : #endif
545 : call rebin( ub_nlevs-2, km, ub_plevse, pint(i,:km+1), table_ox, vmr(i,:km,map(m)) )
546 : #ifdef UB_DEBUG
547 : write(iulog,*) 'set_fstrat_vals: ub o3 @ i = ',i
548 : write(iulog,'(1p5g15.7)') vmr(i,:km,map(m))
549 : #endif
550 : end do
551 : deallocate( table_ox )
552 : cycle species_loop
553 : end if
554 : do i = 1,ncol
555 : do k = 1,kmax(i)
556 : pint_vals(1) = mr_ub(i,m,last,kl(i,k),lchnk) &
557 : + delp(i,k) &
558 : * (mr_ub(i,m,last,ku(i,k),lchnk) &
559 : - mr_ub(i,m,last,kl(i,k),lchnk))
560 : pint_vals(2) = mr_ub(i,m,next,kl(i,k),lchnk) &
561 : + delp(i,k) &
562 : * (mr_ub(i,m,next,ku(i,k),lchnk) &
563 : - mr_ub(i,m,next,kl(i,k),lchnk))
564 : if( m /= table_nox_ndx .and. m /= table_h2o_ndx ) then
565 : vmr(i,k,map(m)) = pint_vals(1) &
566 : + dels * (pint_vals(2) - pint_vals(1))
567 : else if( m == table_nox_ndx .and. sim_has_nox ) then
568 :
569 : nox_ubc = pint_vals(1) + dels * (pint_vals(2) - pint_vals(1))
570 : if( no_ndx > 0 ) then
571 : xno = vmr(i,k,no_ndx)
572 : else
573 : xno = 0._r8
574 : end if
575 : if( no2_ndx > 0 ) then
576 : xno2 = vmr(i,k,no2_ndx)
577 : else
578 : xno2 = 0._r8
579 : end if
580 : rno = xno / (xno + xno2)
581 : if( no_ndx > 0 ) then
582 : vmr(i,k,no_ndx) = rno * nox_ubc
583 : end if
584 : if( no2_ndx > 0 ) then
585 : vmr(i,k,no2_ndx) = (1._r8 - rno) * nox_ubc
586 : end if
587 : end if
588 : end do
589 : end do
590 : end if ub_overwrite
591 : end do species_loop
592 :
593 : col_loop2 : do i = 1,ncol
594 : !--------------------------------------------------------
595 : ! ... relax lower stratosphere to extended ubc
596 : ! check to make sure ubc is not being imposed too low
597 : ! levrelax = lowest model level (highest pressure)
598 : ! in which to relax to ubc
599 : !--------------------------------------------------------
600 : levrelax = ltrop(i)
601 : do while( pmid(i,levrelax) > ub_plevs(ub_nlevs) )
602 : levrelax = levrelax - 1
603 : end do
604 : #ifdef DEBUG
605 : if( levrelax /= ltrop(i) ) then
606 : write(iulog,*) 'warning -- raised ubc: ',i, &
607 : ltrop(i)-1,nint(pmid(i,ltrop(i)-1)/100._r8),'mb -->', &
608 : levrelax,nint(pmid(i,levrelax)/100._r8),'mb'
609 : end if
610 : #endif
611 : level_loop2 : do k = kmax(i)+1,levrelax
612 : if( sim_has_nox ) then
613 : if( no_ndx > 0 ) then
614 : xno = vmr(i,k,no_ndx)
615 : else
616 : xno = 0._r8
617 : end if
618 : if( no2_ndx > 0 ) then
619 : xno2 = vmr(i,k,no2_ndx)
620 : else
621 : xno2 = 0._r8
622 : end if
623 : rno = xno / (xno + xno2)
624 : nox_ubc = xno + xno2
625 : end if
626 : do m = 1,ub_nspecies
627 : if( map(m) < 1 ) then
628 : cycle
629 : end if
630 : pint_vals(1) = mr_ub(i,m,last,kl(i,k),lchnk) &
631 : + delp(i,k) &
632 : * (mr_ub(i,m,last,ku(i,k),lchnk) &
633 : - mr_ub(i,m,last,kl(i,k),lchnk))
634 : pint_vals(2) = mr_ub(i,m,next,kl(i,k),lchnk) &
635 : + delp(i,k) &
636 : * (mr_ub(i,m,next,ku(i,k),lchnk) &
637 : - mr_ub(i,m,next,kl(i,k),lchnk))
638 : vmrrelax = pint_vals(1) &
639 : + dels * (pint_vals(2) - pint_vals(1))
640 : if( m /= table_nox_ndx .and. m /= table_h2o_ndx ) then
641 : vmr(i,k,map(m)) = vmr(i,k,map(m)) &
642 : + (vmrrelax - vmr(i,k,map(m))) * facrelax
643 : else if( m == table_nox_ndx .and. sim_has_nox) then
644 :
645 : nox_ubc = nox_ubc + (vmrrelax - nox_ubc) * facrelax
646 : end if
647 : end do
648 : if( sim_has_nox ) then
649 : if( no_ndx > 0 ) then
650 : vmr(i,k,no_ndx) = rno * nox_ubc
651 : end if
652 : if( no2_ndx > 0 ) then
653 : vmr(i,k,no2_ndx) = (1._r8 - rno) * nox_ubc
654 : end if
655 : end if
656 : end do level_loop2
657 :
658 : !--------------------------------------------------------
659 : ! ... set O3S and O3INERT to OX when no synoz
660 : !--------------------------------------------------------
661 : if( ox_ndx > 0 ) then
662 : if( o3s_ndx > 0 ) then
663 : vmr(i,:ltrop(i),o3s_ndx) = vmr(i,:ltrop(i),ox_ndx)
664 : end if
665 : if( o3inert_ndx > 0 ) then
666 : vmr(i,:ltrop(i),o3inert_ndx) = vmr(i,:ltrop(i),ox_ndx)
667 : end if
668 : end if
669 :
670 :
671 : if ( o3a_ndx > 0 ) then
672 : vmr(i,:ltrop(i),o3a_ndx) = (1._r8 - facrelax ) * vmr(i,:ltrop(i),o3a_ndx)
673 : endif
674 : if ( xno_ndx > 0 ) then
675 : vmr(i,:ltrop(i),xno_ndx) = (1._r8 - facrelax ) * vmr(i,:ltrop(i),xno_ndx)
676 : endif
677 : if ( xno2_ndx > 0 ) then
678 : vmr(i,:ltrop(i),xno2_ndx) = (1._r8 - facrelax ) * vmr(i,:ltrop(i),xno2_ndx)
679 : endif
680 :
681 : end do col_loop2
682 :
683 0 : end subroutine set_fstrat_vals
684 :
685 0 : subroutine set_fstrat_h2o( h2o, pmid, ltrop, calday, ncol, lchnk )
686 : !--------------------------------------------------------------------
687 : ! ... set the h2o upper boundary values
688 : !--------------------------------------------------------------------
689 :
690 : !--------------------------------------------------------------------
691 : ! ... dummy args
692 : !--------------------------------------------------------------------
693 : integer, intent(in) :: ncol ! columns in chunk
694 : integer, intent(in) :: ltrop(pcols) ! tropopause vertical index
695 : integer, intent(in) :: lchnk
696 : real(r8), intent(in) :: calday ! day of year including fraction
697 : real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressure (Pa)
698 : real(r8), intent(inout) :: h2o(ncol,pver) ! h2o concentration (mol/mol)
699 :
700 : !--------------------------------------------------------------------
701 : ! ... local variables
702 : !--------------------------------------------------------------------
703 : integer, parameter :: zlower = pver
704 :
705 : integer :: m, last, next, i, k, k1
706 0 : integer :: kmax(ncol)
707 : integer :: levrelax
708 0 : integer :: kl(ncol,zlower)
709 0 : integer :: ku(ncol,zlower)
710 : real(r8) :: vmrrelax
711 : real(r8) :: pinterp
712 : real(r8) :: dels
713 0 : real(r8) :: delp(ncol,zlower)
714 : real(r8) :: pint_vals(2)
715 :
716 0 : h2o_overwrite : if( h2o_ndx > 0 .and. table_h2o_ndx > 0 ) then
717 : !--------------------------------------------------------
718 : ! ... setup the time interpolation
719 : !--------------------------------------------------------
720 0 : if( calday < days(1) ) then
721 0 : next = 1
722 0 : last = 12
723 0 : dels = (365._r8 + calday - days(12)) / (365._r8 + days(1) - days(12))
724 0 : else if( calday >= days(12) ) then
725 0 : next = 1
726 0 : last = 12
727 0 : dels = (calday - days(12)) / (365._r8 + days(1) - days(12))
728 : else
729 0 : do m = 11,1,-1
730 0 : if( calday >= days(m) ) then
731 : exit
732 : end if
733 : end do
734 0 : last = m
735 0 : next = m + 1
736 0 : dels = (calday - days(m)) / (days(m+1) - days(m))
737 : end if
738 0 : dels = max( min( 1._r8,dels ),0._r8 )
739 :
740 : !--------------------------------------------------------
741 : ! ... setup the pressure interpolation
742 : !--------------------------------------------------------
743 0 : do k = 1,zlower
744 0 : do i = 1,ncol
745 0 : if( pmid(i,k) <= ub_plevs(1) ) then
746 0 : kl(i,k) = 1
747 0 : ku(i,k) = 1
748 0 : delp(i,k) = 0._r8
749 0 : else if( pmid(i,k) >= ub_plevs(ub_nlevs) ) then
750 0 : kl(i,k) = ub_nlevs
751 0 : ku(i,k) = ub_nlevs
752 0 : delp(i,k) = 0._r8
753 : else
754 : pinterp = pmid(i,k)
755 0 : do k1 = 2,ub_nlevs
756 0 : if( pinterp <= ub_plevs(k1) ) then
757 0 : ku(i,k) = k1
758 0 : kl(i,k) = k1 - 1
759 0 : delp(i,k) = log( pinterp/ub_plevs(kl(i,k)) ) &
760 0 : / log( ub_plevs(ku(i,k))/ub_plevs(kl(i,k)) )
761 0 : exit
762 : end if
763 : end do
764 : end if
765 : end do
766 : end do
767 :
768 : !--------------------------------------------------------
769 : ! ... Find max level less than 50 mb
770 : ! fix UB vals from top of model to this level
771 : !--------------------------------------------------------
772 0 : do i = 1,ncol
773 0 : do k = 2,pver
774 0 : if( pmid(i,k) > 50.e2_r8 ) then
775 0 : kmax(i) = k
776 0 : exit
777 : end if
778 : end do
779 : end do
780 : !--------------------------------------------------------
781 : ! ... set the mixing ratio at upper boundary
782 : !--------------------------------------------------------
783 0 : m = table_h2o_ndx
784 0 : do i = 1,ncol
785 0 : do k = 1,kmax(i)
786 0 : pint_vals(1) = mr_ub(i,m,last,kl(i,k),lchnk) &
787 : + delp(i,k) &
788 0 : * (mr_ub(i,m,last,ku(i,k),lchnk) &
789 0 : - mr_ub(i,m,last,kl(i,k),lchnk))
790 0 : pint_vals(2) = mr_ub(i,m,next,kl(i,k),lchnk) &
791 : + delp(i,k) &
792 : * (mr_ub(i,m,next,ku(i,k),lchnk) &
793 0 : - mr_ub(i,m,next,kl(i,k),lchnk))
794 : h2o(i,k) = pint_vals(1) &
795 0 : + dels * (pint_vals(2) - pint_vals(1))
796 : end do
797 : end do
798 :
799 0 : col_loop2 : do i = 1,ncol
800 : !--------------------------------------------------------
801 : ! ... relax lower stratosphere to extended ubc
802 : ! check to make sure ubc is not being imposed too low
803 : ! levrelax = lowest model level (highest pressure)
804 : ! in which to relax to ubc
805 : !--------------------------------------------------------
806 0 : levrelax = ltrop(i)
807 0 : do while( pmid(i,levrelax) > ub_plevs(ub_nlevs) )
808 0 : levrelax = levrelax - 1
809 : end do
810 : #ifdef DEBUG
811 0 : if( levrelax /= ltrop(i) ) then
812 0 : write(iulog,*) 'warning -- raised ubc: ',i, &
813 0 : ltrop(i)-1,nint(pmid(i,ltrop(i)-1)/100._r8),'mb -->', &
814 0 : levrelax,nint(pmid(i,levrelax)/100._r8),'mb'
815 : end if
816 : #endif
817 0 : do k = kmax(i)+1,levrelax
818 0 : pint_vals(1) = mr_ub(i,m,last,kl(i,k),lchnk) &
819 : + delp(i,k) &
820 0 : * (mr_ub(i,m,last,ku(i,k),lchnk) &
821 0 : - mr_ub(i,m,last,kl(i,k),lchnk))
822 0 : pint_vals(2) = mr_ub(i,m,next,kl(i,k),lchnk) &
823 : + delp(i,k) &
824 : * (mr_ub(i,m,next,ku(i,k),lchnk) &
825 0 : - mr_ub(i,m,next,kl(i,k),lchnk))
826 : vmrrelax = pint_vals(1) &
827 0 : + dels * (pint_vals(2) - pint_vals(1))
828 0 : h2o(i,k) = h2o(i,k) + (vmrrelax - h2o(i,k)) * facrelax
829 : end do
830 : end do col_loop2
831 : end if h2o_overwrite
832 :
833 0 : end subroutine set_fstrat_h2o
834 :
835 : subroutine rebin( nsrc, ntrg, src_x, trg_x, src, trg )
836 : !---------------------------------------------------------------
837 : ! ... rebin src to trg
838 : !---------------------------------------------------------------
839 :
840 : !---------------------------------------------------------------
841 : ! ... dummy arguments
842 : !---------------------------------------------------------------
843 : integer, intent(in) :: nsrc ! dimension source array
844 : integer, intent(in) :: ntrg ! dimension target array
845 : real(r8), intent(in) :: src_x(nsrc+1) ! source coordinates
846 : real(r8), intent(in) :: trg_x(ntrg+1) ! target coordinates
847 : real(r8), intent(in) :: src(nsrc) ! source array
848 : real(r8), intent(out) :: trg(ntrg) ! target array
849 :
850 : !---------------------------------------------------------------
851 : ! ... local variables
852 : !---------------------------------------------------------------
853 : integer :: i
854 : integer :: si, si1
855 : integer :: sil, siu
856 : real(r8) :: y
857 : real(r8) :: sl, su
858 : real(r8) :: tl, tu
859 :
860 : !---------------------------------------------------------------
861 : ! ... check interval overlap
862 : !---------------------------------------------------------------
863 : ! if( trg_x(1) < src_x(1) .or. trg_x(ntrg+1) > src_x(nsrc+1) ) then
864 : ! write(iulog,*) 'rebin: target grid is outside source grid'
865 : ! write(iulog,*) ' target grid from ',trg_x(1),' to ',trg_x(ntrg+1)
866 : ! write(iulog,*) ' source grid from ',src_x(1),' to ',src_x(nsrc+1)
867 : ! call endrun
868 : ! end if
869 :
870 : do i = 1,ntrg
871 : tl = trg_x(i)
872 : if( tl < src_x(nsrc+1) ) then
873 : do sil = 1,nsrc+1
874 : if( tl <= src_x(sil) ) then
875 : exit
876 : end if
877 : end do
878 : tu = trg_x(i+1)
879 : do siu = 1,nsrc+1
880 : if( tu <= src_x(siu) ) then
881 : exit
882 : end if
883 : end do
884 : y = 0._r8
885 : sil = max( sil,2 )
886 : siu = min( siu,nsrc+1 )
887 : do si = sil,siu
888 : si1 = si - 1
889 : sl = max( tl,src_x(si1) )
890 : su = min( tu,src_x(si) )
891 : y = y + (su - sl)*src(si1)
892 : end do
893 : trg(i) = y/(trg_x(i+1) - trg_x(i))
894 : else
895 : trg(i) = 0._r8
896 : end if
897 : end do
898 :
899 : end subroutine rebin
900 :
901 : end module mo_fstrat
|
