Line data Source code
1 : module upper_bc
2 :
3 : !---------------------------------------------------------------------------------
4 : ! Module to compute the upper boundary conditions for temperature (dry static energy)
5 : ! and trace gases. Uses the MSIS model, and SNOE and TIME GCM and general prescribed UBC data.
6 : !---------------------------------------------------------------------------------
7 :
8 : use shr_kind_mod, only: r8 => shr_kind_r8
9 : use shr_kind_mod, only: cl => shr_kind_cl
10 : use shr_const_mod,only: grav => shr_const_g, & ! gravitational constant (m/s^2)
11 : kboltz => shr_const_boltz, & ! Boltzmann constant
12 : pi => shr_const_pi, & ! pi
13 : rEarth => shr_const_rearth ! Earth radius
14 : use ppgrid, only: pcols, pver, pverp
15 : use constituents, only: pcnst
16 : use cam_logfile, only: iulog
17 : use spmd_utils, only: masterproc
18 : use ref_pres, only: do_molec_diff, ptop_ref
19 : use shr_kind_mod, only: cx=>SHR_KIND_CX
20 : use cam_abortutils,only: endrun
21 : use cam_history, only: addfld, horiz_only, outfld, fieldname_len
22 :
23 : use upper_bc_file, only: upper_bc_file_readnl, upper_bc_file_specified, upper_bc_file_adv, upper_bc_file_get
24 : use infnan, only: nan, assignment(=)
25 :
26 : implicit none
27 : private
28 : save
29 : !
30 : ! Public interfaces
31 : !
32 : public :: ubc_readnl ! read namelist options for UBCs
33 : public :: ubc_init ! global initialization
34 : public :: ubc_timestep_init ! time step initialization
35 : public :: ubc_get_vals ! get ubc values for this step
36 : public :: ubc_get_flxs ! get ub fluxes for this step
37 : public :: ubc_fixed_conc ! returns true for constituents that have fixed UBC
38 : public :: ubc_fixed_temp ! true if temperature at upper boundary is fixed
39 :
40 : character(len=64) :: ubc_specifier(pcnst) = 'NOTSET'
41 :
42 : character(len=16) :: ubc_flds(pcnst) = 'NOTSET'
43 : character(len=16) :: ubc_file_spfr(pcnst) = ' '
44 : character(len=32) :: ubc_source(pcnst) = ' '
45 :
46 : integer :: n_fixed_mmr=0
47 : integer :: n_fixed_vmr=0
48 :
49 : integer, allocatable :: fixed_mmr_ndx(:)
50 : integer, allocatable :: fixed_vmr_ndx(:)
51 : real(r8), allocatable :: fixed_mmr(:)
52 : real(r8), allocatable :: fixed_vmr(:)
53 :
54 : integer :: num_infile = 0
55 : integer :: num_fixed = 0
56 : character(len=2), parameter :: msis_flds(5) = &
57 : (/ 'H ','N ','O ','O2','T ' /)
58 : character(len=2), parameter :: tgcm_flds(1) = &
59 : (/ 'H2' /)
60 : character(len=2), parameter :: snoe_flds(1) = &
61 : (/ 'NO' /)
62 :
63 : logical, protected :: ubc_fixed_temp =.false.
64 : logical :: msis_active =.false.
65 : logical :: tgcm_active =.false.
66 : logical :: snoe_active =.false.
67 :
68 : ! Namelist variables
69 : character(len=cl) :: snoe_ubc_file = 'NONE'
70 : real(r8) :: t_pert_ubc = 0._r8
71 : real(r8) :: no_xfac_ubc = 1._r8
72 :
73 : integer :: h_ndx=-1
74 : integer :: h_msis_ndx=-1, n_msis_ndx=-1, o_msis_ndx=-1, o2_msis_ndx=-1
75 :
76 : character(len=cl) :: tgcm_ubc_file = 'NONE'
77 : integer :: tgcm_ubc_cycle_yr = 0
78 : integer :: tgcm_ubc_fixed_ymd = 0
79 : integer :: tgcm_ubc_fixed_tod = 0
80 : character(len=32) :: tgcm_ubc_data_type = 'CYCLICAL'
81 :
82 : logical :: apply_upper_bc = .false.
83 :
84 : integer, allocatable :: file_spc_ndx(:)
85 : integer, allocatable :: spc_ndx(:)
86 : character(len=fieldname_len), allocatable :: hist_names(:)
87 :
88 : !================================================================================================
89 : contains
90 : !================================================================================================
91 :
92 : !-----------------------------------------------------------------------
93 : !-----------------------------------------------------------------------
94 1536 : subroutine ubc_readnl(nlfile)
95 : use namelist_utils, only : find_group_name
96 : use spmd_utils, only : mpicom, masterprocid, mpi_character, mpi_integer, mpi_real8
97 : use string_utils, only : to_lower
98 :
99 : character(len=*), intent(in) :: nlfile
100 : integer :: unitn, ierr, m, n, ndx_co, ndx_ar
101 :
102 : character(len=*), parameter :: prefix = 'ubc_readnl: '
103 :
104 : namelist /upper_bc_opts/ tgcm_ubc_file,tgcm_ubc_data_type,tgcm_ubc_cycle_yr,tgcm_ubc_fixed_ymd, &
105 : tgcm_ubc_fixed_tod, snoe_ubc_file, no_xfac_ubc, t_pert_ubc
106 : namelist /upper_bc_opts/ ubc_specifier
107 :
108 1536 : if (masterproc) then
109 : ! read namelist
110 2 : open( newunit=unitn, file=trim(nlfile), status='old' )
111 2 : call find_group_name(unitn, 'upper_bc_opts', status=ierr)
112 2 : if (ierr == 0) then
113 0 : read(unitn, upper_bc_opts, iostat=ierr)
114 0 : if (ierr /= 0) then
115 0 : call endrun(prefix//'upper_bc_opts: ERROR reading namelist')
116 : end if
117 : end if
118 2 : close(unitn)
119 :
120 : ! log the UBC options
121 2 : write(iulog,*) prefix//'tgcm_ubc_file = '//trim(tgcm_ubc_file)
122 2 : write(iulog,*) prefix//'tgcm_ubc_data_type = '//trim(tgcm_ubc_data_type)
123 2 : write(iulog,*) prefix//'tgcm_ubc_cycle_yr = ', tgcm_ubc_cycle_yr
124 2 : write(iulog,*) prefix//'tgcm_ubc_fixed_ymd = ', tgcm_ubc_fixed_ymd
125 2 : write(iulog,*) prefix//'tgcm_ubc_fixed_tod = ', tgcm_ubc_fixed_tod
126 2 : write(iulog,*) prefix//'snoe_ubc_file = '//trim(snoe_ubc_file)
127 2 : write(iulog,*) prefix//'t_pert_ubc = ', t_pert_ubc
128 2 : write(iulog,*) prefix//'no_xfac_ubc = ', no_xfac_ubc
129 2 : write(iulog,*) prefix//'ubc_specifier : '
130 :
131 2 : n=1
132 2 : m=1
133 2 : do while(ubc_specifier(n)/='NOTSET')
134 0 : write(iulog,'(i4,a)') n,' '//trim(ubc_specifier(n))
135 :
136 0 : ndx_ar = index(ubc_specifier(n),'->')
137 :
138 0 : if (ndx_ar<1) then
139 0 : call endrun(prefix//'ubc_specifier "'//trim(ubc_specifier(n))//'" must include "->"')
140 : endif
141 :
142 0 : ubc_source(n) = trim(to_lower(adjustl(ubc_specifier(n)(ndx_ar+2:))))
143 :
144 0 : if (trim(ubc_source(n))=='ubc_file') then
145 0 : ubc_file_spfr(m) = trim(ubc_specifier(n)(:ndx_ar-1))
146 0 : m=m+1
147 : endif
148 0 : if (index(ubc_source(n),'mmr')>0) then
149 0 : n_fixed_mmr=n_fixed_mmr+1
150 0 : else if (index(ubc_source(n),'vmr')>0) then
151 0 : n_fixed_vmr=n_fixed_vmr+1
152 : end if
153 :
154 0 : ndx_co = index(ubc_specifier(n),':')
155 :
156 0 : if (ndx_co>0) then
157 0 : ubc_flds(n) = ubc_specifier(n)(:ndx_co-1)
158 : else
159 0 : ubc_flds(n) = ubc_specifier(n)(:ndx_ar-1)
160 : end if
161 :
162 0 : n=n+1
163 : end do
164 2 : num_fixed=n-1
165 2 : num_infile=m-1
166 : end if
167 :
168 :
169 : ! broadcast to all MPI tasks
170 1536 : call mpi_bcast(num_fixed, 1, mpi_integer, masterprocid, mpicom, ierr)
171 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : num_fixed')
172 1536 : call mpi_bcast(num_infile, 1, mpi_integer, masterprocid, mpicom, ierr)
173 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : num_infile')
174 1536 : call mpi_bcast(n_fixed_mmr, 1, mpi_integer, masterprocid, mpicom, ierr)
175 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : n_fixed_mmr')
176 1536 : call mpi_bcast(n_fixed_vmr, 1, mpi_integer, masterprocid, mpicom, ierr)
177 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : n_fixed_vmr')
178 1536 : call mpi_bcast(tgcm_ubc_file, len(tgcm_ubc_file), mpi_character, masterprocid, mpicom, ierr)
179 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : tgcm_ubc_file')
180 1536 : call mpi_bcast(tgcm_ubc_data_type, len(tgcm_ubc_data_type),mpi_character, masterprocid, mpicom, ierr)
181 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : tgcm_ubc_data_type')
182 1536 : call mpi_bcast(tgcm_ubc_cycle_yr, 1, mpi_integer, masterprocid, mpicom, ierr)
183 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : tgcm_ubc_cycle_yr')
184 1536 : call mpi_bcast(tgcm_ubc_fixed_ymd, 1, mpi_integer, masterprocid, mpicom, ierr)
185 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : tgcm_ubc_fixed_ymd')
186 1536 : call mpi_bcast(tgcm_ubc_fixed_tod, 1, mpi_integer, masterprocid, mpicom, ierr)
187 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : tgcm_ubc_fixed_tod')
188 1536 : call mpi_bcast(snoe_ubc_file, len(snoe_ubc_file), mpi_character, masterprocid, mpicom, ierr)
189 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : snoe_ubc_file')
190 1536 : call mpi_bcast(t_pert_ubc, 1, mpi_real8, masterprocid, mpicom, ierr)
191 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : t_pert_ubc')
192 1536 : call mpi_bcast(no_xfac_ubc,1, mpi_real8, masterprocid, mpicom, ierr)
193 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : no_xfac_ubc')
194 1536 : call mpi_bcast(ubc_specifier, pcnst*len(ubc_specifier(1)), mpi_character,masterprocid, mpicom, ierr)
195 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : ubc_specifier')
196 1536 : call mpi_bcast(ubc_flds, pcnst*len(ubc_flds(1)), mpi_character,masterprocid, mpicom, ierr)
197 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : ubc_flds')
198 1536 : call mpi_bcast(ubc_file_spfr, pcnst*len(ubc_file_spfr(1)), mpi_character,masterprocid, mpicom, ierr)
199 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : ubc_file_spfr')
200 1536 : call mpi_bcast(ubc_source, pcnst*len(ubc_source(1)), mpi_character,masterprocid, mpicom, ierr)
201 1536 : if (ierr /= 0) call endrun(prefix//'mpi_bcast error : ubc_source')
202 :
203 1536 : apply_upper_bc = num_fixed>0
204 :
205 1536 : if (apply_upper_bc) then
206 0 : call upper_bc_file_readnl(nlfile)
207 : end if
208 :
209 1536 : end subroutine ubc_readnl
210 :
211 : !===============================================================================
212 :
213 1536 : subroutine ubc_init()
214 : !-----------------------------------------------------------------------
215 : ! Initialization of time independent fields for the upper boundary condition
216 : ! Calls initialization routine for MSIS, TGCM and SNOE
217 : !-----------------------------------------------------------------------
218 : use mo_tgcm_ubc, only: tgcm_ubc_inti
219 : use mo_snoe, only: snoe_inti
220 : use mo_msis_ubc, only: msis_ubc_inti
221 : use constituents,only: cnst_get_ind
222 : use upper_bc_file, only: upper_bc_file_init
223 :
224 : !---------------------------Local workspace-----------------------------
225 : logical, parameter :: zonal_avg = .false.
226 : integer :: m, mm, ierr
227 : integer :: mmrndx, vmrndx, m_mmr, m_vmr
228 :
229 : real(r8) :: val
230 : character(len=32) :: str
231 : character(len=*), parameter :: prefix = 'ubc_init: '
232 : !-----------------------------------------------------------------------
233 :
234 1536 : call cnst_get_ind('H', h_ndx, abort=.false.) ! for H fluxes UBC (WACCMX)
235 :
236 1536 : if (.not.apply_upper_bc) return
237 :
238 0 : if (num_infile>0) then
239 0 : call upper_bc_file_init( ubc_file_spfr(:num_infile) )
240 : endif
241 :
242 : ! possible MSIS, TGCM, SNOE and ubc_file inputs
243 :
244 0 : mm=1
245 :
246 0 : allocate(hist_names(num_fixed), stat=ierr)
247 0 : if (ierr /= 0) call endrun(prefix//'allocate error : hist_names')
248 0 : allocate(spc_ndx(num_fixed), stat=ierr)
249 0 : if (ierr /= 0) call endrun(prefix//'allocate error : spc_ndx')
250 0 : spc_ndx=-1
251 0 : if (num_infile>0) then
252 0 : allocate(file_spc_ndx(num_infile), stat=ierr)
253 0 : if (ierr /= 0) call endrun(prefix//'allocate error : file_spc_ndx')
254 0 : file_spc_ndx=-1
255 : end if
256 0 : if (n_fixed_mmr>0) then
257 0 : allocate(fixed_mmr_ndx(n_fixed_mmr), stat=ierr)
258 0 : if (ierr /= 0) call endrun(prefix//'allocate error : fixed_mmr_ndx')
259 0 : allocate(fixed_mmr(n_fixed_mmr), stat=ierr)
260 0 : if (ierr /= 0) call endrun(prefix//'allocate error : fixed_mmr')
261 0 : fixed_mmr_ndx=-1
262 0 : fixed_mmr = nan
263 : end if
264 0 : if (n_fixed_vmr>0) then
265 0 : allocate(fixed_vmr_ndx(n_fixed_vmr), stat=ierr)
266 0 : if (ierr /= 0) call endrun(prefix//'allocate error : fixed_vmr_ndx')
267 0 : allocate(fixed_vmr(n_fixed_vmr), stat=ierr)
268 0 : if (ierr /= 0) call endrun(prefix//'allocate error : fixed_vmr')
269 0 : fixed_vmr_ndx=-1
270 0 : fixed_vmr = nan
271 : end if
272 :
273 0 : m_mmr = 0
274 0 : m_vmr = 0
275 :
276 0 : do m = 1,num_fixed
277 0 : hist_names(m) = trim(ubc_flds(m))//'_UBC'
278 0 : if (ubc_flds(m)=='T') then
279 0 : ubc_fixed_temp=.true.
280 0 : spc_ndx(m) = -1
281 0 : call addfld(hist_names(m), horiz_only, 'I', 'K', trim(ubc_flds(m))//' at upper boundary' )
282 : else
283 0 : call cnst_get_ind(ubc_flds(m), spc_ndx(m), abort=.true.)
284 0 : call addfld(hist_names(m), horiz_only, 'I', 'kg/kg', trim(ubc_flds(m))//' at upper boundary' )
285 : end if
286 :
287 0 : if (trim(ubc_source(m))=='msis') then
288 0 : if (do_molec_diff .and. any(msis_flds==ubc_flds(m))) then
289 0 : msis_active = .true.
290 0 : if (trim(ubc_flds(m))=='H') h_msis_ndx=spc_ndx(m)
291 0 : if (trim(ubc_flds(m))=='N') n_msis_ndx=spc_ndx(m)
292 0 : if (trim(ubc_flds(m))=='O') o_msis_ndx=spc_ndx(m)
293 0 : if (trim(ubc_flds(m))=='O2') o2_msis_ndx=spc_ndx(m)
294 : else
295 0 : call endrun(prefix//'MSIS is not allowed in this configuration')
296 : end if
297 0 : else if (trim(ubc_source(m))=='tgcm') then
298 0 : if (do_molec_diff .and. any(tgcm_flds==ubc_flds(m))) then
299 0 : tgcm_active = .true.
300 : else
301 0 : call endrun(prefix//'TGCM is not allowed in this configuration')
302 : end if
303 0 : else if (trim(ubc_source(m))=='snoe') then
304 0 : if (do_molec_diff .and. any(snoe_flds==ubc_flds(m))) then
305 0 : snoe_active = .true.
306 : else
307 0 : call endrun(prefix//'SNOE is not allowed in this configuration')
308 : end if
309 0 : else if (trim(ubc_source(m))=='ubc_file') then
310 0 : file_spc_ndx(mm) = spc_ndx(m)
311 0 : mm = mm+1
312 : else
313 0 : mmrndx = index(trim(ubc_source(m)),'mmr')
314 0 : vmrndx = index(trim(ubc_source(m)),'vmr')
315 0 : if (mmrndx>0 .and. vmrndx>0) then
316 0 : call endrun(prefix//'incorrect units in UBC source: '//trim(ubc_source(m)))
317 : end if
318 0 : if (mmrndx>0) then
319 0 : str = ubc_source(m)(:mmrndx-1)
320 0 : read(str,*) val
321 0 : m_mmr = m_mmr + 1
322 0 : fixed_mmr(m_mmr) = val
323 0 : fixed_mmr_ndx(m_mmr) = spc_ndx(m)
324 0 : else if (vmrndx>0) then
325 0 : str = ubc_source(m)(:vmrndx-1)
326 0 : read(str,*) val
327 0 : m_vmr = m_vmr + 1
328 0 : fixed_vmr(m_vmr) = val
329 0 : fixed_vmr_ndx(m_vmr) = spc_ndx(m)
330 : else
331 0 : call endrun(prefix//'unrecognized UBC source: '//trim(ubc_source(m)))
332 : end if
333 : end if
334 : end do
335 :
336 0 : if (tgcm_active) then
337 : !-----------------------------------------------------------------------
338 : ! ... initialize the tgcm upper boundary module
339 : !-----------------------------------------------------------------------
340 : call tgcm_ubc_inti( tgcm_ubc_file, tgcm_ubc_data_type, tgcm_ubc_cycle_yr, &
341 0 : tgcm_ubc_fixed_ymd, tgcm_ubc_fixed_tod)
342 0 : if (masterproc) write(iulog,*) 'ubc_init: after tgcm_ubc_inti'
343 : endif
344 :
345 0 : if (snoe_active) then
346 : !-----------------------------------------------------------------------
347 : ! ... initialize the snoe module
348 : !-----------------------------------------------------------------------
349 0 : call snoe_inti(snoe_ubc_file)
350 0 : if (masterproc) write(iulog,*) 'ubc_init: after snoe_inti'
351 : endif
352 :
353 0 : if (msis_active) then
354 : !-----------------------------------------------------------------------
355 : ! ... initialize the msis module
356 : !-----------------------------------------------------------------------
357 0 : call msis_ubc_inti( zonal_avg, n_msis_ndx,h_msis_ndx,o_msis_ndx,o2_msis_ndx )
358 0 : if (masterproc) write(iulog,*) 'ubc_init: after msis_ubc_inti'
359 : endif
360 :
361 1536 : end subroutine ubc_init
362 :
363 : !===============================================================================
364 : !===============================================================================
365 :
366 0 : pure logical function ubc_fixed_conc(name)
367 :
368 : character(len=*), intent(in) :: name
369 :
370 : integer :: m
371 :
372 0 : ubc_fixed_conc = .false.
373 :
374 0 : do m = 1,num_fixed
375 0 : if ( trim(ubc_flds(m)) == trim(name) ) then
376 0 : ubc_fixed_conc = .true.
377 0 : return
378 : endif
379 : end do
380 :
381 1536 : end function ubc_fixed_conc
382 :
383 : !===============================================================================
384 :
385 741888 : subroutine ubc_timestep_init(pbuf2d, state)
386 : !-----------------------------------------------------------------------
387 : ! timestep dependent setting
388 : !-----------------------------------------------------------------------
389 :
390 : use solar_parms_data, only: kp=>solar_parms_kp, ap=>solar_parms_ap, f107=>solar_parms_f107
391 : use solar_parms_data, only: f107a=>solar_parms_f107a, f107p=>solar_parms_f107p
392 : use mo_msis_ubc, only: msis_timestep_init
393 : use mo_tgcm_ubc, only: tgcm_timestep_init
394 : use mo_snoe, only: snoe_timestep_init
395 : use physics_types, only: physics_state
396 : use ppgrid, only: begchunk, endchunk
397 : use physics_buffer, only: physics_buffer_desc
398 :
399 : type(physics_state), intent(in) :: state(begchunk:endchunk)
400 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
401 :
402 370944 : if (.not.apply_upper_bc) return
403 :
404 0 : if (num_infile>0) then
405 0 : call upper_bc_file_adv( pbuf2d, state )
406 : end if
407 0 : if (msis_active) then
408 0 : call msis_timestep_init( ap, f107p, f107a )
409 : end if
410 0 : if (tgcm_active) then
411 0 : call tgcm_timestep_init( pbuf2d, state )
412 : end if
413 0 : if (snoe_active) then
414 0 : call snoe_timestep_init( kp, f107 )
415 : end if
416 :
417 370944 : end subroutine ubc_timestep_init
418 :
419 : !===============================================================================
420 :
421 1489176 : subroutine ubc_get_vals (lchnk, ncol, pint, zi, ubc_temp, ubc_mmr)
422 :
423 : !-----------------------------------------------------------------------
424 : ! interface routine for vertical diffusion and pbl scheme
425 : !-----------------------------------------------------------------------
426 370944 : use mo_msis_ubc, only: get_msis_ubc
427 : use mo_snoe, only: set_no_ubc, ndx_no
428 : use mo_tgcm_ubc, only: set_tgcm_ubc
429 : use cam_abortutils, only: endrun
430 : use air_composition, only: rairv, mbarv ! gas constant, mean mass
431 : use constituents, only: cnst_mw ! Needed for ubc_flux
432 :
433 : !------------------------------Arguments--------------------------------
434 : integer, intent(in) :: lchnk ! chunk identifier
435 : integer, intent(in) :: ncol ! number of atmospheric columns
436 : real(r8), intent(in) :: pint(pcols,pverp) ! interface pressures
437 : real(r8), intent(in) :: zi(pcols,pverp) ! interface geoptl height above sfc
438 :
439 : real(r8), intent(out) :: ubc_temp(pcols) ! upper bndy temperature (K)
440 : real(r8), intent(out) :: ubc_mmr(pcols,pcnst) ! upper bndy mixing ratios (kg/kg)
441 :
442 : !---------------------------Local storage-------------------------------
443 : integer :: m ! constituent index
444 : real(r8) :: rho_top(pcols) ! density at top interface
445 : real(r8) :: z_top(pcols) ! height of top interface (km)
446 2978352 : real(r8) :: vals(pcols,num_infile)
447 :
448 : real(r8), parameter :: m2km = 1.e-3_r8 ! meter to km
449 :
450 : !-----------------------------------------------------------------------
451 :
452 1489176 : ubc_mmr(:,:) = 0._r8
453 1489176 : ubc_temp(:) = nan
454 :
455 1489176 : if (.not. apply_upper_bc) return
456 :
457 : ! UBC_FILE
458 0 : if (num_infile>0) then
459 0 : call upper_bc_file_get(lchnk, ncol, vals)
460 0 : do m = 1,num_infile
461 0 : if (file_spc_ndx(m)>0) then
462 0 : ubc_mmr(:ncol,file_spc_ndx(m)) = vals(:ncol,m)
463 : else
464 0 : ubc_temp(:ncol) = vals(:ncol,m)
465 : end if
466 :
467 : end do
468 :
469 : endif
470 :
471 : ! MSIS
472 0 : if (msis_active) then
473 0 : call get_msis_ubc( lchnk, ncol, ubc_temp, ubc_mmr )
474 0 : if( t_pert_ubc /= 0._r8 ) then
475 0 : ubc_temp(:ncol) = ubc_temp(:ncol) + t_pert_ubc
476 0 : if( any( ubc_temp(:ncol) < 0._r8 ) ) then
477 0 : write(iulog,*) 'ubc_get_vals: msis temp < 0 after applying offset = ',t_pert_ubc
478 0 : call endrun('ubc_get_vals: msis temp < 0 after applying t_pert_ubc')
479 : end if
480 : end if
481 : end if
482 :
483 : ! SNOE
484 0 : if (snoe_active) then
485 :
486 0 : rho_top(:ncol) = pint(:ncol,1) / (rairv(:ncol,1,lchnk)*ubc_temp(:ncol))
487 0 : z_top(:ncol) = m2km * zi(:ncol,1)
488 :
489 0 : call set_no_ubc ( lchnk, ncol, z_top, ubc_mmr, rho_top )
490 0 : if( ndx_no > 0 .and. no_xfac_ubc /= 1._r8 ) then
491 0 : ubc_mmr(:ncol,ndx_no) = no_xfac_ubc * ubc_mmr(:ncol,ndx_no)
492 : end if
493 :
494 : endif
495 :
496 : ! TIE-GCM
497 0 : if (tgcm_active) then
498 0 : call set_tgcm_ubc( lchnk, ncol, ubc_mmr )
499 : endif
500 :
501 : ! fixed values
502 0 : do m = 1,n_fixed_mmr
503 0 : ubc_mmr(:ncol,fixed_mmr_ndx(m)) = fixed_mmr(m)
504 : end do
505 0 : do m = 1,n_fixed_vmr
506 0 : ubc_mmr(:ncol,fixed_vmr_ndx(m)) = cnst_mw(fixed_vmr_ndx(m))*fixed_vmr(m)/mbarv(:ncol,1,lchnk)
507 : end do
508 :
509 : ! diagnostic output
510 0 : do m = 1,num_fixed
511 0 : if (ubc_flds(m)=='T') then
512 0 : call outfld(hist_names(m),ubc_temp(:ncol),ncol,lchnk)
513 : else
514 0 : call outfld(hist_names(m),ubc_mmr(:ncol,spc_ndx(m)),ncol,lchnk)
515 : end if
516 : end do
517 :
518 1489176 : end subroutine ubc_get_vals
519 :
520 : !===============================================================================
521 :
522 0 : subroutine ubc_get_flxs (lchnk, ncol, pint, zi, t, q, phis, ubc_flux)
523 :
524 1489176 : use physconst, only: avogad, rga
525 : use air_composition, only: rairv
526 : use constituents, only: cnst_mw
527 : !------------------------------Arguments--------------------------------
528 : integer, intent(in) :: lchnk ! chunk identifier
529 : integer, intent(in) :: ncol ! number of atmospheric columns
530 : real(r8), intent(in) :: pint(pcols,pverp) ! interface pressures
531 : real(r8), intent(in) :: zi(pcols,pverp) ! interface geoptl height above sfc
532 : real(r8), intent(in) :: t(pcols,pver) ! midpoint temperature
533 : real(r8), intent(in),target :: q(pcols,pver,pcnst) ! contituent mixing ratios (kg/kg)
534 : real(r8), intent(in) :: phis(pcols) ! Surface geopotential (m2/s2)
535 :
536 : real(r8), intent(out) :: ubc_flux(pcols,pcnst) ! upper bndy flux (kg/s/m^2)
537 :
538 : !---------------------------Local storage-------------------------------
539 : integer :: iCol ! column loop counter
540 :
541 : real(r8), parameter :: h_escape_flx_factor = 2.03e-13_r8 ! for hydrogen escape flux due to charge exchange
542 : ! adopted from TIME-GCM (R. G. Roble, pp. 1-21, AGU Geophys. Monogr. Ser 87, 1995) following
543 : ! Liu, S.C., and T. M. Donahue, Mesospheric hydrogen related to exospheric escape mechanisms, J. Atmos. Sci.,
544 : ! 31, 1466-1470, 1974. (Equation 4 there). DOI: 10.1175/1520-0469(1974)031<1466:Mhrtee>2.0.Co;2
545 : ! https://journals.ametsoc.org/view/journals/atsc/31/5/1520-0469_1974_031_1466_mhrtee_2_0_co_2.xml
546 :
547 : real(r8), parameter :: hfluxlimitfac = 0.72_r8 ! Hydrogen upper boundary flux limiting factor
548 :
549 : real(r8) :: nmbartop ! Top level density (rho)
550 : real(r8) :: zkt ! Factor for H Jean's escape flux calculation
551 :
552 0 : real(r8), pointer :: qh_top(:) ! Top level hydrogen mixing ratio (kg/kg)
553 :
554 0 : ubc_flux(:,:) = nan
555 :
556 0 : qh_top => q(:,1,h_ndx)
557 :
558 0 : do iCol = 1, ncol
559 : !--------------------------------------------------
560 : ! Get total density (rho) at top level
561 : !--------------------------------------------------
562 0 : nmbartop = 0.5_r8 * (pint(iCol,1) + pint(iCol,2)) / ( rairv(iCol,1,lchnk) * t(iCol,1) )
563 :
564 : !---------------------------------------------------------------------
565 : ! Calculate factor for Jean's escape flux once here, used twice below
566 : !---------------------------------------------------------------------
567 : zkt = (rEarth + ( 0.5_r8 * ( zi(iCol,1) + zi(iCol,2) ) + rga * phis(iCol) ) ) * &
568 0 : cnst_mw(h_ndx) / avogad * grav / ( kboltz * t(iCol,1) )
569 :
570 : ubc_flux(iCol,h_ndx) = hfluxlimitfac * SQRT(kboltz/(2.0_r8 * pi * cnst_mw(h_ndx) / avogad)) * &
571 0 : qh_top(iCol) * nmbartop * &
572 0 : SQRT(t(iCol,1)) * (1._r8 + zkt) * EXP(-zkt)
573 :
574 : ubc_flux(iCol,h_ndx) = ubc_flux(iCol,h_ndx) * &
575 0 : (h_escape_flx_factor * qh_top(iCol) * nmbartop / (cnst_mw(h_ndx) / avogad) * t(iCol,1))
576 :
577 : enddo
578 :
579 0 : end subroutine ubc_get_flxs
580 :
581 : end module upper_bc
|
