diff --git a/Externals.cfg b/Externals.cfg index 13ddc8f76e..dc66b3c485 100644 --- a/Externals.cfg +++ b/Externals.cfg @@ -21,7 +21,7 @@ externals = Externals.cfg required = True [cmeps] -tag = cmeps0.14.18 +tag = cmeps0.14.24 protocol = git repo_url = https://github.com/ESCOMP/CMEPS.git local_path = components/cmeps diff --git a/bld/build-namelist b/bld/build-namelist index 20cdfebe6c..c63c36cbbc 100755 --- a/bld/build-namelist +++ b/bld/build-namelist @@ -3475,14 +3475,16 @@ if ( length($nl->get_value('soil_erod_file'))>0 ) { else { if ($chem =~ /trop_strat/ or $chem =~ /waccm_ma/ or $chem =~ /waccm_tsmlt/ or $chem =~ /trop_mozart/) { add_default($nl, 'dust_emis_fact', 'ver'=>'chem'); - # set scaling of lightning NOx production - add_default($nl, 'lght_no_prd_factor' ); } else { add_default($nl, 'dust_emis_fact'); } } } +if (chem_has_species($cfg, 'NO')) { + # set scaling of lightning NOx production + add_default($nl, 'lght_no_prd_factor' ); +} # Seasalt emissions tuning factor if ($chem =~ /_mam(\d)/) { @@ -4111,6 +4113,21 @@ add_default($nl, 'cam_snapshot_before_num'); add_default($nl, 'cam_snapshot_after_num'); check_snapshot_settings(); +if ($opts{'cmeps'}) { + # advertise the nature of ozone data passed to surface models + if ($rad_prog_ozone) { + add_default($nl, 'atm_ozone_frequency', 'val'=>'subdaily'); + } else { + add_default($nl, 'atm_ozone_frequency', 'val'=>'multiday_average'); + } + # for lightning flash freq to CTSM + if ($simple_phys or $aqua_mode) { + add_default($nl, 'atm_provides_lightning', 'val'=>'.false.'); + } else { + add_default($nl, 'atm_provides_lightning', 'val'=>'.true.'); + } +} + #----------------------------------------------------------------------------------------------- # Write output files @@ -4127,16 +4144,8 @@ my %nl_group = (); foreach my $name (@nl_groups) { $nl_group{$name} = ''; } # Dry deposition, MEGAN VOC emis and ozone namelists -@comp_groups = qw(drydep_inparm megan_emis_nl fire_emis_nl carma_inparm ndep_inparm ozone_coupling_nl); +@comp_groups = qw(drydep_inparm megan_emis_nl fire_emis_nl carma_inparm ndep_inparm ozone_coupling_nl lightning_coupling_nl); -# nature of ozone data passed to surface models -- only if cmeps (nuopc) coupling is used -if ($opts{'cmeps'}) { - if ($rad_prog_ozone) { - add_default($nl, 'atm_ozone_frequency', 'val'=>'subdaily'); - } else { - add_default($nl, 'atm_ozone_frequency', 'val'=>'multiday_average'); - } -} $outfile = "$opts{'dir'}/drv_flds_in"; $nl->write($outfile, 'groups'=>\@comp_groups); if ($print>=1) { diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index 19fe6883b1..0457b3e499 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -6197,14 +6197,8 @@ List of species that are constrained in the stratosphere. Default: set by build-namelist. - -Full pathname of dataset for land mask applied to the lighting NOx production -Default: set by build-namelist. - - + group="lightning_nl" valid_values="" > Multiplication factor applied to the lighting NOx production Default: 1.0. @@ -7312,6 +7306,12 @@ coarser temporal resolution. Default: set by build-namelist. + +If TRUE atmosphere model will provide prognosed lightning flash frequency. +Default: FALSE + + 0 .and. cldbot_ndx>0 + + if (.not.calc_lightning) return + + calc_nox_prod = lght_no_prd_factor>0._r8 + + if (calc_nox_prod) then + + if (masterproc) write(iulog,*) prefix,'lightning no production scaling factor = ',factor + + !---------------------------------------------------------------------- + ! ... vdist(kk,itype) = % of lightning nox between (kk-1) and (kk) + ! km for profile itype + !---------------------------------------------------------------------- + allocate(vdist(16,3),stat=astat) + if( astat /= 0 ) then + write(iulog,*) prefix,'failed to allocate vdist; error = ',astat + call endrun(prefix//'failed to allocate vdist') + end if + vdist(:,1) = (/ 3.0_r8, 3.0_r8, 3.0_r8, 3.0_r8, 3.4_r8, 3.5_r8, 3.6_r8, 4.0_r8, & ! midlat cont + 5.0_r8, 7.0_r8, 9.0_r8, 14.0_r8, 16.0_r8, 14.0_r8, 8.0_r8, 0.5_r8 /) + vdist(:,2) = (/ 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 6.1_r8, & ! trop marine + 17.0_r8, 15.4_r8, 14.5_r8, 13.0_r8, 12.5_r8, 2.8_r8, 0.9_r8, 0.3_r8 /) + vdist(:,3) = (/ 2.0_r8, 2.0_r8, 2.0_r8, 1.5_r8, 1.5_r8, 1.5_r8, 3.0_r8, 5.8_r8, & ! trop cont + 7.6_r8, 9.6_r8, 11.0_r8, 14.0_r8, 14.0_r8, 14.0_r8, 8.2_r8, 2.3_r8 /) + + allocate( prod_no(pcols,pver,begchunk:endchunk),stat=astat ) + if( astat /= 0 ) then + write(iulog,*) prefix, 'failed to allocate prod_no; error = ',astat + call endrun(prefix//'failed to allocate prod_no') + end if + geo_factor = ngcols_p/(4._r8*pi) + + call addfld( 'LNO_COL_PROD', horiz_only, 'I', 'Tg N yr-1', 'lightning column NO source' ) + call addfld( 'LNO_PROD', (/ 'lev' /), 'I', 'molecules/cm3/s', 'lightning insitu NO source' ) + call addfld( 'FLASHENGY', horiz_only, 'I', 'J', 'lightning flash energy' ) ! flash energy + + call phys_getopts( history_cesm_forcing_out = history_cesm_forcing ) + if ( history_cesm_forcing ) then + call add_default('LNO_COL_PROD',1,' ') + endif + + if (is_first_step()) then + call pbuf_set_field(pbuf2d, flsh_frq_ndx, 0.0_r8) + endif - call phys_getopts( history_cesm_forcing_out = history_cesm_forcing ) - - no_ndx = get_spc_ndx('NO') - xno_ndx = get_spc_ndx('XNO') - - has_no_lightning_prod = no_ndx>0 .or. xno_ndx>0 - if (.not.has_no_lightning_prod) return - - - if( lght_no_prd_factor /= 1._r8 ) then - factor = factor*lght_no_prd_factor - end if - - - if (masterproc) write(iulog,*) 'lght_inti: lightning no production scaling factor = ',factor - - !---------------------------------------------------------------------- - ! ... vdist(kk,itype) = % of lightning nox between (kk-1) and (kk) - ! km for profile itype - !---------------------------------------------------------------------- - vdist(:,1) = (/ 3.0_r8, 3.0_r8, 3.0_r8, 3.0_r8, 3.4_r8, 3.5_r8, 3.6_r8, 4.0_r8, & ! midlat cont - 5.0_r8, 7.0_r8, 9.0_r8, 14.0_r8, 16.0_r8, 14.0_r8, 8.0_r8, 0.5_r8 /) - vdist(:,2) = (/ 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 2.5_r8, 6.1_r8, & ! trop marine - 17.0_r8, 15.4_r8, 14.5_r8, 13.0_r8, 12.5_r8, 2.8_r8, 0.9_r8, 0.3_r8 /) - vdist(:,3) = (/ 2.0_r8, 2.0_r8, 2.0_r8, 1.5_r8, 1.5_r8, 1.5_r8, 3.0_r8, 5.8_r8, & ! trop cont - 7.6_r8, 9.6_r8, 11.0_r8, 14.0_r8, 14.0_r8, 14.0_r8, 8.2_r8, 2.3_r8 /) - - allocate( prod_no(pcols,pver,begchunk:endchunk),stat=astat ) - if( astat /= 0 ) then - write(iulog,*) 'lght_inti: failed to allocate prod_no; error = ',astat - call endrun - end if - allocate( flash_freq(pcols,begchunk:endchunk),stat=astat ) - if( astat /= 0 ) then - write(iulog,*) 'lght_inti: failed to allocate flash_freq; error = ',astat - call endrun - end if - allocate( glob_prod_no_col(pcols,begchunk:endchunk),stat=astat ) - if( astat /= 0 ) then - write(iulog,*) 'lght_inti: failed to allocate glob_prod_no_col; error = ',astat - call endrun - end if - prod_no(:,:,:) = 0._r8 - flash_freq(:,:) = 0._r8 - geo_factor = ngcols_p/(4._r8*pi) - - - call addfld( 'LNO_COL_PROD', horiz_only, 'I', 'TG N/YR', 'lighting column NO source' ) - call addfld( 'LNO_PROD', (/ 'lev' /), 'I', '/cm3/s', 'lighting insitu NO source' ) - call addfld( 'FLASHFRQ', horiz_only, 'I', '1/MIN', 'lighting flash rate' ) ! flash frequency in grid box per minute (PPP) - call addfld( 'FLASHENGY', horiz_only, 'I', ' ', 'lighting flash rate' ) ! flash frequency in grid box per minute (PPP) - call addfld( 'CLDHGT', horiz_only, 'I', 'KM', 'cloud top height' ) ! cloud top height - call addfld( 'DCHGZONE', horiz_only, 'I', 'KM', 'depth of discharge zone' ) ! depth of discharge zone - call addfld( 'CGIC', horiz_only, 'I', 'RATIO', 'ratio of cloud-ground/intracloud discharges' ) ! ratio of cloud-ground/intracloud discharges - - if ( history_cesm_forcing ) then - call add_default('LNO_COL_PROD',1,' ') endif - end subroutine lightning_inti + call addfld( 'FLASHFRQ', horiz_only, 'I', 'min-1', 'lightning flash rate' ) ! flash frequency in grid box per minute (PPP) + call addfld( 'CLDHGT', horiz_only, 'I', 'km', 'cloud top height' ) ! cloud top height + call addfld( 'DCHGZONE', horiz_only, 'I', 'km', 'depth of discharge zone' ) ! depth of discharge zone + call addfld( 'CGIC', horiz_only, 'I', '1', 'ratio of cloud-ground/intracloud discharges' ) ! ratio of cloud-ground/intracloud discharges + call addfld( 'LGHTNG_CLD2GRND', horiz_only, 'I', 'min-1', 'clound-to-ground lightning flash rate') ! clound to ground flash frequency + end subroutine lightning_init + + !------------------------------------------------------------------------- + !------------------------------------------------------------------------- subroutine lightning_no_prod( state, pbuf2d, cam_in ) !---------------------------------------------------------------------- ! ... set no production from lightning !---------------------------------------------------------------------- use physics_types, only : physics_state - - use physics_buffer, only : pbuf_get_index, physics_buffer_desc, pbuf_get_field, pbuf_get_chunk use physconst, only : rga use phys_grid, only : get_rlat_all_p, get_wght_all_p use cam_history, only : outfld use camsrfexch, only : cam_in_t use shr_reprosum_mod, only : shr_reprosum_calc - use mo_constants, only : rearth, d2r - implicit none + use mo_constants, only : rearth, d2r !---------------------------------------------------------------------- ! ... dummy args !---------------------------------------------------------------------- type(physics_state), intent(in) :: state(begchunk:endchunk) ! physics state - type(physics_buffer_desc), pointer :: pbuf2d(:,:) type(cam_in_t), intent(in) :: cam_in(begchunk:endchunk) ! physics state !---------------------------------------------------------------------- ! ... local variables !---------------------------------------------------------------------- - real(r8), parameter :: land = 1._r8 - real(r8), parameter :: secpyr = 365._r8 * 8.64e4_r8 + real(r8), parameter :: land = 1._r8 + real(r8), parameter :: secpyr = 365._r8 * 8.64e4_r8 - integer :: i, c integer :: cldtind ! level index for cloud top integer :: cldbind ! level index for cloud base > 273k integer :: k, kk, zlow_ind, zhigh_ind, itype @@ -162,16 +227,20 @@ subroutine lightning_no_prod( state, pbuf2d, cam_in ) real(r8) :: flash_energy(pcols,begchunk:endchunk) ! energy of flashes per second real(r8) :: prod_no_col(pcols,begchunk:endchunk) ! global no production rate for diagnostics real(r8) :: wrk, wrk1, wrk2(1) + integer :: icol ! column index integer :: ncol ! columns per chunk - integer :: lchnk ! columns per chunk + integer :: lchnk ! chunk index real(r8),pointer :: cldtop(:) ! cloud top level index real(r8),pointer :: cldbot(:) ! cloud bottom level index real(r8) :: zmid(pcols,pver) ! geopot height above surface at midpoints (m) real(r8) :: zint(pcols,pver+1,begchunk:endchunk) ! geopot height above surface at interfaces (m) real(r8) :: zsurf(pcols) ! geopot height above surface at interfaces (m) - real(r8) :: rlats(pcols,begchunk:endchunk) ! column latitudes in chunks + real(r8) :: rlats(pcols) ! column latitudes in chunks real(r8) :: wght(pcols) + real(r8) :: glob_prod_no_col(pcols,begchunk:endchunk) + real(r8) :: flash_freq(pcols,begchunk:endchunk) + !---------------------------------------------------------------------- ! ... parameters to determine cg/ic ratio [price and rind, 1993] !---------------------------------------------------------------------- @@ -184,26 +253,29 @@ subroutine lightning_no_prod( state, pbuf2d, cam_in ) real(r8), parameter :: m2km = 1.e-3_r8 real(r8), parameter :: km2cm = 1.e5_r8 real(r8), parameter :: lat25 = 25._r8*d2r ! 25 degrees latitude in radians - integer :: cldtop_ndx, cldbot_ndx + real(r8) :: flash_freq_land, flash_freq_ocn + real(r8), pointer :: cld2grnd_flash_freq(:) + + if (.not.calc_lightning) return - if (.not.has_no_lightning_prod) return + nullify(cld2grnd_flash_freq) !---------------------------------------------------------------------- ! ... initialization !---------------------------------------------------------------------- flash_freq(:,:) = 0._r8 - prod_no(:,:,:) = 0._r8 - prod_no_col(:,:) = 0._r8 cldhgt(:,:) = 0._r8 dchgzone(:,:) = 0._r8 cgic(:,:) = 0._r8 flash_energy(:,:) = 0._r8 - glob_prod_no_col(:,:) = 0._r8 - cldtop_ndx = pbuf_get_index('CLDTOP') - cldbot_ndx = pbuf_get_index('CLDBOT') + if (calc_nox_prod) then + prod_no(:,:,:) = 0._r8 + prod_no_col(:,:) = 0._r8 + glob_prod_no_col(:,:) = 0._r8 + end if !-------------------------------------------------------------------------------- ! ... estimate flash frequency and resulting no emissions @@ -223,29 +295,30 @@ subroutine lightning_no_prod( state, pbuf2d, cam_in ) ! with 1e17 n atoms per j. the total number of n atoms is then distributed ! over the complete column of grid boxes. !-------------------------------------------------------------------------------- - Chunk_loop : do c = begchunk,endchunk - ncol = state(c)%ncol - lchnk = state(c)%lchnk + Chunk_loop : do lchnk = begchunk,endchunk + ncol = state(lchnk)%ncol + call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk), flsh_frq_ndx, cld2grnd_flash_freq ) call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk), cldtop_ndx, cldtop ) call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk), cldbot_ndx, cldbot ) - zsurf(:ncol) = state(c)%phis(:ncol)*rga - call get_rlat_all_p(c, ncol, rlats(1,c) ) - call get_wght_all_p(c, ncol, wght) + zsurf(:ncol) = state(lchnk)%phis(:ncol)*rga + call get_wght_all_p(lchnk, pcols, wght) do k = 1,pver - zmid(:ncol,k) = state(c)%zm(:ncol,k) + zsurf(:ncol) - zint(:ncol,k,c) = state(c)%zi(:ncol,k) + zsurf(:ncol) + zmid(:ncol,k) = state(lchnk)%zm(:ncol,k) + zsurf(:ncol) + zint(:ncol,k,lchnk) = state(lchnk)%zi(:ncol,k) + zsurf(:ncol) end do - zint(:ncol,pver+1,c) = state(c)%zi(:ncol,pver+1) + zsurf(:ncol) + zint(:ncol,pver+1,lchnk) = state(lchnk)%zi(:ncol,pver+1) + zsurf(:ncol) + + cld2grnd_flash_freq(:) = 0.0_r8 - col_loop : do i = 1,ncol + col_loop : do icol = 1,ncol !-------------------------------------------------------------------------------- ! ... find cloud top and bottom level above 273k !-------------------------------------------------------------------------------- - cldtind = nint( cldtop(i) ) - cldbind = nint( cldbot(i) ) + cldtind = nint( cldtop(icol) ) + cldbind = nint( cldbot(icol) ) do - if( cldbind <= cldtind .or. state(c)%t(i,cldbind) < t0 ) then + if( cldbind <= cldtind .or. state(lchnk)%t(icol,cldbind) < t0 ) then exit end if cldbind = cldbind - 1 @@ -254,58 +327,77 @@ subroutine lightning_no_prod( state, pbuf2d, cam_in ) !-------------------------------------------------------------------------------- ! ... compute cloud top height and depth of charging zone !-------------------------------------------------------------------------------- - cldhgt(i,c) = m2km*max( 0._r8,zint(i,cldtind,c) ) - dchgz = cldhgt(i,c) - m2km*zmid(i,cldbind) - dchgzone(i,c) = dchgz + cldhgt(icol,lchnk) = m2km*max( 0._r8,zint(icol,cldtind,lchnk) ) + dchgz = cldhgt(icol,lchnk) - m2km*zmid(icol,cldbind) + dchgzone(icol,lchnk) = dchgz !-------------------------------------------------------------------------------- ! ... compute flash frequency for given cloud top height ! (flashes storm^-1 min^-1) !-------------------------------------------------------------------------------- - flash_freq_land = 3.44e-5_r8 * cldhgt(i,c)**4.9_r8 - flash_freq_ocn = 6.40e-4_r8 * cldhgt(i,c)**1.7_r8 - flash_freq(i,c) = cam_in(c)%landfrac(i)*flash_freq_land + & - cam_in(c)%ocnfrac(i) *flash_freq_ocn + flash_freq_land = 3.44e-5_r8 * cldhgt(icol,lchnk)**4.9_r8 + flash_freq_ocn = 6.40e-4_r8 * cldhgt(icol,lchnk)**1.7_r8 + flash_freq(icol,lchnk) = cam_in(lchnk)%landfrac(icol)*flash_freq_land + & + cam_in(lchnk)%ocnfrac(icol) *flash_freq_ocn !-------------------------------------------------------------------------------- - ! ... compute cg/ic ratio - ! cgic = proportion of cg flashes (=pg from ppp paper) + ! cgic = proportion of cloud-to-ground flashes + ! NOx from lightning 1. Global distribution based on lightning physics, C Price et al + ! JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 102, NO. D5, PAGES 5929-5941, MARCH 20, 1997 + ! (https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/96JD03504) + ! eq 14 !-------------------------------------------------------------------------------- - cgic(i,c) = 1._r8/((((ca*dchgz + cb)*dchgz + cc) *dchgz + cd)*dchgz + ce) + cgic(icol,lchnk) = 1._r8/((((ca*dchgz + cb)*dchgz + cc) *dchgz + cd)*dchgz + ce) if( dchgz < 5.5_r8 ) then - cgic(i,c) = 0._r8 + cgic(icol,lchnk) = 0._r8 else if( dchgz > 14._r8 ) then - cgic(i,c) = .02_r8 + cgic(icol,lchnk) = .02_r8 end if - !-------------------------------------------------------------------------------- - ! ... compute flash energy (cg*6.7e9 + ic*6.7e8) - ! and convert to total energy per second - ! set ic = cg - !-------------------------------------------------------------------------------- - flash_energy(i,c) = 6.7e9_r8 * flash_freq(i,c)/60._r8 - !-------------------------------------------------------------------------------- - ! ... LKE Aug 23, 2005: scale production to account for different grid - ! box sizes. This requires a reduction in the overall fudge factor - ! (e.g., from 1.2 to 0.5) - !-------------------------------------------------------------------------------- - flash_energy(i,c) = flash_energy(i,c) * wght(i) * geo_factor - !-------------------------------------------------------------------------------- - ! ... compute number of n atoms produced per second - ! and convert to n atoms per second per cm2 and apply fudge factor - !-------------------------------------------------------------------------------- - prod_no_col(i,c) = 1.e17_r8*flash_energy(i,c)/(1.e4_r8*rearth*rearth*wght(i)) * factor - - !-------------------------------------------------------------------------------- - ! ... compute global no production rate in tgn/yr: - ! tgn per second: * 14.00674 * 1.65979e-24 * 1.e-12 - ! nb: 1.65979e-24 = 1/avo - ! tgn per year: * secpyr - !-------------------------------------------------------------------------------- - glob_prod_no_col(i,c) = 1.e17_r8*flash_energy(i,c) & - * 14.00674_r8 * 1.65979e-24_r8 * 1.e-12_r8 * secpyr * factor + cld2grnd_flash_freq(icol) = cam_in(lchnk)%landfrac(icol)*flash_freq_land*cgic(icol,lchnk) ! cld-to-grnd flash frq (per min) + + if (calc_nox_prod) then + !-------------------------------------------------------------------------------- + ! ... compute flash energy (cg*6.7e9 + ic*6.7e8) + ! and convert to total energy per second + ! set ic = cg + !-------------------------------------------------------------------------------- + flash_energy(icol,lchnk) = 6.7e9_r8 * flash_freq(icol,lchnk)/60._r8 + !-------------------------------------------------------------------------------- + ! ... LKE Aug 23, 2005: scale production to account for different grid + ! box sizes. This requires a reduction in the overall fudge factor + ! (e.g., from 1.2 to 0.5) + !-------------------------------------------------------------------------------- + flash_energy(icol,lchnk) = flash_energy(icol,lchnk) * wght(icol) * geo_factor + !-------------------------------------------------------------------------------- + ! ... compute number of n atoms produced per second + ! and convert to n atoms per second per cm2 and apply fudge factor + !-------------------------------------------------------------------------------- + prod_no_col(icol,lchnk) = 1.e17_r8*flash_energy(icol,lchnk)/(1.e4_r8*rearth*rearth*wght(icol)) * factor + + !-------------------------------------------------------------------------------- + ! ... compute global no production rate in tgn/yr: + ! tgn per second: * 14.00674 * 1.65979e-24 * 1.e-12 + ! nb: 1.65979e-24 = 1/avo + ! tgn per year: * secpyr + !-------------------------------------------------------------------------------- + glob_prod_no_col(icol,lchnk) = 1.e17_r8*flash_energy(icol,lchnk) & + * 14.00674_r8 * 1.65979e-24_r8 * 1.e-12_r8 * secpyr * factor + end if end if cloud_layer end do Col_loop + + call outfld( 'LGHTNG_CLD2GRND', cld2grnd_flash_freq, pcols, lchnk ) end do Chunk_loop + + do lchnk = begchunk,endchunk + call outfld( 'FLASHFRQ', flash_freq(:,lchnk), pcols, lchnk ) + call outfld( 'CGIC', cgic(:,lchnk), pcols, lchnk ) + call outfld( 'CLDHGT', cldhgt(:,lchnk), pcols, lchnk ) + call outfld( 'DCHGZONE', dchgzone(:,lchnk), pcols, lchnk ) + enddo + + if (.not.calc_nox_prod) return + !-------------------------------------------------------------------------------- ! ... Accumulate global total, convert to flashes per second ! ... Accumulate global NO production rate @@ -325,29 +417,29 @@ subroutine lightning_no_prod( state, pbuf2d, cam_in ) !-------------------------------------------------------------------------------- ! ... Distribute production up to cloud top [Pickering et al., 1998 (JGR)] !-------------------------------------------------------------------------------- - do c = begchunk,endchunk - ncol = state(c)%ncol - lchnk = state(c)%lchnk + do lchnk = begchunk,endchunk + call get_rlat_all_p(lchnk, pcols, rlats) + ncol = state(lchnk)%ncol call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk), cldtop_ndx, cldtop ) - do i = 1,ncol - cldtind = nint( cldtop(i) ) - if( prod_no_col(i,c) > 0._r8 ) then - if( cldhgt(i,c) > 0._r8 ) then - if( abs( rlats(i,c) ) > lat25 ) then - itype = 1 ! midlatitude continental - else if( nint( cam_in(c)%landfrac(i) ) == land ) then - itype = 3 ! tropical continental + do icol = 1,ncol + cldtind = nint( cldtop(icol) ) + if( prod_no_col(icol,lchnk) > 0._r8 ) then + if( cldhgt(icol,lchnk) > 0._r8 ) then + if( abs( rlats(icol) ) > lat25 ) then + itype = 1 ! midlatitude continental + else if( nint( cam_in(lchnk)%landfrac(icol) ) == land ) then + itype = 3 ! tropical continental else - itype = 2 ! topical marine + itype = 2 ! topical marine end if frac_sum = 0._r8 do k = cldtind,pver - zlow = zint(i,k+1,c) * m2km ! lower interface height (km) - zlow_scal = zlow * 16._r8/cldhgt(i,c) ! scale to 16 km convection height - zlow_ind = max( 1,INT(zlow_scal)+1 ) ! lowest vdist index to include in layer - zhigh = zint(i,k,c) * m2km ! upper interface height (km) - zhigh_scal = zhigh * 16._r8/cldhgt(i,c) ! height (km) scaled to 16km convection height - zhigh_ind = max( 1,MIN( 16,INT(zhigh_scal)+1 ) ) ! highest vdist index to include in layer + zlow = zint(icol,k+1,lchnk) * m2km ! lower interface height (km) + zlow_scal = zlow * 16._r8/cldhgt(icol,lchnk) ! scale to 16 km convection height + zlow_ind = max( 1,INT(zlow_scal)+1 ) ! lowest vdist index to include in layer + zhigh = zint(icol,k,lchnk) * m2km ! upper interface height (km) + zhigh_scal = zhigh * 16._r8/cldhgt(icol,lchnk) ! height (km) scaled to 16km convection height + zhigh_ind = max( 1,MIN( 16,INT(zhigh_scal)+1 ) ) ! highest vdist index to include in layer do kk = zlow_ind,zhigh_ind wrk = kk wrk1 = kk - 1 @@ -355,11 +447,11 @@ subroutine lightning_no_prod( state, pbuf2d, cam_in ) - max( zlow_scal,wrk1 ) fraction = max( 0._r8, min( 1._r8,fraction ) ) frac_sum = frac_sum + fraction*vdist(kk,itype) - prod_no(i,k,c) = prod_no(i,k,c) & ! sum the fraction of column NOx in layer k + prod_no(icol,k,lchnk) = prod_no(icol,k,lchnk) & ! sum the fraction of column NOx in layer k + fraction*vdist(kk,itype)*.01_r8 end do - prod_no(i,k,c) = prod_no_col(i,c) * prod_no(i,k,c) & ! multiply fraction by column amount - / (km2cm*(zhigh - zlow)) ! and convert to atom N cm^-3 s^-1 + prod_no(icol,k,lchnk) = prod_no_col(icol,lchnk) * prod_no(icol,k,lchnk) & ! multiply fraction by column amount + / (km2cm*(zhigh - zlow)) ! and convert to atom N cm^-3 s^-1 end do end if end if @@ -370,15 +462,10 @@ subroutine lightning_no_prod( state, pbuf2d, cam_in ) !-------------------------------------------------------------------------------- ! ... output lightning no production to history file !-------------------------------------------------------------------------------- - do c = begchunk,endchunk - lchnk = state(c)%lchnk - call outfld( 'LNO_PROD', prod_no(:,:,c), pcols, lchnk ) - call outfld( 'LNO_COL_PROD', glob_prod_no_col(:,c), pcols, lchnk ) - call outfld( 'FLASHFRQ', flash_freq(:,c), pcols, lchnk ) - call outfld( 'FLASHENGY', flash_energy(:,c), pcols, lchnk ) - call outfld( 'CLDHGT', cldhgt(:,c), pcols, lchnk ) - call outfld( 'DCHGZONE', dchgzone(:,c), pcols, lchnk ) - call outfld( 'CGIC', cgic(:,c), pcols, lchnk ) + do lchnk = begchunk,endchunk + call outfld( 'LNO_PROD', prod_no(:,:,lchnk), pcols, lchnk ) + call outfld( 'LNO_COL_PROD', glob_prod_no_col(:,lchnk), pcols, lchnk ) + call outfld( 'FLASHENGY', flash_energy(:,lchnk), pcols, lchnk ) enddo end subroutine lightning_no_prod diff --git a/src/chemistry/mozart/ocean_emis.F90 b/src/chemistry/mozart/ocean_emis.F90 index 26819fd846..289cafeb77 100644 --- a/src/chemistry/mozart/ocean_emis.F90 +++ b/src/chemistry/mozart/ocean_emis.F90 @@ -3,23 +3,23 @@ ! Ref: Carpenter et al Chem Soc Rev (2012); Johnson, Ocean sci (2010) ! ------------------------------------------------------------------------------------ ! Required inputs for the air-sea flux module: -! - Seawater concentration (nanomoles per liter) and Sea surface salinity +! - Seawater concentration (nanomoles per liter) and Sea surface salinity ! (parts per thousand) read from namelist (netCDF) ! - Concentration in the gas-phase (pptv), air temperature (K), 10m windspeed (m/s), ! surface pressure (atm), sea surface temperature (K): all from other modules ! ------------------------------------------------------------------------------------ ! Key subroutines: -! ocean_emis_readnl(..): Read salinity from namelist (user_nl_cam). +! ocean_emis_readnl(..): Read salinity from namelist (user_nl_cam). ! Salinity not time-dependent. Flux depends very weakly on it -! ocean_emis_init(...): Interpolate salinity, initialize the library for the flux +! ocean_emis_init(...): Interpolate salinity, initialize the library for the flux ! reading time-dependent seawater conc. from user_nl_cam ! ocean_emis_advance(...): process the seawater concentration -! ocean_emis_getflux(...): calculate the air-sea flux (upward or downward), +! ocean_emis_getflux(...): calculate the air-sea flux (upward or downward), ! then add to total surface flux (sflx) ! ------------------------------------------------------------------------------------ ! Last built: 9 March 2018. ! Written by: Siyuan Wang (ACOM/NCAR) siyuan@ucar.edu -! Acknowledgement: Francis Vitt (NCAR). and of course Dr. Peppurr too +! Acknowledgement: Francis Vitt (NCAR). and of course Dr. Peppurr too ! ==================================================================================== module ocean_emis @@ -33,7 +33,7 @@ module ocean_emis use tracer_data, only : trfld,trfile use chem_mods, only : gas_pcnst use cam_logfile, only : iulog - use ioFileMod, only : getfil + use ioFileMod, only : getfil implicit none @@ -57,9 +57,9 @@ module ocean_emis logical :: switch_bubble type(Csw), allocatable :: Csw_nM(:) - integer :: n_Csw_files + integer :: n_Csw_files - real(r8), allocatable :: salinity(:,:) + real(r8), allocatable :: salinity(:,:) ! ================ ! Air-sea exchange @@ -69,32 +69,32 @@ module ocean_emis Integer, Parameter :: HowManySalts = 5 ! Change this number if you wanna add more salts Integer, Parameter :: HowManySaltProperties = 7 ! Don't touch this (unless you wanna add more fields) - Type GasLib + Type GasLib Character(16) :: CmpdName Real(r8), Dimension(HowManyProperties) :: CmpdProperties End Type GasLib - Type SaltLib + Type SaltLib Character(16) :: SaltName - Real(r8), Dimension(HowManySaltProperties) :: SaltProperties + Real(r8), Dimension(HowManySaltProperties) :: SaltProperties End Type SaltLib Type(GasLib), Dimension(HowManyMolecules) :: GasList ! Library for the trace gas properties Type(SaltLib), Dimension(HowManySalts) :: SaltList ! Library for the salt properties - ! =========================== + ! =========================== ! seawater concentration: ! =========================== - character(len=cl) :: csw_specifier(gas_pcnst) = '' + character(len=cl) :: csw_specifier(gas_pcnst) = '' character(len=24) :: csw_time_type = 'CYCLICAL' ! 'CYCLICAL' | 'SERIAL' | 'INTERP_MISSING_MONTHS' integer :: csw_cycle_yr = 0 - logical :: bubble_mediated_transfer = .false. + logical :: bubble_mediated_transfer = .false. character(len=cl) :: ocean_salinity_file = 'NONE' contains -!-------------------------------------------------------------------------------- -!-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- subroutine ocean_emis_readnl(nlfile) use namelist_utils, only : find_group_name @@ -105,7 +105,7 @@ subroutine ocean_emis_readnl(nlfile) integer :: unitn, ierr character(len=*), parameter :: subname = 'ocean_emis_readnl' - ! =================== + ! =================== ! Namelist definition ! =================== namelist /ocean_emis_nl/ ocean_salinity_file @@ -125,7 +125,7 @@ subroutine ocean_emis_readnl(nlfile) end if close(unitn) end if - + ! ============================ ! Broadcast namelist variables ! ============================ @@ -151,7 +151,7 @@ subroutine ocean_emis_init() use pio, only : file_desc_t, pio_inq_dimid, pio_inq_dimlen, pio_inq_varid, pio_get_var use pio, only : PIO_NOWRITE, PIO_NOERR use pio, only : pio_seterrorhandling, PIO_BCAST_ERROR, pio_closefile - use phys_grid, only : get_ncols_p, get_rlon_all_p, get_rlat_all_p + use phys_grid, only : get_ncols_p, get_rlon_all_p, get_rlat_all_p use interpolate_data, only : lininterp_init, lininterp, interp_type, lininterp_finish use mo_constants, only : pi @@ -162,19 +162,19 @@ subroutine ocean_emis_init() real(r8), allocatable :: file_lats(:), file_lons(:) real(r8), allocatable :: wrk2d(:,:) real(r8) :: to_lats(pcols), to_lons(pcols) - type(interp_type) :: lon_wgts, lat_wgts + type(interp_type) :: lon_wgts, lat_wgts real(r8), parameter :: zero=0_r8, twopi=2_r8*pi, degs2rads = pi/180._r8 character(len=*), parameter :: subname = 'ocean_emis_init' - + if (trim(ocean_salinity_file) == 'NONE') return call getfil( ocean_salinity_file, filen, 0 ) call cam_pio_openfile( fid, filen, PIO_NOWRITE) - + call pio_seterrorhandling(fid, PIO_BCAST_ERROR) - + ierr = pio_inq_dimid( fid, 'lon', dimid ) if (ierr /= PIO_NOERR) then call endrun(subname//': pio_inq_dimid lon FAILED') @@ -225,6 +225,7 @@ subroutine ocean_emis_init() endif allocate(salinity(pcols,begchunk:endchunk)) + salinity = 0._r8 do c=begchunk,endchunk @@ -235,17 +236,22 @@ subroutine ocean_emis_init() call lininterp_init(file_lons, file_nlon, to_lons, ncols, 2, lon_wgts, zero, twopi) call lininterp_init(file_lats, file_nlat, to_lats, ncols, 1, lat_wgts) - call lininterp(wrk2d, file_nlon, file_nlat, salinity(1:ncols,c), ncols, lon_wgts, lat_wgts) + call lininterp(wrk2d, file_nlon, file_nlat, salinity(1:ncols,c), ncols, lon_wgts, lat_wgts) call lininterp_finish(lon_wgts) call lininterp_finish(lat_wgts) end do + ! fill in missing values with climatology for modern-day + where(salinity < 0._r8) + salinity = 33.0_r8 + end where + deallocate( file_lons, file_lats ) deallocate( wrk2d ) - call addfld('OCN_SALINITY', horiz_only, 'A', 'parts per thousands', 'ocean salinity' ) + call addfld('OCN_SALINITY', horiz_only, 'A', 'parts per thousands', 'ocean salinity' ) ! ====================================================== ! initializing the libraries for the air-sea flux module @@ -253,17 +259,17 @@ subroutine ocean_emis_init() Call CmpLibInitialization() Call SaltLibInitialization() - ! --------------------------------------------- + ! --------------------------------------------- ! Read seawater concentration: WSY ! --------------------------------------------- call cseawater_ini() call pio_closefile (fid) - + end subroutine ocean_emis_init -!-------------------------------------------------------------------------------- -!-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- subroutine ocean_emis_advance( pbuf2d, state ) ! ------------------------------- ! check serial case for time span @@ -274,7 +280,7 @@ subroutine ocean_emis_advance( pbuf2d, state ) use tracer_data, only : advance_trcdata use physics_buffer, only : physics_buffer_desc - type(physics_state), intent(in) :: state(begchunk:endchunk) + type(physics_state), intent(in) :: state(begchunk:endchunk) type(physics_buffer_desc), pointer :: pbuf2d(:,:) integer :: m @@ -286,12 +292,12 @@ subroutine ocean_emis_advance( pbuf2d, state ) end subroutine ocean_emis_advance -!-------------------------------------------------------------------------------- -!-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- subroutine ocean_emis_getflux(lchnk, ncol, state, u10, sst, ocnfrac, icefrac, sflx) use physics_types, only : physics_state - use ppgrid, only : pver + use ppgrid, only : pver integer, intent(in) :: lchnk, ncol type(physics_state), target, intent(in) :: state ! Physics state variables @@ -301,13 +307,13 @@ subroutine ocean_emis_getflux(lchnk, ncol, state, u10, sst, ocnfrac, icefrac, sf real(r8), intent(in) :: icefrac(:) ! Ice fraction real(r8), intent(inout) :: sflx(:,:) ! Surface emissions (kg/m^2/s) - integer :: m, isec, SpeciesID - real(r8) :: Csw_col(ncol) - real(r8) :: MW_species - real(r8) :: oceanflux_kg_m2_s(ncol) + integer :: i, m, isec, SpeciesID + real(r8) :: Csw_col(ncol) + real(r8) :: MW_species + real(r8) :: oceanflux_kg_m2_s(ncol) if (trim(ocean_salinity_file) == 'NONE') return - + ! ================================================== ! Get seawater concentrations and calculate the flux ! ================================================== @@ -317,28 +323,30 @@ subroutine ocean_emis_getflux(lchnk, ncol, state, u10, sst, ocnfrac, icefrac, sf isec = 1 Csw_col(:ncol) = Csw_nM(m)%scalefactor*Csw_nM(m)%fields(isec)%data(:ncol,1,lchnk) - MW_species = MolecularWeight(SpeciesIndex( Csw_nM(m)%species )) + MW_species = MolecularWeight(SpeciesIndex( Csw_nM(m)%species )) call cnst_get_ind( trim(Csw_nM(m)%species), SpeciesID, abort=.true. ) oceanflux_kg_m2_s = 0.0_r8 - where (ocnfrac(:ncol) >= 0.2_r8 .and. Csw_col(:ncol) >= 0._r8) ! calculate flux only for ocean - oceanflux_kg_m2_s(:ncol) = Flux_kg_m2_s( & - Csw_nM(m)%species, & ! name of species - state%q(:ncol,pver,SpeciesID) * (28.97_r8/MW_species) * 1.0e+12_r8, & ! air concentration (ppt) - Csw_col(:ncol), & ! sea water concentration (nM) - state%t(:ncol,pver), & ! air temperature (K) - u10(:ncol), & ! wind speed at 10m (m/s) <- should use this - state%ps(:ncol) / 101325.0_r8, & ! surface pressure (atm) - sst(:ncol), & ! sea surface temperautre (K) - salinity(:ncol,lchnk), & ! ocean salinity (parts per thousands) - switch_bubble, & ! bubble-mediated transfer: on or off - ncol ) - end where + do i = 1,ncol + if (ocnfrac(i) >= 0.2_r8 .and. Csw_col(i) >= 0._r8) then + ! calculate flux only for ocean + oceanflux_kg_m2_s(i) = Flux_kg_m2_s( & + Csw_nM(m)%species, & ! name of species + state%q(i,pver,SpeciesID) * (28.97_r8/MW_species) * 1.0e+12_r8, & ! air concentration (ppt) + Csw_col(i), & ! sea water concentration (nM) + state%t(i,pver), & ! air temperature (K) + u10(i), & ! wind speed at 10m (m/s) <- should use this + state%ps(i) / 101325.0_r8, & ! surface pressure (atm) + sst(i), & ! sea surface temperautre (K) + salinity(i,lchnk), & ! ocean salinity (parts per thousands) + switch_bubble ) ! bubble-mediated transfer: on or off + end if + end do ! =========================================================================== - ! Add the ocean flux to the other fluxes + ! Add the ocean flux to the other fluxes ! Make sure this ocean module is called after other surface emissions are set ! =========================================================================== sflx(:ncol,SpeciesID) = sflx(:ncol,SpeciesID) + oceanflux_kg_m2_s(:ncol) * ocnfrac(:ncol) @@ -355,10 +363,8 @@ subroutine ocean_emis_getflux(lchnk, ncol, state, u10, sst, ocnfrac, icefrac, sf end subroutine ocean_emis_getflux - -!-------------------------------------------------------------------------------- -!-------------------------------------------------------------------------------- - + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- Subroutine CmpLibInitialization() ! ===================================================================================== ! This is the lookup table for molecular weight, Vb, and Henry's law constant @@ -377,7 +383,7 @@ Subroutine CmpLibInitialization() 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, & 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 /)) 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, & - 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 /)) + 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 /)) 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, & 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/)) 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, & @@ -409,10 +415,12 @@ Subroutine CmpLibInitialization() ! -------------------------------------------------------------------------------- End Subroutine CmpLibInitialization + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- Subroutine SaltLibInitialization() ! ================================================================================ - ! This is the lookup table for common solutes in seawater and the parameters to - ! calculate the dynamic viscosity of seawater. + ! This is the lookup table for common solutes in seawater and the parameters to + ! calculate the dynamic viscosity of seawater. ! You may add other solutes or change the mass fractions. ! -------------------------------------------------------------------------------- ! Col 1: mass fraction of solute @@ -431,6 +439,8 @@ Subroutine SaltLibInitialization() ! --------------------------------------------- End Subroutine SaltLibInitialization + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- Function SpeciesIndex(SpeciesName) ! ============================================== ! This function is to look for the species index @@ -439,7 +449,7 @@ Function SpeciesIndex(SpeciesName) Character(Len=16) :: SpeciesName SpeciesIndex = -1 ! return -1 if species is not found - + Do i = 1, HowManyMolecules If (trim(SpeciesName) == trim(GasList(i)%CmpdName)) Then SpeciesIndex = i @@ -448,13 +458,15 @@ Function SpeciesIndex(SpeciesName) End Do End Function SpeciesIndex - Function Flux_kg_m2_s(SpeciesName,Cgas_ppt,Cwater_nM,T_air_K,u10_m_s,P_atm,T_water_K,& - Salinity_PartsPerThousand,switch_bubble,ncol) + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + Real(r8) Function Flux_kg_m2_s(SpeciesName,Cgas_ppt,Cwater_nM,T_air_K,u10_m_s,P_atm,T_water_K,& + Salinity_PartsPerThousand,switch_bubble) ! =========================================================================== ! This is the main module function. Input variables: ! --------------------------------------------------------------------------- ! - SpeciesName: name of species - ! - Cgas_ppt: mixing ratio (parts per trillion) of trace gas of interest + ! - Cgas_ppt: mixing ratio (parts per trillion) of trace gas of interest ! in the gas-phase (lowest modeling layer) ! - Cwater_nM: concentration of trace gas of interest in the surface ocean ! - T_air_K: temperature in the lowest modeling layer @@ -463,52 +475,51 @@ Function Flux_kg_m2_s(SpeciesName,Cgas_ppt,Cwater_nM,T_air_K,u10_m_s,P_atm,T_wat ! - T_water_K: sea surface temperature ! - Salinity_PartsPerThousand: surface ocean salinity ! - switch_bubble: bubble-mediated transfer switch - ! All must be 1D arrays with same dimension(ncol, so CESM-compatible) ! =========================================================================== - Integer :: ncol, SpeciesID - Character(16) :: SpeciesName - Real(r8), Dimension(ncol) :: Flux_kg_m2_s - Real(r8), Dimension(ncol) :: Cgas_ppt, Cwater_nM, T_air_K, u10_m_s, P_atm, T_water_K, Salinity_PartsPerThousand - Real(r8), Dimension(ncol) :: H_gas_over_liquid_dimless, kt_m_s - Logical :: switch_bubble + Character(16),intent(in) :: SpeciesName + Real(r8),intent(in) :: Cgas_ppt, Cwater_nM, T_air_K, u10_m_s, P_atm, T_water_K, Salinity_PartsPerThousand + Logical ,intent(in) :: switch_bubble - where(Salinity_PartsPerThousand .lt. 0.0_r8) Salinity_PartsPerThousand = 33.0_r8 + Integer :: SpeciesID + Real(r8) :: H_gas_over_liquid_dimless, kt_m_s - SpeciesID = SpeciesIndex(SpeciesName) - H_gas_over_liquid_dimless = 1.0_r8/(Henry_M_atm(SpeciesID,T_water_K,Salinity_PartsPerThousand,ncol)*& + SpeciesID = SpeciesIndex(SpeciesName) + H_gas_over_liquid_dimless = 1.0_r8/(Henry_M_atm(SpeciesID,T_water_K,Salinity_PartsPerThousand)*& 0.082_r8*T_water_K) If (switch_bubble) then ! -------------------------------------------------------- ! k_water parameterization with bubble-induced enhancement ! -------------------------------------------------------- kt_m_s = (1.0_r8/k_water_m_s_bubble(SpeciesID, T_water_K, Salinity_PartsPerThousand, & - u10_m_s, Cgas_ppt, P_atm, T_air_K, ncol) & - + 1.0_r8/k_air_m_s(SpeciesID, u10_m_s, T_air_K, P_atm, ncol)& + u10_m_s, Cgas_ppt, P_atm, T_air_K) & + + 1.0_r8/k_air_m_s(SpeciesID, u10_m_s, T_air_K, P_atm)& /H_gas_over_liquid_dimless)**(-1.0_r8) else ! ------------------------------------------------ ! Original k_water parameterization, scaled to CO2 ! ------------------------------------------------ - kt_m_s = (1.0_r8/k_water_m_s(SpeciesID, T_water_K, Salinity_PartsPerThousand, u10_m_s, ncol) & - + 1.0_r8/k_air_m_s(SpeciesID, u10_m_s, T_air_K, P_atm, ncol)/H_gas_over_liquid_dimless)**(-1.0_r8) + kt_m_s = (1.0_r8/k_water_m_s(SpeciesID, T_water_K, Salinity_PartsPerThousand, u10_m_s) & + + 1.0_r8/k_air_m_s(SpeciesID, u10_m_s, T_air_K, P_atm)/H_gas_over_liquid_dimless)**(-1.0_r8) endif Flux_kg_m2_s = kt_m_s * (Cwater_nM*1E-9_r8*1000.0_r8 & - Cgas_ppt*1E-12_r8*(101325.0_r8*P_atm)/8.314_r8/T_air_K/H_gas_over_liquid_dimless) & ! g/m2/s * MolecularWeight(SpeciesIndex(SpeciesName)) / 1000.0_r8 ! convert to kg/m2/s End Function Flux_kg_m2_s - - Function k_air_m_s(SpeciesIndex, u10_m_s, T_air_K, P_atm, ncol) + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + Real(r8) Function k_air_m_s(SpeciesIndex, u10_m_s, T_air_K, P_atm) use shr_const_mod, only: vonKarman=>SHR_CONST_KARMAN ! ============================================================================= - ! Air-side transfer velocity. Slightly modified NOAA COARE (Fairall et al 2003; - ! Feffery et al 2010), as recommended by Johnson Ocean Sci. 2010. + ! Air-side transfer velocity. Slightly modified NOAA COARE (Fairall et al 2003; + ! Feffery et al 2010), as recommended by Johnson Ocean Sci. 2010. ! Dynamic viscosity of air: Tsilingiris 2008 ! ============================================================================= - Integer :: ncol, SpeciesIndex - Real(r8), Dimension(ncol) :: k_air_m_s - Real(r8), Dimension(ncol) :: u10_m_s, T_air_K, P_atm, ustar_m_s, DragCoeff - Real(r8), Dimension(ncol) :: DynamicViscosityAir_kg_m_s, DensityAir_kg_m3, DiffusivityInAir, SchmidtNumberInAir + Integer ,intent(in) :: SpeciesIndex + Real(r8),intent(in) :: u10_m_s, T_air_K, P_atm + + Real(r8) :: ustar_m_s, DragCoeff + Real(r8) :: DynamicViscosityAir_kg_m_s, DensityAir_kg_m3, DiffusivityInAir, SchmidtNumberInAir ! WSY: If local friction velocity is available from the model, might as well use that? ustar_m_s = u10_m_s * sqrt(6.1E-4_r8 + 6.3E-5_r8 * u10_m_s) @@ -516,53 +527,53 @@ Function k_air_m_s(SpeciesIndex, u10_m_s, T_air_K, P_atm, ncol) DynamicViscosityAir_kg_m_s = 1.715747771E-5_r8 + 4.722402075E-8_r8 * (T_air_K-273.15_r8) & - 3.663027156E-10_r8 * ((T_air_K-273.15_r8)**2.0_r8) & + 1.873236686E-12_r8 * ((T_air_K-273.15_r8)**3.0_r8) & - - 8.050218737E-14_r8 * ((T_air_K-273.15_r8)**4.0_r8) + - 8.050218737E-14_r8 * ((T_air_K-273.15_r8)**4.0_r8) DensityAir_kg_m3 = 1.293393662_r8 - 5.538444326e-3_r8 * (T_air_K-273.15_r8) & + 3.860201577e-5_r8 * (T_air_K-273.15_r8)**2.0_r8 & - 5.2536065e-7_r8 * (T_air_K-273.15_r8)**3.0_r8 - DiffusivityInAir = DiffusivityInAir_cm2_s(SpeciesIndex, T_air_K, P_atm, ncol) - SchmidtNumberInAir = DynamicViscosityAir_kg_m_s / DensityAir_kg_m3 / (DiffusivityInAir/10000.0_r8) + DiffusivityInAir = DiffusivityInAir_cm2_s(SpeciesIndex, T_air_K, P_atm) + SchmidtNumberInAir = DynamicViscosityAir_kg_m_s / DensityAir_kg_m3 / (DiffusivityInAir/10000.0_r8) k_air_m_s = 1E-3_r8 + ustar_m_s / (13.3_r8*(SchmidtNumberInAir**0.5_r8)+(DragCoeff**(-0.5_r8))-& 5.0_r8+log(SchmidtNumberInAir)/2.0_r8/vonKarman) End Function k_air_m_s - - - - Function k_water_m_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, u10_m_s, ncol) + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + Real(r8) Function k_water_m_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, u10_m_s) ! ================================================================================ ! Water-side transfer velocity. Ref: Nightingale et al (2000). Salinity considered ! ================================================================================ - Integer :: ncol, SpeciesIndex - Real(r8), Dimension(ncol) :: k_water_m_s - Real(r8), Dimension(ncol) :: T_water_K, Salinity_PartsPerThousand, u10_m_s - Real(r8), Dimension(ncol) :: DiffusivityInWater, SchmidtNumberInWater - Real(r8) :: SchmidtNumberInWater_CO2ref + Integer ,intent(in) :: SpeciesIndex + Real(r8),intent(in) :: T_water_K, Salinity_PartsPerThousand, u10_m_s + + Real(r8) :: DiffusivityInWater, SchmidtNumberInWater + Real(r8) :: SchmidtNumberInWater_CO2ref + SchmidtNumberInWater_CO2ref = 660.0_r8 ! this is the Schmidt number of CO2 at 20 degC in fresh water - DiffusivityInWater = DiffusivityInWater_cm2_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, ncol) - SchmidtNumberInWater = DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand, ncol) / 1000.0_r8 & - / DensityWater_kg_m3(T_water_K,Salinity_PartsPerThousand,ncol)/(DiffusivityInWater/10000.0_r8) + DiffusivityInWater = DiffusivityInWater_cm2_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand) + SchmidtNumberInWater = DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand) / 1000.0_r8 & + / DensityWater_kg_m3(T_water_K,Salinity_PartsPerThousand)/(DiffusivityInWater/10000.0_r8) k_water_m_s = ((0.222_r8*(u10_m_s**2.0_r8)+0.333_r8*u10_m_s)*& ((SchmidtNumberInWater/SchmidtNumberInWater_CO2ref)**(-0.5_r8)))/360000.0_r8 End Function k_water_m_s - - - - Function k_water_m_s_bubble(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, u10_m_s, Cgas_ppt, P_atm, T_air_K, ncol) + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + Real(r8) Function k_water_m_s_bubble(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, u10_m_s, Cgas_ppt, P_atm, T_air_K) ! ============================================================== ! Water-side transfer velocity. Ref: Asher and Wanninkhof (1998). ! ============================================================== - Integer :: ncol, SpeciesIndex - Real(r8), Dimension(ncol) :: k_water_m_s_bubble - Real(r8), Dimension(ncol) :: T_water_K, Salinity_PartsPerThousand, u10_m_s, Cgas_ppt, P_atm, T_air_K - Real(r8), Dimension(ncol) :: DiffusivityInWater, SchmidtNumberInWater - Real(r8), Dimension(ncol) :: FracCoverage_WhiteCaps, OstwaldSolubilityCoefficient - DiffusivityInWater = DiffusivityInWater_cm2_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, ncol) - SchmidtNumberInWater = DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand, ncol) / 1000.0_r8 & - / DensityWater_kg_m3(T_water_K,Salinity_PartsPerThousand,ncol)/(DiffusivityInWater/10000.0_r8) - FracCoverage_WhiteCaps = 2.56e-6_r8 * (u10_m_s - 1.77_r8)**3.0_r8 - OstwaldSolubilityCoefficient = Henry_M_atm(SpeciesIndex,T_water_K,Salinity_PartsPerThousand,ncol) ! just Henry's law (M/atm) + Integer, intent(in) :: SpeciesIndex + Real(r8),intent(in) :: T_water_K, Salinity_PartsPerThousand, u10_m_s, Cgas_ppt, P_atm, T_air_K + + Real(r8) :: DiffusivityInWater, SchmidtNumberInWater + Real(r8) :: FracCoverage_WhiteCaps, OstwaldSolubilityCoefficient + + DiffusivityInWater = DiffusivityInWater_cm2_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand) + SchmidtNumberInWater = DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand) / 1000.0_r8 & + / DensityWater_kg_m3(T_water_K,Salinity_PartsPerThousand)/(DiffusivityInWater/10000.0_r8) + FracCoverage_WhiteCaps = 2.56e-6_r8 * (u10_m_s - 1.77_r8)**3.0_r8 + OstwaldSolubilityCoefficient = Henry_M_atm(SpeciesIndex,T_water_K,Salinity_PartsPerThousand) ! just Henry's law (M/atm) OstwaldSolubilityCoefficient = OstwaldSolubilityCoefficient * (Cgas_ppt*1.0E-12_r8*P_atm) ! mol / L OstwaldSolubilityCoefficient = OstwaldSolubilityCoefficient * 0.082_r8 * T_air_K / P_atm ! L / L k_water_m_s_bubble = ((47.0_r8*u10_m_s + FracCoverage_WhiteCaps*(115200.0_r8 - 47.0_r8* u10_m_s)) & @@ -570,40 +581,46 @@ Function k_water_m_s_bubble(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, + FracCoverage_WhiteCaps * (-37.0_r8/OstwaldSolubilityCoefficient & + 6120.0_r8*(OstwaldSolubilityCoefficient**(-0.37_r8)) *(SchmidtNumberInWater**(-0.18_r8)))) & * 2.8e-6_r8 - End Function k_water_m_s_bubble - - + End Function k_water_m_s_bubble - Function DiffusivityInAir_cm2_s(SpeciesIndex, T_air_K, P_atm, ncol) + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + Real(r8) Function DiffusivityInAir_cm2_s(SpeciesIndex, T_air_K, P_atm) ! ============================ ! Ref: Johnson Ocean Sci. 2010 ! ============================ - Integer :: ncol, SpeciesIndex - Real(r8), Dimension(ncol) :: DiffusivityInAir_cm2_s, T_air_K, P_atm + Integer ,intent(in) :: SpeciesIndex + Real(r8),intent(in) :: T_air_K, P_atm + Real(r8), parameter :: MW_air = 28.97_r8 ! molecular weight for air Real(r8), parameter :: Va = 20.1_r8 ! molar volume for air Real(r8) :: Vb, MW_species + Vb = LiquidMolarVolume_cm3_mol(SpeciesIndex) MW_species = MolecularWeight(SpeciesIndex) DiffusivityInAir_cm2_s = 0.001_r8 * (T_air_K**1.75_r8) & ! oh f* me * (((MW_air + MW_species)/(MW_air*MW_species))**0.5_r8) & / ((P_atm*(Va**(1.0_r8/3.0_r8)+Vb**(1.0_r8/3.0_r8)))**2.0_r8) - End Function DiffusivityInAir_cm2_s - + End Function DiffusivityInAir_cm2_s - Function DiffusivityInWater_cm2_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, ncol) + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + Real(r8) Function DiffusivityInWater_cm2_s(SpeciesIndex, T_water_K, Salinity_PartsPerThousand) ! ================================================= ! Ref: Johnson Ocean Sci. 2010. Salinity considered ! ================================================= - Integer :: ncol, SpeciesIndex - Real(r8), Dimension(ncol) :: DiffusivityInWater_cm2_s, DynamicViscosityWater, T_water_K, Salinity_PartsPerThousand + Integer, intent(in) :: SpeciesIndex + Real(r8),intent(in) :: T_water_K, Salinity_PartsPerThousand + Real(r8), parameter :: AssociationFactor = 2.6_r8 ! ... for water - Real(r8) :: Vb, MW_species + Real(r8) :: DynamicViscosityWater, Vb, MW_species + Vb = LiquidMolarVolume_cm3_mol(SpeciesIndex) MW_species = MolecularWeight(SpeciesIndex) - DynamicViscosityWater = DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand, ncol) + + DynamicViscosityWater = DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand) ! ------------------------------------------------- ! Wilke and Chang 1955: this seems to be a bit high ! ------------------------------------------------- @@ -617,47 +634,51 @@ Function DiffusivityInWater_cm2_s(SpeciesIndex, T_water_K, Salinity_PartsPerThou End Function DiffusivityInWater_cm2_s - - Function DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand, ncol) + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + Real(r8) Function DynamicViscosityWater_g_m_s(T_water_K, Salinity_PartsPerThousand) ! ================================================= ! Ref: Johnson Ocean Sci. 2010. Salinity considered ! ================================================= - Integer :: ncol - Real(r8), Dimension(ncol) :: DynamicViscosityWater_g_m_s, T_water_K, Salinity_PartsPerThousand - Real(r8), Dimension(ncol) :: MassFrac_water, DynamicViscosityPureWater_g_m_s, SaltViscosity, sum_w_ln_SaltViscosity - Integer :: j, n + Real(r8),intent(in) :: T_water_K, Salinity_PartsPerThousand + + Real(r8) :: MassFrac_water, DynamicViscosityPureWater_g_m_s, SaltViscosity, sum_w_ln_SaltViscosity + Integer :: n + sum_w_ln_SaltViscosity = 0.0_r8 MassFrac_water = 1.0_r8 - Salinity_PartsPerThousand / 1000.0_r8 DynamicViscosityPureWater_g_m_s = ((T_water_K-273.15_r8)+246.0_r8) & - / (0.05594_r8*(T_water_K-273.15_r8)**2.0_r8+5.2842_r8*(T_water_K-273.15_r8)+137.37_r8) - Do j = 1, ncol - If (Salinity_PartsPerThousand(j) == 0.0_r8) Then ! pure water - DynamicViscosityWater_g_m_s(j) = DynamicViscosityPureWater_g_m_s(j) + / (0.05594_r8*(T_water_K-273.15_r8)**2.0_r8+5.2842_r8*(T_water_K-273.15_r8)+137.37_r8) + + If (Salinity_PartsPerThousand == 0.0_r8) Then ! pure water + DynamicViscosityWater_g_m_s = DynamicViscosityPureWater_g_m_s Else ! salty water Do n = 1, HowManySalts - SaltViscosity(j) = exp((SaltList(n)%SaltProperties(2) * & - (Salinity_PartsPerThousand(j)/1000.0_r8)**SaltList(n)%SaltProperties(3) & + SaltViscosity = exp((SaltList(n)%SaltProperties(2) * & + (Salinity_PartsPerThousand/1000.0_r8)**SaltList(n)%SaltProperties(3) & + SaltList(n)%SaltProperties(4)) & - / (SaltList(n)%SaltProperties(5)*(T_water_K(j)-273.15_r8) + 1.0_r8)) & - / (SaltList(n)%SaltProperties(6) * (Salinity_PartsPerThousand(j) / & + / (SaltList(n)%SaltProperties(5)*(T_water_K-273.15_r8) + 1.0_r8)) & + / (SaltList(n)%SaltProperties(6) * (Salinity_PartsPerThousand / & 1000.0_r8)**SaltList(n)%SaltProperties(7) + 1.0_r8) - sum_w_ln_SaltViscosity(j) = sum_w_ln_SaltViscosity(j) + (Salinity_PartsPerThousand(j)/1000.0_r8) & - * SaltList(n)%SaltProperties(1) * log(SaltViscosity(j)) + sum_w_ln_SaltViscosity = sum_w_ln_SaltViscosity + (Salinity_PartsPerThousand/1000.0_r8) & + * SaltList(n)%SaltProperties(1) * log(SaltViscosity) End Do - DynamicViscosityWater_g_m_s(j) = exp(MassFrac_water(j) & - * log(DynamicViscosityPureWater_g_m_s(j)) + sum_w_ln_SaltViscosity(j)) + DynamicViscosityWater_g_m_s = exp(MassFrac_water & + * log(DynamicViscosityPureWater_g_m_s) + sum_w_ln_SaltViscosity) Endif - End Do - End Function DynamicViscosityWater_g_m_s + End Function DynamicViscosityWater_g_m_s - Function DensityWater_kg_m3(T_water_K, Salinity_PartsPerThousand, ncol) + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + Real(r8) Function DensityWater_kg_m3(T_water_K, Salinity_PartsPerThousand) ! ==================================================== ! Ref: Millero and Poisson (1981). Salinity considered ! ==================================================== - Integer :: ncol - Real(r8), Dimension(ncol) :: DensityWater_kg_m3, T_water_K, Salinity_PartsPerThousand - Real(r8), Dimension(ncol) :: DensityPureWater_kg_m3, FactorA, FactorB, FactorC + Real(r8), intent(in) :: T_water_K, Salinity_PartsPerThousand + + Real(r8) :: DensityPureWater_kg_m3, FactorA, FactorB, FactorC + DensityPureWater_kg_m3 = 999.842594_r8 + 0.06793952_r8*(T_water_K-273.15_r8) & - 0.00909529_r8*((T_water_K-273.15_r8)**2.0_r8) & + 0.0001001685_r8*((T_water_K-273.15_r8)**3.0_r8) & @@ -669,41 +690,46 @@ Function DensityWater_kg_m3(T_water_K, Salinity_PartsPerThousand, ncol) FactorC = 0.00048314_r8 DensityWater_kg_m3 = DensityPureWater_kg_m3 + FactorA*Salinity_PartsPerThousand & + FactorB*(Salinity_PartsPerThousand**(2.0_r8/3.0_r8)) + FactorC*Salinity_PartsPerThousand - End Function DensityWater_kg_m3 + End Function DensityWater_kg_m3 - Function Henry_M_atm(SpeciesIndex, T_water_K, Salinity_PartsPerThousand, ncol) + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- + Real(r8) Function Henry_M_atm(SpeciesIndex, T_water_K, Salinity_PartsPerThousand) ! ========================================================================================= ! Ref: Sander compilation 2015. Salt-in or salt-out estimated based on Setschenow constants ! ========================================================================================= - Integer :: ncol, j - Integer :: SpeciesIndex - Real(r8), Dimension(ncol) :: Henry_M_atm, T_water_K, Salinity_PartsPerThousand - Real(r8), Dimension(ncol) :: Heff_M_atm_PureWater, Setschenow, Heff_M_atm_SaltyWater + Integer, intent(in) :: SpeciesIndex + Real(r8), intent(in) :: T_water_K, Salinity_PartsPerThousand + + Real(r8) :: Heff_M_atm_PureWater, Setschenow, Heff_M_atm_SaltyWater + Heff_M_atm_PureWater = GasList(SpeciesIndex)%CmpdProperties(15) * & exp(GasList(SpeciesIndex)%CmpdProperties(16) * (1.0_r8/T_water_K - 1.0_r8/298.0_r8)) - Do j = 1, ncol - If (Salinity_PartsPerThousand(j)==0.0_r8) Then - Henry_M_atm(j) = Heff_M_atm_PureWater(j) - Else - Setschenow(j) = log(LiquidMolarVolume_cm3_mol(SpeciesIndex)) * & - (7.33532E-4_r8 + 3.39615E-5_r8 * log(Heff_M_atm_PureWater(j)) & - - 2.40888E-6_r8 * ((log(Heff_M_atm_PureWater(j)))**2.0_r8) & - + 1.57114E-7_r8 * ((log(Heff_M_atm_PureWater(j)))**3.0_r8)) - Heff_M_atm_SaltyWater(j) = Heff_M_atm_PureWater(j) * 10.0_r8**(Setschenow(j)*Salinity_PartsPerThousand(j)) - Henry_M_atm(j) = Heff_M_atm_SaltyWater(j) - Endif - End Do - End Function Henry_M_atm + If (Salinity_PartsPerThousand==0.0_r8) Then + Henry_M_atm = Heff_M_atm_PureWater + Else + Setschenow = log(LiquidMolarVolume_cm3_mol(SpeciesIndex)) * & + (7.33532E-4_r8 + 3.39615E-5_r8 * log(Heff_M_atm_PureWater) & + - 2.40888E-6_r8 * ((log(Heff_M_atm_PureWater))**2.0_r8) & + + 1.57114E-7_r8 * ((log(Heff_M_atm_PureWater))**3.0_r8)) + Heff_M_atm_SaltyWater = Heff_M_atm_PureWater * 10.0_r8**(Setschenow*Salinity_PartsPerThousand) + Henry_M_atm = Heff_M_atm_SaltyWater + Endif + + End Function Henry_M_atm + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- Function MolecularWeight(SpeciesIndex) Real(r8) :: MolecularWeight Integer :: SpeciesIndex MolecularWeight = GasList(SpeciesIndex)%CmpdProperties(1) End Function MolecularWeight - + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- Function LiquidMolarVolume_cm3_mol(SpeciesIndex) ! =========================================================================== ! If no measurements available, i.e. GasList(SpeciesIndex)%CmpdProperties(14) @@ -712,7 +738,7 @@ Function LiquidMolarVolume_cm3_mol(SpeciesIndex) Real(r8) :: LiquidMolarVolume_cm3_mol Integer :: SpeciesIndex - If (GasList(SpeciesIndex)%CmpdProperties(14)/=0.0_r8) Then + If (GasList(SpeciesIndex)%CmpdProperties(14)/=0.0_r8) Then LiquidMolarVolume_cm3_mol = GasList(SpeciesIndex)%CmpdProperties(14) Else LiquidMolarVolume_cm3_mol = 7.0_r8*GasList(SpeciesIndex)%CmpdProperties(2) ! C @@ -731,18 +757,20 @@ Function LiquidMolarVolume_cm3_mol(SpeciesIndex) End Function LiquidMolarVolume_cm3_mol + !-------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------- subroutine cseawater_ini() - use mo_chem_utls, only : get_spc_ndx - use tracer_data, only : trcdata_init - use cam_pio_utils, only : cam_pio_openfile + use mo_chem_utls, only : get_spc_ndx + use tracer_data, only : trcdata_init + use cam_pio_utils, only : cam_pio_openfile use pio, only : pio_inquire, pio_nowrite, pio_closefile, pio_inq_varndims use pio, only : pio_inq_varname, file_desc_t, pio_get_att, PIO_NOERR, PIO_GLOBAL - use pio, only : pio_seterrorhandling, PIO_BCAST_ERROR - use string_utils, only : GLC + use pio, only : pio_seterrorhandling, PIO_BCAST_ERROR + use string_utils, only : GLC integer :: i, j, l, m, n, nn, astat, vid, ierr, nvars, isec - integer :: indx(gas_pcnst) + integer :: indx(gas_pcnst) type(file_desc_t) :: ncid character(len=16) :: csw_species(gas_pcnst) character(len=256) :: csw_filenam(gas_pcnst) @@ -766,7 +794,7 @@ subroutine cseawater_ini() character(len=*), parameter :: subname = 'cseawater_ini' - ! ======================================================== + ! ======================================================== ! Read sea water concentration specifier from the namelist ! ======================================================== @@ -827,7 +855,7 @@ subroutine cseawater_ini() ! ------------------------------------------- ! Setup the seawater concentration type array ! ------------------------------------------- - do m=1,n_Csw_files + do m=1,n_Csw_files Csw_nM(m)%spc_ndx = csw_indexes(indx(m)) Csw_nM(m)%units = 'nM' Csw_nM(m)%species = csw_species(indx(m)) @@ -898,9 +926,9 @@ subroutine cseawater_ini() deallocate(vndims) ! Global attribute 'input_method' overrides the srf_emis_type namelist setting on - ! a file-by-file basis. If the emis file does not contain the 'input_method' + ! a file-by-file basis. If the emis file does not contain the 'input_method' ! attribute then the srf_emis_type namelist setting is used. - ierr = pio_get_att(ncid, PIO_GLOBAL, 'input_method', file_interp_type) + ierr = pio_get_att(ncid, PIO_GLOBAL, 'input_method', file_interp_type) if ( ierr == PIO_NOERR) then l = GLC(file_interp_type) csw_time_type(1:l) = file_interp_type(1:l) @@ -932,5 +960,4 @@ subroutine cseawater_ini() end subroutine cseawater_ini - end module ocean_emis diff --git a/src/control/camsrfexch.F90 b/src/control/camsrfexch.F90 index f978e4923c..de1ea4ce6e 100644 --- a/src/control/camsrfexch.F90 +++ b/src/control/camsrfexch.F90 @@ -61,6 +61,7 @@ module camsrfexch real(r8) :: co2prog(pcols) ! prognostic co2 real(r8) :: co2diag(pcols) ! diagnostic co2 real(r8) :: ozone(pcols) ! surface ozone concentration (mole/mole) + real(r8) :: lightning_flash_freq(pcols) ! cloud-to-ground lightning flash frequency (/min) real(r8) :: psl(pcols) real(r8) :: bcphiwet(pcols) ! wet deposition of hydrophilic black carbon real(r8) :: bcphidry(pcols) ! dry deposition of hydrophilic black carbon @@ -302,6 +303,7 @@ subroutine atm2hub_alloc( cam_out ) cam_out(c)%co2prog(:) = 0._r8 cam_out(c)%co2diag(:) = 0._r8 cam_out(c)%ozone(:) = 0._r8 + cam_out(c)%lightning_flash_freq(:) = 0._r8 cam_out(c)%psl(:) = 0._r8 cam_out(c)%bcphidry(:) = 0._r8 cam_out(c)%bcphodry(:) = 0._r8 @@ -423,7 +425,7 @@ subroutine cam_export(state,cam_out,pbuf) integer :: psl_idx integer :: prec_dp_idx, snow_dp_idx, prec_sh_idx, snow_sh_idx integer :: prec_sed_idx,snow_sed_idx,prec_pcw_idx,snow_pcw_idx - integer :: srf_ozone_idx + integer :: srf_ozone_idx, lightning_idx real(r8), pointer :: psl(:) @@ -436,6 +438,7 @@ subroutine cam_export(state,cam_out,pbuf) real(r8), pointer :: prec_pcw(:) ! total precipitation from Hack convection real(r8), pointer :: snow_pcw(:) ! snow from Hack convection real(r8), pointer :: o3_ptr(:,:), srf_o3_ptr(:) + real(r8), pointer :: lightning_ptr(:) !----------------------------------------------------------------------- lchnk = state%lchnk @@ -453,6 +456,7 @@ subroutine cam_export(state,cam_out,pbuf) prec_pcw_idx = pbuf_get_index('PREC_PCW', errcode=i) snow_pcw_idx = pbuf_get_index('SNOW_PCW', errcode=i) srf_ozone_idx = pbuf_get_index('SRFOZONE', errcode=i) + lightning_idx = pbuf_get_index('LGHT_FLASH_FREQ', errcode=i) if (prec_dp_idx > 0) then call pbuf_get_field(pbuf, prec_dp_idx, prec_dp) @@ -512,6 +516,12 @@ subroutine cam_export(state,cam_out,pbuf) cam_out%ozone(:ncol) = o3_ptr(:ncol,pver) * mwdry/mwo3 ! mole/mole endif + ! get cloud to ground lightning flash freq (/min) to export to surface models + if (lightning_idx>0) then + call pbuf_get_field(pbuf, lightning_idx, lightning_ptr) + cam_out%lightning_flash_freq(:ncol) = lightning_ptr(:ncol) + end if + ! ! Precipation and snow rates from shallow convection, deep convection and stratiform processes. ! Compute total convective and stratiform precipitation and snow rates diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90 index 356ec0f6d4..f8f182c30b 100644 --- a/src/control/runtime_opts.F90 +++ b/src/control/runtime_opts.F90 @@ -98,6 +98,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) use lunar_tides, only: lunar_tides_readnl use upper_bc, only: ubc_readnl use phys_grid_ctem, only: phys_grid_ctem_readnl + use mo_lightning, only: lightning_readnl !---------------------------Arguments----------------------------------- @@ -165,6 +166,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) call rad_data_readnl(nlfilename) call modal_aer_opt_readnl(nlfilename) call chem_readnl(nlfilename) + call lightning_readnl(nlfilename) call prescribed_volcaero_readnl(nlfilename) call prescribed_strataero_readnl(nlfilename) call solar_data_readnl(nlfilename) diff --git a/src/cpl/nuopc/atm_import_export.F90 b/src/cpl/nuopc/atm_import_export.F90 index 3c0cbba542..8c28b120fa 100644 --- a/src/cpl/nuopc/atm_import_export.F90 +++ b/src/cpl/nuopc/atm_import_export.F90 @@ -60,7 +60,8 @@ module atm_import_export integer :: drydep_nflds = -huge(1) ! number of dry deposition velocity fields lnd-> atm integer :: megan_nflds = -huge(1) ! number of MEGAN voc fields from lnd-> atm integer :: emis_nflds = -huge(1) ! number of fire emission fields from lnd-> atm - integer, public :: ndep_nflds = -huge(1) ! number of nitrogen deposition fields from atm->lnd/ocn + integer, public :: ndep_nflds = -huge(1) ! number of nitrogen deposition fields from atm->lnd/ocn + logical :: atm_provides_lightning = .false. ! cld to grnd lightning flash freq (min-1) character(*),parameter :: F01 = "('(cam_import_export) ',a,i8,2x,i8,2x,d21.14)" character(*),parameter :: F02 = "('(cam_import_export) ',a,i8,2x,i8,2x,i8,2x,d21.14)" character(*),parameter :: u_FILE_u = __FILE__ @@ -79,6 +80,7 @@ subroutine read_surface_fields_namelists() use shr_fire_emis_mod , only : shr_fire_emis_readnl use shr_carma_mod , only : shr_carma_readnl use shr_ndep_mod , only : shr_ndep_readnl + use shr_lightning_coupling_mod, only : shr_lightning_coupling_readnl character(len=*), parameter :: nl_file_name = 'drv_flds_in' @@ -88,6 +90,7 @@ subroutine read_surface_fields_namelists() call shr_megan_readnl(nl_file_name, megan_nflds) call shr_fire_emis_readnl(nl_file_name, emis_nflds) call shr_carma_readnl(nl_file_name, carma_fields) + call shr_lightning_coupling_readnl(nl_file_name, atm_provides_lightning) end subroutine read_surface_fields_namelists @@ -203,6 +206,11 @@ subroutine advertise_fields(gcomp, flds_scalar_name, rc) ! Assume that 2 fields are always sent as part of Faxa_ndep call fldlist_add(fldsFrAtm_num, fldsFrAtm, 'Faxa_ndep', ungridded_lbound=1, ungridded_ubound=2) + ! lightning flash freq + if (atm_provides_lightning) then + call fldlist_add(fldsFrAtm_num, fldsFrAtm, 'Sa_lightning') + end if + ! Now advertise above export fields if (masterproc) write(iulog,*) trim(subname)//' advertise export fields' do n = 1,fldsFrAtm_num @@ -917,6 +925,7 @@ subroutine export_fields( gcomp, model_mesh, model_clock, cam_out, rc) real(r8), pointer :: fldptr_ptem(:) , fldptr_pslv(:) real(r8), pointer :: fldptr_co2prog(:) , fldptr_co2diag(:) real(r8), pointer :: fldptr_ozone(:) + real(r8), pointer :: fldptr_lght(:) character(len=*), parameter :: subname='(atm_import_export:export_fields)' !--------------------------------------------------------------------------- @@ -1046,6 +1055,18 @@ subroutine export_fields( gcomp, model_mesh, model_clock, cam_out, rc) end do end if + call state_getfldptr(exportState, 'Sa_lightning', fldptr=fldptr_lght, exists=exists, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (exists) then + g = 1 + do c = begchunk,endchunk + do i = 1,get_ncols_p(c) + fldptr_lght(g) = cam_out(c)%lightning_flash_freq(i) ! cloud-to-ground lightning flash frequency (/min) + g = g + 1 + end do + end do + end if + call state_getfldptr(exportState, 'Sa_co2prog', fldptr=fldptr_co2prog, exists=exists, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (exists) then diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 061591f9ad..244ac339f6 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -114,6 +114,7 @@ subroutine phys_register use cam_control_mod, only: moist_physics use chemistry, only: chem_register + use mo_lightning, only: lightning_register use cloud_fraction, only: cldfrc_register use rk_stratiform, only: rk_stratiform_register use microp_driver, only: microp_driver_register @@ -267,6 +268,9 @@ subroutine phys_register ! register chemical constituents including aerosols ... call chem_register() + ! add prognostic lightning flash freq pbuf fld + call lightning_register() + ! co2 constituents call co2_register() @@ -716,6 +720,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use cam_control_mod, only: initial_run use check_energy, only: check_energy_init use chemistry, only: chem_init + use mo_lightning, only: lightning_init use prescribed_ozone, only: prescribed_ozone_init use prescribed_ghg, only: prescribed_ghg_init use prescribed_aero, only: prescribed_aero_init @@ -858,6 +863,9 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) ! Prognostic chemistry. call chem_init(phys_state,pbuf2d) + ! Lightning flash frq and NOx prod + call lightning_init( pbuf2d ) + ! Prescribed tracers call prescribed_ozone_init() call prescribed_ghg_init() @@ -1251,9 +1259,9 @@ subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d, cam_out, & ! call get_met_srf2( cam_in ) #endif - ! Set lightning production of NO + ! lightning flash freq and prod rate of NOx call t_startf ('lightning_no_prod') - call lightning_no_prod( phys_state, pbuf2d, cam_in ) + call lightning_no_prod( phys_state, pbuf2d, cam_in ) call t_stopf ('lightning_no_prod') call t_barrierf('sync_ac_physics', mpicom) diff --git a/src/physics/cam_dev/physpkg.F90 b/src/physics/cam_dev/physpkg.F90 index 968a1339d0..7ca40f8cf3 100644 --- a/src/physics/cam_dev/physpkg.F90 +++ b/src/physics/cam_dev/physpkg.F90 @@ -110,6 +110,7 @@ subroutine phys_register use cam_control_mod, only: moist_physics use chemistry, only: chem_register + use mo_lightning, only: lightning_register use cloud_fraction, only: cldfrc_register use microp_driver, only: microp_driver_register use microp_aero, only: microp_aero_register @@ -253,6 +254,9 @@ subroutine phys_register ! register chemical constituents including aerosols ... call chem_register() + ! add prognostic lightning flash freq pbuf fld + call lightning_register() + ! co2 constituents call co2_register() @@ -699,6 +703,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use cam_control_mod, only: initial_run use check_energy, only: check_energy_init use chemistry, only: chem_init + use mo_lightning, only: lightning_init use prescribed_ozone, only: prescribed_ozone_init use prescribed_ghg, only: prescribed_ghg_init use prescribed_aero, only: prescribed_aero_init @@ -832,6 +837,9 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) ! Prognostic chemistry. call chem_init(phys_state,pbuf2d) + ! Lightning flash frq and NOx prod + call lightning_init( pbuf2d ) + ! Prescribed tracers call prescribed_ozone_init() call prescribed_ghg_init() @@ -1196,9 +1204,9 @@ subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d, cam_out, & ! call get_met_srf2( cam_in ) #endif - ! Set lightning production of NO + ! lightning flash freq and prod rate of NOx call t_startf ('lightning_no_prod') - call lightning_no_prod( phys_state, pbuf2d, cam_in ) + call lightning_no_prod( phys_state, pbuf2d, cam_in ) call t_stopf ('lightning_no_prod') call t_barrierf('sync_ac_physics', mpicom)