Skip to content

+Rescaled flux arguments to tracer_vertdiff() #17

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
marshallward merged 2 commits into from
Dec 5, 2021
Merged

+Rescaled flux arguments to tracer_vertdiff() #17

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
13a1cde
+Rescaled flux arguments to tracer_vertdiff()
Hallberg-NOAA Nov 29, 2021
2dbd055
Merge branch 'dev/gfdl' into tracer_sfc_flux_rescale
marshallward Dec 4, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 12 additions & 13 deletions src/core/MOM_forcing_type.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,8 +189,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 [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].
cfc11_flux => NULL(), & !< flux of cfc_11 into the ocean [CU R Z T-1 kg m-3 ~> mol m-2 s-1]
cfc12_flux => NULL(), & !< flux of cfc_12 into the ocean [CU R Z T-1 kg 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]

Expand Down Expand Up @@ -1089,20 +1089,19 @@ subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift)
call hchksum(fluxes%seaice_melt_heat, mesg//" fluxes%seaice_melt_heat", G%HI, &
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)
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%L_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)
call hchksum(fluxes%cfc11_flux, mesg//" fluxes%cfc11_flux", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s)
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)
call hchksum(fluxes%cfc12_flux, mesg//" fluxes%cfc12_flux", G%HI, haloshift=hshift, scale=US%RZ_T_to_kg_m2s)
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)) &
call hchksum(fluxes%TKE_tidal, mesg//" fluxes%TKE_tidal", G%HI, haloshift=hshift, &
scale=US%RZ3_T3_to_W_m2)
call hchksum(fluxes%TKE_tidal, mesg//" fluxes%TKE_tidal", G%HI, haloshift=hshift, scale=US%RZ3_T3_to_W_m2)
if (associated(fluxes%ustar_tidal)) &
call hchksum(fluxes%ustar_tidal, mesg//" fluxes%ustar_tidal", G%HI, haloshift=hshift, scale=US%Z_to_m*US%s_to_T)
if (associated(fluxes%lrunoff)) &
Expand Down Expand Up @@ -1301,22 +1300,22 @@ 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
! units for cfc11_flux and cfc12_flux are [Conc R Z T-1 ~> 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%Z_to_m*US%s_to_T,&
'Gas exchange flux of CFC11 into the ocean ', &
'mol m-2 s-1', conversion=US%RZ_T_to_kg_m2s, &
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%Z_to_m*US%s_to_T,&
'Gas exchange flux of CFC12 into the ocean ', &
'mol m-2 s-1', conversion=US%RZ_T_to_kg_m2s, &
cmor_field_name='fgcfc12', &
cmor_long_name='Surface Downward CFC12 Flux', &
cmor_standard_name='surface_downward_cfc12_flux')
Expand Down
63 changes: 26 additions & 37 deletions src/tracer/MOM_CFC_cap.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -304,46 +304,37 @@ subroutine CFC_cap_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, C
! 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 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

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
! 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, &
call applyTracerBoundaryFluxesInOut(G, GV, CS%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)
call tracer_vertdiff(h_work, ea, eb, dt, CS%CFC11, G, GV, sfc_flux=fluxes%cfc11_flux)

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, &
call applyTracerBoundaryFluxesInOut(G, GV, CS%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)
call tracer_vertdiff(h_work, ea, eb, dt, CS%CFC12, G, GV, sfc_flux=fluxes%cfc12_flux)
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)
call tracer_vertdiff(h_old, ea, eb, dt, CS%CFC11, G, GV, sfc_flux=fluxes%cfc11_flux)
call tracer_vertdiff(h_old, ea, eb, dt, CS%CFC12, G, GV, sfc_flux=fluxes%cfc12_flux)
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, (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)
if (CS%id_cfc11_cmor > 0) call post_data(CS%id_cfc11_cmor, (GV%Rho0*US%R_to_kg_m3)*CS%CFC11, CS%diag)
if (CS%id_cfc12_cmor > 0) call post_data(CS%id_cfc12_cmor, (GV%Rho0*US%R_to_kg_m3)*CS%CFC12, CS%diag)

end subroutine CFC_cap_column_physics

Expand Down Expand Up @@ -449,11 +440,10 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id
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_11, sc_12 ! The Schmidt numbers of CFC 11 and CFC 12 [nondim].
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
real, parameter :: pa_to_atm = 9.8692316931427e-6 ! factor for converting from Pa to atm [atm Pa-1].
real :: press_to_atm ! converts from model pressure units to atm [atm T2 R-1 L-2 ~> atm Pa-1]
integer :: i, j, m, is, ie, js, je

is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
Expand Down Expand Up @@ -488,7 +478,6 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id
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
Expand All @@ -506,14 +495,14 @@ subroutine CFC_cap_fluxes(fluxes, sfc_state, G, US, Rho0, Time, id_cfc11_atm, id
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: CU Z T-1 kg m-3 = mol kg-1 Z T-1 kg m-3 ~> mol m-2 s-1
! CFC flux units: CU R Z T-1 = mol kg-1 R Z T-1 ~> 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
fluxes%cfc11_flux(i,j) = kw(i,j) * (cair(i,j) - sfc_state%sfc_CFC11(i,j)) * Rho0

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
fluxes%cfc12_flux(i,j) = kw(i,j) * (cair(i,j) - sfc_state%sfc_CFC12(i,j)) * Rho0
enddo ; enddo

end subroutine CFC_cap_fluxes
Expand All @@ -524,7 +513,7 @@ subroutine get_solubility(alpha_11, alpha_12, ta, sal , mask)
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(in ) :: mask !< ocean mask [nondim]

! Local variables

Expand Down Expand Up @@ -574,17 +563,17 @@ subroutine comp_CFC_schmidt(sst_in, cfc11_sc, cfc12_sc)
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
real , parameter :: a_11 = 3579.2 ! CFC11 Schmidt number fit coefficient [nondim]
real , parameter :: b_11 = -222.63 ! CFC11 Schmidt number fit coefficient [degC-1]
real , parameter :: c_11 = 7.5749 ! CFC11 Schmidt number fit coefficient [degC-2]
real , parameter :: d_11 = -0.14595 ! CFC11 Schmidt number fit coefficient [degC-3]
real , parameter :: e_11 = 0.0011874 ! CFC11 Schmidt number fit coefficient [degC-4]
real , parameter :: a_12 = 3828.1 ! CFC12 Schmidt number fit coefficient [nondim]
real , parameter :: b_12 = -249.86 ! CFC12 Schmidt number fit coefficient [degC-1]
real , parameter :: c_12 = 8.7603 ! CFC12 Schmidt number fit coefficient [degC-2]
real , parameter :: d_12 = -0.1716 ! CFC12 Schmidt number fit coefficient [degC-3]
real , parameter :: e_12 = 0.001408 ! CFC12 Schmidt number fit coefficient [degC-4]
real :: sst ! A range-limited sea surface temperature [degC]


! clip SST to avoid bad values
Expand Down
42 changes: 15 additions & 27 deletions src/tracer/MOM_OCMIP2_CFC.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,9 +31,6 @@ module MOM_OCMIP2_CFC
public OCMIP2_CFC_column_physics, OCMIP2_CFC_surface_state
public OCMIP2_CFC_stock, OCMIP2_CFC_end


integer, parameter :: NTR = 2 !< the number of tracers in this module.

!> The control structure for the OCMPI2_CFC tracer package
type, public :: OCMIP2_CFC_CS ; private
character(len=200) :: IC_file !< The file in which the CFC initial values can
Expand Down Expand Up @@ -96,18 +93,16 @@ function register_OCMIP2_CFC(HI, GV, param_file, CS, tr_Reg, restart_CS)
type(tracer_registry_type), &
pointer :: tr_Reg !< A pointer to the tracer registry.
type(MOM_restart_CS), target, intent(inout) :: restart_CS !< MOM restart control struct
! This subroutine is used to register tracer fields and subroutines
! to be used with MOM.

! Local variables
character(len=40) :: mdl = "MOM_OCMIP2_CFC" ! 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"
# 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.
real :: d11_dflt(4), d12_dflt(4) ! in the expressions for the solubility and
real :: e11_dflt(3), e12_dflt(3) ! Schmidt numbers [various units by element].
character(len=48) :: flux_units ! The units for tracer fluxes.
logical :: register_OCMIP2_CFC
integer :: isd, ied, jsd, jed, nz, m
Expand Down Expand Up @@ -330,10 +325,6 @@ subroutine initialize_OCMIP2_CFC(restart, day, G, GV, US, h, diag, OBC, CS, &
type(sponge_CS), pointer :: sponge_CSp !< A pointer to the control structure for
!! the sponges, if they are in use.
!! Otherwise this may be unassociated.
! This subroutine initializes the NTR tracer fields in tr(:,:,:,:)
! and it sets up the tracer output.

logical :: from_file = .false.

if (.not.associated(CS)) return

Expand Down Expand Up @@ -441,9 +432,8 @@ subroutine OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US

! Local variables
real, dimension(SZI_(G),SZJ_(G)) :: &
CFC11_flux, & ! The fluxes of CFC11 and CFC12 into the ocean, in the
CFC12_flux ! units of CFC concentrations times meters per second.
real, pointer, dimension(:,:,:) :: CFC11 => NULL(), CFC12 => NULL()
CFC11_flux, & ! The fluxes of CFC11 and CFC12 into the ocean, in unscaled units of
CFC12_flux ! CFC concentrations times meters per second [CU R Z T-1 ~> CU kg m-2 s-1]
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work ! Used so that h can be modified [H ~> m or kg m-2]
integer :: i, j, k, m, is, ie, js, je, nz, idim(4), jdim(4)

Expand All @@ -452,35 +442,33 @@ subroutine OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US

if (.not.associated(CS)) return

CFC11 => CS%CFC11 ; CFC12 => CS%CFC12

! These two calls unpack the fluxes from the input arrays.
! The -GV%Rho0 changes the sign convention of the flux and changes the units
! of the flux from [Conc. m s-1] to [Conc. kg m-2 T-1].
! The -GV%Rho0 changes the sign convention of the flux and with the scaling factors changes
! the units of the flux from [Conc. m s-1] to [Conc. R Z T-1 ~> Conc. kg m-2 s-1].
call extract_coupler_type_data(fluxes%tr_fluxes, CS%ind_cfc_11_flux, CFC11_flux, &
scale_factor=-GV%Rho0*US%R_to_kg_m3*US%T_to_s, idim=idim, jdim=jdim)
scale_factor=-GV%Rho0*US%m_to_Z*US%T_to_s, idim=idim, jdim=jdim)
call extract_coupler_type_data(fluxes%tr_fluxes, CS%ind_cfc_12_flux, CFC12_flux, &
scale_factor=-GV%Rho0*US%R_to_kg_m3*US%T_to_s, idim=idim, jdim=jdim)
scale_factor=-GV%Rho0*US%m_to_Z*US%T_to_s, idim=idim, jdim=jdim)

! 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, &
call applyTracerBoundaryFluxesInOut(G, GV, CS%CFC11, dt, fluxes, h_work, &
evap_CFL_limit, minimum_forcing_depth)
call tracer_vertdiff(h_work, ea, eb, dt, CFC11, G, GV, sfc_flux=CFC11_flux)
call tracer_vertdiff(h_work, ea, eb, dt, CS%CFC11, G, GV, sfc_flux=CFC11_flux)

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, &
call applyTracerBoundaryFluxesInOut(G, GV, CS%CFC12, dt, fluxes, h_work, &
evap_CFL_limit, minimum_forcing_depth)
call tracer_vertdiff(h_work, ea, eb, dt, CFC12, G, GV, sfc_flux=CFC12_flux)
call tracer_vertdiff(h_work, ea, eb, dt, CS%CFC12, G, GV, sfc_flux=CFC12_flux)
else
call tracer_vertdiff(h_old, ea, eb, dt, CFC11, G, GV, sfc_flux=CFC11_flux)
call tracer_vertdiff(h_old, ea, eb, dt, CFC12, G, GV, sfc_flux=CFC12_flux)
call tracer_vertdiff(h_old, ea, eb, dt, CS%CFC11, G, GV, sfc_flux=CFC11_flux)
call tracer_vertdiff(h_old, ea, eb, dt, CS%CFC12, G, GV, sfc_flux=CFC12_flux)
endif

! Write out any desired diagnostics from tracer sources & sinks here.
Expand Down
Loading