Line data Source code
1 : ! ====================================================================================
2 : ! Air-sea exchange module for trace gases
3 : ! Ref: Carpenter et al Chem Soc Rev (2012); Johnson, Ocean sci (2010)
4 : ! ------------------------------------------------------------------------------------
5 : ! Required inputs for the air-sea flux module:
6 : ! - Seawater concentration (nanomoles per liter) and Sea surface salinity
7 : ! (parts per thousand) read from namelist (netCDF)
8 : ! - Concentration in the gas-phase (pptv), air temperature (K), 10m windspeed (m/s),
9 : ! surface pressure (atm), sea surface temperature (K): all from other modules
10 : ! ------------------------------------------------------------------------------------
11 : ! Key subroutines:
12 : ! ocean_emis_readnl(..): Read salinity from namelist (user_nl_cam).
13 : ! Salinity not time-dependent. Flux depends very weakly on it
14 : ! ocean_emis_init(...): Interpolate salinity, initialize the library for the flux
15 : ! reading time-dependent seawater conc. from user_nl_cam
16 : ! ocean_emis_advance(...): process the seawater concentration
17 : ! ocean_emis_getflux(...): calculate the air-sea flux (upward or downward),
18 : ! then add to total surface flux (sflx)
19 : ! ------------------------------------------------------------------------------------
20 : ! Last built: 9 March 2018.
21 : ! Written by: Siyuan Wang (ACOM/NCAR) siyuan@ucar.edu
22 : ! Acknowledgement: Francis Vitt (NCAR). and of course Dr. Peppurr too
23 : ! ====================================================================================
24 :
25 : module ocean_emis
26 :
27 : use shr_kind_mod, only : r8=>shr_kind_r8, cl=>shr_kind_cl
28 : use ppgrid, only : pcols, begchunk,endchunk
29 : use spmd_utils, only : masterproc
30 : use cam_abortutils, only : endrun
31 : use cam_history, only : addfld, add_default, horiz_only, outfld
32 : use constituents, only : cnst_get_ind
33 : use tracer_data, only : trfld,trfile
34 : use chem_mods, only : gas_pcnst
35 : use cam_logfile, only : iulog
36 : use ioFileMod, only : getfil
37 :
38 : implicit none
39 :
40 : private
41 : public :: ocean_emis_readnl
42 : public :: ocean_emis_init
43 : public :: ocean_emis_getflux
44 : public :: ocean_emis_advance
45 : public :: ocean_emis_species
46 :
47 : type :: Csw
48 : integer :: spc_ndx
49 : real(r8) :: scalefactor
50 : character(len=256) :: filename
51 : character(len=16) :: species
52 : character(len=8) :: units
53 : integer :: nsectors
54 : character(len=32),pointer :: sectors(:)
55 : type(trfld), pointer :: fields(:)
56 : type(trfile) :: file
57 : end type Csw
58 :
59 : logical :: switch_bubble
60 : type(Csw), allocatable :: Csw_nM(:)
61 : integer :: n_Csw_files = 0
62 :
63 : real(r8), allocatable :: salinity(:,:)
64 :
65 : ! ================
66 : ! Air-sea exchange
67 : ! ================
68 : Integer, Parameter :: HowManyMolecules = 16 ! Change this number if you wanna add more species
69 : Integer, Parameter :: HowManyProperties = 16 ! Don't touch this (unless you wanna add more fields)
70 : Integer, Parameter :: HowManySalts = 5 ! Change this number if you wanna add more salts
71 : Integer, Parameter :: HowManySaltProperties = 7 ! Don't touch this (unless you wanna add more fields)
72 :
73 : Type GasLib
74 : Character(16) :: CmpdName
75 : Real(r8), Dimension(HowManyProperties) :: CmpdProperties
76 : End Type GasLib
77 :
78 : Type SaltLib
79 : Character(16) :: SaltName
80 : Real(r8), Dimension(HowManySaltProperties) :: SaltProperties
81 : End Type SaltLib
82 :
83 : Type(GasLib), Dimension(HowManyMolecules) :: GasList ! Library for the trace gas properties
84 : Type(SaltLib), Dimension(HowManySalts) :: SaltList ! Library for the salt properties
85 :
86 : ! ===========================
87 : ! seawater concentration:
88 : ! ===========================
89 : character(len=cl) :: csw_specifier(gas_pcnst) = ''
90 : character(len=24) :: csw_time_type = 'CYCLICAL' ! 'CYCLICAL' | 'SERIAL' | 'INTERP_MISSING_MONTHS'
91 : integer :: csw_cycle_yr = 0
92 : logical :: bubble_mediated_transfer = .false.
93 : character(len=cl) :: ocean_salinity_file = 'NONE'
94 :
95 : contains
96 :
97 : !--------------------------------------------------------------------------------
98 : !--------------------------------------------------------------------------------
99 1536 : subroutine ocean_emis_readnl(nlfile)
100 :
101 : use namelist_utils, only : find_group_name
102 : use spmd_utils, only : mpicom, masterprocid, mpi_character, mpi_integer, mpi_logical
103 :
104 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
105 :
106 : integer :: unitn, ierr
107 : character(len=*), parameter :: subname = 'ocean_emis_readnl'
108 :
109 : ! ===================
110 : ! Namelist definition
111 : ! ===================
112 : namelist /ocean_emis_nl/ ocean_salinity_file
113 : namelist /ocean_emis_nl/ csw_specifier, csw_time_type, csw_cycle_yr, bubble_mediated_transfer
114 :
115 : ! =============
116 : ! Read namelist
117 : ! =============
118 1536 : if (masterproc) then
119 2 : open( newunit=unitn, file=trim(nlfile), status='old' )
120 2 : call find_group_name(unitn, 'ocean_emis_nl', status=ierr)
121 2 : if (ierr == 0) then
122 2 : read(unitn, ocean_emis_nl, iostat=ierr)
123 2 : if (ierr /= 0) then
124 0 : call endrun(subname // ':: ERROR reading namelist')
125 : end if
126 : end if
127 2 : close(unitn)
128 : end if
129 :
130 : ! ============================
131 : ! Broadcast namelist variables
132 : ! ============================
133 1536 : call mpi_bcast(ocean_salinity_file, len(ocean_salinity_file), mpi_character, masterprocid, mpicom, ierr)
134 1536 : call mpi_bcast(csw_specifier, cl*gas_pcnst, mpi_character, masterprocid, mpicom, ierr)
135 1536 : call mpi_bcast(csw_time_type, len(csw_time_type), mpi_character, masterprocid, mpicom, ierr)
136 1536 : call mpi_bcast(csw_cycle_yr, 1, mpi_integer, masterprocid, mpicom, ierr)
137 1536 : call mpi_bcast(bubble_mediated_transfer, 1, mpi_logical, masterprocid, mpicom, ierr)
138 :
139 1536 : if (masterproc) then
140 2 : write(iulog,*) 'ocean_emis_readnl: ocean_salinity_file = '//trim(ocean_salinity_file)
141 2 : write(iulog,*) 'ocean_emis_readnl: bubble_mediated_transfer = ',bubble_mediated_transfer
142 : end if
143 :
144 1536 : end subroutine ocean_emis_readnl
145 :
146 : !--------------------------------------------------------------------------------
147 : !--------------------------------------------------------------------------------
148 1536 : subroutine ocean_emis_init()
149 :
150 : use ioFileMod, only : getfil
151 : use cam_pio_utils, only : cam_pio_openfile
152 : use pio, only : file_desc_t, pio_inq_dimid, pio_inq_dimlen, pio_inq_varid, pio_get_var
153 : use pio, only : PIO_NOWRITE, PIO_NOERR
154 : use pio, only : pio_seterrorhandling, PIO_BCAST_ERROR, pio_closefile
155 : use phys_grid, only : get_ncols_p, get_rlon_all_p, get_rlat_all_p
156 : use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish
157 : use mo_constants, only : pi
158 :
159 : integer :: file_nlat, file_nlon
160 : integer :: dimid, vid, ierr, c, ncols
161 : type(file_desc_t) :: fid
162 : character(len=cl) :: filen
163 1536 : real(r8), allocatable :: file_lats(:), file_lons(:)
164 1536 : real(r8), allocatable :: wrk2d(:,:)
165 : real(r8) :: to_lats(pcols), to_lons(pcols)
166 : type(interp_type) :: lon_wgts, lat_wgts
167 :
168 : real(r8), parameter :: zero=0_r8, twopi=2_r8*pi, degs2rads = pi/180._r8
169 :
170 : character(len=*), parameter :: subname = 'ocean_emis_init'
171 :
172 1536 : if (trim(ocean_salinity_file) == 'NONE') return
173 :
174 1536 : call getfil( ocean_salinity_file, filen, 0 )
175 1536 : call cam_pio_openfile( fid, filen, PIO_NOWRITE)
176 :
177 1536 : call pio_seterrorhandling(fid, PIO_BCAST_ERROR)
178 :
179 1536 : ierr = pio_inq_dimid( fid, 'lon', dimid )
180 1536 : if (ierr /= PIO_NOERR) then
181 0 : call endrun(subname//': pio_inq_dimid lon FAILED')
182 : endif
183 1536 : ierr = pio_inq_dimlen( fid, dimid, file_nlon )
184 1536 : if (ierr /= PIO_NOERR) then
185 0 : call endrun(subname//': pio_inq_dimlen file_nlon FAILED')
186 : endif
187 :
188 4608 : allocate( file_lons(file_nlon) )
189 1536 : ierr = pio_inq_varid( fid, 'lon', vid )
190 1536 : if (ierr /= PIO_NOERR) then
191 0 : call endrun(subname//': pio_inq_varid lon FAILED')
192 : endif
193 1536 : ierr = pio_get_var( fid, vid, file_lons )
194 1536 : if (ierr /= PIO_NOERR) then
195 0 : call endrun(subname//': pio_get_var file_lons FAILED')
196 : endif
197 554496 : file_lons = file_lons * degs2rads
198 :
199 1536 : ierr = pio_inq_dimid( fid, 'lat', dimid )
200 1536 : if (ierr /= PIO_NOERR) then
201 0 : call endrun(subname//': pio_inq_dimid lat FAILED')
202 : endif
203 1536 : ierr = pio_inq_dimlen( fid, dimid, file_nlat )
204 1536 : if (ierr /= PIO_NOERR) then
205 0 : call endrun(subname//': pio_inq_dimlen file_nlat FAILED')
206 : endif
207 4608 : allocate( file_lats(file_nlat) )
208 1536 : ierr = pio_inq_varid( fid, 'lat', vid )
209 1536 : if (ierr /= PIO_NOERR) then
210 0 : call endrun(subname//': pio_inq_varid lat FAILED')
211 : endif
212 1536 : ierr = pio_get_var( fid, vid, file_lats )
213 1536 : if (ierr /= PIO_NOERR) then
214 0 : call endrun(subname//': pio_inq_varid SSS FAILED')
215 : endif
216 278016 : file_lats = file_lats * degs2rads
217 :
218 6144 : allocate( wrk2d( file_nlon, file_nlat ) )
219 1536 : ierr = pio_inq_varid( fid, 'SSS', vid )
220 1536 : if (ierr /= PIO_NOERR) then
221 0 : call endrun(subname//': pio_inq_varid SSS FAILED')
222 : endif
223 1536 : ierr = pio_get_var( fid, vid, wrk2d )
224 1536 : if (ierr /= PIO_NOERR) then
225 0 : call endrun(subname//': pio_get_var SSS FAILED')
226 : endif
227 :
228 4608 : allocate(salinity(pcols,begchunk:endchunk))
229 106800 : salinity = 0._r8
230 :
231 7728 : do c=begchunk,endchunk
232 :
233 6192 : ncols = get_ncols_p(c)
234 6192 : call get_rlat_all_p(c, pcols, to_lats)
235 6192 : call get_rlon_all_p(c, pcols, to_lons)
236 :
237 6192 : call lininterp_init(file_lons, file_nlon, to_lons, ncols, 2, lon_wgts, zero, twopi)
238 6192 : call lininterp_init(file_lats, file_nlat, to_lats, ncols, 1, lat_wgts)
239 :
240 6192 : call lininterp(wrk2d, file_nlon, file_nlat, salinity(1:ncols,c), ncols, lon_wgts, lat_wgts)
241 :
242 6192 : call lininterp_finish(lon_wgts)
243 7728 : call lininterp_finish(lat_wgts)
244 :
245 : end do
246 :
247 : ! fill in missing values with climatology for modern-day
248 106800 : where(salinity < 0._r8)
249 : salinity = 33.0_r8
250 : end where
251 :
252 1536 : deallocate( file_lons, file_lats )
253 1536 : deallocate( wrk2d )
254 :
255 1536 : call addfld('OCN_SALINITY', horiz_only, 'A', 'parts per thousands', 'ocean salinity' )
256 :
257 : ! ======================================================
258 : ! initializing the libraries for the air-sea flux module
259 : ! ======================================================
260 1536 : Call CmpLibInitialization()
261 1536 : Call SaltLibInitialization()
262 :
263 : ! ---------------------------------------------
264 : ! Read seawater concentration: WSY
265 : ! ---------------------------------------------
266 1536 : call cseawater_ini()
267 :
268 1536 : call pio_closefile (fid)
269 :
270 3072 : end subroutine ocean_emis_init
271 :
272 : !--------------------------------------------------------------------------------
273 : !--------------------------------------------------------------------------------
274 741888 : subroutine ocean_emis_advance( pbuf2d, state )
275 : ! -------------------------------
276 : ! check serial case for time span
277 : ! -------------------------------
278 :
279 1536 : use physics_types, only : physics_state
280 : use ppgrid, only : begchunk, endchunk
281 : use tracer_data, only : advance_trcdata
282 : use physics_buffer, only : physics_buffer_desc
283 :
284 : type(physics_state), intent(in) :: state(begchunk:endchunk)
285 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
286 : integer :: m
287 :
288 370944 : if (trim(ocean_salinity_file) == 'NONE') return
289 :
290 741888 : do m = 1,n_Csw_files
291 741888 : call advance_trcdata( Csw_nM(m)%fields, Csw_nM(m)%file, state, pbuf2d )
292 : end do
293 :
294 370944 : end subroutine ocean_emis_advance
295 :
296 : !--------------------------------------------------------------------------------
297 : !--------------------------------------------------------------------------------
298 1489176 : subroutine ocean_emis_getflux(lchnk, ncol, state, u10, sst, ocnfrac, icefrac, sflx)
299 :
300 370944 : use physics_types, only : physics_state
301 : use ppgrid, only : pver
302 :
303 : integer, intent(in) :: lchnk, ncol
304 : type(physics_state), target, intent(in) :: state ! Physics state variables
305 : real(r8), intent(in) :: u10(:) ! Wind speed at 10m
306 : real(r8), intent(in) :: sst(:) ! Sea surface temperature
307 : real(r8), intent(in) :: ocnfrac(:) ! Ocean fraction
308 : real(r8), intent(in) :: icefrac(:) ! Ice fraction
309 : real(r8), intent(inout) :: sflx(:,:) ! Surface emissions (kg/m^2/s)
310 :
311 : integer :: i, m, isec, SpeciesID
312 2978352 : real(r8) :: Csw_col(ncol)
313 : real(r8) :: MW_species
314 2978352 : real(r8) :: oceanflux_kg_m2_s(ncol)
315 :
316 1489176 : if (trim(ocean_salinity_file) == 'NONE') return
317 :
318 : ! ==================================================
319 : ! Get seawater concentrations and calculate the flux
320 : ! ==================================================
321 2978352 : Csw_loop : do m = 1, n_Csw_files
322 :
323 24865776 : Csw_col(:) = 0._r8
324 1489176 : isec = 1
325 24865776 : Csw_col(:ncol) = Csw_nM(m)%scalefactor*Csw_nM(m)%fields(isec)%data(:ncol,1,lchnk)
326 :
327 1489176 : MW_species = MolecularWeight(SpeciesIndex( Csw_nM(m)%species ))
328 :
329 1489176 : call cnst_get_ind( trim(Csw_nM(m)%species), SpeciesID, abort=.true. )
330 :
331 24865776 : oceanflux_kg_m2_s = 0.0_r8
332 :
333 24865776 : do i = 1,ncol
334 24865776 : if (ocnfrac(i) >= 0.2_r8 .and. Csw_col(i) >= 0._r8) then
335 : ! calculate flux only for ocean
336 : oceanflux_kg_m2_s(i) = Flux_kg_m2_s( &
337 0 : Csw_nM(m)%species, & ! name of species
338 0 : state%q(i,pver,SpeciesID) * (28.97_r8/MW_species) * 1.0e+12_r8, & ! air concentration (ppt)
339 : Csw_col(i), & ! sea water concentration (nM)
340 0 : state%t(i,pver), & ! air temperature (K)
341 15708747 : u10(i), & ! wind speed at 10m (m/s) <- should use this
342 0 : state%ps(i) / 101325.0_r8, & ! surface pressure (atm)
343 0 : sst(i), & ! sea surface temperautre (K)
344 0 : salinity(i,lchnk), & ! ocean salinity (parts per thousands)
345 31417494 : switch_bubble ) ! bubble-mediated transfer: on or off
346 : end if
347 : end do
348 :
349 : ! ===========================================================================
350 : ! Add the ocean flux to the other fluxes
351 : ! Make sure this ocean module is called after other surface emissions are set
352 : ! ===========================================================================
353 24865776 : sflx(:ncol,SpeciesID) = sflx(:ncol,SpeciesID) + oceanflux_kg_m2_s(:ncol) * ocnfrac(:ncol)
354 :
355 : ! ====================================
356 : ! Write stuff into the historial files
357 : ! ====================================
358 1489176 : call outfld('OCN_FLUX_' // trim(Csw_nM(m)%species), oceanflux_kg_m2_s(:ncol), ncol, lchnk)
359 2978352 : call outfld('Csw_' // trim(Csw_nM(m)%species), Csw_col(:ncol), ncol, lchnk)
360 :
361 : end do Csw_loop
362 :
363 1489176 : call outfld('OCN_SALINITY', salinity(:ncol,lchnk), ncol, lchnk)
364 :
365 1489176 : end subroutine ocean_emis_getflux
366 :
367 : !--------------------------------------------------------------------------------
368 : !--------------------------------------------------------------------------------
369 1536 : Subroutine CmpLibInitialization()
370 : ! =====================================================================================
371 : ! This is the lookup table for molecular weight, Vb, and Henry's law constant
372 : ! NOTE: IF YOU ADDED NEW SPECIES, REMEMBER TO UPDATE THIS -> HowManyMolecules
373 : ! -------------------------------------------------------------------------------------
374 : ! Col 1: molecular weight (g/mol)
375 : ! Col 2-10: number of C, H, N, O, S, Br, Cl, F, and I atoms in the formula, respectively
376 : ! Col 11-13: number of double bound, triple bond, and rings in the molecule, respectively
377 : ! Col 14: Vb (cm^3/mol) values if available.
378 : ! Col 15-16: Henry's law constant at 298 K (M/atm) and temperature dependency (in K^-1)
379 : ! -------------------------------------------------------------------------------------
380 : ! Ref: Henry's law constants from Sander ACP 2014 (compilation)
381 : ! =====================================================================================
382 : GasList(1) = GasLib('CH3OH', (/ 32.05_r8, 1.0_r8, 4.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, &
383 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 42.7_r8, 203.0_r8, 5645.0_r8 /))
384 : GasList(2) = GasLib('C2H5OH', (/ 46.07_r8, 2.0_r8, 6.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, &
385 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 190.0_r8, 6500.0_r8 /))
386 : GasList(3) = GasLib('CH2O', (/ 30.03_r8, 1.0_r8, 2.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, &
387 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 3230.0_r8, 7100.0_r8 /))
388 : GasList(4) = GasLib('CH3CHO', (/ 44.05_r8, 2.0_r8, 4.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, &
389 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 12.9_r8, 5890.0_r8/))
390 : GasList(5) = GasLib('PROPANAL', (/ 58.08_r8, 3.0_r8, 6.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, &
391 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 10.0_r8, 4330.0_r8 /))
392 : GasList(6) = GasLib('CH3COCH3', (/ 58.08_r8, 3.0_r8, 6.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, &
393 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 77.6_r8, 27.8_r8, 5530.0_r8 /))
394 : GasList(7) = GasLib('NO', (/ 30.01_r8, 0.0_r8, 0.0_r8, 1.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, &
395 26112 : 0.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 1.92e-3_r8, 1762.0_r8/)) ! NO
396 : GasList(8) = GasLib('HNO2', (/ 47.01_r8, 0.0_r8, 1.0_r8, 1.0_r8, 2.0_r8, 0.0_r8, 0.0_r8, &
397 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 49.0_r8, 4900.0_r8/))
398 : GasList(9) = GasLib('HCN', (/ 27.03_r8, 1.0_r8, 1.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, &
399 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 9.02_r8, 8258.0_r8/))
400 : GasList(10) = GasLib('CH3CN', (/ 41.05_r8, 2.0_r8, 3.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, &
401 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 55.5_r8, 52.8_r8, 3970.0_r8/))
402 : GasList(11) = GasLib('PAN', (/ 121.05_r8, 2.0_r8, 3.0_r8, 1.0_r8, 5.0_r8, 0.0_r8, 0.0_r8, &
403 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 2.8_r8, 5730.0_r8/))
404 : GasList(12) = GasLib('CH3BR', (/ 94.94_r8, 1.0_r8, 3.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 1.0_r8, &
405 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.16_r8, 3100.0_r8/))
406 : GasList(13) = GasLib('CH2BR2', (/ 173.88_r8, 1.0_r8, 2.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 2.0_r8, &
407 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 1.1_r8, 4000.0_r8/))
408 : GasList(14) = GasLib('CHBR3', (/ 252.73_r8, 1.0_r8, 1.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 3.0_r8, &
409 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 1.4_r8, 5000.0_r8/))
410 : GasList(15) = GasLib('DMS', (/ 62.13_r8, 2.0_r8, 6.0_r8, 0.0_r8, 0.0_r8, 1.0_r8, 0.0_r8, &
411 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.54_r8, 3460.0_r8/))
412 : GasList(16) = GasLib('CH3NO3', (/ 74.04_r8, 1.0_r8, 3.0_r8, 1.0_r8, 3.0_r8, 0.0_r8, 0.0_r8, &
413 26112 : 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 0.0_r8, 2.0_r8, 4700.0_r8/))
414 : ! --------------------------------------------------------------------------------
415 : ! ADD NEW SPECIES HERE (molecular weight, formula info, and Henry's law constants)
416 : ! --------------------------------------------------------------------------------
417 1489176 : End Subroutine CmpLibInitialization
418 :
419 : !--------------------------------------------------------------------------------
420 : !--------------------------------------------------------------------------------
421 1536 : Subroutine SaltLibInitialization()
422 : ! ================================================================================
423 : ! This is the lookup table for common solutes in seawater and the parameters to
424 : ! calculate the dynamic viscosity of seawater.
425 : ! You may add other solutes or change the mass fractions.
426 : ! --------------------------------------------------------------------------------
427 : ! Col 1: mass fraction of solute
428 : ! Col 2-7: some parameter
429 : ! --------------------------------------------------------------------------------
430 : ! Ref: typical mass fractions of solute: Laliberte (2007)
431 : ! parameterizations: Laliberte (2007) and Millero et al (2008)
432 : ! ================================================================================
433 12288 : SaltList(1) = SaltLib('NaCl', (/ 0.798_r8, 16.22_r8, 1.3229_r8, 1.4849_r8, 0.0074691_r8, 30.78_r8, 2.0583_r8 /))
434 12288 : SaltList(2) = SaltLib('KCl', (/ 0.022_r8, 6.4883_r8, 1.3175_r8, -0.7785_r8, 0.09272_r8, -1.3_r8, 2.0811_r8 /))
435 12288 : SaltList(3) = SaltLib('CaCl2', (/ 0.033_r8, 32.028_r8, 0.78792_r8, -1.1495_r8, 0.0026995_r8, 780860.0_r8, 5.8442_r8 /))
436 12288 : SaltList(4) = SaltLib('MgCl2', (/ 0.047_r8, 24.032_r8, 2.2694_r8, 3.7108_r8, 0.021853_r8, -1.1236_r8, 0.14474_r8/))
437 12288 : SaltList(5) = SaltLib('MgSO4', (/ 0.100_r8, 72.269_r8, 2.2238_r8, 6.6037_r8, 0.0079004_r8, 3340.1_r8, 6.1304_r8 /))
438 : ! ---------------------------------------------
439 : ! ADD NEW SALT HERE (mass fraction, and others)
440 : ! ---------------------------------------------
441 1536 : End Subroutine SaltLibInitialization
442 :
443 : !--------------------------------------------------------------------------------
444 : !--------------------------------------------------------------------------------
445 32906670 : Function SpeciesIndex(SpeciesName)
446 : ! ==============================================
447 : ! This function is to look for the species index
448 : ! ==============================================
449 : Integer :: SpeciesIndex, i
450 : Character(Len=16) :: SpeciesName
451 :
452 493600050 : SpeciesIndex = -1 ! return -1 if species is not found
453 :
454 493600050 : Do i = 1, HowManyMolecules
455 493600050 : If (trim(SpeciesName) == trim(GasList(i)%CmpdName)) Then
456 32906670 : SpeciesIndex = i
457 32906670 : Exit
458 : Endif
459 : End Do
460 32906670 : End Function SpeciesIndex
461 :
462 : !--------------------------------------------------------------------------------
463 : !--------------------------------------------------------------------------------
464 15708747 : Real(r8) Function Flux_kg_m2_s(SpeciesName,Cgas_ppt,Cwater_nM,T_air_K,u10_m_s,P_atm,T_water_K,&
465 : Salinity_PartsPerThousand,switch_bubble)
466 : ! ===========================================================================
467 : ! This is the main module function. Input variables:
468 : ! ---------------------------------------------------------------------------
469 : ! - SpeciesName: name of species
470 : ! - Cgas_ppt: mixing ratio (parts per trillion) of trace gas of interest
471 : ! in the gas-phase (lowest modeling layer)
472 : ! - Cwater_nM: concentration of trace gas of interest in the surface ocean
473 : ! - T_air_K: temperature in the lowest modeling layer
474 : ! - u10_m_s: wind speed at 10 m above sea surface level
475 : ! - P_atm: air pressure in atm at sea surface level
476 : ! - T_water_K: sea surface temperature
477 : ! - Salinity_PartsPerThousand: surface ocean salinity
478 : ! - switch_bubble: bubble-mediated transfer switch
479 : ! ===========================================================================
480 : Character(16),intent(in) :: SpeciesName
481 : Real(r8),intent(in) :: Cgas_ppt, Cwater_nM, T_air_K, u10_m_s, P_atm, T_water_K, Salinity_PartsPerThousand
482 : Logical ,intent(in) :: switch_bubble
483 :
484 : Integer :: SpeciesID
485 : Real(r8) :: H_gas_over_liquid_dimless, kt_m_s
486 :
487 15708747 : SpeciesID = SpeciesIndex(SpeciesName)
488 : H_gas_over_liquid_dimless = 1.0_r8/(Henry_M_atm(SpeciesID,T_water_K,Salinity_PartsPerThousand)*&
489 15708747 : 0.082_r8*T_water_K)
490 15708747 : If (switch_bubble) then
491 : ! --------------------------------------------------------
492 : ! k_water parameterization with bubble-induced enhancement
493 : ! --------------------------------------------------------
494 : kt_m_s = (1.0_r8/k_water_m_s_bubble(SpeciesID, T_water_K, Salinity_PartsPerThousand, &
495 : u10_m_s, Cgas_ppt, P_atm, T_air_K) &
496 : + 1.0_r8/k_air_m_s(SpeciesID, u10_m_s, T_air_K, P_atm)&
497 0 : /H_gas_over_liquid_dimless)**(-1.0_r8)
498 : else
499 : ! ------------------------------------------------
500 : ! Original k_water parameterization, scaled to CO2
501 : ! ------------------------------------------------
502 : kt_m_s = (1.0_r8/k_water_m_s(SpeciesID, T_water_K, Salinity_PartsPerThousand, u10_m_s) &
503 15708747 : + 1.0_r8/k_air_m_s(SpeciesID, u10_m_s, T_air_K, P_atm)/H_gas_over_liquid_dimless)**(-1.0_r8)
504 : endif
505 : Flux_kg_m2_s = kt_m_s * (Cwater_nM*1E-9_r8*1000.0_r8 &
506 : - Cgas_ppt*1E-12_r8*(101325.0_r8*P_atm)/8.314_r8/T_air_K/H_gas_over_liquid_dimless) & ! g/m2/s
507 15708747 : * MolecularWeight(SpeciesIndex(SpeciesName)) / 1000.0_r8 ! convert to kg/m2/s
508 15708747 : End Function Flux_kg_m2_s
509 :
510 : !--------------------------------------------------------------------------------
511 : !--------------------------------------------------------------------------------
512 15708747 : Real(r8) Function k_air_m_s(SpeciesIndex, u10_m_s, T_air_K, P_atm)
513 : use shr_const_mod, only: vonKarman=>SHR_CONST_KARMAN
514 : ! =============================================================================
515 : ! Air-side transfer velocity. Slightly modified NOAA COARE (Fairall et al 2003;
516 : ! Feffery et al 2010), as recommended by Johnson Ocean Sci. 2010.
517 : ! Dynamic viscosity of air: Tsilingiris 2008
518 : ! =============================================================================
519 : Integer ,intent(in) :: SpeciesIndex
520 : Real(r8),intent(in) :: u10_m_s, T_air_K, P_atm
521 :
522 : Real(r8) :: ustar_m_s, DragCoeff
523 : Real(r8) :: DynamicViscosityAir_kg_m_s, DensityAir_kg_m3, DiffusivityInAir, SchmidtNumberInAir
524 :
525 : ! WSY: If local friction velocity is available from the model, might as well use that?
526 15708747 : ustar_m_s = u10_m_s * sqrt(6.1E-4_r8 + 6.3E-5_r8 * u10_m_s)
527 15708747 : DragCoeff = (ustar_m_s / u10_m_s)**2.0_r8
528 : DynamicViscosityAir_kg_m_s = 1.715747771E-5_r8 + 4.722402075E-8_r8 * (T_air_K-273.15_r8) &
529 : - 3.663027156E-10_r8 * ((T_air_K-273.15_r8)**2.0_r8) &
530 : + 1.873236686E-12_r8 * ((T_air_K-273.15_r8)**3.0_r8) &
531 15708747 : - 8.050218737E-14_r8 * ((T_air_K-273.15_r8)**4.0_r8)
532 : DensityAir_kg_m3 = 1.293393662_r8 - 5.538444326e-3_r8 * (T_air_K-273.15_r8) &
533 : + 3.860201577e-5_r8 * (T_air_K-273.15_r8)**2.0_r8 &
534 15708747 : - 5.2536065e-7_r8 * (T_air_K-273.15_r8)**3.0_r8
535 15708747 : DiffusivityInAir = DiffusivityInAir_cm2_s(SpeciesIndex, T_air_K, P_atm)
536 15708747 : SchmidtNumberInAir = DynamicViscosityAir_kg_m_s / DensityAir_kg_m3 / (DiffusivityInAir/10000.0_r8)
537 : k_air_m_s = 1E-3_r8 + ustar_m_s / (13.3_r8*(SchmidtNumberInAir**0.5_r8)+(DragCoeff**(-0.5_r8))-&
538 15708747 : 5.0_r8+log(SchmidtNumberInAir)/2.0_r8/vonKarman)
539 15708747 : End Function k_air_m_s
540 :
541 : !--------------------------------------------------------------------------------
542 : !--------------------------------------------------------------------------------
543 15708747 : Real(r8) Function k_water_m_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, u10_m_s)
544 : ! ================================================================================
545 : ! Water-side transfer velocity. Ref: Nightingale et al (2000). Salinity considered
546 : ! ================================================================================
547 : Integer ,intent(in) :: SpeciesIndex
548 : Real(r8),intent(in) :: T_water_K, Salinity_PartsPerThousand, u10_m_s
549 :
550 : Real(r8) :: DiffusivityInWater, SchmidtNumberInWater
551 : Real(r8) :: SchmidtNumberInWater_CO2ref
552 :
553 15708747 : SchmidtNumberInWater_CO2ref = 660.0_r8 ! this is the Schmidt number of CO2 at 20 degC in fresh water
554 15708747 : DiffusivityInWater = DiffusivityInWater_cm2_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand)
555 : SchmidtNumberInWater = DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand) / 1000.0_r8 &
556 15708747 : / DensityWater_kg_m3(T_water_K,Salinity_PartsPerThousand)/(DiffusivityInWater/10000.0_r8)
557 : k_water_m_s = ((0.222_r8*(u10_m_s**2.0_r8)+0.333_r8*u10_m_s)*&
558 15708747 : ((SchmidtNumberInWater/SchmidtNumberInWater_CO2ref)**(-0.5_r8)))/360000.0_r8
559 15708747 : End Function k_water_m_s
560 :
561 : !--------------------------------------------------------------------------------
562 : !--------------------------------------------------------------------------------
563 0 : Real(r8) Function k_water_m_s_bubble(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, u10_m_s, Cgas_ppt, P_atm, T_air_K)
564 : ! ==============================================================
565 : ! Water-side transfer velocity. Ref: Asher and Wanninkhof (1998).
566 : ! ==============================================================
567 : Integer, intent(in) :: SpeciesIndex
568 : Real(r8),intent(in) :: T_water_K, Salinity_PartsPerThousand, u10_m_s, Cgas_ppt, P_atm, T_air_K
569 :
570 : Real(r8) :: DiffusivityInWater, SchmidtNumberInWater
571 : Real(r8) :: FracCoverage_WhiteCaps, OstwaldSolubilityCoefficient
572 :
573 0 : DiffusivityInWater = DiffusivityInWater_cm2_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand)
574 : SchmidtNumberInWater = DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand) / 1000.0_r8 &
575 0 : / DensityWater_kg_m3(T_water_K,Salinity_PartsPerThousand)/(DiffusivityInWater/10000.0_r8)
576 0 : FracCoverage_WhiteCaps = 2.56e-6_r8 * (u10_m_s - 1.77_r8)**3.0_r8
577 0 : OstwaldSolubilityCoefficient = Henry_M_atm(SpeciesIndex,T_water_K,Salinity_PartsPerThousand) ! just Henry's law (M/atm)
578 0 : OstwaldSolubilityCoefficient = OstwaldSolubilityCoefficient * (Cgas_ppt*1.0E-12_r8*P_atm) ! mol / L
579 0 : OstwaldSolubilityCoefficient = OstwaldSolubilityCoefficient * 0.082_r8 * T_air_K / P_atm ! L / L
580 : k_water_m_s_bubble = ((47.0_r8*u10_m_s + FracCoverage_WhiteCaps*(115200.0_r8 - 47.0_r8* u10_m_s)) &
581 : *(SchmidtNumberInWater**(-0.5_r8)) &
582 : + FracCoverage_WhiteCaps * (-37.0_r8/OstwaldSolubilityCoefficient &
583 : + 6120.0_r8*(OstwaldSolubilityCoefficient**(-0.37_r8)) *(SchmidtNumberInWater**(-0.18_r8)))) &
584 0 : * 2.8e-6_r8
585 :
586 0 : End Function k_water_m_s_bubble
587 :
588 : !--------------------------------------------------------------------------------
589 : !--------------------------------------------------------------------------------
590 15708747 : Real(r8) Function DiffusivityInAir_cm2_s(SpeciesIndex, T_air_K, P_atm)
591 : ! ============================
592 : ! Ref: Johnson Ocean Sci. 2010
593 : ! ============================
594 : Integer ,intent(in) :: SpeciesIndex
595 : Real(r8),intent(in) :: T_air_K, P_atm
596 :
597 : Real(r8), parameter :: MW_air = 28.97_r8 ! molecular weight for air
598 : Real(r8), parameter :: Va = 20.1_r8 ! molar volume for air
599 : Real(r8) :: Vb, MW_species
600 :
601 15708747 : Vb = LiquidMolarVolume_cm3_mol(SpeciesIndex)
602 15708747 : MW_species = MolecularWeight(SpeciesIndex)
603 : DiffusivityInAir_cm2_s = 0.001_r8 * (T_air_K**1.75_r8) & ! oh f* me
604 : * (((MW_air + MW_species)/(MW_air*MW_species))**0.5_r8) &
605 15708747 : / ((P_atm*(Va**(1.0_r8/3.0_r8)+Vb**(1.0_r8/3.0_r8)))**2.0_r8)
606 :
607 15708747 : End Function DiffusivityInAir_cm2_s
608 :
609 : !--------------------------------------------------------------------------------
610 : !--------------------------------------------------------------------------------
611 15708747 : Real(r8) Function DiffusivityInWater_cm2_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand)
612 : ! =================================================
613 : ! Ref: Johnson Ocean Sci. 2010. Salinity considered
614 : ! =================================================
615 : Integer, intent(in) :: SpeciesIndex
616 : Real(r8),intent(in) :: T_water_K, Salinity_PartsPerThousand
617 :
618 : Real(r8), parameter :: AssociationFactor = 2.6_r8 ! ... for water
619 : Real(r8) :: DynamicViscosityWater, Vb, MW_species
620 :
621 15708747 : Vb = LiquidMolarVolume_cm3_mol(SpeciesIndex)
622 15708747 : MW_species = MolecularWeight(SpeciesIndex)
623 :
624 15708747 : DynamicViscosityWater = DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand)
625 : ! -------------------------------------------------
626 : ! Wilke and Chang 1955: this seems to be a bit high
627 : ! -------------------------------------------------
628 : ! DiffusivityInWater_cm2_s = 7.4E-8_r8 * T_water_K * sqrt(AssociationFactor*MW_species) /
629 : ! DynamicViscosityWater / (Vb**0.6_r8)
630 : ! ----------------------
631 : ! Hayduk and Minhas 1982
632 : ! ----------------------
633 : DiffusivityInWater_cm2_s = 1.25E-8_r8 * (T_water_K**1.52_r8) * (DynamicViscosityWater**(9.58_r8/Vb - 1.12_r8)) &
634 15708747 : * (Vb**(-0.19_r8)-0.292_r8)
635 :
636 15708747 : End Function DiffusivityInWater_cm2_s
637 :
638 : !--------------------------------------------------------------------------------
639 : !--------------------------------------------------------------------------------
640 31417494 : Real(r8) Function DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand)
641 : ! =================================================
642 : ! Ref: Johnson Ocean Sci. 2010. Salinity considered
643 : ! =================================================
644 : Real(r8),intent(in) :: T_water_K, Salinity_PartsPerThousand
645 :
646 : Real(r8) :: MassFrac_water, DynamicViscosityPureWater_g_m_s, SaltViscosity, sum_w_ln_SaltViscosity
647 : Integer :: n
648 :
649 31417494 : sum_w_ln_SaltViscosity = 0.0_r8
650 31417494 : MassFrac_water = 1.0_r8 - Salinity_PartsPerThousand / 1000.0_r8
651 : DynamicViscosityPureWater_g_m_s = ((T_water_K-273.15_r8)+246.0_r8) &
652 31417494 : / (0.05594_r8*(T_water_K-273.15_r8)**2.0_r8+5.2842_r8*(T_water_K-273.15_r8)+137.37_r8)
653 :
654 31417494 : If (Salinity_PartsPerThousand == 0.0_r8) Then ! pure water
655 : DynamicViscosityWater_g_m_s = DynamicViscosityPureWater_g_m_s
656 : Else ! salty water
657 188504964 : Do n = 1, HowManySalts
658 157087470 : SaltViscosity = exp((SaltList(n)%SaltProperties(2) * &
659 : (Salinity_PartsPerThousand/1000.0_r8)**SaltList(n)%SaltProperties(3) &
660 : + SaltList(n)%SaltProperties(4)) &
661 : / (SaltList(n)%SaltProperties(5)*(T_water_K-273.15_r8) + 1.0_r8)) &
662 : / (SaltList(n)%SaltProperties(6) * (Salinity_PartsPerThousand / &
663 157087470 : 1000.0_r8)**SaltList(n)%SaltProperties(7) + 1.0_r8)
664 : sum_w_ln_SaltViscosity = sum_w_ln_SaltViscosity + (Salinity_PartsPerThousand/1000.0_r8) &
665 188504964 : * SaltList(n)%SaltProperties(1) * log(SaltViscosity)
666 : End Do
667 : DynamicViscosityWater_g_m_s = exp(MassFrac_water &
668 31417494 : * log(DynamicViscosityPureWater_g_m_s) + sum_w_ln_SaltViscosity)
669 : Endif
670 :
671 31417494 : End Function DynamicViscosityWater_g_m_s
672 :
673 : !--------------------------------------------------------------------------------
674 : !--------------------------------------------------------------------------------
675 15708747 : Real(r8) Function DensityWater_kg_m3(T_water_K, Salinity_PartsPerThousand)
676 : ! ====================================================
677 : ! Ref: Millero and Poisson (1981). Salinity considered
678 : ! ====================================================
679 : Real(r8), intent(in) :: T_water_K, Salinity_PartsPerThousand
680 :
681 : Real(r8) :: DensityPureWater_kg_m3, FactorA, FactorB, FactorC
682 :
683 : DensityPureWater_kg_m3 = 999.842594_r8 + 0.06793952_r8*(T_water_K-273.15_r8) &
684 : - 0.00909529_r8*((T_water_K-273.15_r8)**2.0_r8) &
685 : + 0.0001001685_r8*((T_water_K-273.15_r8)**3.0_r8) &
686 : - 0.000001120083_r8*((T_water_K-273.15_r8)**4.0_r8) &
687 15708747 : + 0.000000006536332_r8*((T_water_K-273.15_r8)**5.0_r8)
688 : FactorA = 0.824493_r8 - 0.0040899_r8*(T_water_K-273.15_r8) + 0.000076438_r8*((T_water_K-273.15_r8)**2.0_r8) &
689 15708747 : - 0.00000082467_r8*((T_water_K-273.15_r8)**3.0_r8) + 0.0000000053875_r8*((T_water_K-273.15_r8)**4.0_r8)
690 15708747 : FactorB = - 0.00572466_r8 + 0.00010277_r8*(T_water_K-273.15_r8) - 0.0000016546_r8*((T_water_K-273.15_r8)**2.0_r8)
691 15708747 : FactorC = 0.00048314_r8
692 : DensityWater_kg_m3 = DensityPureWater_kg_m3 + FactorA*Salinity_PartsPerThousand &
693 15708747 : + FactorB*(Salinity_PartsPerThousand**(2.0_r8/3.0_r8)) + FactorC*Salinity_PartsPerThousand
694 :
695 15708747 : End Function DensityWater_kg_m3
696 :
697 : !--------------------------------------------------------------------------------
698 : !--------------------------------------------------------------------------------
699 15708747 : Real(r8) Function Henry_M_atm(SpeciesIndex, T_water_K, Salinity_PartsPerThousand)
700 : ! =========================================================================================
701 : ! Ref: Sander compilation 2015. Salt-in or salt-out estimated based on Setschenow constants
702 : ! =========================================================================================
703 : Integer, intent(in) :: SpeciesIndex
704 : Real(r8), intent(in) :: T_water_K, Salinity_PartsPerThousand
705 :
706 : Real(r8) :: Heff_M_atm_PureWater, Setschenow, Heff_M_atm_SaltyWater
707 :
708 15708747 : Heff_M_atm_PureWater = GasList(SpeciesIndex)%CmpdProperties(15) * &
709 15708747 : exp(GasList(SpeciesIndex)%CmpdProperties(16) * (1.0_r8/T_water_K - 1.0_r8/298.0_r8))
710 :
711 15708747 : If (Salinity_PartsPerThousand==0.0_r8) Then
712 : Henry_M_atm = Heff_M_atm_PureWater
713 : Else
714 : Setschenow = log(LiquidMolarVolume_cm3_mol(SpeciesIndex)) * &
715 : (7.33532E-4_r8 + 3.39615E-5_r8 * log(Heff_M_atm_PureWater) &
716 : - 2.40888E-6_r8 * ((log(Heff_M_atm_PureWater))**2.0_r8) &
717 15708747 : + 1.57114E-7_r8 * ((log(Heff_M_atm_PureWater))**3.0_r8))
718 15708747 : Heff_M_atm_SaltyWater = Heff_M_atm_PureWater * 10.0_r8**(Setschenow*Salinity_PartsPerThousand)
719 15708747 : Henry_M_atm = Heff_M_atm_SaltyWater
720 : Endif
721 :
722 15708747 : End Function Henry_M_atm
723 :
724 : !--------------------------------------------------------------------------------
725 : !--------------------------------------------------------------------------------
726 48615417 : Function MolecularWeight(SpeciesIndex)
727 : Real(r8) :: MolecularWeight
728 : Integer :: SpeciesIndex
729 48615417 : MolecularWeight = GasList(SpeciesIndex)%CmpdProperties(1)
730 48615417 : End Function MolecularWeight
731 :
732 : !--------------------------------------------------------------------------------
733 : !--------------------------------------------------------------------------------
734 47126241 : Function LiquidMolarVolume_cm3_mol(SpeciesIndex)
735 : ! ===========================================================================
736 : ! If no measurements available, i.e. GasList(SpeciesIndex)%CmpdProperties(14)
737 : ! is zero, Schroeder's group contribution method is used
738 : ! ===========================================================================
739 : Real(r8) :: LiquidMolarVolume_cm3_mol
740 : Integer :: SpeciesIndex
741 :
742 47126241 : If (GasList(SpeciesIndex)%CmpdProperties(14)/=0.0_r8) Then
743 : LiquidMolarVolume_cm3_mol = GasList(SpeciesIndex)%CmpdProperties(14)
744 : Else
745 47126241 : LiquidMolarVolume_cm3_mol = 7.0_r8*GasList(SpeciesIndex)%CmpdProperties(2) ! C
746 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol + 7.0_r8*GasList(SpeciesIndex)%CmpdProperties(3) ! H
747 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol + 7.0_r8*GasList(SpeciesIndex)%CmpdProperties(4) ! N
748 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol + 7.0_r8*GasList(SpeciesIndex)%CmpdProperties(5) ! O
749 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol + 21.0_r8*GasList(SpeciesIndex)%CmpdProperties(6) ! S
750 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol + 31.5_r8*GasList(SpeciesIndex)%CmpdProperties(7) ! Br
751 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol + 24.5_r8*GasList(SpeciesIndex)%CmpdProperties(8) ! Cl
752 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol + 10.5_r8*GasList(SpeciesIndex)%CmpdProperties(9) ! F
753 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol + 38.5_r8*GasList(SpeciesIndex)%CmpdProperties(10) ! I
754 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol - 7.0_r8*GasList(SpeciesIndex)%CmpdProperties(13) ! Ring
755 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol + 7.0_r8*GasList(SpeciesIndex)%CmpdProperties(11) ! Double bond
756 47126241 : LiquidMolarVolume_cm3_mol = LiquidMolarVolume_cm3_mol + 14.0_r8*GasList(SpeciesIndex)%CmpdProperties(12) ! Triple bond
757 : Endif
758 :
759 47126241 : End Function LiquidMolarVolume_cm3_mol
760 :
761 : !--------------------------------------------------------------------------------
762 : !--------------------------------------------------------------------------------
763 1536 : subroutine cseawater_ini()
764 :
765 : use mo_chem_utls, only : get_spc_ndx
766 : use tracer_data, only : trcdata_init
767 : use cam_pio_utils, only : cam_pio_openfile
768 : use pio, only : pio_inquire, pio_nowrite, pio_closefile, pio_inq_varndims
769 : use pio, only : pio_inq_varname, file_desc_t, pio_get_att, PIO_NOERR, PIO_GLOBAL
770 : use pio, only : pio_seterrorhandling, PIO_BCAST_ERROR
771 : use string_utils, only : GLC
772 : use phys_control, only : phys_getopts
773 :
774 : integer :: i, j, l, m, n, nn, astat, vid, ierr, nvars, isec
775 : integer :: indx(gas_pcnst)
776 : type(file_desc_t) :: ncid
777 : character(len=16) :: csw_species(gas_pcnst)
778 : character(len=256) :: csw_filenam(gas_pcnst)
779 : integer :: csw_indexes(gas_pcnst)
780 : real(r8) :: csw_scalefactor(gas_pcnst)
781 :
782 : character(len=80) :: file_interp_type = ' '
783 :
784 : character(len=16) :: spc_name
785 : character(len=256) :: filename
786 : character(len=256) :: tmp_string = ' '
787 : character(len=32) :: xchr = ' '
788 : real(r8) :: xdbl
789 :
790 : character(len=32) :: varname
791 : character(len=256) :: locfn
792 1536 : integer, allocatable :: vndims(:)
793 : character(len=1), parameter :: filelist = ''
794 : character(len=1), parameter :: datapath = ''
795 : logical, parameter :: rmv_file = .false.
796 :
797 : character(len=*), parameter :: subname = 'cseawater_ini'
798 :
799 : logical :: history_chemistry
800 1536 : call phys_getopts(history_chemistry_out=history_chemistry)
801 :
802 : ! ========================================================
803 : ! Read sea water concentration specifier from the namelist
804 : ! ========================================================
805 :
806 1536 : nn = 0
807 1536 : indx(:) = 0
808 :
809 1536 : switch_bubble = bubble_mediated_transfer
810 :
811 3072 : count_Csw: do n=1,gas_pcnst
812 3072 : if ( len_trim(csw_specifier(n) ) == 0 ) then
813 : exit count_Csw
814 : endif
815 :
816 1536 : i = scan(csw_specifier(n),'->')
817 1536 : spc_name = trim(adjustl(csw_specifier(n)(:i-1)))
818 :
819 : ! --------------------------
820 : ! get the scaling factor
821 : ! --------------------------
822 1536 : tmp_string = adjustl(csw_specifier(n)(i+2:))
823 1536 : j = scan( tmp_string, '*' )
824 1536 : if (j>0) then
825 0 : xchr = tmp_string(1:j-1) ! get the multipler (left of the '*')
826 0 : read( xchr, * ) xdbl ! convert the string to a real
827 0 : tmp_string = adjustl(tmp_string(j+1:)) ! get the filepath name (right of the '*')
828 : else
829 1536 : xdbl = 1._r8
830 : endif
831 1536 : filename = trim(tmp_string)
832 :
833 1536 : m = get_spc_ndx(spc_name)
834 :
835 1536 : if ( m<1 ) then
836 0 : if (masterproc) write(iulog,*) 'cseawater_inti: spc_name ',spc_name,' is not included in the simulation'
837 0 : call endrun('cseawater_inti: invalid seawater concentration specification')
838 : endif
839 :
840 1536 : nn = nn+1
841 1536 : csw_species(nn) = spc_name
842 1536 : csw_filenam(nn) = filename
843 1536 : csw_indexes(nn) = m
844 1536 : csw_scalefactor(nn) = xdbl
845 :
846 3072 : indx(n)=n
847 :
848 : enddo count_Csw
849 :
850 1536 : n_Csw_files = nn
851 :
852 1536 : if (masterproc) write(iulog,*) subname//': n_Csw_files = ',n_Csw_files
853 :
854 12288 : allocate( Csw_nM(n_Csw_files), stat=astat )
855 1536 : if( astat/= 0 ) then
856 0 : write(iulog,*) subname//': failed to allocate Csw_nM array; error = ',astat
857 0 : call endrun(subname//': failed to allocate Csw_nM array')
858 : end if
859 :
860 : ! -------------------------------------------
861 : ! Setup the seawater concentration type array
862 : ! -------------------------------------------
863 3072 : do m=1,n_Csw_files
864 1536 : Csw_nM(m)%spc_ndx = csw_indexes(indx(m))
865 1536 : Csw_nM(m)%units = 'nM'
866 1536 : Csw_nM(m)%species = csw_species(indx(m))
867 1536 : Csw_nM(m)%filename = csw_filenam(indx(m))
868 3072 : Csw_nM(m)%scalefactor = csw_scalefactor(indx(m))
869 : enddo
870 :
871 : !-----------------------------------------------
872 : ! read emis files to determine number of sectors
873 : !-----------------------------------------------
874 3072 : spc_loop: do m = 1, n_Csw_files
875 :
876 1536 : Csw_nM(m)%nsectors = 0
877 :
878 1536 : call getfil (Csw_nM(m)%filename, locfn, 0)
879 1536 : call cam_pio_openfile ( ncid, trim(locfn), PIO_NOWRITE)
880 :
881 1536 : call pio_seterrorhandling(ncid, PIO_BCAST_ERROR)
882 :
883 1536 : ierr = pio_inquire(ncid, nvariables=nvars)
884 1536 : if (ierr /= PIO_NOERR) then
885 0 : call endrun(subname//': pio_inquire nvars FAILED')
886 : endif
887 :
888 4608 : allocate(vndims(nvars))
889 :
890 9216 : do vid = 1,nvars
891 :
892 7680 : ierr = pio_inq_varndims (ncid, vid, vndims(vid))
893 7680 : if (ierr /= PIO_NOERR) then
894 0 : call endrun(subname//': pio_inq_varndims FAILED')
895 : endif
896 :
897 7680 : if( vndims(vid) < 3 ) then
898 : cycle
899 1536 : elseif( vndims(vid) > 3 ) then
900 0 : ierr = pio_inq_varname (ncid, vid, varname)
901 0 : if (ierr /= PIO_NOERR) then
902 0 : call endrun(subname//': pio_inq_varname varname FAILED')
903 : endif
904 0 : write(iulog,*) subname//': Skipping variable ', trim(varname),', ndims = ',vndims(vid), &
905 0 : ', species=',trim(Csw_nM(m)%species)
906 0 : cycle
907 : end if
908 :
909 9216 : Csw_nM(m)%nsectors = Csw_nM(m)%nsectors+1
910 :
911 : enddo
912 :
913 4608 : allocate( Csw_nM(m)%sectors(Csw_nM(m)%nsectors), stat=astat )
914 1536 : if( astat/= 0 ) then
915 0 : write(iulog,*) subname//': failed to allocate Csw_nM(m)%sectors array; error = ',astat
916 0 : call endrun(subname//': failed to allocate Csw_nM(m)%sectors array')
917 : end if
918 :
919 1536 : isec = 1
920 :
921 9216 : do vid = 1,nvars
922 9216 : if( vndims(vid) == 3 ) then
923 1536 : ierr = pio_inq_varname(ncid, vid, Csw_nM(m)%sectors(isec))
924 1536 : if (ierr /= PIO_NOERR) then
925 0 : call endrun(subname//': pio_inq_varname Csw_nM(m)%sectors(isec) FAILED')
926 : endif
927 1536 : isec = isec+1
928 : endif
929 :
930 : enddo
931 1536 : deallocate(vndims)
932 :
933 : ! Global attribute 'input_method' overrides the srf_emis_type namelist setting on
934 : ! a file-by-file basis. If the emis file does not contain the 'input_method'
935 : ! attribute then the srf_emis_type namelist setting is used.
936 1536 : ierr = pio_get_att(ncid, PIO_GLOBAL, 'input_method', file_interp_type)
937 1536 : if ( ierr == PIO_NOERR) then
938 0 : l = GLC(file_interp_type)
939 0 : csw_time_type(1:l) = file_interp_type(1:l)
940 0 : csw_time_type(l+1:) = ' '
941 : endif
942 :
943 1536 : call pio_closefile (ncid)
944 :
945 4608 : allocate(Csw_nM(m)%file%in_pbuf(size(Csw_nM(m)%sectors)))
946 3072 : Csw_nM(m)%file%in_pbuf(:) = .false.
947 :
948 0 : call trcdata_init( Csw_nM(m)%sectors, &
949 : Csw_nM(m)%filename, filelist, datapath, &
950 : Csw_nM(m)%fields, &
951 : Csw_nM(m)%file, &
952 3072 : rmv_file, csw_cycle_yr, 0, 0, trim(csw_time_type) )
953 :
954 : enddo spc_loop
955 :
956 : ! ===================================
957 : ! Output stuff to the history files
958 : ! ===================================
959 3072 : do m = 1, n_Csw_files
960 0 : call addfld('OCN_FLUX_' // trim(Csw_nM(m)%species), horiz_only, 'A', 'kg/m2/s', &
961 1536 : 'ocean flux ' // trim(Csw_nM(m)%species) )
962 0 : call addfld('Csw_' // trim(Csw_nM(m)%species), horiz_only, 'A', 'nanomole per liter (nM)', &
963 1536 : 'seeawater concentration ' // trim(Csw_nM(m)%species) )
964 3072 : if (history_chemistry) then
965 1536 : call add_default('OCN_FLUX_' // trim(Csw_nM(m)%species), 1, ' ')
966 : end if
967 : end do
968 :
969 3072 : end subroutine cseawater_ini
970 :
971 : !--------------------------------------------------------------------------------
972 : ! returns TRUE if species has ocean emissions
973 : !--------------------------------------------------------------------------------
974 46080 : pure logical function ocean_emis_species(name)
975 : character(len=*), intent(in) :: name
976 :
977 : integer :: m
978 :
979 46080 : ocean_emis_species = .false.
980 :
981 90624 : spc_loop: do m = 1, n_Csw_files
982 90624 : if (trim(name) == trim(Csw_nM(m)%species)) then
983 1536 : ocean_emis_species = .true.
984 1536 : exit spc_loop
985 : end if
986 : end do spc_loop
987 :
988 1536 : end function ocean_emis_species
989 :
990 0 : end module ocean_emis
|
