From 8a4a5edd39d568e8286a1c1dbf92e72c79fa37a6 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Fri, 13 Apr 2018 11:44:24 -0600 Subject: [PATCH 01/18] distinguish profiles and cvmix schemes --- .../vertical/MOM_tidal_mixing.F90 | 106 +++++++++--------- 1 file changed, 56 insertions(+), 50 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index a59af01afb..c6d522fa93 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -34,6 +34,7 @@ module MOM_tidal_mixing !> Containers for tidal mixing diagnostics type, public :: tidal_mixing_diags + ! TODO: private real, pointer, dimension(:,:,:) :: & Kd_itidal => NULL(),& ! internal tide diffusivity at interfaces (m2/s) Fl_itidal => NULL(),& ! vertical flux of tidal turbulent dissipation (m3/s3) @@ -61,6 +62,7 @@ module MOM_tidal_mixing !> Control structure for tidal mixing module. type, public :: tidal_mixing_cs logical :: debug = .true. + ! TODO: private ! Parameters logical :: int_tide_dissipation ! Internal tide conversion (from barotropic) with @@ -128,9 +130,10 @@ module MOM_tidal_mixing real :: min_thickness ! Minimum thickness allowed [m] ! CVMix-specific parameters + integer :: cvmix_tidal_scheme = -1 ! 1 for Simmons, 2 for Schmittner type(cvmix_tidal_params_type) :: cvmix_tidal_params - type(cvmix_global_params_type) :: cvmix_glb_params ! for Prandtl number only - real :: tidal_max_coef ! maximum allowable tidal diffusivity. [m^2/s] + type(cvmix_global_params_type) :: cvmix_glb_params ! for Prandtl number only + real :: tidal_max_coef ! maximum allowable tidal diffusivity. [m^2/s] ! Data containers real, pointer, dimension(:,:) :: TKE_Niku => NULL() @@ -139,7 +142,7 @@ module MOM_tidal_mixing real, pointer, dimension(:,:) :: mask_itidal => NULL() real, pointer, dimension(:,:) :: h2 => NULL() real, pointer, dimension(:,:) :: tideamp => NULL() ! RMS tidal amplitude (m/s) - real, allocatable,dimension(:,:) :: tidal_qe_2d ! q*E(x,y) + real, allocatable,dimension(:,:) :: tidal_qe_2d ! q*E(x,y) ! TODO: make this E(x,y) only ! Diagnostics type(diag_ctrl), pointer :: diag => NULL() ! structure to regulate diagn output timing @@ -174,12 +177,12 @@ module MOM_tidal_mixing character(len=40) :: mdl = "MOM_tidal_mixing" !< This module's name. character*(20), parameter :: STLAURENT_PROFILE_STRING = "STLAURENT_02" character*(20), parameter :: POLZIN_PROFILE_STRING = "POLZIN_09" -character*(20), parameter :: SIMMONS_PROFILE_STRING = "SIMMONS" -character*(20), parameter :: SCHMITTNER_PROFILE_STRING = "SCHMITTNER" integer, parameter :: STLAURENT_02 = 1 integer, parameter :: POLZIN_09 = 2 -integer, parameter :: SIMMONS_04 = 3 -integer, parameter :: SCHMITTNER = 4 +character*(20), parameter :: SIMMONS_SCHEME_STRING = "SIMMONS" +character*(20), parameter :: SCHMITTNER_SCHEME_STRING = "SCHMITTNER" +integer, parameter :: SIMMONS_04 = 1 +integer, parameter :: SCHMITTNER = 2 contains @@ -197,7 +200,7 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, ! Local variables logical :: read_tideamp character(len=20) :: tmpstr, int_tide_profile_str - character(len=20) :: default_profile_string, tidal_energy_type + character(len=20) :: cvmix_tidal_scheme_str, tidal_energy_type character(len=200) :: filename, h2_file, Niku_TKE_input_file character(len=200) :: tidal_energy_file, tideamp_file type(vardesc) :: vd @@ -239,40 +242,47 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, "drive diapycnal mixing, along the lines of St. Laurent \n"//& "et al. (2002) and Simmons et al. (2004).", default=CS%use_cvmix_tidal) if (CS%int_tide_dissipation) then - default_profile_string = STLAURENT_PROFILE_STRING - if (CS%use_cvmix_tidal) default_profile_string = SIMMONS_PROFILE_STRING - call get_param(param_file, mdl, "INT_TIDE_PROFILE", int_tide_profile_str, & - "INT_TIDE_PROFILE selects the vertical profile of energy \n"//& - "dissipation with INT_TIDE_DISSIPATION. Valid values are:\n"//& - "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//& - "\t decay profile.\n"//& - "\t POLZIN_09 - Use the Polzin WKB-streched algebraic \n"//& - "\t decay profile.", & - default=default_profile_string) - ! TODO: list the newly available profile selections - int_tide_profile_str = uppercase(int_tide_profile_str) - select case (int_tide_profile_str) - case (STLAURENT_PROFILE_STRING) ; CS%int_tide_profile = STLAURENT_02 - case (POLZIN_PROFILE_STRING) ; CS%int_tide_profile = POLZIN_09 - case (SIMMONS_PROFILE_STRING) ; CS%int_tide_profile = SIMMONS_04 - case (SCHMITTNER_PROFILE_STRING) ; CS%int_tide_profile = SCHMITTNER - case default - call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & - "#define INT_TIDE_PROFILE "//trim(int_tide_profile_str)//" found in input file.") - end select - ! Check profile consistency - if (CS%use_cvmix_tidal .and. (CS%int_tide_profile.eq.STLAURENT_02 .or. & - CS%int_tide_profile.eq.POLZIN_09)) then - call MOM_error(FATAL, "tidal_mixing_init: Tidal mixing profile"// & - " "//trim(int_tide_profile_str)//" unavailable in CVMix. Available "//& - "profiles in CVMix are "//trim(SIMMONS_PROFILE_STRING)//" and "//& - trim(SCHMITTNER_PROFILE_STRING)//".") - else if (.not.CS%use_cvmix_tidal .and. (CS%int_tide_profile.eq.SIMMONS_04.or. & - CS%int_tide_profile.eq.SCHMITTNER)) then - call MOM_error(FATAL, "tidal_mixing_init: Tidal mixing profiles "// & - trim(SIMMONS_PROFILE_STRING)//" and "//trim(SCHMITTNER_PROFILE_STRING)//& - " are available only when USE_CVMIX_TIDAL is True.") + ! Read in CVMix tidal scheme if CVMix tidal mixing is on + if (CS%use_cvmix_tidal) then + call get_param(param_file, mdl, "CVMIX_TIDAL_SCHEME", cvmix_tidal_scheme_str, & + "CVMIX_TIDAL_SCHEME selects the CVMix tidal mixing\n"//& + "scheme with INT_TIDE_DISSIPATION. Valid values are:\n"//& + "\t SIMMONS - Use the Simmons et al (2004) tidal \n"//& + "\t mixing scheme.\n"//& + "\t SCHMITTNER - Use the Schmittner et al (2014) tidal \n"//& + "\t mixing scheme.", & + default=SIMMONS_SCHEME_STRING) + cvmix_tidal_scheme_str = uppercase(cvmix_tidal_scheme_str) + + select case (cvmix_tidal_scheme_str) + case (SIMMONS_SCHEME_STRING) ; CS%cvmix_tidal_scheme = SIMMONS_04 + case (SCHMITTNER_SCHEME_STRING) ; CS%cvmix_tidal_scheme = SCHMITTNER + case default + call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & + "#define CVMIX_TIDAL_SCHEME "//trim(cvmix_tidal_scheme_str)//" found in input file.") + end select + endif ! CS%use_cvmix_tidal + + ! Read in vertical profile of tidal energy dissipation + if ( CS%cvmix_tidal_scheme.eq.SCHMITTNER .or. .not. CS%use_cvmix_tidal) then + call get_param(param_file, mdl, "INT_TIDE_PROFILE", int_tide_profile_str, & + "INT_TIDE_PROFILE selects the vertical profile of energy \n"//& + "dissipation with INT_TIDE_DISSIPATION. Valid values are:\n"//& + "\t STLAURENT_02 - Use the St. Laurent et al exponential \n"//& + "\t decay profile.\n"//& + "\t POLZIN_09 - Use the Polzin WKB-streched algebraic \n"//& + "\t decay profile.", & + default=STLAURENT_PROFILE_STRING) + int_tide_profile_str = uppercase(int_tide_profile_str) + + select case (int_tide_profile_str) + case (STLAURENT_PROFILE_STRING) ; CS%int_tide_profile = STLAURENT_02 + case (POLZIN_PROFILE_STRING) ; CS%int_tide_profile = POLZIN_09 + case default + call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & + "#define INT_TIDE_PROFILE "//trim(int_tide_profile_str)//" found in input file.") + end select endif else if (CS%use_cvmix_tidal) then @@ -317,10 +327,6 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, if ((CS%Int_tide_dissipation .and. (CS%int_tide_profile == POLZIN_09)) .or. & (CS%lee_wave_dissipation .and. (CS%lee_wave_profile == POLZIN_09))) then - if (CS%use_cvmix_tidal) then - call MOM_error(FATAL, "tidal_mixing_init: Polzin scheme cannot "// & - "be used when CVMix tidal mixing scheme is active.") - end if call get_param(param_file, mdl, "NU_POLZIN", CS%Nu_Polzin, & "When the Polzin decay profile is used, this is a \n"//& "non-dimensional constant in the expression for the \n"//& @@ -489,14 +495,14 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, do_not_log=.true.) call cvmix_put(CS%cvmix_glb_params,'Prandtl',prandtl_tidal) - int_tide_profile_str = lowercase(int_tide_profile_str) + cvmix_tidal_scheme_str = lowercase(cvmix_tidal_scheme_str) ! TODO: check parameter consistency. (see POP::tidal_mixing.F90::tidal_check) ! Set up CVMix call cvmix_init_tidal(CVmix_tidal_params_user = CS%cvmix_tidal_params, & - mix_scheme = int_tide_profile_str, & + mix_scheme = cvmix_tidal_scheme_str, & efficiency = CS%Mu_itides, & vertical_decay_scale = CS%int_tide_decay_scale, & max_coefficient = CS%tidal_max_coef, & @@ -646,7 +652,7 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) is = G%isc ; ie = G%iec dd => CS%dd - select case (CS%int_tide_profile) + select case (CS%cvmix_tidal_scheme) case (SIMMONS_04) do i=is,ie @@ -713,8 +719,8 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) ! TODO: case (SCHMITTNER) case default - call MOM_error(FATAL, "tidal_mixing_init: The selected"// & - " INT_TIDE_PROFILE is unavailable in CVMix") + call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & + "#define CVMIX_TIDAL_SCHEME found in input file.") end select end subroutine calculate_cvmix_tidal From 20cc62b8ce8787870714a1658efab0df3923bd65 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 17 Apr 2018 13:10:15 -0600 Subject: [PATCH 02/18] add capability to read ER03 energy file --- .../vertical/MOM_tidal_mixing.F90 | 101 ++++++++++++++++-- 1 file changed, 92 insertions(+), 9 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index c6d522fa93..e8b383365d 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -132,8 +132,9 @@ module MOM_tidal_mixing ! CVMix-specific parameters integer :: cvmix_tidal_scheme = -1 ! 1 for Simmons, 2 for Schmittner type(cvmix_tidal_params_type) :: cvmix_tidal_params - type(cvmix_global_params_type) :: cvmix_glb_params ! for Prandtl number only - real :: tidal_max_coef ! maximum allowable tidal diffusivity. [m^2/s] + type(cvmix_global_params_type) :: cvmix_glb_params ! for Prandtl number only + real :: tidal_max_coef ! maximum allowable tidal diffusivity. [m^2/s] + real :: tidal_diss_lim_tc ! dissipation limit for tidal-energy-constituent data ! Data containers real, pointer, dimension(:,:) :: TKE_Niku => NULL() @@ -477,6 +478,10 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, call get_param(param_file, mdl, "TIDAL_MAX_COEF", CS%tidal_max_coef, & "largest acceptable value for tidal diffusivity", & units="m^2/s", default=50e-4) ! the default is 50e-4 in CVMIX, 100e-4 in POP. + call get_param(param_file, mdl, "TIDAL_DISS_LIM_TC", CS%tidal_diss_lim_tc, & + "Min allowable depth for dissipation for tidal-energy-constituent data. \n"//& + "No dissipation contribution is applied above TIDAL_DISS_LIM_TC.", & + units="m", default=0.0) call get_param(param_file, mdl, "TIDAL_ENERGY_FILE",tidal_energy_file, & "The path to the file containing tidal energy \n"//& "dissipation. Used with CVMix tidal mixing schemes.", & @@ -1133,7 +1138,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. & @@ -1293,23 +1298,101 @@ subroutine read_tidal_energy(G, tidal_energy_type, tidal_energy_file, CS) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed if (.not. allocated(CS%tidal_qe_2d)) allocate(CS%tidal_qe_2d(isd:ied,jsd:jed)) - allocate(tidal_energy_flux_2d(isd:ied,jsd:jed)) select case (uppercase(tidal_energy_type(1:4))) - case ('JAYN') ! Jayne 2009 input tidal energy flux + case ('JAYN') ! Jayne 2009 + allocate(tidal_energy_flux_2d(isd:ied,jsd:jed)) call MOM_read_data(tidal_energy_file,'wave_dissipation',tidal_energy_flux_2d, G%domain) CS%tidal_qe_2d = (CS%Gamma_itides) * tidal_energy_flux_2d + deallocate(tidal_energy_flux_2d) + case ('ER03') ! Egbert & Ray 2003 + call read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) case default call MOM_error(FATAL, "read_tidal_energy: Unknown tidal energy file type.") - ! TODO: add more tidal energy file types, e.g., Arbic, ER03, GN13, LGM0, etc. - ! see POP::tidal_mixing.F90 end select - deallocate(tidal_energy_flux_2d) - end subroutine read_tidal_energy +subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + character(len=20), intent(in) :: tidal_energy_type + character(len=200), intent(in) :: tidal_energy_file + type(tidal_mixing_cs), pointer :: CS + + ! local + integer :: k, isd, ied, jsd, jed, nz + real, parameter :: p33 = 1.0/3.0 + real, dimension(SZK_(G)) :: & + z_t, & ! depth from surface to midpoint of input layer + z_w ! depth from surface to top of input layer + real, dimension(SZI_(G),SZJ_(G)) :: & + tidal_qk1, & ! qk1 coefficient used in Schmittner & Egbert + tidal_qo1 ! qo1 coefficient used in Schmittner & Egbert + real, allocatable, dimension(:,:,:) :: & + tc_m2, & ! input lunar semidiurnal tidal energy flux [W/m^2] + tc_s2, & ! input solar semidiurnal tidal energy flux [W/m^2] + tc_k1, & ! input lunar diurnal tidal energy flux [W/m^2] + tc_o1, & ! input lunar diurnal tidal energy flux [W/m^2] + tidal_qe_3d ! sum_tc(q_tc*TC(x,y,z)) = q*E(x,y,z) + + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = G%ke + + allocate(tc_m2(isd:ied,jsd:jed,nz), & + tc_s2(isd:ied,jsd:jed,nz), & + tc_k1(isd:ied,jsd:jed,nz), & + tc_o1(isd:ied,jsd:jed,nz), & + tidal_qe_3d(isd:ied,jsd:jed,nz) ) + + ! read in tidal constituents + ! (NOTE: input z coordinates may differ from the model coordinates, which is fine.) + call MOM_read_data(tidal_energy_file, 'M2', tc_m2, G%domain) + call MOM_read_data(tidal_energy_file, 'S2', tc_s2, G%domain) + call MOM_read_data(tidal_energy_file, 'K1', tc_k1, G%domain) + call MOM_read_data(tidal_energy_file, 'O1', tc_o1, G%domain) + call MOM_read_data(tidal_energy_file, 'z_t', z_t) + call MOM_read_data(tidal_energy_file, 'z_w', z_w) + + ! form tidal_qe_3d from weighted tidal constituents + tidal_qe_3d = 0.0 + + where (abs(G%geoLatT(:,:)) < 30.0) + tidal_qk1(:,:) = p33 + tidal_qo1(:,:) = p33 + elsewhere + tidal_qk1(:,:) = 1.0 + tidal_qo1(:,:) = 1.0 + endwhere + + do k=1,nz + where (z_t(k) <= G%bathyT(:,:) .and. z_w(k) > CS%tidal_diss_lim_tc) + tidal_qe_3d(:,:,k) = p33*tc_m2(:,:,k) + p33*tc_s2(:,:,k) + & + tidal_qk1*tc_k1(:,:,k) + tidal_qo1*tc_o1(:,:,k) + endwhere + enddo + + ! test if qE is positive + if (any(tidal_qe_3d<0)) then + call MOM_error(FATAL, "read_tidal_constituents: Negative tidal_qe_3d terms.") + endif + + ! collapse 3D q*E to 2D q*E + CS%tidal_qe_2d = 0.0 + do k=1,nz + where (z_t(k) <= G%bathyT(:,:)) + CS%tidal_qe_2d(:,:) = CS%tidal_qe_2d(:,:) + tidal_qe_3d(:,:,k) + endwhere + enddo + + deallocate(tc_m2) + deallocate(tc_s2) + deallocate(tc_k1) + deallocate(tc_o1) + deallocate(tidal_qe_3d) + +end subroutine read_tidal_constituents + + !> Clear pointers and deallocate memory subroutine tidal_mixing_end(CS) type(tidal_mixing_cs), pointer :: CS ! This module's control structure From f13294011be66669d5413096057ca0ea67f429be Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 17 Apr 2018 14:52:02 -0600 Subject: [PATCH 03/18] add call to compute_Schmittner_invariant --- .../vertical/MOM_tidal_mixing.F90 | 54 +++++++++++++++++-- 1 file changed, 49 insertions(+), 5 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index e8b383365d..2d8415d162 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -19,6 +19,7 @@ module MOM_tidal_mixing use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc use cvmix_tidal, only : cvmix_init_tidal, cvmix_compute_Simmons_invariant use cvmix_tidal, only : cvmix_coeffs_tidal, cvmix_tidal_params_type +use cvmix_tidal, only : cvmix_compute_Schmittner_invariant use cvmix_kinds_and_types, only : cvmix_global_params_type use cvmix_put_get, only : cvmix_put @@ -182,7 +183,7 @@ module MOM_tidal_mixing integer, parameter :: POLZIN_09 = 2 character*(20), parameter :: SIMMONS_SCHEME_STRING = "SIMMONS" character*(20), parameter :: SCHMITTNER_SCHEME_STRING = "SCHMITTNER" -integer, parameter :: SIMMONS_04 = 1 +integer, parameter :: SIMMONS = 1 integer, parameter :: SCHMITTNER = 2 contains @@ -257,7 +258,7 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, cvmix_tidal_scheme_str = uppercase(cvmix_tidal_scheme_str) select case (cvmix_tidal_scheme_str) - case (SIMMONS_SCHEME_STRING) ; CS%cvmix_tidal_scheme = SIMMONS_04 + case (SIMMONS_SCHEME_STRING) ; CS%cvmix_tidal_scheme = SIMMONS case (SCHMITTNER_SCHEME_STRING) ; CS%cvmix_tidal_scheme = SCHMITTNER case default call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & @@ -646,9 +647,12 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) ! local real, dimension(SZK_(G)+1) :: Kd_tidal !< tidal diffusivity [m2/s] real, dimension(SZK_(G)+1) :: Kv_tidal !< tidal viscosity [m2/s] - real, dimension(SZK_(G)+1) :: vert_dep !< vertical deposition needed for Simmons tidal mixing. + real, dimension(SZK_(G)+1) :: vert_dep !< vertical deposition real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers (m) + + real, allocatable, dimension(:,:) :: exp_hab_zetar + integer :: i, k, is, ie real :: dh, hcorr, Simmons_coeff real, parameter :: rho_fw = 1000.0 ! fresh water density [kg/m^3] ! TODO: when coupled, get this from CESM (SHR_CONST_RHOFW) @@ -658,7 +662,7 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) dd => CS%dd select case (CS%cvmix_tidal_scheme) - case (SIMMONS_04) + case (SIMMONS) do i=is,ie if (G%mask2dT(i,j)<1) cycle @@ -722,7 +726,47 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) enddo ! i=is,ie - ! TODO: case (SCHMITTNER) + case (SCHMITTNER) + + allocate(exp_hab_zetar(G%ke+1,G%ke+1)) + + do i=is,ie + + if (G%mask2dT(i,j)<1) cycle + + iFaceHeight = 0.0 ! BBL is all relative to the surface + hcorr = 0.0 + do k=1,G%ke + ! cell center and cell bottom in meters (negative values in the ocean) + 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 + + ! form the time-invariant part of Schmittner coefficient term + call cvmix_compute_Schmittner_invariant(nlev = G%ke, & + VertDep = vert_dep, & + rho = rho_fw, & + exp_hab_zetar = exp_hab_zetar, & + zw = iFaceHeight, & + CVmix_tidal_params_user = CS%cvmix_tidal_params) + + ! form the Schmittner coefficient that is based on 3D q*E, which is formed from + ! summing q_i*TidalConstituent_i over the number of constituents. + !call cvmix_compute_SchmittnerCoeff( nlev = G%ke, & + ! energy_flux = , & + ! rho = rho_fw, & + ! SchmittnerCoeff = , & + ! exp_hab_zetar = , & + ! CVmix_tidal_params_user = CS%cvmix_tidal_params) + + enddo ! i=is,ie + + deallocate(exp_hab_zetar) + case default call MOM_error(FATAL, "tidal_mixing_init: Unrecognized setting "// & "#define CVMIX_TIDAL_SCHEME found in input file.") From 0e4dce5cd3feed3298ae85663bdbfe147dd35f0a Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 19 Apr 2018 10:00:09 -0600 Subject: [PATCH 04/18] add cvmix_compute_SchmittnerCoeff --- .../vertical/MOM_tidal_mixing.F90 | 81 ++++++++++--------- 1 file changed, 45 insertions(+), 36 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 2d8415d162..26d0bb3584 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -19,7 +19,7 @@ module MOM_tidal_mixing use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc use cvmix_tidal, only : cvmix_init_tidal, cvmix_compute_Simmons_invariant use cvmix_tidal, only : cvmix_coeffs_tidal, cvmix_tidal_params_type -use cvmix_tidal, only : cvmix_compute_Schmittner_invariant +use cvmix_tidal, only : cvmix_compute_Schmittner_invariant, cvmix_compute_SchmittnerCoeff use cvmix_kinds_and_types, only : cvmix_global_params_type use cvmix_put_get, only : cvmix_put @@ -138,13 +138,14 @@ module MOM_tidal_mixing real :: tidal_diss_lim_tc ! dissipation limit for tidal-energy-constituent data ! Data containers - real, pointer, dimension(:,:) :: TKE_Niku => NULL() - real, pointer, dimension(:,:) :: TKE_itidal => NULL() - real, pointer, dimension(:,:) :: Nb => NULL() - real, pointer, dimension(:,:) :: mask_itidal => NULL() - real, pointer, dimension(:,:) :: h2 => NULL() - real, pointer, dimension(:,:) :: tideamp => NULL() ! RMS tidal amplitude (m/s) - real, allocatable,dimension(:,:) :: tidal_qe_2d ! q*E(x,y) ! TODO: make this E(x,y) only + real, pointer, dimension(:,:) :: TKE_Niku => NULL() + real, pointer, dimension(:,:) :: TKE_itidal => NULL() + real, pointer, dimension(:,:) :: Nb => NULL() + real, pointer, dimension(:,:) :: mask_itidal => NULL() + real, pointer, dimension(:,:) :: h2 => NULL() + real, pointer, dimension(:,:) :: tideamp => NULL() ! RMS tidal amplitude (m/s) + real, allocatable ,dimension(:,:) :: tidal_qe_2d ! q*E(x,y) ! TODO: make this E(x,y) only + real, allocatable ,dimension(:,:,:) :: tidal_qe_3d_in ! q*E(x,y) ! TODO: make this E(x,y) only ! Diagnostics type(diag_ctrl), pointer :: diag => NULL() ! structure to regulate diagn output timing @@ -243,6 +244,11 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, "If true, use an internal tidal dissipation scheme to \n"//& "drive diapycnal mixing, along the lines of St. Laurent \n"//& "et al. (2002) and Simmons et al. (2004).", default=CS%use_cvmix_tidal) + + ! check if tidal mixing is active + tidal_mixing_init = CS%int_tide_dissipation + if (.not. tidal_mixing_init) return + if (CS%int_tide_dissipation) then ! Read in CVMix tidal scheme if CVMix tidal mixing is on @@ -651,6 +657,7 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers (m) + real, allocatable, dimension(:) :: Schmittner_coeff real, allocatable, dimension(:,:) :: exp_hab_zetar integer :: i, k, is, ie @@ -728,7 +735,11 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) case (SCHMITTNER) + ! TODO: correct exp_hab_zetar shapes in cvmix_compute_Schmittner_invariant + ! and cvmix_compute_SchmittnerCoeff low subroutines + allocate(exp_hab_zetar(G%ke+1,G%ke+1)) + allocate(Schmittner_coeff(G%ke)) do i=is,ie @@ -756,12 +767,12 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) ! form the Schmittner coefficient that is based on 3D q*E, which is formed from ! summing q_i*TidalConstituent_i over the number of constituents. - !call cvmix_compute_SchmittnerCoeff( nlev = G%ke, & - ! energy_flux = , & - ! rho = rho_fw, & - ! SchmittnerCoeff = , & - ! exp_hab_zetar = , & - ! CVmix_tidal_params_user = CS%cvmix_tidal_params) + call cvmix_compute_SchmittnerCoeff( nlev = G%ke, & + energy_flux = CS%tidal_qe_3d_in(i,j,:), & ! todo!!!: vertical interpolation + rho = rho_fw, & + SchmittnerCoeff = Schmittner_coeff, & + exp_hab_zetar = exp_hab_zetar, & + CVmix_tidal_params_user = CS%cvmix_tidal_params) enddo ! i=is,ie @@ -1336,20 +1347,20 @@ subroutine read_tidal_energy(G, tidal_energy_type, tidal_energy_file, CS) character(len=200), intent(in) :: tidal_energy_file type(tidal_mixing_cs), pointer :: CS ! local - integer :: isd, ied, jsd, jed + integer :: isd, ied, jsd, jed, nz real, allocatable, dimension(:,:) :: tidal_energy_flux_2d ! input tidal energy flux at T-grid points (W/m^2) - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - - if (.not. allocated(CS%tidal_qe_2d)) allocate(CS%tidal_qe_2d(isd:ied,jsd:jed)) + isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = G%ke select case (uppercase(tidal_energy_type(1:4))) case ('JAYN') ! Jayne 2009 + if (.not. allocated(CS%tidal_qe_2d)) allocate(CS%tidal_qe_2d(isd:ied,jsd:jed)) allocate(tidal_energy_flux_2d(isd:ied,jsd:jed)) call MOM_read_data(tidal_energy_file,'wave_dissipation',tidal_energy_flux_2d, G%domain) CS%tidal_qe_2d = (CS%Gamma_itides) * tidal_energy_flux_2d deallocate(tidal_energy_flux_2d) case ('ER03') ! Egbert & Ray 2003 + if (.not. allocated(CS%tidal_qe_3d_in)) allocate(CS%tidal_qe_3d_in(isd:ied,jsd:jed,nz)) call read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) case default call MOM_error(FATAL, "read_tidal_energy: Unknown tidal energy file type.") @@ -1377,16 +1388,14 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) tc_m2, & ! input lunar semidiurnal tidal energy flux [W/m^2] tc_s2, & ! input solar semidiurnal tidal energy flux [W/m^2] tc_k1, & ! input lunar diurnal tidal energy flux [W/m^2] - tc_o1, & ! input lunar diurnal tidal energy flux [W/m^2] - tidal_qe_3d ! sum_tc(q_tc*TC(x,y,z)) = q*E(x,y,z) + tc_o1 ! input lunar diurnal tidal energy flux [W/m^2] isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = G%ke allocate(tc_m2(isd:ied,jsd:jed,nz), & tc_s2(isd:ied,jsd:jed,nz), & tc_k1(isd:ied,jsd:jed,nz), & - tc_o1(isd:ied,jsd:jed,nz), & - tidal_qe_3d(isd:ied,jsd:jed,nz) ) + tc_o1(isd:ied,jsd:jed,nz) ) ! read in tidal constituents ! (NOTE: input z coordinates may differ from the model coordinates, which is fine.) @@ -1397,8 +1406,8 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) call MOM_read_data(tidal_energy_file, 'z_t', z_t) call MOM_read_data(tidal_energy_file, 'z_w', z_w) - ! form tidal_qe_3d from weighted tidal constituents - tidal_qe_3d = 0.0 + ! form tidal_qe_3d_in from weighted tidal constituents + CS%tidal_qe_3d_in = 0.0 where (abs(G%geoLatT(:,:)) < 30.0) tidal_qk1(:,:) = p33 @@ -1410,29 +1419,28 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) do k=1,nz where (z_t(k) <= G%bathyT(:,:) .and. z_w(k) > CS%tidal_diss_lim_tc) - tidal_qe_3d(:,:,k) = p33*tc_m2(:,:,k) + p33*tc_s2(:,:,k) + & + CS%tidal_qe_3d_in(:,:,k) = p33*tc_m2(:,:,k) + p33*tc_s2(:,:,k) + & tidal_qk1*tc_k1(:,:,k) + tidal_qo1*tc_o1(:,:,k) endwhere enddo ! test if qE is positive - if (any(tidal_qe_3d<0)) then - call MOM_error(FATAL, "read_tidal_constituents: Negative tidal_qe_3d terms.") + if (any(CS%tidal_qe_3d_in<0)) then + call MOM_error(FATAL, "read_tidal_constituents: Negative tidal_qe_3d_in terms.") endif - ! collapse 3D q*E to 2D q*E - CS%tidal_qe_2d = 0.0 - do k=1,nz - where (z_t(k) <= G%bathyT(:,:)) - CS%tidal_qe_2d(:,:) = CS%tidal_qe_2d(:,:) + tidal_qe_3d(:,:,k) - endwhere - enddo + !! collapse 3D q*E to 2D q*E + !CS%tidal_qe_2d = 0.0 + !do k=1,nz + ! where (z_t(k) <= G%bathyT(:,:)) + ! CS%tidal_qe_2d(:,:) = CS%tidal_qe_2d(:,:) + CS%tidal_qe_3d_in(:,:,k) + ! endwhere + !enddo deallocate(tc_m2) deallocate(tc_s2) deallocate(tc_k1) deallocate(tc_o1) - deallocate(tidal_qe_3d) end subroutine read_tidal_constituents @@ -1442,7 +1450,8 @@ subroutine tidal_mixing_end(CS) type(tidal_mixing_cs), pointer :: CS ! This module's control structure !TODO deallocate all the dynamically allocated members here ... - if (allocated(CS%tidal_qe_2d)) deallocate(CS%tidal_qe_2d) + if (allocated(CS%tidal_qe_2d)) deallocate(CS%tidal_qe_2d) + if (allocated(CS%tidal_qe_3d_in)) deallocate(CS%tidal_qe_3d_in) deallocate(CS%dd) deallocate(CS) From 9703e7c55d370bf16defdb9f70b3126ff262f96d Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 19 Apr 2018 18:16:41 -0600 Subject: [PATCH 05/18] remap tidal energy from input coord to model coord --- .../vertical/MOM_tidal_mixing.F90 | 89 ++++++++++++------- 1 file changed, 57 insertions(+), 32 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 4066f71c17..8ca8ecfffb 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -3,25 +3,26 @@ module MOM_tidal_mixing ! 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 : safe_alloc_ptr, post_data -use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag -use MOM_diag_to_Z, only : calc_Zint_diags -use MOM_EOS, only : calculate_density -use MOM_variables, only : thermo_var_ptrs, p3d -use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE -use MOM_debugging, only : hchksum -use MOM_grid, only : ocean_grid_type -use MOM_verticalGrid, only : verticalGrid_type -use MOM_file_parser, only : openParameterBlock, closeParameterBlock -use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_string_functions, only : uppercase, lowercase -use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc -use cvmix_tidal, only : cvmix_init_tidal, cvmix_compute_Simmons_invariant -use cvmix_tidal, only : cvmix_coeffs_tidal, cvmix_tidal_params_type -use cvmix_tidal, only : cvmix_compute_Schmittner_invariant, cvmix_compute_SchmittnerCoeff -use cvmix_kinds_and_types, only : cvmix_global_params_type -use cvmix_put_get, only : cvmix_put +use MOM_diag_mediator, only : diag_ctrl, time_type, register_diag_field +use MOM_diag_mediator, only : safe_alloc_ptr, post_data +use MOM_diag_to_Z, only : diag_to_Z_CS, register_Zint_diag +use MOM_diag_to_Z, only : calc_Zint_diags +use MOM_EOS, only : calculate_density +use MOM_variables, only : thermo_var_ptrs, p3d +use MOM_error_handler, only : MOM_error, is_root_pe, FATAL, WARNING, NOTE +use MOM_debugging, only : hchksum +use MOM_grid, only : ocean_grid_type +use MOM_verticalGrid, only : verticalGrid_type +use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h +use MOM_file_parser, only : openParameterBlock, closeParameterBlock +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type +use MOM_string_functions, only : uppercase, lowercase +use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc +use cvmix_tidal, only : cvmix_init_tidal, cvmix_compute_Simmons_invariant +use cvmix_tidal, only : cvmix_coeffs_tidal, cvmix_tidal_params_type +use cvmix_tidal, only : cvmix_compute_Schmittner_invariant, cvmix_compute_SchmittnerCoeff +use cvmix_kinds_and_types, only : cvmix_global_params_type +use cvmix_put_get, only : cvmix_put implicit none ; private @@ -136,6 +137,7 @@ module MOM_tidal_mixing type(cvmix_global_params_type) :: cvmix_glb_params ! for Prandtl number only real :: tidal_max_coef ! maximum allowable tidal diffusivity. [m^2/s] real :: tidal_diss_lim_tc ! dissipation limit for tidal-energy-constituent data + type(remapping_CS) :: remap_cs ! Data containers real, pointer, dimension(:,:) :: TKE_Niku => NULL() @@ -144,8 +146,9 @@ module MOM_tidal_mixing real, pointer, dimension(:,:) :: mask_itidal => NULL() real, pointer, dimension(:,:) :: h2 => NULL() real, pointer, dimension(:,:) :: tideamp => NULL() ! RMS tidal amplitude (m/s) - real, allocatable ,dimension(:,:) :: tidal_qe_2d ! q*E(x,y) ! TODO: make this E(x,y) only - real, allocatable ,dimension(:,:,:) :: tidal_qe_3d_in ! q*E(x,y) ! TODO: make this E(x,y) only + real, allocatable, dimension(:) :: h_src ! tidal constituent input layer thickness + real, allocatable ,dimension(:,:) :: tidal_qe_2d ! q*E(x,y) ! TODO: make this E(x,y) only + real, allocatable ,dimension(:,:,:) :: tidal_qe_3d_in ! q*E(x,y) ! Diagnostics type(diag_ctrl), pointer :: diag => NULL() ! structure to regulate diagn output timing @@ -660,13 +663,15 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) real, dimension(SZK_(G)+1) :: vert_dep !< vertical deposition real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers (m) - - real, allocatable, dimension(:) :: Schmittner_coeff + real, dimension(SZK_(G)) :: tidal_qe_md !< Tidal dissipation energy interpolated from 3d input to model coordinates + real, dimension(SZK_(G)) :: Schmittner_coeff + real, dimension(SZK_(G)) :: h_m !< Cell thickness [m] real, allocatable, dimension(:,:) :: exp_hab_zetar integer :: i, k, is, ie real :: dh, hcorr, Simmons_coeff real, parameter :: rho_fw = 1000.0 ! fresh water density [kg/m^3] ! TODO: when coupled, get this from CESM (SHR_CONST_RHOFW) + real :: h_neglect, h_neglect_edge type(tidal_mixing_diags), pointer :: dd is = G%isc ; ie = G%iec @@ -743,7 +748,12 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) ! and cvmix_compute_SchmittnerCoeff low subroutines allocate(exp_hab_zetar(G%ke+1,G%ke+1)) - allocate(Schmittner_coeff(G%ke)) + if (GV%Boussinesq) then + h_neglect = GV%m_to_H*1.0e-30 ; h_neglect_edge = GV%m_to_H*1.0e-10 + else + h_neglect = GV%kg_m2_to_H*1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H*1.0e-10 + endif + do i=is,ie @@ -751,9 +761,10 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) iFaceHeight = 0.0 ! BBL is all relative to the surface hcorr = 0.0 + h_m = h(i,j,:)*GV%H_to_m do k=1,G%ke ! cell center and cell bottom in meters (negative values in the ocean) - dh = h(i,j,k) * GV%H_to_m ! Nominal thickness to use for increment + dh = h_m(k) ! 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 @@ -769,10 +780,14 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) zw = iFaceHeight, & CVmix_tidal_params_user = CS%cvmix_tidal_params) + ! remap from input z coordinate to model coordinate: + tidal_qe_md = 0.0 + call remapping_core_h(CS%remap_cs, G%ke, CS%h_src, CS%tidal_qe_3d_in(i,j,:), G%ke, h_m, tidal_qe_md) + ! form the Schmittner coefficient that is based on 3D q*E, which is formed from ! summing q_i*TidalConstituent_i over the number of constituents. call cvmix_compute_SchmittnerCoeff( nlev = G%ke, & - energy_flux = CS%tidal_qe_3d_in(i,j,:), & ! todo!!!: vertical interpolation + energy_flux = tidal_qe_md(:), & rho = rho_fw, & SchmittnerCoeff = Schmittner_coeff, & exp_hab_zetar = exp_hab_zetar, & @@ -1364,7 +1379,6 @@ subroutine read_tidal_energy(G, tidal_energy_type, tidal_energy_file, CS) CS%tidal_qe_2d = (CS%Gamma_itides) * tidal_energy_flux_2d deallocate(tidal_energy_flux_2d) case ('ER03') ! Egbert & Ray 2003 - if (.not. allocated(CS%tidal_qe_3d_in)) allocate(CS%tidal_qe_3d_in(isd:ied,jsd:jed,nz)) call read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) case default call MOM_error(FATAL, "read_tidal_energy: Unknown tidal energy file type.") @@ -1383,8 +1397,8 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) integer :: k, isd, ied, jsd, jed, nz real, parameter :: p33 = 1.0/3.0 real, dimension(SZK_(G)) :: & - z_t, & ! depth from surface to midpoint of input layer - z_w ! depth from surface to top of input layer + z_t, & ! depth from surface to midpoint of input layer [cm] + z_w ! depth from surface to top of input layer [cm] real, dimension(SZI_(G),SZJ_(G)) :: & tidal_qk1, & ! qk1 coefficient used in Schmittner & Egbert tidal_qo1 ! qo1 coefficient used in Schmittner & Egbert @@ -1396,13 +1410,17 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = G%ke + ! allocate CS variables associated with 3d tidal energy dissipation + if (.not. allocated(CS%tidal_qe_3d_in)) allocate(CS%tidal_qe_3d_in(isd:ied,jsd:jed,nz)) + if (.not. allocated(CS%h_src)) allocate(CS%h_src(nz)) + + ! allocate local variables allocate(tc_m2(isd:ied,jsd:jed,nz), & tc_s2(isd:ied,jsd:jed,nz), & tc_k1(isd:ied,jsd:jed,nz), & tc_o1(isd:ied,jsd:jed,nz) ) ! read in tidal constituents - ! (NOTE: input z coordinates may differ from the model coordinates, which is fine.) call MOM_read_data(tidal_energy_file, 'M2', tc_m2, G%domain) call MOM_read_data(tidal_energy_file, 'S2', tc_s2, G%domain) call MOM_read_data(tidal_energy_file, 'K1', tc_k1, G%domain) @@ -1410,8 +1428,6 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) call MOM_read_data(tidal_energy_file, 'z_t', z_t) call MOM_read_data(tidal_energy_file, 'z_w', z_w) - ! form tidal_qe_3d_in from weighted tidal constituents - CS%tidal_qe_3d_in = 0.0 where (abs(G%geoLatT(:,:)) < 30.0) tidal_qk1(:,:) = p33 @@ -1421,7 +1437,11 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) tidal_qo1(:,:) = 1.0 endwhere + CS%tidal_qe_3d_in = 0.0 do k=1,nz + ! input cell thickness + CS%h_src(k) = (z_t(k)-z_w(k))*2.0 *1e-2 + ! form tidal_qe_3d_in from weighted tidal constituents where (z_t(k) <= G%bathyT(:,:) .and. z_w(k) > CS%tidal_diss_lim_tc) CS%tidal_qe_3d_in(:,:,k) = p33*tc_m2(:,:,k) + p33*tc_s2(:,:,k) + & tidal_qk1*tc_k1(:,:,k) + tidal_qo1*tc_o1(:,:,k) @@ -1441,6 +1461,10 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) ! endwhere !enddo + ! initialize input remapping: + call initialize_remapping(CS%remap_cs, remapping_scheme="PPM_IH4", & + boundary_extrapolation=.false., check_remapping=CS%debug) + deallocate(tc_m2) deallocate(tc_s2) deallocate(tc_k1) @@ -1456,6 +1480,7 @@ subroutine tidal_mixing_end(CS) !TODO deallocate all the dynamically allocated members here ... if (allocated(CS%tidal_qe_2d)) deallocate(CS%tidal_qe_2d) if (allocated(CS%tidal_qe_3d_in)) deallocate(CS%tidal_qe_3d_in) + if (allocated(CS%h_src)) deallocate(CS%h_src) deallocate(CS%dd) deallocate(CS) From 93f242f69f1221bd6fb22065f94de036232d41b6 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 24 Apr 2018 13:57:25 -0600 Subject: [PATCH 06/18] call cvmix_coeffs_tidal_schmittner --- .../vertical/MOM_tidal_mixing.F90 | 29 +++++++++++++++---- 1 file changed, 23 insertions(+), 6 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 8ca8ecfffb..ae958a02ed 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -21,6 +21,7 @@ module MOM_tidal_mixing use cvmix_tidal, only : cvmix_init_tidal, cvmix_compute_Simmons_invariant use cvmix_tidal, only : cvmix_coeffs_tidal, cvmix_tidal_params_type use cvmix_tidal, only : cvmix_compute_Schmittner_invariant, cvmix_compute_SchmittnerCoeff +use cvmix_tidal, only : cvmix_coeffs_tidal_schmittner use cvmix_kinds_and_types, only : cvmix_global_params_type use cvmix_put_get, only : cvmix_put @@ -36,7 +37,7 @@ module MOM_tidal_mixing !> Containers for tidal mixing diagnostics type, public :: tidal_mixing_diags - ! TODO: private + private real, pointer, dimension(:,:,:) :: & Kd_itidal => NULL(),& ! internal tide diffusivity at interfaces (m2/s) Fl_itidal => NULL(),& ! vertical flux of tidal turbulent dissipation (m3/s3) @@ -662,6 +663,7 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) real, dimension(SZK_(G)+1) :: Kv_tidal !< tidal viscosity [m2/s] real, dimension(SZK_(G)+1) :: vert_dep !< vertical deposition real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) + real, dimension(SZK_(G)+1) :: SchmittnerSocn real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers (m) real, dimension(SZK_(G)) :: tidal_qe_md !< Tidal dissipation energy interpolated from 3d input to model coordinates real, dimension(SZK_(G)) :: Schmittner_coeff @@ -772,6 +774,8 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) iFaceHeight(k+1) = iFaceHeight(k) - dh enddo + SchmittnerSocn = 0.0 ! TODO: compute this + ! form the time-invariant part of Schmittner coefficient term call cvmix_compute_Schmittner_invariant(nlev = G%ke, & VertDep = vert_dep, & @@ -786,11 +790,24 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) ! form the Schmittner coefficient that is based on 3D q*E, which is formed from ! summing q_i*TidalConstituent_i over the number of constituents. - call cvmix_compute_SchmittnerCoeff( nlev = G%ke, & - energy_flux = tidal_qe_md(:), & - rho = rho_fw, & - SchmittnerCoeff = Schmittner_coeff, & - exp_hab_zetar = exp_hab_zetar, & + call cvmix_compute_SchmittnerCoeff( nlev = G%ke, & + energy_flux = tidal_qe_md(:), & + rho = rho_fw, & + SchmittnerCoeff = Schmittner_coeff, & + exp_hab_zetar = exp_hab_zetar, & + CVmix_tidal_params_user = CS%cvmix_tidal_params) + + + call cvmix_coeffs_tidal_schmittner( Mdiff_out = Kv_tidal, & + Tdiff_out = Kd_tidal, & + Nsqr = N2_int(i,:), & + OceanDepth = -iFaceHeight(G%ke+1), & + vert_dep = vert_dep, & + nlev = G%ke, & + max_nlev = G%ke, & + SchmittnerCoeff = Schmittner_coeff, & + SchmittnerSouthernOcean = SchmittnerSocn, & + CVmix_params = CS%cvmix_glb_params, & CVmix_tidal_params_user = CS%cvmix_tidal_params) enddo ! i=is,ie From c04f455a153947dd4a9a14448e41e3293ee9810e Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 26 Apr 2018 17:03:16 -0600 Subject: [PATCH 07/18] add diagnostics for schmittner --- .../vertical/MOM_tidal_mixing.F90 | 74 ++++++++++++++----- 1 file changed, 54 insertions(+), 20 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index ae958a02ed..c0dd6e1796 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -39,18 +39,19 @@ module MOM_tidal_mixing type, public :: tidal_mixing_diags private real, pointer, dimension(:,:,:) :: & - Kd_itidal => NULL(),& ! internal tide diffusivity at interfaces (m2/s) - Fl_itidal => NULL(),& ! vertical flux of tidal turbulent dissipation (m3/s3) - Kd_lowmode => NULL(),& ! internal tide diffusivity at interfaces - ! due to propagating low modes (m2/s) (BDM) - Fl_lowmode => NULL(),& ! vertical flux of tidal turbulent dissipation - ! due to propagating low modes (m3/s3) (BDM) - Kd_Niku => NULL(),& ! lee-wave diffusivity at interfaces (m2/s) - Kd_Niku_work => NULL(),& ! layer integrated work by lee-wave driven mixing (W/m2) - Kd_Itidal_Work => NULL(),& ! layer integrated work by int tide driven mixing (W/m2) - Kd_Lowmode_Work => NULL(),& ! layer integrated work by low mode driven mixing (W/m2) BDM - N2_int => NULL(),& - vert_dep_3d => NULL() + Kd_itidal => NULL(),& ! internal tide diffusivity at interfaces (m2/s) + Fl_itidal => NULL(),& ! vertical flux of tidal turbulent dissipation (m3/s3) + Kd_lowmode => NULL(),& ! internal tide diffusivity at interfaces + ! due to propagating low modes (m2/s) (BDM) + Fl_lowmode => NULL(),& ! vertical flux of tidal turbulent dissipation + ! due to propagating low modes (m3/s3) (BDM) + Kd_Niku => NULL(),& ! lee-wave diffusivity at interfaces (m2/s) + Kd_Niku_work => NULL(),& ! layer integrated work by lee-wave driven mixing (W/m2) + Kd_Itidal_Work => NULL(),& ! layer integrated work by int tide driven mixing (W/m2) + Kd_Lowmode_Work => NULL(),& ! layer integrated work by low mode driven mixing (W/m2) BDM + N2_int => NULL(),& + vert_dep_3d => NULL(),& + Schmittner_coeff_3d => NULL() real, pointer, dimension(:,:) :: & TKE_itidal_used => NULL(),& ! internal tide TKE input at ocean bottom (W/m2) @@ -177,6 +178,7 @@ module MOM_tidal_mixing integer :: id_Polzin_decay_scale_scaled = -1 integer :: id_N2_int = -1 integer :: id_Simmons_coeff = -1 + integer :: id_Schmittner_coeff = -1 integer :: id_vert_dep = -1 end type tidal_mixing_cs @@ -501,11 +503,6 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, "The path to the file containing tidal energy \n"//& "dissipation. Used with CVMix tidal mixing schemes.", & fail_if_missing=.true.) - tidal_energy_file = trim(CS%inputdir) // trim(tidal_energy_file) - call get_param(param_file, mdl, "TIDAL_ENERGY_TYPE",tidal_energy_type, & - "The type of input tidal energy flux dataset.",& - fail_if_missing=.true.) - ! TODO: list all available tidal energy types here call get_param(param_file, mdl, 'MIN_THICKNESS', CS%min_thickness, default=0.001, & do_not_log=.True.) call get_param(param_file, mdl, "PRANDTL_TIDAL", prandtl_tidal, & @@ -515,11 +512,22 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, do_not_log=.true.) call cvmix_put(CS%cvmix_glb_params,'Prandtl',prandtl_tidal) + tidal_energy_file = trim(CS%inputdir) // trim(tidal_energy_file) + call get_param(param_file, mdl, "TIDAL_ENERGY_TYPE",tidal_energy_type, & + "The type of input tidal energy flux dataset. Valid values are"//& + "\t Jayne\n"//& + "\t ER03 \n",& + fail_if_missing=.true.) + ! Check whether tidal energy input format and CVMix tidal mixing scheme are consistent + if ( .not. ( & + (uppercase(tidal_energy_type(1:4)).eq.'JAYN' .and. CS%cvmix_tidal_scheme.eq.SIMMONS).or. & + (uppercase(tidal_energy_type(1:4)).eq.'ER03' .and. CS%cvmix_tidal_scheme.eq.SCHMITTNER) ) )then + call MOM_error(FATAL, "tidal_mixing_init: Tidal energy file type ("//& + trim(tidal_energy_type)//") is incompatible with CVMix tidal "//& + " mixing scheme: "//trim(cvmix_tidal_scheme_str) ) + endif cvmix_tidal_scheme_str = lowercase(cvmix_tidal_scheme_str) - - ! TODO: check parameter consistency. (see POP::tidal_mixing.F90::tidal_check) - ! Set up CVMix call cvmix_init_tidal(CVmix_tidal_params_user = CS%cvmix_tidal_params, & mix_scheme = cvmix_tidal_scheme_str, & @@ -549,6 +557,8 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, ! TODO: add units CS%id_Simmons_coeff = register_diag_field('ocean_model','Simmons_coeff',diag%axesT1,Time, & 'time-invariant portion of the tidal mixing coefficient using the Simmons', '') + CS%id_Schmittner_coeff = register_diag_field('ocean_model','Schmittner_coeff',diag%axesTi,Time, & + 'time-invariant portion of the tidal mixing coefficient using the Schmittner', '') CS%id_vert_dep = register_diag_field('ocean_model','vert_dep',diag%axesTi,Time, & 'vertical deposition function needed for Simmons et al tidal mixing', '') @@ -810,6 +820,25 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) CVmix_params = CS%cvmix_glb_params, & CVmix_tidal_params_user = CS%cvmix_tidal_params) + do k=1,G%ke + Kd(i,j,k) = Kd(i,j,k) + 0.5*(Kd_tidal(k) + Kd_tidal(k+1) ) + !TODO: Kv(i,j,k) = ???????????? + enddo + + ! diagnostics + ! diagnostics + if (associated(dd%Kd_itidal)) then + dd%Kd_itidal(i,j,:) = Kd_tidal(:) + endif + if (associated(dd%N2_int)) then + dd%N2_int(i,j,:) = N2_int(i,:) + endif + if (associated(dd%Schmittner_coeff_3d)) then + dd%Schmittner_coeff_3d(i,j,:) = Schmittner_coeff(:) + endif + if (associated(dd%vert_dep_3d)) then + dd%vert_dep_3d(i,j,:) = vert_dep(:) + endif enddo ! i=is,ie deallocate(exp_hab_zetar) @@ -1288,6 +1317,9 @@ subroutine setup_tidal_diagnostics(G,CS) if (CS%id_vert_dep > 0) then allocate(dd%vert_dep_3d(isd:ied,jsd:jed,nz+1)) ; dd%vert_dep_3d(:,:,:) = 0.0 endif + if (CS%id_Schmittner_coeff > 0) then + allocate(dd%Schmittner_coeff_3d(isd:ied,jsd:jed,nz)) ; dd%Schmittner_coeff_3d(:,:,:) = 0.0 + endif end subroutine setup_tidal_diagnostics subroutine post_tidal_diagnostics(G,GV,h,CS) @@ -1322,6 +1354,7 @@ subroutine post_tidal_diagnostics(G,GV,h,CS) if (CS%id_N2_int> 0) call post_data(CS%id_N2_int, dd%N2_int, CS%diag) if (CS%id_vert_dep> 0) call post_data(CS%id_vert_dep, dd%vert_dep_3d, CS%diag) if (CS%id_Simmons_coeff> 0) call post_data(CS%id_Simmons_coeff, dd%Simmons_coeff_2d, CS%diag) + if (CS%id_Schmittner_coeff> 0) call post_data(CS%id_Schmittner_coeff, dd%Schmittner_coeff_3d, CS%diag) if (CS%id_Kd_Itidal_Work > 0) & call post_data(CS%id_Kd_Itidal_Work, dd%Kd_Itidal_Work, CS%diag) @@ -1373,6 +1406,7 @@ subroutine post_tidal_diagnostics(G,GV,h,CS) if (associated(dd%N2_int)) deallocate(dd%N2_int) if (associated(dd%vert_dep_3d)) deallocate(dd%vert_dep_3d) if (associated(dd%Simmons_coeff_2d)) deallocate(dd%Simmons_coeff_2d) + if (associated(dd%Schmittner_coeff_3d)) deallocate(dd%Schmittner_coeff_3d) end subroutine post_tidal_diagnostics From 05a2e5985425c7bccdf406786c0094525e0a7f2b Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 2 May 2018 14:53:02 -0600 Subject: [PATCH 08/18] debug 3d tidal energy remapping --- .../vertical/MOM_tidal_mixing.F90 | 101 +++++++++++++----- 1 file changed, 73 insertions(+), 28 deletions(-) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index c0dd6e1796..37b187878d 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -17,7 +17,7 @@ module MOM_tidal_mixing use MOM_file_parser, only : openParameterBlock, closeParameterBlock use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_string_functions, only : uppercase, lowercase -use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc +use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc, field_size use cvmix_tidal, only : cvmix_init_tidal, cvmix_compute_Simmons_invariant use cvmix_tidal, only : cvmix_coeffs_tidal, cvmix_tidal_params_type use cvmix_tidal, only : cvmix_compute_Schmittner_invariant, cvmix_compute_SchmittnerCoeff @@ -51,7 +51,8 @@ module MOM_tidal_mixing Kd_Lowmode_Work => NULL(),& ! layer integrated work by low mode driven mixing (W/m2) BDM N2_int => NULL(),& vert_dep_3d => NULL(),& - Schmittner_coeff_3d => NULL() + Schmittner_coeff_3d => NULL(),& + tidal_qe_md => NULL() real, pointer, dimension(:,:) :: & TKE_itidal_used => NULL(),& ! internal tide TKE input at ocean bottom (W/m2) @@ -147,8 +148,8 @@ module MOM_tidal_mixing real, pointer, dimension(:,:) :: Nb => NULL() real, pointer, dimension(:,:) :: mask_itidal => NULL() real, pointer, dimension(:,:) :: h2 => NULL() - real, pointer, dimension(:,:) :: tideamp => NULL() ! RMS tidal amplitude (m/s) - real, allocatable, dimension(:) :: h_src ! tidal constituent input layer thickness + real, pointer, dimension(:,:) :: tideamp => NULL() !< RMS tidal amplitude [m/s] + real, allocatable, dimension(:) :: h_src !< tidal constituent input layer thickness [m] real, allocatable ,dimension(:,:) :: tidal_qe_2d ! q*E(x,y) ! TODO: make this E(x,y) only real, allocatable ,dimension(:,:,:) :: tidal_qe_3d_in ! q*E(x,y) @@ -179,6 +180,7 @@ module MOM_tidal_mixing integer :: id_N2_int = -1 integer :: id_Simmons_coeff = -1 integer :: id_Schmittner_coeff = -1 + integer :: id_tidal_qe_md = -1 integer :: id_vert_dep = -1 end type tidal_mixing_cs @@ -557,8 +559,10 @@ logical function tidal_mixing_init(Time, G, GV, param_file, diag, diag_to_Z_CSp, ! TODO: add units CS%id_Simmons_coeff = register_diag_field('ocean_model','Simmons_coeff',diag%axesT1,Time, & 'time-invariant portion of the tidal mixing coefficient using the Simmons', '') - CS%id_Schmittner_coeff = register_diag_field('ocean_model','Schmittner_coeff',diag%axesTi,Time, & + CS%id_Schmittner_coeff = register_diag_field('ocean_model','Schmittner_coeff',diag%axesTL,Time, & 'time-invariant portion of the tidal mixing coefficient using the Schmittner', '') + CS%id_tidal_qe_md = register_diag_field('ocean_model','tidal_qe_md',diag%axesTL,Time, & + 'input tidal energy dissipated locally interpolated to model vertical coordinates', '') CS%id_vert_dep = register_diag_field('ocean_model','vert_dep',diag%axesTi,Time, & 'vertical deposition function needed for Simmons et al tidal mixing', '') @@ -673,7 +677,7 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) real, dimension(SZK_(G)+1) :: Kv_tidal !< tidal viscosity [m2/s] real, dimension(SZK_(G)+1) :: vert_dep !< vertical deposition real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m) - real, dimension(SZK_(G)+1) :: SchmittnerSocn + real, dimension(SZK_(G)+1) :: SchmittnerSocn real, dimension(SZK_(G)) :: cellHeight !< Height of cell centers (m) real, dimension(SZK_(G)) :: tidal_qe_md !< Tidal dissipation energy interpolated from 3d input to model coordinates real, dimension(SZK_(G)) :: Schmittner_coeff @@ -796,7 +800,8 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) ! remap from input z coordinate to model coordinate: tidal_qe_md = 0.0 - call remapping_core_h(CS%remap_cs, G%ke, CS%h_src, CS%tidal_qe_3d_in(i,j,:), G%ke, h_m, tidal_qe_md) + call remapping_core_h(CS%remap_cs, size(CS%h_src), CS%h_src, CS%tidal_qe_3d_in(i,j,:), & + G%ke, h_m, tidal_qe_md) ! form the Schmittner coefficient that is based on 3D q*E, which is formed from ! summing q_i*TidalConstituent_i over the number of constituents. @@ -811,7 +816,7 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) call cvmix_coeffs_tidal_schmittner( Mdiff_out = Kv_tidal, & Tdiff_out = Kd_tidal, & Nsqr = N2_int(i,:), & - OceanDepth = -iFaceHeight(G%ke+1), & + OceanDepth = -iFaceHeight(G%ke+1), & vert_dep = vert_dep, & nlev = G%ke, & max_nlev = G%ke, & @@ -825,7 +830,6 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) !TODO: Kv(i,j,k) = ???????????? enddo - ! diagnostics ! diagnostics if (associated(dd%Kd_itidal)) then dd%Kd_itidal(i,j,:) = Kd_tidal(:) @@ -836,6 +840,9 @@ subroutine calculate_cvmix_tidal(h, j, G, GV, CS, N2_int, Kd) if (associated(dd%Schmittner_coeff_3d)) then dd%Schmittner_coeff_3d(i,j,:) = Schmittner_coeff(:) endif + if (associated(dd%tidal_qe_md)) then + dd%tidal_qe_md(i,j,:) = tidal_qe_md(:) + endif if (associated(dd%vert_dep_3d)) then dd%vert_dep_3d(i,j,:) = vert_dep(:) endif @@ -1312,14 +1319,29 @@ subroutine setup_tidal_diagnostics(G,CS) allocate(dd%N2_int(isd:ied,jsd:jed,nz+1)) ; dd%N2_int(:,:,:) = 0.0 endif if (CS%id_Simmons_coeff > 0) then + if (CS%cvmix_tidal_scheme .ne. SIMMONS) then + call MOM_error(FATAL, "setup_tidal_diagnostics: Simmons_coeff diagnostics is available "//& + "only when cvmix_tidal_scheme is Simmons") + endif allocate(dd%Simmons_coeff_2d(isd:ied,jsd:jed)) ; dd%Simmons_coeff_2d(:,:) = 0.0 endif if (CS%id_vert_dep > 0) then allocate(dd%vert_dep_3d(isd:ied,jsd:jed,nz+1)) ; dd%vert_dep_3d(:,:,:) = 0.0 endif if (CS%id_Schmittner_coeff > 0) then + if (CS%cvmix_tidal_scheme .ne. SCHMITTNER) then + call MOM_error(FATAL, "setup_tidal_diagnostics: Schmittner_coeff diagnostics is available "//& + "only when cvmix_tidal_scheme is Schmittner.") + endif allocate(dd%Schmittner_coeff_3d(isd:ied,jsd:jed,nz)) ; dd%Schmittner_coeff_3d(:,:,:) = 0.0 endif + if (CS%id_tidal_qe_md > 0) then + if (CS%cvmix_tidal_scheme .ne. SCHMITTNER) then + call MOM_error(FATAL, "setup_tidal_diagnostics: tidal_qe_md diagnostics is available "//& + "only when cvmix_tidal_scheme is Schmittner.") + endif + allocate(dd%tidal_qe_md(isd:ied,jsd:jed,nz)) ; dd%tidal_qe_md(:,:,:) = 0.0 + endif end subroutine setup_tidal_diagnostics subroutine post_tidal_diagnostics(G,GV,h,CS) @@ -1355,6 +1377,7 @@ subroutine post_tidal_diagnostics(G,GV,h,CS) if (CS%id_vert_dep> 0) call post_data(CS%id_vert_dep, dd%vert_dep_3d, CS%diag) if (CS%id_Simmons_coeff> 0) call post_data(CS%id_Simmons_coeff, dd%Simmons_coeff_2d, CS%diag) if (CS%id_Schmittner_coeff> 0) call post_data(CS%id_Schmittner_coeff, dd%Schmittner_coeff_3d, CS%diag) + if (CS%id_tidal_qe_md> 0) call post_data(CS%id_tidal_qe_md, dd%tidal_qe_md, CS%diag) if (CS%id_Kd_Itidal_Work > 0) & call post_data(CS%id_Kd_Itidal_Work, dd%Kd_Itidal_Work, CS%diag) @@ -1407,6 +1430,7 @@ subroutine post_tidal_diagnostics(G,GV,h,CS) if (associated(dd%vert_dep_3d)) deallocate(dd%vert_dep_3d) if (associated(dd%Simmons_coeff_2d)) deallocate(dd%Simmons_coeff_2d) if (associated(dd%Schmittner_coeff_3d)) deallocate(dd%Schmittner_coeff_3d) + if (associated(dd%tidal_qe_md)) deallocate(dd%tidal_qe_md) end subroutine post_tidal_diagnostics @@ -1445,31 +1469,36 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) type(tidal_mixing_cs), pointer :: CS ! local - integer :: k, isd, ied, jsd, jed, nz - real, parameter :: p33 = 1.0/3.0 - real, dimension(SZK_(G)) :: & - z_t, & ! depth from surface to midpoint of input layer [cm] - z_w ! depth from surface to top of input layer [cm] + integer :: k, isd, ied, jsd, jed, i,j + integer, dimension(4) :: nz_in + real, parameter :: p33 = 1.0/3.0 real, dimension(SZI_(G),SZJ_(G)) :: & tidal_qk1, & ! qk1 coefficient used in Schmittner & Egbert tidal_qo1 ! qo1 coefficient used in Schmittner & Egbert + real, allocatable, dimension(:) :: & + z_t, & ! depth from surface to midpoint of input layer [cm] + z_w ! depth from surface to top of input layer [cm] real, allocatable, dimension(:,:,:) :: & tc_m2, & ! input lunar semidiurnal tidal energy flux [W/m^2] tc_s2, & ! input solar semidiurnal tidal energy flux [W/m^2] tc_k1, & ! input lunar diurnal tidal energy flux [W/m^2] tc_o1 ! input lunar diurnal tidal energy flux [W/m^2] - 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 - ! allocate CS variables associated with 3d tidal energy dissipation - if (.not. allocated(CS%tidal_qe_3d_in)) allocate(CS%tidal_qe_3d_in(isd:ied,jsd:jed,nz)) - if (.not. allocated(CS%h_src)) allocate(CS%h_src(nz)) + ! get number of input levels: + call field_size(tidal_energy_file, 'z_t',nz_in) ! allocate local variables - allocate(tc_m2(isd:ied,jsd:jed,nz), & - tc_s2(isd:ied,jsd:jed,nz), & - tc_k1(isd:ied,jsd:jed,nz), & - tc_o1(isd:ied,jsd:jed,nz) ) + allocate(z_t(nz_in(1)), z_w(nz_in(1)) ) + allocate(tc_m2(isd:ied,jsd:jed,nz_in(1)), & + tc_s2(isd:ied,jsd:jed,nz_in(1)), & + tc_k1(isd:ied,jsd:jed,nz_in(1)), & + tc_o1(isd:ied,jsd:jed,nz_in(1)) ) + + ! allocate CS variables associated with 3d tidal energy dissipation + if (.not. allocated(CS%tidal_qe_3d_in)) allocate(CS%tidal_qe_3d_in(isd:ied,jsd:jed,nz_in(1))) + if (.not. allocated(CS%h_src)) allocate(CS%h_src(nz_in(1))) ! read in tidal constituents call MOM_read_data(tidal_energy_file, 'M2', tc_m2, G%domain) @@ -1479,7 +1508,6 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) call MOM_read_data(tidal_energy_file, 'z_t', z_t) call MOM_read_data(tidal_energy_file, 'z_w', z_w) - where (abs(G%geoLatT(:,:)) < 30.0) tidal_qk1(:,:) = p33 tidal_qo1(:,:) = p33 @@ -1489,37 +1517,54 @@ subroutine read_tidal_constituents(G, tidal_energy_type, tidal_energy_file, CS) endwhere CS%tidal_qe_3d_in = 0.0 - do k=1,nz + do k=1,nz_in(1) ! input cell thickness CS%h_src(k) = (z_t(k)-z_w(k))*2.0 *1e-2 ! form tidal_qe_3d_in from weighted tidal constituents - where (z_t(k) <= G%bathyT(:,:) .and. z_w(k) > CS%tidal_diss_lim_tc) + where ( (z_t(k)*1e-2) <= G%bathyT(:,:) .and. (z_w(k)*1e-2) > CS%tidal_diss_lim_tc) CS%tidal_qe_3d_in(:,:,k) = p33*tc_m2(:,:,k) + p33*tc_s2(:,:,k) + & tidal_qk1*tc_k1(:,:,k) + tidal_qo1*tc_o1(:,:,k) endwhere enddo + !open(unit=1905,file="out_1905.txt",access="APPEND") + !do j=G%jsd,G%jed + ! do i=isd,ied + ! if ( i+G%idg_offset .eq. 90 .and. j+G%jdg_offset .eq. 126) then + ! print *, "-------------------------------------------" + ! do k=50,nz_in(1) + ! write(1905,*) i,j,k + ! write(1905,*) CS%tidal_qe_3d_in(i,j,k), tc_m2(i,j,k) + ! write(1905,*) z_t(k), G%bathyT(i,j), z_w(k),CS%tidal_diss_lim_tc + ! end do + ! endif + ! enddo + !enddo + !close(1905) + ! test if qE is positive - if (any(CS%tidal_qe_3d_in<0)) then + if (any(CS%tidal_qe_3d_in<0.0)) then call MOM_error(FATAL, "read_tidal_constituents: Negative tidal_qe_3d_in terms.") endif !! collapse 3D q*E to 2D q*E !CS%tidal_qe_2d = 0.0 - !do k=1,nz + !do k=1,nz_in(1) ! where (z_t(k) <= G%bathyT(:,:)) ! CS%tidal_qe_2d(:,:) = CS%tidal_qe_2d(:,:) + CS%tidal_qe_3d_in(:,:,k) ! endwhere !enddo ! initialize input remapping: - call initialize_remapping(CS%remap_cs, remapping_scheme="PPM_IH4", & + call initialize_remapping(CS%remap_cs, remapping_scheme="PLM", & boundary_extrapolation=.false., check_remapping=CS%debug) deallocate(tc_m2) deallocate(tc_s2) deallocate(tc_k1) deallocate(tc_o1) + deallocate(z_t) + deallocate(z_w) end subroutine read_tidal_constituents From 9ce97ef0519246ffa826217df223793e69f5d4f4 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 2 May 2018 16:58:47 -0600 Subject: [PATCH 09/18] use unmodified CVMix_compute_Schmittner_invariant interface --- src/parameterizations/vertical/MOM_tidal_mixing.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/parameterizations/vertical/MOM_tidal_mixing.F90 b/src/parameterizations/vertical/MOM_tidal_mixing.F90 index 11effe5dca..9321d3f68c 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -797,10 +797,14 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) ! form the time-invariant part of Schmittner coefficient term call CVMix_compute_Schmittner_invariant(nlev = G%ke, & VertDep = vert_dep, & + efficiency = CS%Mu_itides, & rho = rho_fw, & exp_hab_zetar = exp_hab_zetar, & zw = iFaceHeight, & CVmix_tidal_params_user = CS%CVMix_tidal_params) + !TODO: in above call, there is no need to pass efficiency, since it gets + ! passed via CVMix_init_tidal and stored in CVMix_tidal_params. Change + ! CVMix API to prevent this redundancy. ! remap from input z coordinate to model coordinate: tidal_qe_md = 0.0 From 72774d72b87f8511ef54e26761d5d2cf41a72871 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 8 May 2018 10:32:28 -0600 Subject: [PATCH 10/18] split OBL depth computation and KPP_calculate --- src/parameterizations/vertical/MOM_KPP.F90 | 116 +++++++++++++++------ 1 file changed, 85 insertions(+), 31 deletions(-) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index 697cc26125..48c683f60a 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -406,9 +406,8 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive) end function KPP_init - -!> KPP vertical diffusivity/viscosity and non-local tracer transport -subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & +!> Compute OBL depth +subroutine KPP_compute_OBL(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,& nonLocalTransScalar) @@ -434,21 +433,19 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local transport (m/s) ! Local variables - integer :: i, j, k, km1,kp1 ! Loop indices + integer :: i, j, k, km1 ! Loop indices real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface (m) (negative in ocean) real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface (m) (negative in ocean) real, dimension( G%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces (1/s2) real, dimension( G%ke+1 ) :: N_1d ! Brunt-Vaisala frequency at interfaces (1/s) (floored at 0) real, dimension( G%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars (m/s) - real, dimension( G%ke ) :: Wm_1d ! Profile of vertical velocity scale for momentum (m/s) + !real, dimension( G%ke ) :: Wm_1d ! Profile of vertical velocity scale for momentum (m/s) real, dimension( G%ke ) :: Vt2_1d ! Unresolved velocity for bulk Ri calculation/diagnostic (m2/s2) - real, dimension( G%ke ) :: BulkRi_1d ! Bulk Richardson number for each layer real, dimension( G%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number real, dimension( G%ke ) :: deltaU2 ! square of delta U (shear) in denominator of Bulk Ri (m2/s2) - real, dimension( G%ke+1, 2) :: Kdiffusivity ! Vertical diffusivity at interfaces (m2/s) - real, dimension( G%ke+1 ) :: Kviscosity ! Vertical viscosity at interfaces (m2/s) real, dimension( G%ke+1, 2) :: nonLocalTrans ! Non-local transport for heat/salt at interfaces (non-dimensional) real, dimension( G%ke ) :: surfBuoyFlux2 + real, dimension( G%ke ) :: BulkRi_1d ! Bulk Richardson number for each layer ! for EOS calculation real, dimension( 3*G%ke ) :: rho_1D @@ -457,7 +454,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & real, dimension( 3*G%ke ) :: Salt_1D real :: kOBL, OBLdepth_0d, surfFricVel, surfBuoyFlux, Coriolis - real :: GoRho, pRef, rho1, rhoK, rhoKm1, Uk, Vk, sigma + real :: GoRho, pRef, rho1, rhoK, Uk, Vk, sigma real :: zBottomMinusOffset ! Height of bottom plus a little bit (m) real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth. @@ -471,20 +468,6 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & real :: hcorr ! A cumulative correction arising from inflation of vanished layers (m) integer :: kk, ksfc, ktmp -#ifdef __DO_SAFETY_CHECKS__ - if (CS%debug) then - call hchksum(h, "KPP in: h",G%HI,haloshift=0, scale=GV%H_to_m) - call hchksum(Temp, "KPP in: T",G%HI,haloshift=0) - call hchksum(Salt, "KPP in: S",G%HI,haloshift=0) - call hchksum(u, "KPP in: u",G%HI,haloshift=0) - call hchksum(v, "KPP in: v",G%HI,haloshift=0) - call hchksum(uStar, "KPP in: uStar",G%HI,haloshift=0) - call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0) - call hchksum(Kt, "KPP in: Kt",G%HI,haloshift=0) - call hchksum(Ks, "KPP in: Ks",G%HI,haloshift=0) - endif -#endif - ! some constants GoRho = GV%g_Earth / GV%Rho0 nonLocalTrans(:,:) = 0.0 @@ -492,17 +475,15 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & if (CS%id_Kd_in > 0) call post_data(CS%id_Kd_in, Kt, CS%diag) !$OMP parallel do default(none) shared(G,GV,CS,EOS,uStar,Temp,Salt,u,v,h,GoRho, & -!$OMP buoyFlux, nonLocalTransHeat, & -!$OMP nonLocalTransScalar,Kt,Ks,Kv) & -!$OMP firstprivate(nonLocalTrans) & +!$OMP buoyFlux) & !$OMP private(Coriolis,surfFricVel,SLdepth_0d,hTot,surfTemp, & !$OMP surfHtemp,surfSalt,surfHsalt,surfU, & !$OMP surfHu,surfV,surfHv,iFaceHeight, & !$OMP pRef,km1,cellHeight,Uk,Vk,deltaU2, & -!$OMP rho1,rhoK,rhoKm1,deltaRho,N2_1d,N_1d,delH, & -!$OMP surfBuoyFlux,Ws_1d,Vt2_1d,BulkRi_1d, & -!$OMP OBLdepth_0d,zBottomMinusOffset,Kdiffusivity, & -!$OMP Kviscosity,sigma,kOBL,kk,pres_1D,Temp_1D, & +!$OMP rho1,rhoK,deltaRho,N2_1d,N_1d,delH, & +!$OMP surfBuoyFlux,Ws_1d,BulkRi_1d, & +!$OMP OBLdepth_0d,zBottomMinusOffset, & +!$OMP sigma,kOBL,kk,pres_1D,Temp_1D, & !$OMP Salt_1D,rho_1D,surfBuoyFlux2,ksfc,dh,hcorr) ! loop over horizontal points on processor @@ -746,6 +727,79 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & ! smg: remove code above ! ********************************************************************** + enddo + enddo + +end subroutine + +!> KPP vertical diffusivity/viscosity and non-local tracer transport +subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & + buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,& + nonLocalTransScalar) + + ! Arguments + type(KPP_CS), pointer :: CS !< Control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses (units of H) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: Temp !< potential/cons temp (deg C) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: Salt !< Salinity (ppt) + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Velocity i-component (m/s) + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Velocity j-component (m/s) + type(EOS_type), pointer :: EOS !< Equation of state + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity (m/s) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: buoyFlux !< Surface buoyancy flux (m2/s3) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kt !< (in) Vertical diffusivity of heat w/o KPP (m2/s) + !< (out) Vertical diffusivity including KPP (m2/s) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Ks !< (in) Vertical diffusivity of salt w/o KPP (m2/s) + !< (out) Vertical diffusivity including KPP (m2/s) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kv !< (in) Vertical viscosity w/o KPP (m2/s) + !< (out) Vertical viscosity including KPP (m2/s) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransHeat !< Temp non-local transport (m/s) + real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local transport (m/s) + + ! Local variables + integer :: i, j, k, km1,kp1 ! Loop indices + real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface (m) (negative in ocean) + real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface (m) (negative in ocean) + real, dimension( G%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces (1/s2) + real, dimension( G%ke+1 ) :: N_1d ! Brunt-Vaisala frequency at interfaces (1/s) (floored at 0) + real, dimension( G%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars (m/s) + !real, dimension( G%ke ) :: Wm_1d ! Profile of vertical velocity scale for momentum (m/s) + real, dimension( G%ke ) :: Vt2_1d ! Unresolved velocity for bulk Ri calculation/diagnostic (m2/s2) + real, dimension( G%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number + real, dimension( G%ke ) :: deltaU2 ! square of delta U (shear) in denominator of Bulk Ri (m2/s2) + real, dimension( G%ke+1, 2) :: nonLocalTrans ! Non-local transport for heat/salt at interfaces (non-dimensional) + + real, dimension( G%ke ) :: BulkRi_1d ! Bulk Richardson number for each layer + real, dimension( G%ke+1, 2) :: Kdiffusivity ! Vertical diffusivity at interfaces (m2/s) + real, dimension( G%ke+1 ) :: Kviscosity ! Vertical viscosity at interfaces (m2/s) + + real :: kOBL, OBLdepth_0d, surfFricVel, surfBuoyFlux + real :: sigma + + real :: surfTemp ! Integral and average of temp over the surface layer + real :: surfSalt ! Integral and average of saln over the surface layer + real :: surfU ! Integral and average of u over the surface layer + real :: surfV ! Integral and average of v over the surface layer + +#ifdef __DO_SAFETY_CHECKS__ + if (CS%debug) then + call hchksum(h, "KPP in: h",G%HI,haloshift=0, scale=GV%H_to_m) + call hchksum(Temp, "KPP in: T",G%HI,haloshift=0) + call hchksum(Salt, "KPP in: S",G%HI,haloshift=0) + call hchksum(u, "KPP in: u",G%HI,haloshift=0) + call hchksum(v, "KPP in: v",G%HI,haloshift=0) + call hchksum(uStar, "KPP in: uStar",G%HI,haloshift=0) + call hchksum(buoyFlux, "KPP in: buoyFlux",G%HI,haloshift=0) + call hchksum(Kt, "KPP in: Kt",G%HI,haloshift=0) + call hchksum(Ks, "KPP in: Ks",G%HI,haloshift=0) + endif +#endif + + ! loop over horizontal points on processor + do j = G%jsc, G%jec + do i = G%isc, G%iec ! Call CVMix/KPP to obtain OBL diffusivities, viscosities and non-local transports @@ -880,7 +934,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & if (CS%id_Tsurf > 0) CS%Tsurf(i,j) = surfTemp if (CS%id_Ssurf > 0) CS%Ssurf(i,j) = surfSalt if (CS%id_Usurf > 0) CS%Usurf(i,j) = surfU - if (CS%id_Vsurf > 0) CS%Vsurf(i,j) = surfv + if (CS%id_Vsurf > 0) CS%Vsurf(i,j) = surfV ! Update output of routine From 80d93232cc447c1db452c80388c18076bf81a25e Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 8 May 2018 15:27:38 -0600 Subject: [PATCH 11/18] restore KPP_calculate for now --- src/parameterizations/vertical/MOM_KPP.F90 | 297 ++++++++++++++++++++- 1 file changed, 285 insertions(+), 12 deletions(-) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index 48c683f60a..192ae02389 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -732,6 +732,7 @@ subroutine KPP_compute_OBL(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & end subroutine + !> KPP vertical diffusivity/viscosity and non-local tracer transport subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,& @@ -758,30 +759,43 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransHeat !< Temp non-local transport (m/s) real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local transport (m/s) - ! Local variables +! Local variables integer :: i, j, k, km1,kp1 ! Loop indices real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface (m) (negative in ocean) real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface (m) (negative in ocean) real, dimension( G%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces (1/s2) real, dimension( G%ke+1 ) :: N_1d ! Brunt-Vaisala frequency at interfaces (1/s) (floored at 0) real, dimension( G%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars (m/s) - !real, dimension( G%ke ) :: Wm_1d ! Profile of vertical velocity scale for momentum (m/s) + real, dimension( G%ke ) :: Wm_1d ! Profile of vertical velocity scale for momentum (m/s) real, dimension( G%ke ) :: Vt2_1d ! Unresolved velocity for bulk Ri calculation/diagnostic (m2/s2) + real, dimension( G%ke ) :: BulkRi_1d ! Bulk Richardson number for each layer real, dimension( G%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number real, dimension( G%ke ) :: deltaU2 ! square of delta U (shear) in denominator of Bulk Ri (m2/s2) - real, dimension( G%ke+1, 2) :: nonLocalTrans ! Non-local transport for heat/salt at interfaces (non-dimensional) - - real, dimension( G%ke ) :: BulkRi_1d ! Bulk Richardson number for each layer real, dimension( G%ke+1, 2) :: Kdiffusivity ! Vertical diffusivity at interfaces (m2/s) real, dimension( G%ke+1 ) :: Kviscosity ! Vertical viscosity at interfaces (m2/s) + real, dimension( G%ke+1, 2) :: nonLocalTrans ! Non-local transport for heat/salt at interfaces (non-dimensional) + real, dimension( G%ke ) :: surfBuoyFlux2 - real :: kOBL, OBLdepth_0d, surfFricVel, surfBuoyFlux - real :: sigma + ! for EOS calculation + real, dimension( 3*G%ke ) :: rho_1D + real, dimension( 3*G%ke ) :: pres_1D + real, dimension( 3*G%ke ) :: Temp_1D + real, dimension( 3*G%ke ) :: Salt_1D - real :: surfTemp ! Integral and average of temp over the surface layer - real :: surfSalt ! Integral and average of saln over the surface layer - real :: surfU ! Integral and average of u over the surface layer - real :: surfV ! Integral and average of v over the surface layer + real :: kOBL, OBLdepth_0d, surfFricVel, surfBuoyFlux, Coriolis + real :: GoRho, pRef, rho1, rhoK, rhoKm1, Uk, Vk, sigma + + real :: zBottomMinusOffset ! Height of bottom plus a little bit (m) + real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth. + real :: hTot ! Running sum of thickness used in the surface layer average (m) + real :: delH ! Thickness of a layer (m) + real :: surfHtemp, surfTemp ! Integral and average of temp over the surface layer + real :: surfHsalt, surfSalt ! Integral and average of saln over the surface layer + real :: surfHu, surfU ! Integral and average of u over the surface layer + real :: surfHv, surfV ! Integral and average of v over the surface layer + real :: dh ! The local thickness used for calculating interface positions (m) + real :: hcorr ! A cumulative correction arising from inflation of vanished layers (m) + integer :: kk, ksfc, ktmp #ifdef __DO_SAFETY_CHECKS__ if (CS%debug) then @@ -797,10 +811,268 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & endif #endif + ! some constants + GoRho = GV%g_Earth / GV%Rho0 + nonLocalTrans(:,:) = 0.0 + + if (CS%id_Kd_in > 0) call post_data(CS%id_Kd_in, Kt, CS%diag) + +!$OMP parallel do default(none) shared(G,GV,CS,EOS,uStar,Temp,Salt,u,v,h,GoRho, & +!$OMP buoyFlux, nonLocalTransHeat, & +!$OMP nonLocalTransScalar,Kt,Ks,Kv) & +!$OMP firstprivate(nonLocalTrans) & +!$OMP private(Coriolis,surfFricVel,SLdepth_0d,hTot,surfTemp, & +!$OMP surfHtemp,surfSalt,surfHsalt,surfU, & +!$OMP surfHu,surfV,surfHv,iFaceHeight, & +!$OMP pRef,km1,cellHeight,Uk,Vk,deltaU2, & +!$OMP rho1,rhoK,rhoKm1,deltaRho,N2_1d,N_1d,delH, & +!$OMP surfBuoyFlux,Ws_1d,Vt2_1d,BulkRi_1d, & +!$OMP OBLdepth_0d,zBottomMinusOffset,Kdiffusivity, & +!$OMP Kviscosity,sigma,kOBL,kk,pres_1D,Temp_1D, & +!$OMP Salt_1D,rho_1D,surfBuoyFlux2,ksfc,dh,hcorr) + ! loop over horizontal points on processor do j = G%jsc, G%jec do i = G%isc, G%iec + ! skip calling KPP for land points + if (G%mask2dT(i,j)==0.) cycle + + ! things independent of position within the column + Coriolis = 0.25*( (G%CoriolisBu(i,j) + G%CoriolisBu(i-1,j-1)) & + +(G%CoriolisBu(i-1,j) + G%CoriolisBu(i,j-1)) ) + surfFricVel = uStar(i,j) + + ! Bullk Richardson number computed for each cell in a column, + ! assuming OBLdepth = grid cell depth. After Rib(k) is + ! known for the column, then CVMix interpolates to find + ! the actual OBLdepth. This approach avoids need to iterate + ! on the OBLdepth calculation. It follows that used in MOM5 + ! and POP. + iFaceHeight(1) = 0.0 ! BBL is all relative to the surface + pRef = 0. + hcorr = 0. + do k=1,G%ke + + ! cell center and cell bottom in meters (negative values in the ocean) + 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 + + ! find ksfc for cell where "surface layer" sits + SLdepth_0d = CS%surf_layer_ext*max( max(-cellHeight(k),-iFaceHeight(2) ), CS%minOBLdepth ) + ksfc = k + do ktmp = 1,k + if (-1.0*iFaceHeight(ktmp+1) >= SLdepth_0d) then + ksfc = ktmp + exit + endif + enddo + + ! average temp, saln, u, v over surface layer + ! use C-grid average to get u,v on T-points. + surfHtemp=0.0 + surfHsalt=0.0 + surfHu =0.0 + surfHv =0.0 + hTot =0.0 + do ktmp = 1,ksfc + + ! SLdepth_0d can be between cell interfaces + delH = min( max(0.0, SLdepth_0d - hTot), h(i,j,ktmp)*GV%H_to_m ) + + ! surface layer thickness + hTot = hTot + delH + + ! surface averaged fields + surfHtemp = surfHtemp + Temp(i,j,ktmp) * delH + surfHsalt = surfHsalt + Salt(i,j,ktmp) * delH + surfHu = surfHu + 0.5*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delH + surfHv = surfHv + 0.5*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delH + + enddo + surfTemp = surfHtemp / hTot + surfSalt = surfHsalt / hTot + surfU = surfHu / hTot + surfV = surfHv / hTot + + ! vertical shear between present layer and + ! surface layer averaged surfU,surfV. + ! C-grid average to get Uk and Vk on T-points. + Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU + Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV + deltaU2(k) = Uk**2 + Vk**2 + + ! pressure, temp, and saln for EOS + ! kk+1 = surface fields + ! kk+2 = k fields + ! kk+3 = km1 fields + km1 = max(1, k-1) + kk = 3*(k-1) + pres_1D(kk+1) = pRef + pres_1D(kk+2) = pRef + pres_1D(kk+3) = pRef + Temp_1D(kk+1) = surfTemp + Temp_1D(kk+2) = Temp(i,j,k) + Temp_1D(kk+3) = Temp(i,j,km1) + Salt_1D(kk+1) = surfSalt + Salt_1D(kk+2) = Salt(i,j,k) + Salt_1D(kk+3) = Salt(i,j,km1) + + ! pRef is pressure at interface between k and km1. + ! iterate pRef for next pass through k-loop. + pRef = pRef + GV%H_to_Pa * h(i,j,k) + + ! this difference accounts for penetrating SW + surfBuoyFlux2(k) = buoyFlux(i,j,1) - buoyFlux(i,j,k+1) + + enddo ! k-loop finishes + + ! compute in-situ density + call calculate_density(Temp_1D, Salt_1D, pres_1D, rho_1D, 1, 3*G%ke, EOS) + + ! N2 (can be negative) and N (non-negative) on interfaces. + ! deltaRho is non-local rho difference used for bulk Richardson number. + ! N_1d is local N (with floor) used for unresolved shear calculation. + do k = 1, G%ke + km1 = max(1, k-1) + kk = 3*(k-1) + deltaRho(k) = rho_1D(kk+2) - rho_1D(kk+1) + N2_1d(k) = (GoRho * (rho_1D(kk+2) - rho_1D(kk+3)) ) / & + ((0.5*(h(i,j,km1) + h(i,j,k))+GV%H_subroundoff)*GV%H_to_m) + N_1d(k) = sqrt( max( N2_1d(k), 0.) ) + enddo + N2_1d(G%ke+1 ) = 0.0 + N_1d(G%ke+1 ) = 0.0 + + ! turbulent velocity scales w_s and w_m computed at the cell centers. + ! Note that if sigma > CS%surf_layer_ext, then CVMix_kpp_compute_turbulent_scales + ! computes w_s and w_m velocity scale at sigma=CS%surf_layer_ext. So we only pass + ! sigma=CS%surf_layer_ext for this calculation. + call CVMix_kpp_compute_turbulent_scales( & + CS%surf_layer_ext, & ! (in) Normalized surface layer depth; sigma = CS%surf_layer_ext + -cellHeight, & ! (in) Assume here that OBL depth (m) = -cellHeight(k) + surfBuoyFlux2, & ! (in) Buoyancy flux at surface (m2/s3) + surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) + w_s=Ws_1d, & ! (out) Turbulent velocity scale profile (m/s) + CVMix_kpp_params_user=CS%KPP_params ) + + ! Calculate Bulk Richardson number from eq (21) of LMD94 + BulkRi_1d = CVMix_kpp_compute_bulk_Richardson( & + cellHeight(1:G%ke), & ! Depth of cell center (m) + GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) (1/s) + deltaU2, & ! Square of resolved velocity difference (m2/s2) + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile (m/s) + N_iface=N_1d) ! Buoyancy frequency (1/s) + + + surfBuoyFlux = buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit + ! h to Monin-Obukov (default is false, ie. not used) + + call CVMix_kpp_compute_OBL_depth( & + BulkRi_1d, & ! (in) Bulk Richardson number + iFaceHeight, & ! (in) Height of interfaces (m) + OBLdepth_0d, & ! (out) OBL depth (m) + kOBL, & ! (out) level (+fraction) of OBL extent + zt_cntr=cellHeight, & ! (in) Height of cell centers (m) + surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) + surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) + Coriolis=Coriolis, & ! (in) Coriolis parameter (1/s) + CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters + + ! A hack to avoid KPP reaching the bottom. It was needed during development + ! because KPP was unable to handle vanishingly small layers near the bottom. + if (CS%deepOBLoffset>0.) then + zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1)) + OBLdepth_0d = min( OBLdepth_0d, -zBottomMinusOffset ) + endif + + ! apply some constraints on OBLdepth + if(CS%fixedOBLdepth) OBLdepth_0d = CS%fixedOBLdepth_value + OBLdepth_0d = max( OBLdepth_0d, -iFaceHeight(2) ) ! no shallower than top layer + OBLdepth_0d = min( OBLdepth_0d, -iFaceHeight(G%ke+1) ) ! no deeper than bottom + kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) + +!************************************************************************* +! smg: remove code below + +! Following "correction" step has been found to be unnecessary. +! Code should be removed after further testing. + if (CS%correctSurfLayerAvg) then + + SLdepth_0d = CS%surf_layer_ext * OBLdepth_0d + hTot = h(i,j,1) + surfTemp = Temp(i,j,1) ; surfHtemp = surfTemp * hTot + surfSalt = Salt(i,j,1) ; surfHsalt = surfSalt * hTot + surfU = 0.5*(u(i,j,1)+u(i-1,j,1)) ; surfHu = surfU * hTot + surfV = 0.5*(v(i,j,1)+v(i,j-1,1)) ; surfHv = surfV * hTot + pRef = 0.0 + + do k = 2, G%ke + + ! Recalculate differences with surface layer + Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU + Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV + deltaU2(k) = Uk**2 + Vk**2 + pRef = pRef + GV%H_to_Pa * h(i,j,k) + call calculate_density(surfTemp, surfSalt, pRef, rho1, EOS) + call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS) + deltaRho(k) = rhoK - rho1 + + ! Surface layer averaging (needed for next k+1 iteration of this loop) + if (hTot < SLdepth_0d) then + delH = min( max(0., SLdepth_0d - hTot), h(i,j,k)*GV%H_to_m ) + hTot = hTot + delH + surfHtemp = surfHtemp + Temp(i,j,k) * delH ; surfTemp = surfHtemp / hTot + surfHsalt = surfHsalt + Salt(i,j,k) * delH ; surfSalt = surfHsalt / hTot + surfHu = surfHu + 0.5*(u(i,j,k)+u(i-1,j,k)) * delH ; surfU = surfHu / hTot + surfHv = surfHv + 0.5*(v(i,j,k)+v(i,j-1,k)) * delH ; surfV = surfHv / hTot + endif + + enddo + + BulkRi_1d = CVMix_kpp_compute_bulk_Richardson( & + cellHeight(1:G%ke), & ! Depth of cell center (m) + GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) (1/s) + deltaU2, & ! Square of resolved velocity difference (m2/s2) + ws_cntr=Ws_1d, & ! Turbulent velocity scale profile (m/s) + N_iface=N_1d ) ! Buoyancy frequency (1/s) + + surfBuoyFlux = buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit + ! h to Monin-Obukov (default is false, ie. not used) + + call CVMix_kpp_compute_OBL_depth( & + BulkRi_1d, & ! (in) Bulk Richardson number + iFaceHeight, & ! (in) Height of interfaces (m) + OBLdepth_0d, & ! (out) OBL depth (m) + kOBL, & ! (out) level (+fraction) of OBL extent + zt_cntr=cellHeight, & ! (in) Height of cell centers (m) + surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) + surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) + Coriolis=Coriolis, & ! (in) Coriolis parameter (1/s) + CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters + + if (CS%deepOBLoffset>0.) then + zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1)) + OBLdepth_0d = min( OBLdepth_0d, -zBottomMinusOffset ) + kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) + endif + + ! apply some constraints on OBLdepth + if(CS%fixedOBLdepth) OBLdepth_0d = CS%fixedOBLdepth_value + OBLdepth_0d = max( OBLdepth_0d, -iFaceHeight(2) ) ! no shallower than top layer + OBLdepth_0d = min( OBLdepth_0d, -iFaceHeight(G%ke+1) ) ! no deep than bottom + kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) + + endif ! endif for "correction" step + +! smg: remove code above +! ********************************************************************** + + ! Call CVMix/KPP to obtain OBL diffusivities, viscosities and non-local transports ! Unlike LMD94, we do not match to interior diffusivities. If using the original @@ -934,7 +1206,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & if (CS%id_Tsurf > 0) CS%Tsurf(i,j) = surfTemp if (CS%id_Ssurf > 0) CS%Ssurf(i,j) = surfSalt if (CS%id_Usurf > 0) CS%Usurf(i,j) = surfU - if (CS%id_Vsurf > 0) CS%Vsurf(i,j) = surfV + if (CS%id_Vsurf > 0) CS%Vsurf(i,j) = surfv ! Update output of routine @@ -991,6 +1263,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & end subroutine KPP_calculate + !> Copies KPP surface boundary layer depth into BLD subroutine KPP_get_BLD(CS, BLD, G) type(KPP_CS), pointer :: CS !< Control structure for From 85810d2c2afcc2f19a7709f6623a1dcd75da5103 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 8 May 2018 15:32:51 -0600 Subject: [PATCH 12/18] remove unnecessary KPP_compute_OBL arguments --- src/parameterizations/vertical/MOM_KPP.F90 | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index 192ae02389..a9e9c06f87 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -407,9 +407,7 @@ end function KPP_init !> Compute OBL depth -subroutine KPP_compute_OBL(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & - buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,& - nonLocalTransScalar) +subroutine KPP_compute_OBL(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) ! Arguments type(KPP_CS), pointer :: CS !< Control structure @@ -423,14 +421,6 @@ subroutine KPP_compute_OBL(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & type(EOS_type), pointer :: EOS !< Equation of state real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity (m/s) real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: buoyFlux !< Surface buoyancy flux (m2/s3) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kt !< (in) Vertical diffusivity of heat w/o KPP (m2/s) - !< (out) Vertical diffusivity including KPP (m2/s) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Ks !< (in) Vertical diffusivity of salt w/o KPP (m2/s) - !< (out) Vertical diffusivity including KPP (m2/s) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kv !< (in) Vertical viscosity w/o KPP (m2/s) - !< (out) Vertical viscosity including KPP (m2/s) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransHeat !< Temp non-local transport (m/s) - real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local transport (m/s) ! Local variables integer :: i, j, k, km1 ! Loop indices @@ -472,8 +462,6 @@ subroutine KPP_compute_OBL(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & GoRho = GV%g_Earth / GV%Rho0 nonLocalTrans(:,:) = 0.0 - if (CS%id_Kd_in > 0) call post_data(CS%id_Kd_in, Kt, CS%diag) - !$OMP parallel do default(none) shared(G,GV,CS,EOS,uStar,Temp,Salt,u,v,h,GoRho, & !$OMP buoyFlux) & !$OMP private(Coriolis,surfFricVel,SLdepth_0d,hTot,surfTemp, & @@ -730,7 +718,7 @@ subroutine KPP_compute_OBL(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & enddo enddo -end subroutine +end subroutine KPP_compute_OBL !> KPP vertical diffusivity/viscosity and non-local tracer transport From 82d9228b2927478463834ecd4b1e48d5ff933fe7 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 8 May 2018 16:33:13 -0600 Subject: [PATCH 13/18] rename KPP_compute_OBL as KPP_compute_BLD --- src/parameterizations/vertical/MOM_KPP.F90 | 9 +++++++-- src/parameterizations/vertical/MOM_diabatic_driver.F90 | 6 +++++- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index a9e9c06f87..83f9788418 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -28,6 +28,7 @@ module MOM_KPP #include "MOM_memory.h" public :: KPP_init +public :: KPP_compute_BLD public :: KPP_calculate public :: KPP_end public :: KPP_NonLocalTransport_temp @@ -171,6 +172,10 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive) 'If False, calculates the non-local transport and tendencies but\n'//& 'purely for diagnostic purposes.', & default=.not. CS%passiveMode) + call get_param(paramFile, mdl, 'APPLY_KPP_OBL_FILTER', CS%applyNonLocalTrans, & + 'If True, applies a 1-1-4-1-1 Laplacian filter one time on HBLT.\n'// & + 'computed via CVMix to reduce any horizontal two-grid-point noise.', & + default=.false.) call get_param(paramFile, mdl, 'RI_CRIT', CS%Ri_crit, & 'Critical bulk Richardson number used to define depth of the\n'// & 'surface Ocean Boundary Layer (OBL).', & @@ -407,7 +412,7 @@ end function KPP_init !> Compute OBL depth -subroutine KPP_compute_OBL(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) +subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) ! Arguments type(KPP_CS), pointer :: CS !< Control structure @@ -718,7 +723,7 @@ subroutine KPP_compute_OBL(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) enddo enddo -end subroutine KPP_compute_OBL +end subroutine KPP_compute_BLD !> KPP vertical diffusivity/viscosity and non-local tracer transport diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 02b9896ab7..284b209932 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -49,7 +49,8 @@ module MOM_diabatic_driver use MOM_internal_tides, only : propagate_int_tide use MOM_internal_tides, only : internal_tides_init, internal_tides_end, int_tide_CS use MOM_kappa_shear, only : kappa_shear_is_used -use MOM_KPP, only : KPP_CS, KPP_init, KPP_calculate, KPP_end, KPP_get_BLD +use MOM_KPP, only : KPP_CS, KPP_init, KPP_compute_BLD, KPP_calculate +use MOM_KPP, only : KPP_end, KPP_get_BLD use MOM_KPP, only : KPP_NonLocalTransport_temp, KPP_NonLocalTransport_saln use MOM_opacity, only : opacity_init, set_opacity, opacity_end, opacity_CS use MOM_regularize_layers, only : regularize_layers, regularize_layers_init, regularize_layers_CS @@ -635,6 +636,9 @@ subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, G, G endif !$OMP end parallel + call KPP_compute_BLD(CS%KPP_CSp, G, GV, h, tv%T, tv%S, u, v, tv%eqn_of_state, & + fluxes%ustar, CS%KPP_buoy_flux) + call KPP_calculate(CS%KPP_CSp, G, GV, h, tv%T, tv%S, u, v, tv%eqn_of_state, & fluxes%ustar, CS%KPP_buoy_flux, Kd_heat, Kd_salt, visc%Kv_shear, CS%KPP_NLTheat, CS%KPP_NLTscalar) !$OMP parallel default(none) shared(is,ie,js,je,nz,Kd_salt,Kd_int,visc,CS,Kd_heat) From c2a6ed841167d648985044b0827084cef6699d13 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 9 May 2018 20:55:18 -0600 Subject: [PATCH 14/18] add smoothBLD var --- src/parameterizations/vertical/MOM_KPP.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index 83f9788418..167bd2b589 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -70,6 +70,7 @@ module MOM_KPP character(len=30) :: MatchTechnique !< Method used in CVMix for setting diffusivity and NLT profile functions integer :: NLT_shape !< MOM6 over-ride of CVMix NLT shape function logical :: applyNonLocalTrans !< If True, apply non-local transport to heat and scalars + logical :: smoothBLD !< If True, apply a 1-1-4-1-1 Laplacian filter one time on HBLT. logical :: KPPzeroDiffusivity !< If True, will set diffusivity and viscosity from KPP to zero; for testing purposes. logical :: KPPisAdditive !< If True, will add KPP diffusivity to initial diffusivity. !! If False, will replace initial diffusivity wherever KPP diffusivity is non-zero. @@ -172,7 +173,7 @@ logical function KPP_init(paramFile, G, diag, Time, CS, passive) 'If False, calculates the non-local transport and tendencies but\n'//& 'purely for diagnostic purposes.', & default=.not. CS%passiveMode) - call get_param(paramFile, mdl, 'APPLY_KPP_OBL_FILTER', CS%applyNonLocalTrans, & + call get_param(paramFile, mdl, 'SMOOTH_BLD', CS%smoothBLD, & 'If True, applies a 1-1-4-1-1 Laplacian filter one time on HBLT.\n'// & 'computed via CVMix to reduce any horizontal two-grid-point noise.', & default=.false.) From 34e4cf4fc32dc0dba03c7ffd36d4820332881baf Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 10 May 2018 13:57:32 -0600 Subject: [PATCH 15/18] use 2d CS arrays for OBLdepth and kOBL --- src/parameterizations/vertical/MOM_KPP.F90 | 173 +++++---------------- 1 file changed, 38 insertions(+), 135 deletions(-) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index 167bd2b589..3ae2a1a196 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -106,6 +106,7 @@ module MOM_KPP ! Diagnostics arrays real, allocatable, dimension(:,:) :: OBLdepth !< Depth (positive) of OBL (m) + real, allocatable, dimension(:,:) :: kOBL !< Level (+fraction) of OBL extent real, allocatable, dimension(:,:,:) :: dRho !< Bulk difference in density (kg/m3) real, allocatable, dimension(:,:,:) :: Uz2 !< Square of bulk difference in resolved velocity (m2/s2) real, allocatable, dimension(:,:,:) :: BulkRi !< Bulk Richardson number for each layer (dimensionless) @@ -376,8 +377,11 @@ 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. + allocate( CS%OBLdepth( SZI_(G), SZJ_(G) ) ) + CS%OBLdepth(:,:) = 0. + allocate( CS%kOBL( SZI_(G), SZJ_(G) ) ) + CS%kOBL(:,:) = 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) ) ) @@ -449,7 +453,7 @@ subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) real, dimension( 3*G%ke ) :: Temp_1D real, dimension( 3*G%ke ) :: Salt_1D - real :: kOBL, OBLdepth_0d, surfFricVel, surfBuoyFlux, Coriolis + real :: surfFricVel, surfBuoyFlux, Coriolis real :: GoRho, pRef, rho1, rhoK, Uk, Vk, sigma real :: zBottomMinusOffset ! Height of bottom plus a little bit (m) @@ -476,8 +480,8 @@ subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) !$OMP pRef,km1,cellHeight,Uk,Vk,deltaU2, & !$OMP rho1,rhoK,deltaRho,N2_1d,N_1d,delH, & !$OMP surfBuoyFlux,Ws_1d,BulkRi_1d, & -!$OMP OBLdepth_0d,zBottomMinusOffset, & -!$OMP sigma,kOBL,kk,pres_1D,Temp_1D, & +!$OMP zBottomMinusOffset, & +!$OMP sigma,kk,pres_1D,Temp_1D, & !$OMP Salt_1D,rho_1D,surfBuoyFlux2,ksfc,dh,hcorr) ! loop over horizontal points on processor @@ -624,8 +628,8 @@ subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) call CVMix_kpp_compute_OBL_depth( & BulkRi_1d, & ! (in) Bulk Richardson number iFaceHeight, & ! (in) Height of interfaces (m) - OBLdepth_0d, & ! (out) OBL depth (m) - kOBL, & ! (out) level (+fraction) of OBL extent + CS%OBLdepth(i,j), & ! (out) OBL depth (m) + CS%kOBL(i,j), & ! (out) level (+fraction) of OBL extent zt_cntr=cellHeight, & ! (in) Height of cell centers (m) surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) @@ -636,14 +640,14 @@ subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) ! because KPP was unable to handle vanishingly small layers near the bottom. if (CS%deepOBLoffset>0.) then zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1)) - OBLdepth_0d = min( OBLdepth_0d, -zBottomMinusOffset ) + CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset ) endif ! apply some constraints on OBLdepth - if(CS%fixedOBLdepth) OBLdepth_0d = CS%fixedOBLdepth_value - OBLdepth_0d = max( OBLdepth_0d, -iFaceHeight(2) ) ! no shallower than top layer - OBLdepth_0d = min( OBLdepth_0d, -iFaceHeight(G%ke+1) ) ! no deeper than bottom - kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) + if(CS%fixedOBLdepth) CS%OBLdepth(i,j) = CS%fixedOBLdepth_value + CS%OBLdepth(i,j) = max( CS%OBLdepth(i,j), -iFaceHeight(2) ) ! no shallower than top layer + CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(G%ke+1) ) ! no deeper than bottom + CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) !************************************************************************* ! smg: remove code below @@ -652,7 +656,7 @@ subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) ! Code should be removed after further testing. if (CS%correctSurfLayerAvg) then - SLdepth_0d = CS%surf_layer_ext * OBLdepth_0d + SLdepth_0d = CS%surf_layer_ext * CS%OBLdepth(i,j) hTot = h(i,j,1) surfTemp = Temp(i,j,1) ; surfHtemp = surfTemp * hTot surfSalt = Salt(i,j,1) ; surfHsalt = surfSalt * hTot @@ -696,8 +700,8 @@ subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) call CVMix_kpp_compute_OBL_depth( & BulkRi_1d, & ! (in) Bulk Richardson number iFaceHeight, & ! (in) Height of interfaces (m) - OBLdepth_0d, & ! (out) OBL depth (m) - kOBL, & ! (out) level (+fraction) of OBL extent + CS%OBLdepth(i,j), & ! (out) OBL depth (m) + CS%kOBL(i,j), & ! (out) level (+fraction) of OBL extent zt_cntr=cellHeight, & ! (in) Height of cell centers (m) surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) @@ -706,15 +710,15 @@ subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) if (CS%deepOBLoffset>0.) then zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1)) - OBLdepth_0d = min( OBLdepth_0d, -zBottomMinusOffset ) - kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) + CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset ) + CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) endif ! apply some constraints on OBLdepth - if(CS%fixedOBLdepth) OBLdepth_0d = CS%fixedOBLdepth_value - OBLdepth_0d = max( OBLdepth_0d, -iFaceHeight(2) ) ! no shallower than top layer - OBLdepth_0d = min( OBLdepth_0d, -iFaceHeight(G%ke+1) ) ! no deep than bottom - kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) + if(CS%fixedOBLdepth) CS%OBLdepth(i,j) = CS%fixedOBLdepth_value + CS%OBLdepth(i,j) = max( CS%OBLdepth(i,j), -iFaceHeight(2) ) ! no shallower than top layer + CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(G%ke+1) ) ! no deep than bottom + CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) endif ! endif for "correction" step @@ -776,7 +780,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & real, dimension( 3*G%ke ) :: Temp_1D real, dimension( 3*G%ke ) :: Salt_1D - real :: kOBL, OBLdepth_0d, surfFricVel, surfBuoyFlux, Coriolis + real :: surfFricVel, surfBuoyFlux, Coriolis real :: GoRho, pRef, rho1, rhoK, rhoKm1, Uk, Vk, sigma real :: zBottomMinusOffset ! Height of bottom plus a little bit (m) @@ -821,8 +825,8 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & !$OMP pRef,km1,cellHeight,Uk,Vk,deltaU2, & !$OMP rho1,rhoK,rhoKm1,deltaRho,N2_1d,N_1d,delH, & !$OMP surfBuoyFlux,Ws_1d,Vt2_1d,BulkRi_1d, & -!$OMP OBLdepth_0d,zBottomMinusOffset,Kdiffusivity, & -!$OMP Kviscosity,sigma,kOBL,kk,pres_1D,Temp_1D, & +!$OMP zBottomMinusOffset,Kdiffusivity, & +!$OMP Kviscosity,sigma,kk,pres_1D,Temp_1D, & !$OMP Salt_1D,rho_1D,surfBuoyFlux2,ksfc,dh,hcorr) ! loop over horizontal points on processor @@ -966,106 +970,6 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & surfBuoyFlux = buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit ! h to Monin-Obukov (default is false, ie. not used) - call CVMix_kpp_compute_OBL_depth( & - BulkRi_1d, & ! (in) Bulk Richardson number - iFaceHeight, & ! (in) Height of interfaces (m) - OBLdepth_0d, & ! (out) OBL depth (m) - kOBL, & ! (out) level (+fraction) of OBL extent - zt_cntr=cellHeight, & ! (in) Height of cell centers (m) - surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) - surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) - Coriolis=Coriolis, & ! (in) Coriolis parameter (1/s) - CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters - - ! A hack to avoid KPP reaching the bottom. It was needed during development - ! because KPP was unable to handle vanishingly small layers near the bottom. - if (CS%deepOBLoffset>0.) then - zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1)) - OBLdepth_0d = min( OBLdepth_0d, -zBottomMinusOffset ) - endif - - ! apply some constraints on OBLdepth - if(CS%fixedOBLdepth) OBLdepth_0d = CS%fixedOBLdepth_value - OBLdepth_0d = max( OBLdepth_0d, -iFaceHeight(2) ) ! no shallower than top layer - OBLdepth_0d = min( OBLdepth_0d, -iFaceHeight(G%ke+1) ) ! no deeper than bottom - kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) - -!************************************************************************* -! smg: remove code below - -! Following "correction" step has been found to be unnecessary. -! Code should be removed after further testing. - if (CS%correctSurfLayerAvg) then - - SLdepth_0d = CS%surf_layer_ext * OBLdepth_0d - hTot = h(i,j,1) - surfTemp = Temp(i,j,1) ; surfHtemp = surfTemp * hTot - surfSalt = Salt(i,j,1) ; surfHsalt = surfSalt * hTot - surfU = 0.5*(u(i,j,1)+u(i-1,j,1)) ; surfHu = surfU * hTot - surfV = 0.5*(v(i,j,1)+v(i,j-1,1)) ; surfHv = surfV * hTot - pRef = 0.0 - - do k = 2, G%ke - - ! Recalculate differences with surface layer - Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU - Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV - deltaU2(k) = Uk**2 + Vk**2 - pRef = pRef + GV%H_to_Pa * h(i,j,k) - call calculate_density(surfTemp, surfSalt, pRef, rho1, EOS) - call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS) - deltaRho(k) = rhoK - rho1 - - ! Surface layer averaging (needed for next k+1 iteration of this loop) - if (hTot < SLdepth_0d) then - delH = min( max(0., SLdepth_0d - hTot), h(i,j,k)*GV%H_to_m ) - hTot = hTot + delH - surfHtemp = surfHtemp + Temp(i,j,k) * delH ; surfTemp = surfHtemp / hTot - surfHsalt = surfHsalt + Salt(i,j,k) * delH ; surfSalt = surfHsalt / hTot - surfHu = surfHu + 0.5*(u(i,j,k)+u(i-1,j,k)) * delH ; surfU = surfHu / hTot - surfHv = surfHv + 0.5*(v(i,j,k)+v(i,j-1,k)) * delH ; surfV = surfHv / hTot - endif - - enddo - - BulkRi_1d = CVMix_kpp_compute_bulk_Richardson( & - cellHeight(1:G%ke), & ! Depth of cell center (m) - GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) (1/s) - deltaU2, & ! Square of resolved velocity difference (m2/s2) - ws_cntr=Ws_1d, & ! Turbulent velocity scale profile (m/s) - N_iface=N_1d ) ! Buoyancy frequency (1/s) - - surfBuoyFlux = buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit - ! h to Monin-Obukov (default is false, ie. not used) - - call CVMix_kpp_compute_OBL_depth( & - BulkRi_1d, & ! (in) Bulk Richardson number - iFaceHeight, & ! (in) Height of interfaces (m) - OBLdepth_0d, & ! (out) OBL depth (m) - kOBL, & ! (out) level (+fraction) of OBL extent - zt_cntr=cellHeight, & ! (in) Height of cell centers (m) - surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) - surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) - Coriolis=Coriolis, & ! (in) Coriolis parameter (1/s) - CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters - - if (CS%deepOBLoffset>0.) then - zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1)) - OBLdepth_0d = min( OBLdepth_0d, -zBottomMinusOffset ) - kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) - endif - - ! apply some constraints on OBLdepth - if(CS%fixedOBLdepth) OBLdepth_0d = CS%fixedOBLdepth_value - OBLdepth_0d = max( OBLdepth_0d, -iFaceHeight(2) ) ! no shallower than top layer - OBLdepth_0d = min( OBLdepth_0d, -iFaceHeight(G%ke+1) ) ! no deep than bottom - kOBL = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, OBLdepth_0d ) - - endif ! endif for "correction" step - -! smg: remove code above -! ********************************************************************** - ! Call CVMix/KPP to obtain OBL diffusivities, viscosities and non-local transports @@ -1076,7 +980,7 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & if (CS%SW_METHOD .eq. SW_METHOD_ALL_SW) then surfBuoyFlux = buoyFlux(i,j,1) elseif (CS%SW_METHOD .eq. SW_METHOD_MXL_SW) then - surfBuoyFlux = buoyFlux(i,j,1) - buoyFlux(i,j,int(kOBL)+1) ! We know the actual buoyancy flux into the OBL + surfBuoyFlux = buoyFlux(i,j,1) - buoyFlux(i,j,int(CS%kOBL(i,j))+1) ! We know the actual buoyancy flux into the OBL elseif (CS%SW_METHOD .eq. SW_METHOD_LV1_SW) then surfBuoyFlux = buoyFlux(i,j,1) - buoyFlux(i,j,2) endif @@ -1099,8 +1003,8 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & Kviscosity, & ! (in) Original viscosity (m2/s) Kdiffusivity(:,1), & ! (in) Original heat diffusivity (m2/s) Kdiffusivity(:,2), & ! (in) Original salt diffusivity (m2/s) - OBLdepth_0d, & ! (in) OBL depth (m) - kOBL, & ! (in) level (+fraction) of OBL extent + CS%OBLdepth(i,j), & ! (in) OBL depth (m) + CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent nonLocalTrans(:,1),& ! (out) Non-local heat transport (non-dimensional) nonLocalTrans(:,2),& ! (out) Non-local salt transport (non-dimensional) surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) @@ -1122,26 +1026,26 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & if (surfBuoyFlux < 0.0) then if (CS%NLT_shape == NLT_SHAPE_CUBIC) then do k = 2, G%ke - sigma = min(1.0,-iFaceHeight(k)/OBLdepth_0d) + sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) nonLocalTrans(k,1) = (1.0 - sigma)**2 * (1.0 + 2.0*sigma) !* nonLocalTrans(k,2) = nonLocalTrans(k,1) enddo elseif (CS%NLT_shape == NLT_SHAPE_PARABOLIC) then do k = 2, G%ke - sigma = min(1.0,-iFaceHeight(k)/OBLdepth_0d) + sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) nonLocalTrans(k,1) = (1.0 - sigma)**2 !*CS%CS2 nonLocalTrans(k,2) = nonLocalTrans(k,1) enddo elseif (CS%NLT_shape == NLT_SHAPE_LINEAR) then do k = 2, G%ke - sigma = min(1.0,-iFaceHeight(k)/OBLdepth_0d) + sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) nonLocalTrans(k,1) = (1.0 - sigma)!*CS%CS2 nonLocalTrans(k,2) = nonLocalTrans(k,1) enddo elseif (CS%NLT_shape == NLT_SHAPE_CUBIC_LMD) then ! Sanity check (should agree with CVMix result using simple matching) do k = 2, G%ke - sigma = min(1.0,-iFaceHeight(k)/OBLdepth_0d) + sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j)) nonLocalTrans(k,1) = CS%CS2 * sigma*(1.0 -sigma)**2 nonLocalTrans(k,2) = nonLocalTrans(k,1) enddo @@ -1163,8 +1067,8 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & ! recompute wscale for diagnostics, now that we in fact know boundary layer depth if (CS%id_Ws > 0) then call CVMix_kpp_compute_turbulent_scales( & - -CellHeight/OBLdepth_0d, & ! (in) Normalized boundary layer coordinate - OBLdepth_0d, & ! (in) OBL depth (m) + -CellHeight/CS%OBLdepth(i,j), & ! (in) Normalized boundary layer coordinate + CS%OBLdepth(i,j), & ! (in) OBL depth (m) surfBuoyFlux, & ! (in) Buoyancy flux at surface (m2/s3) surfFricVel, & ! (in) Turbulent friction velocity at surface (m/s) w_s=Ws_1d, & ! (out) Turbulent velocity scale profile (m/s) @@ -1184,13 +1088,12 @@ 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 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(:) if (CS%id_sigma > 0) then CS%sigma(i,j,:) = 0. - if (OBLdepth_0d>0.) CS%sigma(i,j,:) = -iFaceHeight/OBLdepth_0d + if (CS%OBLdepth(i,j)>0.) CS%sigma(i,j,:) = -iFaceHeight/CS%OBLdepth(i,j) endif if (CS%id_N > 0) CS%N(i,j,:) = N_1d(:) if (CS%id_N2 > 0) CS%N2(i,j,:) = N2_1d(:) From 29a458ac92590330618d161af08c8a25bd17ca5c Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 10 May 2018 14:18:48 -0600 Subject: [PATCH 16/18] remove unnecessary vars in KPP_calculate --- src/parameterizations/vertical/MOM_KPP.F90 | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index 3ae2a1a196..eb5d27d7bf 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -780,10 +780,9 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & real, dimension( 3*G%ke ) :: Temp_1D real, dimension( 3*G%ke ) :: Salt_1D - real :: surfFricVel, surfBuoyFlux, Coriolis + real :: surfFricVel, surfBuoyFlux real :: GoRho, pRef, rho1, rhoK, rhoKm1, Uk, Vk, sigma - real :: zBottomMinusOffset ! Height of bottom plus a little bit (m) real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth. real :: hTot ! Running sum of thickness used in the surface layer average (m) real :: delH ! Thickness of a layer (m) @@ -819,13 +818,13 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & !$OMP buoyFlux, nonLocalTransHeat, & !$OMP nonLocalTransScalar,Kt,Ks,Kv) & !$OMP firstprivate(nonLocalTrans) & -!$OMP private(Coriolis,surfFricVel,SLdepth_0d,hTot,surfTemp, & +!$OMP private(surfFricVel,SLdepth_0d,hTot,surfTemp, & !$OMP surfHtemp,surfSalt,surfHsalt,surfU, & !$OMP surfHu,surfV,surfHv,iFaceHeight, & !$OMP pRef,km1,cellHeight,Uk,Vk,deltaU2, & !$OMP rho1,rhoK,rhoKm1,deltaRho,N2_1d,N_1d,delH, & !$OMP surfBuoyFlux,Ws_1d,Vt2_1d,BulkRi_1d, & -!$OMP zBottomMinusOffset,Kdiffusivity, & +!$OMP Kdiffusivity, & !$OMP Kviscosity,sigma,kk,pres_1D,Temp_1D, & !$OMP Salt_1D,rho_1D,surfBuoyFlux2,ksfc,dh,hcorr) @@ -837,8 +836,6 @@ subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & if (G%mask2dT(i,j)==0.) cycle ! things independent of position within the column - Coriolis = 0.25*( (G%CoriolisBu(i,j) + G%CoriolisBu(i-1,j-1)) & - +(G%CoriolisBu(i-1,j) + G%CoriolisBu(i,j-1)) ) surfFricVel = uStar(i,j) ! Bullk Richardson number computed for each cell in a column, From 8f730569e1d091ec3219bb73db95da545434b1fc Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 10 May 2018 15:52:01 -0600 Subject: [PATCH 17/18] implement smoothing on OBL depth --- src/parameterizations/vertical/MOM_KPP.F90 | 68 ++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/src/parameterizations/vertical/MOM_KPP.F90 b/src/parameterizations/vertical/MOM_KPP.F90 index eb5d27d7bf..e479460ebe 100644 --- a/src/parameterizations/vertical/MOM_KPP.F90 +++ b/src/parameterizations/vertical/MOM_KPP.F90 @@ -728,9 +728,77 @@ subroutine KPP_compute_BLD(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, buoyFlux) enddo enddo + if (CS%smoothBLD) call KPP_smooth_BLD(CS,G,GV,h) + end subroutine KPP_compute_BLD +!> Apply a 1-1-4-1-1 Laplacian filter one time on BLD to reduce any horizontal two-grid-point noise +subroutine KPP_smooth_BLD(CS,G,GV,h) + ! Arguments + type(KPP_CS), pointer :: CS !< Control structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid + type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid + real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses (units of H) + + ! local + real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface (m) (negative in ocean) + real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface (m) (negative in ocean) + integer :: i, j, k + real :: wc, ww, we, wn, ws ! averaging weights for smoothing + real :: dh ! The local thickness used for calculating interface positions (m) + real :: hcorr ! A cumulative correction arising from inflation of vanished layers (m) + + ! apply smoothing on OBL depth + do j = G%jsc, G%jec + do i = G%isc, G%iec + + ! skip land points + if (G%mask2dT(i,j)==0.) cycle + + ! compute weights + ww = 0.125 * G%mask2dT(i-1,j) + we = 0.125 * G%mask2dT(i+1,j) + ws = 0.125 * G%mask2dT(i,j-1) + wn = 0.125 * G%mask2dT(i,j+1) + wc = 1.0 - (ww+we+wn+ws) + + CS%OBLdepth(i,j) = wc * CS%OBLdepth(i,j) & + + ww * CS%OBLdepth(i-1,j) & + + we * CS%OBLdepth(i+1,j) & + + ws * CS%OBLdepth(i,j-1) & + + wn * CS%OBLdepth(i,j+1) + enddo + enddo + + ! Update kOBL for smoothed OBL depths + do j = G%jsc, G%jec + do i = G%isc, G%iec + + ! skip land points + if (G%mask2dT(i,j)==0.) cycle + + iFaceHeight(1) = 0.0 ! BBL is all relative to the surface + hcorr = 0. + do k=1,G%ke + + ! cell center and cell bottom in meters (negative values in the ocean) + 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 + + CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) ) + + enddo + enddo + +end subroutine KPP_smooth_BLD + + !> KPP vertical diffusivity/viscosity and non-local tracer transport subroutine KPP_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, & buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,& From bde57f32885d5286f4c33a7b500de4cc043db46e Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Tue, 15 May 2018 09:54:56 -0600 Subject: [PATCH 18/18] rm trailing whitespaces --- 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 9321d3f68c..cb868d9e95 100644 --- a/src/parameterizations/vertical/MOM_tidal_mixing.F90 +++ b/src/parameterizations/vertical/MOM_tidal_mixing.F90 @@ -802,7 +802,7 @@ subroutine calculate_CVMix_tidal(h, j, G, GV, CS, N2_int, Kd) exp_hab_zetar = exp_hab_zetar, & zw = iFaceHeight, & CVmix_tidal_params_user = CS%CVMix_tidal_params) - !TODO: in above call, there is no need to pass efficiency, since it gets + !TODO: in above call, there is no need to pass efficiency, since it gets ! passed via CVMix_init_tidal and stored in CVMix_tidal_params. Change ! CVMix API to prevent this redundancy.