Skip to content

Commit

Permalink
Pseudo salt and salt have very small differences. Suspect difference …
Browse files Browse the repository at this point in the history 
…because of freshwater flux during ice formation/melting
  • Loading branch information
Andrew Shao committed Sep 13, 2016
1 parent 36acb8a commit 982c373
Show file tree
Hide file tree
Showing 4 changed files with 70 additions and 29 deletions.
42 changes: 30 additions & 12 deletions src/parameterizations/vertical/MOM_diabatic_driver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -137,6 +137,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
type(diag_ctrl), pointer :: diag !< structure used to regulate timing of diagnostic output
real :: MLDdensityDifference !< Density difference used to determine MLD_user
integer :: nsw !< SW_NBANDS
Expand Down Expand Up @@ -920,8 +922,13 @@ 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
print *, "Using tracer_vertdiff for T/S"
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)

Expand Down Expand Up @@ -1000,7 +1007,13 @@ 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
print *, "Using tracer_vertdiff for T/S"
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
Expand Down Expand Up @@ -1121,13 +1134,13 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
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, eatr, ebtr, fluxes, dt, G, GV, tv, &
CS%optics, CS%tracer_flow_CSp, &
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%optics, CS%tracer_flow_CSp, CS%debug)
endif

elseif (associated(visc%Kd_extra_S)) then ! extra diffusivity for passive tracers
Expand All @@ -1151,25 +1164,25 @@ subroutine diabatic(u, v, h, tv, fluxes, visc, ADp, CDp, dt, G, GV, CS)
enddo ; enddo ; enddo
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, &
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%optics, CS%tracer_flow_CSp, CS%debug)
endif

else
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, ea, eb, fluxes, dt, G, GV, tv, &
CS%optics, CS%tracer_flow_CSp, &
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%optics, CS%tracer_flow_CSp, CS%debug)
endif

endif ! (CS%mix_boundary_tracers)
Expand Down Expand Up @@ -1457,7 +1470,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

Expand Down Expand Up @@ -1887,6 +1900,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"
Expand Down
23 changes: 13 additions & 10 deletions src/tracer/MOM_tracer_diabatic.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
call applyTracerBoundaryFluxesInOut(G, GV, tr, dt, sfc_src, fluxes, h_work, &
evap_CFL_limit, minimum_forcing_depth)
! Zero out surface fluxes because they were applied in the call to applyTracerBoundaryFluxesInOut
do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; enddo ; enddo
do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; enddo ; enddo
endif

if (present(sink_rate)) then
Expand Down Expand Up @@ -156,7 +156,7 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
d1(i) = b_denom_1 * b1(i)
h_tr = h_work(i,j,1) + h_neglect
tr(i,j,1) = b1(i)*(h_tr*tr(i,j,1) + sfc_src(i,j))
tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j)
endif ; enddo
do k=2,nz-1 ; do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
c1(i,k) = eb(i,j,k-1) * b1(i)
Expand Down Expand Up @@ -189,30 +189,31 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
else
!$OMP do
do j=js,je
do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
do i=is,ie ; if (G%mask2dT(i,j) > -0.5) then
h_tr = h_work(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))
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
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_work(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
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_work(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)) + &
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
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
Expand Down Expand Up @@ -292,6 +293,8 @@ subroutine applyTracerBoundaryFluxesInOut(G, GV, Tr, dt, sfc_src, fluxes, h, &
! 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)
Expand Down
14 changes: 10 additions & 4 deletions src/tracer/MOM_tracer_flow_control.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -355,7 +355,7 @@ 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, &
evap_CFL_limit, minimum_forcing_depth)
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
Expand All @@ -364,9 +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.

Expand All @@ -388,6 +389,11 @@ 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.")
Expand Down Expand Up @@ -443,7 +449,7 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, dt, G, GV, tv, o
#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, &
G, GV, CS%pseudo_salt_tracer_CSp, tv, debug,&
evap_CFL_limit=evap_CFL_limit, &
minimum_forcing_depth=minimum_forcing_depth)

Expand Down Expand Up @@ -480,7 +486,7 @@ subroutine call_tracer_column_fns(h_old, h_new, ea, eb, fluxes, dt, G, GV, tv, o
#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)
G, GV, CS%pseudo_salt_tracer_CSp, tv, debug)


endif
Expand Down
20 changes: 17 additions & 3 deletions src/tracer/pseudo_salt_tracer.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,7 @@ module pseudo_salt_tracer
!* *
!********+*********+*********+*********+*********+*********+*********+**

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
Expand Down Expand Up @@ -247,7 +248,7 @@ subroutine initialize_pseudo_salt_tracer(restart, day, G, GV, h, diag, OBC, CS,
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=js,je ; do i=is,ie
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
Expand Down Expand Up @@ -300,7 +301,7 @@ subroutine initialize_pseudo_salt_tracer(restart, day, G, GV, h, diag, OBC, CS,

end subroutine initialize_pseudo_salt_tracer

subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, CS, tv, &
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
Expand All @@ -309,8 +310,10 @@ subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G
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.
Expand All @@ -330,6 +333,12 @@ subroutine pseudo_salt_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_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]
Expand Down Expand Up @@ -360,8 +369,13 @@ subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G
do m=1,CS%ntr
call tracer_vertdiff(h_old, ea, eb, dt, CS%tr(:,:,:,m), G, GV)
enddo
endif
endif

! if(debug) then
call hchksum(tv%S,"S post pseudo-salt vertdiff", G%HI)
call hchksum(CS%tr(:,:,:,m),"pseudo_salt post pseudo-salt vertdiff", G%HI)
! endif

allocate(local_tr(G%isd:G%ied,G%jsd:G%jed,nz))
do m=1,CS%ntr
if (CS%id_tracer(m)>0) then
Expand Down

0 comments on commit 982c373

Please sign in to comment.