From 217816cd1d916fdeb4f52e8ddfeb191aacf3830f Mon Sep 17 00:00:00 2001 From: Gustavo Marques Date: Tue, 7 Jan 2025 16:29:59 -0700 Subject: [PATCH] Add option to read MLE Frontal-Length Scale via file Introduced the `MLE_FRONT_LENGTH_FROM_FILE` parameter to allow specifying a spatially and temporally varying MLE frontal-length scale via netCDF file. Also added a new diagnostic, `mle_fl`, to save the frontal-length scale used in the mixed-layer restratification parameterization. --- .../lateral/MOM_mixed_layer_restrat.F90 | 73 +++++++++++++++++-- 1 file changed, 68 insertions(+), 5 deletions(-) diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 4e733f4ab9..732ec32693 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -9,12 +9,15 @@ module MOM_mixed_layer_restrat use MOM_diag_mediator, only : diag_update_remap_grids use MOM_domains, only : pass_var, To_West, To_South, Omit_Corners use MOM_error_handler, only : MOM_error, FATAL, WARNING -use MOM_file_parser, only : get_param, log_version, param_file_type +use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_file_parser, only : openParameterBlock, closeParameterBlock use MOM_forcing_type, only : mech_forcing, find_ustar use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type +use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init +use MOM_interpolate, only : external_field use MOM_intrinsic_functions, only : cuberoot +use MOM_io, only : slasher use MOM_lateral_mixing_coeffs, only : VarMix_CS use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS use MOM_unit_scaling, only : unit_scale_type @@ -48,6 +51,7 @@ module MOM_mixed_layer_restrat !! upscaling of buoyancy gradients that is otherwise represented !! by the parameter FOX_KEMPER_ML_RESTRAT_COEF. If MLE_FRONT_LENGTH is !! non-zero, it is recommended to set FOX_KEMPER_ML_RESTRAT_COEF=1.0. + logical :: fl_from_file !< If true, read the MLE front-length scale from a netCDF file. logical :: MLE_use_PBL_MLD !< If true, use the MLD provided by the PBL parameterization. !! if false, MLE will calculate a MLD based on a density difference !! based on the parameter MLE_DENSITY_DIFF. @@ -94,6 +98,9 @@ module MOM_mixed_layer_restrat type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the !! timing of diagnostic output. + type(external_field) :: sbc_fl !< A handle used in time interpolation of + !! front-lenght scales read from a file. + type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock. logical :: use_Stanley_ML !< If true, use the Stanley parameterization of SGS T variance real :: ustar_min !< A minimum value of ustar in thickness units to avoid numerical !! problems [H T-1 ~> m s-1 or kg m-2 s-1] @@ -123,6 +130,7 @@ module MOM_mixed_layer_restrat integer :: id_ustar = -1 integer :: id_bflux = -1 integer :: id_lfbod = -1 + integer :: id_mle_fl = -1 !>@} end type mixedlayer_restrat_CS @@ -207,7 +215,8 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix, htot_fast, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2] Rml_av_fast, & ! Negative g_Rho0 times the average mixed layer density or G_Earth ! times the average specific volume [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2] - MLD_slow, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2] + MLD_slow, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2] + mle_fl_2d, & ! MLE frontal lengh-scale [H ~> m or kg m-2] htot_slow, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2] Rml_av_slow ! Negative g_Rho0 times the average mixed layer density or G_Earth ! times the average specific volume [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2] @@ -257,6 +266,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix, real :: I_LFront ! The inverse of the frontal length scale [L-1 ~> m-1] real :: vonKar_x_pi2 ! A scaling constant that is approximately the von Karman constant times ! pi squared [nondim] + character(len=128) :: mesg logical :: line_is_empty, keep_going, res_upscale integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz @@ -338,7 +348,21 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix, h_neglect = GV%H_subroundoff if (CS%front_length>0.) then res_upscale = .true. - I_LFront = 1. / CS%front_length + do j=js-1,je+1 ; do i=is-1,ie+1 + mle_fl_2d(i,j) = CS%front_length + enddo ; enddo + elseif (CS%front_length == 0. .and. CS%fl_from_file) then + res_upscale = .true. + call time_interp_external(CS%sbc_fl, CS%Time, mle_fl_2d, turns=G%HI%turns) + call pass_var(mle_fl_2d, G%domain, halo=1) + do j=js,je ; do i=is,ie + if ((G%mask2dT(i,j) > 0.0) .and. (mle_fl_2d(i,j) < 0.0)) then + write(mesg,'(" Time_interp negative MLE frontal-lenght scale of ",(1pe12.4)," at i,j = ",& + & 2(i3), "lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') & + mle_fl_2d(i,j), i, j, G%geoLonT(i,j), G%geoLatT(i,j) + call MOM_error(FATAL, "MOM_mixed_layer_restrat mixedlayer_restrat_OM4: "//trim(mesg)) + endif + enddo ; enddo else res_upscale = .false. endif @@ -458,6 +482,8 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix, u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j))) absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) + ! Compute I_LFront = 1 / (frontal length scale) [m-1] + I_LFront = 2. / ( mle_fl_2d(i,j) + mle_fl_2d(i+1,j) ) ! If needed, res_scaling_fac = min( ds, L_d ) / l_f if (res_upscale) res_scaling_fac = & ( sqrt( 0.5 * ( (G%dxCu(I,j)**2) + (G%dyCu(I,j)**2) ) ) * I_LFront ) & @@ -543,7 +569,8 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix, !$OMP do do J=js-1,je ; do i=is,ie u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1))) - + ! Compute I_LFront = 1 / (frontal length scale) [m-1] + I_LFront = 2. / ( mle_fl_2d(i,j) + mle_fl_2d(i,j+1) ) absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) ! If needed, res_scaling_fac = min( ds, L_d ) / l_f if (res_upscale) res_scaling_fac = & @@ -652,6 +679,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix, if (CS%id_Rml > 0) call post_data(CS%id_Rml, Rml_av_fast, CS%diag) if (CS%id_uDml > 0) call post_data(CS%id_uDml, uDml_diag, CS%diag) if (CS%id_vDml > 0) call post_data(CS%id_vDml, vDml_diag, CS%diag) + if (CS%id_mle_fl > 0) call post_data(CS%id_mle_fl, mle_fl_2d, CS%diag) if (CS%id_uml > 0) then do J=js,je ; do i=is-1,ie @@ -1572,7 +1600,7 @@ end function growth_time !> Initialize the mixed layer restratification module logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, restart_CS) - type(time_type), intent(in) :: Time !< Current model time + type(time_type), target, intent(in) :: Time !< Current model time type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type @@ -1589,6 +1617,12 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, ! temperature variance [nondim] integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags ! This include declares and sets the variable "version". + character(len=200) :: inputdir ! The directory where NetCDF input files + character(len=240) :: mle_fl_filename ! A file from which chl_a concentrations are to be read. + character(len=128) :: mle_fl_file ! Data containing MLE front-lenght scale. Used + ! when reading from file. + character(len=32) :: fl_varname ! Name of front-lenght scale variable in mle_fl_file. + # include "version_variable.h" integer :: i, j @@ -1604,6 +1638,7 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, if (.not. mixedlayer_restrat_init) return CS%initialized = .true. + CS%Time => Time ! Nonsense values to cause problems when these parameters are not used CS%MLE_MLD_decay_time = -9.e9*US%s_to_T @@ -1613,6 +1648,7 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, CS%MLE_MLD_stretch = -9.e9 CS%use_Stanley_ML = .false. CS%use_Bodner = .false. + CS%fl_from_file = .false. call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.) call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & @@ -1735,6 +1771,29 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, "by the parameter FOX_KEMPER_ML_RESTRAT_COEF. If MLE_FRONT_LENGTH is "//& "non-zero, it is recommended to set FOX_KEMPER_ML_RESTRAT_COEF=1.0.",& units="m", default=0.0, scale=US%m_to_L) + call get_param(param_file, mdl, "MLE_FRONT_LENGTH_FROM_FILE", CS%fl_from_file, & + "If true, the MLE front-lenght scale is read from a file.", default=.false.) + if (CS%fl_from_file) then + call time_interp_external_init() + + call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") + call get_param(param_file, mdl, "MLE_FL_FILE", mle_fl_file, & + "MLE_FL_FILE is the file containing MLE front-lenght scale. "//& + "It is used when MLE_FRONT_LENGTH_FROM_FILE is true.", fail_if_missing=.true.) + mle_fl_filename = trim(slasher(inputdir))//trim(mle_fl_file) + call log_param(param_file, mdl, "INPUTDIR/MLE_FL_FILE", mle_fl_filename) + call get_param(param_file, mdl, "FL_VARNAME", fl_varname, & + "Name of MLE front-lenght scale variable in MLE_FL_FILE.", default='mle_fl') + if (modulo(G%Domain%turns, 4) /= 0) then + CS%sbc_fl = init_external_field(mle_fl_filename, trim(fl_varname), MOM_domain=G%Domain%domain_in) + else + CS%sbc_fl = init_external_field(mle_fl_filename, trim(fl_varname), MOM_domain=G%Domain) + endif + endif + if (CS%fl_from_file .and. CS%front_length>0.0) call MOM_error(FATAL, "mixedlayer_restrat_init: "// & + "MLE_FRONT_LENGTH_FROM_FILE cannot be true when MLE_FRONT_LENGTH > 0.0. "// & + "If you want to use MLE_FRONT_LENGTH, set MLE_FRONT_LENGTH_FROM_FILE to false." // & + "If you want to use MLE_FRONT_LENGTH_FROM_FILE, set MLE_FRONT_LENGTH to 0.0.") call get_param(param_file, mdl, "MLE_USE_PBL_MLD", CS%MLE_use_PBL_MLD, & "If true, the MLE parameterization will use the mixed-layer "//& "depth provided by the active PBL parameterization. If false, "//& @@ -1838,6 +1897,10 @@ logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, CS%id_lfbod = register_diag_field('ocean_model', 'lf_bodner', diag%axesT1, Time, & 'Front length in Bodner mixed layer restratificiation parameterization', & 'm', conversion=US%L_to_m) + else + CS%id_mle_fl = register_diag_field('ocean_model', 'mle_fl', diag%axesT1, Time, & + 'Frontal length scale used in the mixed layer restratificiation parameterization', & + 'm', conversion=US%L_to_m) endif ! If MLD_filtered is being used, we need to update halo regions after a restart