LCOV - code coverage report
Current view: | top level - chemistry/mozart - ocean_emis.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 280 | 332 | 84.3 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 19 | 24 | 79.2 %

              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

Generated by: LCOV version 1.14