From 276954f3cf5005858a713e324c348b7fd677aac4 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 14 Jun 2021 08:40:19 -0600 Subject: [PATCH 1/2] Fixed a misplaced parentheses This commit fixes misplaced parentheses in the tridiagonal solvers used in subroutines tracer_vertdiff_Eulerian and tracer_vertdiff. This bug has been reported in https://github.com/NOAA-GFDL/MOM6/issues/1415 This commit also changes the mask condition used in these subroutines. --- src/tracer/MOM_tracer_diabatic.F90 | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index 9be4af08dc..cbfe2d7f0a 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -152,7 +152,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) d1(i) = b_denom_1 * b1(i) h_tr = h_old(i,j,1) + h_neglect - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) endif ; enddo do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = eb(i,j,k-1) * b1(i) @@ -185,14 +185,14 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & else !$OMP do do j=js,je - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then h_tr = h_old(i,j,1) + h_neglect b_denom_1 = h_tr + ea(i,j,1) b1(i) = 1.0 / (b_denom_1 + eb(i,j,1)) d1(i) = h_tr * b1(i) - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) endif ; enddo - do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = eb(i,j,k-1) * b1(i) h_tr = h_old(i,j,k) + h_neglect b_denom_1 = h_tr + d1(i) * ea(i,j,k) @@ -200,7 +200,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & d1(i) = b_denom_1 * b1(i) tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ea(i,j,k) * tr(i,j,k-1)) endif ; enddo ; enddo - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,nz) = eb(i,j,nz-1) * b1(i) h_tr = h_old(i,j,nz) + h_neglect b_denom_1 = h_tr + d1(i)*ea(i,j,nz) @@ -208,7 +208,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + & ea(i,j,nz) * tr(i,j,nz-1)) endif ; enddo - do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) endif ; enddo ; enddo enddo @@ -267,8 +267,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & !! ensure positive definiteness [H ~> m or kg m-2]. real :: h_neglect !< A thickness that is so small it is usually lost !! in roundoff and can be neglected [H ~> m or kg m-2]. - logical :: convert_flux = .true. - + logical :: convert_flux integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -279,6 +278,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & return endif + convert_flux = .true. if (present(convert_flux_in)) convert_flux = convert_flux_in h_neglect = GV%H_subroundoff sink_dist = 0.0 @@ -351,7 +351,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) d1(i) = b_denom_1 * b1(i) h_tr = h_old(i,j,1) + h_neglect - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) endif ; enddo do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = ent(i,j,K) * b1(i) @@ -384,14 +384,14 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & else !$OMP do do j=js,je - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then h_tr = h_old(i,j,1) + h_neglect b_denom_1 = h_tr + ent(i,j,1) b1(i) = 1.0 / (b_denom_1 + ent(i,j,2)) d1(i) = h_tr * b1(i) - tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j) + tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j)) endif ; enddo - do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,k) = ent(i,j,K) * b1(i) h_tr = h_old(i,j,k) + h_neglect b_denom_1 = h_tr + d1(i) * ent(i,j,K) @@ -399,7 +399,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & d1(i) = b_denom_1 * b1(i) tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ent(i,j,K) * tr(i,j,k-1)) endif ; enddo ; enddo - do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then c1(i,nz) = ent(i,j,nz) * b1(i) h_tr = h_old(i,j,nz) + h_neglect b_denom_1 = h_tr + d1(i)*ent(i,j,nz) @@ -407,7 +407,7 @@ subroutine tracer_vertdiff_Eulerian(h_old, ent, dt, tr, G, GV, & tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + & ent(i,j,nz) * tr(i,j,nz-1)) endif ; enddo - do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then + do k=nz-1,1,-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1) endif ; enddo ; enddo enddo From 3aade32dd6b8033603d3e1b21678dfc26c8d3281 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 14 Jun 2021 12:30:37 -0600 Subject: [PATCH 2/2] CFCs implementation via NUOPC cap This commit adds a new module (MOM_CFC_cap) that can used to simulate chlorofluorocarbons (CFCs) tracers. It also includes the modifications needed to 1) allocate arrays in the fluxes and surface types, and 2) activate the MOM_CFC_cap module via NUOPC cap. --- config_src/drivers/nuopc_cap/mom_cap.F90 | 12 +- .../drivers/nuopc_cap/mom_cap_methods.F90 | 17 +- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 14 +- .../nuopc_cap/mom_surface_forcing_nuopc.F90 | 66 +- src/core/MOM_forcing_type.F90 | 74 +- src/core/MOM_unit_tests.F90 | 4 +- src/core/MOM_variables.F90 | 17 +- src/tracer/MOM_CFC_cap.F90 | 704 ++++++++++++++++++ src/tracer/MOM_tracer_flow_control.F90 | 33 +- 9 files changed, 922 insertions(+), 19 deletions(-) create mode 100644 src/tracer/MOM_CFC_cap.F90 diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 2d79674606..42cc579546 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -31,7 +31,7 @@ module MOM_cap_mod use MOM_ocean_model_nuopc, only: ice_ocean_boundary_type use MOM_grid, only: ocean_grid_type, get_global_grid_size use MOM_ocean_model_nuopc, only: ocean_model_restart, ocean_public_type, ocean_state_type -use MOM_ocean_model_nuopc, only: ocean_model_init_sfc +use MOM_ocean_model_nuopc, only: ocean_model_init_sfc, ocean_model_flux_init use MOM_ocean_model_nuopc, only: ocean_model_init, update_ocean_model, ocean_model_end use MOM_ocean_model_nuopc, only: get_ocean_grid, get_eps_omesh use MOM_cap_time, only: AlarmInit @@ -649,6 +649,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ocean_public%is_ocean_pe = .true. call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(restartfiles)) + ! GMM, this call is not needed for NCAR. Check with EMC. + ! If this can be deleted, perhaps we should also delete ocean_model_flux_init + call ocean_model_flux_init(ocean_state) + call ocean_model_init_sfc(ocean_state, ocean_public) call mpp_get_compute_domain(ocean_public%domain, isc, iec, jsc, jec) @@ -668,6 +672,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary% seaice_melt_heat (isc:iec,jsc:jec),& Ice_ocean_boundary% seaice_melt (isc:iec,jsc:jec), & Ice_ocean_boundary% mi (isc:iec,jsc:jec), & + Ice_ocean_boundary% ice_fraction (isc:iec,jsc:jec), & + Ice_ocean_boundary% u10_sqr (isc:iec,jsc:jec), & Ice_ocean_boundary% p (isc:iec,jsc:jec), & Ice_ocean_boundary% lrunoff_hflx (isc:iec,jsc:jec), & Ice_ocean_boundary% frunoff_hflx (isc:iec,jsc:jec), & @@ -689,6 +695,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) Ice_ocean_boundary%seaice_melt = 0.0 Ice_ocean_boundary%seaice_melt_heat= 0.0 Ice_ocean_boundary%mi = 0.0 + Ice_ocean_boundary%ice_fraction = 0.0 + Ice_ocean_boundary%u10_sqr = 0.0 Ice_ocean_boundary%p = 0.0 Ice_ocean_boundary%lrunoff_hflx = 0.0 Ice_ocean_boundary%frunoff_hflx = 0.0 @@ -747,6 +755,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call fld_list_add(fldsToOcn_num, fldsToOcn, "inst_pres_height_surface" , "will provide") call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofl" , "will provide") !-> liquid runoff call fld_list_add(fldsToOcn_num, fldsToOcn, "Foxx_rofi" , "will provide") !-> ice runoff + call fld_list_add(fldsToOcn_num, fldsToOcn, "Si_ifrac" , "will provide") !-> ice fraction + call fld_list_add(fldsToOcn_num, fldsToOcn, "So_duu10n" , "will provide") !-> wind^2 at 10m call fld_list_add(fldsToOcn_num, fldsToOcn, "mean_fresh_water_to_ocean_rate", "will provide") call fld_list_add(fldsToOcn_num, fldsToOcn, "net_heat_flx_to_ocn" , "will provide") !These are not currently used and changing requires a nuopc dictionary change diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 0a4e12c647..b7e3b13c1e 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -250,12 +250,27 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, !---- ! Note - preset values to 0, if field does not exist in importState, then will simply return ! and preset value will be used - ice_ocean_boundary%mi(:,:) = 0._ESMF_KIND_R8 call state_getimport(importState, 'mass_of_overlying_ice', & isc, iec, jsc, jec, ice_ocean_boundary%mi,rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + !---- + ! sea-ice fraction + !---- + ice_ocean_boundary%ice_fraction(:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'Si_ifrac', & + isc, iec, jsc, jec, ice_ocean_boundary%ice_fraction, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + !---- + ! 10m wind squared + !---- + ice_ocean_boundary%u10_sqr(:,:) = 0._ESMF_KIND_R8 + call state_getimport(importState, 'So_duu10n', & + isc, iec, jsc, jec, ice_ocean_boundary%u10_sqr, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + !---- ! Partitioned Stokes Drift Components !---- diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index 493762f4bc..9d5f81f079 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -40,7 +40,6 @@ module MOM_ocean_model_nuopc use MOM_time_manager, only : operator(/=), operator(<=), operator(>=) use MOM_time_manager, only : operator(<), real_to_time_type, time_type_to_real use time_interp_external_mod,only : time_interp_external_init -use MOM_tracer_flow_control, only : call_tracer_register, tracer_flow_control_init use MOM_tracer_flow_control, only : call_tracer_flux_init use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : surface @@ -240,6 +239,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! tracer fluxes, and can be used to spawn related !! internal variables in the ice model. character(len=*), optional, intent(in) :: input_restart_file !< If present, name of restart file to read + ! Local variables real :: Rho0 ! The Boussinesq ocean density, in kg m-3. real :: G_Earth ! The gravitational acceleration in m s-2. @@ -247,7 +247,9 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! The actual depth over which melt potential is computed will !! min(HFrz, OBLD), where OBLD is the boundary layer depth. !! If HFrz <= 0 (default), melt potential will not be computed. - logical :: use_melt_pot!< If true, allocate melt_potential array + logical :: use_melt_pot !< If true, allocate melt_potential array + logical :: use_CFC !< If true, allocated arrays for surface CFCs. + ! This include declares and sets the variable "version". #include "version_variable.h" @@ -366,10 +368,14 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i use_melt_pot=.false. endif + call get_param(param_file, mdl, "USE_CFC_CAP", use_CFC, & + default=.false., do_not_log=.true.) + ! Consider using a run-time flag to determine whether to do the diagnostic ! vertical integrals, since the related 3-d sums are not negligible in cost. call allocate_surface_state(OS%sfc_state, OS%grid, use_temperature, & - do_integrals=.true., gas_fields_ocn=gas_fields_ocn, use_meltpot=use_melt_pot) + do_integrals=.true., gas_fields_ocn=gas_fields_ocn, & + use_meltpot=use_melt_pot, use_cfcs=use_CFC) call surface_forcing_init(Time_in, OS%grid, OS%US, param_file, OS%diag, & OS%forcing_CSp, OS%restore_salinity, OS%restore_temp) @@ -416,6 +422,8 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i endif + call extract_surface_state(OS%MOM_CSp, OS%sfc_state) + call close_param_file(param_file) call diag_mediator_close_registration(OS%diag) diff --git a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 index 7168823fbc..1054234e6b 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -21,6 +21,7 @@ module MOM_surface_forcing_nuopc use MOM_forcing_type, only : allocate_mech_forcing, deallocate_mech_forcing use MOM_get_input, only : Get_MOM_Input, directories use MOM_grid, only : ocean_grid_type +use MOM_CFC_cap, only : CFC_cap_fluxes use MOM_io, only : slasher, write_version_number, MOM_read_data use MOM_restart, only : register_restart_field, restart_init, MOM_restart_CS use MOM_restart, only : restart_init_end, save_restart, restore_state @@ -79,7 +80,7 @@ module MOM_surface_forcing_nuopc !! the correction for the atmospheric (and sea-ice) !! pressure limited by max_p_surf instead of the !! full atmospheric pressure. The default is true. - + logical :: use_CFC !< enables the MOM_CFC_cap tracer package. real :: gust_const !< constant unresolved background gustiness for ustar [R L Z T-1 ~> Pa] logical :: read_gust_2d !< If true, use a 2-dimensional gustiness supplied !! from an input file. @@ -128,6 +129,7 @@ module MOM_surface_forcing_nuopc type(diag_ctrl), pointer :: diag !< structure to regulate diagnostic output timing character(len=200) :: inputdir !< directory where NetCDF input files are + character(len=200) :: CFC_BC_file !< filename with cfc11 and cfc12 data character(len=200) :: salt_restore_file !< filename for salt restoring data character(len=30) :: salt_restore_var_name !< name of surface salinity in salt_restore_file logical :: mask_srestore !< if true, apply a 2-dimensional mask to the surface @@ -141,9 +143,13 @@ module MOM_surface_forcing_nuopc !! temperature restoring fluxes. The masking file should be !! in inputdir/temp_restore_mask.nc and the field should !! be named 'mask' + character(len=30) :: cfc11_var_name !< name of cfc11 in CFC_BC_file + character(len=30) :: cfc12_var_name !< name of cfc11 in CFC_BC_file real, pointer, dimension(:,:) :: trestore_mask => NULL() !< mask for SST restoring - integer :: id_srestore = -1 !< id number for time_interp_external. - integer :: id_trestore = -1 !< id number for time_interp_external. + integer :: id_srestore = -1 !< id number for time_interp_external. + integer :: id_trestore = -1 !< id number for time_interp_external. + integer :: id_cfc11_atm = -1 !< id number for time_interp_external. + integer :: id_cfc12_atm = -1 !< id number for time_interp_external. ! Diagnostics handles type(forcing_diags), public :: handles @@ -178,6 +184,8 @@ module MOM_surface_forcing_nuopc real, pointer, dimension(:,:) :: frunoff_hflx =>NULL() !< heat content of frozen runoff [W/m2] real, pointer, dimension(:,:) :: p =>NULL() !< pressure of overlying ice and atmosphere !< on ocean surface [Pa] + real, pointer, dimension(:,:) :: ice_fraction =>NULL() !< mass of ice [nondim] + real, pointer, dimension(:,:) :: u10_sqr =>NULL() !< wind speed squared at 10m [m2/s2] real, pointer, dimension(:,:) :: mi =>NULL() !< mass of ice [kg/m2] real, pointer, dimension(:,:) :: ice_rigidity =>NULL() !< rigidity of the sea ice, sea-ice and !! ice-shelves, expressed as a coefficient @@ -234,6 +242,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! local variables real, dimension(SZI_(G),SZJ_(G)) :: & + cfc11_atm, & !< CFC11 concentration in the atmopshere [???????] + cfc12_atm, & !< CFC11 concentration in the atmopshere [???????] data_restore, & !< The surface value toward which to restore [g/kg or degC] SST_anom, & !< Instantaneous sea surface temperature anomalies from a target value [deg C] SSS_anom, & !< Instantaneous sea surface salinity anomalies from a target value [g/kg] @@ -293,7 +303,8 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! flux type has been used. if (fluxes%dt_buoy_accum < 0) then call allocate_forcing_type(G, fluxes, water=.true., heat=.true., ustar=.true., & - press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug) + press=.true., fix_accum_bug=CS%fix_ustar_gustless_bug, & + cfc=CS%use_CFC) call safe_alloc_ptr(fluxes%sw_vis_dir,isd,ied,jsd,jed) call safe_alloc_ptr(fluxes%sw_vis_dif,isd,ied,jsd,jed) @@ -433,6 +444,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, enddo ; enddo endif + ! Check that liquid runoff has a place to go if (CS%liquid_runoff_from_data .and. .not. associated(IOB%lrunoff)) then call MOM_error(FATAL, "liquid runoff is being added via data_override but "// & @@ -496,6 +508,14 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, if (associated(IOB%seaice_melt)) & fluxes%seaice_melt(i,j) = kg_m2_s_conversion * G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0) + ! sea ice fraction [nondim] + if (associated(IOB%ice_fraction)) & + fluxes%ice_fraction(i,j) = G%mask2dT(i,j) * IOB%ice_fraction(i-i0,j-j0) + + ! 10-m wind speed squared [m2/s2] + if (associated(IOB%u10_sqr)) & + fluxes%u10_sqr(i,j) = US%m_to_L**2 * US%T_to_s**2 * G%mask2dT(i,j) * IOB%u10_sqr(i-i0,j-j0) + fluxes%latent(i,j) = 0.0 ! notice minus sign since fprec is positive into the ocean if (associated(IOB%fprec)) then @@ -550,6 +570,11 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, fluxes%accumulate_p_surf = .true. ! Multiple components may contribute to surface pressure. endif + ! CFCs + if (CS%use_CFC) then + call CFC_cap_fluxes(fluxes, sfc_state, G, CS%Rho0, Time, CS%id_cfc11_atm, CS%id_cfc11_atm) + endif + if (associated(IOB%salt_flux)) then do j=js,je ; do i=is,ie fluxes%salt_flux(i,j) = G%mask2dT(i,j)*(fluxes%salt_flux(i,j) + kg_m2_s_conversion*IOB%salt_flux(i-i0,j-j0)) @@ -1161,6 +1186,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "coupler. This is used for testing and should be =1.0 for any "//& "production runs.", default=1.0) + call get_param(param_file, mdl, "USE_CFC_CAP", CS%use_CFC, & + default=.false., do_not_log=.true.) + if (restore_salt) then call get_param(param_file, mdl, "FLUXCONST", CS%Flux_const, & "The constant that relates the restoring surface fluxes to the relative "//& @@ -1173,7 +1201,6 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "The name of the surface salinity variable to read from "//& "SALT_RESTORE_FILE for restoring salinity.", & default="salt") - call get_param(param_file, mdl, "SRESTORE_AS_SFLUX", CS%salt_restore_as_sflux, & "If true, the restoring of salinity is applied as a salt "//& "flux instead of as a freshwater flux.", default=.false.) @@ -1269,7 +1296,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, enddo ; enddo endif - call time_interp_external_init + ! initialize time interpolator module + call time_interp_external_init() ! Optionally read a x-y gustiness field in place of a global ! constant. @@ -1320,8 +1348,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, call get_param(param_file, mdl, "ALLOW_ICEBERG_FLUX_DIAGNOSTICS", iceberg_flux_diags, & "If true, makes available diagnostics of fluxes from icebergs "//& "as seen by MOM6.", default=.false.) + call register_forcing_type_diags(Time, diag, US, CS%use_temperature, CS%handles, & - use_berg_fluxes=iceberg_flux_diags) + use_berg_fluxes=iceberg_flux_diags, use_cfcs=CS%use_CFC) call get_param(param_file, mdl, "ALLOW_FLUX_ADJUSTMENTS", CS%allow_flux_adjustments, & "If true, allows flux adjustments to specified via the "//& @@ -1355,6 +1384,27 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, endif endif ; endif + ! Do not log these params here since they are logged in the CFC cap module + if (CS%use_CFC) then + call get_param(param_file, mdl, "CFC_BC_FILE", CS%CFC_BC_file, & + "The file in which the CFC-11 and CFC-12 atm concentrations can be "//& + "found (units must be parts per trillion), or an empty string for "//& + "internal BC generation (TODO).", default=" ", do_not_log=.true.) + if ((len_trim(CS%CFC_BC_file) > 0) .and. (scan(CS%CFC_BC_file,'/') == 0)) then + ! Add the directory if CFC_BC_file is not already a complete path. + CS%CFC_BC_file = trim(slasher(CS%inputdir))//trim(CS%CFC_BC_file) + call get_param(param_file, mdl, "CFC11_VARIABLE", CS%cfc11_var_name, & + "The name of the variable representing CFC-11 in "//& + "CFC_BC_FILE.", default="CFC_11", do_not_log=.true.) + call get_param(param_file, mdl, "CFC12_VARIABLE", CS%cfc12_var_name, & + "The name of the variable representing CFC-12 in "//& + "CFC_BC_FILE.", default="CFC_12", do_not_log=.true.) + + CS%id_cfc11_atm = init_external_field(CS%CFC_BC_file, CS%cfc11_var_name, domain=G%Domain%mpp_domain) + CS%id_cfc12_atm = init_external_field(CS%CFC_BC_file, CS%cfc12_var_name, domain=G%Domain%mpp_domain) + endif + endif + ! Set up any restart fields associated with the forcing. call restart_init(param_file, CS%restart_CSp, "MOM_forcing.res") call restart_init_end(CS%restart_CSp) @@ -1424,6 +1474,8 @@ subroutine ice_ocn_bnd_type_chksum(id, timestep, iobt) write(outunit,100) 'iobt%lrunoff ' , mpp_chksum( iobt%lrunoff ) write(outunit,100) 'iobt%frunoff ' , mpp_chksum( iobt%frunoff ) write(outunit,100) 'iobt%p ' , mpp_chksum( iobt%p ) + write(outunit,100) 'iobt%ice_fraction ' , mpp_chksum( iobt%ice_fraction ) + write(outunit,100) 'iobt%u10_sqr ' , mpp_chksum( iobt%u10_sqr ) if (associated(iobt%ustar_berg)) & write(outunit,100) 'iobt%ustar_berg ' , mpp_chksum( iobt%ustar_berg ) if (associated(iobt%area_berg)) & diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 682ad03397..7d8375a8af 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -183,6 +183,13 @@ module MOM_forcing_type real :: C_p !< heat capacity of seawater [Q degC-1 ~> J kg-1 degC-1]. !! C_p is is the same value as in thermovar_ptrs_type. + ! CFC-related arrays needed in the MOM_CFC_cap module + real, pointer, dimension(:,:) :: & + cfc11_flux => NULL(), & !< flux of cfc_11 into the ocean [mol m-2 s-1]. + cfc12_flux => NULL(), & !< flux of cfc_12 into the ocean [mol m-2 s-1]. + ice_fraction => NULL(), & !< fraction of sea ice coverage at h-cells, from 0 to 1 [nondim]. + u10_sqr => NULL() !< wind magnitude at 10 m squared [L2 T-2 ~> m2 s-2] + ! passive tracer surface fluxes type(coupler_2d_bc_type) :: tr_fluxes !< This structure contains arrays of !! of named fields used for passive tracer fluxes. @@ -348,6 +355,12 @@ module MOM_forcing_type integer :: id_TKE_tidal = -1 integer :: id_buoy = -1 + ! cfc-related diagnostics handles + integer :: id_cfc11 = -1 + integer :: id_cfc12 = -1 + integer :: id_ice_fraction = -1 + integer :: id_u10_sqr = -1 + ! iceberg diagnostic handles integer :: id_ustar_berg = -1 integer :: id_area_berg = -1 @@ -1087,6 +1100,10 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) haloshift=hshift, scale=US%QRZ_T_to_W_m2) if (associated(fluxes%p_surf)) & call hchksum(fluxes%p_surf, mesg//" fluxes%p_surf", G%HI, haloshift=hshift , scale=US%RL2_T2_to_Pa) + if (associated(fluxes%u10_sqr)) & + call hchksum(fluxes%u10_sqr, mesg//" fluxes%u10_sqr", G%HI, haloshift=hshift , scale=US%Z_to_m**2*US%s_to_T**2) + if (associated(fluxes%ice_fraction)) & + call hchksum(fluxes%ice_fraction, mesg//" fluxes%ice_fraction", G%HI, haloshift=hshift) if (associated(fluxes%salt_flux)) & call hchksum(fluxes%salt_flux, mesg//" fluxes%salt_flux", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%TKE_tidal)) & @@ -1239,13 +1256,14 @@ end subroutine forcing_SinglePointPrint !> Register members of the forcing type for diagnostics -subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes) +subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, use_berg_fluxes, use_cfcs) type(time_type), intent(in) :: Time !< time type type(diag_ctrl), intent(inout) :: diag !< diagnostic control type type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type logical, intent(in) :: use_temperature !< True if T/S are in use type(forcing_diags), intent(inout) :: handles !< handles for diagnostics logical, optional, intent(in) :: use_berg_fluxes !< If true, allow iceberg flux diagnostics + logical, optional, intent(in) :: use_cfcs !< If true, allow cfc related diagnostics ! Clock for forcing diagnostics handles%id_clock_forcing=cpu_clock_id('(Ocean forcing diagnostics)', grain=CLOCK_ROUTINE) @@ -1288,6 +1306,34 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, endif endif + ! units for cfc11_flux and cfc12_flux are mol m-2 s-1 + ! See: + ! http://clipc-services.ceda.ac.uk/dreq/u/0940cbee6105037e4b7aa5579004f124.html + ! http://clipc-services.ceda.ac.uk/dreq/u/e9e21426e4810d0bb2d3dddb24dbf4dc.html + if (present(use_cfcs)) then + if (use_cfcs) then + handles%id_cfc11 = register_diag_field('ocean_model', 'cfc11_flux', diag%axesT1, Time, & + 'Gas exchange flux of CFC11 into the ocean ', 'mol m-2 s-1', & + conversion= US%m_to_Z**2*US%s_to_T,& + cmor_field_name='fgcfc11', & + cmor_long_name='Surface Downward CFC11 Flux', & + cmor_standard_name='surface_downward_cfc11_flux') + + handles%id_cfc12 = register_diag_field('ocean_model', 'cfc12_flux', diag%axesT1, Time, & + 'Gas exchange flux of CFC12 into the ocean ', 'mol m-2 s-1', & + conversion= US%m_to_Z**2*US%s_to_T,& + cmor_field_name='fgcfc12', & + cmor_long_name='Surface Downward CFC12 Flux', & + cmor_standard_name='surface_downward_cfc12_flux') + + handles%id_ice_fraction = register_diag_field('ocean_model', 'ice_fraction', diag%axesT1, Time, & + 'Fraction of cell area covered by sea ice', 'm2 m-2') + + handles%id_u10_sqr = register_diag_field('ocean_model', 'u10_sqr', diag%axesT1, Time, & + 'Wind magnitude at 10m, squared', 'm2 s-2', conversion=US%Z_to_m**2*US%s_to_T**2) + endif + endif + handles%id_psurf = register_diag_field('ocean_model', 'p_surf', diag%axesT1, Time, & 'Pressure at ice-ocean or atmosphere-ocean interface', & 'Pa', conversion=US%RL2_T2_to_Pa, cmor_field_name='pso', & @@ -2834,6 +2880,19 @@ subroutine forcing_diagnostics(fluxes_in, sfc_state, G_in, US, time_end, diag, h if (handles%id_netFWGlobalScl > 0) & call post_data(handles%id_netFWGlobalScl, fluxes%netFWGlobalScl, diag) + ! post diagnostics related to cfcs ==================================== + + if ((handles%id_cfc11 > 0) .and. associated(fluxes%cfc11_flux)) & + call post_data(handles%id_cfc11, fluxes%cfc11_flux, diag) + + if ((handles%id_cfc11 > 0) .and. associated(fluxes%cfc12_flux)) & + call post_data(handles%id_cfc12, fluxes%cfc12_flux, diag) + + if ((handles%id_ice_fraction > 0) .and. associated(fluxes%ice_fraction)) & + call post_data(handles%id_ice_fraction, fluxes%ice_fraction, diag) + + if ((handles%id_u10_sqr > 0) .and. associated(fluxes%u10_sqr)) & + call post_data(handles%id_u10_sqr, fluxes%u10_sqr, diag) ! remaining boundary terms ================================================== @@ -2872,7 +2931,7 @@ end subroutine forcing_diagnostics !> Conditionally allocate fields within the forcing type subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & - shelf, iceberg, salt, fix_accum_bug) + shelf, iceberg, salt, fix_accum_bug, cfc) type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(forcing), intent(inout) :: fluxes !< A structure containing thermodynamic forcing fields logical, optional, intent(in) :: water !< If present and true, allocate water fluxes @@ -2884,6 +2943,7 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & logical, optional, intent(in) :: salt !< If present and true, allocate salt fluxes logical, optional, intent(in) :: fix_accum_bug !< If present and true, avoid using a bug in !! accumulation of ustar_gustless + logical, optional, intent(in) :: cfc !< If present and true, allocate cfc fluxes ! Local variables integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB @@ -2939,6 +2999,12 @@ subroutine allocate_forcing_by_group(G, fluxes, water, heat, ustar, press, & call myAlloc(fluxes%area_berg,isd,ied,jsd,jed, iceberg) call myAlloc(fluxes%mass_berg,isd,ied,jsd,jed, iceberg) + !These fields should only on allocated when USE_CFC_CAP is activated. + call myAlloc(fluxes%cfc11_flux,isd,ied,jsd,jed, cfc) + call myAlloc(fluxes%cfc12_flux,isd,ied,jsd,jed, cfc) + call myAlloc(fluxes%ice_fraction,isd,ied,jsd,jed, cfc) + call myAlloc(fluxes%u10_sqr,isd,ied,jsd,jed, cfc) + if (present(fix_accum_bug)) fluxes%gustless_accum_bug = .not.fix_accum_bug end subroutine allocate_forcing_by_group @@ -3185,6 +3251,10 @@ subroutine deallocate_forcing_type(fluxes) if (associated(fluxes%ustar_berg)) deallocate(fluxes%ustar_berg) if (associated(fluxes%area_berg)) deallocate(fluxes%area_berg) if (associated(fluxes%mass_berg)) deallocate(fluxes%mass_berg) + if (associated(fluxes%ice_fraction)) deallocate(fluxes%ice_fraction) + if (associated(fluxes%u10_sqr)) deallocate(fluxes%u10_sqr) + if (associated(fluxes%cfc11_flux)) deallocate(fluxes%cfc11_flux) + if (associated(fluxes%cfc12_flux)) deallocate(fluxes%cfc12_flux) call coupler_type_destructor(fluxes%tr_fluxes) diff --git a/src/core/MOM_unit_tests.F90 b/src/core/MOM_unit_tests.F90 index 86493aad93..08f8dea634 100644 --- a/src/core/MOM_unit_tests.F90 +++ b/src/core/MOM_unit_tests.F90 @@ -11,7 +11,7 @@ module MOM_unit_tests use MOM_diag_vkernels, only : diag_vkernels_unit_tests use MOM_random, only : random_unit_tests use MOM_lateral_boundary_diffusion, only : near_boundary_unit_tests - +use MOM_CFC_cap, only : CFC_cap_unit_tests implicit none ; private public unit_tests @@ -41,6 +41,8 @@ subroutine unit_tests(verbosity) "MOM_unit_tests: random_unit_tests FAILED") if (near_boundary_unit_tests(verbose)) call MOM_error(FATAL, & "MOM_unit_tests: near_boundary_unit_tests FAILED") + if (CFC_cap_unit_tests(verbose)) call MOM_error(FATAL, & + "MOM_unit_tests: CFC_cap_unit_tests FAILED") endif end subroutine unit_tests diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 886ee77510..2ea8943e67 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -42,6 +42,8 @@ module MOM_variables SST, & !< The sea surface temperature [degC]. SSS, & !< The sea surface salinity [ppt ~> psu or gSalt/kg]. sfc_density, & !< The mixed layer density [R ~> kg m-3]. + sfc_cfc11, & !< Sea surface concentration of CFC11 [mol kg-1]. + sfc_cfc12, & !< Sea surface concentration of CFC12 [mol kg-1]. Hml, & !< The mixed layer depth [Z ~> m]. u, & !< The mixed layer zonal velocity [L T-1 ~> m s-1]. v, & !< The mixed layer meridional velocity [L T-1 ~> m s-1]. @@ -300,7 +302,8 @@ module MOM_variables !> Allocates the fields for the surface (return) properties of !! the ocean model. Unused fields are unallocated. subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & - gas_fields_ocn, use_meltpot, use_iceshelves, omit_frazil) + gas_fields_ocn, use_meltpot, use_iceshelves, & + omit_frazil, use_cfcs) type(ocean_grid_type), intent(in) :: G !< ocean grid structure type(surface), intent(inout) :: sfc_state !< ocean surface state type to be allocated. logical, optional, intent(in) :: use_temperature !< If true, allocate the space for thermodynamic variables. @@ -313,13 +316,14 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & !! tracer fluxes, and can be used to spawn related !! internal variables in the ice model. logical, optional, intent(in) :: use_meltpot !< If true, allocate the space for melt potential + logical, optional, intent(in) :: use_cfcs !< If true, allocate the space for cfcs logical, optional, intent(in) :: use_iceshelves !< If true, allocate the space for the stresses !! under ice shelves. logical, optional, intent(in) :: omit_frazil !< If present and false, do not allocate the space to !! pass frazil fluxes to the coupler ! local variables - logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil + logical :: use_temp, alloc_integ, use_melt_potential, alloc_iceshelves, alloc_frazil, alloc_cfcs integer :: is, ie, js, je, isd, ied, jsd, jed integer :: isdB, iedB, jsdB, jedB @@ -330,6 +334,7 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & use_temp = .true. ; if (present(use_temperature)) use_temp = use_temperature alloc_integ = .true. ; if (present(do_integrals)) alloc_integ = do_integrals use_melt_potential = .false. ; if (present(use_meltpot)) use_melt_potential = use_meltpot + alloc_cfcs = .false. ; if (present(use_cfcs)) alloc_cfcs = use_cfcs alloc_iceshelves = .false. ; if (present(use_iceshelves)) alloc_iceshelves = use_iceshelves alloc_frazil = .true. ; if (present(omit_frazil)) alloc_frazil = .not.omit_frazil @@ -353,6 +358,11 @@ subroutine allocate_surface_state(sfc_state, G, use_temperature, do_integrals, & allocate(sfc_state%melt_potential(isd:ied,jsd:jed)) ; sfc_state%melt_potential(:,:) = 0.0 endif + if (alloc_cfcs) then + allocate(sfc_state%sfc_cfc11(isd:ied,jsd:jed)) ; sfc_state%sfc_cfc11(:,:) = 0.0 + allocate(sfc_state%sfc_cfc12(isd:ied,jsd:jed)) ; sfc_state%sfc_cfc12(:,:) = 0.0 + endif + if (alloc_integ) then ! Allocate structures for the vertically integrated ocean_mass, ocean_heat, and ocean_salt. allocate(sfc_state%ocean_mass(isd:ied,jsd:jed)) ; sfc_state%ocean_mass(:,:) = 0.0 @@ -396,7 +406,8 @@ subroutine deallocate_surface_state(sfc_state) if (allocated(sfc_state%ocean_heat)) deallocate(sfc_state%ocean_heat) if (allocated(sfc_state%ocean_salt)) deallocate(sfc_state%ocean_salt) if (allocated(sfc_state%salt_deficit)) deallocate(sfc_state%salt_deficit) - + if (allocated(sfc_state%sfc_cfc11)) deallocate(sfc_state%sfc_cfc11) + if (allocated(sfc_state%sfc_cfc12)) deallocate(sfc_state%sfc_cfc12) call coupler_type_destructor(sfc_state%tr_fields) sfc_state%arrays_allocated = .false. diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 new file mode 100644 index 0000000000..58ba419c62 --- /dev/null +++ b/src/tracer/MOM_CFC_cap.F90 @@ -0,0 +1,704 @@ +!> Simulates CFCs using atmospheric pressure, wind speed and sea ice cover +!! provided via cap (only NUOPC cap is implemented so far). +module MOM_CFC_cap + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl, register_diag_field, post_data +use MOM_error_handler, only : MOM_error, FATAL, WARNING +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_forcing_type, only : forcing +use MOM_hor_index, only : hor_index_type +use MOM_grid, only : ocean_grid_type +use MOM_io, only : file_exists, MOM_read_data, slasher +use MOM_io, only : vardesc, var_desc, query_vardesc, stdout +use MOM_open_boundary, only : ocean_OBC_type +use MOM_restart, only : query_initialized, MOM_restart_CS +use MOM_time_manager, only : time_type +use time_interp_external_mod, only : init_external_field, time_interp_external +use MOM_tracer_registry, only : register_tracer, tracer_registry_type +use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut +use MOM_tracer_Z_init, only : tracer_Z_init +use MOM_unit_scaling, only : unit_scale_type +use MOM_variables, only : surface +use MOM_verticalGrid, only : verticalGrid_type + +implicit none ; private + +#include + +public register_CFC_cap, initialize_CFC_cap, CFC_cap_unit_tests +public CFC_cap_column_physics, CFC_cap_surface_state, CFC_cap_fluxes +public CFC_cap_stock, CFC_cap_end + +integer, parameter :: NTR = 2 !< the number of tracers in this module. + +!> The control structure for the CFC_cap tracer package +type, public :: CFC_cap_CS ; private + character(len=200) :: IC_file !< The file in which the CFC initial values can + !! be found, or an empty string for internal initilaization. + logical :: Z_IC_file !< If true, the IC_file is in Z-space. The default is false. + type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. + type(tracer_registry_type), pointer :: tr_Reg => NULL() !< A pointer to the MOM6 tracer registry + real, pointer, dimension(:,:,:) :: & + CFC11 => NULL(), & !< The CFC11 concentration [mol kg-1]. + CFC12 => NULL() !< The CFC12 concentration [mol kg-1]. + ! In the following variables a suffix of _11 refers to CFC11 and _12 to CFC12. + real :: CFC11_IC_val = 0.0 !< The initial value assigned to CFC11 [mol kg-1]. + real :: CFC12_IC_val = 0.0 !< The initial value assigned to CFC12 [mol kg-1]. + real :: CFC11_land_val = -1.0 !< The value of CFC11 used where land is masked out [mol kg-1]. + real :: CFC12_land_val = -1.0 !< The value of CFC12 used where land is masked out [mol kg-1]. + logical :: tracers_may_reinit !< If true, tracers may be reset via the initialization code + !! if they are not found in the restart files. + character(len=16) :: CFC11_name !< CFC11 variable name + character(len=16) :: CFC12_name !< CFC12 variable name + type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to regulate + !! the timing of diagnostic output. + type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< Model restart control structure + + ! The following vardesc types contain a package of metadata about each tracer. + type(vardesc) :: CFC11_desc !< A set of metadata for the CFC11 tracer + type(vardesc) :: CFC12_desc !< A set of metadata for the CFC12 tracer + !>@{ Diagnostic IDs + integer :: id_cfc11_cmor = -1, id_cfc12_cmor = -1 + !>@} +end type CFC_cap_CS + +contains + +!> Register the CFCs to be used with MOM and read the parameters +!! that are used with this tracer package +function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) + type(hor_index_type), intent(in) :: HI !< A horizontal index type structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters. + type(CFC_cap_CS), pointer :: CS !< A pointer that is set to point to the control + !! structure for this module. + type(tracer_registry_type), & + pointer :: tr_Reg !< A pointer to the tracer registry. + type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure. + + ! Local variables + character(len=40) :: mdl = "MOM_CFC_cap" ! This module's name. + character(len=200) :: inputdir ! The directory where NetCDF input files are. + ! This include declares and sets the variable "version". +#include "version_variable.h" + real, dimension(:,:,:), pointer :: tr_ptr => NULL() + real :: a11_dflt(4), a12_dflt(4) ! Default values of the various coefficients + real :: d11_dflt(4), d12_dflt(4) ! In the expressions for the solubility and + real :: e11_dflt(3), e12_dflt(3) ! Schmidt numbers. + character(len=48) :: flux_units ! The units for tracer fluxes. + character(len=48) :: dummy ! Dummy variable to store params that need to be logged here. + logical :: register_CFC_cap + integer :: isd, ied, jsd, jed, nz, m + + isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke + + if (associated(CS)) then + call MOM_error(WARNING, "register_CFC_cap called with an "// & + "associated control structure.") + return + endif + allocate(CS) + + ! Read all relevant parameters and write them to the model log. + call log_version(param_file, mdl, version, "") + call get_param(param_file, mdl, "CFC_IC_FILE", CS%IC_file, & + "The file in which the CFC initial values can be "//& + "found, or an empty string for internal initialization.", & + default=" ") + if ((len_trim(CS%IC_file) > 0) .and. (scan(CS%IC_file,'/') == 0)) then + ! Add the directory if CS%IC_file is not already a complete path. + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + CS%IC_file = trim(slasher(inputdir))//trim(CS%IC_file) + call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", CS%IC_file) + endif + call get_param(param_file, mdl, "CFC_IC_FILE_IS_Z", CS%Z_IC_file, & + "If true, CFC_IC_FILE is in depth space, not layer space", & + default=.false.) + call get_param(param_file, mdl, "TRACERS_MAY_REINIT", CS%tracers_may_reinit, & + "If true, tracers may go through the initialization code "//& + "if they are not found in the restart files. Otherwise "//& + "it is a fatal error if tracers are not found in the "//& + "restart files of a restarted run.", default=.false.) + + ! the following params are not used in this module. Instead, they are used in + ! the cap but are logged here to keep all the CFC cap params together. + call get_param(param_file, mdl, "CFC_BC_FILE", dummy, & + "The file in which the CFC-11 and CFC-12 atm concentrations can be "//& + "found (units must be parts per trillion), or an empty string for "//& + "internal BC generation (TODO).", default=" ") + if ((len_trim(dummy) > 0) .and. (scan(dummy,'/') == 0)) then + call get_param(param_file, mdl, "CFC11_VARIABLE", dummy, & + "The name of the variable representing CFC-11 in "//& + "CFC_BC_FILE.", default="CFC_11") + call get_param(param_file, mdl, "CFC12_VARIABLE", dummy, & + "The name of the variable representing CFC-12 in "//& + "CFC_BC_FILE.", default="CFC_12") + endif + + ! The following vardesc types contain a package of metadata about each tracer, + ! including, the name; units; longname; and grid information. + CS%CFC11_name = "CFC_11" ; CS%CFC12_name = "CFC_12" + CS%CFC11_desc = var_desc(CS%CFC11_name,"mol kg-1","Moles Per Unit Mass of CFC-11 in sea water", caller=mdl) + CS%CFC12_desc = var_desc(CS%CFC12_name,"mol kg-1","Moles Per Unit Mass of CFC-12 in sea water", caller=mdl) + if (GV%Boussinesq) then ; flux_units = "mol s-1" + else ; flux_units = "mol m-3 kg s-1" ; endif + + allocate(CS%CFC11(isd:ied,jsd:jed,nz)) ; CS%CFC11(:,:,:) = 0.0 + allocate(CS%CFC12(isd:ied,jsd:jed,nz)) ; CS%CFC12(:,:,:) = 0.0 + + ! This pointer assignment is needed to force the compiler not to do a copy in + ! the registration calls. Curses on the designers and implementers of F90. + tr_ptr => CS%CFC11 + ! Register CFC11 for horizontal advection, diffusion, and restarts. + call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, & + tr_desc=CS%CFC11_desc, registry_diags=.true., & + flux_units=flux_units, & + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + ! Do the same for CFC12 + tr_ptr => CS%CFC12 + call register_tracer(tr_ptr, Tr_Reg, param_file, HI, GV, & + tr_desc=CS%CFC12_desc, registry_diags=.true., & + flux_units=flux_units, & + restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit) + + CS%tr_Reg => tr_Reg + CS%restart_CSp => restart_CS + register_CFC_cap = .true. + +end function register_CFC_cap + +!> Initialize the CFC tracer fields and set up the tracer output. +subroutine initialize_CFC_cap(restart, day, G, GV, US, h, diag, OBC, CS) + logical, intent(in) :: restart !< .true. if the fields have already been + !! read from a restart file. + type(time_type), target, intent(in) :: day !< Time of the start of the run. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + type(diag_ctrl), target, intent(in) :: diag !< A structure that is used to regulate + !! diagnostic output. + type(ocean_OBC_type), pointer :: OBC !< This open boundary condition type + !! specifies whether, where, and what + !! open boundary conditions are used. + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + + ! local variables + logical :: from_file = .false. + + if (.not.associated(CS)) return + + CS%Time => day + CS%diag => diag + + if (.not.restart .or. (CS%tracers_may_reinit .and. & + .not.query_initialized(CS%CFC11, CS%CFC11_name, CS%restart_CSp))) & + call init_tracer_CFC(h, CS%CFC11, CS%CFC11_name, CS%CFC11_land_val, & + CS%CFC11_IC_val, G, GV, US, CS) + + if (.not.restart .or. (CS%tracers_may_reinit .and. & + .not.query_initialized(CS%CFC12, CS%CFC12_name, CS%restart_CSp))) & + call init_tracer_CFC(h, CS%CFC12, CS%CFC12_name, CS%CFC12_land_val, & + CS%CFC12_IC_val, G, GV, US, CS) + + + ! cmor diagnostics + ! CFC11 cmor conventions: http://clipc-services.ceda.ac.uk/dreq/u/42625c97b8fe75124a345962c4430982.html + CS%id_cfc11_cmor = register_diag_field('ocean_model', 'cfc11', diag%axesTL, day, & + 'Mole Concentration of CFC11 in Sea Water', 'mol m-3') + ! CFC12 cmor conventions: http://clipc-services.ceda.ac.uk/dreq/u/3ab8e10027d7014f18f9391890369235.html + CS%id_cfc12_cmor = register_diag_field('ocean_model', 'cfc12', diag%axesTL, day, & + 'Mole Concentration of CFC12 in Sea Water', 'mol m-3') + + if (associated(OBC)) then + ! Steal from updated DOME in the fullness of time. + ! GMM: TODO this must be coded if we intend to use this module in regional applications + endif + +end subroutine initialize_CFC_cap + +!>This subroutine initializes a tracer array. +subroutine init_tracer_CFC(h, tr, name, land_val, IC_val, G, GV, US, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: tr !< The tracer concentration array + character(len=*), intent(in) :: name !< The tracer name + real, intent(in) :: land_val !< A value the tracer takes over land + real, intent(in) :: IC_val !< The initial condition value for the tracer + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + + ! local variables + logical :: OK + integer :: i, j, k, is, ie, js, je, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + if (len_trim(CS%IC_file) > 0) then + ! Read the tracer concentrations from a netcdf file. + if (.not.file_exists(CS%IC_file, G%Domain)) & + call MOM_error(FATAL, "initialize_CFC_cap: Unable to open "//CS%IC_file) + if (CS%Z_IC_file) then + OK = tracer_Z_init(tr, h, CS%IC_file, name, G, GV, US) + if (.not.OK) then + OK = tracer_Z_init(tr, h, CS%IC_file, trim(name), G, GV, US) + if (.not.OK) call MOM_error(FATAL,"initialize_CFC_cap: "//& + "Unable to read "//trim(name)//" from "//& + trim(CS%IC_file)//".") + endif + else + call MOM_read_data(CS%IC_file, trim(name), tr, G%Domain) + endif + else + do k=1,nz ; do j=js,je ; do i=is,ie + if (G%mask2dT(i,j) < 0.5) then + tr(i,j,k) = land_val + else + tr(i,j,k) = IC_val + endif + enddo ; enddo ; enddo + endif + +end subroutine init_tracer_CFC + +!> Applies diapycnal diffusion, souces and sinks and any other column +!! tracer physics to the CFC cap tracers. CFCs are relatively simple, +!! as they are passive tracers with only a surface flux as a source. +subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, & + evap_CFL_limit, minimum_forcing_depth) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: ea !< an array to which the amount of fluid entrained + !! from the layer above during this call will be + !! added [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: eb !< an array to which the amount of fluid entrained + !! from the layer below during this call will be + !! added [H ~> m or kg m-2]. + type(forcing), intent(in) :: fluxes!< A structure containing pointers to thermodynamic + !! and tracer forcing fields. Unused fields have NULL ptrs. + real, intent(in) :: dt !< The amount of time covered by this call [T ~> s] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + real, optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can + !! be fluxed out of the top layer in a timestep [nondim] + real, optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which + !! fluxes can be applied [H ~> m or kg m-2] + + ! The arguments to this subroutine are redundant in that + ! h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1) + + ! Local variables + real, pointer, dimension(:,:,:) :: CFC11 => NULL(), CFC12 => NULL() + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2] + real :: scale_factor ! convert from [Conc. m s-1] to [Conc. kg m-2 T-1] + integer :: i, j, k, m, is, ie, js, je, nz + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + ! CFC_flux **unscaled** units are [mol m-2 s-1] which is the same as [CU kg m-2 s-1], + ! where CU = mol kg-1. However, CFC_flux has been scaled already. + ! CFC_flux **scaled** units are [CU L T-1 R]. We want [CU kg m-2 T-1] because + ! these units are what tracer_vertdiff needs. + ! Therefore, we need to convert [L R] to [kg m-2] using scale_factor. + scale_factor = US%L_to_Z * US%RZ_to_kg_m2 ! this will give [CU kg m-2 T-1], which is what verdiff wants + + if (.not.associated(CS)) return + + CFC11 => CS%CFC11 ; CFC12 => CS%CFC12 + + ! Use a tridiagonal solver to determine the concentrations after the + ! surface source is applied and diapycnal advection and diffusion occurs. + if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then + do k=1,nz ;do j=js,je ; do i=is,ie + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, CFC11, dt, fluxes, h_work, & + evap_CFL_limit, minimum_forcing_depth) + call tracer_vertdiff(h_work, ea, eb, dt, CFC11, G, GV, sfc_flux=fluxes%cfc11_flux*scale_factor) + + do k=1,nz ;do j=js,je ; do i=is,ie + h_work(i,j,k) = h_old(i,j,k) + enddo ; enddo ; enddo + call applyTracerBoundaryFluxesInOut(G, GV, CFC12, dt, fluxes, h_work, & + evap_CFL_limit, minimum_forcing_depth) + call tracer_vertdiff(h_work, ea, eb, dt, CFC12, G, GV, sfc_flux=fluxes%cfc12_flux*scale_factor) + else + call tracer_vertdiff(h_old, ea, eb, dt, CFC11, G, GV, sfc_flux=fluxes%cfc11_flux*scale_factor) + call tracer_vertdiff(h_old, ea, eb, dt, CFC12, G, GV, sfc_flux=fluxes%cfc12_flux*scale_factor) + endif + + ! If needed, write out any desired diagnostics from tracer sources & sinks here. + if (CS%id_cfc11_cmor > 0) call post_data(CS%id_cfc11_cmor, CFC11*GV%Rho0, CS%diag) + if (CS%id_cfc12_cmor > 0) call post_data(CS%id_cfc12_cmor, CFC12*GV%Rho0, CS%diag) + +end subroutine CFC_cap_column_physics + +!> Calculates the mass-weighted integral of all tracer stocks, +!! returning the number of stocks it has calculated. If the stock_index +!! is present, only the stock corresponding to that coded index is returned. +function CFC_cap_stock(h, stocks, G, GV, CS, names, units, stock_index) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + real, dimension(:), intent(out) :: stocks !< the mass-weighted integrated amount of each + !! tracer, in kg times concentration units [kg conc]. + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + character(len=*), dimension(:), intent(out) :: names !< The names of the stocks calculated. + character(len=*), dimension(:), intent(out) :: units !< The units of the stocks calculated. + integer, optional, intent(in) :: stock_index !< The coded index of a specific + !! stock being sought. + integer :: CFC_cap_stock !< The number of stocks calculated here. + + ! Local variables + real :: stock_scale ! The dimensional scaling factor to convert stocks to kg [kg H-1 L-2 ~> kg m-3 or nondim] + real :: mass ! The cell volume or mass [H L2 ~> m3 or kg] + integer :: i, j, k, is, ie, js, je, nz + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke + + CFC_cap_stock = 0 + if (.not.associated(CS)) return + + if (present(stock_index)) then ; if (stock_index > 0) then + ! Check whether this stock is available from this routine. + + ! No stocks from this routine are being checked yet. Return 0. + return + endif ; endif + + call query_vardesc(CS%CFC11_desc, name=names(1), units=units(1), caller="CFC_cap_stock") + call query_vardesc(CS%CFC12_desc, name=names(2), units=units(2), caller="CFC_cap_stock") + units(1) = trim(units(1))//" kg" ; units(2) = trim(units(2))//" kg" + + stock_scale = G%US%L_to_m**2 * GV%H_to_kg_m2 + stocks(1) = 0.0 ; stocks(2) = 0.0 + do k=1,nz ; do j=js,je ; do i=is,ie + mass = G%mask2dT(i,j) * G%areaT(i,j) * h(i,j,k) + stocks(1) = stocks(1) + CS%CFC11(i,j,k) * mass + stocks(2) = stocks(2) + CS%CFC12(i,j,k) * mass + enddo ; enddo ; enddo + stocks(1) = stock_scale * stocks(1) + stocks(2) = stock_scale * stocks(2) + + CFC_cap_stock = 2 + +end function CFC_cap_stock + +!> Extracts the ocean surface CFC concentrations and copies them to sfc_state. +subroutine CFC_cap_surface_state(sfc_state, G, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(surface), intent(inout) :: sfc_state !< A structure containing fields that + !! describe the surface state of the ocean. + type(CFC_cap_CS), pointer :: CS!< The control structure returned by a previous + !! call to register_CFC_cap. + + ! Local variables + integer :: i, j, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + if (.not.associated(CS)) return + + do j=js,je ; do i=is,ie + sfc_state%sfc_cfc11(i,j) = CS%CFC11(i,j,1) + sfc_state%sfc_cfc12(i,j) = CS%CFC12(i,j,1) + enddo ; enddo + +end subroutine CFC_cap_surface_state + +!> Orchestrates the calculation of the CFC fluxes [mol m-2 s-1], including getting the ATM +!! concentration, and calculating the solubility, Schmidt number, and gas exchange. +subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc12_atm) + type(ocean_grid_type), intent(in ) :: G !< The ocean's grid structure. + type(surface), intent(in ) :: sfc_state !< A structure containing fields + !! that describe the surface state of the ocean. + type(forcing), intent(inout) :: fluxes !< A structure containing pointers + !! to thermodynamic and tracer forcing fields. Unused fields + !! have NULL ptrs. + real, intent(in ) :: Rho0 !< The mean ocean density [kg m-3] + type(time_type), intent(in ) :: Time !< The time of the fluxes, used for interpolating the + !! CFC's concentration in the atmosphere. + integer, optional, intent(inout):: id_cfc11_atm !< id number for time_interp_external. + integer, optional, intent(inout):: id_cfc12_atm !< id number for time_interp_external. + + ! Local variables + real, dimension(SZI_(G),SZJ_(G)) :: & + CFC11_Csurf, & ! The CFC-11 surface concentrations times the Schmidt number term [mol kg-1]. + CFC12_Csurf, & ! The CFC-12 surface concentrations times the Schmidt number term [mol kg-1]. + CFC11_alpha, & ! The CFC-11 solubility [mol kg-1 atm-1]. + CFC12_alpha, & ! The CFC-12 solubility [mol kg-1 atm-1]. + kw_wo_sc_no_term, & ! gas transfer velocity, without the Schmidt number term [m s-1]. + cair, & ! The surface gas concentration in equilibrium with the atmosphere (saturation concentration) + ! [mol kg-1]. + cfc11_atm, & !< CFC11 concentration in the atmopshere [pico mol/mol] + cfc12_atm !< CFC11 concentration in the atmopshere [pico mol/mol] + real :: ta ! Absolute sea surface temperature [hectoKelvin] + real :: sal ! Surface salinity [PSU]. + real :: alpha_11 ! The solubility of CFC 11 [mol kg-1 atm-1]. + real :: alpha_12 ! The solubility of CFC 12 [mol kg-1 atm-1]. + real :: sc_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12. + real :: sc_no_term ! A term related to the Schmidt number. + real :: kw_coeff ! A coefficient used to scale the piston velocity [L T-1 ~> m s-1] + real, parameter :: pa_to_atm = 9.8692316931427e-6 ! factor for converting from Pa to atm. + integer :: i, j, m, is, ie, js, je + + is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec + + ! CFC11 ATM concentration + if (present(id_cfc11_atm) .and. (id_cfc11_atm /= -1)) then + call time_interp_external(id_cfc11_atm, Time, cfc11_atm) + ! convert from ppt (pico mol/mol) to mol/mol + cfc11_atm = cfc11_atm * 1.0e-12 + else + ! TODO: create cfc11_atm internally + call MOM_error(FATAL, "CFC_cap_fluxes: option to create cfc11_atm internally" //& + "has not been implemented yet.") + endif + + ! CFC12 ATM concentration + if (present(id_cfc12_atm) .and. (id_cfc12_atm /= -1)) then + call time_interp_external(id_cfc12_atm, Time, cfc12_atm) + ! convert from ppt (pico mol/mol) to mol/mol + cfc12_atm = cfc12_atm * 1.0e-12 + else + ! TODO: create cfc11_atm internally + call MOM_error(FATAL, "CFC_cap_fluxes: option to create cfc12_atm internally" //& + "has not been implemented yet.") + endif + + do j=js,je ; do i=is,ie + ! ta in hectoKelvin + ta = max(0.01, (sfc_state%SST(i,j) + 273.15) * 0.01) + sal = sfc_state%SSS(i,j) + + ! Calculate solubilities + call get_solubility(alpha_11, alpha_12, ta, sal , G%mask2dT(i,j)) + + ! Calculate Schmidt numbers using coefficients given by + ! Wanninkhof (2014); doi:10.4319/lom.2014.12.351. + call comp_CFC_schmidt(sfc_state%SST(i,j), sc_11, sc_12) + sc_no_term = sqrt(660.0 / sc_11) + + CFC11_alpha(i,j) = alpha_11 * sc_no_term + CFC11_Csurf(i,j) = sfc_state%sfc_CFC11(i,j) * sc_no_term + sc_no_term = sqrt(660.0 / sc_12) + CFC12_alpha(i,j) = alpha_12 * sc_no_term + CFC12_Csurf(i,j) = sfc_state%sfc_CFC12(i,j) * sc_no_term + + !--------------------------------------------------------------------- + ! Gas exchange/piston velocity parameter + !--------------------------------------------------------------------- + ! From a = 0.251 cm/hr s^2/m^2 in Wannikhof 2014 + ! 6.97e-07 is used to convert from cm/hr to m/s, kw = m/s [L/T] + kw_coeff = 6.97e-07 * G%US%m_to_L * G%US%T_to_s + + kw_wo_sc_no_term(i,j) = kw_coeff * ((1.0 - fluxes%ice_fraction(i,j))*fluxes%u10_sqr(i,j)) + + ! air concentrations and cfcs BC's fluxes + ! CFC flux units: mol kg-1 s-1 kg m-2 + cair(i,j) = pa_to_atm * CFC11_alpha(i,j) * cfc11_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc11_flux(i,j) = kw_wo_sc_no_term(i,j) * (cair(i,j) - CFC11_Csurf(i,j)) * Rho0 + cair(i,j) = pa_to_atm * CFC12_alpha(i,j) * cfc12_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc12_flux(i,j) = kw_wo_sc_no_term(i,j) * (cair(i,j) - CFC12_Csurf(i,j)) * Rho0 + enddo ; enddo + +end subroutine CFC_cap_fluxes + +!> Calculates the CFC's solubility function following Warner and Weiss (1985) DSR, vol 32. +subroutine get_solubility(alpha_11, alpha_12, ta, sal , mask) + real, intent(inout) :: alpha_11 !< The solubility of CFC 11 [mol kg-1 atm-1] + real, intent(inout) :: alpha_12 !< The solubility of CFC 12 [mol kg-1 atm-1] + real, intent(in ) :: ta !< Absolute sea surface temperature [hectoKelvin] + real, intent(in ) :: sal !< Surface salinity [PSU]. + real, intent(in ) :: mask !< ocean mask + + ! Local variables + real :: d11_dflt(4), d12_dflt(4) ! values of the various coefficients + real :: e11_dflt(3), e12_dflt(3) ! in the expressions for the solubility + real :: d1_11, d1_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [nondim] + real :: d2_11, d2_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [hectoKelvin-1] + real :: d3_11, d3_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [log(hectoKelvin)-1] + real :: d4_11, d4_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [hectoKelvin-2] + real :: e1_11, e1_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-1] + real :: e2_11, e2_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-1 hectoKelvin-1] + real :: e3_11, e3_12 ! Coefficients for calculating CFC11 and CFC12 solubilities [PSU-2 hectoKelvin-2] + real :: factor ! factor to use in the solubility conversion + + !----------------------------------------------------------------------- + ! Solubility coefficients for alpha in mol/(kg atm) for CFC11 (_11) and CFC12 (_12) + ! from Table 5 in Warner and Weiss (1985) DSR, vol 32. + !----------------------------------------------------------------------- + d11_dflt(:) = (/ -232.0411, 322.5546, 120.4956, -1.39165 /) + e11_dflt(:) = (/ -0.146531, 0.093621, -0.0160693 /) + d12_dflt(:) = (/ -220.2120, 301.8695, 114.8533, -1.39165 /) + e12_dflt(:) = (/ -0.147718, 0.093175, -0.0157340 /) + + d1_11 = d11_dflt(1) + d2_11 = d11_dflt(2) + d3_11 = d11_dflt(3) + d4_11 = d11_dflt(4) + + e1_11 = e11_dflt(1) + e2_11 = e11_dflt(2) + e3_11 = e11_dflt(3) + + d1_12 = d12_dflt(1) + d2_12 = d12_dflt(2) + d3_12 = d12_dflt(3) + d4_12 = d12_dflt(4) + + e1_12 = e12_dflt(1) + e2_12 = e12_dflt(2) + e3_12 = e12_dflt(3) + + ! Calculate solubilities using Warner and Weiss (1985) DSR, vol 32. + ! The following is from Eq. 9 in Warner and Weiss (1985) + ! The factor 1.0e+03 is for the conversion from mol/(l * atm) to mol/(m3 * atm) + ! The factor 1.e-09 converts from mol/(l * atm) to mol/(m3 * pptv) + factor = 1.0 + alpha_11 = exp(d1_11 + d2_11/ta + d3_11*log(ta) + d4_11*ta**2 +& + sal * ((e3_11 * ta + e2_11) * ta + e1_11)) * & + factor * mask + alpha_12 = exp(d1_12 + d2_12/ta + d3_12*log(ta) + d4_12*ta**2 +& + sal * ((e3_12 * ta + e2_12) * ta + e1_12)) * & + factor * mask + +end subroutine get_solubility + + +!> Compute Schmidt numbers of CFCs following Wanninkhof (2014); doi:10.4319/lom.2014.12.351 +!! Range of validity of fit is -2:40. +subroutine comp_CFC_schmidt(sst_in, cfc11_sc, cfc12_sc) + real, intent(in) :: sst_in !< The sea surface temperature [degC]. + real, intent(inout) :: cfc11_sc !< Schmidt number of CFC11 [nondim]. + real, intent(inout) :: cfc12_sc !< Schmidt number of CFC12 [nondim]. + + !local variables + real , parameter :: a_11 = 3579.2 + real , parameter :: b_11 = -222.63 + real , parameter :: c_11 = 7.5749 + real , parameter :: d_11 = -0.14595 + real , parameter :: e_11 = 0.0011874 + real , parameter :: a_12 = 3828.1 + real , parameter :: b_12 = -249.86 + real , parameter :: c_12 = 8.7603 + real , parameter :: d_12 = -0.1716 + real , parameter :: e_12 = 0.001408 + real :: sst + + + ! clip SST to avoid bad values + sst = MAX(-2.0, MIN(40.0, sst_in)) + cfc11_sc = a_11 + sst * (b_11 + sst * (c_11 + sst * (d_11 + sst * e_11))) + cfc12_sc = a_12 + sst * (b_12 + sst * (c_12 + sst * (d_12 + sst * e_12))) + +end subroutine comp_CFC_schmidt + +!> Deallocate any memory associated with the CFC cap tracer package +subroutine CFC_cap_end(CS) + type(CFC_cap_CS), pointer :: CS !< The control structure returned by a + !! previous call to register_CFC_cap. + ! local variables + integer :: m + + if (associated(CS)) then + if (associated(CS%CFC11)) deallocate(CS%CFC11) + if (associated(CS%CFC12)) deallocate(CS%CFC12) + + deallocate(CS) + endif +end subroutine CFC_cap_end + +!> Unit tests for the CFC cap module. +logical function CFC_cap_unit_tests(verbose) + logical, intent(in) :: verbose !< If true, output additional + !! information for debugging unit tests + + ! Local variables + real :: dummy1, dummy2, ta, sal + character(len=120) :: test_name ! Title of the unit test + + CFC_cap_unit_tests = .false. + write(stdout,*) '==== MOM_CFC_cap =======================' + + ! test comp_CFC_schmidt, Table 1 in Wanninkhof (2014); doi:10.4319/lom.2014.12.351 + test_name = 'Schmidt number calculation' + call comp_CFC_schmidt(20.0, dummy1, dummy2) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy1, 1179.0, 0.5) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy2, 1188.0, 0.5) + + if (.not. CFC_cap_unit_tests) write(stdout,'(2x,a)') "Passed "//test_name + + test_name = 'Solubility function, SST = 1.0 C, and SSS = 10 psu' + ta = max(0.01, (1.0 + 273.15) * 0.01); sal = 10. + ! cfc1 = 3.238 10-2 mol kg-1 atm-1 + ! cfc2 = 7.943 10-3 mol kg-1 atm-1 + call get_solubility(dummy1, dummy2, ta, sal , 1.0) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy1, 3.238e-2, 5.0e-6) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy2, 7.943e-3, 5.0e-6) + + if (.not. CFC_cap_unit_tests) write(stdout,'(2x,a)')"Passed "//test_name + + test_name = 'Solubility function, SST = 20.0 C, and SSS = 35 psu' + ta = max(0.01, (20.0 + 273.15) * 0.01); sal = 35. + ! cfc1 = 0.881 10-2 mol kg-1 atm-1 + ! cfc2 = 2.446 10-3 mol kg-1 atm-1 + call get_solubility(dummy1, dummy2, ta, sal , 1.0) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy1, 8.8145e-3, 5.0e-8) + CFC_cap_unit_tests = CFC_cap_unit_tests .or. & + compare_values(verbose, test_name, dummy2, 2.4462e-3, 5.0e-8) + if (.not. CFC_cap_unit_tests) write(stdout,'(2x,a)')"Passed "//test_name + +end function CFC_cap_unit_tests + +!> Test that ans and calc are approximately equal by computing the difference +!! and comparing it against limit. +logical function compare_values(verbose, test_name, calc, ans, limit) + logical, intent(in) :: verbose !< If true, write results to stdout + character(len=80), intent(in) :: test_name !< Brief description of the unit test + real, intent(in) :: calc !< computed value + real, intent(in) :: ans !< correct value + real, intent(in) :: limit !< value above which test fails + + ! Local variables + real :: diff + + diff = ans - calc + + compare_values = .false. + if (diff > limit ) then + compare_values = .true. + write(stdout,*) "CFC_cap_unit_tests, UNIT TEST FAILED: ", test_name + write(stdout,10) calc, ans + elseif (verbose) then + write(stdout,10) calc, ans + endif + +10 format("calc=",f20.16," ans",f20.16) +end function compare_values + +!> \namespace mom_CFC_cap +!! +!! This module contains the code that is needed to simulate +!! CFC-11 and CFC-12 using atmospheric and sea ice variables +!! provided via cap (only NUOPC cap is implemented so far). + +end module MOM_CFC_cap diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90 index 26ef197ae2..439e9b5396 100644 --- a/src/tracer/MOM_tracer_flow_control.F90 +++ b/src/tracer/MOM_tracer_flow_control.F90 @@ -42,6 +42,9 @@ module MOM_tracer_flow_control use MOM_OCMIP2_CFC, only : register_OCMIP2_CFC, initialize_OCMIP2_CFC, flux_init_OCMIP2_CFC use MOM_OCMIP2_CFC, only : OCMIP2_CFC_column_physics, OCMIP2_CFC_surface_state use MOM_OCMIP2_CFC, only : OCMIP2_CFC_stock, OCMIP2_CFC_end, OCMIP2_CFC_CS +use MOM_CFC_cap, only : register_CFC_cap, initialize_CFC_cap +use MOM_CFC_cap, only : CFC_cap_column_physics, CFC_cap_surface_state +use MOM_CFC_cap, only : CFC_cap_stock, CFC_cap_end, CFC_cap_CS use oil_tracer, only : register_oil_tracer, initialize_oil_tracer use oil_tracer, only : oil_tracer_column_physics, oil_tracer_surface_state use oil_tracer, only : oil_stock, oil_tracer_end, oil_tracer_CS @@ -80,6 +83,7 @@ module MOM_tracer_flow_control logical :: use_oil = .false. !< If true, use the oil tracer package logical :: use_advection_test_tracer = .false. !< If true, use the advection_test_tracer package logical :: use_OCMIP2_CFC = .false. !< If true, use the OCMIP2_CFC tracer package + logical :: use_CFC_cap = .false. !< If true, use the CFC_cap tracer package logical :: use_MOM_generic_tracer = .false. !< If true, use the MOM_generic_tracer packages logical :: use_pseudo_salt_tracer = .false. !< If true, use the psuedo_salt tracer package logical :: use_boundary_impulse_tracer = .false. !< If true, use the boundary impulse tracer package @@ -94,6 +98,7 @@ module MOM_tracer_flow_control type(oil_tracer_CS), pointer :: oil_tracer_CSp => NULL() type(advection_test_tracer_CS), pointer :: advection_test_tracer_CSp => NULL() type(OCMIP2_CFC_CS), pointer :: OCMIP2_CFC_CSp => NULL() + type(CFC_cap_CS), pointer :: CFC_cap_CSp => NULL() type(MOM_generic_tracer_CS), pointer :: MOM_generic_tracer_CSp => NULL() type(pseudo_salt_tracer_CS), pointer :: pseudo_salt_tracer_CSp => NULL() type(boundary_impulse_tracer_CS), pointer :: boundary_impulse_tracer_CSp => NULL() @@ -114,7 +119,7 @@ subroutine call_tracer_flux_init(verbosity) type(param_file_type) :: param_file ! A structure to parse for run-time parameters character(len=40) :: mdl = "call_tracer_flux_init" ! This module's name. - logical :: use_OCMIP_CFCs, use_MOM_generic_tracer + logical :: use_OCMIP_CFCs, use_MOM_generic_tracer, use_CFC_caps ! Determine which tracer routines with tracer fluxes are to be called. Note ! that not every tracer package is required to have a flux_init call. @@ -193,6 +198,9 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) call get_param(param_file, mdl, "USE_OCMIP2_CFC", CS%use_OCMIP2_CFC, & "If true, use the MOM_OCMIP2_CFC tracer package.", & default=.false.) + call get_param(param_file, mdl, "USE_CFC_CAP", CS%use_CFC_cap, & + "If true, use the MOM_CFC_cap tracer package.", & + default=.false.) call get_param(param_file, mdl, "USE_generic_tracer", CS%use_MOM_generic_tracer, & "If true and _USE_GENERIC_TRACER is defined as a "//& "preprocessor macro, use the MOM_generic_tracer packages.", & @@ -237,6 +245,9 @@ subroutine call_tracer_register(HI, GV, US, param_file, CS, tr_Reg, restart_CS) if (CS%use_OCMIP2_CFC) CS%use_OCMIP2_CFC = & register_OCMIP2_CFC(HI, GV, param_file, CS%OCMIP2_CFC_CSp, & tr_Reg, restart_CS) + if (CS%use_CFC_cap) CS%use_CFC_cap = & + register_CFC_cap(HI, GV, param_file, CS%CFC_cap_CSp, & + tr_Reg, restart_CS) if (CS%use_MOM_generic_tracer) CS%use_MOM_generic_tracer = & register_MOM_generic_tracer(HI, GV, param_file, CS%MOM_generic_tracer_CSp, & tr_Reg, restart_CS) @@ -317,6 +328,9 @@ subroutine tracer_flow_control_init(restart, day, G, GV, US, h, param_file, diag if (CS%use_OCMIP2_CFC) & call initialize_OCMIP2_CFC(restart, day, G, GV, US, h, diag, OBC, CS%OCMIP2_CFC_CSp, & sponge_CSp) + if (CS%use_CFC_cap) & + call initialize_CFC_cap(restart, day, G, GV, US, h, diag, OBC, CS%CFC_cap_CSp) + if (CS%use_MOM_generic_tracer) & call initialize_MOM_generic_tracer(restart, day, G, GV, US, h, param_file, diag, OBC, & CS%MOM_generic_tracer_CSp, sponge_CSp, ALE_sponge_CSp) @@ -462,6 +476,11 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, G, GV, US, CS%OCMIP2_CFC_CSp, & evap_CFL_limit=evap_CFL_limit, & minimum_forcing_depth=minimum_forcing_depth) + if (CS%use_CFC_cap) & + call CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, CS%CFC_cap_CSp, & + evap_CFL_limit=evap_CFL_limit, & + minimum_forcing_depth=minimum_forcing_depth) if (CS%use_MOM_generic_tracer) then if (US%QRZ_T_to_W_m2 /= 1.0) call MOM_error(FATAL, "MOM_generic_tracer_column_physics "//& "has not been written to permit dimensionsal rescaling. Set all 4 of the "//& @@ -516,6 +535,9 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, if (CS%use_OCMIP2_CFC) & call OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, & G, GV, US, CS%OCMIP2_CFC_CSp) + if (CS%use_CFC_cap) & + call CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, & + G, GV, US, CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) then if (US%QRZ_T_to_W_m2 /= 1.0) call MOM_error(FATAL, "MOM_generic_tracer_column_physics "//& "has not been written to permit dimensionsal rescaling. Set all 4 of the "//& @@ -624,6 +646,12 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, CS, stock_names, stock_uni set_pkg_name, max_ns, ns_tot, stock_names, stock_units) endif + if (CS%use_CFC_cap) then + ns = CFC_cap_stock(h, values, G, GV, CS%CFC_cap_CSp, names, units, stock_index) + call store_stocks("MOM_CFC_cap", ns, names, units, values, index, stock_values, & + set_pkg_name, max_ns, ns_tot, stock_names, stock_units) + endif + if (CS%use_advection_test_tracer) then ns = advection_test_stock( h, values, G, GV, CS%advection_test_tracer_CSp, & names, units, stock_index ) @@ -753,6 +781,8 @@ subroutine call_tracer_surface_state(sfc_state, h, G, GV, CS) call advection_test_tracer_surface_state(sfc_state, h, G, GV, CS%advection_test_tracer_CSp) if (CS%use_OCMIP2_CFC) & call OCMIP2_CFC_surface_state(sfc_state, h, G, GV, CS%OCMIP2_CFC_CSp) + if (CS%use_CFC_cap) & + call CFC_cap_surface_state(sfc_state, G, CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) & call MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS%MOM_generic_tracer_CSp) @@ -772,6 +802,7 @@ subroutine tracer_flow_control_end(CS) if (CS%use_oil) call oil_tracer_end(CS%oil_tracer_CSp) if (CS%use_advection_test_tracer) call advection_test_tracer_end(CS%advection_test_tracer_CSp) if (CS%use_OCMIP2_CFC) call OCMIP2_CFC_end(CS%OCMIP2_CFC_CSp) + if (CS%use_CFC_cap) call CFC_cap_end(CS%CFC_cap_CSp) if (CS%use_MOM_generic_tracer) call end_MOM_generic_tracer(CS%MOM_generic_tracer_CSp) if (CS%use_pseudo_salt_tracer) call pseudo_salt_tracer_end(CS%pseudo_salt_tracer_CSp) if (CS%use_boundary_impulse_tracer) call boundary_impulse_tracer_end(CS%boundary_impulse_tracer_CSp)