diff --git a/src/core/MOM.F90 b/src/core/MOM.F90
index d0387f5011..985b8ee96a 100644
--- a/src/core/MOM.F90
+++ b/src/core/MOM.F90
@@ -2073,7 +2073,7 @@ subroutine initialize_MOM(Time, param_file, dirs, CS, Time_in)
(LEN_TRIM(dirs%input_filename) == 1))
call tracer_flow_control_init(.not.new_sim, Time, G, GV, CS%h, param_file, &
CS%diag, CS%OBC, CS%tracer_flow_CSp, CS%sponge_CSp, &
- CS%ALE_sponge_CSp, CS%diag_to_Z_CSp)
+ CS%ALE_sponge_CSp, CS%diag_to_Z_CSp, CS%tv)
call cpu_clock_begin(id_clock_pass_init)
call do_group_pass(CS%pass_uv_T_S_h, G%Domain)
diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90
index 390d2379f4..651eb3ee48 100644
--- a/src/core/MOM_forcing_type.F90
+++ b/src/core/MOM_forcing_type.F90
@@ -80,7 +80,10 @@ module MOM_forcing_type
vprec => NULL(), & !< virtual liquid precip associated w/ SSS restoring ( kg/(m^2 s) )
lrunoff => NULL(), & !< liquid river runoff entering ocean ( kg/(m^2 s) )
frunoff => NULL(), & !< frozen river runoff (calving) entering ocean ( kg/(m^2 s) )
- seaice_melt => NULL() !< seaice melt (positive) or formation (negative) ( kg/(m^2 s) )
+ seaice_melt => NULL(), & !< seaice melt (positive) or formation (negative) ( kg/(m^2 s) )
+ netMassIn => NULL(), & !< Sum of water mass flux out of the ocean ( kg/(m^2 s) )
+ netMassOut => NULL(), & !< Net water mass flux into of the ocean ( kg/(m^2 s) )
+ netSalt => NULL() !< Net salt entering the ocean
! heat associated with water crossing ocean surface
real, pointer, dimension(:,:) :: &
@@ -410,7 +413,6 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
netMassInOut(i) = GV%kg_m2_to_H * netMassInOut(i)
netMassOut(i) = GV%kg_m2_to_H * netMassOut(i)
-
! surface heat fluxes from radiation and turbulent fluxes (K * H)
! (H=m for Bouss, H=kg/m2 for non-Bouss)
net_heat(i) = scale * dt * J_m2_to_H * &
@@ -488,9 +490,10 @@ subroutine extractFluxes1d(G, GV, fluxes, optics, nsw, j, dt,
! Convert salt_flux from kg (salt)/(m^2 * s) to
! Boussinesq: (ppt * m)
! non-Bouss: (g/m^2)
- if (ASSOCIATED(fluxes%salt_flux)) &
+ if (ASSOCIATED(fluxes%salt_flux)) then
Net_salt(i) = (scale * dt * (1000.0 * fluxes%salt_flux(i,j))) * GV%kg_m2_to_H
-
+ fluxes%netSalt(i,j) = Net_salt(i)
+ endif
! Diagnostics follow...
! Initialize heat_content_massin that is diagnosed in mixedlayer_convection or
@@ -2186,6 +2189,9 @@ subroutine allocate_forcing_type(G, fluxes, stress, ustar, water, heat, shelf, p
call myAlloc(fluxes%lrunoff,isd,ied,jsd,jed, water)
call myAlloc(fluxes%frunoff,isd,ied,jsd,jed, water)
call myAlloc(fluxes%seaice_melt,isd,ied,jsd,jed, water)
+ call myAlloc(fluxes%netMassOut,isd,ied,jsd,jed, water)
+ call myAlloc(fluxes%netMassIn,isd,ied,jsd,jed, water)
+ call myAlloc(fluxes%netSalt,isd,ied,jsd,jed, water)
call myAlloc(fluxes%sw,isd,ied,jsd,jed, heat)
call myAlloc(fluxes%lw,isd,ied,jsd,jed, heat)
diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90
index 587e55f6ce..906c59cd9b 100644
--- a/src/parameterizations/vertical/MOM_diabatic_aux.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90
@@ -100,9 +100,9 @@ module MOM_diabatic_aux
!! at the river mouths to "rivermix_depth" meters
real :: rivermix_depth = 0.0 !< The depth to which rivers are mixed if
!! do_rivermix = T, in m.
- real :: minimum_forcing_depth = 0.001 !< The smallest depth over which forcing is
+ real, public :: minimum_forcing_depth = 0.001 !< The smallest depth over which forcing is
!! applied, in m.
- real :: evap_CFL_limit = 0.8 !< The largest fraction of a layer that can be
+ real, public :: evap_CFL_limit = 0.8 !< The largest fraction of a layer that can be
!! evaporated in one time-step (non-dim).
logical :: reclaim_frazil !< If true, try to use any frazil heat deficit to
@@ -574,7 +574,6 @@ subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S)
real :: c1(SZIB_(G),SZK_(G)) ! tridiagonal solver.
real :: h_tr, b_denom_1
integer :: i, j, k
-
!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,hold,eb,T,S,ea) &
!$OMP private(h_tr,b1,d1,c1,b_denom_1)
do j=js,je
@@ -803,15 +802,14 @@ end subroutine diagnoseMLDbyDensityDifference
!> Update the thickness, temperature, and salinity due to thermodynamic
!! boundary forcing (contained in fluxes type) applied to h, tv%T and tv%S,
!! and calculate the TKE implications of this heating.
-subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, ea, h, tv, &
+subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, &
aggregate_FW_forcing, cTKE, dSV_dT, dSV_dS)
type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux
type(ocean_grid_type), intent(in) :: G !< Grid structure
- type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
+ type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
real, intent(in) :: dt !< Time-step over which forcing is applied (s)
type(forcing), intent(inout) :: fluxes !< Surface fluxes container
type(optics_type), pointer :: optics !< Optical properties container
- real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: ea !< The entrainment distance at interfaces (H units)
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness in H units
type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics container
!> If False, treat in/out fluxes separately.
@@ -945,13 +943,15 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, ea, h, tv, &
! ea is for passive tracers
do i=is,ie
- ea(i,j,1) = netMassInOut(i)
+ ! ea(i,j,1) = netMassInOut(i)
if (aggregate_FW_forcing) then
netMassOut(i) = netMassInOut(i)
netMassIn(i) = 0.
else
netMassIn(i) = netMassInOut(i) - netMassOut(i)
endif
+ fluxes%netMassOut(i,j) = netMassOut(i)
+ fluxes%netMassIn(i,j) = netMassIn(i)
enddo
! Apply the surface boundary fluxes in three steps:
@@ -973,7 +973,6 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, ea, h, tv, &
! Update the forcing by the part to be consumed within the present k-layer.
! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
netMassIn(i) = netMassIn(i) - dThickness
-
! This line accounts for the temperature of the mass exchange
Temp_in = T2d(i,k)
Salin_in = 0.0
diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90
index b6024426dc..092e0cc37b 100644
--- a/src/parameterizations/vertical/MOM_diabatic_driver.F90
+++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90
@@ -55,6 +55,7 @@ module MOM_diabatic_driver
use MOM_ALE_sponge, only : apply_ALE_sponge, ALE_sponge_CS
use MOM_time_manager, only : operator(<=), time_type ! for testing itides (BDM)
use MOM_tracer_flow_control, only : call_tracer_column_fns, tracer_flow_control_CS
+use MOM_tracer_diabatic, only : tracer_vertdiff
use MOM_variables, only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs
use MOM_variables, only : cont_diag_ptrs, MOM_thermovar_chksum, p3d
use MOM_verticalGrid, only : verticalGrid_type
@@ -141,6 +142,8 @@ module MOM_diabatic_driver
!! is statically unstable.
logical :: debug !< If true, write verbose checksums for debugging purposes.
logical :: debugConservation !< If true, monitor conservation and extrema.
+ logical :: tracer_tridiag !< If true, use tracer_vertdiff instead of tridiagTS for
+ !< vertical diffusion of T and S
logical :: debug_energy_req ! If true, test the mixing energy requirement code.
type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output
real :: MLDdensityDifference !< Density difference used to determine MLD_user
@@ -244,6 +247,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
! one time step (m for Bouss, kg/m^2 for non-Bouss)
Kd, & ! diapycnal diffusivity of layers (m^2/sec)
h_orig, & ! initial layer thicknesses (m for Bouss, kg/m^2 for non-Bouss)
+ h_prebound, & ! initial layer thicknesses (m for Bouss, kg/m^2 for non-Bouss)
hold, & ! layer thickness before diapycnal entrainment, and later
! the initial layer thicknesses (if a mixed layer is used),
! (m for Bouss, kg/m^2 for non-Bouss)
@@ -411,6 +415,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
if (CS%debugConservation) call MOM_state_stats('geothermal', u, v, h, tv%T, tv%S, G)
endif
+
! Whenever thickness changes let the diag manager know, target grids
! for vertical remapping may need to be regenerated.
call diag_update_target_grids(CS%diag)
@@ -719,7 +724,7 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
if (CS%useALEalgorithm) then
call cpu_clock_begin(id_clock_remap)
- ! Changes made to following fields: ea(:,:,1), h, tv%T and tv%S.
+ ! Changes made to following fields: h, tv%T and tv%S.
! save prior values for diagnostics
if(CS%boundary_forcing_tendency_diag) then
@@ -730,10 +735,14 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
enddo ; enddo ; enddo
endif
+ do k=1,nz ; do j=js,je ; do i=is,ie
+ h_prebound(i,j,k) = h(i,j,k)
+ enddo ; enddo ; enddo
+!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_orig,h,eaml,ebml)
if (CS%use_energetic_PBL) then
-
call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, &
- ea, h, tv, CS%aggregate_FW_forcing, cTKE, dSV_dT, dSV_dS)
+ h, tv, CS%aggregate_FW_forcing, cTKE, dSV_dT, dSV_dS)
+
if (CS%debug) then
call hchksum(ea, "after applyBoundaryFluxes ea",G%HI,haloshift=0)
@@ -782,9 +791,8 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
endif
else
-
call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, &
- ea, h, tv, CS%aggregate_FW_forcing)
+ h, tv, CS%aggregate_FW_forcing)
endif ! endif for CS%use_energetic_PBL
@@ -804,7 +812,6 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
endif ! endif for (CS%useALEalgorithm)
-
! Update h according to divergence of the difference between
! ea and eb. We keep a record of the original h in hold.
! In the following, the checks for negative values are to guard
@@ -926,7 +933,12 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
! This simpler form allows T & S to be too dense for the layers
! between the buffer layers and the interior.
! Changes: T, S
- call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
+ if (CS%tracer_tridiag) then
+ call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV)
+ call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV)
+ else
+ call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
+ endif
endif ! massless_match_targets
call cpu_clock_end(id_clock_tridiag)
@@ -1005,7 +1017,12 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
endif
! Changes T and S via the tridiagonal solver; no change to h
- call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
+ if(CS%tracer_tridiag) then
+ call tracer_vertdiff(hold, ea, eb, dt, tv%T, G, GV)
+ call tracer_vertdiff(hold, ea, eb, dt, tv%S, G, GV)
+ else
+ call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
+ endif
! diagnose temperature, salinity, heat, and salt tendencies
if(CS%diabatic_diff_tendency_diag) then
@@ -1071,7 +1088,6 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
enddo ; enddo ; enddo
endif
-
! mixing of passive tracers from massless boundary layers to interior
call cpu_clock_begin(id_clock_tracers)
if (CS%mix_boundary_tracers) then
@@ -1122,8 +1138,19 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
enddo ; enddo
do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ; enddo
enddo
- call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, dt, G, GV, tv, &
- CS%optics, CS%tracer_flow_CSp)
+
+
+ if (CS%useALEalgorithm) then
+ ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
+ ! so hold should be h_orig
+ call call_tracer_column_fns(h_prebound, h, ea, eb, fluxes, dt, G, GV, tv, &
+ CS%optics, CS%tracer_flow_CSp, CS%debug, &
+ evap_CFL_limit = CS%diabatic_aux_CSp%evap_CFL_limit, &
+ minimum_forcing_depth = CS%diabatic_aux_CSp%minimum_forcing_depth)
+ else
+ call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, dt, G, GV, tv, &
+ CS%optics, CS%tracer_flow_CSp, CS%debug)
+ endif
elseif (associated(visc%Kd_extra_S)) then ! extra diffusivity for passive tracers
@@ -1144,12 +1171,28 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
eatr(i,j,k) = ea(i,j,k) + add_ent
enddo ; enddo ; enddo
- call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, dt, G, GV, tv, &
- CS%optics, CS%tracer_flow_CSp)
+ if (CS%useALEalgorithm) then
+ ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
+ call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, dt, G, GV, tv, &
+ CS%optics, CS%tracer_flow_CSp, CS%debug,&
+ evap_CFL_limit = CS%diabatic_aux_CSp%evap_CFL_limit, &
+ minimum_forcing_depth = CS%diabatic_aux_CSp%minimum_forcing_depth)
+ else
+ call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, dt, G, GV, tv, &
+ CS%optics, CS%tracer_flow_CSp, CS%debug)
+ endif
else
- call call_tracer_column_fns(hold, h, ea, eb, fluxes, dt, G, GV, tv, &
- CS%optics, CS%tracer_flow_CSp)
+ if (CS%useALEalgorithm) then
+ ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
+ call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, dt, G, GV, tv, &
+ CS%optics, CS%tracer_flow_CSp, CS%debug, &
+ evap_CFL_limit = CS%diabatic_aux_CSp%evap_CFL_limit, &
+ minimum_forcing_depth = CS%diabatic_aux_CSp%minimum_forcing_depth)
+ else
+ call call_tracer_column_fns(hold, h, ea, eb, fluxes, dt, G, GV, tv, &
+ CS%optics, CS%tracer_flow_CSp, CS%debug)
+ endif
endif ! (CS%mix_boundary_tracers)
@@ -1436,7 +1479,7 @@ subroutine adiabatic(h, tv, fluxes, dt, G, GV, CS)
zeros(:,:,:) = 0.0
call call_tracer_column_fns(h, h, zeros, zeros, fluxes, dt, G, GV, tv, &
- CS%optics, CS%tracer_flow_CSp)
+ CS%optics, CS%tracer_flow_CSp, CS%debug)
end subroutine adiabatic
@@ -1870,6 +1913,11 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag,
"entrainment at the bottom is at least sqrt(Kd_BBL_tr*dt) \n"//&
"over the same distance.", units="m2 s-1", default=0.)
endif
+
+ call get_param(param_file, mod, "TRACER_TRIDIAG", CS%tracer_tridiag, &
+ "If true, use the passive tracer tridiagonal solver for T and S\n", &
+ default=.false.)
+
! Register all available diagnostics for this module.
if (GV%Boussinesq) then ; thickness_units = "meter"
diff --git a/src/tracer/DOME_tracer.F90 b/src/tracer/DOME_tracer.F90
index 993b67eea9..4af723e5cb 100644
--- a/src/tracer/DOME_tracer.F90
+++ b/src/tracer/DOME_tracer.F90
@@ -68,7 +68,7 @@ module DOME_tracer
use MOM_time_manager, only : time_type, get_time
use MOM_tracer_registry, only : register_tracer, tracer_registry_type
use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values
-use MOM_tracer_registry, only : tracer_vertdiff
+use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
use MOM_variables, only : surface
use MOM_verticalGrid, only : verticalGrid_type
@@ -415,13 +415,16 @@ subroutine initialize_DOME_tracer(restart, day, G, GV, h, diag, OBC, CS, &
end subroutine initialize_DOME_tracer
-subroutine DOME_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS)
+subroutine DOME_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, &
+ evap_CFL_limit, minimum_forcing_depth)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old, h_new, ea, eb
type(forcing), intent(in) :: fluxes
real, intent(in) :: dt
type(DOME_tracer_CS), pointer :: CS
+ real, optional,intent(in) :: evap_CFL_limit
+ real, optional,intent(in) :: minimum_forcing_depth
! This subroutine applies diapycnal diffusion and any other column
! tracer physics or chemistry to the tracers from this file.
! This is a simple example of a set of advected passive tracers.
@@ -447,14 +450,26 @@ subroutine DOME_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV,
real :: b1(SZI_(G)) ! b1 and c1 are variables used by the
real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver.
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
integer :: i, j, k, is, ie, js, je, nz, m
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
if (.not.associated(CS)) return
- do m=1,NTR
- call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
- enddo
+ if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
+ do m=1,NTR
+ 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, CS%tr(:,:,:,m) , dt, fluxes, h_work, &
+ evap_CFL_limit, minimum_forcing_depth)
+ call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ else
+ do m=1,NTR
+ call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ endif
if (CS%mask_tracers) then
do m = 1,NTR ; if (CS%id_tracer(m) > 0) then
diff --git a/src/tracer/ISOMIP_tracer.F90 b/src/tracer/ISOMIP_tracer.F90
index 6b826f40d5..9715b8a3dc 100644
--- a/src/tracer/ISOMIP_tracer.F90
+++ b/src/tracer/ISOMIP_tracer.F90
@@ -29,7 +29,7 @@ module ISOMIP_tracer
use MOM_time_manager, only : time_type, get_time
use MOM_tracer_registry, only : register_tracer, tracer_registry_type
use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values
-use MOM_tracer_registry, only : tracer_vertdiff
+use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
use MOM_variables, only : surface
use MOM_open_boundary, only : ocean_OBC_type
use MOM_verticalGrid, only : verticalGrid_type
@@ -302,13 +302,16 @@ end subroutine initialize_ISOMIP_tracer
!> This subroutine applies diapycnal diffusion and any other column
! tracer physics or chemistry to the tracers from this file.
! This is a simple example of a set of advected passive tracers.
-subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS)
+subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, &
+ evap_CFL_limit, minimum_forcing_depth)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old, h_new, ea, eb
type(forcing), intent(in) :: fluxes
real, intent(in) :: dt
type(ISOMIP_tracer_CS), pointer :: CS
+ real, optional,intent(in) :: evap_CFL_limit
+ real, optional,intent(in) :: minimum_forcing_depth
! Arguments: h_old - Layer thickness before entrainment, in m or kg m-2.
! (in) h_new - Layer thickness after entrainment, in m or kg m-2.
@@ -331,6 +334,7 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G
real :: b1(SZI_(G)) ! b1 and c1 are variables used by the
real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver.
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
integer :: i, j, k, is, ie, js, je, nz, m
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
@@ -349,9 +353,20 @@ subroutine ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, G
enddo ; enddo
enddo
- do m=1,NTR
- call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
- enddo
+ if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
+ do m=1,NTR
+ 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, CS%tr(:,:,:,m) , dt, fluxes, h_work, &
+ evap_CFL_limit, minimum_forcing_depth)
+ call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ else
+ do m=1,NTR
+ call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ endif
if (CS%mask_tracers) then
do m = 1,NTR ; if (CS%id_tracer(m) > 0) then
diff --git a/src/tracer/MOM_OCMIP2_CFC.F90 b/src/tracer/MOM_OCMIP2_CFC.F90
index 8ad94af4c0..cba46ef234 100644
--- a/src/tracer/MOM_OCMIP2_CFC.F90
+++ b/src/tracer/MOM_OCMIP2_CFC.F90
@@ -79,7 +79,7 @@ module MOM_OCMIP2_CFC
use MOM_time_manager, only : time_type, get_time
use MOM_tracer_registry, only : register_tracer, tracer_registry_type
use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values
-use MOM_tracer_registry, only : tracer_vertdiff
+use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
use MOM_tracer_Z_init, only : tracer_Z_init
use MOM_variables, only : surface
use MOM_verticalGrid, only : verticalGrid_type
@@ -525,13 +525,16 @@ subroutine init_tracer_CFC(h, tr, name, land_val, IC_val, G, CS)
end subroutine init_tracer_CFC
-subroutine OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS)
+subroutine OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, &
+ evap_CFL_limit, minimum_forcing_depth)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old, h_new, ea, eb
type(forcing), intent(in) :: fluxes
real, intent(in) :: dt
type(OCMIP2_CFC_CS), pointer :: CS
+ real, optional,intent(in) :: evap_CFL_limit
+ real, optional,intent(in) :: minimum_forcing_depth
! This subroutine applies diapycnal diffusion and any other column
! tracer physics or chemistry to the tracers from this file.
! CFCs are relatively simple, as they are passive tracers. with only a surface
@@ -562,6 +565,7 @@ subroutine OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS
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, CFC12
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
integer :: i, j, k, is, ie, js, je, nz, m
is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
@@ -580,8 +584,24 @@ subroutine OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS
! Use a tridiagonal solver to determine the concentrations after the
! surface source is applied and diapycnal advection and diffusion occurs.
- 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)
+ 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=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, &
+ evap_CFL_limit, minimum_forcing_depth)
+ call tracer_vertdiff(h_work, ea, eb, dt, 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)
+ endif
! Write out any desired diagnostics.
if (CS%mask_tracers) then
diff --git a/src/tracer/MOM_generic_tracer.F90 b/src/tracer/MOM_generic_tracer.F90
index d45cd78ca5..b36e1bf59d 100644
--- a/src/tracer/MOM_generic_tracer.F90
+++ b/src/tracer/MOM_generic_tracer.F90
@@ -61,13 +61,12 @@ module MOM_generic_tracer
use MOM_hor_index, only : hor_index_type
use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc
use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS
- use MOM_spatial_means, only : global_area_mean
use MOM_sponge, only : set_up_sponge_field, sponge_CS
use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS
use MOM_time_manager, only : time_type, get_time, set_time
+ use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
use MOM_tracer_registry, only : register_tracer, tracer_registry_type
use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values
- use MOM_tracer_registry, only : tracer_vertdiff
use MOM_tracer_Z_init, only : tracer_Z_init
use MOM_tracer_initialization_from_Z, only : MOM_initialize_tracer_from_Z
use MOM_variables, only : surface, thermo_var_ptrs
@@ -87,7 +86,7 @@ module MOM_generic_tracer
public MOM_generic_tracer_fluxes_accumulate
type, public :: MOM_generic_tracer_CS ; private
- character(len=200) :: IC_file ! The file in which the generic tracer initial values can
+ character(len = 200) :: IC_file ! The file in which the generic tracer initial values can
! be found, or an empty string for internal initialization.
logical :: Z_IC_file ! If true, the generic_tracer IC_file is in Z-space. The default is false.
real :: tracer_IC_val = 0.0 ! The initial value assigned to tracers.
@@ -125,25 +124,25 @@ module MOM_generic_tracer
! Register these tracers for restart
!
!
- ! call register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
+ ! call register_MOM_generic_tracer(G, param_file, CS, diag, tr_Reg, restart_CS)
!
!
function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
- type(hor_index_type), intent(in) :: HI
- type(verticalGrid_type), intent(in) :: GV
- type(param_file_type), intent(in) :: param_file
- type(MOM_generic_tracer_CS), pointer :: CS
- type(tracer_registry_type), pointer :: tr_Reg
- type(MOM_restart_CS), pointer :: restart_CS
+ type(hor_index_type), intent(in) :: HI
+ type(verticalGrid_type), intent(in) :: GV
+ type(param_file_type), intent(in) :: param_file
+ type(MOM_generic_tracer_CS), pointer :: CS
+ type(tracer_registry_type), pointer :: tr_Reg
+ type(MOM_restart_CS), pointer :: restart_CS
! This subroutine is used to register tracer fields and subroutines
! to be used with MOM.
- ! Arguments: HI - A horizontal index type structure.
- ! (in) GV - The ocean's vertical grid structure.
+ ! Arguments: G - The ocean's grid structure.
! (in) param_file - A structure indicating the open file to parse for
! model parameter values.
! (in/out) CS - A pointer that is set to point to the control structure
! for this module
+ ! (in) diag - A structure that is used to regulate diagnostic output.
! (in/out) tr_Reg - A pointer that is set to point to the control structure
! for the tracer advection and diffusion module.
! (in) restart_CS - A pointer to the restart control structure.
@@ -159,8 +158,8 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
character(len=fm_string_len) :: g_tracer_name,longname,units
real, dimension(:,:,:,:), pointer :: tr_field
real, dimension(:,:,:), pointer :: tr_ptr
- real, dimension(HI%isd:HI%ied, HI%jsd:HI%jed, GV%ke) :: grid_tmask
- integer, dimension(HI%isd:HI%ied, HI%jsd:HI%jed) :: grid_kmt
+ real, dimension(HI%isd:HI%ied, HI%jsd:HI%jed,GV%ke) :: grid_tmask
+ integer, dimension(HI%isd:HI%ied, HI%jsd:HI%jed) :: grid_kmt
type(vardesc) :: vdesc
register_MOM_generic_tracer = .false.
@@ -171,6 +170,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
endif
allocate(CS)
+
!Register all the generic tracers used and create the list of them.
!This can be called by ALL PE's. No array fields allocated.
if (.not. g_registered) then
@@ -178,6 +178,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
g_registered = .true.
endif
+
! Read all relevant parameters and write them to the model log.
call log_version(param_file, sub_name, version, "")
call get_param(param_file, sub_name, "GENERIC_TRACER_IC_FILE", CS%IC_file, &
@@ -204,7 +205,8 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
ntau=1 ! MOM needs the fields at only one time step
- ! At this point HI%mask2dT and CS%diag%axesTL are not allocated.
+
+ ! At this point G%mask2dT and CS%diag%axesTL are not allocated.
! postpone diag_registeration to initialize_MOM_generic_tracer
!Fields cannot be diag registered as they are allocated and have to registered later.
@@ -216,7 +218,7 @@ function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
! Initialize all generic tracers
!
call generic_tracer_init(HI%isc,HI%iec,HI%jsc,HI%jec,HI%isd,HI%ied,HI%jsd,HI%jed,&
- GV%ke,ntau,axes,grid_tmask,grid_kmt,set_time(0,0))
+ GV%ke,ntau,axes,grid_tmask,grid_kmt,set_time(0,0))
!
@@ -288,7 +290,7 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h
type(param_file_type), intent(in) :: param_file
- type(diag_ctrl), target, intent(in) :: diag
+ type(diag_ctrl), target, intent(in) :: diag
type(ocean_OBC_type), pointer :: OBC
type(MOM_generic_tracer_CS), pointer :: CS
type(sponge_CS), pointer :: sponge_CSp
@@ -303,7 +305,6 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia
! (in) G - The ocean's grid structure.
! (in) GV - The ocean's vertical grid structure.
! (in) h - Layer thickness, in m or kg m-2.
- ! (in) diag - A structure that is used to regulate diagnostic output.
! (in) OBC - This open boundary condition type specifies whether, where,
! and what open boundary conditions are used.
! (in/out) CS - The control structure returned by a previous call to
@@ -323,24 +324,22 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia
character(len=fm_string_len) :: g_tracer_name, longname, units
real, dimension(:,:,:,:), pointer :: tr_field
real, dimension(:,:,:), pointer :: tr_ptr
- real, dimension(G%isd:G%ied, G%jsd:G%jed,1:GV%ke) :: grid_tmask
- integer, dimension(G%isd:G%ied, G%jsd:G%jed) :: grid_kmt
+ real, dimension(G%isd:G%ied, G%jsd:G%jed,1:G%ke) :: grid_tmask
+ integer, dimension(G%isd:G%ied, G%jsd:G%jed) :: grid_kmt
!! 2010/02/04 Add code to re-initialize Generic Tracers if needed during a model simulation
!! By default, restart cpio should not contain a Generic Tracer IC file and step below will be skipped.
!! Ideally, the generic tracer IC file should have the tracers on Z levels.
- isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke
+ isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = G%ke
+ CS%diag=>diag
!Get the tracer list
if(.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//&
": No tracer in the list.")
!For each tracer name get its fields
g_tracer=>CS%g_tracer_list
- ! diag is an opaque type that contains information about the diagnostics.
- CS%diag => diag
-
do
if(INDEX(CS%IC_file, '_NULL_') .ne. 0) then
call MOM_error(WARNING,"The name of the IC_file "//trim(CS%IC_file)//&
@@ -433,15 +432,14 @@ subroutine initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, dia
do j = G%jsd, G%jed ; do i = G%isd, G%ied
if (G%mask2dT(i,j) .gt. 0) then
grid_tmask(i,j,:) = 1.0
- grid_kmt(i,j) = GV%ke ! Tell the code that a layer thicker than 1m is the bottom layer.
+ grid_kmt(i,j) = G%ke ! Tell the code that a layer thicker than 1m is the bottom layer.
endif
enddo ; enddo
-
call g_tracer_set_common(G%isc,G%iec,G%jsc,G%jec,G%isd,G%ied,G%jsd,G%jed,&
GV%ke,1,CS%diag%axesTL%handles,grid_tmask,grid_kmt,day)
! Register generic tracer modules diagnostics
-
+
#ifdef _USE_MOM6_DIAG
call g_tracer_set_csdiag(CS%diag)
#endif
@@ -510,7 +508,8 @@ end subroutine initialize_MOM_generic_tracer
!
!
- subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, tv, optics)
+ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, tv, optics, &
+ evap_CFL_limit, minimum_forcing_depth)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old, h_new, ea, eb
@@ -519,6 +518,8 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G
type(MOM_generic_tracer_CS), pointer :: CS
type(thermo_var_ptrs), intent(in) :: tv
type(optics_type), intent(in) :: optics
+ real, optional,intent(in) :: evap_CFL_limit
+ real, optional,intent(in) :: minimum_forcing_depth
! This subroutine applies diapycnal diffusion and any other column
! tracer physics or chemistry to the tracers from this file.
! CFCs are relatively simple, as they are passive tracers. with only a surface
@@ -539,6 +540,10 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G
! (in) GV - The ocean's vertical grid structure.
! (in) CS - The control structure returned by a previous call to
! register_MOM_generic_tracer.
+ ! (in) evap_CFL_limit - Limits how much water can be fluxed out of the top layer
+ ! Stored previously in diabatic CS.
+ ! (in) minimum_forcing_depth - The smallest depth over which fluxes can be applied
+ ! Stored previously in diabatic CS.
!
! 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]
@@ -548,14 +553,12 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G
character(len=fm_string_len) :: g_tracer_name
real, dimension(:,:), pointer :: stf_array,trunoff_array,runoff_tracer_flux_array
- real :: surface_field(SZI_(G),SZJ_(G))
- real :: sosga
-
- real, dimension(G%isd:G%ied,G%jsd:G%jed,GV%ke) :: rho_dzt, dzt
- real, dimension(G%isd:G%ied,G%jsd:G%jed) :: hblt_depth
+ real, dimension(G%isd:G%ied,G%jsd:G%jed,G%ke) :: rho_dzt, dzt
+ real, dimension(G%isd:G%ied,G%jsd:G%jed) :: hblt_depth
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work
integer :: i, j, k, isc, iec, jsc, jec, nk
- isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke
+ isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = G%ke
!Get the tracer list
if(.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL,&
@@ -620,27 +623,44 @@ subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G
enddo; enddo ; enddo
- do j=jsc,jec ; do i=isc,iec
- surface_field(i,j) = tv%S(i,j,1)
- enddo ; enddo
- sosga = global_area_mean(surface_field, G)
-
!
!Calculate tendencies (i.e., field changes at dt) from the sources / sinks
!
- call generic_tracer_source(tv%T,tv%S,sosga,rho_dzt,dzt,hblt_depth,G%isd,G%jsd,1,dt,&
+ call generic_tracer_source(tv%T,tv%S,rho_dzt,dzt,hblt_depth,G%isd,G%jsd,1,dt,&
G%areaT,get_diag_time_end(CS%diag),&
optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band)
+ ! This uses applyTracerBoundaryFluxesInOut to handle the change in tracer due to freshwater fluxes
+ ! usually in ALE mode
+ if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
+ g_tracer=>CS%g_tracer_list
+ do
+ if (g_tracer_is_prog(g_tracer)) then
+ do k=1,nk ;do j=jsc,jec ; do i=isc,iec
+ h_work(i,j,k) = h_old(i,j,k)
+ enddo ; enddo ; enddo;
+ call applyTracerBoundaryFluxesInOut(G, GV, g_tracer%field(:,:,:,1), dt, fluxes, h_work, &
+ evap_CFL_limit, minimum_forcing_depth)
+ endif
+
+ !traverse the linked list till hit NULL
+ call g_tracer_get_next(g_tracer, g_tracer_next)
+ if(.NOT. associated(g_tracer_next)) exit
+ g_tracer=>g_tracer_next
+ enddo
+ endif
+
!
!Update Tr(n)%field from explicit vertical diffusion
!
-
! Use a tridiagonal solver to determine the concentrations after the
! surface source is applied and diapycnal advection and diffusion occurs.
-
- call generic_tracer_vertdiff_G(h_old, ea, eb, dt, GV%kg_m2_to_H, GV%m_to_H, 1) !Last arg is tau which is always 1 for MOM
+ if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
+ call generic_tracer_vertdiff_G(h_work, ea, eb, dt, GV%kg_m2_to_H, GV%m_to_H, 1) !Last arg is tau which is always 1 for MOM
+ else
+ call generic_tracer_vertdiff_G(h_old, ea, eb, dt, GV%kg_m2_to_H, GV%m_to_H, 1) !Last arg is tau which is always 1 for MOM
+ endif
! Update bottom fields after vertical processes
@@ -701,7 +721,7 @@ function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_inde
character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_stock'
integer :: i, j, k, is, ie, js, je, nz, m
- is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
MOM_generic_tracer_stock = 0
if (.not.associated(CS)) return
@@ -866,8 +886,6 @@ subroutine MOM_generic_tracer_surface_state(state, h, G, CS)
! (in) CS - The control structure returned by a previous call to
! register_MOM_generic_tracer.
- real :: sosga
-
character(len=fm_string_len), parameter :: sub_name = 'MOM_generic_tracer_surface_state'
real, dimension(G%isd:G%ied,G%jsd:G%jed,1:G%ke,1) :: rho0
type(g_tracer_type), pointer :: g_tracer
@@ -876,15 +894,12 @@ subroutine MOM_generic_tracer_surface_state(state, h, G, CS)
!nnz: fake rho0
rho0=1.0
- sosga = global_area_mean(state%SSS, G)
-
call generic_tracer_coupler_set(state%tr_fields,&
ST=state%SST,&
SS=state%SSS,&
- sosga=sosga, &
- rho=rho0,& !nnz: required for MOM5 and previous versions.
+ rho=rho0,& !nnz: required for MOM
ilb=G%isd, jlb=G%jsd,&
- tau=1,model_time=get_diag_time_end(CS%diag))
+ tau=1)
!Output diagnostics via diag_manager for all tracers in this module
! if(.NOT. associated(CS%g_tracer_list)) call mpp_error(FATAL, trim(sub_name)//&
diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90
new file mode 100644
index 0000000000..b11c1b5ec6
--- /dev/null
+++ b/src/tracer/MOM_tracer_diabatic.F90
@@ -0,0 +1,423 @@
+!> This module contains routines that implement physical fluxes of tracers (e.g. due
+!! to surface fluxes or mixing). These are intended to be called from call_tracer_column_fns
+!! in the MOM_tracer_flow_control module.
+module MOM_tracer_diabatic
+
+! This file is part of MOM6. See LICENSE.md for the license.
+use MOM_grid, only : ocean_grid_type
+use MOM_verticalGrid, only : verticalGrid_type
+use MOM_forcing_type, only : forcing
+use MOM_error_handler, only : MOM_error, FATAL, WARNING
+
+implicit none ; private
+
+#include
+public tracer_vertdiff
+public applyTracerBoundaryFluxesInOut
+!> This subroutine solves a tridiagonal equation for the final tracer
+!! concentrations after the dual-entrainments, and possibly sinking or surface
+!! and bottom sources, are applied. The sinking is implemented with an
+!! fully implicit upwind advection scheme.
+
+contains
+
+subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
+ sfc_flux, btm_flux, btm_reservoir, sink_rate, convert_flux_in)
+ type(ocean_grid_type), intent(in) :: G !< ocean grid structure
+ type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
+ real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< layer thickness before entrainment (m or kg m-2)
+ real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: ea !< amount of fluid entrained from the layer above (units of h_work)
+ real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: eb !< amount of fluid entrained from the layer below (units of h_work)
+ real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr !< tracer concentration (in concentration units CU)
+ real, intent(in) :: dt !< amount of time covered by this call (seconds)
+ real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: sfc_flux !< surface flux of the tracer (in CU * kg m-2 s-1)
+ real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: btm_flux !< The (negative upward) bottom flux of the tracer,
+ !! in units of (CU * kg m-2 s-1)
+ real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir (units of CU kg m-2; formerly CU m)
+ real, optional,intent(in) :: sink_rate !< rate at which the tracer sinks, in m s-1
+ logical, optional,intent(in) :: convert_flux_in !< True if the specified sfc_flux needs to be integrated in time
+
+ real :: sink_dist ! The distance the tracer sinks in a time step, in m or kg m-2.
+ real, dimension(SZI_(G),SZJ_(G)) :: &
+ sfc_src, & ! The time-integrated surface source of the tracer, in
+ ! units of m or kg m-2 times a concentration.
+ btm_src ! The time-integrated bottom source of the tracer, in
+ ! units of m or kg m-2 times a concentration.
+ real, dimension(SZI_(G)) :: &
+ b1, & ! b1 is used by the tridiagonal solver, in m-1 or m2 kg-1.
+ d1 ! d1=1-c1 is used by the tridiagonal solver, nondimensional.
+ real :: c1(SZI_(G),SZK_(GV)) ! c1 is used by the tridiagonal solver, ND.
+ real :: h_minus_dsink(SZI_(G),SZK_(GV)) ! The layer thickness minus the
+ ! difference in sinking rates across the layer, in m or kg m-2.
+ ! By construction, 0 <= h_minus_dsink < h_work.
+ real :: sink(SZI_(G),SZK_(GV)+1) ! The tracer's sinking distances at the
+ ! interfaces, limited to prevent characteristics from
+ ! crossing within a single timestep, in m or kg m-2.
+ real :: b_denom_1 ! The first term in the denominator of b1, in m or kg m-2.
+ real :: h_tr ! h_tr is h at tracer points with a h_neglect added to
+ ! ensure positive definiteness, in m or kg m-2.
+ real :: h_neglect ! A thickness that is so small it is usually lost
+ ! in roundoff and can be neglected, in m.
+ logical :: convert_flux = .true.
+
+
+ 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 (present(convert_flux_in)) convert_flux = convert_flux_in
+ h_neglect = GV%H_subroundoff
+ sink_dist = 0.0
+ if (present(sink_rate)) sink_dist = (dt*sink_rate) * GV%m_to_H
+!$OMP parallel default(none) shared(is,ie,js,je,sfc_src,btm_src,sfc_flux,dt,G,GV,btm_flux, &
+!$OMP sink_rate,btm_reservoir,nz,sink_dist,h_work,ea, &
+!$OMP h_neglect,eb,tr) &
+!$OMP private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1)
+!$OMP do
+ do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo; enddo
+ if (present(sfc_flux)) then
+ if(convert_flux) then
+!$OMP do
+ do j = js, je; do i = is,ie
+ sfc_src(i,j) = (sfc_flux(i,j)*dt) * GV%kg_m2_to_H
+ enddo; enddo
+ else
+ do j = js, je; do i = is,ie
+ sfc_src(i,j) = sfc_flux(i,j)
+ enddo; enddo
+ endif
+ endif
+ if (present(btm_flux)) then
+!$OMP do
+ if(convert_flux) then
+ do j = js, je; do i = is,ie
+ btm_src(i,j) = (btm_flux(i,j)*dt) * GV%kg_m2_to_H
+ enddo; enddo
+ else
+ do j = js, je; do i = is,ie
+ btm_src(i,j) = btm_flux(i,j)
+ enddo; enddo
+ endif
+ endif
+
+ if (present(sink_rate)) then
+!$OMP do
+ do j=js,je
+ ! Find the sinking rates at all interfaces, limiting them if necesary
+ ! so that the characteristics do not cross within a timestep.
+ ! If a non-constant sinking rate were used, that would be incorprated
+ ! here.
+ if (present(btm_reservoir)) then
+ do i=is,ie ; sink(i,nz+1) = sink_dist ; enddo
+ do k=2,nz ; do i=is,ie
+ sink(i,K) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k)
+ enddo ; enddo
+ else
+ do i=is,ie ; sink(i,nz+1) = 0.0 ; enddo
+ ! Find the limited sinking distance at the interfaces.
+ do k=nz,2,-1 ; do i=is,ie
+ if (sink(i,K+1) >= sink_dist) then
+ sink(i,K) = sink_dist
+ h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,K+1) - sink(i,K))
+ elseif (sink(i,K+1) + h_old(i,j,k) < sink_dist) then
+ sink(i,K) = sink(i,K+1) + h_old(i,j,k)
+ h_minus_dsink(i,k) = 0.0
+ else
+ sink(i,K) = sink_dist
+ h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,K+1)) - sink(i,K)
+ endif
+ enddo ; enddo
+ endif
+ do i=is,ie
+ sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2))
+ enddo
+
+ ! Now solve the tridiagonal equation for the tracer concentrations.
+ do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+ b_denom_1 = h_minus_dsink(i,1) + ea(i,j,1) + h_neglect
+ 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)
+ 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)
+ b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ea(i,j,k) + sink(i,K)) + &
+ h_neglect
+ b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
+ d1(i) = b_denom_1 * b1(i)
+ h_tr = h_old(i,j,k) + h_neglect
+ tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + &
+ (ea(i,j,k) + sink(i,K)) * tr(i,j,k-1))
+ endif ; enddo ; enddo
+ do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
+ c1(i,nz) = eb(i,j,nz-1) * b1(i)
+ b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ea(i,j,nz) + sink(i,nz)) + &
+ h_neglect
+ b1(i) = 1.0 / (b_denom_1 + eb(i,j,nz))
+ h_tr = h_old(i,j,nz) + h_neglect
+ tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + &
+ (ea(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1))
+ endif ; enddo
+ if (present(btm_reservoir)) then ; do i=is,ie ; if (G%mask2dT(i,j)>0.5) then
+ btm_reservoir(i,j) = btm_reservoir(i,j) + &
+ (sink(i,nz+1)*tr(i,j,nz)) * GV%H_to_kg_m2
+ endif ; enddo ; endif
+
+ 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
+ else
+!$OMP do
+ do j=js,je
+ 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)
+ 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)
+ h_tr = h_old(i,j,k) + h_neglect
+ b_denom_1 = h_tr + d1(i) * ea(i,j,k)
+ b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
+ 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
+ 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)
+ b1(i) = 1.0 / ( b_denom_1 + eb(i,j,nz))
+ 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
+ tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
+ endif ; enddo ; enddo
+ enddo
+ endif
+
+!$OMP end parallel
+
+end subroutine tracer_vertdiff
+
+subroutine applyTracerBoundaryFluxesInOut(G, GV, Tr, dt, fluxes, h, &
+ evap_CFL_limit, minimum_forcing_depth, in_flux_optional, out_flux_optional)
+! This routine is modeled after applyBoundaryFluxesInOut in MOM_diabatic_aux.F90
+! NOTE: Please note that in this routine sfc_flux gets set to zero to ensure that the surface
+! flux of the tracer does not get applied again during a subsequent call to tracer_vertdif
+
+ type(ocean_grid_type), intent(in) :: G !< Grid structure
+ type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: Tr !< Tracer concentration on T-cell
+ real, intent(in) :: dt !< Time-step over which forcing is applied (s)
+ type(forcing), intent(in) :: fluxes !< Surface fluxes container
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness in H units
+ real, intent(in) :: evap_CFL_limit
+ real, intent(in) :: minimum_forcing_depth
+ real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: in_flux_optional ! The total time-integrated amount of tracer!
+ ! that enters with freshwater
+ real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: out_flux_optional ! The total time-integrated amount of tracer!
+ ! that leaves with freshwater
+ integer, parameter :: maxGroundings = 5
+ integer :: numberOfGroundings, iGround(maxGroundings), jGround(maxGroundings)
+ real :: H_limit_fluxes, IforcingDepthScale, Idt
+ real :: dThickness, dTracer
+ real :: fractionOfForcing, hOld, Ithickness
+ real :: RivermixConst ! A constant used in implementing river mixing, in Pa s.
+ real, dimension(SZI_(G)) :: &
+ netMassInOut, & ! surface water fluxes (H units) over time step
+ netMassIn, & ! mass entering ocean surface (H units) over a time step
+ netMassOut ! mass leaving ocean surface (H units) over a time step
+
+ real, dimension(SZI_(G), SZK_(G)) :: h2d, Tr2d
+ real, dimension(SZI_(G),SZJ_(G)) :: in_flux ! The total time-integrated amount of tracer!
+ ! that enters with freshwater
+ real, dimension(SZI_(G),SZJ_(G)) :: out_flux ! The total time-integrated amount of tracer!
+ ! that leaves with freshwater
+ real, dimension(SZI_(G)) :: in_flux_1d, out_flux_1d
+ real :: hGrounding(maxGroundings)
+ real :: Tr_in
+ integer :: i, j, is, ie, js, je, k, nz, n, nsw
+ character(len=45) :: mesg
+
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke
+
+ ! Only apply forcing if fluxes%sw is associated.
+ if (.not.ASSOCIATED(fluxes%sw)) return
+
+ in_flux(:,:) = 0.0 ; out_flux(:,:) = 0.0
+ if(present(in_flux_optional)) then
+ do j=js,je ; do i=is,ie
+ in_flux(i,j) = in_flux_optional(i,j)
+ enddo; enddo
+ endif
+ if(present(out_flux_optional)) then
+ do j=js,je ; do i=is,ie
+ out_flux(i,j) = out_flux_optional(i,j)
+ enddo ; enddo
+ endif
+
+ Idt = 1.0/dt
+ numberOfGroundings = 0
+
+!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,Tr,G,GV,fluxes,dt, &
+!$OMP IforcingDepthScale, &
+!$OMP numberOfGroundings,iGround,jGround, &
+!$OMP hGrounding,CS,Idt,aggregate_FW_forcing) &
+!$OMP private(h2d,Tr2d,netMassInOut,netMassOut, &
+!$OMP sfc_src,fractionOfForcing, &
+!$OMP dThickness,dTracer,hOld,Ithickness, &
+!$OMP netMassIn, Tr_in, Salin_in)
+
+ ! Work in vertical slices for efficiency
+ do j=js,je
+
+ ! Copy state into 2D-slice arrays
+ do k=1,nz ; do i=is,ie
+ h2d(i,k) = h(i,j,k)
+ Tr2d(i,k) = Tr(i,j,k)
+ enddo ; enddo
+
+ do i = is,ie
+ in_flux_1d(i) = in_flux(i,j)
+ out_flux_1d(i) = out_flux(i,j)
+ enddo
+ ! The surface forcing is contained in the fluxes type.
+ ! We aggregate the thermodynamic forcing for a time step into the following:
+ ! These should have been set and stored during a call to applyBoundaryFluxesInOut
+ ! netMassIn = net mass entering at ocean surface over a timestep
+ ! netMassOut = net mass leaving ocean surface (H units) over a time step.
+ ! netMassOut < 0 means mass leaves ocean.
+
+ ! Note here that the aggregateFW flag has already been taken care of in the call to
+ ! applyBoundaryFluxesInOut
+ do i=is,ie
+ netMassOut(i) = fluxes%netMassOut(i,j)
+ netMassIn(i) = fluxes%netMassIn(i,j)
+ enddo
+
+ ! Apply the surface boundary fluxes in three steps:
+ ! A/ update concentration from mass entering the ocean
+ ! B/ update concentration from mass leaving ocean.
+ do i=is,ie
+ if (G%mask2dT(i,j)>0.) then
+
+ ! A/ Update tracer due to incoming mass flux.
+ do k=1,1
+
+ ! Change in state due to forcing
+ dThickness = netMassIn(i) ! Since we are adding mass, we can use all of it
+ dTracer = 0.
+
+ ! Update the forcing by the part to be consumed within the present k-layer.
+ ! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
+ netMassIn(i) = netMassIn(i) - dThickness
+ dTracer = dTracer + in_flux_1d(i)
+ in_flux_1d(i) = 0.0
+
+ ! Update state
+ hOld = h2d(i,k) ! Keep original thickness in hand
+ h2d(i,k) = h2d(i,k) + dThickness ! New thickness
+ if (h2d(i,k) > 0.0) then
+ Ithickness = 1.0/h2d(i,k) ! Inverse new thickness
+ ! The "if"s below avoid changing T/S by roundoff unnecessarily
+ if (dThickness /= 0. .or. dTracer /= 0.) tr2d(i,k) = (hOld*tr2d(i,k)+ dTracer)*Ithickness
+ endif
+
+ enddo ! k=1,1
+
+ ! B/ Update tracer from mass leaving ocean
+ do k=1,nz
+
+ ! Place forcing into this layer if this layer has nontrivial thickness.
+ ! For layers thin relative to 1/IforcingDepthScale, then distribute
+ ! forcing into deeper layers.
+ IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth*GV%m_to_H - netMassOut(i) )
+ ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
+ fractionOfForcing = min(1.0, h2d(i,k)*IforcingDepthScale)
+
+ ! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we
+ ! limit the forcing applied to this cell, leaving the remaining forcing to
+ ! be distributed downwards.
+ if (-fractionOfForcing*netMassOut(i) > evap_CFL_limit*h2d(i,k)) then
+ fractionOfForcing = -evap_CFL_limit*h2d(i,k)/netMassOut(i)
+ endif
+
+ ! Change in state due to forcing
+ dThickness = max( fractionOfForcing*netMassOut(i), -h2d(i,k) )
+ ! Note this is slightly different to how salt is currently treated
+ dTracer = fractionOfForcing*out_flux_1d(i)
+
+ ! Update the forcing by the part to be consumed within the present k-layer.
+ ! If fractionOfForcing = 1, then new netMassOut vanishes.
+ netMassOut(i) = netMassOut(i) - dThickness
+ out_flux_1d(i) = out_flux_1d(i) - dTracer
+
+ ! Update state by the appropriate increment.
+ hOld = h2d(i,k) ! Keep original thickness in hand
+ h2d(i,k) = h2d(i,k) + dThickness ! New thickness
+ if (h2d(i,k) > 0.) then
+ Ithickness = 1.0/h2d(i,k) ! Inverse of new thickness
+ Tr2d(i,k) = (hOld*Tr2d(i,k) + dTracer)*Ithickness
+ elseif (h2d(i,k) < 0.0) then ! h2d==0 is a special limit that needs no extra handling
+ write(0,*) 'applyTracerBoundaryFluxesInOut(): lon,lat=',G%geoLonT(i,j),G%geoLatT(i,j)
+ write(0,*) 'applyTracerBoundaryFluxesInOut(): netTr,netH=',in_flux_1d(i)-out_flux_1d(i),netMassInOut(i)
+ write(0,*) 'applyTracerBoundaryFluxesInOut(): h(n),h(n+1),k=',hOld,h2d(i,k),k
+ call MOM_error(FATAL, "MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//&
+ "Complete mass loss in column!")
+ endif
+
+ enddo ! k
+
+ ! Check if trying to apply fluxes over land points
+ elseif((abs(in_flux_1d(i))+abs(out_flux_1d(i))+abs(netMassIn(i))+abs(netMassOut(i)))>0.) then
+ write(0,*) 'applyTracerBoundaryFluxesInOut(): lon,lat=',G%geoLonT(i,j),G%geoLatT(i,j)
+ write(0,*) 'applyTracerBoundaryFluxesInOut(): in_flux, out_flux, netMassIn,netMassOut=',&
+ in_flux_1d(i), out_flux_1d(i),netMassIn(i),netMassOut(i)
+ call MOM_error(FATAL, "MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//&
+ "Mass loss over land?")
+ endif
+
+ ! If anything remains after the k-loop, then we have grounded out, which is a problem.
+ if (abs(in_flux_1d(i))+abs(out_flux_1d(i)) /= 0.0) then
+!$OMP critical
+ numberOfGroundings = numberOfGroundings +1
+ if (numberOfGroundings<=maxGroundings) then
+ iGround(numberOfGroundings) = i ! Record i,j location of event for
+ jGround(numberOfGroundings) = j ! warning message
+ hGrounding(numberOfGroundings) = abs(in_flux_1d(i))+abs(out_flux_1d(i))
+ endif
+!$OMP end critical
+ endif
+
+ enddo ! i
+
+ ! Step C/ copy updated tracer concentration from the 2d slice now back into model state.
+ do k=1,nz ; do i=is,ie
+ Tr(i,j,k) = Tr2d(i,k)
+ h(i,j,k) = h2d(i,k)
+ enddo ; enddo
+
+ enddo ! j-loop finish
+
+ if (numberOfGroundings>0) then
+ do i = 1, min(numberOfGroundings, maxGroundings)
+ write(mesg(1:45),'(3es15.3)') G%geoLonT( iGround(i), jGround(i) ), &
+ G%geoLatT( iGround(i), jGround(i)) , hGrounding(i)
+ call MOM_error(WARNING, "MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//&
+ "Tracer created. x,y,dh= "//trim(mesg), all_print=.true.)
+ enddo
+
+ if (numberOfGroundings - maxGroundings > 0) then
+ write(mesg, '(i4)') numberOfGroundings - maxGroundings
+ call MOM_error(WARNING, "MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//&
+ trim(mesg) // " groundings remaining")
+ endif
+ endif
+
+end subroutine applyTracerBoundaryFluxesInOut
+end module MOM_tracer_diabatic
diff --git a/src/tracer/MOM_tracer_flow_control.F90 b/src/tracer/MOM_tracer_flow_control.F90
index 3d43765d9d..a07d05e61e 100644
--- a/src/tracer/MOM_tracer_flow_control.F90
+++ b/src/tracer/MOM_tracer_flow_control.F90
@@ -70,13 +70,17 @@ module MOM_tracer_flow_control
use oil_tracer, only : oil_stock, oil_tracer_end, oil_tracer_CS
use advection_test_tracer, only : register_advection_test_tracer, initialize_advection_test_tracer
use advection_test_tracer, only : advection_test_tracer_column_physics, advection_test_tracer_surface_state
-use advection_test_tracer, only : advection_test_tracer_end, advection_test_tracer_CS
+use advection_test_tracer, only : advection_test_stock, advection_test_tracer_end, advection_test_tracer_CS
#ifdef _USE_GENERIC_TRACER
use MOM_generic_tracer, only : register_MOM_generic_tracer, initialize_MOM_generic_tracer
use MOM_generic_tracer, only : MOM_generic_tracer_column_physics, MOM_generic_tracer_surface_state
use MOM_generic_tracer, only : end_MOM_generic_tracer, MOM_generic_tracer_get
use MOM_generic_tracer, only : MOM_generic_tracer_stock, MOM_generic_tracer_min_max, MOM_generic_tracer_CS
#endif
+use pseudo_salt_tracer, only : register_pseudo_salt_tracer, initialize_pseudo_salt_tracer
+use pseudo_salt_tracer, only : pseudo_salt_tracer_column_physics, pseudo_salt_tracer_surface_state
+use pseudo_salt_tracer, only : pseudo_salt_stock, pseudo_salt_tracer_end, pseudo_salt_tracer_CS
+
implicit none ; private
@@ -94,6 +98,7 @@ module MOM_tracer_flow_control
logical :: use_advection_test_tracer = .false.
logical :: use_OCMIP2_CFC = .false.
logical :: use_MOM_generic_tracer = .false.
+ logical :: use_pseudo_salt_tracer = .false.
type(USER_tracer_example_CS), pointer :: USER_tracer_example_CSp => NULL()
type(DOME_tracer_CS), pointer :: DOME_tracer_CSp => NULL()
type(ISOMIP_tracer_CS), pointer :: ISOMIP_tracer_CSp => NULL()
@@ -105,6 +110,7 @@ module MOM_tracer_flow_control
#ifdef _USE_GENERIC_TRACER
type(MOM_generic_tracer_CS), pointer :: MOM_generic_tracer_CSp => NULL()
#endif
+ type(pseudo_salt_tracer_CS), pointer :: pseudo_salt_tracer_CSp => NULL()
end type tracer_flow_control_CS
contains
@@ -172,6 +178,9 @@ subroutine call_tracer_register(HI, GV, param_file, CS, tr_Reg, restart_CS)
"If true and _USE_GENERIC_TRACER is defined as a \n"//&
"preprocessor macro, use the MOM_generic_tracer packages.", &
default=.false.)
+ call get_param(param_file, mod, "USE_PSEUDO_SALT_TRACER", CS%use_pseudo_salt_tracer, &
+ "If true, use the pseudo salt tracer, typically run as a diagnostic.", &
+ default=.false.)
#ifndef _USE_GENERIC_TRACER
if (CS%use_MOM_generic_tracer) call MOM_error(FATAL, &
@@ -210,11 +219,14 @@ subroutine call_tracer_register(HI, GV, param_file, CS, tr_Reg, restart_CS)
register_MOM_generic_tracer(HI, GV, param_file, CS%MOM_generic_tracer_CSp, &
tr_Reg, restart_CS)
#endif
+ if (CS%use_pseudo_salt_tracer) CS%use_pseudo_salt_tracer = &
+ register_pseudo_salt_tracer(HI, GV, param_file, CS%pseudo_salt_tracer_CSp, &
+ tr_Reg, restart_CS)
end subroutine call_tracer_register
subroutine tracer_flow_control_init(restart, day, G, GV, h, param_file, diag, OBC, &
- CS, sponge_CSp, ALE_sponge_CSp, diag_to_Z_CSp)
+ CS, sponge_CSp, ALE_sponge_CSp, diag_to_Z_CSp, tv)
logical, intent(in) :: restart
type(time_type), target, intent(in) :: day
type(ocean_grid_type), intent(inout) :: G
@@ -227,6 +239,7 @@ subroutine tracer_flow_control_init(restart, day, G, GV, h, param_file, diag, OB
type(sponge_CS), pointer :: sponge_CSp
type(ALE_sponge_CS), pointer :: ALE_sponge_CSp
type(diag_to_Z_CS), pointer :: diag_to_Z_CSp
+ type(thermo_var_ptrs), intent(in) :: tv
! This subroutine calls all registered tracer initialization
! subroutines.
@@ -279,6 +292,9 @@ subroutine tracer_flow_control_init(restart, day, G, GV, h, param_file, diag, OB
call initialize_MOM_generic_tracer(restart, day, G, GV, h, param_file, diag, OBC, &
CS%MOM_generic_tracer_CSp, sponge_CSp, ALE_sponge_CSp, diag_to_Z_CSp)
#endif
+ if (CS%use_pseudo_salt_tracer) &
+ call initialize_pseudo_salt_tracer(restart, day, G, GV, h, diag, OBC, CS%pseudo_salt_tracer_CSp, &
+ sponge_CSp, diag_to_Z_CSp, tv)
end subroutine tracer_flow_control_init
@@ -338,7 +354,8 @@ subroutine call_tracer_set_forcing(state, fluxes, day_start, day_interval, G, CS
end subroutine call_tracer_set_forcing
-subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, dt, G, GV, tv, optics, CS)
+subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, dt, G, GV, tv, optics, CS, &
+ debug, evap_CFL_limit, minimum_forcing_depth)
real, dimension(NIMEM_,NJMEM_,NKMEM_), intent(in) :: h_old, h_new, ea, eb
type(forcing), intent(in) :: fluxes
real, intent(in) :: dt
@@ -347,6 +364,10 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, dt, G, GV, tv, o
type(thermo_var_ptrs), intent(in) :: tv
type(optics_type), pointer :: optics
type(tracer_flow_control_CS), pointer :: CS
+ logical, intent(in) :: debug
+ real, optional,intent(in) :: evap_CFL_limit
+ real, optional,intent(in) :: minimum_forcing_depth
+
! This subroutine calls all registered tracer column physics
! subroutines.
@@ -368,39 +389,107 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, dt, G, GV, tv, o
! (in) optics - The structure containing optical properties.
! (in) CS - The control structure returned by a previous call to
! call_tracer_register.
+! (in) evap_CFL_limit - Limits how much water can be fluxed out of the top layer
+! Stored previously in diabatic CS.
+! (in) minimum_forcing_depth - The smallest depth over which fluxes can be applied
+! Stored previously in diabatic CS.
+! (in) debug - Calculates checksums
if (.not. associated(CS)) call MOM_error(FATAL, "call_tracer_column_fns: "// &
"Module must be initialized via call_tracer_register before it is used.")
-! Add calls to tracer column functions here.
- if (CS%use_USER_tracer_example) &
- call tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
- G, GV, CS%USER_tracer_example_CSp)
- if (CS%use_DOME_tracer) &
- call DOME_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
- G, GV, CS%DOME_tracer_CSp)
- if (CS%use_ISOMIP_tracer) &
- call ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
- G, GV, CS%ISOMIP_tracer_CSp)
- if (CS%use_ideal_age) &
- call ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
- G, GV, CS%ideal_age_tracer_CSp)
- if (CS%use_regional_dyes) &
- call dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
- G, GV, CS%dye_tracer_CSp)
- if (CS%use_oil) &
- call oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
- G, GV, CS%oil_tracer_CSp, tv)
- if (CS%use_advection_test_tracer) &
- call advection_test_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
- G, GV, CS%advection_test_tracer_CSp)
- if (CS%use_OCMIP2_CFC) &
- call OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
- G, GV, CS%OCMIP2_CFC_CSp)
+
+ ! Use the applyTracerBoundaryFluxesInOut to handle surface fluxes
+ if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
+ ! Add calls to tracer column functions here.
+ if (CS%use_USER_tracer_example) &
+ call tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%USER_tracer_example_CSp)
+ if (CS%use_DOME_tracer) &
+ call DOME_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%DOME_tracer_CSp, &
+ evap_CFL_limit=evap_CFL_limit, &
+ minimum_forcing_depth=minimum_forcing_depth)
+ if (CS%use_ISOMIP_tracer) &
+ call ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%ISOMIP_tracer_CSp, &
+ evap_CFL_limit=evap_CFL_limit, &
+ minimum_forcing_depth=minimum_forcing_depth)
+ if (CS%use_ideal_age) &
+ call ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%ideal_age_tracer_CSp, &
+ evap_CFL_limit=evap_CFL_limit, &
+ minimum_forcing_depth=minimum_forcing_depth)
+ if (CS%use_regional_dyes) &
+ call dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%dye_tracer_CSp, &
+ evap_CFL_limit=evap_CFL_limit, &
+ minimum_forcing_depth=minimum_forcing_depth)
+ if (CS%use_oil) &
+ call oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%oil_tracer_CSp, tv, &
+ evap_CFL_limit=evap_CFL_limit, &
+ minimum_forcing_depth=minimum_forcing_depth)
+
+ if (CS%use_advection_test_tracer) &
+ call advection_test_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%advection_test_tracer_CSp, &
+ evap_CFL_limit=evap_CFL_limit, &
+ minimum_forcing_depth=minimum_forcing_depth)
+ if (CS%use_OCMIP2_CFC) &
+ call OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%OCMIP2_CFC_CSp, &
+ evap_CFL_limit=evap_CFL_limit, &
+ minimum_forcing_depth=minimum_forcing_depth)
#ifdef _USE_GENERIC_TRACER
- if (CS%use_MOM_generic_tracer) &
- call MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
- G, GV, CS%MOM_generic_tracer_CSp, tv, optics)
+ if (CS%use_MOM_generic_tracer) &
+ call MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%MOM_generic_tracer_CSp, tv, optics, &
+ evap_CFL_limit=evap_CFL_limit, &
+ minimum_forcing_depth=minimum_forcing_depth)
+#endif
+ if (CS%use_pseudo_salt_tracer) &
+ call pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%pseudo_salt_tracer_CSp, tv, debug,&
+ evap_CFL_limit=evap_CFL_limit, &
+ minimum_forcing_depth=minimum_forcing_depth)
+
+
+ else ! Apply tracer surface fluxes using ea on the first layer
+ if (CS%use_USER_tracer_example) &
+ call tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%USER_tracer_example_CSp)
+ if (CS%use_DOME_tracer) &
+ call DOME_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%DOME_tracer_CSp)
+ if (CS%use_ISOMIP_tracer) &
+ call ISOMIP_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%ISOMIP_tracer_CSp)
+ if (CS%use_ideal_age) &
+ call ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%ideal_age_tracer_CSp)
+ if (CS%use_regional_dyes) &
+ call dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%dye_tracer_CSp)
+ if (CS%use_oil) &
+ call oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%oil_tracer_CSp, tv)
+ if (CS%use_advection_test_tracer) &
+ call advection_test_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%advection_test_tracer_CSp)
+ if (CS%use_OCMIP2_CFC) &
+ call OCMIP2_CFC_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%OCMIP2_CFC_CSp)
+#ifdef _USE_GENERIC_TRACER
+ if (CS%use_MOM_generic_tracer) &
+ call MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%MOM_generic_tracer_CSp, tv, optics)
#endif
+ if (CS%use_pseudo_salt_tracer) &
+ call pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, &
+ G, GV, CS%pseudo_salt_tracer_CSp, tv, debug)
+
+
+ endif
end subroutine call_tracer_column_fns
@@ -486,6 +575,14 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, CS, stock_names, stock_uni
call store_stocks("MOM_OCMIP2_CFC", 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 )
+ call store_stocks("advection_test_tracer", ns, names, units, values, index, &
+ stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+ endif
+
#ifdef _USE_GENERIC_TRACER
if (CS%use_MOM_generic_tracer) then
ns = MOM_generic_tracer_stock(h, values, G, GV, CS%MOM_generic_tracer_CSp, &
@@ -498,6 +595,12 @@ subroutine call_tracer_stocks(h, stock_values, G, GV, CS, stock_names, stock_uni
endif
#endif
+ if (CS%use_pseudo_salt_tracer) then
+ ns = pseudo_salt_stock(h, values, G, GV, CS%pseudo_salt_tracer_CSp, &
+ names, units, stock_index)
+ call store_stocks("pseudo_salt_tracer", ns, names, units, values, index, &
+ stock_values, set_pkg_name, max_ns, ns_tot, stock_names, stock_units)
+ endif
if (ns_tot == 0) stock_values(1) = 0.0
@@ -610,6 +713,7 @@ subroutine tracer_flow_control_end(CS)
#ifdef _USE_GENERIC_TRACER
if (CS%use_MOM_generic_tracer) call end_MOM_generic_tracer(CS%MOM_generic_tracer_CSp)
#endif
+ if (CS%use_pseudo_salt_tracer) call pseudo_salt_tracer_end(CS%pseudo_salt_tracer_CSp)
if (associated(CS)) deallocate(CS)
end subroutine tracer_flow_control_end
diff --git a/src/tracer/MOM_tracer_registry.F90 b/src/tracer/MOM_tracer_registry.F90
index a956bf500e..132755def9 100644
--- a/src/tracer/MOM_tracer_registry.F90
+++ b/src/tracer/MOM_tracer_registry.F90
@@ -26,7 +26,6 @@ module MOM_tracer_registry
public add_tracer_diagnostics
public add_tracer_OBC_values
public lock_tracer_registry
-public tracer_vertdiff
public tracer_registry_end
!> The tracer type
@@ -247,180 +246,6 @@ subroutine add_tracer_diagnostics(name, Reg, ad_x, ad_y, df_x, df_y, &
end subroutine add_tracer_diagnostics
-
-!> This subroutine solves a tridiagonal equation for the final tracer
-!! concentrations after the dual-entrainments, and possibly sinking or surface
-!! and bottom sources, are applied. The sinking is implemented with an
-!! fully implicit upwind advection scheme.
-subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
- sfc_flux, btm_flux, btm_reservoir, sink_rate)
- type(ocean_grid_type), intent(in) :: G !< ocean grid structure
- type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
- real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< layer thickness before entrainment (m or kg m-2)
- real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: ea !< amount of fluid entrained from the layer above (units of h_old)
- real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: eb !< amount of fluid entrained from the layer below (units of h_old)
- real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr !< tracer concentration (in concentration units CU)
- real, intent(in) :: dt !< amount of time covered by this call (seconds)
- real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: sfc_flux !< surface flux of the tracer (in CU * kg m-2 s-1)
- real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: btm_flux !< The (negative upward) bottom flux of the tracer,
- !! in units of (CU * kg m-2 s-1)
- real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir (units of CU kg m-2; formerly CU m)
- real, optional,intent(in) :: sink_rate !< rate at which the tracer sinks, in m s-1
-
-
- real :: sink_dist ! The distance the tracer sinks in a time step, in m or kg m-2.
- real, dimension(SZI_(G),SZJ_(G)) :: &
- sfc_src, & ! The time-integrated surface source of the tracer, in
- ! units of m or kg m-2 times a concentration.
- btm_src ! The time-integrated bottom source of the tracer, in
- ! units of m or kg m-2 times a concentration.
- real, dimension(SZI_(G)) :: &
- b1, & ! b1 is used by the tridiagonal solver, in m-1 or m2 kg-1.
- d1 ! d1=1-c1 is used by the tridiagonal solver, nondimensional.
- real :: c1(SZI_(G),SZK_(GV)) ! c1 is used by the tridiagonal solver, ND.
- real :: h_minus_dsink(SZI_(G),SZK_(GV)) ! The layer thickness minus the
- ! difference in sinking rates across the layer, in m or kg m-2.
- ! By construction, 0 <= h_minus_dsink < h_old.
- real :: sink(SZI_(G),SZK_(GV)+1) ! The tracer's sinking distances at the
- ! interfaces, limited to prevent characteristics from
- ! crossing within a single timestep, in m or kg m-2.
- real :: b_denom_1 ! The first term in the denominator of b1, in m or kg m-2.
- real :: h_tr ! h_tr is h at tracer points with a h_neglect added to
- ! ensure positive definiteness, in m or kg m-2.
- real :: h_neglect ! A thickness that is so small it is usually lost
- ! in roundoff and can be neglected, in m.
-
- integer :: i, j, k, is, ie, js, je, nz
- is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
-
- h_neglect = GV%H_subroundoff
- sink_dist = 0.0
- if (present(sink_rate)) sink_dist = (dt*sink_rate) * GV%m_to_H
-!$OMP parallel default(none) shared(is,ie,js,je,sfc_src,btm_src,sfc_flux,dt,G,GV,btm_flux, &
-!$OMP sink_rate,btm_reservoir,nz,sink_dist,h_old,ea, &
-!$OMP h_neglect,eb,tr) &
-!$OMP private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1)
-!$OMP do
- do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo; enddo
- if (present(sfc_flux)) then
-!$OMP do
- do j = js, je; do i = is,ie
- sfc_src(i,j) = (sfc_flux(i,j)*dt) * GV%kg_m2_to_H
- enddo; enddo
- endif
- if (present(btm_flux)) then
-!$OMP do
- do j = js, je; do i = is,ie
- btm_src(i,j) = (btm_flux(i,j)*dt) * GV%kg_m2_to_H
- enddo; enddo
- endif
-
-
- if (present(sink_rate)) then
-!$OMP do
- do j=js,je
- ! Find the sinking rates at all interfaces, limiting them if necesary
- ! so that the characteristics do not cross within a timestep.
- ! If a non-constant sinking rate were used, that would be incorprated
- ! here.
- if (present(btm_reservoir)) then
- do i=is,ie ; sink(i,nz+1) = sink_dist ; enddo
- do k=2,nz ; do i=is,ie
- sink(i,K) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k)
- enddo ; enddo
- else
- do i=is,ie ; sink(i,nz+1) = 0.0 ; enddo
- ! Find the limited sinking distance at the interfaces.
- do k=nz,2,-1 ; do i=is,ie
- if (sink(i,K+1) >= sink_dist) then
- sink(i,K) = sink_dist
- h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,K+1) - sink(i,K))
- elseif (sink(i,K+1) + h_old(i,j,k) < sink_dist) then
- sink(i,K) = sink(i,K+1) + h_old(i,j,k)
- h_minus_dsink(i,k) = 0.0
- else
- sink(i,K) = sink_dist
- h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,K+1)) - sink(i,K)
- endif
- enddo ; enddo
- endif
- do i=is,ie
- sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2))
- enddo
-
- ! Now solve the tridiagonal equation for the tracer concentrations.
- do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
- b_denom_1 = h_minus_dsink(i,1) + ea(i,j,1) + h_neglect
- 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))
- 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)
- b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ea(i,j,k) + sink(i,K)) + &
- h_neglect
- b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
- d1(i) = b_denom_1 * b1(i)
- h_tr = h_old(i,j,k) + h_neglect
- tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + &
- (ea(i,j,k) + sink(i,K)) * tr(i,j,k-1))
- endif ; enddo ; enddo
- do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
- c1(i,nz) = eb(i,j,nz-1) * b1(i)
- b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ea(i,j,nz) + sink(i,nz)) + &
- h_neglect
- b1(i) = 1.0 / (b_denom_1 + eb(i,j,nz))
- h_tr = h_old(i,j,nz) + h_neglect
- tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + &
- (ea(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1))
- endif ; enddo
- if (present(btm_reservoir)) then ; do i=is,ie ; if (G%mask2dT(i,j)>0.5) then
- btm_reservoir(i,j) = btm_reservoir(i,j) + &
- (sink(i,nz+1)*tr(i,j,nz)) * GV%H_to_kg_m2
- endif ; enddo ; endif
-
- 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
- else
-!$OMP do
- do j=js,je
- 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) = b_denom_1 * b1(i)
- 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)
- h_tr = h_old(i,j,k) + h_neglect
- b_denom_1 = h_tr + d1(i) * ea(i,j,k)
- b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
- 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
- c1(i,nz) = eb(i,j,nz-1) * b1(i)
- h_tr = h_old(i,j,nz) + h_neglect
- b1(i) = 1.0 / (h_tr + d1(i) * ea(i,j,nz) + eb(i,j,nz))
- 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
- tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
- endif ; enddo ; enddo
- enddo
- endif
-
-!$OMP end parallel
-
-end subroutine tracer_vertdiff
-
-
!> This subroutine writes out chksums for thermodynamic state variables.
subroutine MOM_tracer_chksum(mesg, Tr, ntr, G)
character(len=*), intent(in) :: mesg !< message that appears on the chksum lines
diff --git a/src/tracer/advection_test_tracer.F90 b/src/tracer/advection_test_tracer.F90
index c5eae04fe2..46c7264757 100644
--- a/src/tracer/advection_test_tracer.F90
+++ b/src/tracer/advection_test_tracer.F90
@@ -68,7 +68,7 @@ module advection_test_tracer
use MOM_time_manager, only : time_type, get_time
use MOM_tracer_registry, only : register_tracer, tracer_registry_type
use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values
-use MOM_tracer_registry, only : tracer_vertdiff
+use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
use MOM_variables, only : surface
use MOM_verticalGrid, only : verticalGrid_type
@@ -81,7 +81,7 @@ module advection_test_tracer
public register_advection_test_tracer, initialize_advection_test_tracer
public advection_test_tracer_surface_state, advection_test_tracer_end
-public advection_test_tracer_column_physics
+public advection_test_tracer_column_physics, advection_test_stock
! ntr is the number of tracers in this module.
integer, parameter :: ntr = 11
@@ -116,6 +116,7 @@ module advection_test_tracer
integer, dimension(NTR) :: ind_tr ! Indices returned by aof_set_coupler_flux
! if it is used and the surface tracer concentrations are to be
! provided to the coupler.
+ integer :: ntr = NTR
type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
! timing of diagnostic output.
@@ -362,13 +363,16 @@ subroutine initialize_advection_test_tracer(restart, day, G, GV, h,diag, OBC, CS
end subroutine initialize_advection_test_tracer
-subroutine advection_test_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS)
+subroutine advection_test_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, &
+ evap_CFL_limit, minimum_forcing_depth)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old, h_new, ea, eb
type(forcing), intent(in) :: fluxes
real, intent(in) :: dt
type(advection_test_tracer_CS), pointer :: CS
+ real, optional,intent(in) :: evap_CFL_limit
+ real, optional,intent(in) :: minimum_forcing_depth
! This subroutine applies diapycnal diffusion and any other column
! tracer physics or chemistry to the tracers from this file.
! This is a simple example of a set of advected passive tracers.
@@ -391,7 +395,7 @@ subroutine advection_test_tracer_column_physics(h_old, h_new, ea, eb, fluxes,
!
! 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]
-
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
real :: b1(SZI_(G)) ! b1 and c1 are variables used by the
real :: c1(SZI_(G),SZK_(G)) ! tridiagonal solver.
integer :: i, j, k, is, ie, js, je, nz, m
@@ -399,9 +403,20 @@ subroutine advection_test_tracer_column_physics(h_old, h_new, ea, eb, fluxes,
if (.not.associated(CS)) return
-! do m=1,NTR
-! call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
-! enddo
+ if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
+ do m=1,CS%ntr
+ 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, CS%tr(:,:,:,m) , dt, fluxes, h_work, &
+ evap_CFL_limit, minimum_forcing_depth)
+ call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ else
+ do m=1,NTR
+ call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ endif
if (CS%mask_tracers) then
do m = 1,NTR ; if (CS%id_tracer(m) > 0) then
@@ -462,6 +477,59 @@ subroutine advection_test_tracer_surface_state(state, h, G, CS)
endif
end subroutine advection_test_tracer_surface_state
+function advection_test_stock(h, stocks, G, GV, CS, names, units, stock_index)
+ type(ocean_grid_type), intent(in) :: G
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h
+ real, dimension(:), intent(out) :: stocks
+ type(verticalGrid_type), intent(in) :: GV
+ type(advection_test_tracer_CS), pointer :: CS
+ character(len=*), dimension(:), intent(out) :: names
+ character(len=*), dimension(:), intent(out) :: units
+ integer, optional, intent(in) :: stock_index
+ integer :: advection_test_stock
+! This function 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.
+
+! Arguments: h - Layer thickness, in m or kg m-2.
+! (out) stocks - the mass-weighted integrated amount of each tracer,
+! in kg times concentration units.
+! (in) G - The ocean's grid structure.
+! (in) GV - The ocean's vertical grid structure.
+! (in) CS - The control structure returned by a previous call to
+! register_ideal_age_tracer.
+! (out) names - the names of the stocks calculated.
+! (out) units - the units of the stocks calculated.
+! (in,opt) stock_index - the coded index of a specific stock being sought.
+! Return value: the number of stocks calculated here.
+
+ integer :: i, j, k, is, ie, js, je, nz, m
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+
+ advection_test_stock = 0
+ if (.not.associated(CS)) return
+ if (CS%ntr < 1) 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
+
+ do m=1,CS%ntr
+ call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="advection_test_stock")
+ stocks(m) = 0.0
+ do k=1,nz ; do j=js,je ; do i=is,ie
+ stocks(m) = stocks(m) + CS%tr(i,j,k,m) * &
+ (G%mask2dT(i,j) * G%areaT(i,j) * h(i,j,k))
+ enddo ; enddo ; enddo
+ stocks(m) = GV%H_to_kg_m2 * stocks(m)
+ enddo
+ advection_test_stock = CS%ntr
+
+end function advection_test_stock
+
subroutine advection_test_tracer_end(CS)
type(advection_test_tracer_CS), pointer :: CS
integer :: m
diff --git a/src/tracer/dye_example.F90 b/src/tracer/dye_example.F90
index 4c4d7cc6bb..c8faba57ca 100644
--- a/src/tracer/dye_example.F90
+++ b/src/tracer/dye_example.F90
@@ -69,7 +69,7 @@ module regional_dyes
use MOM_time_manager, only : time_type, get_time
use MOM_tracer_registry, only : register_tracer, tracer_registry_type
use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values
-use MOM_tracer_registry, only : tracer_vertdiff
+use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
use MOM_tracer_Z_init, only : tracer_Z_init
use MOM_variables, only : surface
use MOM_verticalGrid, only : verticalGrid_type
@@ -369,13 +369,16 @@ subroutine initialize_dye_tracer(restart, day, G, GV, h, diag, OBC, CS, sponge_C
end subroutine initialize_dye_tracer
-subroutine dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS)
+subroutine dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, &
+ evap_CFL_limit, minimum_forcing_depth)
real, dimension(NIMEM_,NJMEM_,NKMEM_), intent(in) :: h_old, h_new, ea, eb
type(forcing), intent(in) :: fluxes
real, intent(in) :: dt
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
- type(dye_tracer_CS), pointer :: CS
+ type(dye_tracer_CS), pointer :: CS
+ real, optional,intent(in) :: evap_CFL_limit
+ real, optional,intent(in) :: minimum_forcing_depth
! This subroutine applies diapycnal diffusion and any other column
! tracer physics or chemistry to the tracers from this file.
! This is a simple example of a set of advected passive tracers.
@@ -398,7 +401,7 @@ subroutine dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS
!
! 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]
-
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
real :: sfc_val ! The surface value for the tracers.
real :: Isecs_per_year ! The number of seconds in a year.
real :: year ! The time in years.
@@ -411,9 +414,20 @@ subroutine dye_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS
if (.not.associated(CS)) return
if (CS%ntr < 1) return
- do m=1,CS%ntr
- call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
- enddo
+ if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
+ do m=1,CS%ntr
+ 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, CS%tr(:,:,:,m) , dt, fluxes, h_work, &
+ evap_CFL_limit, minimum_forcing_depth)
+ call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ else
+ do m=1,CS%ntr
+ call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ endif
do m=1,CS%ntr
do j=G%jsd,G%jed ; do i=G%isd,G%ied
diff --git a/src/tracer/ideal_age_example.F90 b/src/tracer/ideal_age_example.F90
index 04bfa4532d..76d6d6c2f6 100644
--- a/src/tracer/ideal_age_example.F90
+++ b/src/tracer/ideal_age_example.F90
@@ -69,7 +69,7 @@ module ideal_age_example
use MOM_time_manager, only : time_type, get_time
use MOM_tracer_registry, only : register_tracer, tracer_registry_type
use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values
-use MOM_tracer_registry, only : tracer_vertdiff
+use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
use MOM_tracer_Z_init, only : tracer_Z_init
use MOM_variables, only : surface
use MOM_verticalGrid, only : verticalGrid_type
@@ -413,13 +413,16 @@ subroutine initialize_ideal_age_tracer(restart, day, G, GV, h, diag, OBC, CS, &
end subroutine initialize_ideal_age_tracer
-subroutine ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS)
+subroutine ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, &
+ evap_CFL_limit, minimum_forcing_depth)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old, h_new, ea, eb
type(forcing), intent(in) :: fluxes
real, intent(in) :: dt
type(ideal_age_tracer_CS), pointer :: CS
+ real, optional,intent(in) :: evap_CFL_limit
+ real, optional,intent(in) :: minimum_forcing_depth
! This subroutine applies diapycnal diffusion and any other column
! tracer physics or chemistry to the tracers from this file.
! This is a simple example of a set of advected passive tracers.
@@ -442,7 +445,7 @@ subroutine ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G,
!
! 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]
-
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
real :: sfc_val ! The surface value for the tracers.
real :: Isecs_per_year ! The number of seconds in a year.
real :: year ! The time in years.
@@ -453,9 +456,20 @@ subroutine ideal_age_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G,
if (.not.associated(CS)) return
if (CS%ntr < 1) return
- do m=1,CS%ntr
- call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
- enddo
+ if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
+ do m=1,CS%ntr
+ 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, CS%tr(:,:,:,m) , dt, fluxes, h_work, &
+ evap_CFL_limit, minimum_forcing_depth)
+ call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ else
+ do m=1,CS%ntr
+ call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ endif
Isecs_per_year = 1.0 / (365.0*86400.0)
! Set the surface value of tracer 1 to increase exponentially
diff --git a/src/tracer/oil_tracer.F90 b/src/tracer/oil_tracer.F90
index d2725cdb5d..80c5ee66b2 100644
--- a/src/tracer/oil_tracer.F90
+++ b/src/tracer/oil_tracer.F90
@@ -69,7 +69,7 @@ module oil_tracer
use MOM_time_manager, only : time_type, get_time
use MOM_tracer_registry, only : register_tracer, tracer_registry_type
use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values
-use MOM_tracer_registry, only : tracer_vertdiff
+use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
use MOM_tracer_Z_init, only : tracer_Z_init
use MOM_variables, only : surface
use MOM_variables, only : thermo_var_ptrs
@@ -420,7 +420,8 @@ subroutine initialize_oil_tracer(restart, day, G, GV, h, diag, OBC, CS, &
end subroutine initialize_oil_tracer
-subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, tv)
+subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, tv, &
+ evap_CFL_limit, minimum_forcing_depth)
type(ocean_grid_type), intent(in) :: G
type(verticalGrid_type), intent(in) :: GV
real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old, h_new, ea, eb
@@ -428,6 +429,8 @@ subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS
real, intent(in) :: dt
type(oil_tracer_CS), pointer :: CS
type(thermo_var_ptrs), intent(in) :: tv
+ real, optional,intent(in) :: evap_CFL_limit
+ real, optional,intent(in) :: minimum_forcing_depth
! This subroutine applies diapycnal diffusion and any other column
! tracer physics or chemistry to the tracers from this file.
! This is a simple example of a set of advected passive tracers.
@@ -451,6 +454,7 @@ subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS
! 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]
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
real :: Isecs_per_year = 1.0 / (365.0*86400.0)
real :: year, h_total, ldecay
integer :: secs, days
@@ -461,9 +465,20 @@ subroutine oil_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS
if (.not.associated(CS)) return
if (CS%ntr < 1) return
- do m=1,CS%ntr
- call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
- enddo
+ if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
+ do m=1,CS%ntr
+ 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, CS%tr(:,:,:,m) , dt, fluxes, h_work, &
+ evap_CFL_limit, minimum_forcing_depth)
+ call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ else
+ do m=1,CS%ntr
+ call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
+ enddo
+ endif
! Set the surface value of tracer 1 to increase exponentially
! with a 30 year time scale.
diff --git a/src/tracer/pseudo_salt_tracer.F90 b/src/tracer/pseudo_salt_tracer.F90
new file mode 100644
index 0000000000..67526372e3
--- /dev/null
+++ b/src/tracer/pseudo_salt_tracer.F90
@@ -0,0 +1,523 @@
+module pseudo_salt_tracer
+!***********************************************************************
+!* GNU General Public License *
+!* This file is a part of MOM. *
+!* *
+!* MOM is free software; you can redistribute it and/or modify it and *
+!* are expected to follow the terms of the GNU General Public License *
+!* as published by the Free Software Foundation; either version 2 of *
+!* the License, or (at your option) any later version. *
+!* *
+!* MOM is distributed in the hope that it will be useful, but WITHOUT *
+!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
+!* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public *
+!* License for more details. *
+!* *
+!* For the full text of the GNU General Public License, *
+!* write to: Free Software Foundation, Inc., *
+!* 675 Mass Ave, Cambridge, MA 02139, USA. *
+!* or see: http://www.gnu.org/licenses/gpl.html *
+!***********************************************************************
+
+!********+*********+*********+*********+*********+*********+*********+**
+!* *
+!* By Andrew Shao, 2016 *
+!* *
+!* This file contains the routines necessary to model a passive *
+!* tracer that uses the same boundary fluxes as salinity. At the *
+!* beginning of the run, salt is set to the same as tv%S. Any *
+!* deviations between this salt-like tracer and tv%S signifies a *
+!* difference between how active and passive tracers are treated. *
+!* A single subroutine is called from within each file to register *
+!* each of the tracers for reinitialization and advection and to *
+!* register the subroutine that initializes the tracers and set up *
+!* their output and the subroutine that does any tracer physics or *
+!* chemistry along with diapycnal mixing (included here because some *
+!* tracers may float or swim vertically or dye diapycnal processes). *
+!* *
+!* *
+!* Macros written all in capital letters are defined in MOM_memory.h. *
+!* *
+!* A small fragment of the grid is shown below: *
+!* *
+!* j+1 x ^ x ^ x At x: q *
+!* j+1 > o > o > At ^: v *
+!* j x ^ x ^ x At >: u *
+!* j > o > o > At o: h, tr *
+!* j-1 x ^ x ^ x *
+!* i-1 i i+1 At x & ^: *
+!* i i+1 At > & o: *
+!* *
+!* The boundaries always run through q grid points (x). *
+!* *
+!********+*********+*********+*********+*********+*********+*********+**
+
+use MOM_checksums, only : hchksum
+use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
+use MOM_diag_mediator, only : diag_ctrl
+use MOM_diag_to_Z, only : register_Z_tracer, diag_to_Z_CS
+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_grid, only : ocean_grid_type
+use MOM_hor_index, only : hor_index_type
+use MOM_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc
+use MOM_open_boundary, only : ocean_OBC_type
+use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS
+use MOM_sponge, only : set_up_sponge_field, sponge_CS
+use MOM_time_manager, only : time_type, get_time
+use MOM_tracer_registry, only : register_tracer, tracer_registry_type
+use MOM_tracer_registry, only : add_tracer_diagnostics, add_tracer_OBC_values
+use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
+use MOM_tracer_Z_init, only : tracer_Z_init
+use MOM_variables, only : surface
+use MOM_variables, only : thermo_var_ptrs
+use MOM_verticalGrid, only : verticalGrid_type
+use coupler_util, only : set_coupler_values, ind_csurf
+use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux
+
+implicit none ; private
+
+#include
+
+public register_pseudo_salt_tracer, initialize_pseudo_salt_tracer
+public pseudo_salt_tracer_column_physics, pseudo_salt_tracer_surface_state
+public pseudo_salt_stock, pseudo_salt_tracer_end
+
+! NTR_MAX is the maximum number of tracers in this module.
+integer, parameter :: NTR_MAX = 1
+
+type p3d
+ real, dimension(:,:,:), pointer :: p => NULL()
+end type p3d
+
+type, public :: pseudo_salt_tracer_CS ; private
+ integer :: ntr=NTR_MAX ! The number of tracers that are actually used.
+ logical :: coupled_tracers = .false. ! These tracers are not offered to the
+ ! coupler.
+ type(time_type), pointer :: Time ! A pointer to the ocean model's clock.
+ type(tracer_registry_type), pointer :: tr_Reg => NULL()
+ real, pointer :: tr(:,:,:,:) => NULL() ! The array of tracers used in this
+ ! subroutine, in g m-3?
+ real, pointer :: diff(:,:,:,:) => NULL() ! The array of tracers used in this
+ ! subroutine, in g m-3?
+ type(p3d), dimension(NTR_MAX) :: &
+ tr_adx, &! Tracer zonal advective fluxes in g m-3 m3 s-1.An Error Has Occurred
+
+
+ tr_ady, &! Tracer meridional advective fluxes in g m-3 m3 s-1.
+ tr_dfx, &! Tracer zonal diffusive fluxes in g m-3 m3 s-1.
+ tr_dfy ! Tracer meridional diffusive fluxes in g m-3 m3 s-1.
+ logical :: mask_tracers ! If true, pseudo_salt is masked out in massless layers.
+ logical :: pseudo_salt_may_reinit = .true. ! Hard coding since this should not matter
+ integer, dimension(NTR_MAX) :: &
+ ind_tr, & ! Indices returned by aof_set_coupler_flux if it is used and the
+ ! surface tracer concentrations are to be provided to the coupler.
+ id_tracer = -1, id_tr_adx = -1, id_tr_ady = -1, &
+ id_tr_dfx = -1, id_tr_dfy = -1
+ real, dimension(NTR_MAX) :: land_val = -1.0
+
+ type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
+ ! timing of diagnostic output.
+ type(MOM_restart_CS), pointer :: restart_CSp => NULL()
+
+ type(vardesc) :: tr_desc(NTR_MAX)
+end type pseudo_salt_tracer_CS
+
+contains
+
+function register_pseudo_salt_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
+ type(hor_index_type), intent(in) :: HI
+ type(verticalGrid_type), intent(in) :: GV
+ type(param_file_type), intent(in) :: param_file
+ type(pseudo_salt_tracer_CS), pointer :: CS
+ type(tracer_registry_type), pointer :: tr_Reg
+ type(MOM_restart_CS), pointer :: restart_CS
+! This subroutine is used to register tracer fields and subroutines
+! to be used with MOM.
+! Arguments: HI - A horizontal index type structure.
+! (in) GV - The ocean's vertical grid structure.
+! (in) param_file - A structure indicating the open file to parse for
+! model parameter values.
+! (in/out) CS - A pointer that is set to point to the control structure
+! for this module
+! (in/out) tr_Reg - A pointer that is set to point to the control structure
+! for the tracer advection and diffusion module.
+! (in) restart_CS - A pointer to the restart control structure.
+
+! This include declares and sets the variable "version".
+#include "version_variable.h"
+ character(len=40) :: mod = "pseudo_salt_tracer" ! This module's name.
+ character(len=200) :: inputdir ! The directory where the input files are.
+ character(len=48) :: var_name ! The variable's name.
+ character(len=3) :: name_tag ! String for creating identifying pseudo_salt
+ real, pointer :: tr_ptr(:,:,:) => NULL()
+ logical :: register_pseudo_salt_tracer
+ integer :: isd, ied, jsd, jed, nz, m, i, j
+ isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke
+
+ if (associated(CS)) then
+ call MOM_error(WARNING, "register_pseudo_salt_tracer 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, mod, version, "")
+
+ CS%ntr = NTR_MAX
+ allocate(CS%tr(isd:ied,jsd:jed,nz,CS%ntr)) ; CS%tr(:,:,:,:) = 0.0
+ allocate(CS%diff(isd:ied,jsd:jed,nz,CS%ntr)) ; CS%diff(:,:,:,:) = 0.0
+
+ do m=1,CS%ntr
+ ! This is needed to force the compiler not to do a copy in the registration
+ ! calls. Curses on the designers and implementers of Fortran90.
+ CS%tr_desc(m) = var_desc(trim("pseudo_salt_diff"), "kg", &
+ "Difference between pseudo salt passive tracer and salt tracer", caller=mod)
+ tr_ptr => CS%tr(:,:,:,m)
+ call query_vardesc(CS%tr_desc(m), name=var_name, caller="register_pseudo_salt_tracer")
+ ! Register the tracer for the restart file.
+ call register_restart_field(tr_ptr, CS%tr_desc(m), &
+ .not. CS%pseudo_salt_may_reinit, restart_CS)
+ ! Register the tracer for horizontal advection & diffusion.
+ call register_tracer(tr_ptr, CS%tr_desc(m), param_file, HI, GV, tr_Reg, &
+ tr_desc_ptr=CS%tr_desc(m))
+
+ ! Set coupled_tracers to be true (hard-coded above) to provide the surface
+ ! values to the coupler (if any). This is meta-code and its arguments will
+ ! currently (deliberately) give fatal errors if it is used.
+ if (CS%coupled_tracers) &
+ CS%ind_tr(m) = aof_set_coupler_flux(trim(var_name)//'_flux', &
+ flux_type=' ', implementation=' ', caller="register_pseudo_salt_tracer")
+ enddo
+
+ CS%tr_Reg => tr_Reg
+ CS%restart_CSp => restart_CS
+ register_pseudo_salt_tracer = .true.
+
+end function register_pseudo_salt_tracer
+
+subroutine initialize_pseudo_salt_tracer(restart, day, G, GV, h, diag, OBC, CS, &
+ sponge_CSp, diag_to_Z_CSp, tv)
+ logical, intent(in) :: restart
+ type(time_type), target, intent(in) :: day
+ type(ocean_grid_type), intent(in) :: G
+ type(verticalGrid_type), intent(in) :: GV
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h
+ type(diag_ctrl), target, intent(in) :: diag
+ type(ocean_OBC_type), pointer :: OBC
+ type(pseudo_salt_tracer_CS), pointer :: CS
+ type(sponge_CS), pointer :: sponge_CSp
+ type(diag_to_Z_CS), pointer :: diag_to_Z_CSp
+ type(thermo_var_ptrs), intent(in) :: tv
+! This subroutine initializes the CS%ntr tracer fields in tr(:,:,:,:)
+! and it sets up the tracer output.
+
+! Arguments: restart - .true. if the fields have already been read from
+! a restart file.
+! (in) day - Time of the start of the run.
+! (in) G - The ocean's grid structure.
+! (in) GV - The ocean's vertical grid structure.
+! (in) h - Layer thickness, in m or kg m-2.
+! (in) diag - A structure that is used to regulate diagnostic output.
+! (in) OBC - This open boundary condition type specifies whether, where,
+! and what open boundary conditions are used.
+! (in/out) CS - The control structure returned by a previous call to
+! register_pseudo_salt_tracer.
+! (in/out) sponge_CSp - A pointer to the control structure for the sponges, if
+! they are in use. Otherwise this may be unassociated.
+! (in/out) diag_to_Z_Csp - A pointer to the control structure for diagnostics
+! in depth space.
+ character(len=16) :: name ! A variable's name in a NetCDF file.
+ character(len=72) :: longname ! The long name of that variable.
+ character(len=48) :: units ! The dimensions of the variable.
+ character(len=48) :: flux_units ! The units for age tracer fluxes, either
+ ! years m3 s-1 or years kg s-1.
+ logical :: OK
+ integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
+ integer :: IsdB, IedB, JsdB, JedB
+
+ if (.not.associated(CS)) return
+ if (CS%ntr < 1) return
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+ isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
+ IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
+
+ CS%Time => day
+ CS%diag => diag
+ name = "pseudo_salt"
+
+ do m=1,CS%ntr
+ call query_vardesc(CS%tr_desc(m), name=name, caller="initialize_pseudo_salt_tracer")
+ if ((.not.restart) .or. (.not. &
+ query_initialized(CS%tr(:,:,:,m), name, CS%restart_CSp))) then
+ do k=1,nz ; do j=jsd,jed ; do i=isd,ied
+ CS%tr(i,j,k,m) = tv%S(i,j,k)
+ enddo ; enddo ; enddo
+ endif
+ enddo ! Tracer loop
+
+ if (associated(OBC)) then
+ ! All tracers but the first have 0 concentration in their inflows. As this
+ ! is the default value, the following calls are unnecessary.
+ ! do m=1,CS%ntr
+ ! call add_tracer_OBC_values(trim(CS%tr_desc(m)%name), CS%tr_Reg, 0.0)
+ ! enddo
+ endif
+
+ ! This needs to be changed if the units of tracer are changed above.
+ if (GV%Boussinesq) then ; flux_units = "g salt/(m^2 s)"
+ else ; flux_units = "g salt/(m^2 s)" ; endif
+
+ do m=1,CS%ntr
+ ! Register the tracer for the restart file.
+ call query_vardesc(CS%tr_desc(m), name, units=units, longname=longname, &
+ caller="initialize_pseudo_salt_tracer")
+ CS%id_tracer(m) = register_diag_field("ocean_model", trim(name), CS%diag%axesTL, &
+ day, trim(longname) , trim(units))
+ CS%id_tr_adx(m) = register_diag_field("ocean_model", trim(name)//"_adx", &
+ CS%diag%axesCuL, day, trim(longname)//" advective zonal flux" , &
+ trim(flux_units))
+ CS%id_tr_ady(m) = register_diag_field("ocean_model", trim(name)//"_ady", &
+ CS%diag%axesCvL, day, trim(longname)//" advective meridional flux" , &
+ trim(flux_units))
+ CS%id_tr_dfx(m) = register_diag_field("ocean_model", trim(name)//"_dfx", &
+ CS%diag%axesCuL, day, trim(longname)//" diffusive zonal flux" , &
+ trim(flux_units))
+ CS%id_tr_dfy(m) = register_diag_field("ocean_model", trim(name)//"_dfy", &
+ CS%diag%axesCvL, day, trim(longname)//" diffusive zonal flux" , &
+ trim(flux_units))
+ if (CS%id_tr_adx(m) > 0) call safe_alloc_ptr(CS%tr_adx(m)%p,IsdB,IedB,jsd,jed,nz)
+ if (CS%id_tr_ady(m) > 0) call safe_alloc_ptr(CS%tr_ady(m)%p,isd,ied,JsdB,JedB,nz)
+ if (CS%id_tr_dfx(m) > 0) call safe_alloc_ptr(CS%tr_dfx(m)%p,IsdB,IedB,jsd,jed,nz)
+ if (CS%id_tr_dfy(m) > 0) call safe_alloc_ptr(CS%tr_dfy(m)%p,isd,ied,JsdB,JedB,nz)
+
+! Register the tracer for horizontal advection & diffusion.
+ if ((CS%id_tr_adx(m) > 0) .or. (CS%id_tr_ady(m) > 0) .or. &
+ (CS%id_tr_dfx(m) > 0) .or. (CS%id_tr_dfy(m) > 0)) &
+ call add_tracer_diagnostics(name, CS%tr_Reg, CS%tr_adx(m)%p, &
+ CS%tr_ady(m)%p,CS%tr_dfx(m)%p,CS%tr_dfy(m)%p)
+
+ call register_Z_tracer(CS%tr(:,:,:,m), trim(name), longname, units, &
+ day, G, diag_to_Z_CSp)
+ enddo
+
+end subroutine initialize_pseudo_salt_tracer
+
+subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, tv, debug, &
+ evap_CFL_limit, minimum_forcing_depth)
+ type(ocean_grid_type), intent(in) :: G
+ type(verticalGrid_type), intent(in) :: GV
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h_old, h_new, ea, eb
+ type(forcing), intent(in) :: fluxes
+ real, intent(in) :: dt
+ type(pseudo_salt_tracer_CS), pointer :: CS
+ type(thermo_var_ptrs), intent(in) :: tv
+ logical, intent(in) :: debug
+ real, optional,intent(in) :: evap_CFL_limit
+ real, optional,intent(in) :: minimum_forcing_depth
+
+! This subroutine applies diapycnal diffusion and any other column
+! tracer physics or chemistry to the tracers from this file.
+! This is a simple example of a set of advected passive tracers.
+
+! Arguments: h_old - Layer thickness before entrainment, in m or kg m-2.
+! (in) h_new - Layer thickness after entrainment, in m or kg m-2.
+! (in) ea - an array to which the amount of fluid entrained
+! from the layer above during this call will be
+! added, in m or kg m-2.
+! (in) eb - an array to which the amount of fluid entrained
+! from the layer below during this call will be
+! added, in m or kg m-2.
+! (in) fluxes - A structure containing pointers to any possible
+! forcing fields. Unused fields have NULL ptrs.
+! (in) dt - The amount of time covered by this call, in s.
+! (in) G - The ocean's grid structure.
+! (in) GV - The ocean's vertical grid structure.
+! (in) CS - The control structure returned by a previous call to
+! register_pseudo_salt_tracer.
+! (in) tv - Thermodynamic structure with T and S
+! (in) evap_CFL_limit - Limits how much water can be fluxed out of the top layer
+! Stored previously in diabatic CS.
+! (in) minimum_forcing_depth - The smallest depth over which fluxes can be applied
+! Stored previously in diabatic CS.
+! (in) debug - Calculates checksums
+!
+! 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]
+
+ real :: Isecs_per_year = 1.0 / (365.0*86400.0)
+ real :: year, h_total, scale, htot, Ih_limit
+ integer :: secs, days
+ integer :: i, j, k, is, ie, js, je, nz, m, k_max
+ real, allocatable :: local_tr(:,:,:)
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
+ real, dimension(:,:), pointer :: net_salt
+
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+ net_salt=>fluxes%netSalt
+
+ if (.not.associated(CS)) return
+ if (CS%ntr < 1) return
+
+ if (debug) then
+ call hchksum(tv%S,"salt pre pseudo-salt vertdiff", G%HI)
+ call hchksum(CS%tr(:,:,:,1),"pseudo_salt pre pseudo-salt vertdiff", G%HI)
+ endif
+
+ ! This uses applyTracerBoundaryFluxesInOut, usually in ALE mode
+ 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, CS%tr(:,:,:,1), dt, fluxes, h_work, &
+ evap_CFL_limit, minimum_forcing_depth, out_flux_optional=net_salt)
+ call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,1), G, GV)
+ else
+ call tracer_vertdiff(h_work, ea, eb, dt, CS%tr(:,:,:,1), G, GV)
+ endif
+
+ do k=1,nz ; do j=js,je ; do i=is,ie
+ CS%diff(i,j,k,1) = CS%tr(i,j,k,1)-tv%S(i,j,k)
+ enddo ; enddo ; enddo
+
+ if(debug) then
+ call hchksum(tv%S,"salt post pseudo-salt vertdiff", G%HI)
+ call hchksum(CS%tr(:,:,:,1),"pseudo_salt post pseudo-salt vertdiff", G%HI)
+ endif
+
+ allocate(local_tr(G%isd:G%ied,G%jsd:G%jed,nz))
+ do m=1,1
+ if (CS%id_tracer(m)>0) then
+ if (CS%mask_tracers) then
+ do k=1,nz ; do j=js,je ; do i=is,ie
+ if (h_new(i,j,k) < 1.1*GV%Angstrom) then
+ local_tr(i,j,k) = CS%land_val(m)
+ else
+ local_tr(i,j,k) = CS%diff(i,j,k,m)
+ endif
+ enddo ; enddo ; enddo
+ else
+ do k=1,nz ; do j=js,je ; do i=is,ie
+ local_tr(i,j,k) = CS%tr(i,j,k,m)-tv%S(i,j,k)
+ enddo ; enddo ; enddo
+ endif ! CS%mask_tracers
+ call post_data(CS%id_tracer(m),local_tr,CS%diag)
+ endif ! CS%id_tracer(m)>0
+ if (CS%id_tr_adx(m)>0) &
+ call post_data(CS%id_tr_adx(m),CS%tr_adx(m)%p(:,:,:),CS%diag)
+ if (CS%id_tr_ady(m)>0) &
+ call post_data(CS%id_tr_ady(m),CS%tr_ady(m)%p(:,:,:),CS%diag)
+ if (CS%id_tr_dfx(m)>0) &
+ call post_data(CS%id_tr_dfx(m),CS%tr_dfx(m)%p(:,:,:),CS%diag)
+ if (CS%id_tr_dfy(m)>0) &
+ call post_data(CS%id_tr_dfy(m),CS%tr_dfy(m)%p(:,:,:),CS%diag)
+ enddo
+ deallocate(local_tr)
+
+end subroutine pseudo_salt_tracer_column_physics
+
+function pseudo_salt_stock(h, stocks, G, GV, CS, names, units, stock_index)
+ type(ocean_grid_type), intent(in) :: G
+ type(verticalGrid_type), intent(in) :: GV
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h
+ real, dimension(:), intent(out) :: stocks
+ type(pseudo_salt_tracer_CS), pointer :: CS
+ character(len=*), dimension(:), intent(out) :: names
+ character(len=*), dimension(:), intent(out) :: units
+ integer, optional, intent(in) :: stock_index
+ integer :: pseudo_salt_stock
+! This function 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.
+
+! Arguments: h - Layer thickness, in m or kg m-2.
+! (out) stocks - the mass-weighted integrated amount of each tracer,
+! in kg times concentration units.
+! (in) G - The ocean's grid structure.
+! (in) GV - The ocean's vertical grid structure.
+! (in) CS - The control structure returned by a previous call to
+! register_pseudo_salt_tracer.
+! (out) names - the names of the stocks calculated.
+! (out) units - the units of the stocks calculated.
+! (in,opt) stock_index - the coded index of a specific stock being sought.
+! Return value: the number of stocks calculated here.
+
+ integer :: i, j, k, is, ie, js, je, nz, m
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
+
+ pseudo_salt_stock = 0
+ if (.not.associated(CS)) return
+ if (CS%ntr < 1) 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
+
+ do m=1,1
+ call query_vardesc(CS%tr_desc(m), name=names(m), units=units(m), caller="pseudo_salt_stock")
+ units(m) = trim(units(m))//" kg"
+ stocks(m) = 0.0
+ do k=1,nz ; do j=js,je ; do i=is,ie
+ stocks(m) = stocks(m) + CS%diff(i,j,k,m) * &
+ (G%mask2dT(i,j) * G%areaT(i,j) * h(i,j,k))
+ enddo ; enddo ; enddo
+ stocks(m) = GV%H_to_kg_m2 * stocks(m)
+ enddo
+
+ pseudo_salt_stock = CS%ntr
+
+end function pseudo_salt_stock
+
+subroutine pseudo_salt_tracer_surface_state(state, h, G, CS)
+ type(ocean_grid_type), intent(in) :: G
+ type(surface), intent(inout) :: state
+ real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h
+ type(pseudo_salt_tracer_CS), pointer :: CS
+! This particular tracer package does not report anything back to the coupler.
+! The code that is here is just a rough guide for packages that would.
+! Arguments: state - A structure containing fields that describe the
+! surface state of the ocean.
+! (in) h - Layer thickness, in m or kg m-2.
+! (in) G - The ocean's grid structure.
+! (in) CS - The control structure returned by a previous call to
+! register_pseudo_salt_tracer.
+ integer :: m, is, ie, js, je
+ is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
+
+ if (.not.associated(CS)) return
+
+ if (CS%coupled_tracers) then
+ do m=1,CS%ntr
+ ! This call loads the surface vlues into the appropriate array in the
+ ! coupler-type structure.
+ call set_coupler_values(CS%tr(:,:,1,m), state%tr_fields, CS%ind_tr(m), &
+ ind_csurf, is, ie, js, je)
+ enddo
+ endif
+
+end subroutine pseudo_salt_tracer_surface_state
+
+subroutine pseudo_salt_tracer_end(CS)
+ type(pseudo_salt_tracer_CS), pointer :: CS
+ integer :: m
+
+ if (associated(CS)) then
+ if (associated(CS%tr)) deallocate(CS%tr)
+ if (associated(CS%tr)) deallocate(CS%diff)
+ do m=1,CS%ntr
+ if (associated(CS%tr_adx(m)%p)) deallocate(CS%tr_adx(m)%p)
+ if (associated(CS%tr_ady(m)%p)) deallocate(CS%tr_ady(m)%p)
+ if (associated(CS%tr_dfx(m)%p)) deallocate(CS%tr_dfx(m)%p)
+ if (associated(CS%tr_dfy(m)%p)) deallocate(CS%tr_dfy(m)%p)
+ enddo
+
+ deallocate(CS)
+ endif
+end subroutine pseudo_salt_tracer_end
+
+end module pseudo_salt_tracer