From 99ca9a71878c38255dac041d0b833db4e0065c56 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 10 Apr 2018 12:54:25 -0600 Subject: [PATCH 01/39] Adding option to smooth Ri using a 1-2-1 filter --- .../vertical/MOM_cvmix_shear.F90 | 45 +++++++++++++------ 1 file changed, 31 insertions(+), 14 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_shear.F90 b/src/parameterizations/vertical/MOM_cvmix_shear.F90 index 345522126b..062282f596 100644 --- a/src/parameterizations/vertical/MOM_cvmix_shear.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_shear.F90 @@ -30,9 +30,11 @@ module MOM_cvmix_shear !> Control structure including parameters for CVMix interior shear schemes. type, public :: cvmix_shear_cs logical :: use_LMD94, use_PP81 !< Flags for various schemes + logical :: smooth_ri !< If true, smooth Ri using a 1-2-1 filter real :: Ri_zero !< LMD94 critical Richardson number real :: Nu_zero !< LMD94 maximum interior diffusivity - real :: KPP_exp !< + real :: KPP_exp !< Exponent of unitless factor of diff. + !! for KPP internal shear mixing scheme. real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency (1/s2) real, allocatable, dimension(:,:,:) :: S2 !< Squared shear frequency (1/s2) real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number @@ -52,22 +54,22 @@ module MOM_cvmix_shear !> Subroutine for calculating (internal) vertical diffusivities/viscosities subroutine calculate_cvmix_shear(u_H, v_H, h, tv, kd, & kv, G, GV, CS ) - type(ocean_grid_type), intent(in) :: G !< Grid structure. - type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: u_H !< Initial zonal velocity on T points, in m s-1. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: v_H !< Initial meridional velocity on T points, in m s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. - type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kd !< The vertical diffusivity at each interface - !! (not layer!) in m2 s-1. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kv !< The vertical viscosity at each interface - !! (not layer!) in m2 s-1. - type(cvmix_shear_cs), pointer :: CS !< The control structure returned by a previous call to - !! CVMix_shear_init. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kd !< The vertical diffusivity at each interface + !! (not layer!) in m2 s-1. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: kv !< The vertical viscosity at each interface + !! (not layer!) in m2 s-1. + type(cvmix_shear_cs), pointer :: CS !< The control structure returned by a previous call to + !! CVMix_shear_init. ! Local variables integer :: i, j, k, kk, km1 - real :: gorho - real :: pref, DU, DV, DRHO, DZ, N2, S2 + real :: GoRho + real :: pref, DU, DV, DRHO, DZ, N2, S2, dummy real, dimension(2*(G%ke)) :: pres_1d, temp_1d, salt_1d, rho_1d real, dimension(G%ke+1) :: Ri_Grad !< Gradient Richardson number @@ -120,10 +122,21 @@ subroutine calculate_cvmix_shear(u_H, v_H, h, tv, kd, & ! fill 3d arrays, if user asks for diagsnostics if (CS%id_N2 > 0) CS%N2(i,j,k) = N2 if (CS%id_S2 > 0) CS%S2(i,j,k) = S2 - if (CS%id_ri_grad > 0) CS%ri_grad(i,j,k) = Ri_Grad(k) enddo + ! vertically smooth Ri with 1-2-1 weighting + if (CS%smooth_ri) then + dummy = 0.25 * Ri_grad(1) + Ri_grad(G%ke+1) = Ri_grad(G%ke) + do k = 1, G%ke + Ri_Grad(k) = dummy + 0.5 * Ri_Grad(k) + 0.25 * Ri_grad(k+1) + dummy = 0.25 * Ri_grad(k) + enddo + endif + + if (CS%id_ri_grad > 0) CS%ri_grad(i,j,:) = Ri_Grad(:) + ! Call to CVMix wrapper for computing interior mixing coefficients. call cvmix_coeffs_shear(Mdiff_out=kv(i,j,:), & Tdiff_out=kd(i,j,:), & @@ -209,6 +222,10 @@ logical function cvmix_shear_init(Time, G, GV, param_file, diag, CS) "Exponent of unitless factor of diffusivities,"// & " for KPP internal shear mixing scheme." & ,units="nondim", default=3.0) + call get_param(param_file, mdl, "SMOOTH_RI", CS%smooth_ri, & + "If true, vertically smooth the Richardson"// & + "number by applying a 1-2-1 filter once.", & + default = .false.) call cvmix_init_shear(mix_scheme=CS%mix_scheme, & KPP_nu_zero=CS%Nu_Zero, & KPP_Ri_zero=CS%Ri_zero, & From b14213a020730dab9592669e5d681e39f362ebd0 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 12 Apr 2018 16:45:55 -0600 Subject: [PATCH 02/39] Fill Ri_grad in vanished layers with adjacent value, just when Ri smooth is enabled --- .../vertical/MOM_cvmix_shear.F90 | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_shear.F90 b/src/parameterizations/vertical/MOM_cvmix_shear.F90 index 062282f596..dacb02fe59 100644 --- a/src/parameterizations/vertical/MOM_cvmix_shear.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_shear.F90 @@ -38,8 +38,6 @@ module MOM_cvmix_shear real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency (1/s2) real, allocatable, dimension(:,:,:) :: S2 !< Squared shear frequency (1/s2) real, allocatable, dimension(:,:,:) :: ri_grad !< Gradient Richardson number -! real, allocatable, dimension(:,:,:) :: kv !< vertical viscosity at interface (m2/s) -! real, allocatable, dimension(:,:,:) :: kd !< vertical diffusivity at interface (m2/s) character(10) :: Mix_Scheme !< Mixing scheme name (string) ! Daignostic handles and pointers type(diag_ctrl), pointer :: diag => NULL() @@ -72,7 +70,8 @@ subroutine calculate_cvmix_shear(u_H, v_H, h, tv, kd, & real :: pref, DU, DV, DRHO, DZ, N2, S2, dummy real, dimension(2*(G%ke)) :: pres_1d, temp_1d, salt_1d, rho_1d real, dimension(G%ke+1) :: Ri_Grad !< Gradient Richardson number - + real, parameter :: epsln = 1.e-10 !< Threshold to identify + !! vanished layers ! some constants GoRho = GV%g_Earth / GV%Rho0 @@ -125,8 +124,17 @@ subroutine calculate_cvmix_shear(u_H, v_H, h, tv, kd, & enddo - ! vertically smooth Ri with 1-2-1 weighting + Ri_grad(G%ke+1) = Ri_grad(G%ke) + if (CS%smooth_ri) then + ! 1) fill Ri_grad in vanished layers with adjacent value + do k = 2, G%ke + if (h(i,j,k) .le. epsln) Ri_grad(k) = Ri_grad(k-1) + enddo + + Ri_grad(G%ke+1) = Ri_grad(G%ke) + + ! 2) vertically smooth Ri with 1-2-1 filter dummy = 0.25 * Ri_grad(1) Ri_grad(G%ke+1) = Ri_grad(G%ke) do k = 1, G%ke From c4f1f553c06d16172be35a09e02f0c50953aeb5b Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 12 Apr 2018 17:02:31 -0600 Subject: [PATCH 03/39] Update CVMix --- pkg/CVMix-src | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pkg/CVMix-src b/pkg/CVMix-src index 653d7c39f5..534fc38e75 160000 --- a/pkg/CVMix-src +++ b/pkg/CVMix-src @@ -1 +1 @@ -Subproject commit 653d7c39f50047c9d79c1b15caffe5631dad8bbb +Subproject commit 534fc38e759fcb8a2563fa0dc4c0bbf81f758db9 From 0c363ae6d8d7b6fd32bfd0ba128513730d6431e0 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 13 Apr 2018 08:39:58 -0600 Subject: [PATCH 04/39] Always allocate CS%OBLdepth since other modules may need to know OBLdepth --- src/parameterizations/vertical/MOM_KPP.F90 | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index 87ce532a28..baa33e2ffa 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -370,8 +370,9 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive) CS%id_Vsurf = register_diag_field('ocean_model', 'KPP_Vsurf', diag%axesCv1, Time, & 'j-component flow of surface layer (10% of OBL depth) as passed to [CVmix] KPP', 'm/s') - if (CS%id_OBLdepth > 0) allocate( CS%OBLdepth( SZI_(G), SZJ_(G) ) ) - if (CS%id_OBLdepth > 0) CS%OBLdepth(:,:) = 0. + ! CS%OBLdepth should always be allocated, since it may used by other modules + allocate( CS%OBLdepth( SZI_(G), SZJ_(G) ) ); CS%OBLdepth(:,:) = 0. + if (CS%id_BulkDrho > 0) allocate( CS%dRho( SZI_(G), SZJ_(G), SZK_(G) ) ) if (CS%id_BulkDrho > 0) CS%dRho(:,:,:) = 0. if (CS%id_BulkUz2 > 0) allocate( CS%Uz2( SZI_(G), SZJ_(G), SZK_(G) ) ) @@ -864,7 +865,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & endif ! Copy 1d data into 3d diagnostic arrays - if (CS%id_OBLdepth > 0) CS%OBLdepth(i,j) = OBLdepth_0d + CS%OBLdepth(i,j) = OBLdepth_0d if (CS%id_BulkDrho > 0) CS%dRho(i,j,:) = deltaRho(:) if (CS%id_BulkUz2 > 0) CS%Uz2(i,j,:) = deltaU2(:) if (CS%id_BulkRi > 0) CS%BulkRi(i,j,:) = BulkRi_1d(:) @@ -942,7 +943,7 @@ subroutine KPP_get_BLD(CS, BLD, G) type(KPP_CS), pointer :: CS !< Control structure for !! this module type(ocean_grid_type), intent(in) :: G !< Grid structure - real, dimension(SZI_(G),SZJ_(G)), intent(out) :: BLD!< bnd. layer depth (m) + real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD!< bnd. layer depth (m) ! Local variables integer :: i,j do j = G%jsc, G%jec ; do i = G%isc, G%iec From 60f7d665a6952225a074a58347ad9e7c2700742f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 13 Apr 2018 14:32:56 -0600 Subject: [PATCH 05/39] Initialize visc%Kv_slow and update halos --- src/core/MOM.F90 | 3 +++ src/parameterizations/vertical/MOM_set_viscosity.F90 | 4 +++- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 22dbb86b15..cfbb0c101f 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -2292,6 +2292,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & if (associated(CS%visc%Kv_shear)) & call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1) + if (associated(CS%visc%Kv_slow)) & + call pass_var(CS%visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) + call cpu_clock_end(id_clock_pass_init) call register_obsolete_diagnostics(param_file, CS%diag) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 3df3e7b780..18eb80f280 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -1811,7 +1811,9 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) allocate(visc%Kd_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kd_shear(:,:,:) = 0.0 allocate(visc%TKE_turb(isd:ied,jsd:jed,nz+1)) ; visc%TKE_turb(:,:,:) = 0.0 allocate(visc%Kv_shear(isd:ied,jsd:jed,nz+1)) ; visc%Kv_shear(:,:,:) = 0.0 - allocate(visc%Kv_slow(isd:ied,jsd:jed,nz+1)) ; visc%Kv_slow(:,:,:) = 0.0 + + ! MOM_bkgnd_mixing is always used, so always allocate visc%Kv_slow. GMM + allocate(visc%Kv_slow(isd:ied,jsd:jed,nz+1)) ; visc%Kv_slow(:,:,:) = 0.0 vd = var_desc("Kd_shear","m2 s-1","Shear-driven turbulent diffusivity at interfaces", & hor_grid='h', z_grid='i') From d774e5f05ce348ae388c10304757df6150d76314 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 16 Apr 2018 09:10:19 -0600 Subject: [PATCH 06/39] Add "slow" vertical viscosity in vertvisc_coef --- .../vertical/MOM_vert_friction.F90 | 37 +++++++++++++++++++ 1 file changed, 37 insertions(+) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index ff14a698ed..5154e92c33 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1134,6 +1134,43 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m endif endif + ! add "slow" varying vertical viscosity (e.g., from background, tidal etc) + if (associated(visc%Kv_slow)) then + if (work_on_u) then + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + Kv_add(i,K) = Kv_add(i,K) + 0.5 * (visc%Kv_slow(i,j,k) + visc%Kv_slow(i+1,j,k)) + endif ; enddo ; enddo + if (do_OBCs) then + do I=is,ie ; if (do_i(I) .and. (OBC%segnum_u(I,j) /= OBC_NONE)) then + if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + visc%Kv_slow(i,j,k) ; enddo + elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + visc%Kv_slow(i+1,j,k) ; enddo + endif + endif ; enddo + endif + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + a(i,K) = a(i,K) + Kv_add(i,K) + endif ; enddo ; enddo + else + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + Kv_add(i,K) = Kv_add(i,K) + 0.5*(visc%Kv_slow(i,j,k) + visc%Kv_slow(i,j+1,k)) + endif ; enddo ; enddo + if (do_OBCs) then + do i=is,ie ; if (do_i(i) .and. (OBC%segnum_v(i,J) /= OBC_NONE)) then + if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + visc%Kv_slow(i,j,k) ; enddo + elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + visc%Kv_slow(i,j+1,k) ; enddo + endif + endif ; enddo + endif + do K=2,nz ; do i=is,ie ; if (do_i(i)) then + a(i,K) = a(i,K) + Kv_add(i,K) + endif ; enddo ; enddo + endif + endif + do K=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then ! botfn determines when a point is within the influence of the bottom ! boundary layer, going from 1 at the bottom to 0 in the interior. From fd23f9117335d0542312f932e2fc241fbab25058 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 16 Apr 2018 11:28:19 -0600 Subject: [PATCH 07/39] Set visc%Kv_slow to zero in diabatic --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 70412a716b..8e45d26688 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -378,6 +378,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect*h_neglect Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0 + ! visc%Kv_slow must be set to zero + visc%Kv_slow(:,:,:) = 0.0 if (nz == 1) return showCallTree = callTree_showQuery() From 1488b78945b2789b79bb6c7d0589e08674f897bd Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 16 Apr 2018 11:29:00 -0600 Subject: [PATCH 08/39] Add option to save Kv_slow and Kv --- .../vertical/MOM_vert_friction.F90 | 46 ++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 5154e92c33..e391780253 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -2,7 +2,7 @@ module MOM_vert_friction ! This file is part of MOM6. See LICENSE.md for the license. - +use MOM_domains, only : pass_var use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl use MOM_debugging, only : uvchksum, hchksum @@ -116,6 +116,7 @@ module MOM_vert_friction integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1 integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 integer :: id_Ray_u = -1, id_Ray_v = -1, id_taux_bot = -1, id_tauy_bot = -1 + integer :: id_Kv_slow = -1, id_Kv = -1 !>@} type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() @@ -583,6 +584,10 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) ! based on harmonic mean thicknesses, in m or kg m-2. h_ml ! The mixed layer depth, in m or kg m-2. real, allocatable, dimension(:,:) :: hML_u, hML_v + real, allocatable, dimension(:,:,:) :: Kv_h !< Total vertical viscosity at h-points + real :: av_h !< v-drag coefficient at h-points, in m s-1 + real :: au_h !< u-drag coefficient at h-points, in m s-1 + real :: dh !< average thickness between layers k and k+1, in m or kg m-2. real :: zcol(SZI_(G)) ! The height of an interface at h-points, in m or kg m-2. real :: botfn ! A function which goes from 1 at the bottom to 0 much more ! than Hbbl into the interior. @@ -615,6 +620,10 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) I_Hbbl(:) = 1.0 / (CS%Hbbl * GV%m_to_H + h_neglect) I_valBL = 0.0 ; if (CS%harm_BL_val > 0.0) I_valBL = 1.0 / CS%harm_BL_val + if (CS%id_Kv > 0) then + allocate(Kv_h(SZI_(G), SZJ_(G), SZK_(G)+1)) ; Kv_h(:,:,:) = 0.0 + endif + if (CS%debug .or. (CS%id_hML_u > 0)) then allocate(hML_u(G%IsdB:G%IedB,G%jsd:G%jed)) ; hML_u(:,:) = 0.0 endif @@ -955,6 +964,29 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) endif enddo ! end of v-point j loop + ! Total Kv at h points + if (CS%id_Kv > 0) then + do j = js, je + do i = is, ie + ! set surface and bottom values to zero + Kv_h(i,j,1) = 0.0; Kv_h(i,j,nz+1) = 0.0 + do k=2,nz + av_h = 0.5 * (CS%a_v(i,J,k) + CS%a_v(i,J+1,k)) + au_h = 0.5 * (CS%a_u(I,J,k) + CS%a_u(I+1,j,k)) + dh = 0.5 * (h(i,j,K)+h(i,j,K+1)) + if (dh .le. h_neglect) then + Kv_h(i,j,k) = 0.0 + else + Kv_h(i,j,k) = sqrt(av_h**2 + au_h**2) * dh + if (Kv_h(i,j,k) .lt. 0.0) Kv_h(i,j,k) = 0.0 + endif + enddo + enddo + enddo + ! update halos + call pass_var(Kv_h, G%Domain) + endif + if (CS%debug) then call uvchksum("vertvisc_coef h_[uv]", CS%h_u, & CS%h_v, G%HI,haloshift=0, scale=GV%H_to_m) @@ -966,6 +998,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) endif ! Offer diagnostic fields for averaging. + if (CS%id_Kv_slow > 0) call post_data(CS%id_Kv_slow, visc%Kv_slow, CS%diag) + if (CS%id_Kv > 0) call post_data(CS%id_Kv, Kv_h, CS%diag) if (CS%id_au_vv > 0) call post_data(CS%id_au_vv, CS%a_u, CS%diag) if (CS%id_av_vv > 0) call post_data(CS%id_av_vv, CS%a_v, CS%diag) if (CS%id_h_u > 0) call post_data(CS%id_h_u, CS%h_u, CS%diag) @@ -1660,17 +1694,27 @@ subroutine vertvisc_init(MIS, Time, G, GV, param_file, diag, ADp, dirs, & ALLOC_(CS%a_v(isd:ied,JsdB:JedB,nz+1)) ; CS%a_v(:,:,:) = 0.0 ALLOC_(CS%h_v(isd:ied,JsdB:JedB,nz)) ; CS%h_v(:,:,:) = 0.0 + CS%id_Kv_slow = register_diag_field('ocean_model', 'Kv_slow', diag%axesTi, Time, & + 'Slow varying vertical viscosity', 'm2 s-1') + + CS%id_Kv = register_diag_field('ocean_model', 'Kv', diag%axesTi, Time, & + 'Total vertical viscosity', 'm2 s-1') + CS%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, Time, & 'Zonal Viscous Vertical Coupling Coefficient', 'm s-1') + CS%id_av_vv = register_diag_field('ocean_model', 'av_visc', diag%axesCvi, Time, & 'Meridional Viscous Vertical Coupling Coefficient', 'm s-1') CS%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, Time, & 'Thickness at Zonal Velocity Points for Viscosity', thickness_units) + CS%id_h_v = register_diag_field('ocean_model', 'Hv_visc', diag%axesCvL, Time, & 'Thickness at Meridional Velocity Points for Viscosity', thickness_units) + CS%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', diag%axesCu1, Time, & 'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', thickness_units) + CS%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', diag%axesCv1, Time, & 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', thickness_units) From 69c2c96f9812d74d43c8da2421c88e3c5eb9d05f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 16 Apr 2018 11:34:39 -0600 Subject: [PATCH 09/39] Remove trailing space --- src/parameterizations/vertical/MOM_tidal_mixing.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index a59af01afb..2b3ffb00df 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -1127,7 +1127,7 @@ subroutine setup_tidal_diagnostics(G,CS) integer :: isd, ied, jsd, jed, nz type(tidal_mixing_diags), pointer :: dd - isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed; nz = G%ke + isd = G%isd; ied = G%ied; jsd = G%jsd; jed = G%jed; nz = G%ke dd => CS%dd if ((CS%id_Kd_itidal > 0) .or. (CS%id_Kd_itidal_z > 0) .or. & From d5be1f8c9b0ba118ab9a5e41c8e4fd88148c0cce Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 16 Apr 2018 13:15:56 -0600 Subject: [PATCH 10/39] Attempt to add diag. for total vertical visc. * I believe there is something wrong (halo updates?) with the way this is being done. It needs to be fixed! --- .../vertical/MOM_vert_friction.F90 | 48 ++++++++++--------- 1 file changed, 26 insertions(+), 22 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index e391780253..d90748b820 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -2,7 +2,7 @@ module MOM_vert_friction ! This file is part of MOM6. See LICENSE.md for the license. -use MOM_domains, only : pass_var +use MOM_domains, only : pass_var, To_All, Omit_corners use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr use MOM_diag_mediator, only : diag_ctrl use MOM_debugging, only : uvchksum, hchksum @@ -585,8 +585,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) h_ml ! The mixed layer depth, in m or kg m-2. real, allocatable, dimension(:,:) :: hML_u, hML_v real, allocatable, dimension(:,:,:) :: Kv_h !< Total vertical viscosity at h-points - real :: av_h !< v-drag coefficient at h-points, in m s-1 - real :: au_h !< u-drag coefficient at h-points, in m s-1 + real, dimension(SZI_(G),SZJ_(G)) :: av_h, & !< v-drag coefficient at h-points, in m s-1 + au_h !< u-drag coefficient at h-points, in m s-1 real :: dh !< average thickness between layers k and k+1, in m or kg m-2. real :: zcol(SZI_(G)) ! The height of an interface at h-points, in m or kg m-2. real :: botfn ! A function which goes from 1 at the bottom to 0 much more @@ -966,25 +966,29 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) ! Total Kv at h points if (CS%id_Kv > 0) then - do j = js, je - do i = is, ie - ! set surface and bottom values to zero - Kv_h(i,j,1) = 0.0; Kv_h(i,j,nz+1) = 0.0 - do k=2,nz - av_h = 0.5 * (CS%a_v(i,J,k) + CS%a_v(i,J+1,k)) - au_h = 0.5 * (CS%a_u(I,J,k) + CS%a_u(I+1,j,k)) - dh = 0.5 * (h(i,j,K)+h(i,j,K+1)) - if (dh .le. h_neglect) then - Kv_h(i,j,k) = 0.0 - else - Kv_h(i,j,k) = sqrt(av_h**2 + au_h**2) * dh - if (Kv_h(i,j,k) .lt. 0.0) Kv_h(i,j,k) = 0.0 - endif - enddo - enddo - enddo - ! update halos - call pass_var(Kv_h, G%Domain) + !$OMP parallel do default(shared) + do k=2,nz + ! set surface and bottom values to zero + Kv_h(i,j,1) = 0.0; Kv_h(i,j,nz+1) = 0.0 + do j=js,je ; do I=is-1,ie + au_h(I,j) = CS%a_u(I,J,k) + enddo ; enddo + do J=js-1,je ; do i=is,ie + av_h(i,J) = CS%a_v(i,J,k) + enddo ; enddo + do j = js, je; do i = is, ie + dh = 0.5 * (h(i,j,K)+h(i,j,K+1)) + if (dh .le. h_neglect) then + Kv_h(i,j,k) = 0.0 + else + Kv_h(i,j,k) = sqrt((0.5 * (au_h(I,j)+au_h(I-1,j)))**2 + & + (0.5 * (av_h(i,J) + av_h(i,J-1)))**2) * dh + if (Kv_h(i,j,k) .lt. 0.0) Kv_h(i,j,k) = 0.0 + endif + enddo ; enddo + enddo ! k + ! update halos + call pass_var(Kv_h, G%Domain, To_All+Omit_Corners, halo=1) endif if (CS%debug) then From dcfb722cb334d80d47ece603e6725dd3b30f17bd Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 17 Apr 2018 13:18:22 -0600 Subject: [PATCH 11/39] Bug fix in MOM_bkgnd_mixing The background vertical diff was being added to itself, which is wrong and was leading to a increase in kd_bkgnd over time. This commit fixes this problem. --- src/parameterizations/vertical/MOM_bkgnd_mixing.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 index 14c6c3412e..2c55f4b1c5 100644 --- a/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 +++ b/src/parameterizations/vertical/MOM_bkgnd_mixing.F90 @@ -408,7 +408,7 @@ subroutine calculate_bkgnd_mixing(h, tv, N2_lay, kd_lay, kv, j, G, GV, CS) CS%kd_bkgnd(i,j,1) = 0.0; CS%kv_bkgnd(i,j,1) = 0.0 CS%kd_bkgnd(i,j,nz+1) = 0.0; CS%kv_bkgnd(i,j,nz+1) = 0.0 do k=2,nz - CS%kd_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) + 0.5*(kd_lay(i,j,K-1) + kd_lay(i,j,K)) + CS%kd_bkgnd(i,j,k) = 0.5*(kd_lay(i,j,K-1) + kd_lay(i,j,K)) CS%kv_bkgnd(i,j,k) = CS%kd_bkgnd(i,j,k) * CS%prandtl_bkgnd enddo enddo From d09e8099f8bef3d7912e525522e17f43b020261d Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 17 Apr 2018 13:20:48 -0600 Subject: [PATCH 12/39] Add option to diagnose Kv at u and v points Total vertical viscosity can now be diagnosed at u and v points. --- .../vertical/MOM_vert_friction.F90 | 66 +++++++++---------- 1 file changed, 30 insertions(+), 36 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index d90748b820..95773908aa 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -116,7 +116,7 @@ module MOM_vert_friction integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_au_vv = -1, id_av_vv = -1 integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1 integer :: id_Ray_u = -1, id_Ray_v = -1, id_taux_bot = -1, id_tauy_bot = -1 - integer :: id_Kv_slow = -1, id_Kv = -1 + integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1 !>@} type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() @@ -584,10 +584,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) ! based on harmonic mean thicknesses, in m or kg m-2. h_ml ! The mixed layer depth, in m or kg m-2. real, allocatable, dimension(:,:) :: hML_u, hML_v - real, allocatable, dimension(:,:,:) :: Kv_h !< Total vertical viscosity at h-points - real, dimension(SZI_(G),SZJ_(G)) :: av_h, & !< v-drag coefficient at h-points, in m s-1 - au_h !< u-drag coefficient at h-points, in m s-1 - real :: dh !< average thickness between layers k and k+1, in m or kg m-2. + real, allocatable, dimension(:,:,:) :: Kv_v, & !< Total vertical viscosity at u-points + Kv_u !< Total vertical viscosity at v-points real :: zcol(SZI_(G)) ! The height of an interface at h-points, in m or kg m-2. real :: botfn ! A function which goes from 1 at the bottom to 0 much more ! than Hbbl into the interior. @@ -620,8 +618,12 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) I_Hbbl(:) = 1.0 / (CS%Hbbl * GV%m_to_H + h_neglect) I_valBL = 0.0 ; if (CS%harm_BL_val > 0.0) I_valBL = 1.0 / CS%harm_BL_val - if (CS%id_Kv > 0) then - allocate(Kv_h(SZI_(G), SZJ_(G), SZK_(G)+1)) ; Kv_h(:,:,:) = 0.0 + if (CS%id_Kv_u > 0) then + allocate(Kv_u(G%IsdB:G%IedB,G%jsd:G%jed,G%ke)) ; Kv_u(:,:,:) = 0.0 + endif + + if (CS%id_Kv_v > 0) then + allocate(Kv_v(G%isd:G%ied,G%JsdB:G%JedB,G%ke)) ; Kv_v(:,:,:) = 0.0 endif if (CS%debug .or. (CS%id_hML_u > 0)) then @@ -799,6 +801,13 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) CS%h_u(I,j,k) = hvel(I,k) ; enddo ; enddo endif + ! Diagnose total Kv at u-points + if (CS%id_Kv_u > 0) then + do k=1,nz ; do I=Isq,Ieq + if (do_i(I)) Kv_u(I,j,k) = 0.5 * (CS%a_u(I,j,K)+CS%a_u(I,j,K+1)) * CS%h_u(I,j,k) + enddo ; enddo + endif + enddo @@ -962,34 +971,15 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) do K=1,nz+1 ; do i=is,ie ; if (do_i(i)) CS%a_v(i,J,K) = a(i,K) ; enddo ; enddo do k=1,nz ; do i=is,ie ; if (do_i(i)) CS%h_v(i,J,k) = hvel(i,k) ; enddo ; enddo endif - enddo ! end of v-point j loop - ! Total Kv at h points - if (CS%id_Kv > 0) then - !$OMP parallel do default(shared) - do k=2,nz - ! set surface and bottom values to zero - Kv_h(i,j,1) = 0.0; Kv_h(i,j,nz+1) = 0.0 - do j=js,je ; do I=is-1,ie - au_h(I,j) = CS%a_u(I,J,k) - enddo ; enddo - do J=js-1,je ; do i=is,ie - av_h(i,J) = CS%a_v(i,J,k) + ! Diagnose total Kv at v-points + if (CS%id_Kv_v > 0) then + do k=1,nz ; do i=is,ie + if (do_i(I)) Kv_v(i,J,k) = 0.5 * (CS%a_v(i,J,K)+CS%a_v(i,J,K+1)) * CS%h_v(i,J,k) enddo ; enddo - do j = js, je; do i = is, ie - dh = 0.5 * (h(i,j,K)+h(i,j,K+1)) - if (dh .le. h_neglect) then - Kv_h(i,j,k) = 0.0 - else - Kv_h(i,j,k) = sqrt((0.5 * (au_h(I,j)+au_h(I-1,j)))**2 + & - (0.5 * (av_h(i,J) + av_h(i,J-1)))**2) * dh - if (Kv_h(i,j,k) .lt. 0.0) Kv_h(i,j,k) = 0.0 - endif - enddo ; enddo - enddo ! k - ! update halos - call pass_var(Kv_h, G%Domain, To_All+Omit_Corners, halo=1) - endif + endif + + enddo ! end of v-point j loop if (CS%debug) then call uvchksum("vertvisc_coef h_[uv]", CS%h_u, & @@ -1003,7 +993,8 @@ subroutine vertvisc_coef(u, v, h, forces, visc, dt, G, GV, CS, OBC) ! Offer diagnostic fields for averaging. if (CS%id_Kv_slow > 0) call post_data(CS%id_Kv_slow, visc%Kv_slow, CS%diag) - if (CS%id_Kv > 0) call post_data(CS%id_Kv, Kv_h, CS%diag) + if (CS%id_Kv_u > 0) call post_data(CS%id_Kv_u, Kv_u, CS%diag) + if (CS%id_Kv_v > 0) call post_data(CS%id_Kv_v, Kv_v, CS%diag) if (CS%id_au_vv > 0) call post_data(CS%id_au_vv, CS%a_u, CS%diag) if (CS%id_av_vv > 0) call post_data(CS%id_av_vv, CS%a_v, CS%diag) if (CS%id_h_u > 0) call post_data(CS%id_h_u, CS%h_u, CS%diag) @@ -1701,8 +1692,11 @@ subroutine vertvisc_init(MIS, Time, G, GV, param_file, diag, ADp, dirs, & CS%id_Kv_slow = register_diag_field('ocean_model', 'Kv_slow', diag%axesTi, Time, & 'Slow varying vertical viscosity', 'm2 s-1') - CS%id_Kv = register_diag_field('ocean_model', 'Kv', diag%axesTi, Time, & - 'Total vertical viscosity', 'm2 s-1') + CS%id_Kv_u = register_diag_field('ocean_model', 'Kv_u', diag%axesCuL, Time, & + 'Total vertical viscosity at u-points', 'm2 s-1') + + CS%id_Kv_v = register_diag_field('ocean_model', 'Kv_v', diag%axesCvL, Time, & + 'Total vertical viscosity at v-points', 'm2 s-1') CS%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, Time, & 'Zonal Viscous Vertical Coupling Coefficient', 'm s-1') From 901c3016270f60a77b25e0001513e5aea8cd9569 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 17 Apr 2018 16:16:19 -0600 Subject: [PATCH 13/39] Add option to post diags for bkgnd_mixing --- .../vertical/MOM_set_diffusivity.F90 | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index a1f10519e6..7f192d65d9 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -512,12 +512,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif - ! send bkgnd_mixing diagnostics to post_data - if (CS%bkgnd_mixing_csp%id_kd_bkgnd > 0) & - call post_data(CS%bkgnd_mixing_csp%id_kd_bkgnd, CS%bkgnd_mixing_csp%kd_bkgnd, CS%bkgnd_mixing_csp%diag) - if (CS%bkgnd_mixing_csp%id_kv_bkgnd > 0) & - call post_data(CS%bkgnd_mixing_csp%id_kv_bkgnd, CS%bkgnd_mixing_csp%kv_bkgnd, CS%bkgnd_mixing_csp%diag) - if (CS%Kd_add > 0.0) then if (present(Kd_int)) then !$OMP parallel do default(none) shared(is,ie,js,je,nz,Kd_int,CS,Kd) @@ -538,7 +532,12 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & T_f, S_f, dd%Kd_user) endif - ! GMM, post diags... + ! post diagnostics + if (CS%bkgnd_mixing_csp%id_kd_bkgnd > 0) & + call post_data(CS%bkgnd_mixing_csp%id_kd_bkgnd, CS%bkgnd_mixing_csp%kd_bkgnd, CS%bkgnd_mixing_csp%diag) + if (CS%bkgnd_mixing_csp%id_kv_bkgnd > 0) & + call post_data(CS%bkgnd_mixing_csp%id_kv_bkgnd, CS%bkgnd_mixing_csp%kv_bkgnd, CS%bkgnd_mixing_csp%diag) + if (CS%id_Kd_layer > 0) call post_data(CS%id_Kd_layer, Kd, CS%diag) num_z_diags = 0 From a0358fbbaf7044ac5526f6dad8a054fa69e910bb Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 17 Apr 2018 16:52:26 -0600 Subject: [PATCH 14/39] Fix bug in MOM_cvmix_conv --- src/parameterizations/vertical/MOM_cvmix_conv.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 index 55e7d55d6e..956d0c0de3 100644 --- a/src/parameterizations/vertical/MOM_cvmix_conv.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -212,6 +212,7 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, CS, hbl) iFaceHeight(k+1) = iFaceHeight(k) - dh enddo + ! gets index of the level and interface above hbl kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) call cvmix_coeffs_conv(Mdiff_out=CS%kv_conv(i,j,:), & @@ -224,7 +225,7 @@ subroutine calculate_cvmix_conv(h, tv, G, GV, CS, hbl) OBL_ind=kOBL) ! Do not apply mixing due to convection within the boundary layer - do k=1,NINT(hbl(i,j)) + do k=1,kOBL CS%kv_conv(i,j,k) = 0.0 CS%kd_conv(i,j,k) = 0.0 enddo From 9d6cb46774e83f8eaf2616771c9a70ef5713fa88 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 17 Apr 2018 16:54:55 -0600 Subject: [PATCH 15/39] Rename variable in register_diag_field Variable names are now consistent with what is used in other modules. --- src/parameterizations/vertical/MOM_cvmix_conv.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_cvmix_conv.F90 b/src/parameterizations/vertical/MOM_cvmix_conv.F90 index 956d0c0de3..b385880d7a 100644 --- a/src/parameterizations/vertical/MOM_cvmix_conv.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_conv.F90 @@ -128,11 +128,11 @@ logical function cvmix_conv_init(Time, G, GV, param_file, diag, CS) ! Register diagnostics CS%diag => diag - CS%id_N2 = register_diag_field('ocean_model', 'conv_N2', diag%axesTi, Time, & + CS%id_N2 = register_diag_field('ocean_model', 'N2_conv', diag%axesTi, Time, & 'Square of Brunt-Vaisala frequency used by MOM_cvmix_conv module', '1/s2') - CS%id_kd_conv = register_diag_field('ocean_model', 'conv_kd', diag%axesTi, Time, & + CS%id_kd_conv = register_diag_field('ocean_model', 'kd_conv', diag%axesTi, Time, & 'Additional diffusivity added by MOM_cvmix_conv module', 'm2/s') - CS%id_kv_conv = register_diag_field('ocean_model', 'conv_kv', diag%axesTi, Time, & + CS%id_kv_conv = register_diag_field('ocean_model', 'kv_conv', diag%axesTi, Time, & 'Additional viscosity added by MOM_cvmix_conv module', 'm2/s') call cvmix_init_conv(convect_diff=CS%kd_conv_const, & From 78656bbb545f01cd0a2d71af95ec7072ffea89f7 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 19 Apr 2018 08:49:47 -0600 Subject: [PATCH 16/39] Add missing halo update for visc%Kv_slow --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 8e45d26688..eea1eba16a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -1372,6 +1372,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G ! visc%Kv_shear is 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)) & + call pass_var(visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1) + call cpu_clock_end(id_clock_pass) if (.not. CS%useALEalgorithm) then From 3385857b6fe19e7c7e35a0370914cdc339c131fa Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 19 Apr 2018 08:55:22 -0600 Subject: [PATCH 17/39] Initialize Kd, Kd_int and Kv_slow using interior values specified by user --- src/parameterizations/vertical/MOM_set_diffusivity.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 7f192d65d9..b9905977d5 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -276,6 +276,12 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & call MOM_error(FATAL, "set_diffusivity: visc%Kd_extra_T and "//& "visc%Kd_extra_S must be associated when DOUBLE_DIFFUSION is true.") + ! Set Kd, Kd_int and Kv_slow to constant values. + ! If nothing else is specified, this will be the value used. + Kd(:,:,:) = CS%Kd + Kd_int(:,:,:) = CS%Kd + visc%Kv_slow(:,:,:) = CS%Kv + ! Set up arrays for diagnostics. if ((CS%id_N2 > 0) .or. (CS%id_N2_z > 0)) then From 191b5be35ca028a2aa0608f9202048c5a4a22ca3 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 19 Apr 2018 08:56:19 -0600 Subject: [PATCH 18/39] Add a factor of 2 when adding Kv_slow into Kv_add --- .../vertical/MOM_vert_friction.F90 | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index 95773908aa..973f256915 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1165,16 +1165,17 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m ! add "slow" varying vertical viscosity (e.g., from background, tidal etc) if (associated(visc%Kv_slow)) then + ! GMM/ A factor of 2 is also needed here, see comment above from BGR. if (work_on_u) then do K=2,nz ; do i=is,ie ; if (do_i(i)) then - Kv_add(i,K) = Kv_add(i,K) + 0.5 * (visc%Kv_slow(i,j,k) + visc%Kv_slow(i+1,j,k)) + Kv_add(i,K) = Kv_add(i,K) + 1.0 * (visc%Kv_slow(i,j,k) + visc%Kv_slow(i+1,j,k)) endif ; enddo ; enddo if (do_OBCs) then do I=is,ie ; if (do_i(I) .and. (OBC%segnum_u(I,j) /= OBC_NONE)) then if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then - do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + visc%Kv_slow(i,j,k) ; enddo + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + 2. * visc%Kv_slow(i,j,k) ; enddo elseif (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_W) then - do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + visc%Kv_slow(i+1,j,k) ; enddo + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + 2. * visc%Kv_slow(i+1,j,k) ; enddo endif endif ; enddo endif @@ -1183,14 +1184,14 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m endif ; enddo ; enddo else do K=2,nz ; do i=is,ie ; if (do_i(i)) then - Kv_add(i,K) = Kv_add(i,K) + 0.5*(visc%Kv_slow(i,j,k) + visc%Kv_slow(i,j+1,k)) + Kv_add(i,K) = Kv_add(i,K) + 1.0*(visc%Kv_slow(i,j,k) + visc%Kv_slow(i,j+1,k)) endif ; enddo ; enddo if (do_OBCs) then do i=is,ie ; if (do_i(i) .and. (OBC%segnum_v(i,J) /= OBC_NONE)) then if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then - do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + visc%Kv_slow(i,j,k) ; enddo + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + 2. * visc%Kv_slow(i,j,k) ; enddo elseif (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_S) then - do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + visc%Kv_slow(i,j+1,k) ; enddo + do K=2,nz ; Kv_add(i,K) = Kv_add(i,K) + 2. * visc%Kv_slow(i,j+1,k) ; enddo endif endif ; enddo endif From 02acd12398eec40a2f55b73ecfcdbaaa37d8e6b8 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 16 May 2018 16:11:13 -0600 Subject: [PATCH 19/39] Doxygenize subroutine differential_diffuse_T_S --- .../vertical/MOM_diabatic_aux.F90 | 22 ++++++------------- 1 file changed, 7 insertions(+), 15 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 9588ac3a5c..82799d531a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -239,25 +239,17 @@ subroutine make_frazil(h, tv, G, GV, CS, p_surf) end subroutine make_frazil +!> Applies double diffusion to T & S, assuming no diapycal mass +!! fluxes, using a simple triadiagonal solver. subroutine differential_diffuse_T_S(h, tv, visc, dt, G, GV) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) - type(thermo_var_ptrs), intent(inout) :: tv - type(vertvisc_type), intent(in) :: visc - real, intent(in) :: dt - -! This subroutine applies double diffusion to T & S, assuming no diapycal mass -! fluxes, using a simple triadiagonal solver. - -! Arguments: h - Layer thickness, in m or kg m-2. -! (in) tv - A structure containing pointers to any available -! thermodynamic fields. Absent fields have NULL ptrs. -! (in) visc - A structure containing vertical viscosities, bottom boundary -! layer properies, and related fields. -! (in) dt - Time increment, in s. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. + type(thermo_var_ptrs), intent(inout) :: tv !< pointers to any available modynamic fields. + !! Absent fields have NULL ptrs. + type(vertvisc_type), intent(in) :: visc !< structure containing vertical viscosities, + !! layer properies, and related fields. + real, intent(in) :: dt !< Time increment, in s. real, dimension(SZI_(G)) :: & b1_T, b1_S, & ! Variables used by the tridiagonal solvers of T & S, in H. From 507d34e350a3566cfd68bf15899f889f1b955e63 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 17 May 2018 16:35:53 -0600 Subject: [PATCH 20/39] First version of Double-diffusion via CVMix * Delete the old double-diffusion code --- .../vertical/MOM_cvmix_ddiff.F90 | 314 ++++++++++++++++++ .../vertical/MOM_diabatic_driver.F90 | 26 +- .../vertical/MOM_set_diffusivity.F90 | 282 +++------------- .../vertical/MOM_set_viscosity.F90 | 18 +- 4 files changed, 384 insertions(+), 256 deletions(-) create mode 100644 src/parameterizations/vertical/MOM_cvmix_ddiff.F90 diff --git a/src/parameterizations/vertical/MOM_cvmix_ddiff.F90 b/src/parameterizations/vertical/MOM_cvmix_ddiff.F90 new file mode 100644 index 0000000000..8e52c39849 --- /dev/null +++ b/src/parameterizations/vertical/MOM_cvmix_ddiff.F90 @@ -0,0 +1,314 @@ +!> Interface to CVMix double diffusion scheme. +module MOM_CVMix_ddiff + +! This file is part of MOM6. See LICENSE.md for the license. + +use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field +use MOM_diag_mediator, only : post_data +use MOM_EOS, only : calculate_density_derivs +use MOM_variables, only : thermo_var_ptrs +use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_debugging, only : hchksum +use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_file_parser, only : get_param, log_version, param_file_type +use cvmix_ddiff, only : cvmix_init_ddiff, CVMix_coeffs_ddiff +use cvmix_kpp, only : CVmix_kpp_compute_kOBL_depth +implicit none ; private + +#include + +public CVMix_ddiff_init, CVMix_ddiff_end, CVMix_ddiff_is_used, compute_ddiff_coeffs + +!> Control structure including parameters for CVMix double diffusion. +type, public :: CVMix_ddiff_cs + + ! Parameters + real :: strat_param_max !< maximum value for the stratification parameter (nondim) + real :: kappa_ddiff_s !< leading coefficient in formula for salt-fingering regime + !! for salinity diffusion (m^2/s) + real :: ddiff_exp1 !< interior exponent in salt-fingering regime formula (nondim) + real :: ddiff_exp2 !< exterior exponent in salt-fingering regime formula (nondim) + real :: mol_diff !< molecular diffusivity (m^2/s) + real :: kappa_ddiff_param1 !< exterior coefficient in diffusive convection regime (nondim) + real :: kappa_ddiff_param2 !< middle coefficient in diffusive convection regime (nondim) + real :: kappa_ddiff_param3 !< interior coefficient in diffusive convection regime (nondim) + real :: min_thickness !< Minimum thickness allowed (m) + character(len=4) :: diff_conv_type !< type of diffusive convection to use. Options are Marmorino & + !! Caldwell 1976 ("MC76"; default) and Kelley 1988, 1990 ("K90") + logical :: debug !< If true, turn on debugging + + ! Daignostic handles and pointers + type(diag_ctrl), pointer :: diag => NULL() + integer :: id_KT_extra = -1, id_KS_extra = -1, id_R_rho = -1 + + ! Diagnostics arrays + real, allocatable, dimension(:,:,:) :: KT_extra !< double diffusion diffusivity for temp (m2/s) + real, allocatable, dimension(:,:,:) :: KS_extra !< double diffusion diffusivity for salt (m2/s) + real, allocatable, dimension(:,:,:) :: R_rho !< double-diffusion density ratio (nondim) + +end type CVMix_ddiff_cs + +character(len=40) :: mdl = "MOM_CVMix_ddiff" !< This module's name. + +contains + +!> Initialized the CVMix double diffusion module. +logical function CVMix_ddiff_init(Time, G, GV, param_file, diag, CS) + + type(time_type), intent(in) :: Time !< The current time. + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + type(param_file_type), intent(in) :: param_file !< Run-time parameter file handle + type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. + type(CVMix_ddiff_cs), pointer :: CS !< This module's control structure. + + ! Local variables +! real :: prandtl_conv !< Turbulent Prandtl number used in convective instabilities. + +! This include declares and sets the variable "version". +#include "version_variable.h" + + if (associated(CS)) then + call MOM_error(WARNING, "CVMix_ddiff_init called with an associated "// & + "control structure.") + return + endif + allocate(CS) + + ! Read parameters + call log_version(param_file, mdl, version, & + "Parameterization of mixing due to double diffusion processes via CVMix") + call get_param(param_file, mdl, "USE_CVMIX_DDIFF", CVMix_ddiff_init, & + "If true, turns on double diffusive processes via CVMix. \n"// & + "Note that double diffusive processes on viscosity are ignored \n"// & + "in CVMix, see http://cvmix.github.io/ for justification.",& + default=.false.) + + if (.not. CVMix_ddiff_init) return + + call get_param(param_file, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.) + + call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, do_not_log=.True.) + + call openParameterBlock(param_file,'CVMIX_DDIFF') + + call get_param(param_file, mdl, "STRAT_PARAM_MAX", CS%strat_param_max, & + "The maximum value for the double dissusion stratification parameter", & + units="nondim", default=2.55) + + call get_param(param_file, mdl, "KAPPA_DDIFF_S", CS%kappa_ddiff_s, & + "Leading coefficient in formula for salt-fingering regime \n"// & + "for salinity diffusion.", units="m2 s-1", default=1.0e-4) + + call get_param(param_file, mdl, "DDIFF_EXP1", CS%ddiff_exp1, & + "Interior exponent in salt-fingering regime formula.", & + units="nondim", default=1.0) + + call get_param(param_file, mdl, "DDIFF_EXP2", CS%ddiff_exp2, & + "Exterior exponent in salt-fingering regime formula.", & + units="nondim", default=3.0) + + call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM1", CS%kappa_ddiff_param1, & + "Exterior coefficient in diffusive convection regime.", & + units="nondim", default=0.909) + + call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM2", CS%kappa_ddiff_param2, & + "Middle coefficient in diffusive convection regime.", & + units="nondim", default=4.6) + + call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM3", CS%kappa_ddiff_param3, & + "Interior coefficient in diffusive convection regime.", & + units="nondim", default=-0.54) + + call get_param(param_file, mdl, "MOL_DIFF", CS%mol_diff, & + "Molecular diffusivity used in CVMix double diffusion.", & + units="m2 s-1", default=1.5e-6) + + call get_param(param_file, mdl, "DIFF_CONV_TYPE", CS%diff_conv_type, & + "type of diffusive convection to use. Options are Marmorino \n" //& + "and Caldwell 1976 (MC76) and Kelley 1988, 1990 (K90).", & + default="MC76") + + call closeParameterBlock(param_file) + + ! allocate arrays and set them to zero + ! GMM, dont need the following + !allocate(CS%KT_extra(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%KT_extra(:,:,:) = 0. + !allocate(CS%KS_extra(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%KS_extra(:,:,:) = 0. + + ! Register diagnostics + CS%diag => diag + + CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for temperature', 'm2 s-1') + + CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for salinity', 'm2 s-1') + + CS%id_R_rho = register_diag_field('ocean_model','R_rho',diag%axesTi,Time, & + 'Double-diffusion density ratio', 'nondim') + if (CS%id_R_rho > 0) & + allocate(CS%R_rho( SZI_(G), SZJ_(G), SZK_(G)+1)); CS%R_rho(:,:,:) = 0.0 + + call cvmix_init_ddiff(strat_param_max=CS%strat_param_max, & + kappa_ddiff_s=CS%kappa_ddiff_s, & + ddiff_exp1=CS%ddiff_exp1, & + ddiff_exp2=CS%ddiff_exp2, & + mol_diff=CS%mol_diff, & + kappa_ddiff_param1=CS%kappa_ddiff_param1, & + kappa_ddiff_param2=CS%kappa_ddiff_param2, & + kappa_ddiff_param3=CS%kappa_ddiff_param3, & + diff_conv_type=CS%diff_conv_type) + +end function CVMix_ddiff_init + +!> Subroutine for computing vertical diffusion coefficients for the +!! double diffusion mixing parameterization. +subroutine compute_ddiff_coeffs(h, tv, G, GV, j, Kd_T, Kd_S, CS) + + type(ocean_grid_type), intent(in) :: G !< Grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. + type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: Kd_T !< Interface double diffusion diapycnal +! real, dimension(:,:,:), pointer :: Kd_T +! real, dimension(:,:,:), pointer :: Kd_S + !! diffusivity for temp (m2/sec). + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: Kd_S !< Interface double diffusion diapycnal + !! diffusivity for salt (m2/sec). + type(CVMix_ddiff_cs), pointer :: CS !< The control structure returned + !! by a previous call to CVMix_ddiff_init. + integer, intent(in) :: j !< Meridional grid indice. +! real, dimension(:,:), optional, pointer :: hbl !< Depth of ocean boundary layer (m) + + ! local variables + real, dimension(SZK_(G)) :: & + cellHeight, & !< Height of cell centers (m) + dRho_dT, & !< partial derivatives of density wrt temp (kg m-3 degC-1) + dRho_dS, & !< partial derivatives of density wrt saln (kg m-3 PPT-1) + pres_int, & !< pressure at each interface (Pa) + temp_int, & !< temp and at interfaces (degC) + salt_int, & !< salt at at interfaces + alpha_dT, & !< alpha*dT across interfaces + beta_dS, & !< beta*dS across interfaces + dT, & !< temp. difference between adjacent layers (degC) + dS !< salt difference between adjacent layers + + real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) + integer :: kOBL !< level of OBL extent + real :: pref, g_o_rho0, rhok, rhokm1, dz, dh, hcorr + integer :: i, k + + ! initialize dummy variables + pres_int(:) = 0.0; temp_int(:) = 0.0; salt_int(:) = 0.0 + alpha_dT(:) = 0.0; beta_dS(:) = 0.0; dRho_dT(:) = 0.0 + dRho_dS(:) = 0.0; dT(:) = 0.0; dS(:) = 0.0 + + ! set Kd_T and Kd_S to zero to avoid passing values from previous call + Kd_T(:,j,:) = 0.0; Kd_S(:,j,:) = 0.0 + + ! GMM, check this. + !if (.not. associated(hbl)) then + ! allocate(hbl(SZI_(G), SZJ_(G))); + ! hbl(:,:) = 0.0 + !endif + + do i = G%isc, G%iec + + ! skip calling at land points + if (G%mask2dT(i,j) == 0.) cycle + + pRef = 0. + pres_int(1) = pRef + ! we don't have SST and SSS, so let's use values at top-most layer + temp_int(1) = TV%T(i,j,1); salt_int(1) = TV%S(i,j,1) + do k=2,G%ke + ! pressure at interface + pres_int(k) = pRef + GV%H_to_Pa * h(i,j,k-1) + ! temp and salt at interface + ! for temp: (t1*h1 + t2*h2)/(h1+h2) + temp_int(k) = (TV%T(i,j,k-1)*h(i,j,k-1) + TV%T(i,j,k)*h(i,j,k))/(h(i,j,k-1)+h(i,j,k)) + salt_int(k) = (TV%S(i,j,k-1)*h(i,j,k-1) + TV%S(i,j,k)*h(i,j,k))/(h(i,j,k-1)+h(i,j,k)) + ! dT and dS + dT(k) = (TV%T(i,j,k-1)-TV%T(i,j,k)) + dS(k) = (TV%S(i,j,k-1)-TV%S(i,j,k)) + pRef = pRef + GV%H_to_Pa * h(i,j,k-1) + enddo ! k-loop finishes + + call calculate_density_derivs(temp_int(:), salt_int(:), pres_int(:), drho_dT(:), drho_dS(:), 1, G%ke, TV%EQN_OF_STATE) + + ! GMM, explain need for -1 in front of alpha + ! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then ! salt finger case + ! if ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then ! diffusive convection + do k=1,G%ke + alpha_dT(k) = -1.0*drho_dT(k) * dT(k) + beta_dS(k) = drho_dS(k) * dS(k) + enddo + + if (CS%id_R_rho > 0.0) then + do k=1,G%ke + CS%R_rho(i,j,k) = alpha_dT(k)/beta_dS(k) + enddo + endif + + iFaceHeight(1) = 0.0 ! BBL is all relative to the surface + hcorr = 0.0 + ! compute heights at cell center and interfaces + do k=1,G%ke + dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment + dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0) + hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0 + dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness + cellHeight(k) = iFaceHeight(k) - 0.5 * dh + iFaceHeight(k+1) = iFaceHeight(k) - dh + enddo + + ! gets index of the level and interface above hbl + !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j)) + + call CVMix_coeffs_ddiff(Tdiff_out=Kd_T(i,j,:), & + Sdiff_out=Kd_S(i,j,:), & + strat_param_num=alpha_dT(:), & + strat_param_denom=beta_dS(:), & + nlev=G%ke, & + max_nlev=G%ke) + + !if (is_root_pe()) then + ! write(*,*)'drho_dT, alpha_dT', & + ! drho_dT(:), alpha_dT(:) + ! write(*,*)'drho_dS, beta_dS', & + ! drho_dS(:), beta_dS(:) + !endif + + ! Do not apply mixing due to convection within the boundary layer + !do k=1,kOBL + ! Kd_T(i,j,k) = 0.0 + ! Kd_S(i,j,k) = 0.0 + !enddo + + enddo ! i-loop + +end subroutine compute_ddiff_coeffs + +!> Reads the parameter "USE_CVMIX_DDIFF" and returns state. +!! This function allows other modules to know whether this parameterization will +!! be used without needing to duplicate the log entry. +logical function CVMix_ddiff_is_used(param_file) + type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters + call get_param(param_file, mdl, "USE_CVMIX_DDIFF", CVMix_ddiff_is_used, & + default=.false., do_not_log = .true.) + +end function CVMix_ddiff_is_used + +!> Clear pointers and dealocate memory +subroutine CVMix_ddiff_end(CS) + type(CVMix_ddiff_cs), pointer :: CS ! Control structure + + deallocate(CS) + +end subroutine CVMix_ddiff_end + + +end module MOM_CVMix_ddiff diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index eea1eba16a..b109d3642a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -9,7 +9,8 @@ module MOM_diabatic_driver use MOM_checksum_packages, only : MOM_state_chksum, MOM_state_stats use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE -use MOM_CVMix_shear, only : cvmix_shear_is_used +use MOM_CVMix_shear, only : CVMix_shear_is_used +use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used use MOM_diabatic_aux, only : diabatic_aux_init, diabatic_aux_end, diabatic_aux_CS use MOM_diabatic_aux, only : make_frazil, adjust_salt, insert_brine, differential_diffuse_T_S, triDiagTS use MOM_diabatic_aux, only : find_uv_at_h, diagnoseMLDbyDensityDifference, applyBoundaryFluxesInOut @@ -90,8 +91,9 @@ module MOM_diabatic_driver !! in the surface boundary layer. logical :: use_kappa_shear !< If true, use the kappa_shear module to find the !! shear-driven diapycnal diffusivity. - logical :: use_cvmix_shear !< If true, use the CVMix module to find the + logical :: use_CVMix_shear !< If true, use the CVMix module to find the !! shear-driven diapycnal diffusivity. + logical :: use_CVMix_ddiff !< If true, use the CVMix double diffusion module. logical :: use_tidal_mixing !< If true, activate tidal mixing diffusivity. logical :: use_cvmix_conv !< If true, use the CVMix module to get enhanced !! mixing due to convection. @@ -244,7 +246,7 @@ module MOM_diabatic_driver integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap -integer :: id_clock_kpp +integer :: id_clock_kpp, id_clock_CVMix_ddiff contains @@ -721,10 +723,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G ! a diffusivity and happen before KPP. But generally in MOM, we do not match ! KPP boundary layer to interior, so this diffusivity can be computed when convenient. if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T)) then - call cpu_clock_begin(id_clock_differential_diff) + call cpu_clock_begin(id_clock_CVMix_ddiff) call differential_diffuse_T_S(h, tv, visc, dt, G, GV) - call cpu_clock_end(id_clock_differential_diff) + call cpu_clock_end(id_clock_CVMix_ddiff) if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)") if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, G) @@ -737,7 +739,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G enddo ; enddo ; enddo endif - endif @@ -1872,7 +1873,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, real :: Kd integer :: num_mode - logical :: use_temperature, differentialDiffusion + logical :: use_temperature type(vardesc) :: vd ! This "include" declares and sets the variable "version". @@ -1924,11 +1925,10 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "If true, the diffusivity from ePBL is added to all\n"//& "other diffusivities. Otherwise, the larger of kappa-\n"//& "shear and ePBL diffusivities are used.", default=.true.) - call get_param(param_file, mod, "DOUBLE_DIFFUSION", differentialDiffusion, & - "If true, apply parameterization of double-diffusion.", & - default=.false. ) + CS%use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) CS%use_kappa_shear = kappa_shear_is_used(param_file) - CS%use_cvmix_shear = cvmix_shear_is_used(param_file) + CS%use_CVMix_shear = CVMix_shear_is_used(param_file) + if (CS%bulkmixedlayer) then call get_param(param_file, mod, "ML_MIX_FIRST", CS%ML_mix_first, & "The fraction of the mixed layer mixing that is applied \n"//& @@ -2384,8 +2384,8 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=CLOCK_MODULE) id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=CLOCK_ROUTINE) id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=CLOCK_ROUTINE) - id_clock_differential_diff = -1 ; if (differentialDiffusion) & - id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=CLOCK_ROUTINE) + id_clock_CVMix_ddiff = -1 ; if (CS%use_CVMix_ddiff) & + id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion)', grain=CLOCK_ROUTINE) ! initialize the auxiliary diabatic driver module call diabatic_aux_init(Time, G, GV, param_file, diag, CS%diabatic_aux_CSp, & diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index b9905977d5..f33bdea2a8 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -21,8 +21,10 @@ module MOM_set_diffusivity use MOM_intrinsic_functions, only : invcosh use MOM_io, only : slasher, vardesc, var_desc, MOM_read_data use MOM_kappa_shear, only : calculate_kappa_shear, kappa_shear_init, Kappa_shear_CS -use MOM_cvmix_shear, only : calculate_cvmix_shear, cvmix_shear_init, cvmix_shear_cs -use MOM_cvmix_shear, only : cvmix_shear_end +use MOM_CVMix_shear, only : calculate_CVMix_shear, CVMix_shear_init, CVMix_shear_cs +use MOM_CVMix_shear, only : CVMix_shear_end +use MOM_CVMix_ddiff, only : CVMix_ddiff_init, CVMix_ddiff_end, CVMix_ddiff_cs +use MOM_CVMix_ddiff, only : compute_ddiff_coeffs use MOM_bkgnd_mixing, only : calculate_bkgnd_mixing, bkgnd_mixing_init, bkgnd_mixing_cs use MOM_bkgnd_mixing, only : bkgnd_mixing_end, sfc_bkgnd_mixing use MOM_string_functions, only : uppercase @@ -129,18 +131,16 @@ module MOM_set_diffusivity ! shear-driven diapycnal diffusivity. logical :: use_cvmix_shear ! If true, use one of the CVMix modules to find ! shear-driven diapycnal diffusivity. - logical :: double_diffusion ! If true, enable double-diffusive mixing. + logical :: use_CVMix_ddiff ! If true, enable double-diffusive mixing via CVMix. logical :: simple_TKE_to_Kd ! If true, uses a simple estimate of Kd/TKE that ! does not rely on a layer-formulation. - real :: Max_Rrho_salt_fingers ! max density ratio for salt fingering - real :: Max_salt_diff_salt_fingers ! max salt diffusivity for salt fingers (m2/s) - real :: Kv_molecular ! molecular visc for double diff convect (m2/s) character(len=200) :: inputdir type(user_change_diff_CS), pointer :: user_change_diff_CSp => NULL() type(diag_to_Z_CS), pointer :: diag_to_Z_CSp => NULL() type(Kappa_shear_CS), pointer :: kappaShear_CSp => NULL() - type(cvmix_shear_cs), pointer :: cvmix_shear_csp => NULL() + type(CVMix_shear_cs), pointer :: CVMix_shear_csp => NULL() + type(CVMix_ddiff_cs), pointer :: CVMix_ddiff_csp => NULL() type(bkgnd_mixing_cs), pointer :: bkgnd_mixing_csp => NULL() type(int_tide_CS), pointer :: int_tide_CSp => NULL() type(tidal_mixing_cs), pointer :: tm_csp => NULL() @@ -158,11 +158,6 @@ module MOM_set_diffusivity 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 - end type set_diffusivity_CS type diffusivity_diags @@ -172,12 +167,9 @@ module MOM_set_diffusivity Kd_BBL => NULL(),& ! BBL diffusivity at interfaces (m2/s) Kd_work => NULL(),& ! layer integrated work by diapycnal mixing (W/m2) maxTKE => NULL(),& ! energy required to entrain to h_max (m3/s3) - TKE_to_Kd => NULL(),& ! conversion rate (~1.0 / (G_Earth + dRho_lay)) + TKE_to_Kd => NULL() ! conversion rate (~1.0 / (G_Earth + dRho_lay)) ! between TKE dissipated within a layer and Kd ! in that layer, in m2 s-1 / m3 s-3 = s2 m-1 - KT_extra => NULL(),& ! double diffusion diffusivity for temp (m2/s) - KS_extra => NULL() ! double diffusion diffusivity for saln (m2/s) - end type diffusivity_diags ! Clocks @@ -226,17 +218,15 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! massless layers filled vertically by diffusion. real, dimension(SZI_(G),SZK_(G)) :: & - N2_lay, & ! squared buoyancy frequency associated with layers (1/s2) - maxTKE, & ! energy required to entrain to h_max (m3/s3) - TKE_to_Kd ! conversion rate (~1.0 / (G_Earth + dRho_lay)) between - ! TKE dissipated within a layer and Kd in that layer, in - ! m2 s-1 / m3 s-3 = s2 m-1. + N2_lay, & !< squared buoyancy frequency associated with layers (1/s2) + maxTKE, & !< energy required to entrain to h_max (m3/s3) + TKE_to_Kd !< conversion rate (~1.0 / (G_Earth + dRho_lay)) between + !< TKE dissipated within a layer and Kd in that layer, in + !< m2 s-1 / m3 s-3 = s2 m-1. real, dimension(SZI_(G),SZK_(G)+1) :: & - N2_int, & ! squared buoyancy frequency associated at interfaces (1/s2) - dRho_int, & ! locally ref potential density difference across interfaces (in s-2) smg: or kg/m3? - KT_extra, & ! double difusion diffusivity on temperature (m2/sec) - KS_extra ! double difusion diffusivity on salinity (m2/sec) + N2_int, & !< squared buoyancy frequency associated at interfaces (1/s2) + dRho_int !< locally ref potential density difference across interfaces (in s-2) smg: or kg/m3? real :: I_Rho0 ! inverse of Boussinesq density (m3/kg) real :: dissip ! local variable for dissipation calculations (W/m3) @@ -271,10 +261,10 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & use_EOS = associated(tv%eqn_of_state) - if ((CS%double_diffusion) .and. & + if ((CS%use_CVMix_ddiff) .and. & .not.(associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S)) ) & call MOM_error(FATAL, "set_diffusivity: visc%Kd_extra_T and "//& - "visc%Kd_extra_S must be associated when DOUBLE_DIFFUSION is true.") + "visc%Kd_extra_S must be associated when USE_CVMIX_DDIFF is true.") ! Set Kd, Kd_int and Kv_slow to constant values. ! If nothing else is specified, this will be the value used. @@ -299,12 +289,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%id_TKE_to_Kd > 0) then allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz)) ; dd%TKE_to_Kd(:,:,:) = 0.0 endif - if ((CS%id_KT_extra > 0) .or. (CS%id_KT_extra_z > 0)) then - allocate(dd%KT_extra(isd:ied,jsd:jed,nz+1)) ; dd%KT_extra(:,:,:) = 0.0 - endif - if ((CS%id_KS_extra > 0) .or. (CS%id_KS_extra_z > 0)) then - allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1)) ; dd%KS_extra(:,:,:) = 0.0 - endif if ((CS%id_Kd_BBL > 0) .or. (CS%id_Kd_BBL_z > 0)) then allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1)) ; dd%Kd_BBL(:,:,:) = 0.0 endif @@ -376,35 +360,13 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & do K=1,nz+1 ; do i=is,ie ; dd%N2_3d(i,j,K) = N2_int(i,K) ; enddo ; enddo endif - ! add background mixing + ! Add background mixing call calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc%Kv_slow, j, G, GV, CS%bkgnd_mixing_csp) - ! GMM, the following will go into the MOM_cvmix_double_diffusion module - if (CS%double_diffusion) then - call double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, KT_extra, KS_extra) - do K=2,nz ; do i=is,ie - if (KS_extra(i,K) > KT_extra(i,K)) then ! salt fingering - Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KT_extra(i,K) - Kd(i,j,k) = Kd(i,j,k) + 0.5*KT_extra(i,K) - visc%Kd_extra_S(i,j,k) = KS_extra(i,K)-KT_extra(i,K) - visc%Kd_extra_T(i,j,k) = 0.0 - elseif (KT_extra(i,K) > 0.0) then ! double-diffusive convection - Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KS_extra(i,K) - Kd(i,j,k) = Kd(i,j,k) + 0.5*KS_extra(i,K) - visc%Kd_extra_T(i,j,k) = KT_extra(i,K)-KS_extra(i,K) - visc%Kd_extra_S(i,j,k) = 0.0 - else ! There is no double diffusion at this interface. - visc%Kd_extra_T(i,j,k) = 0.0 - visc%Kd_extra_S(i,j,k) = 0.0 - endif - enddo; enddo - if (associated(dd%KT_extra)) then ; do K=1,nz+1 ; do i=is,ie - dd%KT_extra(i,j,K) = KT_extra(i,K) - enddo ; enddo ; endif - - if (associated(dd%KS_extra)) then ; do K=1,nz+1 ; do i=is,ie - dd%KS_extra(i,j,K) = KS_extra(i,K) - enddo ; enddo ; endif + ! Apply double diffusion + ! GMM, we need to pass HBL to compute_ddiff_coeffs, but it is not yet available. + if (CS%use_CVMix_ddiff) then + call compute_ddiff_coeffs(h, tv, G, GV, j, visc%Kd_extra_T, visc%Kd_extra_S, CS%CVMix_ddiff_csp) endif ! Add the input turbulent diffusivity. @@ -502,6 +464,11 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%useKappaShear) call hchksum(visc%Kd_shear,"Turbulent Kd",G%HI,haloshift=0) + if (CS%use_CVMix_ddiff) then + call hchksum(visc%Kd_extra_T, "MOM_set_diffusivity: Kd_extra_T",G%HI,haloshift=0) + call hchksum(visc%Kd_extra_S, "MOM_set_diffusivity: Kd_extra_S",G%HI,haloshift=0) + endif + if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then call uvchksum("BBL Kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, & G%HI, 0, symmetric=.true.) @@ -539,17 +506,27 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif ! post diagnostics + + ! background mixing if (CS%bkgnd_mixing_csp%id_kd_bkgnd > 0) & call post_data(CS%bkgnd_mixing_csp%id_kd_bkgnd, CS%bkgnd_mixing_csp%kd_bkgnd, CS%bkgnd_mixing_csp%diag) if (CS%bkgnd_mixing_csp%id_kv_bkgnd > 0) & call post_data(CS%bkgnd_mixing_csp%id_kv_bkgnd, CS%bkgnd_mixing_csp%kv_bkgnd, CS%bkgnd_mixing_csp%diag) - if (CS%id_Kd_layer > 0) call post_data(CS%id_Kd_layer, Kd, CS%diag) + ! double diffusive mixing + if (CS%CVMix_ddiff_csp%id_KT_extra > 0) & + call post_data(CS%CVMix_ddiff_csp%id_KT_extra, visc%Kd_extra_T, CS%CVMix_ddiff_csp%diag) + if (CS%CVMix_ddiff_csp%id_KS_extra > 0) & + call post_data(CS%CVMix_ddiff_csp%id_KS_extra, visc%Kd_extra_S, CS%CVMix_ddiff_csp%diag) + if (CS%CVMix_ddiff_csp%id_R_rho > 0) & + call post_data(CS%CVMix_ddiff_csp%id_R_rho, CS%CVMix_ddiff_csp%R_rho, CS%CVMix_ddiff_csp%diag) - num_z_diags = 0 + if (CS%id_Kd_layer > 0) call post_data(CS%id_Kd_layer, Kd, CS%diag) + ! tidal mixing call post_tidal_diagnostics(G,GV,h,CS%tm_csp) + num_z_diags = 0 if (CS%tm_csp%Int_tide_dissipation .or. CS%tm_csp%Lee_wave_dissipation .or. & CS%tm_csp%Lowmode_itidal_dissipation) then @@ -573,26 +550,11 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif - if (CS%id_KT_extra > 0) call post_data(CS%id_KT_extra, dd%KT_extra, CS%diag) - if (CS%id_KS_extra > 0) call post_data(CS%id_KS_extra, dd%KS_extra, CS%diag) if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, dd%Kd_BBL, CS%diag) - if (CS%id_KT_extra_z > 0) then - num_z_diags = num_z_diags + 1 - z_ids(num_z_diags) = CS%id_KT_extra_z - z_ptrs(num_z_diags)%p => dd%KT_extra - endif - - if (CS%id_KS_extra_z > 0) then - num_z_diags = num_z_diags + 1 - z_ids(num_z_diags) = CS%id_KS_extra_z - z_ptrs(num_z_diags)%p => dd%KS_extra - endif - if (CS%id_Kd_BBL_z > 0) then num_z_diags = num_z_diags + 1 z_ids(num_z_diags) = CS%id_Kd_BBL_z - z_ptrs(num_z_diags)%p => dd%KS_extra endif if (num_z_diags > 0) & @@ -603,8 +565,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (associated(dd%Kd_user)) deallocate(dd%Kd_user) if (associated(dd%maxTKE)) deallocate(dd%maxTKE) if (associated(dd%TKE_to_Kd)) deallocate(dd%TKE_to_Kd) - if (associated(dd%KT_extra)) deallocate(dd%KT_extra) - if (associated(dd%KS_extra)) deallocate(dd%KS_extra) if (associated(dd%Kd_BBL)) deallocate(dd%Kd_BBL) if (showCallTree) call callTree_leave("set_diffusivity()") @@ -956,119 +916,6 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, & end subroutine find_N2 -! GMM, the following will be moved to a new module - -!> This subroutine sets the additional diffusivities of temperature and -!! salinity due to double diffusion, using the same functional form as is -!! used in MOM4.1, and taken from an NCAR technical note (REF?) that updates -!! what was in Large et al. (1994). All the coefficients here should probably -!! be made run-time variables rather than hard-coded constants. -!! -!! \todo Find reference for NCAR tech note above. -subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, Kd_T_dd, Kd_S_dd) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available - !! thermodynamic fields; absent fields have NULL - !! ptrs. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: T_f !< layer temp in C with the values in massless layers - !! filled vertically by diffusion. - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: S_f !< Layer salinities in PPT with values in massless - !! layers filled vertically by diffusion. - integer, intent(in) :: j !< Meridional index upon which to work. - type(set_diffusivity_CS), pointer :: CS !< Module control structure. - real, dimension(SZI_(G),SZK_(G)+1), & - intent(out) :: Kd_T_dd !< Interface double diffusion diapycnal - !! diffusivity for temp (m2/sec). - real, dimension(SZI_(G),SZK_(G)+1), & - intent(out) :: Kd_S_dd !< Interface double diffusion diapycnal - !! diffusivity for saln (m2/sec). - -! Arguments: -! (in) tv - structure containing pointers to any available -! thermodynamic fields; absent fields have NULL ptrs -! (in) h - layer thickness (m or kg m-2) -! (in) T_f - layer temp in C with the values in massless layers -! filled vertically by diffusion -! (in) S_f - layer salinities in PPT with values in massless layers -! filled vertically by diffusion -! (in) G - ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) CS - module control structure -! (in) j - meridional index upon which to work -! (out) Kd_T_dd - interface double diffusion diapycnal diffusivity for temp (m2/sec) -! (out) Kd_S_dd - interface double diffusion diapycnal diffusivity for saln (m2/sec) - -! This subroutine sets the additional diffusivities of temperature and -! salinity due to double diffusion, using the same functional form as is -! used in MOM4.1, and taken from an NCAR technical note (###REF?) that updates -! what was in Large et al. (1994). All the coefficients here should probably -! be made run-time variables rather than hard-coded constants. - - real, dimension(SZI_(G)) :: & - dRho_dT, & ! partial derivatives of density wrt temp (kg m-3 degC-1) - dRho_dS, & ! partial derivatives of density wrt saln (kg m-3 PPT-1) - pres, & ! pressure at each interface (Pa) - Temp_int, & ! temp and saln at interfaces - Salin_int - - real :: alpha_dT ! density difference between layers due to temp diffs (kg/m3) - real :: beta_dS ! density difference between layers due to saln diffs (kg/m3) - - real :: Rrho ! vertical density ratio - real :: diff_dd ! factor for double-diffusion - real :: prandtl ! flux ratio for diffusive convection regime - - real, parameter :: Rrho0 = 1.9 ! limit for double-diffusive density ratio - real, parameter :: dsfmax = 1.e-4 ! max diffusivity in case of salt fingering - real, parameter :: Kv_molecular = 1.5e-6 ! molecular viscosity (m2/sec) - - integer :: i, k, is, ie, nz - is = G%isc ; ie = G%iec ; nz = G%ke - - if (associated(tv%eqn_of_state)) then - do i=is,ie - pres(i) = 0.0 ; Kd_T_dd(i,1) = 0.0 ; Kd_S_dd(i,1) = 0.0 - Kd_T_dd(i,nz+1) = 0.0 ; Kd_S_dd(i,nz+1) = 0.0 - enddo - do K=2,nz - do i=is,ie - pres(i) = pres(i) + GV%H_to_Pa*h(i,j,k-1) - Temp_Int(i) = 0.5 * (T_f(i,j,k-1) + T_f(i,j,k)) - Salin_Int(i) = 0.5 * (S_f(i,j,k-1) + S_f(i,j,k)) - enddo - call calculate_density_derivs(Temp_int, Salin_int, pres, & - dRho_dT(:), dRho_dS(:), is, ie-is+1, tv%eqn_of_state) - - do i=is,ie - alpha_dT = -1.0*dRho_dT(i) * (T_f(i,j,k-1) - T_f(i,j,k)) - beta_dS = dRho_dS(i) * (S_f(i,j,k-1) - S_f(i,j,k)) - - if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then ! salt finger case - Rrho = min(alpha_dT/beta_dS,Rrho0) - diff_dd = 1.0 - ((RRho-1.0)/(RRho0-1.0)) - diff_dd = dsfmax*diff_dd*diff_dd*diff_dd - Kd_T_dd(i,K) = 0.7*diff_dd - Kd_S_dd(i,K) = diff_dd - elseif ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then ! diffusive convection - Rrho = alpha_dT/beta_dS - diff_dd = Kv_molecular*0.909*exp(4.6*exp(-0.54*(1/Rrho-1))) - prandtl = 0.15*Rrho - if (Rrho > 0.5) prandtl = (1.85-0.85/Rrho)*Rrho - Kd_T_dd(i,K) = diff_dd - Kd_S_dd(i,K) = prandtl*diff_dd - else - Kd_T_dd(i,K) = 0.0 ; Kd_S_dd(i,K) = 0.0 - endif - enddo - enddo - endif - -end subroutine double_diffusion !> This routine adds diffusion sustained by flow energy extracted by bottom drag. subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & maxTKE, kb, G, GV, CS, Kd, Kd_int, Kd_BBL) @@ -2079,45 +1926,6 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp endif endif - - ! GMM, the following should be moved to the DD module - call get_param(param_file, mdl, "DOUBLE_DIFFUSION", CS%double_diffusion, & - "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.", & - default=2.55, units="nondim") - call get_param(param_file, mdl, "MAX_SALT_DIFF_SALT_FINGERS", CS%Max_salt_diff_salt_fingers, & - "Maximum salt diffusivity for salt fingering regime.", & - default=1.e-4, units="m2 s-1") - call get_param(param_file, mdl, "KV_MOLECULAR", CS%Kv_molecular, & - "Molecular viscosity for calculation of fluxes under \n"//& - "double-diffusive convection.", default=1.5e-6, units="m2 s-1") - ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults. - - CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & - 'Double-diffusive diffusivity for temperature', 'm2 s-1') - - CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & - 'Double-diffusive diffusivity for salinity', 'm2 s-1') - - if (associated(diag_to_Z_CSp)) then - vd = var_desc("KT_extra", "m2 s-1", & - "Double-Diffusive Temperature Diffusivity, interpolated to z", & - z_grid='z') - CS%id_KT_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - vd = var_desc("KS_extra", "m2 s-1", & - "Double-Diffusive Salinity Diffusivity, interpolated to z",& - z_grid='z') - CS%id_KS_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - vd = var_desc("Kd_BBL", "m2 s-1", & - "Bottom Boundary Layer Diffusivity", z_grid='z') - CS%id_Kd_BBL_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) - endif - endif - if (CS%user_change_diff) then call user_change_diff_init(Time, G, param_file, diag, CS%user_change_diff_CSp) endif @@ -2131,7 +1939,10 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp id_clock_kappaShear = cpu_clock_id('(Ocean kappa_shear)', grain=CLOCK_MODULE) ! CVMix shear-driven mixing - CS%use_cvmix_shear = cvmix_shear_init(Time, G, GV, param_file, CS%diag, CS%cvmix_shear_csp) + CS%use_CVMix_shear = CVMix_shear_init(Time, G, GV, param_file, CS%diag, CS%CVMix_shear_csp) + + ! CVMix double diffusion mixing + CS%use_CVMix_ddiff = CVMix_ddiff_init(Time, G, GV, param_file, CS%diag, CS%CVMix_ddiff_csp) end subroutine set_diffusivity_init @@ -2146,8 +1957,11 @@ subroutine set_diffusivity_end(CS) if (CS%user_change_diff) & call user_change_diff_end(CS%user_change_diff_CSp) - if (CS%use_cvmix_shear) & - call cvmix_shear_end(CS%cvmix_shear_csp) + if (CS%use_CVMix_shear) & + call CVMix_shear_end(CS%CVMix_shear_csp) + + if (CS%use_CVMix_ddiff) & + call CVMix_ddiff_end(CS%CVMix_ddiff_csp) if (associated(CS)) deallocate(CS) diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index 18eb80f280..ee18094f7e 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -46,6 +46,7 @@ module MOM_set_visc use MOM_kappa_shear, only : kappa_shear_is_used use MOM_cvmix_shear, only : cvmix_shear_is_used use MOM_cvmix_conv, only : cvmix_conv_is_used +use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used use MOM_io, only : vardesc, var_desc use MOM_restart, only : register_restart_field, MOM_restart_CS use MOM_variables, only : thermo_var_ptrs @@ -1791,7 +1792,7 @@ subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS) call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., & do_not_log=.true.) - use_kappa_shear = .false. ; use_cvmix_shear = .false. ; + use_kappa_shear = .false. ; use_cvmix_shear = .false. useKPP = .false. ; useEPBL = .false. ; use_cvmix_conv = .false. ; if (.not.adiabatic) then use_kappa_shear = kappa_shear_is_used(param_file) @@ -1870,7 +1871,8 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) real :: Kv_background real :: omega_frac_dflt integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, i, j, n - logical :: use_kappa_shear, adiabatic, differential_diffusion, use_omega + logical :: use_kappa_shear, adiabatic, use_omega + logical :: use_CVMix_ddiff type(OBC_segment_type), pointer :: segment ! pointer to OBC segment type ! This include declares and sets the variable "version". #include "version_variable.h" @@ -1893,8 +1895,8 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) ! Set default, read and log parameters call log_version(param_file, mdl, version, "") - CS%RiNo_mix = .false. - use_kappa_shear = .false. ; differential_diffusion = .false. !; adiabatic = .false. ! Needed? -AJA + CS%RiNo_mix = .false. ; use_CVMix_ddiff = .false. + use_kappa_shear = .false. !; adiabatic = .false. ! Needed? -AJA 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"//& @@ -1921,11 +1923,9 @@ 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 + call get_param(param_file, mdl, "PRANDTL_TURB", visc%Prandtl_turb, & "The turbulent Prandtl number applied to shear \n"//& "instability.", units="nondim", default=1.0) @@ -2067,7 +2067,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 (differential_diffusion) then + if (use_CVMix_ddiff) 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 81265e64ae209c0685de5d8b1cfbe22049533b61 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 May 2018 08:09:41 -0600 Subject: [PATCH 21/39] Doxygenize MOM_diabatic_aux --- .../vertical/MOM_diabatic_aux.F90 | 132 +++++++----------- 1 file changed, 52 insertions(+), 80 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 82799d531a..89d16f8e87 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -340,27 +340,23 @@ subroutine differential_diffuse_T_S(h, tv, visc, dt, G, GV) end subroutine differential_diffuse_T_S +!> Keep salinity from falling below a small but positive threshold +!! This occurs when the ice model attempts to extract more salt then +!! is actually available to it from the ocean. subroutine adjust_salt(h, tv, G, GV, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) - type(thermo_var_ptrs), intent(inout) :: tv - type(diabatic_aux_CS), intent(in) :: CS - -! Keep salinity from falling below a small but positive threshold -! This occurs when the ice model attempts to extract more salt then -! is actually available to it from the ocean. - -! Arguments: h - Layer thickness, in m. -! (in/out) tv - A structure containing pointers to any available -! thermodynamic fields. Absent fields have NULL ptrs. -! (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 -! diabatic_driver_init. - real :: salt_add_col(SZI_(G),SZJ_(G)) ! The accumulated salt requirement - real :: S_min ! The minimum salinity - real :: mc ! A layer's mass kg m-2 . + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m + !! or kg m-2) + type(thermo_var_ptrs), intent(inout) :: tv !< structure containing pointers to any + !! available thermodynamic fields. + type(diabatic_aux_CS), intent(in) :: CS !< control structure returned by + !! a previous call to diabatic_driver_init. + + ! local variables + real :: salt_add_col(SZI_(G),SZJ_(G)) !< The accumulated salt requirement + real :: S_min !< The minimum salinity + real :: mc !< A layer's mass kg m-2 . integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = G%ke @@ -402,33 +398,29 @@ subroutine adjust_salt(h, tv, G, GV, CS) end subroutine adjust_salt +!> Insert salt from brine rejection into the first layer below +!! the mixed layer which both contains mass and in which the +!! change in layer density remains stable after the addition +!! of salt via brine rejection. subroutine insert_brine(h, tv, G, GV, fluxes, nkmb, CS, dt, id_brine_lay) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) - type(thermo_var_ptrs), intent(inout) :: tv - type(forcing), intent(in) :: fluxes - integer, intent(in) :: nkmb - type(diabatic_aux_CS), intent(in) :: CS - real, intent(in) :: dt + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m + !! or kg m-2) + type(thermo_var_ptrs), intent(inout) :: tv !< structure containing pointers to + !! any available hermodynamic fields. + type(forcing), intent(in) :: fluxes !< tructure containing pointers + !! any possible forcing fields + integer, intent(in) :: nkmb !< number of layers in the mixed and + !! buffer layers + type(diabatic_aux_CS), intent(in) :: CS !< control structure returned by a + !! previous call to diabatic_driver_init. + real, intent(in) :: dt !< time step between calls to this + !! function (s) ?? integer, intent(in) :: id_brine_lay -! Insert salt from brine rejection into the first layer below -! the mixed layer which both contains mass and in which the -! change in layer density remains stable after the addition -! of salt via brine rejection. - -! Arguments: h - Layer thickness, in m. -! (in/out) tv - A structure containing pointers to any available -! thermodynamic fields. Absent fields have NULL ptrs. -! (in) fluxes = A structure containing pointers to any possible -! forcing fields; unused fields have NULL ptrs. -! (in) nkmb - The number of layers in the mixed and buffer layers. -! (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 -! diabatic_driver_init. + ! local variables real :: salt(SZI_(G)) ! The amount of salt rejected from ! sea ice. [grams] real :: dzbr(SZI_(G)) ! cumulative depth over which brine is distributed @@ -531,10 +523,9 @@ subroutine insert_brine(h, tv, G, GV, fluxes, nkmb, CS, dt, id_brine_lay) end subroutine insert_brine +!> Simple tri-diagnonal solver for T and S. +!! "Simple" means it only uses arrays hold, ea and eb. subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S) -! Simple tri-diagnonal solver for T and S -! "Simple" means it only uses arrays hold, ea and eb - ! Arguments type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure integer, intent(in) :: is, ie, js, je @@ -571,35 +562,22 @@ subroutine triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, T, S) enddo end subroutine triDiagTS - +!> Calculates u_h and v_h (velocities at thickness points), +!! optionally using the entrainments (in m) passed in as arguments. subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, ea, eb) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< The zonal velocity, in m s-1 real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< The meridional velocity, in m s-1 real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: u_h, v_h - real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in), optional :: ea, eb -! This subroutine calculates u_h and v_h (velocities at thickness -! points), optionally using the entrainments (in m) passed in as arguments. - -! Arguments: u - Zonal velocity, in m s-1. -! (in) v - Meridional velocity, in m s-1. -! (in) h - Layer thickness, in m or kg m-2. -! (out) u_h - The zonal velocity at thickness points after -! entrainment, in m s-1. -! (out) v_h - The meridional velocity at thickness points after -! entrainment, in m s-1. -! (in) G - The ocean's grid structure. -! (in) GV - The ocean's vertical grid structure. -! (in, opt) ea - The amount of fluid entrained from the layer above within -! this time step, in units of m or kg m-2. Omitting ea is the -! same as setting it to 0. -! (in, opt) eb - The amount of fluid entrained from the layer below within -! this time step, in units of m or kg m-2. Omitting eb is the -! same as setting it to 0. ea and eb must either be both -! present or both absent. - + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: u_h, v_h !< zonal and meridional velocity at thickness + !! points entrainment, in m s-1. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in), optional :: ea, eb !< The amount of fluid entrained + !! from the layer above within this time step + !! , in units of m or kg m-2. Omitting ea is the + !! same as setting it to 0. + + ! local variables real :: b_denom_1 ! The first term in the denominator of b1 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 or kg m-2. @@ -1306,26 +1284,20 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, dt, fluxes, optics, h, tv, & end subroutine applyBoundaryFluxesInOut +!> Initializes this module. subroutine diabatic_aux_init(Time, G, GV, param_file, diag, CS, useALEalgorithm, use_ePBL) type(time_type), intent(in) :: Time type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters - type(diag_ctrl), target, intent(inout) :: diag - type(diabatic_aux_CS), pointer :: CS - logical, intent(in) :: useALEalgorithm - logical, intent(in) :: use_ePBL - -! Arguments: -! (in) Time = current model time -! (in) G = ocean grid structure -! (in) GV - The ocean's vertical grid structure. -! (in) param_file = structure indicating the open file to parse for parameter values -! (in) diag = structure used to regulate diagnostic output -! (in/out) CS = pointer set to point to the control structure for this module -! (in) use_ePBL = If true, use the implicit energetics planetary boundary -! layer scheme to determine the diffusivity in the -! surface boundary layer. + type(diag_ctrl), target, intent(inout) :: diag !< structure used to regulate diagnostic output + type(diabatic_aux_CS), pointer :: CS !< pointer set to point to the ontrol structure for + !! this module + logical, intent(in) :: useALEalgorithm !< If True, uses ALE. + logical, intent(in) :: use_ePBL !< If true, use the implicit energetics + !! planetary boundary layer scheme to determine the + !! diffusivity in the surface boundary layer. + ! local variables type(vardesc) :: vd ! This "include" declares and sets the variable "version". From 05daededb4103ecd3a3c1fbd2fcefd13968f1278 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 May 2018 08:10:36 -0600 Subject: [PATCH 22/39] Fix indentation --- src/parameterizations/vertical/MOM_set_diffusivity.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index f33bdea2a8..37696386be 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -465,8 +465,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%useKappaShear) call hchksum(visc%Kd_shear,"Turbulent Kd",G%HI,haloshift=0) if (CS%use_CVMix_ddiff) then - call hchksum(visc%Kd_extra_T, "MOM_set_diffusivity: Kd_extra_T",G%HI,haloshift=0) - call hchksum(visc%Kd_extra_S, "MOM_set_diffusivity: Kd_extra_S",G%HI,haloshift=0) + call hchksum(visc%Kd_extra_T, "MOM_set_diffusivity: Kd_extra_T",G%HI,haloshift=0) + call hchksum(visc%Kd_extra_S, "MOM_set_diffusivity: Kd_extra_S",G%HI,haloshift=0) endif if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) then From d5ce7a4aa4d94a239e19d5b5330bcc41ac4a0ae3 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 May 2018 13:54:27 -0600 Subject: [PATCH 23/39] Avoid NaNs when computing stratification parameter --- src/parameterizations/vertical/MOM_cvmix_ddiff.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/parameterizations/vertical/MOM_cvmix_ddiff.F90 b/src/parameterizations/vertical/MOM_cvmix_ddiff.F90 index 8e52c39849..da75caf1e3 100644 --- a/src/parameterizations/vertical/MOM_cvmix_ddiff.F90 +++ b/src/parameterizations/vertical/MOM_cvmix_ddiff.F90 @@ -250,6 +250,8 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, j, Kd_T, Kd_S, CS) if (CS%id_R_rho > 0.0) then do k=1,G%ke CS%R_rho(i,j,k) = alpha_dT(k)/beta_dS(k) + ! avoid NaN's + if(CS%R_rho(i,j,k) /= CS%R_rho(i,j,k)) CS%R_rho(i,j,k) = 0.0 enddo endif From abd620c8c730143aae2bcd55c752741a1d602789 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 May 2018 13:55:19 -0600 Subject: [PATCH 24/39] Move description to the end of the module --- .../vertical/MOM_diabatic_aux.F90 | 93 +++++++++---------- 1 file changed, 45 insertions(+), 48 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 89d16f8e87..fa0cca8681 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -2,53 +2,6 @@ module MOM_diabatic_aux ! This file is part of MOM6. See LICENSE.md for the license. -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg, April 1994 - July 2000 * -!* Alistair Adcroft, and Stephen Griffies * -!* * -!* This program contains the subroutine that, along with the * -!* subroutines that it calls, implements diapycnal mass and momentum * -!* fluxes and a bulk mixed layer. The diapycnal diffusion can be * -!* used without the bulk mixed layer. * -!* * -!* diabatic first determines the (diffusive) diapycnal mass fluxes * -!* based on the convergence of the buoyancy fluxes within each layer. * -!* The dual-stream entrainment scheme of MacDougall and Dewar (JPO, * -!* 1997) is used for combined diapycnal advection and diffusion, * -!* calculated implicitly and potentially with the Richardson number * -!* dependent mixing, as described by Hallberg (MWR, 2000). Diapycnal * -!* advection is fundamentally the residual of diapycnal diffusion, * -!* so the fully implicit upwind differencing scheme that is used is * -!* entirely appropriate. The downward buoyancy flux in each layer * -!* is determined from an implicit calculation based on the previously * -!* calculated flux of the layer above and an estimated flux in the * -!* layer below. This flux is subject to the following conditions: * -!* (1) the flux in the top and bottom layers are set by the boundary * -!* conditions, and (2) no layer may be driven below an Angstrom thick-* -!* ness. If there is a bulk mixed layer, the buffer layer is treat- * -!* ed as a fixed density layer with vanishingly small diffusivity. * -!* * -!* diabatic takes 5 arguments: the two velocities (u and v), the * -!* thicknesses (h), a structure containing the forcing fields, and * -!* the length of time over which to act (dt). The velocities and * -!* thickness are taken as inputs and modified within the subroutine. * -!* There is no limit on the time step. * -!* * -!* 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, T, S, buoy, ustar, ea, eb, etc. * -!* 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_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end use MOM_cpu_clock, only : CLOCK_MODULE_DRIVER, CLOCK_MODULE, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr @@ -251,6 +204,7 @@ subroutine differential_diffuse_T_S(h, tv, visc, dt, G, GV) !! layer properies, and related fields. real, intent(in) :: dt !< Time increment, in s. + ! local variables real, dimension(SZI_(G)) :: & b1_T, b1_S, & ! Variables used by the tridiagonal solvers of T & S, in H. d1_T, d1_S ! Variables used by the tridiagonal solvers, nondim. @@ -337,7 +291,6 @@ subroutine differential_diffuse_T_S(h, tv, visc, dt, G, GV) S(i,j,k) = S(i,j,k) + c1_S(i,k+1)*S(i,j,k+1) enddo ; enddo enddo - end subroutine differential_diffuse_T_S !> Keep salinity from falling below a small but positive threshold @@ -1420,4 +1373,48 @@ subroutine diabatic_aux_end(CS) end subroutine diabatic_aux_end +!> \namespace MOM_diabatic_aux +!! +!! This module contains the subroutines that, along with the * +!! subroutines that it calls, implements diapycnal mass and momentum * +!! fluxes and a bulk mixed layer. The diapycnal diffusion can be * +!! used without the bulk mixed layer. * +!! * +!! diabatic first determines the (diffusive) diapycnal mass fluxes * +!! based on the convergence of the buoyancy fluxes within each layer. * +!! The dual-stream entrainment scheme of MacDougall and Dewar (JPO, * +!! 1997) is used for combined diapycnal advection and diffusion, * +!! calculated implicitly and potentially with the Richardson number * +!! dependent mixing, as described by Hallberg (MWR, 2000). Diapycnal * +!! advection is fundamentally the residual of diapycnal diffusion, * +!! so the fully implicit upwind differencing scheme that is used is * +!! entirely appropriate. The downward buoyancy flux in each layer * +!! is determined from an implicit calculation based on the previously * +!! calculated flux of the layer above and an estimated flux in the * +!! layer below. This flux is subject to the following conditions: * +!! (1) the flux in the top and bottom layers are set by the boundary * +!! conditions, and (2) no layer may be driven below an Angstrom thick-* +!! ness. If there is a bulk mixed layer, the buffer layer is treat- * +!! ed as a fixed density layer with vanishingly small diffusivity. * +!! * +!! diabatic takes 5 arguments: the two velocities (u and v), the * +!! thicknesses (h), a structure containing the forcing fields, and * +!! the length of time over which to act (dt). The velocities and * +!! thickness are taken as inputs and modified within the subroutine. * +!! There is no limit on the time step. * +!! * +!! 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, T, S, buoy, ustar, ea, eb, etc. * +!! 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). * +!! * +!!********+*********+*********+*********+*********+*********+*********+** + end module MOM_diabatic_aux From 152a7073db62e9b9b68023d2eb1dd8529414a7fd Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 May 2018 13:56:00 -0600 Subject: [PATCH 25/39] Change cvmix to CVMix --- .../vertical/MOM_diabatic_driver.F90 | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index b109d3642a..abc0d664f5 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -23,8 +23,8 @@ module MOM_diabatic_driver use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag, calc_Zint_diags use MOM_diapyc_energy_req, only : diapyc_energy_req_init, diapyc_energy_req_end use MOM_diapyc_energy_req, only : diapyc_energy_req_calc, diapyc_energy_req_test, diapyc_energy_req_CS -use MOM_cvmix_conv, only : cvmix_conv_init, cvmix_conv_cs -use MOM_cvmix_conv, only : cvmix_conv_end, calculate_cvmix_conv +use MOM_CVMix_conv, only : CVMix_conv_init, CVMix_conv_cs +use MOM_CVMix_conv, only : CVMix_conv_end, calculate_CVMix_conv use MOM_domains, only : pass_var, To_West, To_South, To_All, Omit_Corners use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_tidal_mixing, only : tidal_mixing_init, tidal_mixing_cs @@ -95,7 +95,7 @@ module MOM_diabatic_driver !! shear-driven diapycnal diffusivity. logical :: use_CVMix_ddiff !< If true, use the CVMix double diffusion module. logical :: use_tidal_mixing !< If true, activate tidal mixing diffusivity. - logical :: use_cvmix_conv !< If true, use the CVMix module to get enhanced + logical :: use_CVMix_conv !< If true, use the CVMix module to get enhanced !! mixing due to convection. logical :: use_sponge !< If true, sponges may be applied anywhere in the !! domain. The exact location and properties of @@ -226,7 +226,7 @@ module MOM_diabatic_driver type(diag_to_Z_CS), pointer :: diag_to_Z_CSp => NULL() type(KPP_CS), pointer :: KPP_CSp => NULL() type(tidal_mixing_cs), pointer :: tidal_mixing_csp => NULL() - type(cvmix_conv_cs), pointer :: cvmix_conv_csp => NULL() + type(CVMix_conv_cs), pointer :: CVMix_conv_csp => NULL() type(diapyc_energy_req_CS), pointer :: diapyc_en_rec_CSp => NULL() type(group_pass_type) :: pass_hold_eb_ea !< For group halo pass @@ -530,7 +530,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G if (CS%debug) then call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, haloshift=0) endif - if (CS%use_kappa_shear .or. CS%use_cvmix_shear) then + if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then call find_uv_at_h(u, v, h_orig, u_h, v_h, G, GV, eaml, ebml) if (CS%debug) then @@ -680,13 +680,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G 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) + 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) + 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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -2349,9 +2349,9 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, ! CS%use_tidal_mixing is set to True if an internal tidal dissipation scheme is to be used. CS%use_tidal_mixing = tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, CS%tidal_mixing_CSp) - ! CS%use_cvmix_conv is set to True if CVMix convection will be used, otherwise + ! CS%use_CVMix_conv is set to True if CVMix convection will be used, otherwise ! False. - CS%use_cvmix_conv = cvmix_conv_init(Time, G, GV, param_file, diag, CS%cvmix_conv_csp) + CS%use_CVMix_conv = CVMix_conv_init(Time, G, GV, param_file, diag, CS%CVMix_conv_csp) call entrain_diffusive_init(Time, G, GV, param_file, diag, CS%entrain_diffusive_CSp) @@ -2442,7 +2442,7 @@ subroutine diabatic_driver_end(CS) if (CS%use_tidal_mixing) call tidal_mixing_end(CS%tidal_mixing_CSp) - if (CS%use_cvmix_conv) call cvmix_conv_end(CS%cvmix_conv_csp) + if (CS%use_CVMix_conv) call CVMix_conv_end(CS%CVMix_conv_csp) if (CS%use_energetic_PBL) & call energetic_PBL_end(CS%energetic_PBL_CSp) From cc273629b858ddc197e2f6f43e75731d87cbac3f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Fri, 18 May 2018 15:01:08 -0600 Subject: [PATCH 26/39] Clean up spaces and comments --- .../vertical/MOM_diabatic_driver.F90 | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index abc0d664f5..27ff1de4d9 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -484,13 +484,13 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G endif if (CS%ML_mix_first > 0.0) then -! This subroutine -! (1) Cools the mixed layer. -! (2) Performs convective adjustment by mixed layer entrainment. -! (3) Heats the mixed layer and causes it to detrain to -! Monin-Obukhov depth or minimum mixed layer depth. -! (4) Uses any remaining TKE to drive mixed layer entrainment. -! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) + ! This subroutine: + ! (1) Cools the mixed layer. + ! (2) Performs convective adjustment by mixed layer entrainment. + ! (3) Heats the mixed layer and causes it to detrain to + ! Monin-Obukhov depth or minimum mixed layer depth. + ! (4) Uses any remaining TKE to drive mixed layer entrainment. + ! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) call find_uv_at_h(u, v, h, u_h, v_h, G, GV) call cpu_clock_begin(id_clock_mixedlayer) @@ -525,11 +525,12 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G if (showCallTree) call callTree_waypoint("done with 1st bulkmixedlayer (diabatic)") if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, G) endif - endif + endif ! end CS%bulkmixedlayer if (CS%debug) then call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, haloshift=0) endif + if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then call find_uv_at_h(u, v, h_orig, u_h, v_h, G, GV, eaml, ebml) @@ -586,7 +587,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G CS%int_tide_input%tideamp, CS%int_tide_input%Nb, dt, G, GV, CS%int_tide_CSp) endif if (showCallTree) call callTree_waypoint("done with propagate_int_tide (diabatic)") - endif + endif ! end CS%use_int_tides call cpu_clock_begin(id_clock_set_diffusivity) ! Sets: Kd, Kd_int, visc%Kd_extra_T, visc%Kd_extra_S From 70d88e4e52f9ef500b43583a2888d9d6acd41209 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Mon, 21 May 2018 17:05:14 -0600 Subject: [PATCH 27/39] Add a flag to control if visc%Kv_slow is used This commit adds a flag ADD_KV_SLOW (default is FALSE) that controls if the background vertical viscosity in the interior (i.e., tidal + background + shear + convenction) is addded when computing the coupling coefficient. The purpose of this flag is to be able to recover previous answers and it will likely be removed in the future since this option should always be true. --- src/core/MOM_variables.F90 | 3 + .../vertical/MOM_set_viscosity.F90 | 86 ++++++++++--------- .../vertical/MOM_vert_friction.F90 | 2 +- 3 files changed, 48 insertions(+), 43 deletions(-) diff --git a/src/core/MOM_variables.F90 b/src/core/MOM_variables.F90 index 09305eb9fb..02b0b622a3 100644 --- a/src/core/MOM_variables.F90 +++ b/src/core/MOM_variables.F90 @@ -233,6 +233,9 @@ module MOM_variables !! convection etc). TKE_turb => NULL() !< The turbulent kinetic energy per unit mass defined !! at the interfaces between each layer, in m2 s-2. + logical :: add_Kv_slow !< If True, adds Kv_slow when calculating the + !! 'coupling coefficient' (a[k]) at the interfaces. + !! This is done in find_coupling_coef. end type vertvisc_type !> The BT_cont_type structure contains information about the summed layer diff --git a/src/parameterizations/vertical/MOM_set_viscosity.F90 b/src/parameterizations/vertical/MOM_set_viscosity.F90 index c0b03e5832..ec1b09a5ad 100644 --- a/src/parameterizations/vertical/MOM_set_viscosity.F90 +++ b/src/parameterizations/vertical/MOM_set_viscosity.F90 @@ -2,38 +2,6 @@ module MOM_set_visc ! This file is part of MOM6. See LICENSE.md for the license. -!********+*********+*********+*********+*********+*********+*********+** -!* * -!* By Robert Hallberg, April 1994 - October 2006 * -!* Quadratic Bottom Drag by James Stephens and R. Hallberg. * -!* * -!* This file contains the subroutine that calculates various values * -!* related to the bottom boundary layer, such as the viscosity and * -!* thickness of the BBL (set_viscous_BBL). This would also be the * -!* module in which other viscous quantities that are flow-independent * -!* might be set. This information is transmitted to other modules * -!* via a vertvisc type structure. * -!* * -!* The same code is used for the two velocity components, by * -!* indirectly referencing the velocities and defining a handful of * -!* direction-specific defined variables. * -!* * -!* 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, frhatv, tauy * -!* j x ^ x ^ x At >: u, frhatu, taux * -!* j > o > o > At o: h * -!* 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_debugging, only : uvchksum, hchksum use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr @@ -1859,16 +1827,8 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) type(set_visc_CS), pointer :: CS !< A pointer that is set to point to the control !! structure for this module type(ocean_OBC_type), pointer :: OBC -! Arguments: Time - The current model time. -! (in) G - The ocean's grid 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) diag - A structure that is used to regulate diagnostic output. -! (out) visc - A structure containing vertical viscosities and related -! fields. Allocated here. -! (in/out) CS - A pointer that is set to point to the control structure -! for this module + + ! local variables real :: Csmag_chan_dflt, smag_const1, TKE_decay_dflt, bulk_Ri_ML_dflt real :: Kv_background real :: omega_frac_dflt @@ -2020,6 +1980,15 @@ subroutine set_visc_init(Time, G, GV, param_file, diag, visc, CS, OBC) "The background kinematic viscosity in the interior. \n"//& "The molecular value, ~1e-6 m2 s-1, may be used.", & units="m2 s-1", fail_if_missing=.true.) + + call get_param(param_file, mdl, "ADD_KV_SLOW", visc%add_Kv_slow, & + "If true, the background vertical viscosity in the interior \n"//& + "(i.e., tidal + background + shear + convenction) is addded \n"// & + "when computing the coupling coefficient. The purpose of this \n"// & + "flag is to be able to recover previous answers and it will likely \n"// & + "be removed in the future since this option should always be true.", & + default=.false.) + call get_param(param_file, mdl, "KV_BBL_MIN", CS%KV_BBL_min, & "The minimum viscosities in the bottom boundary layer.", & units="m2 s-1", default=Kv_background) @@ -2117,4 +2086,37 @@ subroutine set_visc_end(visc, CS) deallocate(CS) end subroutine set_visc_end +!> \namespace MOM_set_visc +!!********+*********+*********+*********+*********+*********+*********+** +!!* * +!!* By Robert Hallberg, April 1994 - October 2006 * +!!* Quadratic Bottom Drag by James Stephens and R. Hallberg. * +!!* * +!!* This file contains the subroutine that calculates various values * +!!* related to the bottom boundary layer, such as the viscosity and * +!!* thickness of the BBL (set_viscous_BBL). This would also be the * +!!* module in which other viscous quantities that are flow-independent * +!!* might be set. This information is transmitted to other modules * +!!* via a vertvisc type structure. * +!!* * +!!* The same code is used for the two velocity components, by * +!!* indirectly referencing the velocities and defining a handful of * +!!* direction-specific defined variables. * +!!* * +!!* 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, frhatv, tauy * +!!* j x ^ x ^ x At >: u, frhatu, taux * +!!* j > o > o > At o: h * +!!* 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). * +!!* * +!!********+*********+*********+*********+*********+*********+*********+** + end module MOM_set_visc diff --git a/src/parameterizations/vertical/MOM_vert_friction.F90 b/src/parameterizations/vertical/MOM_vert_friction.F90 index ad0d7fc90d..bafbe5eb59 100644 --- a/src/parameterizations/vertical/MOM_vert_friction.F90 +++ b/src/parameterizations/vertical/MOM_vert_friction.F90 @@ -1195,7 +1195,7 @@ subroutine find_coupling_coef(a, hvel, do_i, h_harm, bbl_thick, kv_bbl, z_i, h_m endif ! add "slow" varying vertical viscosity (e.g., from background, tidal etc) - if (associated(visc%Kv_slow)) then + if (associated(visc%Kv_slow) .and. (visc%add_Kv_slow)) then ! GMM/ A factor of 2 is also needed here, see comment above from BGR. if (work_on_u) then do K=2,nz ; do i=is,ie ; if (do_i(i)) then From 8877c17dccd111e4789f4553eb4e25befee6e091 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 22 May 2018 08:32:43 -0600 Subject: [PATCH 28/39] Rename MOM_cvmix_ddiff.F90 -> MOM_CVMix_ddiff.F90 --- .../vertical/{MOM_cvmix_ddiff.F90 => MOM_CVMix_ddiff.F90} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/parameterizations/vertical/{MOM_cvmix_ddiff.F90 => MOM_CVMix_ddiff.F90} (100%) diff --git a/src/parameterizations/vertical/MOM_cvmix_ddiff.F90 b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 similarity index 100% rename from src/parameterizations/vertical/MOM_cvmix_ddiff.F90 rename to src/parameterizations/vertical/MOM_CVMix_ddiff.F90 From bf6c0030a308c3099fce3674a84a03ad16e358ba Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 22 May 2018 08:41:39 -0600 Subject: [PATCH 29/39] Clean the ddiff code and improve comments --- .../vertical/MOM_CVMix_ddiff.F90 | 27 +++++-------------- 1 file changed, 6 insertions(+), 21 deletions(-) diff --git a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 index da75caf1e3..7137aabfa6 100644 --- a/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 +++ b/src/parameterizations/vertical/MOM_CVMix_ddiff.F90 @@ -64,9 +64,6 @@ logical function CVMix_ddiff_init(Time, G, GV, param_file, diag, CS) type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure. type(CVMix_ddiff_cs), pointer :: CS !< This module's control structure. - ! Local variables -! real :: prandtl_conv !< Turbulent Prandtl number used in convective instabilities. - ! This include declares and sets the variable "version". #include "version_variable.h" @@ -133,11 +130,6 @@ logical function CVMix_ddiff_init(Time, G, GV, param_file, diag, CS) call closeParameterBlock(param_file) - ! allocate arrays and set them to zero - ! GMM, dont need the following - !allocate(CS%KT_extra(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%KT_extra(:,:,:) = 0. - !allocate(CS%KS_extra(SZI_(G), SZJ_(G), SZK_(G)+1)); CS%KS_extra(:,:,:) = 0. - ! Register diagnostics CS%diag => diag @@ -173,8 +165,6 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, j, Kd_T, Kd_S, CS) real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness, in m or kg m-2. type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure. real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: Kd_T !< Interface double diffusion diapycnal -! real, dimension(:,:,:), pointer :: Kd_T -! real, dimension(:,:,:), pointer :: Kd_S !! diffusivity for temp (m2/sec). real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: Kd_S !< Interface double diffusion diapycnal !! diffusivity for salt (m2/sec). @@ -209,7 +199,9 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, j, Kd_T, Kd_S, CS) ! set Kd_T and Kd_S to zero to avoid passing values from previous call Kd_T(:,j,:) = 0.0; Kd_S(:,j,:) = 0.0 - ! GMM, check this. + ! GMM, I am leaving some code commented below. We need to pass BLD to + ! this soubroutine to avoid adding diffusivity above that. This needs + ! to be done once we re-structure the order of the calls. !if (.not. associated(hbl)) then ! allocate(hbl(SZI_(G), SZJ_(G))); ! hbl(:,:) = 0.0 @@ -239,9 +231,9 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, j, Kd_T, Kd_S, CS) call calculate_density_derivs(temp_int(:), salt_int(:), pres_int(:), drho_dT(:), drho_dS(:), 1, G%ke, TV%EQN_OF_STATE) - ! GMM, explain need for -1 in front of alpha - ! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then ! salt finger case - ! if ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then ! diffusive convection + ! The "-1.0" below is needed so that the following criteria is satisfied: + ! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then "salt finger" + ! if ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then "diffusive convection" do k=1,G%ke alpha_dT(k) = -1.0*drho_dT(k) * dT(k) beta_dS(k) = drho_dS(k) * dS(k) @@ -277,13 +269,6 @@ subroutine compute_ddiff_coeffs(h, tv, G, GV, j, Kd_T, Kd_S, CS) nlev=G%ke, & max_nlev=G%ke) - !if (is_root_pe()) then - ! write(*,*)'drho_dT, alpha_dT', & - ! drho_dT(:), alpha_dT(:) - ! write(*,*)'drho_dS, beta_dS', & - ! drho_dS(:), beta_dS(:) - !endif - ! Do not apply mixing due to convection within the boundary layer !do k=1,kOBL ! Kd_T(i,j,k) = 0.0 From 198d75559a5f7ef65c566a3f56f9c641286329e8 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 23 May 2018 11:14:06 -0600 Subject: [PATCH 30/39] Delete visc%kv_slow=0 since this is done in set_diffusivity --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 882ed8cb26..698243a7f6 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -387,9 +387,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect*h_neglect Kd_heat(:,:,:) = 0.0 ; Kd_salt(:,:,:) = 0.0 - ! visc%Kv_slow must be set to zero - visc%Kv_slow(:,:,:) = 0.0 - if (nz == 1) return showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("diabatic(), MOM_diabatic_driver.F90") From 720dbc04da0cd278d6c2fe6f56aa2edc0c14e90b Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 23 May 2018 11:25:31 -0600 Subject: [PATCH 31/39] Doxygenize set_diff + read background kinematic viscosity --- .../vertical/MOM_set_diffusivity.F90 | 201 ++++++++++-------- 1 file changed, 109 insertions(+), 92 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 9835c19912..903868795a 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -45,95 +45,94 @@ module MOM_set_diffusivity public set_diffusivity_end type, public :: set_diffusivity_CS ; private - logical :: debug ! If true, write verbose checksums for debugging. - - logical :: bulkmixedlayer ! If true, a refined bulk mixed layer is used with - ! GV%nk_rho_varies variable density mixed & buffer - ! layers. - real :: FluxRi_max ! The flux Richardson number where the stratification is - ! large enough that N2 > omega2. The full expression for - ! the Flux Richardson number is usually - ! FLUX_RI_MAX*N2/(N2+OMEGA2). The default is 0.2. - logical :: bottomdraglaw ! If true, the bottom stress is calculated with a - ! drag law c_drag*|u|*u. - logical :: BBL_mixing_as_max ! If true, take the maximum of the diffusivity - ! from the BBL mixing and the other diffusivities. - ! Otherwise, diffusivities from the BBL_mixing is - ! added. - logical :: use_LOTW_BBL_diffusivity ! If true, use simpler/less precise, BBL diffusivity. - logical :: LOTW_BBL_use_omega ! If true, use simpler/less precise, BBL diffusivity. - real :: BBL_effic ! efficiency with which the energy extracted - ! by bottom drag drives BBL diffusion (nondim) - real :: cdrag ! quadratic drag coefficient (nondim) - real :: IMax_decay ! inverse of a maximum decay scale for - ! bottom-drag driven turbulence, (1/m) - - real :: Kd ! interior diapycnal diffusivity (m2/s) - real :: Kd_min ! minimum diapycnal diffusivity (m2/s) - real :: Kd_max ! maximum increment for diapycnal diffusivity (m2/s) - ! Set to a negative value to have no limit. - real :: Kd_add ! uniform diffusivity added everywhere without - ! filtering or scaling (m2/s) - real :: Kv ! interior vertical viscosity (m2/s) - real :: Kdml ! mixed layer diapycnal diffusivity (m2/s) - ! when bulkmixedlayer==.false. - real :: Hmix ! mixed layer thickness (meter) when - ! bulkmixedlayer==.false. + logical :: debug !< If true, write verbose checksums for debugging. + + logical :: bulkmixedlayer !< If true, a refined bulk mixed layer is used with + !! GV%nk_rho_varies variable density mixed & buffer + !! layers. + real :: FluxRi_max !< The flux Richardson number where the stratification is + !! large enough that N2 > omega2. The full expression for + !! the Flux Richardson number is usually + !! FLUX_RI_MAX*N2/(N2+OMEGA2). The default is 0.2. + logical :: bottomdraglaw !< If true, the bottom stress is calculated with a + !! drag law c_drag*|u|*u. + logical :: BBL_mixing_as_max !< If true, take the maximum of the diffusivity + !! from the BBL mixing and the other diffusivities. + !! Otherwise, diffusivities from the BBL_mixing is + !! added. + logical :: use_LOTW_BBL_diffusivity !< If true, use simpler/less precise, BBL diffusivity. + logical :: LOTW_BBL_use_omega !< If true, use simpler/less precise, BBL diffusivity. + real :: BBL_effic !< efficiency with which the energy extracted + !! by bottom drag drives BBL diffusion (nondim) + real :: cdrag !< quadratic drag coefficient (nondim) + real :: IMax_decay !< inverse of a maximum decay scale for + !! bottom-drag driven turbulence, (1/m) + real :: Kv !< The interior vertical viscosity (m2/s) + real :: Kd !< interior diapycnal diffusivity (m2/s) + real :: Kd_min !< minimum diapycnal diffusivity (m2/s) + real :: Kd_max !< maximum increment for diapycnal diffusivity (m2/s) + !! Set to a negative value to have no limit. + real :: Kd_add !< uniform diffusivity added everywhere without + !! filtering or scaling (m2/s) + real :: Kdml !< mixed layer diapycnal diffusivity (m2/s) + !! when bulkmixedlayer==.false. + real :: Hmix !< mixed layer thickness (meter) when + !! bulkmixedlayer==.false. type(diag_ctrl), pointer :: diag ! structure to regulate diagn output timing - logical :: limit_dissipation ! If enabled, dissipation is limited to be larger - ! than the following: - real :: dissip_min ! Minimum dissipation (W/m3) - real :: dissip_N0 ! Coefficient a in minimum dissipation = a+b*N (W/m3) - real :: dissip_N1 ! Coefficient b in minimum dissipation = a+b*N (J/m3) - real :: dissip_N2 ! Coefficient c in minimum dissipation = c*N2 (W m-3 s2) - real :: dissip_Kd_min ! Minimum Kd (m2/s) with dissipatio Rho0*Kd_min*N^2 - - real :: TKE_itide_max ! maximum internal tide conversion (W m-2) - ! available to mix above the BBL - real :: omega ! Earth's rotation frequency (s-1) - logical :: ML_radiation ! allow a fraction of TKE available from wind work - ! to penetrate below mixed layer base with a vertical - ! decay scale determined by the minimum of - ! (1) The depth of the mixed layer, or - ! (2) An Ekman length scale. - ! Energy availble to drive mixing below the mixed layer is - ! given by E = ML_RAD_COEFF*MSTAR*USTAR**3. Optionally, if - ! ML_rad_TKE_decay is true, this is further reduced by a factor - ! of exp(-h_ML*Idecay_len_TkE), where Idecay_len_TKE is - ! calculated the same way as in the mixed layer code. - ! The diapycnal diffusivity is KD(k) = E/(N2(k)+OMEGA2), - ! where N2 is the squared buoyancy frequency (s-2) and OMEGA2 - ! is the rotation rate of the earth squared. - real :: ML_rad_kd_max ! Maximum diapycnal diffusivity due to turbulence - ! radiated from the base of the mixed layer (m2/s) - real :: ML_rad_efold_coeff ! non-dim coefficient to scale penetration depth - real :: ML_rad_coeff ! coefficient, which scales MSTAR*USTAR^3 to - ! obtain energy available for mixing below - ! mixed layer base (nondimensional) - logical :: ML_rad_TKE_decay ! If true, apply same exponential decay - ! to ML_rad as applied to the other surface - ! sources of TKE in the mixed layer code. - real :: ustar_min ! A minimum value of ustar to avoid numerical - ! problems (m/s). If the value is small enough, - ! this parameter should not affect the solution. - real :: TKE_decay ! ratio of natural Ekman depth to TKE decay scale (nondim) - real :: mstar ! ratio of friction velocity cubed to - ! TKE input to the mixed layer (nondim) - logical :: ML_use_omega ! If true, use absolute rotation rate instead - ! of the vertical component of rotation when - ! setting the decay scale for mixed layer turbulence. - real :: ML_omega_frac ! When setting the decay scale for turbulence, use - ! this fraction of the absolute rotation rate blended - ! with the local value of f, as f^2 ~= (1-of)*f^2 + of*4*omega^2. - logical :: user_change_diff ! If true, call user-defined code to change diffusivity. - logical :: useKappaShear ! If true, use the kappa_shear module to find the - ! shear-driven diapycnal diffusivity. - logical :: use_CVMix_shear ! If true, use one of the CVMix modules to find - ! shear-driven diapycnal diffusivity. - logical :: use_CVMix_ddiff ! If true, enable double-diffusive mixing via CVMix. - logical :: simple_TKE_to_Kd ! If true, uses a simple estimate of Kd/TKE that - ! does not rely on a layer-formulation. + logical :: limit_dissipation !< If enabled, dissipation is limited to be larger + !! than the following: + real :: dissip_min !< Minimum dissipation (W/m3) + real :: dissip_N0 !< Coefficient a in minimum dissipation = a+b*N (W/m3) + real :: dissip_N1 !< Coefficient b in minimum dissipation = a+b*N (J/m3) + real :: dissip_N2 !< Coefficient c in minimum dissipation = c*N2 (W m-3 s2) + real :: dissip_Kd_min !< Minimum Kd (m2/s) with dissipatio Rho0*Kd_min*N^2 + + real :: TKE_itide_max !< maximum internal tide conversion (W m-2) + !! available to mix above the BBL + real :: omega !< Earth's rotation frequency (s-1) + logical :: ML_radiation !< allow a fraction of TKE available from wind work + !! to penetrate below mixed layer base with a vertical + !! decay scale determined by the minimum of + !! (1) The depth of the mixed layer, or + !! (2) An Ekman length scale. + !! Energy availble to drive mixing below the mixed layer is + !! given by E = ML_RAD_COEFF*MSTAR*USTAR**3. Optionally, if + !! ML_rad_TKE_decay is true, this is further reduced by a factor + !! of exp(-h_ML*Idecay_len_TkE), where Idecay_len_TKE is + !! calculated the same way as in the mixed layer code. + !! The diapycnal diffusivity is KD(k) = E/(N2(k)+OMEGA2), + !! where N2 is the squared buoyancy frequency (s-2) and OMEGA2 + !! is the rotation rate of the earth squared. + real :: ML_rad_kd_max !< Maximum diapycnal diffusivity due to turbulence + !! radiated from the base of the mixed layer (m2/s) + real :: ML_rad_efold_coeff !< non-dim coefficient to scale penetration depth + real :: ML_rad_coeff !< coefficient, which scales MSTAR*USTAR^3 to + !! obtain energy available for mixing below + !! mixed layer base (nondimensional) + logical :: ML_rad_TKE_decay !< If true, apply same exponential decay + !! to ML_rad as applied to the other surface + !! sources of TKE in the mixed layer code. + real :: ustar_min !< A minimum value of ustar to avoid numerical + !! problems (m/s). If the value is small enough, + !! this parameter should not affect the solution. + real :: TKE_decay !< ratio of natural Ekman depth to TKE decay scale (nondim) + real :: mstar !! ratio of friction velocity cubed to + !! TKE input to the mixed layer (nondim) + logical :: ML_use_omega !< If true, use absolute rotation rate instead + !! of the vertical component of rotation when + !! setting the decay scale for mixed layer turbulence. + real :: ML_omega_frac !< When setting the decay scale for turbulence, use + !! this fraction of the absolute rotation rate blended + !! with the local value of f, as f^2 ~= (1-of)*f^2 + of*4*omega^2. + logical :: user_change_diff !< If true, call user-defined code to change diffusivity. + logical :: useKappaShear !< If true, use the kappa_shear module to find the + !! shear-driven diapycnal diffusivity. + logical :: use_CVMix_shear !< If true, use one of the CVMix modules to find + !! shear-driven diapycnal diffusivity. + logical :: use_CVMix_ddiff !< If true, enable double-diffusive mixing via CVMix. + logical :: simple_TKE_to_Kd !< If true, uses a simple estimate of Kd/TKE that + !! does not rely on a layer-formulation. character(len=200) :: inputdir type(user_change_diff_CS), pointer :: user_change_diff_CSp => NULL() @@ -177,6 +176,17 @@ module MOM_set_diffusivity contains +!> Sets the interior vertical diffusion of scalars due to the following processes: +!! 1) Shear-driven mixing: two options, Jackson et at. and KPP interior; +!! 2) Background mixing via CVMix (Bryan-Lewis profile) or the scheme described by +!! Harrison & Hallberg, JPO 2008; +!! 3) Double-diffusion aplpied via CVMix; +!! 4) Tidal mixing: many options available, see MOM_tidal_mixing.F90; +!! In addition, this subroutine has the option to set the interior vertical +!! viscosity associated with processes 2-4 listed above, which is stored in +!! visc%Kv_slow. Vertical viscosity due to shear-driven mixing is passed via +!! visc%Kv_shear +!! GMM, TODO: add contribution from tidal mixing into visc%Kv_slow subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & G, GV, CS, Kd, Kd_int) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. @@ -188,9 +198,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: u_h + intent(in) :: u_h !< zonal thickness transport m^2/s. real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & - intent(in) :: v_h + intent(in) :: v_h !< meridional thickness transport m^2/s. type(thermo_var_ptrs), intent(inout) :: tv !< Structure with pointers to thermodynamic !! fields. Out is for tv%TempxPmE. type(forcing), intent(in) :: fluxes !< Structure of surface fluxes that may be @@ -270,7 +280,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! If nothing else is specified, this will be the value used. Kd(:,:,:) = CS%Kd Kd_int(:,:,:) = CS%Kd - visc%Kv_slow(:,:,:) = CS%Kv + if (associated(visc%Kv_slow)) visc%Kv_slow(:,:,:) = CS%Kv ! Set up arrays for diagnostics. @@ -331,6 +341,10 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & elseif (CS%use_CVMix_shear) then !NOTE{BGR}: this needs to be cleaned up. It works in 1D case, but has not been tested outside. call calculate_CVMix_shear(u_h, v_h, h, tv, visc%Kd_shear, visc%Kv_shear,G,GV,CS%CVMix_shear_csp) + if (CS%debug) then + call hchksum(visc%Kd_shear, "after CVMix_shear visc%Kd_shear",G%HI) + call hchksum(visc%Kv_shear, "after CVMix_shear visc%Kv_shear",G%HI) + endif elseif (associated(visc%Kv_shear)) then visc%Kv_shear(:,:,:) = 0. ! needed if calculate_kappa_shear is not enabled endif @@ -342,8 +356,6 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! set surface diffusivities (CS%bkgnd_mixing_csp%Kd_sfc) call sfc_bkgnd_mixing(G, CS%bkgnd_mixing_csp) -! GMM, fix OMP calls below - !$OMP parallel do default(none) shared(is,ie,js,je,nz,G,GV,CS,h,tv,T_f,S_f,fluxes,dd, & !$OMP Kd,visc, & !$OMP Kd_int,dt,u,v,Omega2) & @@ -1826,6 +1838,11 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp ! set params releted to the background mixing call bkgnd_mixing_init(Time, G, GV, param_file, CS%diag, CS%bkgnd_mixing_csp) + call get_param(param_file, mdl, "KV", CS%Kv, & + "The background kinematic viscosity in the interior. \n"//& + "The molecular value, ~1e-6 m2 s-1, may be used.", & + units="m2 s-1", fail_if_missing=.true.) + call get_param(param_file, mdl, "KD", CS%Kd, & "The background diapycnal diffusivity of density in the \n"//& "interior. Zero or the molecular value, ~1e-7 m2 s-1, \n"//& From 6782ecb9dade745c6c5307cbf69b66b3e3a0dc33 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 23 May 2018 17:13:39 -0600 Subject: [PATCH 32/39] Deleted calls related to layer mode --- .../vertical/MOM_diabatic_driver.F90 | 653 ++++-------------- 1 file changed, 146 insertions(+), 507 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 698243a7f6..a26e44fe48 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -279,6 +279,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & type(diabatic_CS), pointer :: CS !< module control structure type(Wave_parameters_CS), optional, pointer :: Waves !< Surface gravity waves + ! local variables real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: & ea, & ! amount of fluid entrained from the layer above within ! one time step (m for Bouss, kg/m^2 for non-Bouss) @@ -391,6 +392,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & showCallTree = callTree_showQuery() if (showCallTree) call callTree_enter("diabatic(), MOM_diabatic_driver.F90") + if (.not. (CS%useALEalgorithm)) call MOM_error(FATAL, "MOM_diabatic_driver: "// & + "The ALE algorithm must be enabled when using MOM_diabatic_driver.") ! Offer diagnostics of various state varables at the start of diabatic ! these are mostly for debugging purposes. @@ -452,7 +455,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%id_frazil_h > 0) call post_data(CS%id_frazil_h, h, CS%diag) endif call disable_averaging(CS%diag) - endif + endif !associated(tv%T) .AND. associated(tv%frazil) + ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep call enable_averaging(dt, Time_end, CS%diag) if (CS%debugConservation) call MOM_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, G) @@ -482,58 +486,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (associated(CS%optics)) & call set_opacity(CS%optics, fluxes, G, GV, CS%opacity_CSp) - if (CS%bulkmixedlayer) then - if (CS%debug) then - call MOM_forcing_chksum("Before mixedlayer", fluxes, G, haloshift=0) - endif - - if (CS%ML_mix_first > 0.0) then - ! This subroutine: - ! (1) Cools the mixed layer. - ! (2) Performs convective adjustment by mixed layer entrainment. - ! (3) Heats the mixed layer and causes it to detrain to - ! Monin-Obukhov depth or minimum mixed layer depth. - ! (4) Uses any remaining TKE to drive mixed layer entrainment. - ! (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate) - call find_uv_at_h(u, v, h, u_h, v_h, G, GV) - - call cpu_clock_begin(id_clock_mixedlayer) - if (CS%ML_mix_first < 1.0) then - ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???) - call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt*CS%ML_mix_first, & - eaml,ebml, G, GV, CS%bulkmixedlayer_CSp, CS%optics, & - Hml, CS%aggregate_FW_forcing, dt, last_call=.false.) - if (CS%salt_reject_below_ML) & - call insert_brine(h, tv, G, GV, fluxes, nkmb, CS%diabatic_aux_CSp, & - dt*CS%ML_mix_first, CS%id_brine_lay) - else - ! Changes: h, tv%T, tv%S, eaml and ebml (G is also inout???) - call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt, eaml, ebml, & - G, GV, CS%bulkmixedlayer_CSp, CS%optics, & - Hml, CS%aggregate_FW_forcing, dt, last_call=.true.) - endif - - ! 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 - ! from this limitation, in which case we can let salinity=0 and still - ! have salt conserved with SIS2 ice. So for SIS2, we can run with - ! BOUND_SALINITY=False in MOM.F90. - if (associated(tv%S) .and. associated(tv%salt_deficit)) & - call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - call cpu_clock_end(id_clock_mixedlayer) - if (CS%debug) then - call MOM_state_chksum("After mixedlayer ", u, v, h, G, GV, haloshift=0) - call MOM_forcing_chksum("After mixedlayer", fluxes, G, haloshift=0) - endif - if (showCallTree) call callTree_waypoint("done with 1st bulkmixedlayer (diabatic)") - if (CS%debugConservation) call MOM_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, G) - endif - endif ! end CS%bulkmixedlayer - - if (CS%debug) then + if (CS%debug) & call MOM_state_chksum("before find_uv_at_h", u, v, h, G, GV, haloshift=0) - endif if (CS%use_kappa_shear .or. CS%use_CVMix_shear) then if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then @@ -609,7 +563,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call hchksum(Kd_Int, "after set_diffusivity Kd_Int",G%HI,haloshift=0) endif - if (CS%useKPP) then call cpu_clock_begin(id_clock_kpp) ! KPP needs the surface buoyancy flux but does not update state variables. @@ -752,47 +705,22 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! This block sets ea, eb from Kd or Kd_int. - ! If using ALE algorithm, set ea=eb=Kd_int on interfaces for - ! use in the tri-diagonal solver. - ! Otherwise, call entrainment_diffusive() which sets ea and eb - ! based on KD and target densities (ie. does remapping as well). - if (CS%useALEalgorithm) then + ! set ea=eb=Kd_int on interfaces for use in the tri-diagonal solver. - do j=js,je ; do i=is,ie - ea(i,j,1) = 0. - enddo ; enddo + do j=js,je ; do i=is,ie + ea(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 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) - enddo ; enddo ; enddo - do j=js,je ; do i=is,ie - eb(i,j,nz) = 0. - enddo ; enddo - if (showCallTree) call callTree_waypoint("done setting ea,eb from Kd_int (diabatic)") - - else ! .not. CS%useALEalgorithm - ! When not using ALE, calculate layer entrainments/detrainments from - ! diffusivities and differences between layer and target densities - call cpu_clock_begin(id_clock_entrain) - ! Calculate appropriately limited diapycnal mass fluxes to account - ! for diapycnal diffusion and advection. Sets: ea, eb. Changes: kb - call Entrainment_diffusive(u, v, h, tv, fluxes, dt, G, GV, CS%entrain_diffusive_CSp, & - ea, eb, kb, Kd_Lay=Kd, Kd_int=Kd_int) - call cpu_clock_end(id_clock_entrain) - if (showCallTree) call callTree_waypoint("done with Entrainment_diffusive (diabatic)") - - endif ! endif for (CS%useALEalgorithm) - - if (CS%debug) then - call MOM_forcing_chksum("after calc_entrain ", fluxes, G, haloshift=0) - call MOM_thermovar_chksum("after calc_entrain ", tv, G) - call MOM_state_chksum("after calc_entrain ", u, v, h, G, GV, haloshift=0) - call hchksum(ea, "after calc_entrain ea", G%HI, haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after calc_entrain eb", G%HI, haloshift=0, scale=GV%H_to_m) - endif + 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) + enddo ; enddo ; enddo + do j=js,je ; do i=is,ie + eb(i,j,nz) = 0. + enddo ; enddo + if (showCallTree) call callTree_waypoint("done setting ea,eb from Kd_int (diabatic)") ! Save fields before boundary forcing is applied for tendency diagnostics if (CS%boundary_forcing_tendency_diag) then @@ -803,96 +731,90 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo ; enddo ; enddo endif - ! Apply forcing when using the ALE algorithm - if (CS%useALEalgorithm) then - call cpu_clock_begin(id_clock_remap) - - ! Changes made to following fields: h, tv%T and tv%S. - - do k=1,nz ; do j=js,je ; do i=is,ie - h_prebound(i,j,k) = h(i,j,k) - enddo ; enddo ; enddo - if (CS%use_energetic_PBL) then - - skinbuoyflux(:,:) = 0.0 - call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & - h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & - 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(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) - endif + ! Apply forcing + call cpu_clock_begin(id_clock_remap) - call find_uv_at_h(u, v, h, u_h, v_h, G, GV) - call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, & - CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) - - ! If visc%MLD exists, copy the ePBL's MLD into it - if (associated(visc%MLD)) then - call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G) - call pass_var(visc%MLD, G%domain, halo=1) - Hml(:,:) = visc%MLD(:,:) - endif + ! Changes made to following fields: h, tv%T and tv%S. + do k=1,nz ; do j=js,je ; do i=is,ie + h_prebound(i,j,k) = h(i,j,k) + enddo ; enddo ; enddo + if (CS%use_energetic_PBL) then - ! Augment the diffusivities due to those diagnosed in energetic_PBL. - do K=2,nz ; do j=js,je ; do i=is,ie + skinbuoyflux(:,:) = 0.0 + call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & + h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, & + CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux) - if (CS%ePBL_is_additive) then - Kd_add_here = Kd_ePBL(i,j,K) - visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + Kd_ePBL(i,j,K) - else - 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(i,j,k-1) = eb(i,j,k-1) + Ent_int - ea(i,j,k) = ea(i,j,k) + Ent_int - Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here + 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(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) + endif - ! 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) + call find_uv_at_h(u, v, h, u_h, v_h, G, GV) + call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, & + CS%energetic_PBL_CSp, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) - enddo ; enddo ; enddo + ! If visc%MLD exists, copy the ePBL's MLD into it + if (associated(visc%MLD)) then + call energetic_PBL_get_MLD(CS%energetic_PBL_CSp, visc%MLD, G) + call pass_var(visc%MLD, G%domain, halo=1) + Hml(:,:) = visc%MLD(:,:) + endif - 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(Kd_ePBL, "after ePBL Kd_ePBL",G%HI,haloshift=0) + ! Augment the diffusivities 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) + visc%Kv_shear(i,j,K) = visc%Kv_shear(i,j,K) + Kd_ePBL(i,j,K) + else + 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(i,j,k-1) = eb(i,j,k-1) + Ent_int + ea(i,j,k) = ea(i,j,k) + Ent_int + Kd_int(i,j,K) = Kd_int(i,j,K) + Kd_add_here - else - call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & - h, tv, CS%aggregate_FW_forcing, & - CS%evap_CFL_limit, CS%minimum_forcing_depth) - - endif ! endif for CS%use_energetic_PBL - - ! diagnose the tendencies due to boundary forcing - ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme - ! so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards - if (CS%boundary_forcing_tendency_diag) then - call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, G, GV, CS) - if (CS%id_boundary_forcing_h > 0) call post_data(CS%id_boundary_forcing_h, h, CS%diag, alt_h = h_diag) - endif - ! Boundary fluxes may have changed T, S, and h - call diag_update_remap_grids(CS%diag) + ! 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) + + enddo ; enddo ; enddo - call cpu_clock_end(id_clock_remap) if (CS%debug) then - call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, haloshift=0) - call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G) - call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, haloshift=0) + 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(Kd_ePBL, "after ePBL Kd_ePBL",G%HI,haloshift=0) endif - 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) - endif ! endif for (CS%useALEalgorithm) + else + call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, dt, fluxes, CS%optics, & + h, tv, CS%aggregate_FW_forcing, & + CS%evap_CFL_limit, CS%minimum_forcing_depth) + + endif ! endif for CS%use_energetic_PBL + + ! diagnose the tendencies due to boundary forcing + ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme + ! so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards + if (CS%boundary_forcing_tendency_diag) then + call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, G, GV, CS) + if (CS%id_boundary_forcing_h > 0) call post_data(CS%id_boundary_forcing_h, h, CS%diag, alt_h = h_diag) + endif + ! Boundary fluxes may have changed T, S, and h + call diag_update_remap_grids(CS%diag) + call cpu_clock_end(id_clock_remap) + if (CS%debug) then + call MOM_forcing_chksum("after applyBoundaryFluxes ", fluxes, G, haloshift=0) + call MOM_thermovar_chksum("after applyBoundaryFluxes ", tv, G) + call MOM_state_chksum("after applyBoundaryFluxes ", u, v, h, G, GV, haloshift=0) + endif + 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. @@ -902,6 +824,10 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! 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 @@ -937,204 +863,54 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (CS%debugConservation) call MOM_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, G) - ! Here, T and S are updated according to ea and eb. - ! If using the bulk mixed layer, T and S are also updated - ! by surface fluxes (in fluxes%*). - ! This is a very long block. - if (CS%bulkmixedlayer) then + ! calculate change in temperature & salinity due to dia-coordinate surface diffusion + if (associated(tv%T)) then - if (associated(tv%T)) then - call cpu_clock_begin(id_clock_tridiag) - ! Temperature and salinity (as state variables) are treated - ! differently from other tracers to insure massless layers that - ! are lighter than the mixed layer have temperatures and salinities - ! that correspond to their prescribed densities. - if (CS%massless_match_targets) then - !$OMP parallel do default (shared) private(h_tr,b1,d1,c1,b_denom_1) - do j=js,je - do i=is,ie - h_tr = hold(i,j,1) + h_neglect - b1(i) = 1.0 / (h_tr + eb(i,j,1)) - d1(i) = h_tr * b1(i) - tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1)) - tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1)) - enddo - do k=2,nkmb ; do i=is,ie - c1(i,k) = eb(i,j,k-1) * b1(i) - h_tr = hold(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)) - if (k kb(i,j)) then - c1(i,k) = eb(i,j,k-1) * b1(i) - h_tr = hold(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) - tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1)) - tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1)) - elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j)) - ! The bottommost buffer layer might entrain all the mass from some - ! of the interior layers that are thin and lighter in the coordinate - ! density than that buffer layer. The T and S of these newly - ! massless interior layers are unchanged. - tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%T(i,j,k) - tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%S(i,j,k) - endif - enddo ; enddo - - do k=nz-1,nkmb,-1 ; do i=is,ie - if (k >= kb(i,j)) then - tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1) - tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1) - endif - enddo ; enddo - do i=is,ie ; if (kb(i,j) <= nz) then - tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j)) - tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j)) - endif ; enddo - do k=nkmb-1,1,-1 ; do i=is,ie - tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1) - tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1) - enddo ; enddo - enddo ! end of j loop - else ! .not. massless_match_targets - ! This simpler form allows T & S to be too dense for the layers - ! between the buffer layers and the interior. - ! Changes: T, 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) + 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) + endif + call cpu_clock_begin(id_clock_tridiag) - endif ! endif for associated(T) - if (CS%debugConservation) call MOM_state_stats('BML tridiag', u, v, h, tv%T, tv%S, G) + ! 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 + ! from this limitation, in which case we can let salinity=0 and still + ! have salt conserved with SIS2 ice. So for SIS2, we can run with + ! BOUND_SALINITY=False in MOM.F90. + if (associated(tv%S) .and. associated(tv%salt_deficit)) & + call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - if ((CS%ML_mix_first > 0.0) .or. CS%use_geothermal) then - ! The mixed layer code has already been called, but there is some needed - ! bookkeeping. - !$OMP parallel do default(shared) + if (CS%diabatic_diff_tendency_diag) then do k=1,nz ; do j=js,je ; do i=is,ie - hold(i,j,k) = h_orig(i,j,k) - ea(i,j,k) = ea(i,j,k) + eaml(i,j,k) - eb(i,j,k) = eb(i,j,k) + ebml(i,j,k) + temp_diag(i,j,k) = tv%T(i,j,k) + saln_diag(i,j,k) = tv%S(i,j,k) enddo ; enddo ; enddo - if (CS%debug) then - call hchksum(ea, "after ea = ea + eaml",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(eb, "after eb = eb + ebml",G%HI,haloshift=0, scale=GV%H_to_m) - endif endif - if (CS%ML_mix_first < 1.0) then - ! Call the mixed layer code now, perhaps for a second time. - ! This subroutine (1) Cools the mixed layer. - ! (2) Performs convective adjustment by mixed layer entrainment. - ! (3) Heats the mixed layer and causes it to detrain to - ! Monin-Obukhov depth or minimum mixed layer depth. - ! (4) Uses any remaining TKE to drive mixed layer entrainment. - ! (5) Possibly splits the buffer layer into two isopycnal layers. - - call find_uv_at_h(u, v, hold, u_h, v_h, G, GV, ea, eb) - if (CS%debug) call MOM_state_chksum("find_uv_at_h1 ", u, v, h, G, GV, haloshift=0) - - dt_mix = min(dt,dt*(1.0 - CS%ML_mix_first)) - call cpu_clock_begin(id_clock_mixedlayer) - ! Changes: h, tv%T, tv%S, ea and eb (G is also inout???) - call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, & - G, GV, CS%bulkmixedlayer_CSp, CS%optics, & - Hml, CS%aggregate_FW_forcing, dt, last_call=.true.) - - if (CS%salt_reject_below_ML) & - call insert_brine(h, tv, G, GV, fluxes, nkmb, CS%diabatic_aux_CSp, dt_mix, & - CS%id_brine_lay) - - ! 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 - ! from this limitation, in which case we can let salinity=0 and still - ! have salt conserved with SIS2 ice. So for SIS2, we can run with - ! BOUND_SALINITY=False in MOM.F90. - if (associated(tv%S) .and. associated(tv%salt_deficit)) & - call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - - call cpu_clock_end(id_clock_mixedlayer) - if (showCallTree) call callTree_waypoint("done with 2nd bulkmixedlayer (diabatic)") - if (CS%debugConservation) call MOM_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, G) + ! 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 triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, tv%T, tv%S) endif - else ! following block for when NOT using BULKMIXEDLAYER - - - ! calculate change in temperature & salinity due to dia-coordinate surface diffusion - 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) - endif - 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 - ! from this limitation, in which case we can let salinity=0 and still - ! have salt conserved with SIS2 ice. So for SIS2, we can run with - ! BOUND_SALINITY=False in MOM.F90. - if (associated(tv%S) .and. associated(tv%salt_deficit)) & - call adjust_salt(h, tv, G, GV, CS%diabatic_aux_CSp) - - if (CS%diabatic_diff_tendency_diag) then - do k=1,nz ; do j=js,je ; do i=is,ie - temp_diag(i,j,k) = tv%T(i,j,k) - saln_diag(i,j,k) = tv%S(i,j,k) - enddo ; enddo ; enddo - 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 triDiagTS(G, GV, is, ie, js, je, hold, ea, eb, 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 - 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) - 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)) - if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G) + ! 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 + 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) + endif + call cpu_clock_end(id_clock_tridiag) + if (showCallTree) call callTree_waypoint("done with triDiagTS (diabatic)") - endif ! endif for the BULKMIXEDLAYER block + endif ! endif corresponding to if (associated(tv%T)) + if (CS%debugConservation) call MOM_state_stats('triDiagTS', u, v, h, tv%T, tv%S, G) if (CS%debug) then call MOM_state_chksum("after mixed layer ", u, v, h, G, GV, haloshift=0) @@ -1143,14 +919,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call hchksum(eb, "after mixed layer eb", G%HI, scale=GV%H_to_m) endif - if (.not. CS%useALEalgorithm) then - call cpu_clock_begin(id_clock_remap) - call regularize_layers(h, tv, dt, ea, eb, G, GV, CS%regularize_layers_CSp) - call cpu_clock_end(id_clock_remap) - if (showCallTree) call callTree_waypoint("done with regularize_layers (diabatic)") - if (CS%debugConservation) call MOM_state_stats('regularize_layers', u, v, h, tv%T, tv%S, G) - endif - ! Whenever thickness changes let the diag manager know, as the ! target grids for vertical remapping may need to be regenerated. call diag_update_remap_grids(CS%diag) @@ -1187,6 +955,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! mixing of passive tracers from massless boundary layers to interior call cpu_clock_begin(id_clock_tracers) + if (CS%mix_boundary_tracers) then Tr_ea_BBL = sqrt(dt*CS%Kd_BBL_tr) !$OMP parallel do default(shared) private(htot,in_boundary,add_ent) @@ -1235,17 +1004,12 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo - 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, 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) - else - call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug) - endif + call call_tracer_column_fns(h_prebound, h, ea, eb, 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) elseif (associated(visc%Kd_extra_S)) then ! extra diffusivity for passive tracers @@ -1265,56 +1029,31 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & eatr(i,j,k) = ea(i,j,k) + add_ent 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, 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) - else - call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug) - endif + call call_tracer_column_fns(h_prebound, h, eatr, ebtr, 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) 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, eatr, ebtr, 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) - else - call call_tracer_column_fns(hold, h, ea, eb, fluxes, Hml, dt, G, GV, tv, & - CS%optics, CS%tracer_flow_CSp, CS%debug) - endif + call call_tracer_column_fns(h_prebound, h, eatr, ebtr, 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) endif ! (CS%mix_boundary_tracers) - - call cpu_clock_end(id_clock_tracers) - ! sponges if (CS%use_sponge) then call cpu_clock_begin(id_clock_sponge) if (associated(CS%ALE_sponge_CSp)) then ! ALE sponge call apply_ALE_sponge(h, dt, G, CS%ALE_sponge_CSp, CS%Time) - else - ! Layer mode sponge - if (CS%bulkmixedlayer .and. associated(tv%eqn_of_state)) then - do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo - !$OMP parallel do default(shared) - do j=js,je - call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, Rcv_ml(:,j), & - is, ie-is+1, tv%eqn_of_state) - enddo - call apply_sponge(h, dt, G, GV, ea, eb, CS%sponge_CSp, Rcv_ml) - else - call apply_sponge(h, dt, G, GV, ea, eb, CS%sponge_CSp) - endif endif + call cpu_clock_end(id_clock_sponge) if (CS%debug) then call MOM_state_chksum("apply_sponge ", u, v, h, G, GV, haloshift=0) @@ -1337,29 +1076,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & enddo endif -! For momentum, it is only the net flux that homogenizes within -! the mixed layer. Vertical viscosity that is proportional to the -! mixed layer turbulence is applied elsewhere. - if (CS%bulkmixedlayer) then - if (CS%debug) then - call hchksum(ea, "before net flux rearrangement ea",G%HI, scale=GV%H_to_m) - call hchksum(eb, "before net flux rearrangement eb",G%HI, scale=GV%H_to_m) - endif - !$OMP parallel do default(shared) private(net_ent) - do j=js,je - do K=2,GV%nkml ; do i=is,ie - net_ent = ea(i,j,k) - eb(i,j,k-1) - ea(i,j,k) = max(net_ent, 0.0) - eb(i,j,k-1) = max(-net_ent, 0.0) - enddo ; enddo - enddo - if (CS%debug) then - call hchksum(ea, "after net flux rearrangement ea",G%HI, scale=GV%H_to_m) - call hchksum(eb, "after net flux rearrangement eb",G%HI, scale=GV%H_to_m) - endif - endif - -! Initialize halo regions of ea, eb, and hold to default values. + ! 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 @@ -1387,84 +1104,6 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & call cpu_clock_end(id_clock_pass) - if (.not. CS%useALEalgorithm) then - ! Use a tridiagonal solver to determine effect of the diapycnal - ! advection on velocity field. It is assumed that water leaves - ! or enters the ocean with the surface velocity. - if (CS%debug) then - call MOM_state_chksum("before u/v tridiag ", u, v, h, G, GV, haloshift=0) - call hchksum(ea, "before u/v tridiag ea",G%HI, scale=GV%H_to_m) - call hchksum(eb, "before u/v tridiag eb",G%HI, scale=GV%H_to_m) - call hchksum(hold, "before u/v tridiag hold",G%HI, scale=GV%H_to_m) - endif - call cpu_clock_begin(id_clock_tridiag) - !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval) - do j=js,je - do I=Isq,Ieq - if (associated(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,1) = u(I,j,1) - hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect - b1(I) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1))) - d1(I) = hval * b1(I) - u(I,j,1) = b1(I) * (hval * u(I,j,1)) - enddo - do k=2,nz ; do I=Isq,Ieq - if (associated(ADp%du_dt_dia)) ADp%du_dt_dia(I,j,k) = u(I,j,k) - c1(I,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(I) - eaval = ea(i,j,k) + ea(i+1,j,k) - hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect - b1(I) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(I)*eaval)) - d1(I) = (hval + d1(I)*eaval) * b1(I) - u(I,j,k) = (hval*u(I,j,k) + eaval*u(I,j,k-1))*b1(I) - enddo ; enddo - do k=nz-1,1,-1 ; do I=Isq,Ieq - u(I,j,k) = u(I,j,k) + c1(I,k+1)*u(I,j,k+1) - if (associated(ADp%du_dt_dia)) & - ADp%du_dt_dia(I,j,k) = (u(I,j,k) - ADp%du_dt_dia(I,j,k)) * Idt - enddo ; enddo - if (associated(ADp%du_dt_dia)) then - do I=Isq,Ieq - ADp%du_dt_dia(I,j,nz) = (u(I,j,nz)-ADp%du_dt_dia(I,j,nz)) * Idt - enddo - endif - enddo - if (CS%debug) then - call MOM_state_chksum("aft 1st loop tridiag ", u, v, h, G, GV, haloshift=0) - endif - !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval) - do J=Jsq,Jeq - do i=is,ie - if (associated(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,1) = v(i,J,1) - hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect - b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1))) - d1(I) = hval * b1(I) - v(i,J,1) = b1(i) * (hval * v(i,J,1)) - enddo - do k=2,nz ; do i=is,ie - if (associated(ADp%dv_dt_dia)) ADp%dv_dt_dia(i,J,k) = v(i,J,k) - c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i) - eaval = ea(i,j,k) + ea(i,j+1,k) - hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect - b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval)) - d1(i) = (hval + d1(i)*eaval) * b1(i) - v(i,J,k) = (hval*v(i,J,k) + eaval*v(i,J,k-1))*b1(i) - enddo ; enddo - do k=nz-1,1,-1 ; do i=is,ie - v(i,J,k) = v(i,J,k) + c1(i,k+1)*v(i,J,k+1) - if (associated(ADp%dv_dt_dia)) & - ADp%dv_dt_dia(i,J,k) = (v(i,J,k) - ADp%dv_dt_dia(i,J,k)) * Idt - enddo ; enddo - if (associated(ADp%dv_dt_dia)) then - do i=is,ie - ADp%dv_dt_dia(i,J,nz) = (v(i,J,nz)-ADp%dv_dt_dia(i,J,nz)) * Idt - enddo - endif - enddo - call cpu_clock_end(id_clock_tridiag) - if (CS%debug) then - call MOM_state_chksum("after u/v tridiag ", u, v, h, G, GV, haloshift=0) - endif - endif ! useALEalgorithm - call disable_averaging(CS%diag) ! Frazil formation keeps temperature above the freezing point. ! make_frazil is deliberately called at both the beginning and at From 42c4f937538db40aedf600239cbf3bf441a1c92a Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 29 May 2018 09:30:37 -0600 Subject: [PATCH 33/39] Add some comments --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index a26e44fe48..f283f6243a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -893,6 +893,7 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & 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 From 56c461cc8ad7eb31bf695f5f3f4d4ee67c93b7ec Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 29 May 2018 15:51:45 -0600 Subject: [PATCH 34/39] Adds back old double-diffusion code --- .../vertical/MOM_diabatic_driver.F90 | 21 +- .../vertical/MOM_set_diffusivity.F90 | 202 +++++++++++++++++- 2 files changed, 213 insertions(+), 10 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 698243a7f6..d6ca59183e 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -731,11 +731,16 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & ! 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 ! KPP boundary layer to interior, so this diffusivity can be computed when convenient. - if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T)) then - call cpu_clock_begin(id_clock_CVMix_ddiff) + if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T) .and. .not. & + CS%use_CVMix_ddiff) then + ! GMM, fix this + !call cpu_clock_begin(id_clock_CVMix_ddiff) + call cpu_clock_begin(id_clock_differential_diff) call differential_diffuse_T_S(h, tv, visc, dt, G, GV) - call cpu_clock_end(id_clock_CVMix_ddiff) + !call cpu_clock_end(id_clock_CVMix_ddiff) + call cpu_clock_end(id_clock_differential_diff) + if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)") if (CS%debugConservation) call MOM_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, G) @@ -1889,7 +1894,7 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, real :: Kd integer :: num_mode - logical :: use_temperature + logical :: use_temperature, differentialDiffusion type(vardesc) :: vd ! This "include" declares and sets the variable "version". @@ -1940,6 +1945,10 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, "If true, the diffusivity from ePBL is added to all\n"//& "other diffusivities. Otherwise, the larger of kappa-\n"//& "shear and ePBL diffusivities are used.", default=.true.) + call get_param(param_file, mod, "DOUBLE_DIFFUSION", differentialDiffusion, & + "If true, apply parameterization of double-diffusion.", & + default=.false. ) + CS%use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) CS%use_kappa_shear = kappa_shear_is_used(param_file) CS%use_CVMix_shear = CVMix_shear_is_used(param_file) @@ -2402,8 +2411,10 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=CLOCK_MODULE) id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=CLOCK_ROUTINE) id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=CLOCK_ROUTINE) + id_clock_differential_diff = -1 ; if (differentialDiffusion) & + id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=CLOCK_ROUTINE) id_clock_CVMix_ddiff = -1 ; if (CS%use_CVMix_ddiff) & - id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion)', grain=CLOCK_ROUTINE) + id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=CLOCK_ROUTINE) ! initialize the auxiliary diabatic driver module call diabatic_aux_init(Time, G, GV, param_file, diag, CS%diabatic_aux_CSp, & diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 903868795a..b97583aa1c 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -130,9 +130,13 @@ module MOM_set_diffusivity !! shear-driven diapycnal diffusivity. logical :: use_CVMix_shear !< If true, use one of the CVMix modules to find !! shear-driven diapycnal diffusivity. + logical :: double_diffusion !< If true, enable double-diffusive mixing using an old method. logical :: use_CVMix_ddiff !< If true, enable double-diffusive mixing via CVMix. logical :: simple_TKE_to_Kd !< If true, uses a simple estimate of Kd/TKE that !! does not rely on a layer-formulation. + real :: Max_Rrho_salt_fingers !< max density ratio for salt fingering + real :: Max_salt_diff_salt_fingers !< max salt diffusivity for salt fingers (m2/s) + real :: Kv_molecular !< molecular visc for double diff convect (m2/s) character(len=200) :: inputdir type(user_change_diff_CS), pointer :: user_change_diff_CSp => NULL() @@ -156,6 +160,10 @@ module MOM_set_diffusivity 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 end type set_diffusivity_CS @@ -166,9 +174,12 @@ module MOM_set_diffusivity Kd_BBL => NULL(),& ! BBL diffusivity at interfaces (m2/s) Kd_work => NULL(),& ! layer integrated work by diapycnal mixing (W/m2) maxTKE => NULL(),& ! energy required to entrain to h_max (m3/s3) - TKE_to_Kd => NULL() ! conversion rate (~1.0 / (G_Earth + dRho_lay)) + TKE_to_Kd => NULL(),& ! conversion rate (~1.0 / (G_Earth + dRho_lay)) ! between TKE dissipated within a layer and Kd ! in that layer, in m2 s-1 / m3 s-3 = s2 m-1 + KT_extra => NULL(),& ! double diffusion diffusivity for temp (m2/s) + KS_extra => NULL() ! double diffusion diffusivity for saln (m2/s) + end type diffusivity_diags ! Clocks @@ -236,7 +247,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & real, dimension(SZI_(G),SZK_(G)+1) :: & N2_int, & !< squared buoyancy frequency associated at interfaces (1/s2) - dRho_int !< locally ref potential density difference across interfaces (in s-2) smg: or kg/m3? + dRho_int, & !< locally ref potential density difference across interfaces (in s-2) smg: or kg/m3? + KT_extra, & !< double difusion diffusivity on temperature (m2/sec) + KS_extra ! double difusion diffusivity on salinity (m2/sec) real :: I_Rho0 ! inverse of Boussinesq density (m3/kg) real :: dissip ! local variable for dissipation calculations (W/m3) @@ -271,10 +284,10 @@ 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) .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: visc%Kd_extra_T and "//& - "visc%Kd_extra_S must be associated when USE_CVMIX_DDIFF is true.") + "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. ! If nothing else is specified, this will be the value used. @@ -299,6 +312,12 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (CS%id_TKE_to_Kd > 0) then allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz)) ; dd%TKE_to_Kd(:,:,:) = 0.0 endif + if ((CS%id_KT_extra > 0) .or. (CS%id_KT_extra_z > 0)) then + allocate(dd%KT_extra(isd:ied,jsd:jed,nz+1)) ; dd%KT_extra(:,:,:) = 0.0 + endif + if ((CS%id_KS_extra > 0) .or. (CS%id_KS_extra_z > 0)) then + allocate(dd%KS_extra(isd:ied,jsd:jed,nz+1)) ; dd%KS_extra(:,:,:) = 0.0 + endif if ((CS%id_Kd_BBL > 0) .or. (CS%id_Kd_BBL_z > 0)) then allocate(dd%Kd_BBL(isd:ied,jsd:jed,nz+1)) ; dd%Kd_BBL(:,:,:) = 0.0 endif @@ -375,7 +394,35 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! Add background mixing call calculate_bkgnd_mixing(h, tv, N2_lay, Kd, visc%Kv_slow, j, G, GV, CS%bkgnd_mixing_csp) - ! Apply double diffusion + ! Double-diffusion (old method) + if (CS%double_diffusion) then + call double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, KT_extra, KS_extra) + do K=2,nz ; do i=is,ie + if (KS_extra(i,K) > KT_extra(i,K)) then ! salt fingering + Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KT_extra(i,K) + Kd(i,j,k) = Kd(i,j,k) + 0.5*KT_extra(i,K) + visc%Kd_extra_S(i,j,k) = KS_extra(i,K)-KT_extra(i,K) + visc%Kd_extra_T(i,j,k) = 0.0 + elseif (KT_extra(i,K) > 0.0) then ! double-diffusive convection + Kd(i,j,k-1) = Kd(i,j,k-1) + 0.5*KS_extra(i,K) + Kd(i,j,k) = Kd(i,j,k) + 0.5*KS_extra(i,K) + visc%Kd_extra_T(i,j,k) = KT_extra(i,K)-KS_extra(i,K) + visc%Kd_extra_S(i,j,k) = 0.0 + else ! There is no double diffusion at this interface. + visc%Kd_extra_T(i,j,k) = 0.0 + visc%Kd_extra_S(i,j,k) = 0.0 + endif + enddo ; enddo + if (associated(dd%KT_extra)) then ; do K=1,nz+1 ; do i=is,ie + dd%KT_extra(i,j,K) = KT_extra(i,K) + enddo ; enddo ; endif + + if (associated(dd%KS_extra)) then ; do K=1,nz+1 ; do i=is,ie + dd%KS_extra(i,j,K) = KS_extra(i,K) + enddo ; enddo ; endif + endif + + ! Apply double diffusion via CVMix ! GMM, we need to pass HBL to compute_ddiff_coeffs, but it is not yet available. if (CS%use_CVMix_ddiff) then call compute_ddiff_coeffs(h, tv, G, GV, j, visc%Kd_extra_T, visc%Kd_extra_S, CS%CVMix_ddiff_csp) @@ -562,11 +609,26 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & endif + if (CS%id_KT_extra > 0) call post_data(CS%id_KT_extra, dd%KT_extra, CS%diag) + if (CS%id_KS_extra > 0) call post_data(CS%id_KS_extra, dd%KS_extra, CS%diag) if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, dd%Kd_BBL, CS%diag) + if (CS%id_KT_extra_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_KT_extra_z + z_ptrs(num_z_diags)%p => dd%KT_extra + endif + + if (CS%id_KS_extra_z > 0) then + num_z_diags = num_z_diags + 1 + z_ids(num_z_diags) = CS%id_KS_extra_z + z_ptrs(num_z_diags)%p => dd%KS_extra + endif + if (CS%id_Kd_BBL_z > 0) then num_z_diags = num_z_diags + 1 z_ids(num_z_diags) = CS%id_Kd_BBL_z + z_ptrs(num_z_diags)%p => dd%KS_extra endif if (num_z_diags > 0) & @@ -577,6 +639,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & if (associated(dd%Kd_user)) deallocate(dd%Kd_user) if (associated(dd%maxTKE)) deallocate(dd%maxTKE) if (associated(dd%TKE_to_Kd)) deallocate(dd%TKE_to_Kd) + if (associated(dd%KT_extra)) deallocate(dd%KT_extra) + if (associated(dd%KS_extra)) deallocate(dd%KS_extra) if (associated(dd%Kd_BBL)) deallocate(dd%Kd_BBL) if (showCallTree) call callTree_leave("set_diffusivity()") @@ -929,6 +993,97 @@ subroutine find_N2(h, tv, T_f, S_f, fluxes, j, G, GV, CS, dRho_int, & end subroutine find_N2 +!> This subroutine sets the additional diffusivities of temperature and +!! salinity due to double diffusion, using the same functional form as is +!! used in MOM4.1, and taken from an NCAR technical note (REF?) that updates +!! what was in Large et al. (1994). All the coefficients here should probably +!! be made run-time variables rather than hard-coded constants. +!! +!! \todo Find reference for NCAR tech note above. +subroutine double_diffusion(tv, h, T_f, S_f, j, G, GV, CS, Kd_T_dd, Kd_S_dd) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any available + !! thermodynamic fields; absent fields have NULL + !! ptrs. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: h !< Layer thicknesses, in H (usually m or kg m-2). + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: T_f !< layer temp in C with the values in massless layers + !! filled vertically by diffusion. + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & + intent(in) :: S_f !< Layer salinities in PPT with values in massless + !! layers filled vertically by diffusion. + integer, intent(in) :: j !< Meridional index upon which to work. + type(set_diffusivity_CS), pointer :: CS !< Module control structure. + real, dimension(SZI_(G),SZK_(G)+1), & + intent(out) :: Kd_T_dd !< Interface double diffusion diapycnal + !! diffusivity for temp (m2/sec). + real, dimension(SZI_(G),SZK_(G)+1), & + intent(out) :: Kd_S_dd !< Interface double diffusion diapycnal + !! diffusivity for saln (m2/sec). + + real, dimension(SZI_(G)) :: & + dRho_dT, & ! partial derivatives of density wrt temp (kg m-3 degC-1) + dRho_dS, & ! partial derivatives of density wrt saln (kg m-3 PPT-1) + pres, & ! pressure at each interface (Pa) + Temp_int, & ! temp and saln at interfaces + Salin_int + + real :: alpha_dT ! density difference between layers due to temp diffs (kg/m3) + real :: beta_dS ! density difference between layers due to saln diffs (kg/m3) + + real :: Rrho ! vertical density ratio + real :: diff_dd ! factor for double-diffusion + real :: prandtl ! flux ratio for diffusive convection regime + + real, parameter :: Rrho0 = 1.9 ! limit for double-diffusive density ratio + real, parameter :: dsfmax = 1.e-4 ! max diffusivity in case of salt fingering + real, parameter :: Kv_molecular = 1.5e-6 ! molecular viscosity (m2/sec) + + integer :: i, k, is, ie, nz + is = G%isc ; ie = G%iec ; nz = G%ke + + if (associated(tv%eqn_of_state)) then + do i=is,ie + pres(i) = 0.0 ; Kd_T_dd(i,1) = 0.0 ; Kd_S_dd(i,1) = 0.0 + Kd_T_dd(i,nz+1) = 0.0 ; Kd_S_dd(i,nz+1) = 0.0 + enddo + do K=2,nz + do i=is,ie + pres(i) = pres(i) + GV%H_to_Pa*h(i,j,k-1) + Temp_Int(i) = 0.5 * (T_f(i,j,k-1) + T_f(i,j,k)) + Salin_Int(i) = 0.5 * (S_f(i,j,k-1) + S_f(i,j,k)) + enddo + call calculate_density_derivs(Temp_int, Salin_int, pres, & + dRho_dT(:), dRho_dS(:), is, ie-is+1, tv%eqn_of_state) + + do i=is,ie + alpha_dT = -1.0*dRho_dT(i) * (T_f(i,j,k-1) - T_f(i,j,k)) + beta_dS = dRho_dS(i) * (S_f(i,j,k-1) - S_f(i,j,k)) + + if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then ! salt finger case + Rrho = min(alpha_dT/beta_dS,Rrho0) + diff_dd = 1.0 - ((RRho-1.0)/(RRho0-1.0)) + diff_dd = dsfmax*diff_dd*diff_dd*diff_dd + Kd_T_dd(i,K) = 0.7*diff_dd + Kd_S_dd(i,K) = diff_dd + elseif ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then ! diffusive convection + Rrho = alpha_dT/beta_dS + diff_dd = Kv_molecular*0.909*exp(4.6*exp(-0.54*(1/Rrho-1))) + prandtl = 0.15*Rrho + if (Rrho > 0.5) prandtl = (1.85-0.85/Rrho)*Rrho + Kd_T_dd(i,K) = diff_dd + Kd_S_dd(i,K) = prandtl*diff_dd + else + Kd_T_dd(i,K) = 0.0 ; Kd_S_dd(i,K) = 0.0 + endif + enddo + enddo + endif + +end subroutine double_diffusion + !> This routine adds diffusion sustained by flow energy extracted by bottom drag. subroutine add_drag_diffusivity(h, u, v, tv, fluxes, visc, j, TKE_to_Kd, & maxTKE, kb, G, GV, CS, Kd, Kd_int, Kd_BBL) @@ -1945,6 +2100,43 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp endif endif + call get_param(param_file, mdl, "DOUBLE_DIFFUSION", CS%double_diffusion, & + "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.", & + default=2.55, units="nondim") + call get_param(param_file, mdl, "MAX_SALT_DIFF_SALT_FINGERS", CS%Max_salt_diff_salt_fingers, & + "Maximum salt diffusivity for salt fingering regime.", & + default=1.e-4, units="m2 s-1") + call get_param(param_file, mdl, "KV_MOLECULAR", CS%Kv_molecular, & + "Molecular viscosity for calculation of fluxes under \n"//& + "double-diffusive convection.", default=1.5e-6, units="m2 s-1") + ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults. + + CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for temperature', 'm2 s-1') + + CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for salinity', 'm2 s-1') + + if (associated(diag_to_Z_CSp)) then + vd = var_desc("KT_extra", "m2 s-1", & + "Double-Diffusive Temperature Diffusivity, interpolated to z", & + z_grid='z') + CS%id_KT_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + vd = var_desc("KS_extra", "m2 s-1", & + "Double-Diffusive Salinity Diffusivity, interpolated to z",& + z_grid='z') + CS%id_KS_extra_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + vd = var_desc("Kd_BBL", "m2 s-1", & + "Bottom Boundary Layer Diffusivity", z_grid='z') + CS%id_Kd_BBL_z = register_Zint_diag(vd, CS%diag_to_Z_CSp, Time) + endif + endif + if (CS%user_change_diff) then call user_change_diff_init(Time, G, param_file, diag, CS%user_change_diff_CSp) endif From 71e4beab97a43041587b2c668fa9662602cc6e1a Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 29 May 2018 16:06:10 -0600 Subject: [PATCH 35/39] Add fatal error if multiple double-diffusion options are enabled --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index d6ca59183e..f991404149 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -1950,6 +1950,13 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, default=.false. ) CS%use_CVMix_ddiff = CVMix_ddiff_is_used(param_file) + + if (CS%use_CVMix_ddiff .and. differentialDiffusion) then + call MOM_error(FATAL, 'diabatic_driver_init: '// & + 'Multiple double-diffusion options selected (DOUBLE_DIFFUSION and'//& + 'USE_CVMIX_DDIFF), please disable all but one option to proceed.') + endif + CS%use_kappa_shear = kappa_shear_is_used(param_file) CS%use_CVMix_shear = CVMix_shear_is_used(param_file) From 474d664a6f023a601f061bcac9cef0d2189dde5f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 30 May 2018 08:23:04 -0600 Subject: [PATCH 36/39] Set clock for double-diffusion via CVMix --- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 7 +------ src/parameterizations/vertical/MOM_set_diffusivity.F90 | 6 +++++- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index eb43f4b0c2..963547f3c0 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -251,7 +251,7 @@ module MOM_diabatic_driver integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap -integer :: id_clock_kpp, id_clock_CVMix_ddiff +integer :: id_clock_kpp contains @@ -687,11 +687,8 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, & if (associated(visc%Kd_extra_T) .and. associated(visc%Kd_extra_S) .and. associated(tv%T) .and. .not. & CS%use_CVMix_ddiff) then - ! GMM, fix this - !call cpu_clock_begin(id_clock_CVMix_ddiff) call cpu_clock_begin(id_clock_differential_diff) call differential_diffuse_T_S(h, tv, visc, dt, G, GV) - !call cpu_clock_end(id_clock_CVMix_ddiff) call cpu_clock_end(id_clock_differential_diff) if (showCallTree) call callTree_waypoint("done with differential_diffuse_T_S (diabatic)") @@ -2060,8 +2057,6 @@ subroutine diabatic_driver_init(Time, G, GV, param_file, useALEalgorithm, diag, id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=CLOCK_ROUTINE) id_clock_differential_diff = -1 ; if (differentialDiffusion) & id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=CLOCK_ROUTINE) - id_clock_CVMix_ddiff = -1 ; if (CS%use_CVMix_ddiff) & - id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=CLOCK_ROUTINE) ! initialize the auxiliary diabatic driver module call diabatic_aux_init(Time, G, GV, param_file, diag, CS%diabatic_aux_CSp, & diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index b97583aa1c..93018c9dac 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -183,7 +183,7 @@ module MOM_set_diffusivity end type diffusivity_diags ! Clocks -integer :: id_clock_kappaShear +integer :: id_clock_kappaShear, id_clock_CVMix_ddiff contains @@ -425,7 +425,9 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, & ! Apply double diffusion via CVMix ! GMM, we need to pass HBL to compute_ddiff_coeffs, but it is not yet available. if (CS%use_CVMix_ddiff) then + call cpu_clock_begin(id_clock_CVMix_ddiff) call compute_ddiff_coeffs(h, tv, G, GV, j, visc%Kd_extra_T, visc%Kd_extra_S, CS%CVMix_ddiff_csp) + call cpu_clock_end(id_clock_CVMix_ddiff) endif ! Add the input turbulent diffusivity. @@ -2154,6 +2156,8 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp ! CVMix double diffusion mixing CS%use_CVMix_ddiff = CVMix_ddiff_init(Time, G, GV, param_file, CS%diag, CS%CVMix_ddiff_csp) + if (CS%use_CVMix_ddiff) & + id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=CLOCK_MODULE) end subroutine set_diffusivity_init From c527d5a218099dfc72a4c42efcbf5b4650062106 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 30 May 2018 16:09:32 -0600 Subject: [PATCH 37/39] 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 309b4d41c91b1abfe27b62e3e5e8beabd4be1332 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Wed, 30 May 2018 16:49:55 -0600 Subject: [PATCH 38/39] Fix if statement for fatal error when using double diffusion (old and CVMix) --- .../vertical/MOM_set_diffusivity.F90 | 21 +++++++++++-------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 93018c9dac..2e9b7553ab 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -284,9 +284,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 +2106,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.", & @@ -2118,12 +2119,6 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp "double-diffusive convection.", default=1.5e-6, units="m2 s-1") ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults. - CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & - 'Double-diffusive diffusivity for temperature', 'm2 s-1') - - CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & - 'Double-diffusive diffusivity for salinity', 'm2 s-1') - if (associated(diag_to_Z_CSp)) then vd = var_desc("KT_extra", "m2 s-1", & "Double-Diffusive Temperature Diffusivity, interpolated to z", & @@ -2159,6 +2154,14 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp if (CS%use_CVMix_ddiff) & id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=CLOCK_MODULE) + if (CS%use_CVMix_ddiff .or. CS%double_diffusion) then + CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for temperature', 'm2 s-1') + + CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for salinity', 'm2 s-1') + endif + end subroutine set_diffusivity_init !> Clear pointers and dealocate memory From 9e365566cc10ae102cde4c293ff02fface1c80c3 Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Thu, 31 May 2018 08:52:58 -0600 Subject: [PATCH 39/39] Fix a bug This commit fixes a bug in the if statement that checks if visc%Kd_extra_T and visc%Kd_extra_S are associated when either use_CVMix_ddiff or double_diffusion are used. --- .../vertical/MOM_set_diffusivity.F90 | 20 +++++++++---------- 1 file changed, 9 insertions(+), 11 deletions(-) diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 2e9b7553ab..c40b3b0a2b 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -284,9 +284,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: both 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. @@ -2119,6 +2119,12 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp "double-diffusive convection.", default=1.5e-6, units="m2 s-1") ! The default molecular viscosity follows the CCSM4.0 and MOM4p1 defaults. + CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for temperature', 'm2 s-1') + + CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & + 'Double-diffusive diffusivity for salinity', 'm2 s-1') + if (associated(diag_to_Z_CSp)) then vd = var_desc("KT_extra", "m2 s-1", & "Double-Diffusive Temperature Diffusivity, interpolated to z", & @@ -2154,14 +2160,6 @@ subroutine set_diffusivity_init(Time, G, GV, param_file, diag, CS, diag_to_Z_CSp if (CS%use_CVMix_ddiff) & id_clock_CVMix_ddiff = cpu_clock_id('(Double diffusion via CVMix)', grain=CLOCK_MODULE) - if (CS%use_CVMix_ddiff .or. CS%double_diffusion) then - CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, & - 'Double-diffusive diffusivity for temperature', 'm2 s-1') - - CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, & - 'Double-diffusive diffusivity for salinity', 'm2 s-1') - endif - end subroutine set_diffusivity_init !> Clear pointers and dealocate memory