diff --git a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 index e2a557daff..7b32270b4c 100644 --- a/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 +++ b/config_src/drivers/mct_cap/mom_surface_forcing_mct.F90 @@ -1245,7 +1245,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "If true, use a 2-dimensional gustiness supplied from "//& "an input file", default=.false.) call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, & - "The background gustiness in the winds.", units="Pa", default=0.0) + "The background gustiness in the winds.", units="Pa", default=0.0, & + scale=US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z) if (CS%read_gust_2d) then call get_param(param_file, mdl, "GUST_2D_FILE", gust_file, & "The file in which the wind gustiness is found in "//& diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index 0090468a05..ef3d9867b9 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -1113,7 +1113,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) k = k + 1 ! Increment position within gindex if (mask(k) /= 0) then mesh_areas(k) = dataPtr_mesh_areas(k) - model_areas(k) = ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2 + model_areas(k) = ocean_grid%US%L_to_m**2 * ocean_grid%AreaT(i,j) / ocean_grid%Rad_Earth**2 mod2med_areacor(k) = model_areas(k) / mesh_areas(k) med2mod_areacor(k) = mesh_areas(k) / model_areas(k) end if 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 c83e32711c..5d1070954b 100644 --- a/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_surface_forcing_nuopc.F90 @@ -576,7 +576,7 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G, ! 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) + call CFC_cap_fluxes(fluxes, sfc_state, G, US, CS%Rho0, Time, CS%id_cfc11_atm, CS%id_cfc11_atm) endif if (associated(IOB%salt_flux)) then @@ -1396,7 +1396,9 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, restore_salt, "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) + CS%CFC_BC_file = trim(CS%inputdir) // trim(CS%CFC_BC_file) + endif + if (len_trim(CS%CFC_BC_file) > 0) then 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.) diff --git a/config_src/drivers/solo_driver/MESO_surface_forcing.F90 b/config_src/drivers/solo_driver/MESO_surface_forcing.F90 index 679f147797..7c778d9753 100644 --- a/config_src/drivers/solo_driver/MESO_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/MESO_surface_forcing.F90 @@ -243,7 +243,8 @@ subroutine MESO_surface_forcing_init(Time, G, US, param_file, diag, CS) "parameters from vertical units of m to kg m-2.", & units="kg m-3", default=1035.0, scale=US%kg_m3_to_R) call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, & - "The background gustiness in the winds.", units="Pa", default=0.0) + "The background gustiness in the winds.", units="Pa", default=0.0, & + scale=US%kg_m3_to_R*US%m_s_to_L_T**2*US%L_to_Z) call get_param(param_file, mdl, "RESTOREBUOY", CS%restorebuoy, & "If true, the buoyancy fluxes drive the model back "//& diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 7d8375a8af..77eda959ad 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -185,8 +185,8 @@ module MOM_forcing_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]. + cfc11_flux => NULL(), & !< flux of cfc_11 into the ocean [CU Z T-1 kg m-3 = mol Z T-1 m-3 ~> mol m-2 s-1]. + cfc12_flux => NULL(), & !< flux of cfc_12 into the ocean [CU Z T-1 kg m-3 = mol Z T-1 m-3 ~> 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] @@ -1101,9 +1101,13 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) 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) + call hchksum(fluxes%u10_sqr, mesg//" fluxes%u10_sqr", G%HI, haloshift=hshift , scale=US%L_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%cfc11_flux)) & + call hchksum(fluxes%cfc11_flux, mesg//" fluxes%cfc11_flux", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) + if (associated(fluxes%cfc12_flux)) & + call hchksum(fluxes%cfc12_flux, mesg//" fluxes%cfc12_flux", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T) 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)) & @@ -1117,7 +1121,7 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) call hchksum(fluxes%frunoff, mesg//" fluxes%frunoff", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s) if (associated(fluxes%heat_content_lrunoff)) & call hchksum(fluxes%heat_content_lrunoff, mesg//" fluxes%heat_content_lrunoff", G%HI, & - haloshift=hshift, scale=US%RZ_T_to_kg_m2s) + haloshift=hshift, scale=US%QRZ_T_to_W_m2) if (associated(fluxes%heat_content_frunoff)) & call hchksum(fluxes%heat_content_frunoff, mesg//" fluxes%heat_content_frunoff", G%HI, & haloshift=hshift, scale=US%QRZ_T_to_W_m2) @@ -1314,14 +1318,14 @@ subroutine register_forcing_type_diags(Time, diag, US, use_temperature, handles, 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,& + conversion= US%Z_to_m*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,& + conversion= US%Z_to_m*US%s_to_T,& cmor_field_name='fgcfc12', & cmor_long_name='Surface Downward CFC12 Flux', & cmor_standard_name='surface_downward_cfc12_flux') diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 762b2edaea..25946eb631 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -81,7 +81,7 @@ module MOM_MEKE real :: MEKE_advection_factor !< A scaling in front of the advection of MEKE [nondim] real :: MEKE_topographic_beta !< Weight for how much topographic beta is considered !! when computing beta in Rhines scale [nondim] - real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [s-1]. + real :: MEKE_restoring_rate !< Inverse of the timescale used to nudge MEKE toward its equilibrium value [T-1 ~> s-1]. logical :: kh_flux_enabled !< If true, lateral diffusive MEKE flux is enabled. logical :: initialize !< If True, invokes a steady state solver to calculate MEKE. @@ -340,6 +340,10 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h enddo ; enddo endif + if (CS%debug) then + call hchksum(src, "MEKE src", G%HI, haloshift=0, scale=US%L_to_m**2*US%s_to_T**3) + endif + ! Increase EKE by a full time-steps worth of source !$OMP parallel do default(shared) do j=js,je ; do i=is,ie @@ -534,6 +538,10 @@ subroutine step_forward_MEKE(MEKE, h, SN_u, SN_v, visc, dt, G, GV, US, CS, hu, h endif endif ! MEKE_KH>=0 + if (CS%debug) then + call hchksum(MEKE%MEKE, "MEKE post-update MEKE", G%HI, haloshift=0, scale=US%L_T_to_m_s**2) + endif + ! do j=js,je ; do i=is,ie ! MEKE%MEKE(i,j) = MAX(MEKE%MEKE(i,j),0.0) ! enddo ; enddo @@ -688,7 +696,7 @@ subroutine MEKE_equilibrium(CS, MEKE, G, GV, US, SN_u, SN_v, drag_rate_visc, I_m SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) if (CS%MEKE_equilibrium_alt) then - MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 + MEKE%MEKE(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_L*G%bathyT(i,j))**2 / cd2 else FatH = 0.25*((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) + & (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1))) ! Coriolis parameter at h points @@ -824,7 +832,7 @@ subroutine MEKE_equilibrium_restoring(CS, G, US, SN_u, SN_v) ! SN = 0.25*max( (SN_u(I,j) + SN_u(I-1,j)) + (SN_v(i,J) + SN_v(i,J-1)), 0.) ! This avoids extremes values in equilibrium solution due to bad values in SN_u, SN_v SN = min(SN_u(I,j), SN_u(I-1,j), SN_v(i,J), SN_v(i,J-1)) - CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_m*G%bathyT(i,j))**2 / cd2 + CS%equilibrium_value(i,j) = (CS%MEKE_GEOMETRIC_alpha * SN * US%Z_to_L*G%bathyT(i,j))**2 / cd2 enddo ; enddo if (CS%id_MEKE_equilibrium>0) call post_data(CS%id_MEKE_equilibrium, CS%equilibrium_value, CS%diag) @@ -957,7 +965,7 @@ subroutine MEKE_lengthScales_0d(CS, US, area, beta, depth, Rd_dx, SN, EKE, & ! Z Leady = 0. endif if (CS%use_min_lscale) then - LmixScale = 1.e7 + LmixScale = 1.e7*US%m_to_L if (CS%aDeform*Ldeform > 0.) LmixScale = min(LmixScale,CS%aDeform*Ldeform) if (CS%aFrict *Lfrict > 0.) LmixScale = min(LmixScale,CS%aFrict *Lfrict) if (CS%aRhines*Lrhines > 0.) LmixScale = min(LmixScale,CS%aRhines*Lrhines) @@ -1069,7 +1077,7 @@ logical function MEKE_init(Time, G, US, param_file, diag, CS, MEKE, restart_CS) if (CS%MEKE_equilibrium_restoring) then call get_param(param_file, mdl, "MEKE_RESTORING_TIMESCALE", MEKE_restoring_timescale, & "The timescale used to nudge MEKE toward its equilibrium value.", units="s", & - default=1e6, scale=US%T_to_s) + default=1e6, scale=US%s_to_T) CS%MEKE_restoring_rate = 1.0 / MEKE_restoring_timescale endif diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index fb70f5d679..dabeb27159 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -506,10 +506,10 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: slope_x !< Zonal isoneutral slope real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: N2_u !< Buoyancy (Brunt-Vaisala) frequency - !! at u-points [T-2 ~> s-2] + !! at u-points [L2 Z-2 T-2 ~> s-2] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: slope_y !< Meridional isoneutral slope real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: N2_v !< Buoyancy (Brunt-Vaisala) frequency - !! at v-points [T-2 ~> s-2] + !! at v-points [L2 Z-2 T-2 ~> s-2] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(VarMix_CS), pointer :: CS !< Variable mixing coefficients type(ocean_OBC_type), optional, pointer :: OBC !< Open boundaries control structure. @@ -654,9 +654,10 @@ subroutine calc_Visbeck_coeffs_old(h, slope_x, slope_y, N2_u, N2_v, G, GV, US, C endif if (CS%debug) then - call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, haloshift=1) + call uvchksum("calc_Visbeck_coeffs_old slope_[xy]", slope_x, slope_y, G%HI, & + scale=US%Z_to_L, haloshift=1) call uvchksum("calc_Visbeck_coeffs_old N2_u, N2_v", N2_u, N2_v, G%HI, & - scale=US%s_to_T**2, scalar_pair=.true.) + scale=US%L_to_Z**2 * US%s_to_T**2, scalar_pair=.true.) call uvchksum("calc_Visbeck_coeffs_old SN_[uv]", CS%SN_u, CS%SN_v, G%HI, & scale=US%s_to_T, scalar_pair=.true.) endif @@ -1484,7 +1485,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) allocate(CS%Depth_fn_v(isd:ied,JsdB:JedB)) ; CS%Depth_fn_v(:,:) = 0.0 call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_H0", CS%depth_scaled_khth_h0, & "The depth above which KHTH is scaled away.",& - units="m", default=1000.) + units="m", scale=US%m_to_Z, default=1000.) call get_param(param_file, mdl, "DEPTH_SCALED_KHTH_EXP", CS%depth_scaled_khth_exp, & "The exponent used in the depth dependent scaling function for KHTH.",& units="nondim", default=3.0) diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index da62ffc6b7..4d038ed31e 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -433,7 +433,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp call hchksum(e, "thickness_diffuse_1 e", G%HI, haloshift=1, scale=US%Z_to_m) if (use_stored_slopes) then call uvchksum("VarMix%slope_[xy]", VarMix%slope_x, VarMix%slope_y, & - G%HI, haloshift=0) + G%HI, haloshift=0, scale=US%Z_to_L) endif if (associated(tv%eqn_of_state)) then call hchksum(tv%T, "thickness_diffuse T", G%HI, haloshift=1) @@ -505,7 +505,7 @@ subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp enddo do j=js,je ; do i=is,ie - MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(1.0,G%bathyT(i,j)) + MEKE%Kh_diff(i,j) = MEKE%Kh_diff(i,j) / MAX(GV%m_to_H,GV%Z_to_H*G%bathyT(i,j)) enddo ; enddo endif diff --git a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 index 083f5ed000..dec8fad571 100644 --- a/src/parameterizations/vertical/MOM_CVMix_KPP.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_KPP.F90 @@ -644,7 +644,7 @@ subroutine KPP_calculate(CS, G, GV, US, h, uStar, & if (CS%debug) then call hchksum(h, "KPP in: h",G%HI,haloshift=0, scale=GV%H_to_m) call hchksum(uStar, "KPP in: uStar",G%HI,haloshift=0, scale=US%Z_to_m*US%s_to_T) - call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0) + call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0, scale=US%L_to_m**2*US%s_to_T**3) call hchksum(Kt, "KPP in: Kt",G%HI,haloshift=0, scale=US%Z2_T_to_m2_s) call hchksum(Ks, "KPP in: Ks",G%HI,haloshift=0, scale=US%Z2_T_to_m2_s) endif @@ -967,8 +967,8 @@ subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFl if (CS%debug) then call hchksum(Salt, "KPP in: S",G%HI,haloshift=0) call hchksum(Temp, "KPP in: T",G%HI,haloshift=0) - call hchksum(u, "KPP in: u",G%HI,haloshift=0) - call hchksum(v, "KPP in: v",G%HI,haloshift=0) + call hchksum(u, "KPP in: u",G%HI,haloshift=0,scale=US%L_T_to_m_s) + call hchksum(v, "KPP in: v",G%HI,haloshift=0,scale=US%L_T_to_m_s) endif call cpu_clock_begin(id_clock_KPP_compute_BLD) diff --git a/src/parameterizations/vertical/MOM_CVMix_conv.F90 b/src/parameterizations/vertical/MOM_CVMix_conv.F90 index 8b007a0b11..a615c9f40b 100644 --- a/src/parameterizations/vertical/MOM_CVMix_conv.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_conv.F90 @@ -288,8 +288,8 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux) ! call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) ! if (CS%id_kv_conv > 0) & ! call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T) - if (present(Kd)) call hchksum(Kv, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T) - if (present(Kv)) call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T) + if (present(Kd)) call hchksum(Kd, "MOM_CVMix_conv: Kd", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) + if (present(Kv)) call hchksum(Kv, "MOM_CVMix_conv: Kv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s) endif ! send diagnostics to post_data diff --git a/src/tracer/MOM_CFC_cap.F90 b/src/tracer/MOM_CFC_cap.F90 index 58ba419c62..e1770b0d52 100644 --- a/src/tracer/MOM_CFC_cap.F90 +++ b/src/tracer/MOM_CFC_cap.F90 @@ -84,11 +84,7 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) ! 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. + character(len=200) :: dummy ! Dummy variable to store params that need to be logged here. logical :: register_CFC_cap integer :: isd, ied, jsd, jed, nz, m @@ -121,6 +117,12 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) "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.) + call get_param(param_file, mdl, "CFC11_IC_VAL", CS%CFC11_IC_val, & + "Value that CFC_11 is set to when it is not read from a file.", & + units="mol kg-1", default=0.0) + call get_param(param_file, mdl, "CFC12_IC_VAL", CS%CFC12_IC_val, & + "Value that CFC_12 is set to when it is not read from a file.", & + units="mol kg-1", default=0.0) ! 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. @@ -129,6 +131,12 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) "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 + ! Add the directory if dummy is not already a complete path. + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + dummy = trim(slasher(inputdir))//trim(dummy) + call log_param(param_file, mdl, "INPUTDIR/CFC_IC_FILE", dummy) + endif + if (len_trim(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") @@ -142,8 +150,6 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) 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 @@ -154,13 +160,11 @@ function register_CFC_cap(HI, GV, param_file, CS, tr_Reg, restart_CS) ! 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 @@ -302,20 +306,18 @@ subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, C ! 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] + real :: scale_factor ! convert from cfc1[12]_flux to units of sfc_flux in tracer_vertdiff 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 + ! set factor to convert from cfc1[12]_flux units to tracer_vertdiff argument sfc_flux units + ! cfc1[12]_flux units are CU Z T-1 kg m-3 + ! tracer_vertdiff argument sfc_flux units are CU kg m-2 T-1 + scale_factor = US%Z_to_m + CFC11 => CS%CFC11 ; CFC12 => CS%CFC12 ! Use a tridiagonal solver to determine the concentrations after the @@ -340,8 +342,8 @@ subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, C 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) + if (CS%id_cfc11_cmor > 0) call post_data(CS%id_cfc11_cmor, (GV%Rho0*US%R_to_kg_m3)*CFC11, CS%diag) + if (CS%id_cfc12_cmor > 0) call post_data(CS%id_cfc12_cmor, (GV%Rho0*US%R_to_kg_m3)*CFC12, CS%diag) end subroutine CFC_cap_column_physics @@ -421,14 +423,15 @@ 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) +subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id_cfc12_atm) type(ocean_grid_type), intent(in ) :: G !< The ocean's grid structure. + type(unit_scale_type), intent(in ) :: US !< A dimensional unit scaling type 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] + real, intent(in ) :: Rho0 !< The mean ocean density [R ~> 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. @@ -436,11 +439,8 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc ! 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]. + kw_wo_sc_no_term, & ! gas transfer velocity, without the Schmidt number term [Z T-1 ~> m s-1]. + kw, & ! gas transfer velocity [Z T-1 ~> 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] @@ -450,9 +450,10 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc 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 :: kw_coeff ! A coefficient used to compute the piston velocity [Z T-1 T2 L-2 = Z T L-2 ~> s / m] + real :: Rho0_kg_m3 ! Rho0 in non-scaled units [kg m-3] real, parameter :: pa_to_atm = 9.8692316931427e-6 ! factor for converting from Pa to atm. + real :: press_to_atm ! converts from model pressure units to atm integer :: i, j, m, is, ie, js, je is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec @@ -479,6 +480,17 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc "has not been implemented yet.") endif + !--------------------------------------------------------------------- + ! Gas exchange/piston velocity parameter + !--------------------------------------------------------------------- + ! From a = 0.251 cm/hr s^2/m^2 in Wannikhof 2014 + ! = 6.97e-7 m/s s^2/m^2 [Z T-1 T2 L-2 = Z T L-2 ~> s / m] + kw_coeff = (US%m_to_Z*US%s_to_T*US%L_to_m**2) * 6.97e-7 + + ! set unit conversion factors + Rho0_kg_m3 = Rho0 * US%R_to_kg_m3 + press_to_atm = US%R_to_kg_m3*US%L_T_to_m_s**2 * pa_to_atm + do j=js,je ; do i=is,ie ! ta in hectoKelvin ta = max(0.01, (sfc_state%SST(i,j) + 273.15) * 0.01) @@ -490,84 +502,59 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, Rho0, Time, id_cfc11_atm, id_cfc ! 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)) + 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 + ! CFC flux units: CU Z T-1 kg m-3 = mol kg-1 Z T-1 kg m-3 ~> mol m-2 s-1 + kw(i,j) = kw_wo_sc_no_term(i,j) * sqrt(660.0 / sc_11) + cair(i,j) = press_to_atm * alpha_11 * cfc11_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc11_flux(i,j) = kw(i,j) * (cair(i,j) - sfc_state%sfc_CFC11(i,j)) * Rho0_kg_m3 + + kw(i,j) = kw_wo_sc_no_term(i,j) * sqrt(660.0 / sc_12) + cair(i,j) = press_to_atm * alpha_12 * cfc12_atm(i,j) * fluxes%p_surf_full(i,j) + fluxes%cfc12_flux(i,j) = kw(i,j) * (cair(i,j) - sfc_state%sfc_CFC12(i,j)) * Rho0_kg_m3 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 + 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) + + ! Coefficients for calculating CFC11 solubilities ! 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) + + real, parameter :: d1_11 = -232.0411 ! [nondim] + real, parameter :: d2_11 = 322.5546 ! [hectoKelvin-1] + real, parameter :: d3_11 = 120.4956 ! [log(hectoKelvin)-1] + real, parameter :: d4_11 = -1.39165 ! [hectoKelvin-2] + + real, parameter :: e1_11 = -0.146531 ! [PSU-1] + real, parameter :: e2_11 = 0.093621 ! [PSU-1 hectoKelvin-1] + real, parameter :: e3_11 = -0.0160693 ! [PSU-2 hectoKelvin-2] + + ! Coefficients for calculating CFC12 solubilities + ! from Table 5 in Warner and Weiss (1985) DSR, vol 32. + + real, parameter :: d1_12 = -220.2120 ! [nondim] + real, parameter :: d2_12 = 301.8695 ! [hectoKelvin-1] + real, parameter :: d3_12 = 114.8533 ! [log(hectoKelvin)-1] + real, parameter :: d4_12 = -1.39165 ! [hectoKelvin-2] + + real, parameter :: e1_12 = -0.147718 ! [PSU-1] + real, parameter :: e2_12 = 0.093175 ! [PSU-1 hectoKelvin-1] + real, parameter :: e3_12 = -0.0157340 ! [PSU-2 hectoKelvin-2] + + real :: factor ! introduce units to result [mol kg-1 atm-1] + + ! Eq. 9 from Warner and Weiss (1985) DSR, vol 32. 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)) * &