From c266f9f63c51cb4769506e06104c0b82f6774cc6 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 30 May 2018 14:56:10 -0600 Subject: [PATCH 1/7] First step towards using Kd_heat and Kd_salt --- .../vertical/MOM_diabatic_driver.F90 | 90 ++++++++++--------- 1 file changed, 50 insertions(+), 40 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index f283f6243a..843ecf76cd 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -281,9 +281,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & - ea, & ! amount of fluid entrained from the layer above within + ea_s, & ! amount of fluid entrained from the layer above within ! one time step (m for Bouss, kg/m^2 for non-Bouss) - eb, & ! amount of fluid entrained from the layer below within + eb_s, & ! amount of fluid entrained from the layer below within + ! one time step (m for Bouss, kg/m^2 for non-Bouss) + ea_t, & ! amount of fluid entrained from the layer above within + ! one time step (m for Bouss, kg/m^2 for non-Bouss) + eb_t, & ! amount of fluid entrained from the layer below within ! 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) @@ -555,12 +559,34 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call cpu_clock_end(id_clock_set_diffusivity) if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)") + ! Set diffusivities for heat and salt + +!$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat) +!$OMP do + do k=1,nz+1 ; do j=js,je ; do i=is,ie + Kd_salt(i,j,k) = Kd_int(i,j,k) + Kd_heat(i,j,k) = Kd_int(i,j,k) + enddo ; enddo ; enddo + if (associated(visc%Kd_extra_S)) then +!$OMP do + do k=1,nz+1 ; do j=js,je ; do i=is,ie + Kd_salt(i,j,k) = Kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k) + enddo ; enddo ; enddo + endif + if (associated(visc%Kd_extra_T)) then +!$OMP do + do k=1,nz+1 ; do j=js,je ; do i=is,ie + Kd_heat(i,j,k) = Kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k) + enddo ; enddo ; enddo + endif +!$OMP end parallel + if (CS%debug) then call MOM_state_chksum("after set_diffusivity ", u, v, h, G, GV, haloshift=0) call MOM_forcing_chksum("after set_diffusivity ", fluxes, G, haloshift=0) call MOM_thermovar_chksum("after set_diffusivity ", tv, G) - call hchksum(Kd, "after set_diffusivity Kd",G%HI,haloshift=0) - call hchksum(Kd_Int, "after set_diffusivity Kd_Int",G%HI,haloshift=0) + call hchksum(Kd_heat, "after set_diffusivity Kd_heat",G%HI,haloshift=0) + call hchksum(Kd_salt, "after set_diffusivity Kd_salt",G%HI,haloshift=0) endif if (CS%useKPP) then @@ -577,26 +603,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! since the matching to nonzero interior diffusivity can be problematic. ! Changes: Kd_int. Sets: KPP_NLTheat, KPP_NLTscalar -!$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat) -!$OMP do - do k=1,nz+1 ; do j=js,je ; do i=is,ie - Kd_salt(i,j,k) = Kd_int(i,j,k) - Kd_heat(i,j,k) = Kd_int(i,j,k) - enddo ; enddo ; enddo - if (associated(visc%Kd_extra_S)) then -!$OMP do - do k=1,nz+1 ; do j=js,je ; do i=is,ie - Kd_salt(i,j,k) = Kd_salt(i,j,k) + visc%Kd_extra_S(i,j,k) - enddo ; enddo ; enddo - endif - if (associated(visc%Kd_extra_T)) then -!$OMP do - do k=1,nz+1 ; do j=js,je ; do i=is,ie - Kd_heat(i,j,k) = Kd_heat(i,j,k) + visc%Kd_extra_T(i,j,k) - enddo ; enddo ; enddo - endif -!$OMP end parallel - call KPP_compute_BLD(CS%KPP_CSp, G, GV, h, tv%T, tv%S, u, v, tv%eqn_of_state, & fluxes%ustar, CS%KPP_buoy_flux) @@ -610,6 +616,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call pass_var(Hml, G%domain, halo=1) endif +! GMM, I believe the following can be deleted??? if (.not. CS%KPPisPassive) then !$OMP do do k=1,nz+1 ; do j=js,je ; do i=is,ie @@ -628,34 +635,24 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo ; enddo endif endif ! not passive + !$OMP end parallel +!!!!!!! delete above ?? !!!!!!!!!!!!! + call cpu_clock_end(id_clock_kpp) if (showCallTree) call callTree_waypoint("done with KPP_calculate (diabatic)") if (CS%debug) then call MOM_state_chksum("after KPP", u, v, h, G, GV, haloshift=0) call MOM_forcing_chksum("after KPP", fluxes, G, haloshift=0) call MOM_thermovar_chksum("after KPP", tv, G) - call hchksum(Kd, "after KPP Kd",G%HI,haloshift=0) - call hchksum(Kd_Int, "after KPP Kd_Int",G%HI,haloshift=0) + call hchksum(Kd_heat, "after KPP Kd_heat",G%HI,haloshift=0) + call hchksum(Kd_salt, "after KPP Kd_salt",G%HI,haloshift=0) endif endif ! endif for KPP - ! Add vertical diff./visc. due to convection (computed via CVMix) - if (CS%use_CVMix_conv) then - call calculate_CVMix_conv(h, tv, G, GV, CS%CVMix_conv_csp, Hml) - - !!!!!!!! GMM, the following needs to be checked !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - do k=1,nz ; do j=js,je ; do i=is,ie - Kd_int(i,j,k) = Kd_int(i,j,k) + CS%CVMix_conv_csp%kd_conv(i,j,k) - visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + CS%CVMix_conv_csp%kv_conv(i,j,k) - enddo ; enddo ; enddo - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - endif if (CS%useKPP) then - call cpu_clock_begin(id_clock_kpp) if (CS%debug) then call hchksum(CS%KPP_temp_flux, "before KPP_applyNLT netHeat",G%HI,haloshift=0, scale=GV%H_to_m) @@ -676,7 +673,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call MOM_forcing_chksum("after KPP_applyNLT ", fluxes, G, haloshift=0) call MOM_thermovar_chksum("after KPP_applyNLT ", tv, G) endif - endif ! endif for KPP ! Differential diffusion done here. @@ -703,6 +699,20 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif + ! Add vertical diff./visc. due to convection (computed via CVMix) + if (CS%use_CVMix_conv) then + call calculate_CVMix_conv(h, tv, G, GV, CS%CVMix_conv_csp, Hml) + +!$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,visc,CS,Kd_heat) +!$OMP do + do k=1,nz+1 ; do j=js,je ; do i=is,ie + Kd_heat(i,j,k) = Kd_heat(i,j,k) + CS%CVMix_conv_csp%kd_conv(i,j,k) + Kd_salt(i,j,k) = Kd_salt(i,j,k) + CS%CVMix_conv_csp%kd_conv(i,j,k) + visc%Kv_slow(i,j,k) = visc%Kv_slow(i,j,k) + CS%CVMix_conv_csp%kv_conv(i,j,k) + enddo ; enddo ; enddo +!$OMP end parallel + endif + ! This block sets ea, eb from Kd or Kd_int. ! set ea=eb=Kd_int on interfaces for use in the tri-diagonal solver. From 88886d9b66cd7b7e9b7a54778bbd31d1d4b665d6 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 30 May 2018 15:27:22 -0600 Subject: [PATCH 2/7] Doxygenize tracer_vertdiff --- src/tracer/MOM_tracer_diabatic.F90 | 45 +++++++++++++++--------------- 1 file changed, 23 insertions(+), 22 deletions(-) diff --git a/src/tracer/MOM_tracer_diabatic.F90 b/src/tracer/MOM_tracer_diabatic.F90 index 0bdd327033..c8ce7700db 100644 --- a/src/tracer/MOM_tracer_diabatic.F90 +++ b/src/tracer/MOM_tracer_diabatic.F90 @@ -15,13 +15,13 @@ module MOM_tracer_diabatic #include public tracer_vertdiff public applyTracerBoundaryFluxesInOut + +contains + !> 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 @@ -43,27 +43,28 @@ subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, & 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. + ! local variables + 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. + 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. + 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. From e15fe54cb87681d2902e41fb34df54ca9c220772 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 31 May 2018 09:04:51 -0600 Subject: [PATCH 3/7] Add missing code relate to old double-diffusion method --- src/parameterizations/vertical/MOM_set_viscosity.F90 | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index ec1b09a5ad..33a7fbaa4f 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1834,7 +1834,7 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) real :: omega_frac_dflt integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, i, j, n logical :: use_kappa_shear, adiabatic, use_omega - logical :: use_CVMix_ddiff + logical :: use_CVMix_ddiff, differential_diffusion type(OBC_segment_type), pointer :: segment ! pointer to OBC segment type ! This include declares and sets the variable "version". #include "version_variable.h" @@ -1859,6 +1859,7 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) call log_version(param_file, mdl, version, "") CS%RiNo_mix = .false. ; use_CVMix_ddiff = .false. use_kappa_shear = .false. !; adiabatic = .false. ! Needed? -AJA + differential_diffusion = .false. call get_param(param_file, mdl, "BOTTOMDRAGLAW", CS%bottomdraglaw, & "If true, the bottom stress is calculated with a drag \n"//& "law of the form c_drag*|u|*u. The velocity magnitude \n"//& @@ -1885,6 +1886,10 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) if (.not.adiabatic) then use_kappa_shear = kappa_shear_is_used(param_file) CS%RiNo_mix = use_kappa_shear + call get_param(param_file, mdl, "DOUBLE_DIFFUSION", differential_diffusion, & + "If true, increase diffusivitives for temperature or salt \n"//& + "based on double-diffusive paramaterization from MOM4/KPP.", & + default=.false.) use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) endif @@ -2038,7 +2043,7 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) Time, 'Rayleigh drag velocity at v points', 'm s-1') endif - if (use_CVMix_ddiff) then + if (use_CVMix_ddiff .or. differential_diffusion) then allocate(visc%Kd_extra_T(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_T = 0.0 allocate(visc%Kd_extra_S(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_S = 0.0 endif From fcdd55deebe96018ec3c33c590ede0db3024597e Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 31 May 2018 09:05:09 -0600 Subject: [PATCH 4/7] Add missing code relate to old double-diffusion method --- .../vertical/MOM_set_diffusivity.F90 | 30 +++++++------------ 1 file changed, 10 insertions(+), 20 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 93018c9dac..fb6f25226d 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -148,22 +148,11 @@ module MOM_set_diffusivity type(int_tide_CS), pointer :: int_tide_CSp => NULL() type(tidal_mixing_cs), pointer :: tm_csp => NULL() - integer :: id_maxTKE = -1 - integer :: id_TKE_to_Kd = -1 - - integer :: id_Kd_user = -1 - integer :: id_Kd_layer = -1 - integer :: id_Kd_BBL = -1 - integer :: id_Kd_BBL_z = -1 - integer :: id_Kd_user_z = -1 - integer :: id_Kd_Work = -1 - - integer :: id_N2 = -1 - integer :: id_N2_z = -1 - integer :: id_KT_extra = -1 - integer :: id_KS_extra = -1 - integer :: id_KT_extra_z = -1 - integer :: id_KS_extra_z = -1 + integer :: id_maxTKE = -1, id_TKE_to_Kd = -1, id_Kd_user = -1 + integer :: id_Kd_layer = -1, id_Kd_BBL = -1, id_Kd_BBL_z = -1 + integer :: id_Kd_user_z = -1, id_N2 = -1, id_N2_z = -1 + integer :: id_Kd_Work = -1, id_KT_extra = -1, id_KS_extra = -1 + integer :: id_KT_extra_z = -1, id_KS_extra_z = -1 end type set_diffusivity_CS @@ -284,9 +273,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & use_EOS = associated(tv%eqn_of_state) - if ((CS%use_CVMix_ddiff) .or. CS%double_diffusion .and. & - .not.(associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S)) ) & - call MOM_error(FATAL, "set_diffusivity: visc%Kd_extra_T and "//& + if ((CS%use_CVMix_ddiff .or. CS%double_diffusion) .and. .not. & + (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S))) & + call MOM_error(FATAL, "set_diffusivity: both visc%Kd_extra_T and "//& "visc%Kd_extra_S must be associated when USE_CVMIX_DDIFF or DOUBLE_DIFFUSION are true.") ! Set Kd, Kd_int and Kv_slow to constant values. @@ -2106,6 +2095,7 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp "If true, increase diffusivitives for temperature or salt \n"//& "based on double-diffusive paramaterization from MOM4/KPP.", & default=.false.) + if (CS%double_diffusion) then call get_param(param_file, mdl, "MAX_RRHO_SALT_FINGERS", CS%Max_Rrho_salt_fingers, & "Maximum density ratio for salt fingering regime.", & @@ -2137,7 +2127,7 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp "Bottom Boundary Layer Diffusivity", z_grid='z') CS%id_Kd_BBL_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) endif - endif + endif ! old double-diffusion if (CS%user_change_diff) then call user_change_diff_init(Time, G, param_file, diag, CS%user_change_diff_CSp) From 4c36f8f9faec86a43278b09e9e2f31fdc78d7a26 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 31 May 2018 11:26:58 -0600 Subject: [PATCH 5/7] Computes diffusivity for salt and heat separetely This commit adds the capability to compute vertical diffusivities for salt and heat separetely (via tracer_vertdiff). --- .../vertical/MOM_diabatic_driver.F90 | 183 ++++++++++-------- 1 file changed, 106 insertions(+), 77 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index c391ceb1c2..962594ae00 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -172,8 +172,9 @@ module MOM_diabatic_driver integer :: id_cg1 = -1 ! diag handle for mode-1 speed (BDM) integer, allocatable, dimension(:) :: id_cn ! diag handle for all mode speeds (BDM) - integer :: id_dudt_dia = -1, id_dvdt_dia = -1, id_wd = -1 - integer :: id_ea = -1, id_eb = -1, id_Kd_z = -1 + integer :: id_wd = -1, id_ea = -1, id_eb = -1 ! used by layer diabatic + integer :: id_dudt_dia = -1, id_dvdt_dia = -1, id_ea_s = -1, id_eb_s = -1 + integer :: id_ea_t = -1, id_eb_t = -1, id_Kd_z = -1 integer :: id_Kd_heat = -1, id_Kd_salt = -1, id_Kd_interface = -1, id_Kd_ePBL = -1 integer :: id_Tdif_z = -1, id_Tadv_z = -1, id_Sdif_z = -1, id_Sadv_z = -1 integer :: id_Tdif = -1, id_Tadv = -1, id_Sdif = -1, id_Sadv = -1 @@ -675,7 +676,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif endif ! endif for KPP - ! Differential diffusion done here. + ! GMM, this is the "old" method for applying differential diffusion. + ! TODO\ the following will not work with KPP. We need to add a FATAL + ! error to avoid that. + ! Changes: tv%T, tv%S ! If using matching within the KPP scheme, then this step needs to provide ! a diffusivity and happen before KPP. But generally in MOM, we do not match @@ -693,18 +697,21 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! increment heat and salt diffusivity. ! CS%useKPP==.true. already has extra_T and extra_S included if (.not. CS%useKPP) then +!$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,visc,Kd_heat) +!$OMP do do K=2,nz ; do j=js,je ; do i=is,ie Kd_heat(i,j,K) = Kd_heat(i,j,K) + visc%Kd_extra_T(i,j,K) Kd_salt(i,j,K) = Kd_salt(i,j,K) + visc%Kd_extra_S(i,j,K) enddo ; enddo ; enddo +!$OMP end parallel endif endif - ! Add vertical diff./visc. due to convection (computed via CVMix) + ! Calculate vertical mixing due to convection (computed via CVMix) if (CS%use_CVMix_conv) then call calculate_CVMix_conv(h, tv, G, GV, CS%CVMix_conv_csp, Hml) - + ! Increment vertical diffusion and viscosity due to convection !$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,visc,CS,Kd_heat) !$OMP do do k=1,nz+1 ; do j=js,je ; do i=is,ie @@ -716,23 +723,26 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif - ! This block sets ea, eb from Kd or Kd_int. - ! set ea=eb=Kd_int on interfaces for use in the tri-diagonal solver. + ! set ea_t=eb_t=Kd_heat and ea_s=eb_s=Kd_salt on interfaces for use in the + ! tri-diagonal solver. do j=js,je ; do i=is,ie - ea(i,j,1) = 0. + ea_t(i,j,1) = 0.; ea_s(i,j,1) = 0. enddo ; enddo -!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_neglect,h,ea,GV,dt,Kd_int,eb) & + +!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_neglect,h,ea_t,ea_s,GV,dt,Kd_salt,Kd_heat,eb_t,eb_s) & !$OMP private(hval) do k=2,nz ; do j=js,je ; do i=is,ie hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) - ea(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_int(i,j,k) - eb(i,j,k-1) = ea(i,j,k) + ea_t(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_heat(i,j,k) + eb_t(i,j,k-1) = ea_t(i,j,k) + ea_s(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_salt(i,j,k) + eb_s(i,j,k-1) = ea_s(i,j,k) enddo ; enddo ; enddo do j=js,je ; do i=is,ie - eb(i,j,nz) = 0. + eb_t(i,j,nz) = 0.; eb_s(i,j,nz) = 0. enddo ; enddo - if (showCallTree) call callTree_waypoint("done setting ea,eb from Kd_int (diabatic)") + if (showCallTree) call callTree_waypoint("done setting ea_t,ea_s,eb_t,eb_s from Kd_heat and Kd_salt (diabatic)") ! Save fields before boundary forcing is applied for tendency diagnostics if (CS%boundary_forcing_tendency_diag) then @@ -758,8 +768,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) if (CS%debug) then - call hchksum(ea, "after applyBoundaryFluxes ea",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after applyBoundaryFluxes eb",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(ea_t, "after applyBoundaryFluxes ea_t",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb_t, "after applyBoundaryFluxes eb_t",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(ea_s, "after applyBoundaryFluxes ea_s",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb_s, "after applyBoundaryFluxes eb_s",G%HI,haloshift=0, scale=GV%H_to_m) call hchksum(cTKE, "after applyBoundaryFluxes cTKE",G%HI,haloshift=0) call hchksum(dSV_dT, "after applyBoundaryFluxes dSV_dT",G%HI,haloshift=0) call hchksum(dSV_dS, "after applyBoundaryFluxes dSV_dS",G%HI,haloshift=0) @@ -787,8 +799,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif Ent_int = Kd_add_here * (GV%m_to_H**2 * dt) / & (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect) - eb(i,j,k-1) = eb(i,j,k-1) + Ent_int - ea(i,j,k) = ea(i,j,k) + Ent_int + eb_t(i,j,k-1) = eb_t(i,j,k-1) + Ent_int + ea_t(i,j,k) = ea_t(i,j,k) + Ent_int + eb_s(i,j,k-1) = eb_s(i,j,k-1) + Ent_int + ea_s(i,j,k) = ea_s(i,j,k) + Ent_int Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here ! for diagnostics @@ -798,8 +812,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo ; enddo if (CS%debug) then - call hchksum(ea, "after ePBL ea",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after ePBL eb",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(ea_t, "after ePBL ea_t",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb_t, "after ePBL eb_t",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(ea_s, "after ePBL ea_s",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb_s, "after ePBL eb_s",G%HI,haloshift=0, scale=GV%H_to_m) call hchksum(Kd_ePBL, "after ePBL Kd_ePBL",G%HI,haloshift=0) endif @@ -844,9 +860,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & do j=js,je do i=is,ie hold(i,j,1) = h(i,j,1) - h(i,j,1) = h(i,j,1) + (eb(i,j,1) - ea(i,j,2)) + h(i,j,1) = h(i,j,1) + (eb_t(i,j,1) - ea_t(i,j,2)) hold(i,j,nz) = h(i,j,nz) - h(i,j,nz) = h(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1)) + h(i,j,nz) = h(i,j,nz) + (ea_t(i,j,nz) - eb_t(i,j,nz-1)) if (h(i,j,1) <= 0.0) then h(i,j,1) = GV%Angstrom endif @@ -856,8 +872,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo do k=2,nz-1 ; do i=is,ie hold(i,j,k) = h(i,j,k) - h(i,j,k) = h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + & - (eb(i,j,k) - ea(i,j,k+1))) + h(i,j,k) = h(i,j,k) + ((ea_t(i,j,k) - eb_t(i,j,k-1)) + & + (eb_t(i,j,k) - ea_t(i,j,k+1))) if (h(i,j,k) <= 0.0) then h(i,j,k) = GV%Angstrom endif @@ -871,6 +887,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call MOM_forcing_chksum("after negative check ", fluxes, G, haloshift=0) call MOM_thermovar_chksum("after negative check ", tv, G) endif + if (showCallTree) call callTree_waypoint("done with h=ea-eb (diabatic)") if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G) @@ -879,11 +896,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (associated(tv%T)) then if (CS%debug) then - call hchksum(ea, "before triDiagTS ea ",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "before triDiagTS eb ",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(ea_t, "before triDiagTS ea_t ",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb_t, "before triDiagTS eb_t ",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(ea_s, "before triDiagTS ea_s ",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(eb_s, "before triDiagTS eb_s ",G%HI,haloshift=0, scale=GV%H_to_m) endif - call cpu_clock_begin(id_clock_tridiag) + call cpu_clock_begin(id_clock_tridiag) ! Keep salinity from falling below a small but positive threshold. ! This constraint is needed for SIS1 ice model, which can extract ! more salt than is present in the ocean. SIS2 does not suffer @@ -901,13 +920,20 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif ! Changes T and S via the tridiagonal solver; no change to h - 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 tracer_vertdiff(hold, ea_t, eb_t, dt, tv%T, G, GV) + call tracer_vertdiff(hold, ea_s, eb_s, dt, tv%S, G, GV) - call triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S) - endif + ! GMM, with the new approach of having ea,eb for temp and salt + ! only tracer_vertdiff can be used at this time. Therefore, + ! I am commenting the following code. Should this be deleted? + + !if (CS%tracer_tridiag) then + ! call tracer_vertdiff(hold, ea_t, eb_t, dt, tv%T, G, GV) + ! call tracer_vertdiff(hold, ea_s, eb_s, dt, tv%S, G, GV) + !else + ! + ! call triDiagTS(G, GV, is, ie, js, je, hold, ea_t, eb_t, tv%T, tv%S) + !endif ! diagnose temperature, salinity, heat, and salt tendencies ! Note: hold here refers to the thicknesses from before the dual-entraintment when using @@ -917,8 +943,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, G, GV, CS) if (CS%id_diabatic_diff_h > 0) call post_data(CS%id_diabatic_diff_h, hold, CS%diag, alt_h = hold) endif - call cpu_clock_end(id_clock_tridiag) + if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)") endif ! endif corresponding to if (associated(tv%T)) @@ -928,8 +954,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%debug) then call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, haloshift=0) call MOM_thermovar_chksum("after mixed layer ", tv, G) - call hchksum(ea, "after mixed layer ea", G%HI, scale=GV%H_to_m) - call hchksum(eb, "after mixed layer eb", G%HI, scale=GV%H_to_m) endif ! Whenever thickness changes let the diag manager know, as the @@ -945,9 +969,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo !$OMP parallel do default(shared) do K=2,nz ; do j=js,je ; do i=is,ie - Tdif_flx(i,j,K) = (Idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * & + Tdif_flx(i,j,K) = (Idt * 0.5*(ea_t(i,j,k) + eb_t(i,j,k-1))) * & (tv%T(i,j,k-1) - tv%T(i,j,k)) - Tadv_flx(i,j,K) = (Idt * (ea(i,j,k) - eb(i,j,k-1))) * & + Tadv_flx(i,j,K) = (Idt * (ea_t(i,j,k) - eb_t(i,j,k-1))) * & 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k)) enddo ; enddo ; enddo endif @@ -959,9 +983,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo !$OMP parallel do default(shared) do K=2,nz ; do j=js,je ; do i=is,ie - Sdif_flx(i,j,K) = (Idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * & + Sdif_flx(i,j,K) = (Idt * 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))) * & (tv%S(i,j,k-1) - tv%S(i,j,k)) - Sadv_flx(i,j,K) = (Idt * (ea(i,j,k) - eb(i,j,k-1))) * & + Sadv_flx(i,j,K) = (Idt * (ea_s(i,j,k) - eb_s(i,j,k-1))) * & 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k)) enddo ; enddo ; enddo endif @@ -974,7 +998,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & !$OMP parallel do default(shared) private(htot,in_boundary,add_ent) do j=js,je do i=is,ie - ebtr(i,j,nz) = eb(i,j,nz) + ebtr(i,j,nz) = eb_s(i,j,nz) htot(i) = 0.0 in_boundary(i) = (G%mask2dT(i,j) > 0.0) enddo @@ -992,19 +1016,20 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & add_ent = ((dt * CS%Kd_min_tr) * GV%m_to_H**2) * & ((h(i,j,k-1)+h(i,j,k)+h_neglect) / & (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - & - 0.5*(ea(i,j,k) + eb(i,j,k-1)) + 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1)) if (htot(i) < Tr_ea_BBL) then add_ent = max(0.0, add_ent, & - (Tr_ea_BBL - htot(i)) - min(ea(i,j,k),eb(i,j,k-1))) + (Tr_ea_BBL - htot(i)) - min(ea_s(i,j,k),eb_s(i,j,k-1))) elseif (add_ent < 0.0) then add_ent = 0.0 ; in_boundary(i) = .false. endif - ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent - eatr(i,j,k) = ea(i,j,k) + add_ent + ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent + eatr(i,j,k) = ea_s(i,j,k) + add_ent else - ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k) + ebtr(i,j,k-1) = eb_s(i,j,k-1) ; eatr(i,j,k) = ea_s(i,j,k) endif + if (associated(visc%Kd_extra_S)) then ; if (visc%Kd_extra_S(i,j,k) > 0.0) then add_ent = ((dt * visc%Kd_extra_S(i,j,k)) * GV%m_to_H**2) / & (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + & @@ -1013,13 +1038,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & eatr(i,j,k) = eatr(i,j,k) + add_ent endif ; endif enddo ; enddo - do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ; enddo + do i=is,ie ; eatr(i,j,1) = ea_s(i,j,1) ; enddo enddo ! 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, Hml, dt, G, GV, tv, & + call call_tracer_column_fns(h_prebound, h, ea_s, eb_s, fluxes, Hml, dt, G, GV, tv, & CS%optics, CS%tracer_flow_CSp, CS%debug, & evap_CFL_limit = CS%evap_CFL_limit, & minimum_forcing_depth = CS%minimum_forcing_depth) @@ -1027,7 +1052,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & elseif (associated(visc%Kd_extra_S)) then ! extra diffusivity for passive tracers do j=js,je ; do i=is,ie - ebtr(i,j,nz) = eb(i,j,nz) ; eatr(i,j,1) = ea(i,j,1) + ebtr(i,j,nz) = eb_s(i,j,nz) ; eatr(i,j,1) = ea_s(i,j,1) enddo ; enddo !$OMP parallel do default(shared) private(add_ent) do k=nz,2,-1 ; do j=js,je ; do i=is,ie @@ -1038,8 +1063,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & else add_ent = 0.0 endif - ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent - eatr(i,j,k) = ea(i,j,k) + add_ent + ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent + eatr(i,j,k) = ea_s(i,j,k) + add_ent enddo ; enddo ; enddo ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied @@ -1075,30 +1100,24 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif ! CS%use_sponge -! Save the diapycnal mass fluxes as a diagnostic field. - if (associated(CDp%diapyc_vel)) then - !$OMP parallel do default(shared) - do j=js,je - do K=2,nz ; do i=is,ie - CDp%diapyc_vel(i,j,K) = Idt * (GV%H_to_m * (ea(i,j,k) - eb(i,j,k-1))) - enddo ; enddo - do i=is,ie - CDp%diapyc_vel(i,j,1) = 0.0 - CDp%diapyc_vel(i,j,nz+1) = 0.0 - enddo - enddo - endif - ! Initialize halo regions of ea, eb, and hold to default values. !$OMP parallel do default(shared) do k=1,nz do i=is-1,ie+1 - hold(i,js-1,k) = GV%Angstrom ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0 - hold(i,je+1,k) = GV%Angstrom ; ea(i,je+1,k) = 0.0 ; eb(i,je+1,k) = 0.0 + hold(i,js-1,k) = GV%Angstrom + ea_t(i,js-1,k) = 0.0 ; eb_t(i,js-1,k) = 0.0 + ea_s(i,js-1,k) = 0.0 ; eb_s(i,js-1,k) = 0.0 + hold(i,je+1,k) = GV%Angstrom + ea_t(i,je+1,k) = 0.0 ; eb_t(i,je+1,k) = 0.0 + ea_s(i,je+1,k) = 0.0 ; eb_s(i,je+1,k) = 0.0 enddo do j=js,je - hold(is-1,j,k) = GV%Angstrom ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0 - hold(ie+1,j,k) = GV%Angstrom ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0 + hold(is-1,j,k) = GV%Angstrom + ea_t(is-1,j,k) = 0.0 ; eb_t(is-1,j,k) = 0.0 + ea_s(is-1,j,k) = 0.0 ; eb_s(is-1,j,k) = 0.0 + hold(ie+1,j,k) = GV%Angstrom + ea_t(ie+1,j,k) = 0.0 ; eb_t(ie+1,j,k) = 0.0 + ea_s(ie+1,j,k) = 0.0 ; eb_s(ie+1,j,k) = 0.0 enddo enddo @@ -1106,10 +1125,12 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (G%symmetric) then ; dir_flag = To_All+Omit_Corners else ; dir_flag = To_West+To_South+Omit_Corners ; endif call create_group_pass(CS%pass_hold_eb_ea, hold, G%Domain, dir_flag, halo=1) - call create_group_pass(CS%pass_hold_eb_ea, eb, G%Domain, dir_flag, halo=1) - call create_group_pass(CS%pass_hold_eb_ea, ea, G%Domain, dir_flag, halo=1) + call create_group_pass(CS%pass_hold_eb_ea, eb_t, G%Domain, dir_flag, halo=1) + call create_group_pass(CS%pass_hold_eb_ea, eb_s, G%Domain, dir_flag, halo=1) + call create_group_pass(CS%pass_hold_eb_ea, ea_t, G%Domain, dir_flag, halo=1) + call create_group_pass(CS%pass_hold_eb_ea, ea_s, G%Domain, dir_flag, halo=1) call do_group_pass(CS%pass_hold_eb_ea, G%Domain) - ! visc%Kv_shear is not in the group pass because it has larger vertical extent. + ! visc%Kv_shear and visc%Kv_slow are not in the group pass because it has larger vertical extent. if (associated(visc%Kv_shear)) & call pass_var(visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) if (associated(visc%Kv_slow)) & @@ -1154,12 +1175,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%id_Kd_salt > 0) call post_data(CS%id_Kd_salt, Kd_salt, CS%diag) if (CS%id_Kd_ePBL > 0) call post_data(CS%id_Kd_ePBL, Kd_ePBL, CS%diag) - if (CS%id_ea > 0) call post_data(CS%id_ea, ea, CS%diag) - if (CS%id_eb > 0) call post_data(CS%id_eb, eb, CS%diag) + if (CS%id_ea_t > 0) call post_data(CS%id_ea_t, ea_t, CS%diag) + if (CS%id_eb_t > 0) call post_data(CS%id_eb_t, eb_t, CS%diag) + if (CS%id_ea_s > 0) call post_data(CS%id_ea_s, ea_s, CS%diag) + if (CS%id_eb_s > 0) call post_data(CS%id_eb_s, eb_s, CS%diag) if (CS%id_dudt_dia > 0) call post_data(CS%id_dudt_dia, ADp%du_dt_dia, CS%diag) if (CS%id_dvdt_dia > 0) call post_data(CS%id_dvdt_dia, ADp%dv_dt_dia, CS%diag) - if (CS%id_wd > 0) call post_data(CS%id_wd, CDp%diapyc_vel, CS%diag) if (CS%id_MLD_003 > 0 .or. CS%id_subMLN2 > 0 .or. CS%id_mlotstsq > 0) then call diagnoseMLDbyDensityDifference(CS%id_MLD_003, h, tv, 0.03, G, GV, CS%diag, & @@ -1720,16 +1742,24 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, if (GV%Boussinesq) then ; thickness_units = "m" else ; thickness_units = "kg m-2" ; endif + ! used by layer diabatic CS%id_ea = register_diag_field('ocean_model','ea',diag%axesTL,Time, & 'Layer entrainment from above per timestep','m') CS%id_eb = register_diag_field('ocean_model','eb',diag%axesTL,Time, & 'Layer entrainment from below per timestep', 'm') + + CS%id_ea_t = register_diag_field('ocean_model','ea_t',diag%axesTL,Time, & + 'Layer (heat) entrainment from above per timestep','m') + CS%id_eb_t = register_diag_field('ocean_model','eb_t',diag%axesTL,Time, & + 'Layer (heat) entrainment from below per timestep', 'm') + CS%id_ea_s = register_diag_field('ocean_model','ea_s',diag%axesTL,Time, & + 'Layer (salt) entrainment from above per timestep','m') + CS%id_eb_s = register_diag_field('ocean_model','eb_s',diag%axesTL,Time, & + 'Layer (salt) entrainment from below per timestep', 'm') CS%id_dudt_dia = register_diag_field('ocean_model','dudt_dia',diag%axesCuL,Time, & 'Zonal Acceleration from Diapycnal Mixing', 'm s-2') CS%id_dvdt_dia = register_diag_field('ocean_model','dvdt_dia',diag%axesCvL,Time, & 'Meridional Acceleration from Diapycnal Mixing', 'm s-2') - CS%id_wd = register_diag_field('ocean_model','wd',diag%axesTi,Time, & - 'Diapycnal Velocity', 'm s-1') if (CS%use_int_tides) then CS%id_cg1 = register_diag_field('ocean_model','cn1', diag%axesT1, & Time, 'First baroclinic mode (eigen) speed', 'm s-1') @@ -1799,7 +1829,6 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, if (CS%id_dudt_dia > 0) call safe_alloc_ptr(ADp%du_dt_dia,IsdB,IedB,jsd,jed,nz) if (CS%id_dvdt_dia > 0) call safe_alloc_ptr(ADp%dv_dt_dia,isd,ied,JsdB,JedB,nz) - if (CS%id_wd > 0) call safe_alloc_ptr(CDp%diapyc_vel,isd,ied,jsd,jed,nz+1) ! diagnostics for values prior to diabatic and prior to ALE CS%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, Time, & From 139f6afdd400ab29ac658b119e5d5bc172e543a4 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 1 Jun 2018 08:54:40 -0600 Subject: [PATCH 6/7] Add comments --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 962594ae00..1f0b640a2c 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -560,7 +560,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call cpu_clock_end(id_clock_set_diffusivity) if (showCallTree) call callTree_waypoint("done with set_diffusivity (diabatic)") - ! Set diffusivities for heat and salt + ! Set diffusivities for heat and salt separately !$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat) !$OMP do @@ -568,6 +568,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & Kd_salt(i,j,k) = Kd_int(i,j,k) Kd_heat(i,j,k) = Kd_int(i,j,k) enddo ; enddo ; enddo + ! Add contribution from double diffusion if (associated(visc%Kd_extra_S)) then !$OMP do do k=1,nz+1 ; do j=js,je ; do i=is,ie From b01dda016ce4761193b5aa8c8f41cd4bfe3c9351 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 1 Jun 2018 14:02:24 -0600 Subject: [PATCH 7/7] Clean unecessary layer-related code --- .../vertical/MOM_diabatic_driver.F90 | 180 ++++++------------ 1 file changed, 53 insertions(+), 127 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 1f0b640a2c..0f66625f49 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -293,7 +293,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & 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 +! 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) dSV_dT, & ! The partial derivatives of specific volume with temperature @@ -723,28 +723,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & !$OMP end parallel endif - - ! set ea_t=eb_t=Kd_heat and ea_s=eb_s=Kd_salt on interfaces for use in the - ! tri-diagonal solver. - - do j=js,je ; do i=is,ie - ea_t(i,j,1) = 0.; ea_s(i,j,1) = 0. - enddo ; enddo - -!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_neglect,h,ea_t,ea_s,GV,dt,Kd_salt,Kd_heat,eb_t,eb_s) & -!$OMP private(hval) - do k=2,nz ; do j=js,je ; do i=is,ie - hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) - ea_t(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_heat(i,j,k) - eb_t(i,j,k-1) = ea_t(i,j,k) - ea_s(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_salt(i,j,k) - eb_s(i,j,k-1) = ea_s(i,j,k) - enddo ; enddo ; enddo - do j=js,je ; do i=is,ie - eb_t(i,j,nz) = 0.; eb_s(i,j,nz) = 0. - enddo ; enddo - if (showCallTree) call callTree_waypoint("done setting ea_t,ea_s,eb_t,eb_s from Kd_heat and Kd_salt (diabatic)") - ! Save fields before boundary forcing is applied for tendency diagnostics if (CS%boundary_forcing_tendency_diag) then do k=1,nz ; do j=js,je ; do i=is,ie @@ -789,7 +767,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & Hml(:,:) = visc%MLD(:,:) endif - ! Augment the diffusivities due to those diagnosed in energetic_PBL. + ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL. do K=2,nz ; do j=js,je ; do i=is,ie if (CS%ePBL_is_additive) then Kd_add_here = Kd_ePBL(i,j,K) @@ -798,17 +776,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & Kd_add_here = max(Kd_ePBL(i,j,K) - visc%Kd_shear(i,j,K), 0.0) visc%Kv_shear(i,j,K) = max(visc%Kv_shear(i,j,K), Kd_ePBL(i,j,K)) endif - Ent_int = Kd_add_here * (GV%m_to_H**2 * dt) / & - (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect) - eb_t(i,j,k-1) = eb_t(i,j,k-1) + Ent_int - ea_t(i,j,k) = ea_t(i,j,k) + Ent_int - eb_s(i,j,k-1) = eb_s(i,j,k-1) + Ent_int - ea_s(i,j,k) = ea_s(i,j,k) + Ent_int - Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here - - ! for diagnostics - Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_int(i,j,K) - Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_int(i,j,K) + + Kd_heat(i,j,K) = Kd_heat(i,j,K) + Kd_add_here + Kd_salt(i,j,K) = Kd_salt(i,j,K) + Kd_add_here enddo ; enddo ; enddo @@ -845,54 +815,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (showCallTree) call callTree_waypoint("done with applyBoundaryFluxes (diabatic)") if (CS%debugConservation) call MOM_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, G) - ! 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 - ! against instances where entrainment drives a layer to - ! negative thickness. This situation will never happen if - ! enough iterations are permitted in Calculate_Entrainment. - ! Even if too few iterations are allowed, it is still guarded - ! against. In other words the checks are probably unnecessary. - - ! GMM, should the code below be deleted? eb(i,j,k-1) = ea(i,j,k), - ! see above, so h should not change. - - !$OMP parallel do default(shared) - do j=js,je - do i=is,ie - hold(i,j,1) = h(i,j,1) - h(i,j,1) = h(i,j,1) + (eb_t(i,j,1) - ea_t(i,j,2)) - hold(i,j,nz) = h(i,j,nz) - h(i,j,nz) = h(i,j,nz) + (ea_t(i,j,nz) - eb_t(i,j,nz-1)) - if (h(i,j,1) <= 0.0) then - h(i,j,1) = GV%Angstrom - endif - if (h(i,j,nz) <= 0.0) then - h(i,j,nz) = GV%Angstrom - endif - enddo - do k=2,nz-1 ; do i=is,ie - hold(i,j,k) = h(i,j,k) - h(i,j,k) = h(i,j,k) + ((ea_t(i,j,k) - eb_t(i,j,k-1)) + & - (eb_t(i,j,k) - ea_t(i,j,k+1))) - if (h(i,j,k) <= 0.0) then - h(i,j,k) = GV%Angstrom - endif - enddo ; enddo - enddo - ! Checks for negative thickness may have changed layer thicknesses - call diag_update_remap_grids(CS%diag) - - if (CS%debug) then - call MOM_state_chksum("after negative check ", u, v, h, G, GV, haloshift=0) - call MOM_forcing_chksum("after negative check ", fluxes, G, haloshift=0) - call MOM_thermovar_chksum("after negative check ", tv, G) - endif - if (showCallTree) call callTree_waypoint("done with h=ea-eb (diabatic)") if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G) - ! calculate change in temperature & salinity due to dia-coordinate surface diffusion if (associated(tv%T)) then @@ -920,29 +845,53 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo ; enddo endif - ! Changes T and S via the tridiagonal solver; no change to h - call tracer_vertdiff(hold, ea_t, eb_t, dt, tv%T, G, GV) - call tracer_vertdiff(hold, ea_s, eb_s, dt, tv%S, G, GV) - - ! GMM, with the new approach of having ea,eb for temp and salt - ! only tracer_vertdiff can be used at this time. Therefore, - ! I am commenting the following code. Should this be deleted? - - !if (CS%tracer_tridiag) then - ! call tracer_vertdiff(hold, ea_t, eb_t, dt, tv%T, G, GV) - ! call tracer_vertdiff(hold, ea_s, eb_s, dt, tv%S, G, GV) - !else - ! - ! call triDiagTS(G, GV, is, ie, js, je, hold, ea_t, eb_t, tv%T, tv%S) - !endif - - ! diagnose temperature, salinity, heat, and salt tendencies - ! Note: hold here refers to the thicknesses from before the dual-entraintment when using - ! the bulk mixed layer scheme. Otherwise in ALE-mode, layer thicknesses will have changed - ! In either case, tendencies should be posted on hold + ! set ea_t=eb_t=Kd_heat and ea_s=eb_s=Kd_salt on interfaces for use in the + ! tri-diagonal solver. + + do j=js,je ; do i=is,ie + ea_t(i,j,1) = 0.; ea_s(i,j,1) = 0. + enddo ; enddo + +!$OMP parallel do default(none) shared(is,ie,js,je,nz,h_neglect,h,ea_t,ea_s,GV,dt,Kd_salt,Kd_heat,eb_t,eb_s) & +!$OMP private(hval) + do k=2,nz ; do j=js,je ; do i=is,ie + hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k))) + ea_t(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_heat(i,j,k) + eb_t(i,j,k-1) = ea_t(i,j,k) + ea_s(i,j,k) = (GV%m_to_H**2) * dt * hval * Kd_salt(i,j,k) + eb_s(i,j,k-1) = ea_s(i,j,k) + enddo ; enddo ; enddo + do j=js,je ; do i=is,ie + eb_t(i,j,nz) = 0.; eb_s(i,j,nz) = 0. + enddo ; enddo + if (showCallTree) call callTree_waypoint("done setting ea_t,ea_s,eb_t,eb_s from Kd_heat" //& + "and Kd_salt (diabatic)") + + ! Initialize halo regions of ea, eb, and hold to default values. + !$OMP parallel do default(shared) + do k=1,nz + do i=is-1,ie+1 + ea_t(i,js-1,k) = 0.0 ; eb_t(i,js-1,k) = 0.0 + ea_s(i,js-1,k) = 0.0 ; eb_s(i,js-1,k) = 0.0 + ea_t(i,je+1,k) = 0.0 ; eb_t(i,je+1,k) = 0.0 + ea_s(i,je+1,k) = 0.0 ; eb_s(i,je+1,k) = 0.0 + enddo + do j=js,je + ea_t(is-1,j,k) = 0.0 ; eb_t(is-1,j,k) = 0.0 + ea_s(is-1,j,k) = 0.0 ; eb_s(is-1,j,k) = 0.0 + ea_t(ie+1,j,k) = 0.0 ; eb_t(ie+1,j,k) = 0.0 + ea_s(ie+1,j,k) = 0.0 ; eb_s(ie+1,j,k) = 0.0 + enddo + enddo + + ! Changes T and S via the tridiagonal solver; no change to h + call tracer_vertdiff(h, ea_t, eb_t, dt, tv%T, G, GV) + call tracer_vertdiff(h, ea_s, eb_s, dt, tv%S, G, GV) + + + ! In ALE-mode, layer thicknesses do not change. Therefore, we can use h below if (CS%diabatic_diff_tendency_diag) then - call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, G, GV, CS) - if (CS%id_diabatic_diff_h > 0) call post_data(CS%id_diabatic_diff_h, hold, CS%diag, alt_h = hold) + call diagnose_diabatic_diff_tendency(tv, h, temp_diag, saln_diag, dt, G, GV, CS) endif call cpu_clock_end(id_clock_tridiag) @@ -1033,7 +982,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (associated(visc%Kd_extra_S)) then ; if (visc%Kd_extra_S(i,j,k) > 0.0) then add_ent = ((dt * visc%Kd_extra_S(i,j,k)) * GV%m_to_H**2) / & - (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + & + (0.5 * (h(i,j,k-1) + h(i,j,k)) + & h_neglect) ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent eatr(i,j,k) = eatr(i,j,k) + add_ent @@ -1059,7 +1008,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & do k=nz,2,-1 ; do j=js,je ; do i=is,ie if (visc%Kd_extra_S(i,j,k) > 0.0) then add_ent = ((dt * visc%Kd_extra_S(i,j,k)) * GV%m_to_H**2) / & - (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + & + (0.5 * (h(i,j,k-1) + h(i,j,k)) + & h_neglect) else add_ent = 0.0 @@ -1100,32 +1049,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & endif endif ! CS%use_sponge - - ! Initialize halo regions of ea, eb, and hold to default values. - !$OMP parallel do default(shared) - do k=1,nz - do i=is-1,ie+1 - hold(i,js-1,k) = GV%Angstrom - ea_t(i,js-1,k) = 0.0 ; eb_t(i,js-1,k) = 0.0 - ea_s(i,js-1,k) = 0.0 ; eb_s(i,js-1,k) = 0.0 - hold(i,je+1,k) = GV%Angstrom - ea_t(i,je+1,k) = 0.0 ; eb_t(i,je+1,k) = 0.0 - ea_s(i,je+1,k) = 0.0 ; eb_s(i,je+1,k) = 0.0 - enddo - do j=js,je - hold(is-1,j,k) = GV%Angstrom - ea_t(is-1,j,k) = 0.0 ; eb_t(is-1,j,k) = 0.0 - ea_s(is-1,j,k) = 0.0 ; eb_s(is-1,j,k) = 0.0 - hold(ie+1,j,k) = GV%Angstrom - ea_t(ie+1,j,k) = 0.0 ; eb_t(ie+1,j,k) = 0.0 - ea_s(ie+1,j,k) = 0.0 ; eb_s(ie+1,j,k) = 0.0 - enddo - enddo - call cpu_clock_begin(id_clock_pass) if (G%symmetric) then ; dir_flag = To_All+Omit_Corners else ; dir_flag = To_West+To_South+Omit_Corners ; endif - call create_group_pass(CS%pass_hold_eb_ea, hold, G%Domain, dir_flag, halo=1) call create_group_pass(CS%pass_hold_eb_ea, eb_t, G%Domain, dir_flag, halo=1) call create_group_pass(CS%pass_hold_eb_ea, eb_s, G%Domain, dir_flag, halo=1) call create_group_pass(CS%pass_hold_eb_ea, ea_t, G%Domain, dir_flag, halo=1)