diff --git a/src/core/MOM_continuity_PPM.F90 b/src/core/MOM_continuity_PPM.F90 index 49a1a6e674..de0285c2a0 100644 --- a/src/core/MOM_continuity_PPM.F90 +++ b/src/core/MOM_continuity_PPM.F90 @@ -1017,10 +1017,12 @@ subroutine set_zonal_BT_cont(u, h_in, hL, hR, BT_cont, uh_tot_0, duhdu_tot_0, & do k=1,nz ; do I=ish-1,ieh ; if (do_i(I)) then visc_rem_lim = max(visc_rem(I,k), min_visc_rem*visc_rem_max(I)) - if (u(I,j,k) + duR(I)*visc_rem_lim > -du_CFL(I)*visc_rem(I,k)) & - duR(I) = -(u(I,j,k) + du_CFL(I)*visc_rem(I,k)) / visc_rem_lim - if (u(I,j,k) + duL(I)*visc_rem_lim < du_CFL(I)*visc_rem(I,k)) & - duL(I) = -(u(I,j,k) - du_CFL(I)*visc_rem(I,k)) / visc_rem_lim + if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points. + if (u(I,j,k) + duR(I)*visc_rem_lim > -du_CFL(I)*visc_rem(I,k)) & + duR(I) = -(u(I,j,k) + du_CFL(I)*visc_rem(I,k)) / visc_rem_lim + if (u(I,j,k) + duL(I)*visc_rem_lim < du_CFL(I)*visc_rem(I,k)) & + duL(I) = -(u(I,j,k) - du_CFL(I)*visc_rem(I,k)) / visc_rem_lim + endif endif ; enddo ; enddo do k=1,nz @@ -1770,10 +1772,12 @@ subroutine set_merid_BT_cont(v, h_in, hL, hR, BT_cont, vh_tot_0, dvhdv_tot_0, & do k=1,nz ; do i=ish,ieh ; if (do_i(i)) then visc_rem_lim = max(visc_rem(i,k), min_visc_rem*visc_rem_max(i)) - if (v(i,J,k) + dvR(i)*visc_rem_lim > -dv_CFL(i)*visc_rem(i,k)) & - dvR(i) = -(v(i,J,k) + dv_CFL(i)*visc_rem(i,k)) / visc_rem_lim - if (v(i,J,k) + dvL(i)*visc_rem_lim < dv_CFL(i)*visc_rem(i,k)) & - dvL(i) = -(v(i,J,k) - dv_CFL(i)*visc_rem(i,k)) / visc_rem_lim + if (visc_rem_lim > 0.0) then ! This is almost always true for ocean points. + if (v(i,J,k) + dvR(i)*visc_rem_lim > -dv_CFL(i)*visc_rem(i,k)) & + dvR(i) = -(v(i,J,k) + dv_CFL(i)*visc_rem(i,k)) / visc_rem_lim + if (v(i,J,k) + dvL(i)*visc_rem_lim < dv_CFL(i)*visc_rem(i,k)) & + dvL(i) = -(v(i,J,k) - dv_CFL(i)*visc_rem(i,k)) / visc_rem_lim + endif endif ; enddo ; enddo do k=1,nz do i=ish,ieh ; if (do_i(i)) then diff --git a/src/diagnostics/MOM_diag_to_Z.F90 b/src/diagnostics/MOM_diag_to_Z.F90 index d19fd82dba..8eb34d9604 100644 --- a/src/diagnostics/MOM_diag_to_Z.F90 +++ b/src/diagnostics/MOM_diag_to_Z.F90 @@ -850,6 +850,8 @@ subroutine register_Z_tracer(tr_ptr, name, long_name, units, Time, G, CS, standa character(len=256) :: posted_cmor_standard_name character(len=256) :: posted_cmor_long_name + if (CS%nk_zspace<1) return + if (present(standard_name)) then posted_standard_name = standard_name else @@ -858,7 +860,7 @@ subroutine register_Z_tracer(tr_ptr, name, long_name, units, Time, G, CS, standa call register_Z_tracer_low(tr_ptr, name, long_name, units, trim(posted_standard_name), Time, G, CS) - if (present(cmor_field_name)) then + if (present(cmor_field_name)) then ! Fallback values for strings set to "NULL" posted_cmor_units = "not provided" ! posted_cmor_standard_name = "not provided" ! values might be replaced with a CS%missing field? @@ -966,9 +968,10 @@ subroutine MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS) character(len=40) :: mod = "MOM_diag_to_Z" ! module name character(len=200) :: in_dir, zgrid_file ! strings for directory/file - character(len=48) :: flux_units + character(len=48) :: flux_units, string integer :: z_axis, zint_axis integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nk, id_test + logical :: diag_mediator_is_using_z isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nk = G%ke IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB @@ -991,7 +994,19 @@ subroutine MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS) "The file that specifies the vertical grid for \n"//& "depth-space diagnostics, or blank to disable \n"//& "depth-space output.", default="") + + ! Check that the diag_mediator z-sapce remapping is not using the same module name + string = '' + call get_param(param_file, mod, "DIAG_REMAP_Z_MODULE_SUFFIX", string, & + default='_z_new', do_not_log=.true.) + diag_mediator_is_using_z = .false. + if (trim(string) == '_z') diag_mediator_is_using_z = .true. + if (len_trim(zgrid_file) > 0) then + if (diag_mediator_is_using_z) call MOM_error(FATAL, "MOM_diag_to_Z_init:"// & + "Z_OUTPUT_GRID_FILE can not be used when DIAG_REMAP_Z_MODULE_SUFFIX='_z'." // & + " Z_OUTPUT_GRID_FILE='"//trim(zgrid_file)//"'") + call get_param(param_file, mod, "INPUTDIR", in_dir, & "The directory in which input files are found.", default=".") in_dir = slasher(in_dir) @@ -1004,7 +1019,7 @@ subroutine MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS) "from the size of the variable zw in the output grid file.", & units="nondim") else - in_dir = "" ; CS%nk_zspace = -1 + CS%nk_zspace = -1 endif if (CS%nk_zspace > 0) then @@ -1041,7 +1056,7 @@ subroutine MOM_diag_to_Z_init(Time, G, GV, param_file, diag, CS) missing_value=CS%missing_trans) if (CS%id_vh_z>0) call safe_alloc_ptr(CS%vh_z,isd,ied,JsdB,JedB,CS%nk_zspace) - else + elseif (.not. diag_mediator_is_using_z) then ! Check whether diag-table is requesting any z-space files; issue a warning if it is. @@ -1229,6 +1244,7 @@ function ocean_register_diag_with_z(tr_ptr, vardesc_tr, G, Time, CS) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nk = G%ke if (.not.associated(CS)) call MOM_error(FATAL, & "register_Z_tracer: Module must be initialized before it is used.") + if (CS%nk_zspace<1) return if (CS%num_tr_used >= MAX_FIELDS_) then call MOM_error(WARNING,"ocean_register_diag_with_z: Attempted to register and use "//& diff --git a/src/framework/MOM_diag_mediator.F90 b/src/framework/MOM_diag_mediator.F90 index ebfd8e2569..be8d5cc104 100644 --- a/src/framework/MOM_diag_mediator.F90 +++ b/src/framework/MOM_diag_mediator.F90 @@ -28,6 +28,8 @@ module MOM_diag_mediator !********+*********+*********+*********+*********+*********+*********+** use MOM_coms, only : PE_here +use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end +use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE use MOM_error_handler, only : MOM_error, FATAL, is_root_pe use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type @@ -162,13 +164,16 @@ module MOM_diag_mediator ! Number of z levels used for remapping integer :: nz_remap - ! Define z star on u, v, T grids, these are the interface positions + ! Output grid thicknesses real, dimension(:,:,:), allocatable :: h_zoutput ! Keep track of which remapping is needed for diagnostic output logical :: do_z_remapping_on_u, do_z_remapping_on_v, do_z_remapping_on_T logical :: remapping_initialized + !> String appended to module name for z*-remapped diagnostics + character(len=6) :: z_remap_suffix = '_z_new' + ! Pointer to H and G for remapping real, dimension(:,:,:), pointer :: h => null() type(ocean_grid_type), pointer :: G => null() @@ -181,6 +186,9 @@ module MOM_diag_mediator end type diag_ctrl +! CPU clocks +integer :: id_clock_diag_mediator, id_clock_diag_z_remap, id_clock_diag_grid_updates + contains !> Sets up diagnostics axes @@ -312,6 +320,17 @@ subroutine set_axes_info(G, GV, param_file, diag_cs, set_vertical) "Depth of cell center", direction=-1) id_zzi = diag_axis_init('zi_remap', diag_cs%zi_remap, "meters", "z", & 'Depth of interfaces', direction=-1) + call get_param(param_file, mod, "DIAG_REMAP_Z_MODULE_SUFFIX", diag_cs%z_remap_suffix, & + 'This is the string attached to the end of "ocean_model"\n'// & + 'for use in the model column of the diag_table to indicate\n'// & + 'a diagnostic should be remapped to z*-coordinates.', & + default='_z_new') + if (trim(diag_cs%z_remap_suffix) == '_z') then + ! This will conflict with the older MOM_diag_to_Z module for z-output + call get_param(param_file, mod, "Z_OUTPUT_GRID_FILE", string, default="", do_not_log=.true.) + if (len(trim(string))>0) call MOM_error(FATAL,"MOM_diag_mediator, set_axes_info: "// & + "Z_OUTPUT_GRID_FILE must be blank to use DIAG_REMAP_Z_MODULE_SUFFIX='_z'") + endif else ! In this case the axes associated with these will never be used, however ! they need to be positive otherwise FMS complains. @@ -471,6 +490,7 @@ subroutine post_data_0d(diag_field_id, field, diag_cs, is_static, mask) logical :: used, is_stat type(diag_type), pointer :: diag => null() + if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator) is_stat = .false. ; if (present(is_static)) is_stat = is_static ! Iterate over list of diag 'variants', e.g. CMOR aliases, call send_data @@ -487,6 +507,7 @@ subroutine post_data_0d(diag_field_id, field, diag_cs, is_static, mask) diag => diag%next enddo + if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator) end subroutine post_data_0d subroutine post_data_1d_k(diag_field_id, field, diag_cs, is_static) @@ -507,6 +528,7 @@ subroutine post_data_1d_k(diag_field_id, field, diag_cs, is_static) integer :: isv, iev, jsv, jev type(diag_type), pointer :: diag => null() + if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator) is_stat = .false. ; if (present(is_static)) is_stat = is_static ! Iterate over list of diag 'variants', e.g. CMOR aliases. @@ -522,6 +544,7 @@ subroutine post_data_1d_k(diag_field_id, field, diag_cs, is_static) diag => diag%next enddo + if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator) end subroutine post_data_1d_k subroutine post_data_2d(diag_field_id, field, diag_cs, is_static, mask) @@ -541,6 +564,7 @@ subroutine post_data_2d(diag_field_id, field, diag_cs, is_static, mask) type(diag_type), pointer :: diag => null() + if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator) ! Iterate over list of diag 'variants' (e.g. CMOR aliases) and post each. call assert(diag_field_id < diag_cs%next_free_diag_id, & 'post_data_2d: Unregistered diagnostic id') @@ -550,6 +574,7 @@ subroutine post_data_2d(diag_field_id, field, diag_cs, is_static, mask) diag => diag%next enddo + if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator) end subroutine post_data_2d subroutine post_data_2d_low(diag, field, diag_cs, is_static, mask) @@ -657,6 +682,7 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) type(diag_type), pointer :: diag => null() real, allocatable :: remapped_field(:,:,:) + if (id_clock_diag_mediator>0) call cpu_clock_begin(id_clock_diag_mediator) ! Iterate over list of diag 'variants', e.g. CMOR aliases, different vertical ! grids, and post each. call assert(diag_field_id < diag_cs%next_free_diag_id, & @@ -671,8 +697,10 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) call MOM_error(FATAL,"post_data_3d: no mask for regridded field.") endif + if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) allocate(remapped_field(DIM_I(field),DIM_J(field), diag_cs%nz_remap)) call remap_diag_to_z(field, diag, diag_cs, remapped_field) + if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) if (associated(diag%mask3d)) then ! Since 3d masks do not vary in the vertical, just use as much as is ! needed. @@ -681,12 +709,15 @@ subroutine post_data_3d(diag_field_id, field, diag_cs, is_static, mask) else call post_data_3d_low(diag, remapped_field, diag_cs, is_static) endif + if (id_clock_diag_z_remap>0) call cpu_clock_begin(id_clock_diag_z_remap) deallocate(remapped_field) + if (id_clock_diag_z_remap>0) call cpu_clock_end(id_clock_diag_z_remap) else call post_data_3d_low(diag, field, diag_cs, is_static, mask) endif diag => diag%next enddo + if (id_clock_diag_mediator>0) call cpu_clock_end(id_clock_diag_mediator) end subroutine post_data_3d @@ -817,6 +848,7 @@ subroutine diag_update_target_grids(diag_cs) if (.not. allocated(diag_cs%zi_remap)) then return endif + if (id_clock_diag_grid_updates>0) call cpu_clock_begin(id_clock_diag_grid_updates) if (.not. diag_cs%remapping_initialized) then call assert(allocated(diag_cs%zi_remap), & @@ -854,6 +886,7 @@ subroutine diag_update_target_grids(diag_cs) ! when doing remapping. diag_cs%h_old(:,:,:) = diag_cs%h(:,:,:) #endif + if (id_clock_diag_grid_updates>0) call cpu_clock_end(id_clock_diag_grid_updates) end subroutine diag_update_target_grids @@ -1127,7 +1160,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & ! Remap to z vertical coordinate, note that only diagnostics on layers ! (not interfaces) are supported, also B axes are not supported yet if (is_layer_axes(axes, diag_cs) .and. (.not. is_B_axes(axes, diag_cs)) .and. axes%rank == 3) then - if (get_diag_field_id_fms(trim(module_name)//'_z_new', field_name) /= DIAG_FIELD_NOT_FOUND) then + if (get_diag_field_id_fms(trim(module_name)//trim(diag_cs%z_remap_suffix), field_name) /= DIAG_FIELD_NOT_FOUND) then if (.not. allocated(diag_cs%zi_remap)) then call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// & 'destination grid spec provided, see param DIAG_REMAP_Z_GRID_DEF') @@ -1139,7 +1172,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & call set_diag_mask(z_remap_diag, diag_cs, axes) call set_diag_remap_axes(z_remap_diag, diag_cs, axes) call assert(associated(z_remap_diag%remap_axes), 'register_diag_field: remap axes not set') - fms_id = register_diag_field_fms(module_name//'_z_new', field_name, & + fms_id = register_diag_field_fms(module_name//trim(diag_cs%z_remap_suffix), field_name, & z_remap_diag%remap_axes%handles, & init_time, long_name=long_name, units=units, missing_value=MOM_missing_value, & range=range, mask_variant=mask_variant, standard_name=standard_name, & @@ -1160,13 +1193,13 @@ function register_diag_field(module_name, field_name, axes, init_time, & if (is_root_pe() .and. diag_CS%doc_unit > 0) then msg = '' if (present(cmor_field_name)) msg = 'CMOR equivalent is "'//trim(cmor_field_name)//'"' - call log_available_diag(associated(z_remap_diag), module_name//'_z_new', field_name, & + call log_available_diag(associated(z_remap_diag), module_name//trim(diag_cs%z_remap_suffix), field_name, & cm_string, msg, diag_CS, long_name, units, standard_name) endif ! Remap to z vertical coordinate with CMOR names and attributes if (present(cmor_field_name)) then - if (get_diag_field_id_fms(module_name//'_z_new', cmor_field_name) /= DIAG_FIELD_NOT_FOUND) then + if (get_diag_field_id_fms(module_name//trim(diag_cs%z_remap_suffix), cmor_field_name) /= DIAG_FIELD_NOT_FOUND) then if (.not. allocated(diag_cs%zi_remap)) then call MOM_error(FATAL, 'register_diag_field: Request to regrid but no '// & 'destination grid spec provided, see param DIAG_REMAP_Z_GRID_DEF') @@ -1178,7 +1211,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & call set_diag_mask(cmor_z_remap_diag, diag_cs, axes) call set_diag_remap_axes(cmor_z_remap_diag, diag_cs, axes) call assert(associated(cmor_z_remap_diag%remap_axes), 'register_diag_field: remap axes not set') - fms_id = register_diag_field_fms(module_name//'_z_new', cmor_field_name, & + fms_id = register_diag_field_fms(module_name//trim(diag_cs%z_remap_suffix), cmor_field_name, & cmor_z_remap_diag%remap_axes%handles, & init_time, long_name=trim(posted_cmor_long_name), units=trim(posted_cmor_units), missing_value=MOM_missing_value, & range=range, mask_variant=mask_variant, standard_name=trim(posted_cmor_standard_name), & @@ -1198,7 +1231,7 @@ function register_diag_field(module_name, field_name, axes, init_time, & endif if (is_root_pe() .and. diag_CS%doc_unit > 0) then msg = 'native name is "'//trim(field_name)//'"' - call log_available_diag(associated(cmor_z_remap_diag), module_name//'_z_new', cmor_field_name, & + call log_available_diag(associated(cmor_z_remap_diag), module_name//trim(diag_cs%z_remap_suffix), cmor_field_name, & cm_string, msg, diag_CS, posted_cmor_long_name, posted_cmor_units, & posted_cmor_standard_name) endif @@ -1629,6 +1662,10 @@ subroutine diag_mediator_init(G, param_file, diag_cs, err_msg, doc_file_dir) character(len=240) :: doc_file, doc_file_dflt, doc_path character(len=40) :: mod = "MOM_diag_mediator" ! This module's name. + id_clock_diag_mediator = cpu_clock_id('(Ocean diagnostics framework)', grain=CLOCK_MODULE) + id_clock_diag_z_remap = cpu_clock_id('(Ocean diagnostics remapping)', grain=CLOCK_ROUTINE) + id_clock_diag_grid_updates = cpu_clock_id('(Ocean diagnostics grid updates)', grain=CLOCK_ROUTINE) + call diag_manager_init(err_msg=err_msg) ! Allocate and initialise list of all diagnostics (and variants) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index d0963411fe..08825eaff5 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -854,7 +854,7 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h) ! of salinity and temperature within each layer. character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path real :: scale_factor, min_thickness - integer :: i, j + integer :: i, j, k call get_param(PF, mod, "SURFACE_PRESSURE_FILE", p_surf_file, & "The initial condition file for the surface height.", & @@ -883,7 +883,11 @@ subroutine trim_for_ice(PF, G, GV, ALE_CSp, tv, h) ! call pressure_gradient_ppm(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h) ! endif else - call MOM_error(FATAL, "trim_for_ice: Does not work without ALE mode") +! call MOM_error(FATAL, "trim_for_ice: Does not work without ALE mode") + do k=1,G%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec + T_t(i,j,k) = tv%T(i,j,k) ; T_b(i,j,k) = tv%T(i,j,k) + S_t(i,j,k) = tv%S(i,j,k) ; S_b(i,j,k) = tv%S(i,j,k) + enddo ; enddo ; enddo endif do j=G%jsc,G%jec ; do i=G%isc,G%iec @@ -1236,6 +1240,7 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, param_file, eqn_of_state, P_Ref real :: drho_dT(SZK_(G)) ! Derivative of density with temperature in kg m-3 K-1. ! real :: drho_dS(SZK_(G)) ! Derivative of density with salinity in kg m-3 PSU-1. ! real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 in kg m-3. + logical :: fit_salin ! If true, accept the prescribed temperature and fit the salinity. character(len=40) :: mod = "initialize_temp_salt_fit" ! This subroutine's name. integer :: i, j, k, itt, nz nz = G%ke @@ -1248,27 +1253,44 @@ subroutine initialize_temp_salt_fit(T, S, G, GV, param_file, eqn_of_state, P_Ref call get_param(param_file, mod, "S_REF", S_Ref, & "A reference salinity used in initialization.", units="PSU", & default=35.0) + call get_param(param_file, mod, "FIT_SALINITY", fit_salin, & + "If true, accept the prescribed temperature and fit the \n"//& + "salinity; otherwise take salinity and fit temperature.", & + default=.false.) do k=1,nz pres(k) = P_Ref ; S0(k) = S_Ref + T0(k) = T_Ref enddo - T0(1) = T_Ref call calculate_density(T0(1),S0(1),pres(1),rho_guess(1),eqn_of_state) call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,1,eqn_of_state) -! A first guess of the layers' temperatures. ! - do k=nz,1,-1 - T0(k) = T0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dT(1) - enddo - -! Refine the guesses for each layer. ! - do itt=1,6 - call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state) - call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state) - do k=1,nz - T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) + if (fit_salin) then +! A first guess of the layers' temperatures. + do k=nz,1,-1 + S0(k) = max(0.0, S0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dS(1)) enddo - enddo +! Refine the guesses for each layer. + do itt=1,6 + call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state) + call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state) + do k=1,nz + S0(k) = max(0.0, S0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dS(k)) + enddo + enddo + else +! A first guess of the layers' temperatures. + do k=nz,1,-1 + T0(k) = T0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dT(1) + enddo + do itt=1,6 + call calculate_density(T0,S0,pres,rho_guess,1,nz,eqn_of_state) + call calculate_density_derivs(T0,S0,pres,drho_dT,drho_dS,1,nz,eqn_of_state) + do k=1,nz + T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k) + enddo + enddo + endif do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied T(i,j,k) = T0(k) ; S(i,j,k) = S0(k) diff --git a/src/initialization/MOM_tracer_initialization_from_Z.F90 b/src/initialization/MOM_tracer_initialization_from_Z.F90 index ccd216e9ba..2dcbe8511f 100644 --- a/src/initialization/MOM_tracer_initialization_from_Z.F90 +++ b/src/initialization/MOM_tracer_initialization_from_Z.F90 @@ -23,7 +23,7 @@ module MOM_tracer_initialization_from_Z use MOM_verticalGrid, only : setVerticalGridAxes use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type use MOM_EOS, only : int_specific_vol_dp -use MOM_ALE, only : ALE_initRegridding, ALE_CS, ALE_initThicknessToCoord +use MOM_ALE, only : ALE_initRegridding, ALE_CS, ALE_initThicknessToCoord, ALE_remap_scalar use MOM_regridding, only : regridding_CS use MOM_remapping, only : remapping_CS, initialize_remapping use MOM_remapping, only : remapping_core_h @@ -106,10 +106,9 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, PF, src_file, src_var_nam, real, dimension(:), allocatable :: h1, h2, hTarget, deltaE, tmpT1d real, dimension(:), allocatable :: tmpT1dIn real :: zTopOfCell, zBottomOfCell - type(regridding_CS) :: regridCS ! Regridding parameters and work arrays type(remapping_CS) :: remapCS ! Remapping parameters and work arrays - real, dimension(:,:,:), allocatable :: tmp1 + real, dimension(:,:,:), allocatable :: hSrc real :: tempAvg, missing_value integer :: nPoints, ans @@ -167,6 +166,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, PF, src_file, src_var_nam, call cpu_clock_begin(id_clock_ALE) ! First we reserve a work space for reconstructions of the source data allocate( h1(kd) ) + allocate( hSrc(isd:ied,jsd:jed,kd) ) allocate( tmpT1dIn(kd) ) call initialize_remapping( remapCS, remapScheme, boundary_extrapolation=.false. ) ! Data for reconstructions ! Next we initialize the regridding package so that it knows about the target grid @@ -174,9 +174,7 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, PF, src_file, src_var_nam, allocate( h2(nz) ) allocate( tmpT1d(nz) ) allocate( deltaE(nz+1) ) - ! This call can be more general but is hard-coded for z* coordinates... ???? - call ALE_initRegridding( G, GV, PF, mod, regridCS, hTarget ) ! sets regridCS and hTarget(1:nz) - ! For each column ... + do j = js, je ; do i = is, ie if (G%mask2dT(i,j)>0.) then ! Build the source grid @@ -196,22 +194,15 @@ subroutine MOM_initialize_tracer_from_Z(h, tr, G, GV, PF, src_file, src_var_nam, zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k enddo h1(kd) = h1(kd) + ( zTopOfCell + G%bathyT(i,j) ) ! In case data is deeper than model - ! Build the target grid combining hTarget and topography - zTopOfCell = 0. ; zBottomOfCell = 0. - do k = 1, nz - zBottomOfCell = max( zTopOfCell - hTarget(k), -G%bathyT(i,j) ) - h2(k) = zTopOfCell - zBottomOfCell - zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k - enddo - ! Now remap from h1 to h2 - call remapping_core_h( nPoints, h1, tmpT1dIn, nz, h2, tmpT1d, remapCS ) ! sets tmpT1d -!!!MJH h(i,j,:) = h2(:) - tr(i,j,:) = tmpT1d(:) else tr(i,j,:) = 0. -!!!MJH h(i,j,:) = 0. endif ! mask2dT + hSrc(i,j,:) = h1(:) enddo ; enddo + + call ALE_remap_scalar(remapCS, G, kd, hSrc, tr_z, h, tr, all_cells=.true. ) + + deallocate( hSrc ) deallocate( h1 ) deallocate( h2 ) deallocate( hTarget ) diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 26ae3bdfb2..0231035d0d 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -267,55 +267,52 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, fluxes, dt, MLD, G, GV, ! U - Component !$OMP do - do j=js,je - do i=is-1,ie ; utimescale_diag(i,j) = 0.0 ; enddo - do i=is-1,ie ; vtimescale_diag(i,j) = 0.0 ; enddo - do I=is-1,ie - h_vel = 0.5*((htot(i,j) + htot(i+1,j)) + h_neglect) * GV%H_to_m + do j=js,je ; do I=is-1,ie + h_vel = 0.5*((htot(i,j) + htot(i+1,j)) + h_neglect) * GV%H_to_m - u_star = 0.5*(fluxes%ustar(i,j) + fluxes%ustar(i+1,j)) - absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) - ! peak ML visc: u_star * 0.41 * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) - ! momentum mixing rate: pi^2*visc/h_ml^2 - ! 0.41 is the von Karmen constant, 9.8696 = pi^2. - mom_mixrate = (0.41*9.8696)*u_star**2 / & - (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star) - timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) + u_star = 0.5*(fluxes%ustar(i,j) + fluxes%ustar(i+1,j)) + absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) + ! peak ML visc: u_star * 0.41 * (h_ml*u_star)/(absf*h_ml + 4.0*u_star) + ! momentum mixing rate: pi^2*visc/h_ml^2 + ! 0.41 is the von Karmen constant, 9.8696 = pi^2. + mom_mixrate = (0.41*9.8696)*u_star**2 / & + (absf*h_vel**2 + 4.0*(h_vel+dz_neglect)*u_star) + timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) - timescale = timescale * CS%ml_restrat_coef -! timescale = timescale*(2?)*(L_def/L_MLI)*min(EKE/MKE,1.0 + G%dyCv(i,j)**2/L_def**2)) + timescale = timescale * CS%ml_restrat_coef +! timescale = timescale*(2?)*(L_def/L_MLI)*min(EKE/MKE,1.0 + G%dyCv(i,j)**2/L_def**2)) - utimescale_diag(I,j) = timescale - uDml(I) = timescale * G%mask2dCu(I,j)*G%dyCu(I,j)* & - G%IdxCu(I,j)*(Rml_av(i+1,j)-Rml_av(i,j)) * (h_vel**2 * GV%m_to_H) + uDml(I) = timescale * G%mask2dCu(I,j)*G%dyCu(I,j)* & + G%IdxCu(I,j)*(Rml_av(i+1,j)-Rml_av(i,j)) * (h_vel**2 * GV%m_to_H) - if (uDml(i) == 0) then - do k=1,nz ; uhml(I,j,k) = 0.0 ; enddo - else - IhTot = 2.0 / ((htot(i,j) + htot(i+1,j)) + h_neglect) - zIHbelowVel = 0.0 - ! a(k) relates the sublayer transport to uDml with a linear profile. - ! The sum of a(k) through the mixed layers must be 0. - do k=1,nz - hAtVel = 0.5*(h(i,j,k) + h(i+1,j,k)) - zIHaboveVel = zIHbelowVel ! z/H for upper interface - zIHbelowVel = zIHbelowVel - (hAtVel * IhTot) ! z/H for lower interface - a(k) = PSI( zIHaboveVel ) - PSI( zIHbelowVel ) - if (a(k)*uDml(I) > 0.0) then - if (a(k)*uDml(I) > h_avail(i,j,k)) uDml(I) = h_avail(i,j,k) / a(k) - elseif (a(k)*uDml(I) < 0.0) then - if (-a(k)*uDml(I) > h_avail(i+1,j,k)) uDml(I) = -h_avail(i+1,j,k) / a(k) - endif - enddo - do k=1,nz - uhml(I,j,k) = a(k)*uDml(I) - uhtr(I,j,k) = uhtr(I,j,k) + uhml(I,j,k)*dt - enddo - endif - enddo - uDml_diag(is:ie,j) = uDml(is:ie) - enddo + if (uDml(i) == 0) then + do k=1,nz ; uhml(I,j,k) = 0.0 ; enddo + else + IhTot = 2.0 / ((htot(i,j) + htot(i+1,j)) + h_neglect) + zIHbelowVel = 0.0 + ! a(k) relates the sublayer transport to uDml with a linear profile. + ! The sum of a(k) through the mixed layers must be 0. + do k=1,nz + hAtVel = 0.5*(h(i,j,k) + h(i+1,j,k)) + zIHaboveVel = zIHbelowVel ! z/H for upper interface + zIHbelowVel = zIHbelowVel - (hAtVel * IhTot) ! z/H for lower interface + a(k) = PSI( zIHaboveVel ) - PSI( zIHbelowVel ) + if (a(k)*uDml(I) > 0.0) then + if (a(k)*uDml(I) > h_avail(i,j,k)) uDml(I) = h_avail(i,j,k) / a(k) + elseif (a(k)*uDml(I) < 0.0) then + if (-a(k)*uDml(I) > h_avail(i+1,j,k)) uDml(I) = -h_avail(i+1,j,k) / a(k) + endif + enddo + do k=1,nz + uhml(I,j,k) = a(k)*uDml(I) + uhtr(I,j,k) = uhtr(I,j,k) + uhml(I,j,k)*dt + enddo + endif + + utimescale_diag(I,j) = timescale + uDml_diag(I,j) = uDml(I) + enddo ; enddo ! V- component !$OMP do @@ -332,9 +329,7 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, fluxes, dt, MLD, G, GV, timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2) timescale = timescale * CS%ml_restrat_coef -! timescale = timescale*(2?)*(L_def/L_MLI)*min(EKE/MKE,1.0 + G%dyCv(i,j)**2/L_def**2)) - - vtimescale_diag(i,J) = timescale +! timescale = timescale*(2?)*(L_def/L_MLI)*min(EKE/MKE,1.0 + G%dyCv(i,j)**2/L_def**2)) vDml(i) = timescale * G%mask2dCv(i,J)*G%dxCv(i,J)* & G%IdyCv(i,J)*(Rml_av(i,j+1)-Rml_av(i,j)) * (h_vel**2 * GV%m_to_H) @@ -361,9 +356,10 @@ subroutine mixedlayer_restrat_general(h, uhtr, vhtr, tv, fluxes, dt, MLD, G, GV, vhtr(i,J,k) = vhtr(i,J,k) + vhml(i,J,k)*dt enddo endif - enddo - vDml_diag(is:ie,j) = vDml(is:ie) - enddo + + vtimescale_diag(i,J) = timescale + vDml_diag(i,J) = vDml(i) + enddo ; enddo !$OMP do do j=js,je ; do k=1,nz ; do i=is,ie diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90 index 4bf5e2dc0a..91acfe91c5 100644 --- a/src/user/Rossby_front_2d_initialization.F90 +++ b/src/user/Rossby_front_2d_initialization.F90 @@ -139,8 +139,8 @@ end subroutine Rossby_front_initialize_temperature_salinity subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, param_file) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure - real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: u !< i-component of velocity [m/s] - real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: v !< j-component of velocity [m/s] + real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: u !< i-component of velocity [m/s] + real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: v !< j-component of velocity [m/s] real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Thickness [H] type(param_file_type), intent(in) :: param_file !< Parameter file handle @@ -163,7 +163,8 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, param_file) do j = G%jsc,G%jec ; do I = G%isc-1,G%iec+1 f = 0.5*( G%CoriolisBu(I,j) + G%CoriolisBu(I,j-1) ) - dUdT = ( G%g_Earth * dRho_dT ) / ( f * GV%Rho0 ) + dUdT = 0.0 ; if (abs(f) > 0.0) & + dUdT = ( G%g_Earth * dRho_dT ) / ( f * GV%Rho0 ) Dml = Hml( G, G%geoLatT(i,j) ) Ty = dTdy( G, T_range, G%geoLatT(i,j) ) zi = 0. @@ -174,7 +175,7 @@ subroutine Rossby_front_initialize_velocity(u, v, h, G, GV, param_file) zm = max( zc + Dml, 0. ) ! Height above bottom of mixed layer u(I,j,k) = dUdT * Ty * zm ! Thermal wind starting at base of ML enddo - end do ; end do + enddo ; enddo end subroutine Rossby_front_initialize_velocity