Skip to content

Commit

Permalink
Merge pull request #17 from Hallberg-NOAA/tracer_sfc_flux_rescale
Browse files Browse the repository at this point in the history 
+Rescaled flux arguments to tracer_vertdiff()
  • Loading branch information
@marshallward
marshallward authored Dec 5, 2021
2 parents 585cc70 + 2dbd055 commit a44e460
Show file tree
Hide file tree
Showing 4 changed files with 67 additions and 93 deletions.
25 changes: 12 additions & 13 deletions src/core/MOM_forcing_type.F90
View file
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
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
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

0 comments on commit a44e460

Please sign in to comment.