diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000000..09130b35b9 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "generic_tracers"] + path = src/access/generic_tracers + url = https://github.com/ACCESS-NRI/GFDL-generic-tracers.git diff --git a/exp/MOM_compile.csh b/exp/MOM_compile.csh index 4dc3361f8e..3145b6b927 100755 --- a/exp/MOM_compile.csh +++ b/exp/MOM_compile.csh @@ -57,8 +57,10 @@ if ( $help ) then echo " EBM : ocean-seaice-land-atmosphere coupled model with energy balance atmosphere" echo " ACCESS-CM : ocean component of ACCESS-CM model." echo " ACCESS-OM : ocean component of ACCESS-OM model." - echo " ACCESS-ESM : ocean component of ACCESS-ESM model with CSIRO BGC (Wombat)." - echo " ACCESS-OM-BGC: ocean component of ACCESS-OM model with CSIRO BGC (Wombat)." + echo " ACCESS-ESM : ocean component of ACCESS-ESM model with support for generic tracer WOMBATlite." + echo " ACCESS-OM-BGC: ocean component of ACCESS-OM model with CSIRO BGC (Wombat). Wombat has now been" + echo " implemented as a generic tracer and is available in the ACCESS-OM model type." + echo " ACCESS-OM-BGC is retained only for legacy." echo echo "--platform followed by the platform name that has a corresponding environ file in the ../bin dir, default is gfortran" echo @@ -104,13 +106,13 @@ endif if ( $type == EBM ) then set cppDefs = ( "-Duse_netCDF -Duse_netCDF3 -Duse_libMPI -DLAND_BND_TRACERS -DOVERLOAD_C8 -DOVERLOAD_C4 -DOVERLOAD_R4" ) else if( $type == ACCESS-OM ) then - set cppDefs = ( "-Duse_netCDF -Duse_libMPI -DACCESS_OM" ) + set cppDefs = ( "-Duse_netCDF -Duse_libMPI -DACCESS_OM -DUSE_OCEAN_BGC" ) else if( $type == ACCESS-OM-BGC ) then set cppDefs = ( "-Duse_netCDF -Duse_libMPI -DACCESS_OM -DCSIRO_BGC" ) else if( $type == ACCESS-CM ) then set cppDefs = ( "-Duse_netCDF -Duse_libMPI -DACCESS_CM" ) else if( $type == ACCESS-ESM ) then - set cppDefs = ( "-Duse_netCDF -Duse_libMPI -DACCESS_CM -DCSIRO_BGC" ) + set cppDefs = ( "-Duse_netCDF -Duse_libMPI -DACCESS_CM -DUSE_OCEAN_BGC" ) endif if ( $unit_testing ) then @@ -221,12 +223,12 @@ if( $type == MOM_solo ) then set srcList = ( mom5/drivers ) set libs = "$executable:h:h/lib_ocean/lib_ocean.a $executable:h:h/lib_FMS/lib_FMS.a" else if( $type == ACCESS-CM || $type == ACCESS-ESM) then - set srcList = ( accesscm_coupler ) + set srcList = ( access/accesscm_coupler ) set includes = "-I$executable:h:h/lib_FMS -I$executable:h:h/$type/lib_ocean" set libs = "$executable:h:h/$type/lib_ocean/lib_ocean.a" setenv OASIS true else if( $type == ACCESS-OM || $type == ACCESS-OM-BGC) then - set srcList = ( accessom_coupler ) + set srcList = ( access/accessom_coupler ) set includes = "-I$executable:h:h/lib_FMS -I$executable:h:h/$type/lib_ocean" set libs = "$executable:h:h/$type/lib_ocean/lib_ocean.a $executable:h:h/lib_FMS/lib_FMS.a" else if( $type == MOM_SIS ) then @@ -254,6 +256,10 @@ else exit 1 endif +if( $type == ACCESS-OM || $type == ACCESS-ESM) then + set srcList = ( $srcList access/shared ) +endif + # Always include FMS set libs = "$libs $executable:h:h/lib_FMS/lib_FMS.a" diff --git a/exp/ocean_compile.csh b/exp/ocean_compile.csh index c547360ce8..fbfcd09ad4 100644 --- a/exp/ocean_compile.csh +++ b/exp/ocean_compile.csh @@ -6,8 +6,10 @@ set lib_name = "lib_ocean" if( $type == ACCESS-OM || $type == ACCESS-CM || $type == ACCESS-OM-BGC || $type == ACCESS-ESM) then set srcList = ( $srcList mom5/ocean_access ) - if( $type == ACCESS-OM-BGC || $type == ACCESS-ESM) then + if( $type == ACCESS-OM-BGC ) then set srcList = ( $srcList mom5/ocean_csiro_bgc ) + else if ( $type == ACCESS-OM || $type == ACCESS-ESM ) then + set srcList = ( $srcList mom5/ocean_bgc access/generic_tracers/generic_tracers access/generic_tracers/mocsy/src ) endif mkdir -p $executable:h:h/$type/$lib_name cd $executable:h:h/$type/$lib_name diff --git a/src/access/README.md b/src/access/README.md new file mode 100644 index 0000000000..837f117693 --- /dev/null +++ b/src/access/README.md @@ -0,0 +1,9 @@ +# ACCESS-specific code + +This directory includes code that is specific to ACCESS-models. + +[ACCESS-specific generic tracers](https://github.com/ACCESS-NRI/GFDL-generic-tracers) are included as a git submodule. To retrieve the code for the submodule run: + +```bash +git submodule update --init --recursive +``` diff --git a/src/accesscm_coupler/cpl_netcdf_setup.F90 b/src/access/accesscm_coupler/cpl_netcdf_setup.F90 similarity index 100% rename from src/accesscm_coupler/cpl_netcdf_setup.F90 rename to src/access/accesscm_coupler/cpl_netcdf_setup.F90 diff --git a/src/accesscm_coupler/mom_oasis3_interface.F90 b/src/access/accesscm_coupler/mom_oasis3_interface.F90 similarity index 91% rename from src/accesscm_coupler/mom_oasis3_interface.F90 rename to src/access/accesscm_coupler/mom_oasis3_interface.F90 index d34798554e..06dccc794f 100644 --- a/src/accesscm_coupler/mom_oasis3_interface.F90 +++ b/src/access/accesscm_coupler/mom_oasis3_interface.F90 @@ -94,6 +94,9 @@ module mom_oasis3_interface_mod ocean_public_type, & ocean_domain_type use time_manager_mod, only: time_type +use gtracer_flux_mod, only: set_coupler_type_data, extract_coupler_type_data +use coupler_types_mod, only: coupler_2d_bc_type, ind_pcair, ind_u10, ind_psurf, ind_csurf, ind_flux +use constants_mod, only: WTMCO2, hlv ! Timing @@ -355,7 +358,7 @@ subroutine coupler_init(Dom, Time, Time_step_coupled, Run_len, dt_cpld) mom_name_write(6)='frazil' mom_name_write(7)='dssldx' mom_name_write(8)='dssldy' - mom_name_write(9)='co2_o' + mom_name_write(9)='co2_o' ! Ocean surface pCO2 is not used by any other models mom_name_write(10)='co2fx_o' @@ -486,7 +489,7 @@ subroutine coupler_init(Dom, Time, Time_step_coupled, Run_len, dt_cpld) end subroutine coupler_init !======================================================================= -subroutine into_coupler(step, Ocean_sfc, Time, before_ocean_update) +subroutine into_coupler(step, Ocean_sfc, Ice_ocean_boundary, Time, before_ocean_update) !------------------------------------------! use ocean_operators_mod, only : GRAD_BAROTROPIC_P !GRAD_SURF_sealev @@ -496,6 +499,7 @@ subroutine into_coupler(step, Ocean_sfc, Time, before_ocean_update) implicit none type (ocean_public_type) :: Ocean_sfc +type (ice_ocean_boundary_type) :: Ice_ocean_boundary type (time_type),optional :: Time integer, intent(in) :: step @@ -561,9 +565,17 @@ subroutine into_coupler(step, Ocean_sfc, Time, before_ocean_update) case('dssldy') vtmp(iisd:iied,jjsd:jjed) = Ocean_sfc%gradient(iisd:iied,jjsd:jjed,2) case('co2_o') + ! Note, this is not actually used by the other models + ! If this is needed in the future with generic WOMBATlite, it can be calculated from the csurf + ! and alpha fields in the Ocean_sfc%fields coupler_bc_type "co2_flux" boundary condition: + ! pco2 [ppmv] = 1e6 * (co2_csurf / co2_alpha) vtmp(iisd:iied,jjsd:jjed) = Ocean_sfc%co2(iisd:iied,jjsd:jjed) case('co2fx_o') - vtmp(iisd:iied,jjsd:jjed) = Ocean_sfc%co2flux(iisd:iied,jjsd:jjed) + ! vtmp(iisd:iied,jjsd:jjed) = Ocean_sfc%co2flux(iisd:iied,jjsd:jjed) + ! Extract the flux field in the IOB%fluxes coupler_bc_type "co2_flux" boundary condition, + ! converting from [mol/m^2/s] to [kg(CO2)/m^2/s] and to positive downwards + call extract_coupler_type_data(Ice_ocean_boundary%fluxes, "co2_flux", ind_flux, vtmp, & + scale_factor=-1.e-3*WTMCO2, idim=(/iisc,iisc,iiec,iiec/), jdim=(/jjsc,jjsc,jjec,jjec/)) case DEFAULT call mpp_error(FATAL,& '==>Error from into_coupler: Unknown quantity.') @@ -611,16 +623,16 @@ subroutine into_coupler(step, Ocean_sfc, Time, before_ocean_update) end subroutine into_coupler !----------------------------------------------------------------------------------- -subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Time) +subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Atm_fields, Time) ! This is all highly user dependent. -use constants_mod, only : hlv ! 2.500e6 J/kg use auscom_ice_mod, only : chk_i2o_fields, chk_fields_period, chk_fields_start_time implicit none type (ocean_public_type) :: Ocean_sfc type (ice_ocean_boundary_type) :: Ice_ocean_boundary +type (coupler_2d_bc_type) :: Atm_fields type (time_type),optional :: Time real, dimension(isg:ieg,jsg:jeg) :: gtmp @@ -631,7 +643,7 @@ subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Time) frac_nir_dir=0.5*0.57, frac_nir_dif=0.5*0.57 ! shortwave partitioning character*80 :: fname = 'fields_i2o_in_ocn.nc' - integer :: ncid,currstep,ll,ilout + integer :: ncid,currstep,ll,ilout,n data currstep/0/ save currstep @@ -723,7 +735,11 @@ subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Time) case('wfiform') Ice_ocean_boundary%wfiform(iisc:iiec,jjsc:jjec) = vwork(iisc:iiec,jjsc:jjec) case('co2_io') - Ice_ocean_boundary%co2(iisc:iiec,jjsc:jjec) = vwork(iisc:iiec,jjsc:jjec) + ! Ice_ocean_boundary%co2(iisc:iiec,jjsc:jjec) = vwork(iisc:iiec,jjsc:jjec) + ! Set the pcair field in the Atm_fields coupler_bc_type "co2_flux" boundary condition, + ! converting from [ppmv] to [mol/mol] + call set_coupler_type_data(vwork, "co2_flux", ind_pcair, Atm_fields, & + scale_factor=1.e-6, idim=(/iisc,iisc,iiec,iiec/), jdim=(/jjsc,jjsc,jjec,jjec/)) case('wnd_io') Ice_ocean_boundary%wnd(iisc:iiec,jjsc:jjec) = vwork(iisc:iiec,jjsc:jjec) !20171024: 2 more i2o fields: water and heat fluxes due to land ice discharge into ocean @@ -752,6 +768,18 @@ subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Time) if(jf .ne. 1) call mpp_clock_end(id_oasis_recv1) enddo !jf + + ! Set the u10 and psurf fields in the Atm_fields coupler_bc_types + do n = 1, Atm_fields%num_bcs + if ((Atm_fields%bc(n)%flux_type .eq. 'air_sea_gas_flux_generic') .or. & + (Atm_fields%bc(n)%flux_type .eq. 'air_sea_gas_flux')) then + call set_coupler_type_data(ice_ocean_boundary%wnd, Atm_fields%bc(n)%name, ind_u10, & + Atm_fields, idim=(/iisc,iisc,iiec,iiec/), jdim=(/jjsc,jjsc,jjec,jjec/)) + call set_coupler_type_data(ice_ocean_boundary%p, Atm_fields%bc(n)%name, ind_psurf, & + Atm_fields, idim=(/iisc,iisc,iiec,iiec/), jdim=(/jjsc,jjsc,jjec,jjec/)) + endif + enddo + call mpp_clock_end(id_oasis_recv) if (chk_i2o_fields .and. (mod(step, chk_fields_period) == 0) .and. (step >= chk_fields_start_time) .and. (mpp_pe() == mpp_root_pe())) then @@ -761,13 +789,14 @@ subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Time) end subroutine from_coupler !----------------------------------------------------------------------------------- -subroutine write_coupler_restart(step,Ocean_sfc,write_restart) +subroutine write_coupler_restart(step,Ocean_sfc,Ice_ocean_boundary,write_restart) use auscom_ice_mod, only : auscom_ice_heatflux_new logical, intent(in) :: write_restart integer, intent(in) :: step type (ocean_public_type) :: Ocean_sfc +type (ice_ocean_boundary_type) :: Ice_ocean_boundary integer :: ncid,ll,ilout real, dimension(iisd:iied,jjsd:jjed) :: vtmp @@ -794,7 +823,13 @@ subroutine write_coupler_restart(step,Ocean_sfc,write_restart) case('dssldy'); vtmp = Ocean_sfc%gradient(iisd:iied,jjsd:jjed,2); fld_ice='ssly_i' case('frazil'); vtmp = Ocean_sfc%frazil(iisd:iied,jjsd:jjed); fld_ice='pfmice_i' case('co2_o'); vtmp = Ocean_sfc%co2(iisd:iied,jjsd:jjed); fld_ice='co2_oi' - case('co2fx_o'); vtmp = Ocean_sfc%co2flux(iisd:iied,jjsd:jjed); fld_ice='co2fx_oi' + case('co2fx_o') + ! vtmp = Ocean_sfc%co2flux(iisd:iied,jjsd:jjed); fld_ice='co2fx_oi' + ! Extract the flux field in the IOB%fluxes coupler_bc_type "co2_flux" boundary condition, + ! converting from [mol/m^2/s] to [kg(CO2)/m^2/s] and to positive downwards + call extract_coupler_type_data(Ice_ocean_boundary%fluxes, "co2_flux", ind_flux, vtmp, & + scale_factor=-1.e-3*WTMCO2, idim=(/iisc,iisc,iiec,iiec/), jdim=(/jjsc,jjsc,jjec,jjec/)) + fld_ice='co2fx_oi' end select if (parallel_coupling) then diff --git a/src/accesscm_coupler/ocean_solo.F90 b/src/access/accesscm_coupler/ocean_solo.F90 similarity index 89% rename from src/accesscm_coupler/ocean_solo.F90 rename to src/access/accesscm_coupler/ocean_solo.F90 index 9944be7f30..0fd42996ac 100644 --- a/src/accesscm_coupler/ocean_solo.F90 +++ b/src/access/accesscm_coupler/ocean_solo.F90 @@ -154,11 +154,16 @@ program main use auscom_ice_parameters_mod, only: redsea_gulfbay_sfix, do_sfix_now, sfix_hours, int_sec + use coupler_types_mod, only: coupler_2d_bc_type, coupler_type_data_override, coupler_type_send_data + use gtracer_flux_mod, only: flux_exchange_init, atmos_ocean_fluxes_calc + use gtracer_flux_mod, only: gas_fields_restore, gas_fields_restart + implicit none type (ocean_public_type) :: Ocean_sfc type (ocean_state_type), pointer :: Ocean_state - type(ice_ocean_boundary_type), target :: Ice_ocean_boundary + type(ice_ocean_boundary_type), target :: Ice_ocean_boundary + type(coupler_2d_bc_type), target :: Atm_fields ! define some time types type(time_type) :: Time_init ! initial time for experiment @@ -171,7 +176,8 @@ program main type(time_type) :: Time_restart type(time_type) :: Time_restart_current type(time_type) :: Time_last_sfix - type(time_type) :: Time_sfix + type(time_type) :: Time_sfix + type(time_type) :: Time_next integer :: sfix_seconds character(len=17) :: calendar = 'julian' @@ -393,56 +399,13 @@ program main call data_override_init(Ocean_domain_in = Ocean_sfc%domain) override_clock = mpp_clock_id('Override', flags=flags,grain=CLOCK_COMPONENT) - - call mpp_get_compute_domain(Ocean_sfc%domain, isc, iec, jsc, jec) - - allocate ( Ice_ocean_boundary% u_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% v_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% t_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% q_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% salt_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% lw_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_vis_dir (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_vis_dif (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_nir_dir (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_nir_dif (isc:iec,jsc:jec), & - Ice_ocean_boundary% lprec (isc:iec,jsc:jec), & - Ice_ocean_boundary% fprec (isc:iec,jsc:jec), & - Ice_ocean_boundary% runoff (isc:iec,jsc:jec), & - Ice_ocean_boundary% calving (isc:iec,jsc:jec), & - Ice_ocean_boundary% p (isc:iec,jsc:jec), & - Ice_ocean_boundary% aice(isc:iec,jsc:jec), & - Ice_ocean_boundary% mh_flux(isc:iec,jsc:jec), & - Ice_ocean_boundary% wfimelt(isc:iec,jsc:jec), & - Ice_ocean_boundary% wfiform(isc:iec,jsc:jec)) - allocate ( Ice_ocean_boundary%co2(isc:iec,jsc:jec), & - Ice_ocean_boundary%wnd(isc:iec,jsc:jec), & - Ice_ocean_boundary%licefw(isc:iec,jsc:jec), & - Ice_ocean_boundary%liceht(isc:iec,jsc:jec)) - - Ice_ocean_boundary%u_flux = 0.0 - Ice_ocean_boundary%v_flux = 0.0 - Ice_ocean_boundary%t_flux = 0.0 - Ice_ocean_boundary%q_flux = 0.0 - Ice_ocean_boundary%salt_flux = 0.0 - Ice_ocean_boundary%lw_flux = 0.0 - Ice_ocean_boundary%sw_flux_vis_dir = 0.0 - Ice_ocean_boundary%sw_flux_vis_dif = 0.0 - Ice_ocean_boundary%sw_flux_nir_dir = 0.0 - Ice_ocean_boundary%sw_flux_nir_dif = 0.0 - Ice_ocean_boundary%lprec = 0.0 - Ice_ocean_boundary%fprec = 0.0 - Ice_ocean_boundary%runoff = 0.0 - Ice_ocean_boundary%calving = 0.0 - Ice_ocean_boundary%p = 0.0 - Ice_ocean_boundary%aice = 0.0 - Ice_ocean_boundary%mh_flux = 0.0 - Ice_ocean_boundary% wfimelt = 0.0 - Ice_ocean_boundary% wfiform = 0.0 - Ice_ocean_boundary%co2 = 0.0 - Ice_ocean_boundary%wnd = 0.0 - Ice_ocean_boundary%licefw = 0.0 - Ice_ocean_boundary%liceht = 0.0 + + ! Initialise the boundary values, including initialising and setting boundary values + ! in FMS coupler types + call flux_exchange_init(Time, Ocean_sfc, Ocean_state, Ice_ocean_boundary, Atm_fields) + + ! Restore ocean FMS coupler type fields from restart file + call gas_fields_restore(Ocean_sfc) coupler_init_clock = mpp_clock_id('OASIS init', grain=CLOCK_COMPONENT) call mpp_clock_begin(coupler_init_clock) @@ -452,12 +415,25 @@ program main ! loop over the coupled calls do nc=1, num_cpld_calls - call external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, nc, dt_cpld ) + call external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, Atm_fields, nc, dt_cpld ) + ! Potentially override fields from the data_table call mpp_clock_begin(override_clock) - call ice_ocn_bnd_from_data(Ice_ocean_boundary) + Time_next = Time + Time_step_coupled + call coupler_type_data_override('OCN', Atm_fields, Time_next) + call ice_ocn_bnd_from_data(Ice_ocean_boundary, Time_next) call mpp_clock_end(override_clock) + ! Calculate the extra tracer fluxes + call mpp_get_compute_domain(Ocean_sfc%domain, isc, iec, jsc, jec) + call atmos_ocean_fluxes_calc(Atm_fields, Ocean_sfc%fields, Ice_ocean_boundary%fluxes, & + Ice_ocean_boundary%aice, isc, iec, jsc, jec) + + ! Send FMS coupler type diagnostics + call coupler_type_send_data(Ice_ocean_boundary%fluxes, Time_next) + call coupler_type_send_data(Ocean_sfc%fields, Time_next) + call coupler_type_send_data(Atm_fields, Time_next) + if (debug_this_module) then call write_boundary_chksums(Ice_ocean_boundary) endif @@ -473,24 +449,25 @@ program main call update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, Time, Time_step_coupled) - Time = Time + Time_step_coupled + Time = Time_next if( Time >= Time_restart ) then - Time_restart_current = Time - Time_restart = increment_date(Time, restart_interval(1), restart_interval(2), & + Time_restart_current = Time + Time_restart = increment_date(Time, restart_interval(1), restart_interval(2), & restart_interval(3), restart_interval(4), restart_interval(5), restart_interval(6) ) - timestamp = date_to_string(time_restart_current) + timestamp = date_to_string(time_restart_current) write(stdoutunit,*) '=> NOTE from program ocean_solo: intermediate restart file is written and ', & trim(timestamp),' is appended as prefix to each restart file name' call ocean_model_restart(Ocean_state, timestamp) call ocean_solo_restart(Time, Time_restart_current, timestamp) + call gas_fields_restart(Ocean_sfc, timestamp) end if call external_coupler_sbc_after(Ice_ocean_boundary, Ocean_sfc, nc, dt_cpld ) enddo - call external_coupler_restart( dt_cpld, num_cpld_calls, Ocean_sfc) + call external_coupler_restart( dt_cpld, num_cpld_calls, Ocean_sfc, Ice_ocean_boundary) ! close some of the main components call ocean_model_end(Ocean_sfc, Ocean_state, Time) @@ -500,8 +477,9 @@ program main ! need to reset pelist before calling mpp_clock_end ! call mpp_set_current_pelist() - ! write restart file + ! write restart files call ocean_solo_restart(Time_end, Time_restart_current) + call gas_fields_restart(Ocean_sfc) call fms_io_exit @@ -558,13 +536,14 @@ subroutine ocean_solo_restart(Time_run, Time_res, time_stamp) end subroutine ocean_solo_restart !==================================================================== -! get forcing data from data_overide -subroutine ice_ocn_bnd_from_data(x) +! get forcing data from data_override +subroutine ice_ocn_bnd_from_data(x, Time_next) type (ice_ocean_boundary_type) :: x type(time_type) :: Time_next - Time_next = Time + Time_step_coupled + integer :: m, n + call data_override('OCN', 't_flux', x%t_flux , Time_next) call data_override('OCN', 'u_flux', x%u_flux , Time_next) call data_override('OCN', 'v_flux', x%v_flux , Time_next) @@ -582,12 +561,21 @@ subroutine ice_ocn_bnd_from_data(x) call data_override('OCN', 'p', x%p , Time_next) call data_override('OCN', 'aice', x%aice , Time_next) call data_override('OCN', 'mh_flux', x%mh_flux , Time_next) + + ! Overriding ice_ocean_boundary%fluxes here avoids unnecessary calculation + ! of overridden fluxes. However, we cannot use coupler_type_data_override + ! here since it does not set the override flag on overridden fields + do n = 1, x%fluxes%num_bcs + do m = 1, x%fluxes%bc(n)%num_fields + call data_override('OCN', x%fluxes%bc(n)%field(m)%name, & + x%fluxes%bc(n)%field(m)%values, Time_next, & + override=x%fluxes%bc(n)%field(m)%override) + enddo + enddo call mpp_sync() end subroutine ice_ocn_bnd_from_data - - ! Here we provide some hooks for calling an interface between the OASIS3 coupler and MOM. ! The mom_oasis3_interface module is NOT general and it is expected that the user will ! heavily modify it depending on the coupling strategy. @@ -619,7 +607,7 @@ subroutine external_coupler_sbc_init(Dom, dt_cpld, Run_len) call coupler_init(Dom, dt_cpld=dt_cpld, Run_len=Run_len) end subroutine external_coupler_sbc_init -subroutine external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, nsteps, dt_cpld ) +subroutine external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, Atm_fields, nsteps, dt_cpld ) ! Perform transfers before ocean time stepping ! May need special tratment on first call. @@ -628,6 +616,7 @@ subroutine external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, nsteps, dt implicit none type (ice_ocean_boundary_type), intent(INOUT) :: Ice_ocean_boundary type (ocean_public_type) , intent(INOUT) :: Ocean_sfc + type (coupler_2d_bc_type), intent(INOUT) :: Atm_fields integer , intent(IN) :: nsteps, dt_cpld integer :: rtimestep ! Receive timestep @@ -635,8 +624,8 @@ subroutine external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, nsteps, dt rtimestep = (nsteps-1) * dt_cpld ! runtime in this run segment! stimestep = rtimestep - call from_coupler( rtimestep, Ocean_sfc, Ice_ocean_boundary ) - call into_coupler( stimestep, Ocean_sfc, before_ocean_update = .true.) + call from_coupler( rtimestep, Ocean_sfc, Ice_ocean_boundary, Atm_fields ) + call into_coupler( stimestep, Ocean_sfc, Ice_ocean_boundary, before_ocean_update = .true.) end subroutine external_coupler_sbc_before subroutine external_coupler_sbc_after(Ice_ocean_boundary, Ocean_sfc, nsteps, dt_cpld ) @@ -652,19 +641,21 @@ subroutine external_coupler_sbc_after(Ice_ocean_boundary, Ocean_sfc, nsteps, dt_ integer :: stimestep ! Send timestep stimestep = nsteps * dt_cpld ! runtime in this run segment! - if (stimestep < num_cpld_calls*dt_cpld) call into_coupler(stimestep, Ocean_sfc, before_ocean_update = .false.) + if (stimestep < num_cpld_calls*dt_cpld) call into_coupler(stimestep, Ocean_sfc, & + Ice_ocean_boundary, before_ocean_update = .false.) end subroutine external_coupler_sbc_after -subroutine external_coupler_restart( dt_cpld, num_cpld_calls, Ocean_sfc) +subroutine external_coupler_restart( dt_cpld, num_cpld_calls, Ocean_sfc, Ice_ocean_boundary) !Clean up as appropriate and write a restart use mom_oasis3_interface_mod, only : write_coupler_restart implicit none integer, intent(in) :: dt_cpld, num_cpld_calls integer :: timestep type (ocean_public_type) :: Ocean_sfc + type (ice_ocean_boundary_type) :: Ice_ocean_boundary timestep = num_cpld_calls * dt_cpld - call write_coupler_restart(timestep, Ocean_sfc, write_restart=.true.) + call write_coupler_restart(timestep, Ocean_sfc, Ice_ocean_boundary, write_restart=.true.) end subroutine external_coupler_restart subroutine external_coupler_exit diff --git a/src/accessom_coupler/cpl_netcdf_setup.F90 b/src/access/accessom_coupler/cpl_netcdf_setup.F90 similarity index 100% rename from src/accessom_coupler/cpl_netcdf_setup.F90 rename to src/access/accessom_coupler/cpl_netcdf_setup.F90 diff --git a/src/accessom_coupler/mom_oasis3_interface.F90 b/src/access/accessom_coupler/mom_oasis3_interface.F90 similarity index 97% rename from src/accessom_coupler/mom_oasis3_interface.F90 rename to src/access/accessom_coupler/mom_oasis3_interface.F90 index 939ecf8319..058434930a 100644 --- a/src/accessom_coupler/mom_oasis3_interface.F90 +++ b/src/access/accessom_coupler/mom_oasis3_interface.F90 @@ -106,6 +106,7 @@ module mom_oasis3_interface_mod ocean_public_type, & ocean_domain_type use time_manager_mod, only: time_type, get_time +use gtracer_flux_mod, only: set_coupler_type_data ! Timing @@ -657,16 +658,18 @@ subroutine into_coupler(step, Ocean_sfc, Time, before_ocean_update) end subroutine into_coupler !----------------------------------------------------------------------------------- -subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Time) +subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Atm_fields, Time) ! This is all highly user dependent. use constants_mod, only : hlv ! 2.500e6 J/kg use auscom_ice_mod, only : chk_i2o_fields, chk_fields_period, chk_fields_start_time +use coupler_types_mod, only: coupler_2d_bc_type, ind_u10, ind_psurf implicit none type (ocean_public_type) :: Ocean_sfc type (ice_ocean_boundary_type) :: Ice_ocean_boundary +type (coupler_2d_bc_type) :: Atm_fields type (time_type),optional :: Time real, dimension(isg:ieg,jsg:jeg) :: gtmp @@ -674,7 +677,7 @@ subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Time) integer, intent(in) :: step character*80 :: fname = 'fields_i2o_in_ocn.nc' - integer :: ncid,currstep,ll,ilout + integer :: ncid,currstep,ll,ilout,n data currstep/0/ save currstep @@ -798,6 +801,18 @@ subroutine from_coupler(step,Ocean_sfc,Ice_ocean_boundary, Time) if(jf .ne. 1) call mpp_clock_end(id_oasis_recv1) enddo !jf + + ! Set the u10 and psurf fields in the Atm_fields coupler_bc_types + do n = 1, Atm_fields%num_bcs + if ((Atm_fields%bc(n)%flux_type .eq. 'air_sea_gas_flux_generic') .or. & + (Atm_fields%bc(n)%flux_type .eq. 'air_sea_gas_flux')) then + call set_coupler_type_data(ice_ocean_boundary%wnd, Atm_fields%bc(n)%name, ind_u10, & + Atm_fields, idim=(/iisc,iisc,iiec,iiec/), jdim=(/jjsc,jjsc,jjec,jjec/)) + call set_coupler_type_data(ice_ocean_boundary%p, Atm_fields%bc(n)%name, ind_psurf, & + Atm_fields, idim=(/iisc,iisc,iiec,iiec/), jdim=(/jjsc,jjsc,jjec,jjec/)) + endif + enddo + call mpp_clock_end(id_oasis_recv) if (chk_i2o_fields .and. (mod(step, chk_fields_period) == 0) .and. (step >= chk_fields_start_time) .and. (mpp_pe() == mpp_root_pe())) then diff --git a/src/accessom_coupler/ocean_solo.F90 b/src/access/accessom_coupler/ocean_solo.F90 similarity index 89% rename from src/accessom_coupler/ocean_solo.F90 rename to src/access/accessom_coupler/ocean_solo.F90 index 838759cde0..5baf2de772 100644 --- a/src/accessom_coupler/ocean_solo.F90 +++ b/src/access/accessom_coupler/ocean_solo.F90 @@ -111,11 +111,16 @@ program main use accessom2_mod, only : accessom2_type => accessom2 use coupler_mod, only : coupler_type => coupler + use coupler_types_mod, only: coupler_2d_bc_type, coupler_type_data_override, coupler_type_send_data + use gtracer_flux_mod, only: flux_exchange_init, atmos_ocean_fluxes_calc + use gtracer_flux_mod, only: gas_fields_restore, gas_fields_restart + implicit none type (ocean_public_type) :: Ocean_sfc type (ocean_state_type), pointer :: Ocean_state type(ice_ocean_boundary_type), target :: Ice_ocean_boundary + type(coupler_2d_bc_type), target :: Atm_fields type(accessom2_type) :: accessom2 type(coupler_type) :: coupler @@ -130,7 +135,8 @@ program main type(time_type) :: Time_restart type(time_type) :: Time_restart_current type(time_type) :: Time_last_sfix - type(time_type) :: Time_sfix + type(time_type) :: Time_sfix + type(time_type) :: Time_next integer :: sfix_seconds character(len=17) :: calendar = 'julian' @@ -388,61 +394,14 @@ program main call data_override_init(Ocean_domain_in = Ocean_sfc%domain) override_clock = mpp_clock_id('Override', flags=flags,grain=CLOCK_COMPONENT) + + ! Initialise the boundary values, including initialising and setting boundary values + ! in FMS coupler types + call flux_exchange_init(Time, Ocean_sfc, Ocean_state, Ice_ocean_boundary, Atm_fields) + + ! Restore ocean FMS coupler type fields from restart file + call gas_fields_restore(Ocean_sfc) - call mpp_get_compute_domain(Ocean_sfc%domain, isc, iec, jsc, jec) - - allocate ( Ice_ocean_boundary% u_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% v_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% t_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% q_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% salt_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% lw_flux (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_vis_dir (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_vis_dif (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_nir_dir (isc:iec,jsc:jec), & - Ice_ocean_boundary% sw_flux_nir_dif (isc:iec,jsc:jec), & - Ice_ocean_boundary% lprec (isc:iec,jsc:jec), & - Ice_ocean_boundary% fprec (isc:iec,jsc:jec), & - Ice_ocean_boundary% runoff (isc:iec,jsc:jec), & - Ice_ocean_boundary% calving (isc:iec,jsc:jec), & - Ice_ocean_boundary% p (isc:iec,jsc:jec), & - Ice_ocean_boundary% aice(isc:iec,jsc:jec), & - Ice_ocean_boundary% mh_flux(isc:iec,jsc:jec), & - Ice_ocean_boundary% wfimelt(isc:iec,jsc:jec), & - Ice_ocean_boundary% wfiform(isc:iec,jsc:jec), & - Ice_ocean_boundary% licefw(isc:iec,jsc:jec), & - Ice_ocean_boundary% liceht(isc:iec,jsc:jec), & - Ice_ocean_boundary%wnd(isc:iec,jsc:jec)) -#if defined(ACCESS_OM) && defined(CSIRO_BGC) - allocate ( Ice_ocean_boundary%iof_nit(isc:iec,jsc:jec), & - Ice_ocean_boundary%iof_alg(isc:iec,jsc:jec)) -#endif - Ice_ocean_boundary%u_flux = 0.0 - Ice_ocean_boundary%v_flux = 0.0 - Ice_ocean_boundary%t_flux = 0.0 - Ice_ocean_boundary%q_flux = 0.0 - Ice_ocean_boundary%salt_flux = 0.0 - Ice_ocean_boundary%lw_flux = 0.0 - Ice_ocean_boundary%sw_flux_vis_dir = 0.0 - Ice_ocean_boundary%sw_flux_vis_dif = 0.0 - Ice_ocean_boundary%sw_flux_nir_dir = 0.0 - Ice_ocean_boundary%sw_flux_nir_dif = 0.0 - Ice_ocean_boundary%lprec = 0.0 - Ice_ocean_boundary%fprec = 0.0 - Ice_ocean_boundary%runoff = 0.0 - Ice_ocean_boundary%calving = 0.0 - Ice_ocean_boundary%p = 0.0 - Ice_ocean_boundary%aice = 0.0 - Ice_ocean_boundary%mh_flux = 0.0 - Ice_ocean_boundary% wfimelt = 0.0 - Ice_ocean_boundary% wfiform = 0.0 - Ice_ocean_boundary%licefw = 0.0 - Ice_ocean_boundary%liceht = 0.0 - Ice_ocean_boundary%wnd = 0.0 -#if defined(ACCESS_OM) && defined(CSIRO_BGC) - Ice_ocean_boundary%iof_nit = 0.0 - Ice_ocean_boundary%iof_alg = 0.0 -#endif coupler_init_clock = mpp_clock_id('OASIS init', grain=CLOCK_COMPONENT) call mpp_clock_begin(coupler_init_clock) call external_coupler_sbc_init(Ocean_sfc%domain, dt_cpld, Run_len, & @@ -454,12 +413,25 @@ program main call mpp_clock_begin(main_clock) do nc=1, num_cpld_calls - call external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, nc, dt_cpld) + call external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, Atm_fields, nc, dt_cpld) + ! Potentially override fields from the data_table call mpp_clock_begin(override_clock) - call ice_ocn_bnd_from_data(Ice_ocean_boundary) + Time_next = Time + Time_step_coupled + call coupler_type_data_override('OCN', Atm_fields, Time_next) + call ice_ocn_bnd_from_data(Ice_ocean_boundary, Time_next) call mpp_clock_end(override_clock) + ! Calculate the extra tracer fluxes + call mpp_get_compute_domain(Ocean_sfc%domain, isc, iec, jsc, jec) + call atmos_ocean_fluxes_calc(Atm_fields, Ocean_sfc%fields, Ice_ocean_boundary%fluxes, & + Ice_ocean_boundary%aice, isc, iec, jsc, jec) + + ! Send FMS coupler type diagnostics + call coupler_type_send_data(Ice_ocean_boundary%fluxes, Time_next) + call coupler_type_send_data(Ocean_sfc%fields, Time_next) + call coupler_type_send_data(Atm_fields, Time_next) + if (debug_this_module) then call write_boundary_chksums(Ice_ocean_boundary) endif @@ -475,20 +447,21 @@ program main call update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, Time, Time_step_coupled) - Time = Time + Time_step_coupled + Time = Time_next if ( mpp_pe() == mpp_root_pe() ) then call accessom2%progress_date(int(dt_cpld)) endif if( Time >= Time_restart ) then - Time_restart_current = Time - Time_restart = increment_date(Time, restart_interval(1), restart_interval(2), & + Time_restart_current = Time + Time_restart = increment_date(Time, restart_interval(1), restart_interval(2), & restart_interval(3), restart_interval(4), restart_interval(5), restart_interval(6) ) - timestamp = date_to_string(time_restart_current) + timestamp = date_to_string(time_restart_current) write(stdoutunit,*) '=> NOTE from program ocean_solo: intermediate restart file is written and ', & trim(timestamp),' is appended as prefix to each restart file name' call ocean_model_restart(Ocean_state, timestamp) call ocean_solo_restart(Time, Time_restart_current, timestamp) + call gas_fields_restart(Ocean_sfc, timestamp) end if call external_coupler_sbc_after(Ice_ocean_boundary, Ocean_sfc, nc, dt_cpld ) @@ -508,8 +481,9 @@ program main ! need to reset pelist before calling mpp_clock_end ! call mpp_set_current_pelist() - ! write restart file + ! write restart files call ocean_solo_restart(Time_end, Time_restart_current) + call gas_fields_restart(Ocean_sfc) call fms_io_exit @@ -574,12 +548,13 @@ end subroutine ocean_solo_restart !==================================================================== ! get forcing data from data_overide -subroutine ice_ocn_bnd_from_data(x) +subroutine ice_ocn_bnd_from_data(x, Time_next) type (ice_ocean_boundary_type) :: x type(time_type) :: Time_next - Time_next = Time + Time_step_coupled + integer :: m, n + call data_override('OCN', 't_flux', x%t_flux , Time_next) call data_override('OCN', 'u_flux', x%u_flux , Time_next) call data_override('OCN', 'v_flux', x%v_flux , Time_next) @@ -597,12 +572,21 @@ subroutine ice_ocn_bnd_from_data(x) call data_override('OCN', 'p', x%p , Time_next) call data_override('OCN', 'aice', x%aice , Time_next) call data_override('OCN', 'mh_flux', x%mh_flux , Time_next) + + ! Overriding ice_ocean_boundary%fluxes here avoids unnecessary calculation + ! of overridden fluxes. However, we cannot use coupler_type_data_override + ! here since it does not set the override flag on overridden fields + do n = 1, x%fluxes%num_bcs + do m = 1, x%fluxes%bc(n)%num_fields + call data_override('OCN', x%fluxes%bc(n)%field(m)%name, & + x%fluxes%bc(n)%field(m)%values, Time_next, & + override=x%fluxes%bc(n)%field(m)%override) + enddo + enddo call mpp_sync() end subroutine ice_ocn_bnd_from_data - - ! Here we provide some hooks for calling an interface between the OASIS3 coupler and MOM. ! The mom_oasis3_interface module is NOT general and it is expected that the user will ! heavily modify it depending on the coupling strategy. @@ -626,7 +610,7 @@ subroutine external_coupler_sbc_init(Dom, dt_cpld, Run_len, & coupling_field_timesteps=coupling_field_timesteps) end subroutine external_coupler_sbc_init -subroutine external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, nsteps, dt_cpld ) +subroutine external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, Atm_fields, nsteps, dt_cpld ) ! Perform transfers before ocean time stepping ! May need special tratment on first call. @@ -635,6 +619,7 @@ subroutine external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, nsteps, dt implicit none type (ice_ocean_boundary_type), intent(INOUT) :: Ice_ocean_boundary type (ocean_public_type) , intent(INOUT) :: Ocean_sfc + type (coupler_2d_bc_type), intent(INOUT) :: Atm_fields integer , intent(IN) :: nsteps, dt_cpld integer :: rtimestep ! Receive timestep @@ -642,7 +627,7 @@ subroutine external_coupler_sbc_before(Ice_ocean_boundary, Ocean_sfc, nsteps, dt rtimestep = (nsteps-1) * dt_cpld ! runtime in this run segment! stimestep = rtimestep - call from_coupler( rtimestep, Ocean_sfc, Ice_ocean_boundary ) + call from_coupler( rtimestep, Ocean_sfc, Ice_ocean_boundary, Atm_fields ) call into_coupler( stimestep, Ocean_sfc, before_ocean_update = .true.) end subroutine external_coupler_sbc_before diff --git a/src/access/generic_tracers b/src/access/generic_tracers new file mode 160000 index 0000000000..6dcadbbaef --- /dev/null +++ b/src/access/generic_tracers @@ -0,0 +1 @@ +Subproject commit 6dcadbbaef22c8c7b20f7ed5cf75df37dbade3b8 diff --git a/src/access/shared/gtracer_flux.F90 b/src/access/shared/gtracer_flux.F90 new file mode 100644 index 0000000000..4e264ac344 --- /dev/null +++ b/src/access/shared/gtracer_flux.F90 @@ -0,0 +1,738 @@ +!> Contains routines for handling FMS coupler_bc_type tracer flux structures +!! and calculating fluxes when using generic tracers + +module gtracer_flux_mod + + use time_manager_mod, only: time_type + use coupler_types_mod, only: coupler_1d_bc_type, coupler_2d_bc_type + use coupler_types_mod, only: coupler_type_spawn, coupler_type_set_diags + use coupler_types_mod, only: coupler_type_set_data, coupler_type_extract_data + use coupler_types_mod, only: coupler_type_register_restarts, coupler_type_restore_state + use coupler_types_mod, only: ind_flux, ind_deltap, ind_kw + use coupler_types_mod, only: ind_pcair, ind_u10, ind_psurf + use coupler_types_mod, only: ind_alpha, ind_csurf, ind_sc_no + use coupler_types_mod, only: ind_runoff, ind_deposition + use coupler_types_mod, only: atmos_ocean_type_fluxes_init => coupler_types_init + use ocean_model_mod, only: ocean_public_type, ocean_state_type + use ocean_model_mod, only: ocean_model_init_sfc, ocean_model_flux_init + use ocean_types_mod, only: ice_ocean_boundary_type + use atmos_ocean_fluxes_mod, only: atmos_ocean_fluxes_init + use mpp_domains_mod, only: mpp_get_compute_domain + use field_manager_mod, only: fm_field_name_len, fm_type_name_len, fm_loop_over_list, fm_change_list + use fm_util_mod, only: fm_util_get_real_array + use fms_io_mod, only: restart_file_type, save_restart + use constants_mod, only: wtmair, rdgas, vonkarm + use mpp_mod, only: mpp_error, FATAL + + implicit none; private + + ! Public member functions + public :: flux_exchange_init + public :: gas_fields_restore + public :: gas_fields_restart + public :: set_coupler_type_data + public :: extract_coupler_type_data + public :: atmos_ocean_fluxes_calc + + character(len=*), parameter :: mod_name = 'gtracer_flux_mod' + real, parameter :: epsln=1.0e-30 + + contains + + !> \brief Initialise the boundary values, including initialising and setting boundary + !! values in FMS coupler types. Based loosely on the flux_exchange_init subroutine in + !! flux_exchange_mod (coupler/flux_exchange.F90) + subroutine flux_exchange_init (Time, Ocean, Ocean_state, Ice_ocean_boundary, atm_fields) + type(time_type), intent(in) :: Time + type(ocean_public_type), intent(inout) :: Ocean + type(ocean_state_type), pointer :: Ocean_state + type(ice_ocean_boundary_type), intent(inout) :: Ice_ocean_boundary + type(coupler_2d_bc_type), intent(inout) :: atm_fields + + ! local variables + integer :: isc,iec,jsc,jec + type(coupler_1d_bc_type) :: gas_fields_atm ! tracer fields in atm + type(coupler_1d_bc_type) :: gas_fields_ocn ! tracer fields atop the ocean + type(coupler_1d_bc_type) :: gas_fluxes ! tracer fluxes between the atm and ocean + + ! Initialise FMS coupler types + call atmos_ocean_type_fluxes_init( ) + call ocean_model_flux_init(Ocean_state) + call atmos_ocean_fluxes_init(gas_fluxes, gas_fields_atm, gas_fields_ocn) + + call mpp_get_compute_domain(Ocean%domain, isc, iec, jsc, jec) + + ! Allocate Ice_ocean_boundary fields + allocate( & + Ice_ocean_boundary%u_flux(isc:iec,jsc:jec), & + Ice_ocean_boundary%v_flux(isc:iec,jsc:jec), & + Ice_ocean_boundary%t_flux(isc:iec,jsc:jec), & + Ice_ocean_boundary%q_flux(isc:iec,jsc:jec), & + Ice_ocean_boundary%salt_flux(isc:iec,jsc:jec), & + Ice_ocean_boundary%lw_flux(isc:iec,jsc:jec), & + Ice_ocean_boundary%sw_flux_vis_dir(isc:iec,jsc:jec), & + Ice_ocean_boundary%sw_flux_vis_dif(isc:iec,jsc:jec), & + Ice_ocean_boundary%sw_flux_nir_dir(isc:iec,jsc:jec), & + Ice_ocean_boundary%sw_flux_nir_dif(isc:iec,jsc:jec), & + Ice_ocean_boundary%lprec(isc:iec,jsc:jec), & + Ice_ocean_boundary%fprec(isc:iec,jsc:jec), & + Ice_ocean_boundary%runoff(isc:iec,jsc:jec), & + Ice_ocean_boundary%calving(isc:iec,jsc:jec), & + Ice_ocean_boundary%p(isc:iec,jsc:jec), & + Ice_ocean_boundary%aice(isc:iec,jsc:jec), & + Ice_ocean_boundary%mh_flux(isc:iec,jsc:jec), & + Ice_ocean_boundary%wfimelt(isc:iec,jsc:jec), & + Ice_ocean_boundary%wfiform(isc:iec,jsc:jec), & + Ice_ocean_boundary%licefw(isc:iec,jsc:jec), & + Ice_ocean_boundary%liceht(isc:iec,jsc:jec), & + Ice_ocean_boundary%wnd(isc:iec,jsc:jec)) +#if defined(ACCESS_CM) + allocate(Ice_ocean_boundary%co2(isc:iec,jsc:jec)) +#endif + + Ice_ocean_boundary%u_flux = 0.0 + Ice_ocean_boundary%v_flux = 0.0 + Ice_ocean_boundary%t_flux = 0.0 + Ice_ocean_boundary%q_flux = 0.0 + Ice_ocean_boundary%salt_flux = 0.0 + Ice_ocean_boundary%lw_flux = 0.0 + Ice_ocean_boundary%sw_flux_vis_dir = 0.0 + Ice_ocean_boundary%sw_flux_vis_dif = 0.0 + Ice_ocean_boundary%sw_flux_nir_dir = 0.0 + Ice_ocean_boundary%sw_flux_nir_dif = 0.0 + Ice_ocean_boundary%lprec = 0.0 + Ice_ocean_boundary%fprec = 0.0 + Ice_ocean_boundary%runoff = 0.0 + Ice_ocean_boundary%calving = 0.0 + Ice_ocean_boundary%p = 0.0 + Ice_ocean_boundary%aice = 0.0 + Ice_ocean_boundary%mh_flux = 0.0 + Ice_ocean_boundary% wfimelt = 0.0 + Ice_ocean_boundary% wfiform = 0.0 + Ice_ocean_boundary%licefw = 0.0 + Ice_ocean_boundary%liceht = 0.0 + Ice_ocean_boundary%wnd = 0.0 +#if defined(ACCESS_CM) + Ice_ocean_boundary%co2 = 0.0 +#endif + + ! Spawn 2D Ocean%fields FMS coupler type from 1D gas_fields_ocn + call coupler_type_spawn(gas_fields_ocn, Ocean%fields, (/isc,isc,iec,iec/), & + (/jsc,jsc,jec,jec/), suffix = '_ocn') + call coupler_type_set_diags(Ocean%fields, "ocean_sfc", Ocean%axes(1:2), Time) + + ! This sets the boundary values in Ocean%fields + call ocean_model_init_sfc(Ocean_state, Ocean) + + ! Spawn 2D Ice_ocean_boundary%fluxes FMS coupler type from 1D gas_fluxes + ! Annoyingly, spawning doesn't copy param array, so add manually + call coupler_type_spawn(gas_fluxes, Ice_ocean_boundary%fluxes, (/isc,isc,iec,iec/), & + (/jsc,jsc,jec,jec/), suffix='_ice_ocn') + call add_gas_fluxes_param(Ice_ocean_boundary%fluxes) + call coupler_type_set_diags(Ice_ocean_boundary%fluxes, "ocean_flux", Ocean%axes(1:2), Time) + + ! Spawn 2D atm_fields FMS coupler type from 1D gas_fields_atm + call coupler_type_spawn(gas_fields_atm, atm_fields, (/isc,isc,iec,iec/), & + (/jsc,jsc,jec,jec/), suffix='_atm') + call coupler_type_set_diags(atm_fields, "atmos_sfc", Ocean%axes(1:2), Time) + + end subroutine flux_exchange_init + + !> \brief Restore FMS coupler_bc_type state from ocean restart file. Based loosely on + !! code from the coupler_init subroutine in coupler_main (coupler/coupler_main.F90 and + !! FMScoupler) + subroutine gas_fields_restore(Ocean) + type(ocean_public_type), intent(inout) :: Ocean + + ! local variables + type(restart_file_type), dimension(:), pointer :: ocn_bc_restart => NULL() + integer :: num_ocn_bc_restart + + call coupler_type_register_restarts(Ocean%fields, ocn_bc_restart, num_ocn_bc_restart, & + Ocean%domain, ocean_restart=.true.) + + call coupler_type_restore_state(Ocean%fields, test_by_field=.true.) + + end subroutine gas_fields_restore + + !> \brief Write ocean restart file for FMS coupler_bc_type state. Based loosely on + !! code from the coupler_restart subroutine in coupler_main (coupler/coupler_main.F90) + subroutine gas_fields_restart(Ocean, timestamp) + type(ocean_public_type), intent(inout) :: Ocean + character(len=64), optional :: timestamp + + ! local variables + type(restart_file_type), dimension(:), pointer :: ocn_bc_restart => NULL() + integer :: num_ocn_bc_restart + integer :: n + + call coupler_type_register_restarts(Ocean%fields, ocn_bc_restart, num_ocn_bc_restart, & + Ocean%domain, ocean_restart=.true.) + + do n = 1, num_ocn_bc_restart + call save_restart(ocn_bc_restart(n), timestamp) + enddo + + end subroutine gas_fields_restart + + !> Retrieve param array from field_manager and add to FMS coupler_bc_type. This is + !! needed because the coupler_type_spawn routine does not copy the param array into + !! the spawned type. This routine is based on the atmos_ocean_fluxes_init subroutine + !! in atmos_ocean_fluxes_mod (shared/coupler/atmos_ocean_fluxes.F90) + !! https://github.com/NOAA-GFDL/FMS/blob/7f585284f1487c0629f8075be350385e6e75dd2f/coupler/atmos_ocean_fluxes.F90#L448 + subroutine add_gas_fluxes_param(gas_fluxes) + type(coupler_2d_bc_type), intent(inout) :: gas_fluxes !< FMS coupler_bc_type to add param to + + ! local variables + integer :: n + character(len=fm_field_name_len) :: name + character(len=fm_type_name_len) :: typ + integer :: ind + + character(len=*), parameter :: sub_name = 'add_gas_fluxes_param' + character(len=*), parameter :: error_header =& + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + + n = 0 + do while (fm_loop_over_list('/coupler_mod/fluxes', name, typ, ind)) + if (typ .ne. 'list') then + call mpp_error(FATAL, trim(error_header) // ' ' // trim(name) // ' is not a list') + endif + + n = n + 1 + + if (.not. fm_change_list('/coupler_mod/fluxes/' // trim(name))) then + call mpp_error(FATAL, trim(error_header) // ' Problem changing to ' // trim(name)) + endif + + if (gas_fluxes%bc(n)%name .eq. name) then + gas_fluxes%bc(n)%param => fm_util_get_real_array('param') + else + call mpp_error(FATAL, trim(error_header) // ' Problem setting param array pointer') + endif + enddo + end subroutine add_gas_fluxes_param + + !> Set single field by name in a 2D FMS coupler type from a two-dimensional array. + !! This is a convenience wrapper on coupler_type_set_data routine in + !! coupler_types_mod (shared/coupler/coupler_types.F90) + subroutine set_coupler_type_data(array_in, name, field_index, var, scale_factor, halo_size, idim, jdim, override) + real, dimension(1:,1:), intent(in) :: array_in !< The source array for the field; its size must match the + !! size of the data being copied unless idim and jdim are + !! supplied. + character(*), intent(in) :: name !< The name of the boundary condition to set + integer, intent(in) :: field_index !< The index of the field in the boundary condition that + !! is being set. + type(coupler_2d_bc_type), intent(inout) :: var !< BC_type structure with the data to set + real, optional, intent(in) :: scale_factor !< A scaling factor for the data that is being added + integer, optional, intent(in) :: halo_size !< The extent of the halo to copy; 0 by default + integer, dimension(4), optional, intent(in) :: idim !< The data and computational domain extents of the first + !! dimension of the output array in a non-decreasing list + integer, dimension(4), optional, intent(in) :: jdim !< The data and computational domain extents of the second + !! dimension of the output array in a non-decreasing list + logical, optional, intent(in) :: override !< If true, set the override flag for the field being set + + integer :: n + logical :: override_flag = .false. + logical :: fmatch = .false. + + character(len=*), parameter :: sub_name = 'set_coupler_type_data' + character(len=*), parameter :: error_header =& + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + + if (present(override)) override_flag = override + + do n = 1, var%num_bcs + if ( var%bc(n)%name .eq. trim(name) ) then + call coupler_type_set_data(array_in, n, field_index, var, & + scale_factor=scale_factor, halo_size=halo_size, idim=idim, jdim=jdim) + if ( override_flag ) then + var%bc(n)%field(field_index)%override = override_flag + endif + fmatch = .true. + exit + endif + enddo + if (.not. fmatch) then + call mpp_error(FATAL, trim(error_header) // ' No coupler type BC called '//trim(name)) + endif + end subroutine set_coupler_type_data + + !> Extract single field by name from a 2D FMS coupler type into a two-dimensional array. + !! This is a convenience wrapper on coupler_type_extract_data routine in + !! coupler_types_mod (shared/coupler/coupler_types.F90) + subroutine extract_coupler_type_data(var_in, name, field_index, array_out, scale_factor, halo_size, idim, jdim) + type(coupler_2d_bc_type), intent(in) :: var_in !< BC_type structure with the data to extract + character(*), intent(in) :: name !< The name of the boundary condition to extract + integer, intent(in) :: field_index !< The index of the field in the boundary condition that + !! is being extracted. + real, dimension(1:,1:), intent(out) :: array_out !< The recipient array for the field; its size + !! must match the size of the data being copied + !! unless idim and jdim are supplied. + real, optional, intent(in) :: scale_factor !< A scaling factor for the data that is being extracted + integer, optional, intent(in) :: halo_size !< The extent of the halo to copy; 0 by default + integer, dimension(4), optional, intent(in) :: idim !< The data and computational domain extents of the first + !! dimension of the output array in a non-decreasing list + integer, dimension(4), optional, intent(in) :: jdim !< The data and computational domain extents of the second + !! dimension of the output array in a non-decreasing list + + integer :: n + logical :: fmatch = .false. + + character(len=*), parameter :: sub_name = 'extract_coupler_type_data' + character(len=*), parameter :: error_header =& + '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + + do n = 1, var_in%num_bcs + if ( var_in%bc(n)%name .eq. trim(name) ) then + call coupler_type_extract_data(var_in, n, field_index, array_out, & + scale_factor=scale_factor, halo_size=halo_size, idim=idim, jdim=jdim) + fmatch = .true. + exit + endif + enddo + if (.not. fmatch) then + call mpp_error(FATAL, trim(error_header) // ' No coupler type BC called '//trim(name)) + endif + end subroutine extract_coupler_type_data + + !> \brief Calculate the FMS coupler_bc_type ocean tracer fluxes. Units should be mol/m^2/s. + !! Upward flux is positive. + !! + !! This routine is based on the atmos_ocean_fluxes_calc subroutine in atmos_ocean_fluxes_mod + !! (shared/coupler/atmos_ocean_fluxes.F90 and FMScoupler) modified in the following ways: + !! - Operate on 2D inputs, rather than 1D + !! - Add calculation for 'air_sea_deposition' from atmos_ocean_dep_fluxes_calc subroutine + !! - Multiply fluxes by ice_fraction input, rather than masking based on seawater input + !! - Make tsurf input optional, as it is only used by a few implementations + !! - Use ind_runoff rather than ind_deposition in runoff flux calculation (note, their + !! values are equal) + !! - Rename gas_fields_ice to gas_fields_ocn + subroutine atmos_ocean_fluxes_calc(gas_fields_atm, gas_fields_ocn, gas_fluxes,& + ice_fraction, isc, iec, jsc, jec, tsurf, ustar, cd_m) + type(coupler_2d_bc_type), intent(in) :: gas_fields_atm ! fields in atm + !< Structure containing atmospheric surface variables that are used in the calculation + !! of the atmosphere-ocean tracer fluxes. + type(coupler_2d_bc_type), intent(in) :: gas_fields_ocn ! fields atop the ocean + !< Structure containing ocean surface variables that are used in the calculation of the + !! atmosphere-ocean tracer fluxes. + type(coupler_2d_bc_type), intent(inout) :: gas_fluxes ! fluxes between the atm and ocean + !< Structure containing the gas fluxes between the atmosphere and the ocean and + !! parameters related to the calculation of these fluxes. + real, intent(in) :: ice_fraction(isc:iec,jsc:jec) !< sea ice fraction + integer, intent(in) :: isc !< The start i-index of cell centers within + !! the computational domain + integer, intent(in) :: iec !< The end i-index of cell centers within the + !! computational domain + integer, intent(in) :: jsc !< The start j-index of cell centers within + !! the computational domain + integer, intent(in) :: jec !< The end j-index of cell centers within the + !! computational domain + real, intent(in), optional :: tsurf(isc:iec,jsc:jec) !< surface temperature + real, intent(in), optional :: ustar(isc:iec,jsc:jec) !< friction velocity, not + !! used + real, intent(in), optional :: cd_m (isc:iec,jsc:jec) !< drag coefficient, not + !! used + + ! local variables + character(len=*), parameter :: sub_name = 'atmos_ocean_fluxes_calc' + character(len=*), parameter :: error_header =& + & '==>Error from ' // trim(mod_name) // '(' // trim(sub_name) // '):' + real, parameter :: permeg=1.0e-6 + + integer :: n + integer :: i + integer :: j + real, dimension(:,:), allocatable :: kw + real, dimension(:,:), allocatable :: cair + character(len=128) :: error_string + + ! Return if no fluxes to be calculated + if (gas_fluxes%num_bcs .le. 0) return + + if (.not. associated(gas_fluxes%bc)) then + if (gas_fluxes%num_bcs .ne. 0) then + call mpp_error(FATAL, trim(error_header) // ' Number of gas fluxes not zero') + else + return + endif + endif + + do n = 1, gas_fluxes%num_bcs + ! only do calculations if the flux has not been overridden + if ( .not. gas_fluxes%bc(n)%field(ind_flux)%override) then + if (gas_fluxes%bc(n)%flux_type .eq. 'air_sea_gas_flux_generic') then + if (.not. allocated(kw)) then + allocate( kw(isc:iec,jsc:jec) ) + allocate ( cair(isc:iec,jsc:jec) ) + elseif ((size(kw(:,:), dim=1) .ne. iec-isc+1) .or. (size(kw(:,:), dim=2) .ne. jec-jsc+1)) then + call mpp_error(FATAL, trim(error_header) // ' Sizes of flux fields do not match') + endif + + if (gas_fluxes%bc(n)%implementation .eq. 'ocmip2') then + do j = jsc,jec + do i = isc,iec + gas_fluxes%bc(n)%field(ind_kw)%values(i,j) =& + & (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) * & + & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j)**2 + cair(i,j) = & + gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * & + gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) * & + gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(2) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) =& + & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *& + & sqrt(660. / (gas_fields_ocn%bc(n)%field(ind_sc_no)%values(i,j) + epsln)) *& + & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) + gas_fluxes%bc(n)%field(ind_deltap)%values(i,j) =& + & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) / & + (gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * permeg + epsln) + enddo + enddo + elseif (gas_fluxes%bc(n)%implementation .eq. 'duce') then + if (.not. present(tsurf)) then + call mpp_error(FATAL, trim(error_header) // ' Implementation ' //& + trim(gas_fluxes%bc(n)%implementation) // ' for ' // trim(gas_fluxes%bc(n)%name) //& + ' requires input tsurf') + endif + do j = jsc,jec + do i = isc,iec + gas_fluxes%bc(n)%field(ind_kw)%values(i,j) = & + & (1 - ice_fraction(i,j)) * gas_fields_atm%bc(n)%field(ind_u10)%values(i,j) /& + & (770.+45.*gas_fluxes%bc(n)%param(1)**(1./3.)) *& + & 101325./(rdgas*wtmair*1e-3*tsurf(i,j) *& + & max(gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j),epsln)) + !alpha: mol/m3/atm + cair(i,j) = & + gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * & + gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) * & + gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * 9.86923e-6 + cair(i,j) = max(cair(i,j),0.) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) =& + & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *& + & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j)) + gas_fluxes%bc(n)%field(ind_deltap)%values(i,j) =& + & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j)) /& + & (gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * permeg + epsln) + enddo + enddo + elseif (gas_fluxes%bc(n)%implementation .eq. 'johnson') then + if (.not. present(tsurf)) then + call mpp_error(FATAL, trim(error_header) // ' Implementation ' //& + trim(gas_fluxes%bc(n)%implementation) // ' for ' // trim(gas_fluxes%bc(n)%name) //& + ' requires input tsurf') + endif + !f1p: not sure how to pass salinity. For now, just force at 35. + do j = jsc,jec + do i = isc,iec + !calc_kw(tk,p,u10,h,vb,mw,sc_w,ustar,cd_m) + gas_fluxes%bc(n)%field(ind_kw)%values(i,j) =& + & (1 - ice_fraction(i,j)) * calc_kw(tsurf(i,j),& + & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j),& + & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j),& + & 101325./(rdgas*wtmair*1e-3*tsurf(i,j)*max(gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j),epsln)),& + & gas_fluxes%bc(n)%param(2),& + & gas_fluxes%bc(n)%param(1),& + & gas_fields_ocn%bc(n)%field(ind_sc_no)%values(i,j)) + cair(i,j) =& + & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * 9.86923e-6 + cair(i,j) = max(cair(i,j),0.) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) =& + & gas_fluxes%bc(n)%field(ind_kw)%values(i,j) *& + & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j)) + gas_fluxes%bc(n)%field(ind_deltap)%values(i,j) =& + & (max(gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j),0.) - cair(i,j)) /& + & (gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) * permeg + epsln) + enddo + enddo + else + call mpp_error(FATAL, ' Unknown implementation (' //& + & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name)) + endif + elseif (gas_fluxes%bc(n)%flux_type .eq. 'air_sea_gas_flux') then + if (.not. allocated(kw)) then + allocate( kw(isc:iec,jsc:jec) ) + allocate ( cair(isc:iec,jsc:jec) ) + elseif ((size(kw(:,:), dim=1) .ne. iec-isc+1) .or. (size(kw(:,:), dim=2) .ne. jec-jsc+1)) then + call mpp_error(FATAL, trim(error_header) // ' Sizes of flux fields do not match') + endif + + if (gas_fluxes%bc(n)%implementation .eq. 'ocmip2_data') then + do j = jsc,jec + do i = isc,iec + kw(i,j) = (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) *& + & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j) + cair(i,j) =& + & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(2) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = kw(i,j) *& + & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) + enddo + enddo + elseif (gas_fluxes%bc(n)%implementation .eq. 'ocmip2') then + do j = jsc,jec + do i = isc,iec + kw(i,j) = (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) *& + & gas_fields_atm%bc(n)%field(ind_u10)%values(i,j)**2 + cair(i,j) =& + & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(2) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = kw(i,j) *& + & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) + enddo + enddo + elseif (gas_fluxes%bc(n)%implementation .eq. 'linear') then + do j = jsc,jec + do i = isc,iec + kw(i,j) = (1 - ice_fraction(i,j)) * gas_fluxes%bc(n)%param(1) *& + & max(0.0, gas_fields_atm%bc(n)%field(ind_u10)%values(i,j) - gas_fluxes%bc(n)%param(2)) + cair(i,j) =& + & gas_fields_ocn%bc(n)%field(ind_alpha)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_pcair)%values(i,j) *& + & gas_fields_atm%bc(n)%field(ind_psurf)%values(i,j) * gas_fluxes%bc(n)%param(3) + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = kw(i,j) *& + & (gas_fields_ocn%bc(n)%field(ind_csurf)%values(i,j) - cair(i,j)) + enddo + enddo + else + call mpp_error(FATAL, ' Unknown implementation (' //& + & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name)) + endif + elseif (gas_fluxes%bc(n)%flux_type .eq. 'air_sea_deposition') then + if (gas_fluxes%bc(n)%param(1) .le. 0.0) then + write (error_string, '(1pe10.3)') gas_fluxes%bc(n)%param(1) + call mpp_error(FATAL, 'Bad parameter (' // trim(error_string) //& + & ') for air_sea_deposition for ' // trim(gas_fluxes%bc(n)%name)) + endif + + if (gas_fluxes%bc(n)%implementation .eq. 'dry') then + do j = jsc,jec + do i = isc,iec + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = (1 - ice_fraction(i,j)) *& + gas_fields_atm%bc(n)%field(ind_deposition)%values(i,j) / gas_fluxes%bc(n)%param(1) + enddo + enddo + elseif (gas_fluxes%bc(n)%implementation .eq. 'wet') then + do j = jsc,jec + do i = isc,iec + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = (1 - ice_fraction(i,j)) *& + gas_fields_atm%bc(n)%field(ind_deposition)%values(i,j) / gas_fluxes%bc(n)%param(1) + enddo + enddo + else + call mpp_error(FATAL, 'Unknown implementation (' //& + & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name)) + endif + elseif (gas_fluxes%bc(n)%flux_type .eq. 'land_sea_runoff') then + if (gas_fluxes%bc(n)%param(1) .le. 0.0) then + write (error_string, '(1pe10.3)') gas_fluxes%bc(n)%param(1) + call mpp_error(FATAL, ' Bad parameter (' // trim(error_string) //& + & ') for land_sea_runoff for ' // trim(gas_fluxes%bc(n)%name)) + endif + + if (gas_fluxes%bc(n)%implementation .eq. 'river') then + do j = jsc,jec + do i = isc,iec + gas_fluxes%bc(n)%field(ind_flux)%values(i,j) = (1 - ice_fraction(i,j)) *& + & gas_fields_atm%bc(n)%field(ind_runoff)%values(i,j) /& + & gas_fluxes%bc(n)%param(1) + enddo + enddo + else + call mpp_error(FATAL, ' Unknown implementation (' //& + & trim(gas_fluxes%bc(n)%implementation) // ') for ' // trim(gas_fluxes%bc(n)%name)) + endif + else + call mpp_error(FATAL, ' Unknown flux_type (' // trim(gas_fluxes%bc(n)%flux_type) //& + & ') for ' // trim(gas_fluxes%bc(n)%name)) + endif + endif + enddo + + if (allocated(kw)) then + deallocate(kw) + deallocate(cair) + endif + end subroutine atmos_ocean_fluxes_calc + + !> Calculate \f$k_w\f$ + !! + !! Taken from Johnson, Ocean Science, 2010. (http://doi.org/10.5194/os-6-913-2010) + !! + !! Uses equations defined in Liss[1974], + !! \f[ + !! F = K_g(c_g - H C_l) = K_l(c_g/H - C_l) + !! \f] + !! where \f$c_g\f$ and \f$C_l\f$ are the bulk gas and liquid concentrations, \f$H\f$ + !! is the Henry's law constant (\f$H = c_{sg}/C_{sl}\f$, where \f$c_{sg}\f$ is the + !! equilibrium concentration in gas phase (\f$g/cm^3\f$ of air) and \f$C_{sl}\f$ is the + !! equilibrium concentration of unionised dissolved gas in liquid phase (\f$g/cm^3\f$ + !! of water)), + !! \f[ + !! 1/K_g = 1/k_g + H/k_l + !! \f] + !! and + !! \f[ + !! 1/K_l = 1/k_l + 1/{Hk_g} + !! \f] + !! where \f$k_g\f$ and \f$k_l\f$ are the exchange constants for the gas and liquid + !! phases, respectively. + !! + !! This routine was copied from FMScoupler at + !! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d387153064644325c96a5e9e2935139d5e3c/full/atmos_ocean_fluxes_calc.F90 + real function calc_kw(tk, p, u10, h, vb, mw, sc_w, ustar, cd_m) + real, intent(in) :: tk !< temperature at surface in kelvin + real, intent(in) :: p !< pressure at surface in pa + real, intent(in) :: u10 !< wind speed at 10m above the surface in m/s + real, intent(in) :: h !< Henry's law constant (\f$H=c_sg/C_sl\f$) (unitless) + real, intent(in) :: vb !< Molar volume + real, intent(in) :: mw !< molecular weight (g/mol) + real, intent(in) :: sc_w + real, intent(in), optional :: ustar !< Friction velocity (m/s). If not provided, + !! ustar = \f$u_{10} \sqrt{C_D}\f$. + real, intent(in), optional :: cd_m !< Drag coefficient (\f$C_D\f$). Used only if + !! ustar is provided. + !! If ustar is not provided, + !! cd_m = \f$6.1 \times 10^{-4} + 0.63 \times 10^{-4} *u_10\f$ + + real :: ra,rl,tc + + tc = tk-273.15 + ra = 1./max(h*calc_ka(tc,p,mw,vb,u10,ustar,cd_m),epsln) + rl = 1./max(calc_kl(tc,u10,sc_w),epsln) + calc_kw = 1./max(ra+rl,epsln) + end function calc_kw + + !> Calculate \f$k_a\f$ + !! + !! See calc_kw + !! + !! This routine was copied from FMScoupler at + !! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d387153064644325c96a5e9e2935139d5e3c/full/atmos_ocean_fluxes_calc.F90 + real function calc_ka(t, p, mw, vb, u10, ustar, cd_m) + real, intent(in) :: t !< temperature at surface in C + real, intent(in) :: p !< pressure at surface in pa + real, intent(in) :: mw !< molecular weight (g/mol) + real, intent(in) :: vb !< molar volume + real, intent(in) :: u10 !< wind speed at 10m above the surface in m/s + real, intent(in), optional :: ustar !< Friction velocity (m/s). If not provided, + !! ustar = \f$u_{10} \sqrt{C_D}\f$. + real, intent(in), optional :: cd_m !< Drag coefficient (\f$C_D\f$). Used only if + !! ustar is provided. + !! If ustar is not provided, + !! cd_m = \f$6.1 \times 10^{-4} + 0.63 \times 10^{-4} *u_10\f$ + + real :: sc + real :: ustar_t, cd_m_t + + if (.not. present(ustar)) then + !drag coefficient + cd_m_t = 6.1e-4 +0.63e-4*u10 + !friction velocity + ustar_t = u10*sqrt(cd_m_t) + else + cd_m_t = cd_m + ustar_t = ustar + end if + sc = schmidt_g(t,p,mw,vb) + calc_ka = 1e-3+ustar_t/(13.3*sqrt(sc)+1/sqrt(cd_m_t)-5.+log(sc)/(2.*vonkarm)) + end function calc_ka + + !> Calculate \f$k_l\f$ + !! + !! See calc_kw, and Nightingale, Global Biogeochemical Cycles, 2000 + !! (https://doi.org/10.1029/1999GB900091) + !! + !! This routine was copied from FMScoupler at + !! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d387153064644325c96a5e9e2935139d5e3c/full/atmos_ocean_fluxes_calc.F90 + real function calc_kl(t, v, sc) + real, intent(in) :: t !< temperature at surface in C + real, intent(in) :: v !< wind speed at surface in m/s + real, intent(in) :: sc + + calc_kl = (((0.222*v**2)+0.333*v)*(max(sc,epsln)/600.)**(-0.5))/(100.*3600.) + end function calc_kl + + !> Schmidt number of the gas in air + !! + !! This routine was copied from FMScoupler at + !! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d387153064644325c96a5e9e2935139d5e3c/full/atmos_ocean_fluxes_calc.F90 + real function schmidt_g(t, p, mw, vb) + real, intent(in) :: t !< temperature at surface in C + real, intent(in) :: p !< pressure at surface in pa + real, intent(in) :: mw !< molecular weight (g/mol) + real, intent(in) :: vb !< molar volume + + real :: d,v + + d = d_air(t,p,mw,vb) + v = v_air(t) + schmidt_g = v / d + end function schmidt_g + + !> From Fuller, Industrial & Engineering Chemistry (https://doi.org/10.1021/ie50677a007) + !! + !! This routine was copied from FMScoupler at + !! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d387153064644325c96a5e9e2935139d5e3c/full/atmos_ocean_fluxes_calc.F90 + real function d_air(t, p, mw, vb) + real, intent(in) :: t !< temperature in c + real, intent(in) :: p !< pressure in pa + real, intent(in) :: mw !< molecular weight (g/mol) + real, intent(in) :: vb !< diffusion coefficient (\f$cm3/mol\f$) + + real, parameter :: ma = 28.97d0 !< molecular weight air in g/mol + real, parameter :: va = 20.1d0 !< diffusion volume for air (\f$cm^3/mol\f$) + + real :: pa + + ! convert p to atm + pa = 9.8692d-6*p + d_air = 1d-3 *& + & (t+273.15d0)**(1.75d0)*sqrt(1d0/ma + 1d0/mw)/(pa*(va**(1d0/3d0)+vb**(1d0/3d0))**2d0) + ! d_air is in cm2/s convert to m2/s + d_air = d_air * 1d-4 + end function d_air + + !> kinematic viscosity in air + !! + !! This routine was copied from FMScoupler at + !! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d387153064644325c96a5e9e2935139d5e3c/full/atmos_ocean_fluxes_calc.F90 + real function p_air(t) + real, intent(in) :: t + + real, parameter :: sd_0 = 1.293393662d0,& + & sd_1 = -5.538444326d-3,& + & sd_2 = 3.860201577d-5,& + & sd_3 = -5.2536065d-7 + p_air = sd_0+(sd_1*t)+(sd_2*t**2)+(sd_3*t**3) + end function p_air + + !> Kinematic viscosity in air (\f$m^2/s\f$ + !! + !! This routine was copied from FMScoupler at + !! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d387153064644325c96a5e9e2935139d5e3c/full/atmos_ocean_fluxes_calc.F90 + real function v_air(t) + real, intent(in) :: t !< temperature in C + v_air = n_air(t)/p_air(t) + end function v_air + + !> dynamic viscosity in air + !! + !! This routine was copied from FMScoupler at + !! https://github.com/NOAA-GFDL/FMScoupler/blob/6442d387153064644325c96a5e9e2935139d5e3c/full/atmos_ocean_fluxes_calc.F90 + real function n_air(t) + real, intent(in) :: t !< temperature in C + + real, parameter :: sv_0 = 1.715747771d-5,& + & sv_1 = 4.722402075d-8,& + & sv_2 = -3.663027156d-10,& + & sv_3 = 1.873236686d-12,& + & sv_4 = -8.050218737d-14 + ! in n.s/m^2 (pa.s) + n_air = sv_0+(sv_1*t)+(sv_2*t**2)+(sv_3*t**3)+(sv_4*t**4) + end function n_air + +end module gtracer_flux_mod diff --git a/src/mom5/ocean_bgc/ocean_generic_tracer.F90 b/src/mom5/ocean_bgc/ocean_generic_tracer.F90 index 737345bd23..e4cf03afa6 100644 --- a/src/mom5/ocean_bgc/ocean_generic_tracer.F90 +++ b/src/mom5/ocean_bgc/ocean_generic_tracer.F90 @@ -27,10 +27,10 @@ module ocean_generic_mod use field_manager_mod, only: fm_get_index,fm_string_len, fm_new_value use generic_tracer, only: generic_tracer_init, generic_tracer_source, generic_tracer_update_from_bottom - use generic_tracer, only: generic_tracer_coupler_get, generic_tracer_coupler_set, generic_tracer_register_diag + use generic_tracer, only: generic_tracer_coupler_accumulate, generic_tracer_coupler_set, generic_tracer_register_diag use generic_tracer, only: generic_tracer_end, generic_tracer_get_list, do_generic_tracer, generic_tracer_register use generic_tracer, only: generic_tracer_coupler_zero, generic_tracer_vertdiff_G, generic_tracer_vertdiff_M - use generic_tracer, only: generic_tracer_diag + use generic_tracer, only: generic_tracer_diag, generic_tracer_update_from_coupler use g_tracer_utils, only: g_tracer_get_name,g_tracer_get_alias,g_tracer_set_values,g_tracer_get_common use g_tracer_utils, only: g_tracer_get_next,g_tracer_type,g_tracer_is_prog,g_tracer_flux_init @@ -51,6 +51,7 @@ module ocean_generic_mod public ocean_generic_sum_sfc public ocean_generic_zero_sfc public ocean_generic_sbc + public ocean_generic_sbc_adjust public ocean_generic_init public ocean_generic_flux_init public ocean_generic_column_physics @@ -336,9 +337,14 @@ subroutine ocean_generic_sbc(Ice_ocean_boundary_fluxes,Disd,Djsd, T_prog, runoff character(len=fm_string_len) :: g_tracer_name character(len=fm_string_len), parameter :: sub_name = 'update_generic_tracer_sbc' - !Extract the tracer surface fields from coupler - call generic_tracer_coupler_get(Ice_ocean_boundary_fluxes) + !Extract the tracer surface fields from coupler + ! dts: change to use generic_tracer_coupler_accumulate with a weight of 1. for consistency + ! with MOM6. Note this means that the generic_?_update_from_coupler method is no longer + ! called here. Instead this will be called later in ocean_generic_sbc_adjust when flux + ! adjustments are available + call generic_tracer_coupler_accumulate(Ice_ocean_boundary_fluxes, 1.) + ! dts: Is this still needed here, since this is now also done in ocean_generic_sbc_adjust? !Update T_prog fields from generic tracer fields ! !Get the tracer list @@ -394,6 +400,87 @@ subroutine ocean_generic_sbc(Ice_ocean_boundary_fluxes,Disd,Djsd, T_prog, runoff enddo end subroutine ocean_generic_sbc + ! + ! + ! Adjust tracer surface boundary conditions + ! + ! + ! This subroutine applies any adjustments to the generic tracers surface boundary conditions, + ! for example, adjustment of surface tracer fluxes due to salinity restoring. + ! + ! + ! + subroutine ocean_generic_sbc_adjust(Disd, Djsd, T_prog, salt_flux_added, runoff) + integer, intent(in) :: Disd, Djsd + type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog + real, intent(in), dimension(Disd:,Djsd:) :: salt_flux_added + real, intent(in), dimension(Disd:,Djsd:) :: runoff + + type(g_tracer_type), pointer :: g_tracer_list,g_tracer,g_tracer_next + integer :: g_tracer_index + character(len=fm_string_len) :: g_tracer_name + character(len=fm_string_len), parameter :: sub_name = 'ocean_generic_sbc_adjust' + + ! Adjust tracer fields via the generic_?_update_from_coupler method + call generic_tracer_update_from_coupler(Disd, Djsd, salt_flux_added) + + !Update T_prog fields from generic tracer fields + ! + !Get the tracer list + call generic_tracer_get_list(g_tracer_list) + if(.NOT. associated(g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//& + ": No tracer in the list.") + !For each tracer name get its T_prog index and get its flux fields + g_tracer=>g_tracer_list + do + if(g_tracer_is_prog(g_tracer)) then + call g_tracer_get_alias(g_tracer,g_tracer_name) + g_tracer_index = fm_get_index(trim('/ocean_mod/prog_tracers/'//g_tracer_name)) + if (g_tracer_index .le. 0) & + call mpp_error(FATAL,trim(sub_name) // ' Could not get the index for '//g_tracer_name) + + if (_ALLOCATED(g_tracer%stf) )& + call g_tracer_get_values(g_tracer,g_tracer_name,'stf', T_prog(g_tracer_index)%stf, Disd,Djsd) + + if (_ALLOCATED(g_tracer%btf) )& + call g_tracer_get_values(g_tracer,g_tracer_name,'btf', T_prog(g_tracer_index)%btf, Disd,Djsd) + + !If the tracer has runoff fill in the T_prog(n)%trunoff and T_prog(n)%runoff_tracer_flux + if(_ALLOCATED(g_tracer%trunoff)) then + !Fill in T_prog(n)%trunoff + + call g_tracer_get_values(g_tracer,g_tracer_name,'trunoff',T_prog(g_tracer_index)%trunoff,Disd,Djsd) + + !Fill in T_prog(n)%runoff_tracer_flux + T_prog(g_tracer_index)%runoff_tracer_flux = T_prog(g_tracer_index)%trunoff * runoff + + !Set g_tracer%runoff_tracer_flux + call g_tracer_set_values(g_tracer,g_tracer_name,'runoff_tracer_flux',T_prog(g_tracer_index)%runoff_tracer_flux,Disd,Djsd) + ! + !Fill in T_prog(n)%triver in MOM + !Note: This is done so that MOM can apply the river fluxes through setting either + ! the runoff and calving fluxes (when discharge_combine_runoff_calve=.false.) + ! or + ! the total river concentration (when discharge_combine_runoff_calve=.true.) + ! + !Assume zero calving flux for the generic tracers. + !T_prog(g_tracer_index)%tcalving = 0 !MOM default + T_prog(g_tracer_index)%triver = T_prog(g_tracer_index)%trunoff + + endif + + endif + + !traverse the linked list till hit NULL + call g_tracer_get_next(g_tracer, g_tracer_next) + if(.NOT. associated(g_tracer_next)) exit + g_tracer=>g_tracer_next + + enddo + end subroutine ocean_generic_sbc_adjust + ! ! ! Column physics for generic tracers. diff --git a/src/mom5/ocean_core/ocean_model.F90 b/src/mom5/ocean_core/ocean_model.F90 index 1f98edfa25..176773177e 100644 --- a/src/mom5/ocean_core/ocean_model.F90 +++ b/src/mom5/ocean_core/ocean_model.F90 @@ -1605,7 +1605,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, Ocean_state, Ocean_sfc, & ! compute "flux adjustments" (e.g., surface tracer restoring, flux correction) call mpp_clock_begin(id_flux_adjust) call flux_adjust(Time, T_diag(1:num_diag_tracers), Dens, Ext_mode, & - T_prog(1:num_prog_tracers), Velocity, river, melt, pme) + T_prog(1:num_prog_tracers), Velocity, river, melt, pme, runoff) call mpp_clock_end(id_flux_adjust) diff --git a/src/mom5/ocean_core/ocean_sbc.F90 b/src/mom5/ocean_core/ocean_sbc.F90 index d0ac92f0a5..4d40cecd3e 100644 --- a/src/mom5/ocean_core/ocean_sbc.F90 +++ b/src/mom5/ocean_core/ocean_sbc.F90 @@ -522,7 +522,7 @@ module ocean_sbc_mod use ocean_parameters_mod, only: MOM_BGRID, MOM_CGRID use ocean_riverspread_mod, only: spread_river_horz use ocean_tempsalt_mod, only: pottemp_from_contemp -use ocean_tpm_mod, only: ocean_tpm_sum_sfc, ocean_tpm_avg_sfc, ocean_tpm_sbc +use ocean_tpm_mod, only: ocean_tpm_sum_sfc, ocean_tpm_avg_sfc, ocean_tpm_sbc, ocean_tpm_sbc_adjust use ocean_tpm_mod, only: ocean_tpm_zero_sfc, ocean_tpm_sfc_end use ocean_types_mod, only: ocean_grid_type, ocean_domain_type, ocean_public_type use ocean_types_mod, only: ocean_time_type, ocean_thickness_type @@ -4466,7 +4466,7 @@ end subroutine get_ocean_sbc ! ! -subroutine flux_adjust(Time, T_diag, Dens, Ext_mode, T_prog, Velocity, river, melt, pme) +subroutine flux_adjust(Time, T_diag, Dens, Ext_mode, T_prog, Velocity, river, melt, pme, runoff) #if defined(ACCESS_CM) || defined(ACCESS_OM) use auscom_ice_parameters_mod, only : use_ioaice, aice_cutoff @@ -4482,6 +4482,7 @@ subroutine flux_adjust(Time, T_diag, Dens, Ext_mode, T_prog, Velocity, river, me real, dimension(isd:,jsd:), intent(in) :: river real, dimension(isd:,jsd:), intent(in) :: melt real, dimension(isd:,jsd:), intent(inout) :: pme + real, dimension(isd:,jsd:), intent(inout) :: runoff real, dimension(isd:ied,jsd:jed) :: open_ocean_mask real, dimension(isd:ied,jsd:jed) :: pme_restore, flx_restore @@ -4766,6 +4767,9 @@ subroutine flux_adjust(Time, T_diag, Dens, Ext_mode, T_prog, Velocity, river, me endif ! endif for if (id_correction(index_salt) > 0 ) + ! apply adjustments to tracer surface boundary conditions due to salt restoring/correction + + call ocean_tpm_sbc_adjust(Dom, T_prog, flx_restore(:,:)+flx_correct(:,:), runoff) ! diagnostics for salinity or pme restoring and correction diff --git a/src/mom5/ocean_param/sources/ocean_shortwave.F90 b/src/mom5/ocean_param/sources/ocean_shortwave.F90 index 6b5140486f..5cf141da43 100644 --- a/src/mom5/ocean_param/sources/ocean_shortwave.F90 +++ b/src/mom5/ocean_param/sources/ocean_shortwave.F90 @@ -294,11 +294,11 @@ subroutine sw_source (Time, Thickness, Dens, T_prog, T_diag, swflx, swflx_vis, T if(use_shortwave_gfdl) then call sw_source_gfdl (Time, Thickness, T_diag(:), swflx, swflx_vis, index_irr, Temp, sw_frac_zt, opacity) elseif(use_shortwave_csiro) then - call sw_source_csiro (Time, Thickness, T_diag(:), swflx, index_irr, Temp, sw_frac_zt) + call sw_source_csiro (Time, Thickness, T_diag(:), swflx, index_irr, Temp, sw_frac_zt, opacity) elseif(use_shortwave_jerlov) then call sw_source_jerlov (Thickness, T_diag(:), swflx, swflx_vis, index_irr, Temp, sw_frac_zt, opacity) elseif(use_shortwave_ext) then - call sw_source_ext(Time, Thickness, T_diag(:), swflx, Temp, sw_frac_zt) + call sw_source_ext(Time, Thickness, T_diag(:), swflx, Temp, sw_frac_zt, opacity) endif ! add heating rate to thickness*density weighted temperature @@ -364,7 +364,7 @@ end subroutine sw_source ! ! -subroutine sw_source_ext (Time, Thickness, T_diag, swflx, Temp, sw_frac_zt) +subroutine sw_source_ext (Time, Thickness, T_diag, swflx, Temp, sw_frac_zt, opacity) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness @@ -372,6 +372,7 @@ subroutine sw_source_ext (Time, Thickness, T_diag, swflx, Temp, sw_frac_zt) real, dimension(isd:,jsd:) , intent(in) :: swflx type(ocean_prog_tracer_type), intent(inout) :: Temp real, dimension(isd:,jsd:,:), intent(inout) :: sw_frac_zt + real, dimension(isd:,jsd:,:), intent(inout) :: opacity integer :: k, j, i real :: sw_heat_rate(isd:ied,jsd:jed,1:nk) @@ -381,7 +382,7 @@ subroutine sw_source_ext (Time, Thickness, T_diag, swflx, Temp, sw_frac_zt) ! ! for our testing purposes we'll just invoke one of our methods ! Note here that index_irr is "global" - call sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_frac_zt) + call sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_frac_zt, opacity) ! the heating rate is stored in sw_heat_rate = Temp%wrk1 diff --git a/src/mom5/ocean_param/sources/ocean_shortwave_csiro.F90 b/src/mom5/ocean_param/sources/ocean_shortwave_csiro.F90 index 43f4072734..a960d2cef4 100644 --- a/src/mom5/ocean_param/sources/ocean_shortwave_csiro.F90 +++ b/src/mom5/ocean_param/sources/ocean_shortwave_csiro.F90 @@ -100,6 +100,7 @@ module ocean_shortwave_csiro_mod use mpp_mod, only: mpp_clock_id, mpp_clock_begin, mpp_clock_end use mpp_mod, only: CLOCK_ROUTINE use time_interp_external_mod, only: time_interp_external, init_external_field +use constants_mod, only: epsln use ocean_domains_mod, only: get_local_indices use ocean_types_mod, only: ocean_time_type, ocean_domain_type, ocean_grid_type @@ -306,7 +307,7 @@ end subroutine ocean_shortwave_csiro_init ! surface tracer flux "stf(i,j,temp)" ! ! -subroutine sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_frac_zt) +subroutine sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_frac_zt, opacity) type(ocean_time_type), intent(in) :: Time type(ocean_thickness_type), intent(in) :: Thickness @@ -315,6 +316,7 @@ subroutine sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_ integer, intent(in) :: index_irr type(ocean_prog_tracer_type), intent(inout) :: Temp real, dimension(isd:,jsd:,:), intent(inout) :: sw_frac_zt + real, dimension(isd:,jsd:,:), intent(inout) :: opacity real, dimension(isd:ied,jsd:jed,0:nk) :: sw_frac_zw real, dimension(isd:ied,jsd:jed) :: zt_sw @@ -422,6 +424,29 @@ subroutine sw_source_csiro (Time, Thickness, T_diag, swflx, index_irr, Temp, sw_ enddo enddo + ! Compute opacity, which is used by biogeochemistry. + ! We split off the k=1 level, since sw_frac_zw(k=0)=sw_frac_top, + ! which is typically set to sw_frac_top=0.0 for purposes of + ! accounting (as swflx is also in stf(index_temp). A value + ! of sw_frac_zw(k=0)=0.0 to compute opacity would result in a + ! negative opacity at k=1, which is not physical. Instead, for + ! purposes of opacity calculation, we need sw_frac_zw(k=0)=1.0. + k=1 + do j=jsc,jec + do i=isc,iec + opacity(i,j,k) = -log( sw_frac_zw(i,j,k) / (F_vis + epsln) + epsln) & + / (Thickness%dzt(i,j,k) + epsln) + enddo + enddo + do k=2,nk-1 + do j=jsc,jec + do i=isc,iec + opacity(i,j,k) = -log( sw_frac_zw(i,j,k) / (sw_frac_zw(i,j,k-1) + epsln) + epsln) & + / (Thickness%dzt(i,j,k) + epsln) + enddo + enddo + enddo + end subroutine sw_source_csiro ! NAME="sw_source_csiro" diff --git a/src/mom5/ocean_tracers/ocean_tpm.F90 b/src/mom5/ocean_tracers/ocean_tpm.F90 index 576dde39fb..b8cd30a1ad 100644 --- a/src/mom5/ocean_tracers/ocean_tpm.F90 +++ b/src/mom5/ocean_tracers/ocean_tpm.F90 @@ -46,6 +46,10 @@ module ocean_tpm_mod !{ ! this functionality may be moved into a new, generalized ! boundary condition manager. ! +! ocean_tpm_sbc_adjust: Calls specified routines to adjust +! surface boundary condition, e.g. due to virtual fluxes +! that arise from salinity restoring +! ! ocean_tpm_bbc: Calls specified routines to handle bottom ! coundary condition calculations. ! @@ -221,6 +225,7 @@ module ocean_tpm_mod !{ use ocean_generic_mod, only: ocean_generic_sum_sfc use ocean_generic_mod, only: ocean_generic_zero_sfc use ocean_generic_mod, only: ocean_generic_sbc +use ocean_generic_mod, only: ocean_generic_sbc_adjust use ocean_generic_mod, only: ocean_generic_init use ocean_generic_mod, only: ocean_generic_column_physics use ocean_generic_mod, only: ocean_generic_end @@ -267,6 +272,7 @@ module ocean_tpm_mod !{ public ocean_tpm_init public ocean_tpm_flux_init public ocean_tpm_sbc +public ocean_tpm_sbc_adjust public ocean_tpm_source public ocean_tpm_start public ocean_tpm_tracer @@ -1336,6 +1342,54 @@ subroutine ocean_tpm_sbc(Domain, Grid, T_prog, Time, Ice_ocean_boundary_fluxes, end subroutine ocean_tpm_sbc !} ! NAME="ocean_tpm_sbc" +!####################################################################### +! +! +! +! call subroutines to adjust surface boundary condition, e.g. due to +! virtual fluxes that arise from salinity restoring +! +! + +subroutine ocean_tpm_sbc_adjust(Domain, T_prog, salt_flux_added, runoff) + + +implicit none + +! +!----------------------------------------------------------------------- +! Arguments +!----------------------------------------------------------------------- +! + +type(ocean_domain_type), intent(in) :: Domain +type(ocean_prog_tracer_type), dimension(:), intent(inout) :: T_prog +real, dimension(Domain%isd:,Domain%jsd:), intent(in) :: salt_flux_added +real, dimension(Domain%isd:,Domain%jsd:), intent(in) :: runoff + +! +!----------------------------------------------------------------------- +! local parameters +!----------------------------------------------------------------------- +! + +! +!----------------------------------------------------------------------- +! local variables +!----------------------------------------------------------------------- +! + + +#ifdef USE_OCEAN_BGC + +if (do_generic_tracer) call ocean_generic_sbc_adjust(Domain%isd, Domain%jsd, T_prog, salt_flux_added, runoff) + +#endif + +return + +end subroutine ocean_tpm_sbc_adjust !} +! NAME="ocean_tpm_sbc_adjust" !####################################################################### ! diff --git a/src/ocean_shared/generic_tracers/generic_tracer.F90 b/src/ocean_shared/generic_tracers/generic_tracer.F90 index 68c3ac74cf..c358f3f5e5 100644 --- a/src/ocean_shared/generic_tracers/generic_tracer.F90 +++ b/src/ocean_shared/generic_tracers/generic_tracer.F90 @@ -86,7 +86,9 @@ module generic_tracer public generic_tracer_source public generic_tracer_diag public generic_tracer_update_from_bottom + public generic_tracer_update_from_coupler public generic_tracer_coupler_get + public generic_tracer_coupler_accumulate public generic_tracer_coupler_set public generic_tracer_coupler_zero public generic_tracer_end @@ -304,6 +306,28 @@ subroutine generic_tracer_coupler_get(IOB_struc) end subroutine generic_tracer_coupler_get + ! + ! + ! Accumulate the boundary values (%stf and %triver) from coupler fluxes + ! + ! + ! This subroutine accumulates coupler values for those generic tracers that have flux + ! exchange with atmosphere. This routine simply wraps g_tracer_coupler_get for MOM5, + ! but the API is provided here for consistency with https://github.com/NOAA-GFDL/ocean_BGC + ! + ! + ! + subroutine generic_tracer_coupler_accumulate(IOB_struc, weight) + type(coupler_2d_bc_type), intent(in) :: IOB_struc + real, intent(in) :: weight + + !All generic tracers + !Update tracer boundary values (%stf and %triver) from coupler fluxes foreach tracer in the prog_tracer_list + call g_tracer_coupler_get(tracer_list, IOB_struc) + + end subroutine generic_tracer_coupler_accumulate ! ! @@ -481,6 +505,44 @@ subroutine generic_tracer_update_from_bottom(dt, tau, model_time) end subroutine generic_tracer_update_from_bottom + ! + ! + ! Modify the values obtained from the coupler + ! + ! + ! Calls the corresponding generic_X_update_from_coupler routine for each package X. + ! + ! + ! + ! Lower bounds of x and y extents of input arrays on data domain + ! + ! + ! Surface salt flux into ocean from restoring or flux adjustment [g/m^2/sec] + ! + ! + + subroutine generic_tracer_update_from_coupler(ilb, jlb, salt_flux_added) + integer, intent(in) :: ilb, jlb + real, dimension(ilb:,jlb:), intent(in) :: salt_flux_added + + character(len=fm_string_len), parameter :: sub_name = 'generic_tracer_update_from_coupler' + + !Specific tracers + ! if(do_generic_CFC) call generic_CFC_update_from_coupler(tracer_list) !Nothing to do + + if(do_generic_TOPAZ) call generic_TOPAZ_update_from_coupler(tracer_list) + + if(do_generic_BLING) call generic_BLING_update_from_coupler(tracer_list) + + if(do_generic_miniBLING) call generic_miniBLING_update_from_coupler(tracer_list) + + if(do_generic_COBALT) call generic_COBALT_update_from_coupler(tracer_list) + + return + + end subroutine generic_tracer_update_from_coupler ! !