diff --git a/src/ALE/MOM_ALE.F90 b/src/ALE/MOM_ALE.F90 index 137f6cee9b..c300c8bd2f 100644 --- a/src/ALE/MOM_ALE.F90 +++ b/src/ALE/MOM_ALE.F90 @@ -1463,17 +1463,23 @@ subroutine ALE_writeCoordinateFile( CS, GV, directory ) end subroutine ALE_writeCoordinateFile !> Set h to coordinate values for fixed coordinate systems -subroutine ALE_initThicknessToCoord( CS, G, GV, h ) +subroutine ALE_initThicknessToCoord( CS, G, GV, h, height_units ) type(ALE_CS), intent(inout) :: CS !< module control structure type(ocean_grid_type), intent(in) :: G !< module grid structure type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< layer thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< layer thickness in thickness units + !! [H ~> m or kg m-2] or height units [Z ~> m] + logical, optional, intent(in) :: height_units !< If present and true, the + !! thicknesses are in height units ! Local variables + real :: scale ! A scaling value for the thicknesses [nondim] or [H Z-1 ~> nondim or kg m-3] integer :: i, j + scale = GV%Z_to_H + if (present(height_units)) then ; if (height_units) scale = 1.0 ; endif do j = G%jsd,G%jed ; do i = G%isd,G%ied - h(i,j,:) = GV%Z_to_H * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j)+G%Z_ref ) + h(i,j,:) = scale * getStaticThickness( CS%regridCS, 0., G%bathyT(i,j)+G%Z_ref ) enddo ; enddo end subroutine ALE_initThicknessToCoord diff --git a/src/diagnostics/MOM_obsolete_params.F90 b/src/diagnostics/MOM_obsolete_params.F90 index 7564137de8..21a09dfdbb 100644 --- a/src/diagnostics/MOM_obsolete_params.F90 +++ b/src/diagnostics/MOM_obsolete_params.F90 @@ -56,6 +56,7 @@ subroutine find_obsolete_params(param_file) hint="Instead use OBC_SEGMENT_xxx_VELOCITY_NUDGING_TIMESCALES.") enddo + call obsolete_logical(param_file, "CONVERT_THICKNESS_UNITS", .true.) call obsolete_logical(param_file, "MASK_MASSLESS_TRACERS", .false.) call obsolete_logical(param_file, "SALT_REJECT_BELOW_ML", .false.) diff --git a/src/initialization/MOM_state_initialization.F90 b/src/initialization/MOM_state_initialization.F90 index bd0931c694..dafdb553f3 100644 --- a/src/initialization/MOM_state_initialization.F90 +++ b/src/initialization/MOM_state_initialization.F90 @@ -150,7 +150,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & optional, intent(in) :: mass_shelf !< The mass per unit area of the overlying !! ice shelf [ R Z ~> kg m-2 ] ! Local variables - real :: depth_tot(SZI_(G),SZJ_(G)) ! The nominal total depth of the ocean [Z ~> m] + real :: depth_tot(SZI_(G),SZJ_(G)) ! The nominal total depth of the ocean [Z ~> m] + real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! The layer thicknesses in geopotential (z) units [Z ~> m] character(len=200) :: inputdir ! The directory where NetCDF input files are. character(len=200) :: config real :: H_rescale ! A rescaling factor for thicknesses from the representation in @@ -226,6 +227,9 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & !do k=1,nz ; do j=js,je ; do i=is,ie ! h(i,j,k) = 0. !enddo + + ! Initialize the layer thicknesses. + dz(:,:,:) = 0.0 endif ! Set the nominal depth of the ocean, which might be different from the bathymetric @@ -250,6 +254,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & "salinities from a Z-space file on a latitude-longitude grid.", & default=.false., do_not_log=just_read) + convert = new_sim ! Thicknesses are initialized in height units in most cases. if (from_Z_file) then ! Initialize thickness and T/S from z-coordinate data in a file. if (.NOT.use_temperature) call MOM_error(FATAL,"MOM_initialize_state : "//& @@ -257,14 +262,18 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & call MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, & just_read=just_read, frac_shelf_h=frac_shelf_h) + convert = .false. else ! Initialize thickness, h. call get_param(PF, mdl, "THICKNESS_CONFIG", config, & "A string that determines how the initial layer "//& "thicknesses are specified for a new run: \n"//& " \t file - read interface heights from the file specified \n"//& + " \t\t by (THICKNESS_FILE).\n"//& " \t thickness_file - read thicknesses from the file specified \n"//& " \t\t by (THICKNESS_FILE).\n"//& + " \t mass_file - read thicknesses in units of mass per unit area from the file \n"//& + " \t\t specified by (THICKNESS_FILE).\n"//& " \t coord - determined by ALE coordinate.\n"//& " \t uniform - uniform thickness layers evenly distributed \n"//& " \t\t between the surface and MAXIMUM_DEPTH. \n"//& @@ -289,51 +298,57 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & default="uniform", do_not_log=just_read) select case (trim(config)) case ("file") - call initialize_thickness_from_file(h, depth_tot, G, GV, US, PF, .false., just_read=just_read) + call initialize_thickness_from_file(dz, depth_tot, G, GV, US, PF, file_has_thickness=.false., & + mass_file=.false., just_read=just_read) case ("thickness_file") - call initialize_thickness_from_file(h, depth_tot, G, GV, US, PF, .true., just_read=just_read) + call initialize_thickness_from_file(dz, depth_tot, G, GV, US, PF, file_has_thickness=.true., & + mass_file=.false., just_read=just_read) + case ("mass_file") + call initialize_thickness_from_file(h, depth_tot, G, GV, US, PF, file_has_thickness=.true., & + mass_file=.true., just_read=just_read) + convert = .false. case ("coord") if (new_sim .and. useALE) then - call ALE_initThicknessToCoord( ALE_CSp, G, GV, h ) + call ALE_initThicknessToCoord( ALE_CSp, G, GV, dz, height_units=.true. ) elseif (new_sim) then call MOM_error(FATAL, "MOM_initialize_state: USE_REGRIDDING must be True "//& "for THICKNESS_CONFIG of 'coord'") endif - case ("uniform"); call initialize_thickness_uniform(h, depth_tot, G, GV, PF, & + case ("uniform"); call initialize_thickness_uniform(dz, depth_tot, G, GV, PF, & just_read=just_read) - case ("list"); call initialize_thickness_list(h, depth_tot, G, GV, US, PF, & + case ("list"); call initialize_thickness_list(dz, depth_tot, G, GV, US, PF, & just_read=just_read) - case ("DOME"); call DOME_initialize_thickness(h, depth_tot, G, GV, PF, & + case ("DOME"); call DOME_initialize_thickness(dz, depth_tot, G, GV, PF, & just_read=just_read) - case ("ISOMIP"); call ISOMIP_initialize_thickness(h, depth_tot, G, GV, US, PF, tv, & + case ("ISOMIP"); call ISOMIP_initialize_thickness(dz, depth_tot, G, GV, US, PF, tv, & just_read=just_read) - case ("benchmark"); call benchmark_initialize_thickness(h, depth_tot, G, GV, US, PF, & + case ("benchmark"); call benchmark_initialize_thickness(dz, depth_tot, G, GV, US, PF, & tv%eqn_of_state, tv%P_Ref, just_read=just_read) - case ("Neverworld","Neverland"); call Neverworld_initialize_thickness(h, depth_tot, & + case ("Neverworld","Neverland"); call Neverworld_initialize_thickness(dz, depth_tot, & G, GV, US, PF, tv%P_Ref) case ("search"); call initialize_thickness_search() - case ("circle_obcs"); call circle_obcs_initialize_thickness(h, depth_tot, G, GV, PF, & + case ("circle_obcs"); call circle_obcs_initialize_thickness(dz, depth_tot, G, GV, US, PF, & just_read=just_read) - case ("lock_exchange"); call lock_exchange_initialize_thickness(h, G, GV, US, & + case ("lock_exchange"); call lock_exchange_initialize_thickness(dz, G, GV, US, & PF, just_read=just_read) - case ("external_gwave"); call external_gwave_initialize_thickness(h, G, GV, US, & + case ("external_gwave"); call external_gwave_initialize_thickness(dz, G, GV, US, & PF, just_read=just_read) - case ("DOME2D"); call DOME2d_initialize_thickness(h, depth_tot, G, GV, US, PF, & + case ("DOME2D"); call DOME2d_initialize_thickness(dz, depth_tot, G, GV, US, PF, & just_read=just_read) - case ("adjustment2d"); call adjustment_initialize_thickness(h, G, GV, US, & + case ("adjustment2d"); call adjustment_initialize_thickness(dz, G, GV, US, & PF, just_read=just_read) - case ("sloshing"); call sloshing_initialize_thickness(h, depth_tot, G, GV, US, PF, & + case ("sloshing"); call sloshing_initialize_thickness(dz, depth_tot, G, GV, US, PF, & just_read=just_read) - case ("seamount"); call seamount_initialize_thickness(h, depth_tot, G, GV, US, PF, & + case ("seamount"); call seamount_initialize_thickness(dz, depth_tot, G, GV, US, PF, & just_read=just_read) - case ("dumbbell"); call dumbbell_initialize_thickness(h, depth_tot, G, GV, US, PF, & + case ("dumbbell"); call dumbbell_initialize_thickness(dz, depth_tot, G, GV, US, PF, & just_read=just_read) - case ("soliton"); call soliton_initialize_thickness(h, depth_tot, G, GV, US) - case ("phillips"); call Phillips_initialize_thickness(h, depth_tot, G, GV, US, PF, & + case ("soliton"); call soliton_initialize_thickness(dz, depth_tot, G, GV, US) + case ("phillips"); call Phillips_initialize_thickness(dz, depth_tot, G, GV, US, PF, & just_read=just_read) - case ("rossby_front"); call Rossby_front_initialize_thickness(h, G, GV, US, & + case ("rossby_front"); call Rossby_front_initialize_thickness(dz, G, GV, US, & PF, just_read=just_read) - case ("USER"); call user_initialize_thickness(h, G, GV, PF, & + case ("USER"); call user_initialize_thickness(dz, G, GV, PF, & just_read=just_read) case default ; call MOM_error(FATAL, "MOM_initialize_state: "//& "Unrecognized layer thickness configuration "//trim(config)) @@ -374,26 +389,26 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & G, GV, US, PF, just_read=just_read) case ("linear"); call initialize_temp_salt_linear(tv%T, tv%S, G, GV, US, PF, & just_read=just_read) - case ("DOME2D"); call DOME2d_initialize_temperature_salinity (tv%T, tv%S, h, & + case ("DOME2D"); call DOME2d_initialize_temperature_salinity (tv%T, tv%S, dz, & G, GV, US, PF, just_read=just_read) - case ("ISOMIP"); call ISOMIP_initialize_temperature_salinity (tv%T, tv%S, h, & + case ("ISOMIP"); call ISOMIP_initialize_temperature_salinity (tv%T, tv%S, dz, & depth_tot, G, GV, US, PF, eos, just_read=just_read) case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, & - tv%S, h, depth_tot, G, GV, US, PF, just_read=just_read) + tv%S, dz, depth_tot, G, GV, US, PF, just_read=just_read) case ("baroclinic_zone"); call baroclinic_zone_init_temperature_salinity( tv%T, & - tv%S, h, depth_tot, G, GV, US, PF, just_read=just_read) + tv%S, dz, depth_tot, G, GV, US, PF, just_read=just_read) case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, & - tv%S, h, G, GV, US, PF, just_read=just_read) + tv%S, dz, G, GV, US, PF, just_read=just_read) case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, & - tv%S, h, G, GV, US, PF, just_read=just_read) + tv%S, dz, G, GV, US, PF, just_read=just_read) case ("dumbbell"); call dumbbell_initialize_temperature_salinity(tv%T, & - tv%S, h, G, GV, US, PF, just_read=just_read) + tv%S, dz, G, GV, US, PF, just_read=just_read) case ("rossby_front"); call Rossby_front_initialize_temperature_salinity ( tv%T, & - tv%S, h, G, GV, US, PF, just_read=just_read) - case ("SCM_CVMix_tests"); call SCM_CVMix_tests_TS_init(tv%T, tv%S, h, & + tv%S, dz, G, GV, US, PF, just_read=just_read) + case ("SCM_CVMix_tests"); call SCM_CVMix_tests_TS_init(tv%T, tv%S, dz, & G, GV, US, PF, just_read=just_read) case ("dense"); call dense_water_initialize_TS(G, GV, US, PF, tv%T, tv%S, & - h, just_read=just_read) + dz, just_read=just_read) case ("USER"); call user_init_temperature_salinity(tv%T, tv%S, G, GV, PF, & just_read=just_read) case default ; call MOM_error(FATAL, "MOM_initialize_state: "//& @@ -404,8 +419,10 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & if (use_temperature .and. use_OBC) & call fill_temp_salt_segments(G, GV, US, OBC, tv) - ! Calculate the initial surface displacement under ice shelf + ! Convert thicknesses from geometric distances in depth units to thickness units or mass-per-unit-area. + if (new_sim .and. convert) call convert_thickness(dz, h, G, GV, US, tv) + ! Handle the initial surface displacement under ice shelf call get_param(PF, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, & "If true, depress the initial surface to avoid huge "//& "tsunamis when a large surface pressure is applied.", & @@ -415,10 +432,43 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & "at the depth where the hydrostatic pressure matches the imposed "//& "surface pressure which is read from file.", default=.false., & do_not_log=just_read) + if (depress_sfc .and. trim_ic_for_p_surf) call MOM_error(FATAL, "MOM_initialize_state: "//& + "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True") - if (new_sim) then - if (use_ice_shelf .and. present(mass_shelf) .and. .not. (trim_ic_for_p_surf .or. depress_sfc)) & - call calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) + if (new_sim .and. debug .and. (depress_sfc .or. trim_ic_for_p_surf)) & + call hchksum(h, "Pre-depress: h ", G%HI, haloshift=1, scale=GV%H_to_MKS) + + ! Remove the mass that would be displaced by an ice shelf or inverse barometer. + if (depress_sfc) then + call depress_surface(h, G, GV, US, PF, tv, just_read=just_read) + elseif (trim_ic_for_p_surf) then + call trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read=just_read) + elseif (new_sim .and. use_ice_shelf .and. present(mass_shelf)) then + call calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h) + endif + + ! Perhaps we want to run the regridding coordinate generator for multiple + ! iterations here so the initial grid is consistent with the coordinate + if (useALE) then + call get_param(PF, mdl, "REGRID_ACCELERATE_INIT", regrid_accelerate, & + "If true, runs REGRID_ACCELERATE_ITERATIONS iterations of the regridding "//& + "algorithm to push the initial grid to be consistent with the initial "//& + "condition. Useful only for state-based and iterative coordinates.", & + default=.false., do_not_log=just_read) + if (regrid_accelerate) then + call get_param(PF, mdl, "REGRID_ACCELERATE_ITERATIONS", regrid_iterations, & + "The number of regridding iterations to perform to generate "//& + "an initial grid that is consistent with the initial conditions.", & + default=1, do_not_log=just_read) + + call get_param(PF, mdl, "DT", dt, "Timestep", & + units="s", scale=US%s_to_T, fail_if_missing=.true.) + + if (new_sim .and. debug) & + call hchksum(h, "Pre-ALE_regrid: h ", G%HI, haloshift=1, scale=GV%H_to_MKS) + call ALE_regrid_accelerated(ALE_CSp, G, GV, h, tv, regrid_iterations, u, v, OBC, tracer_Reg, & + dt=dt, initial=.true.) + endif endif ! The thicknesses in halo points might be needed to initialize the velocities. @@ -438,21 +488,15 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & " \t USER - call a user modified routine.", default="zero", & do_not_log=just_read) select case (trim(config)) - case ("file"); call initialize_velocity_from_file(u, v, G, GV, US, PF, & - just_read=just_read) - case ("zero"); call initialize_velocity_zero(u, v, G, GV, PF, & - just_read=just_read) - case ("uniform"); call initialize_velocity_uniform(u, v, G, GV, US, PF, & - just_read=just_read) - case ("circular"); call initialize_velocity_circular(u, v, G, GV, US, PF, & - just_read=just_read) - case ("phillips"); call Phillips_initialize_velocity(u, v, G, GV, US, PF, & - just_read=just_read) + case ("file"); call initialize_velocity_from_file(u, v, G, GV, US, PF, just_read) + case ("zero"); call initialize_velocity_zero(u, v, G, GV, PF, just_read) + case ("uniform"); call initialize_velocity_uniform(u, v, G, GV, US, PF, just_read) + case ("circular"); call initialize_velocity_circular(u, v, G, GV, US, PF, just_read) + case ("phillips"); call Phillips_initialize_velocity(u, v, G, GV, US, PF, just_read) case ("rossby_front"); call Rossby_front_initialize_velocity(u, v, h, & - G, GV, US, PF, just_read=just_read) - case ("soliton"); call soliton_initialize_velocity(u, v, h, G, GV, US) - case ("USER"); call user_initialize_velocity(u, v, G, GV, US, PF, & - just_read=just_read) + G, GV, US, PF, just_read) + case ("soliton"); call soliton_initialize_velocity(u, v, G, GV, US) + case ("USER"); call user_initialize_velocity(u, v, G, GV, US, PF, just_read) case default ; call MOM_error(FATAL, "MOM_initialize_state: "//& "Unrecognized velocity configuration "//trim(config)) end select @@ -462,49 +506,8 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & call uvchksum("MOM_initialize_state [uv]", u, v, G%HI, haloshift=1, scale=US%L_T_to_m_s) endif - ! Optionally convert the thicknesses from m to kg m-2. This is particularly - ! useful in a non-Boussinesq model. - call get_param(PF, mdl, "CONVERT_THICKNESS_UNITS", convert, & - "If true, convert the thickness initial conditions from "//& - "units of m to kg m-2 or vice versa, depending on whether "//& - "BOUSSINESQ is defined. This does not apply if a restart "//& - "file is read.", default=.not.GV%Boussinesq, do_not_log=just_read) - - if (new_sim .and. convert .and. .not.GV%Boussinesq) & - ! Convert thicknesses from geometric distances to mass-per-unit-area. - call convert_thickness(h, G, GV, US, tv) - - ! Remove the mass that would be displaced by an ice shelf or inverse barometer. - if (depress_sfc .and. trim_ic_for_p_surf) call MOM_error(FATAL, "MOM_initialize_state: "//& - "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True") - if (new_sim .and. debug .and. (depress_sfc .or. trim_ic_for_p_surf)) & - call hchksum(h, "Pre-depress: h ", G%HI, haloshift=1, scale=GV%H_to_m) - if (depress_sfc) call depress_surface(h, G, GV, US, PF, tv, just_read=just_read) - if (trim_ic_for_p_surf) call trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read=just_read) - - ! Perhaps we want to run the regridding coordinate generator for multiple - ! iterations here so the initial grid is consistent with the coordinate - if (useALE) then - call get_param(PF, mdl, "REGRID_ACCELERATE_INIT", regrid_accelerate, & - "If true, runs REGRID_ACCELERATE_ITERATIONS iterations of the regridding "//& - "algorithm to push the initial grid to be consistent with the initial "//& - "condition. Useful only for state-based and iterative coordinates.", & - default=.false., do_not_log=just_read) - if (regrid_accelerate) then - call get_param(PF, mdl, "REGRID_ACCELERATE_ITERATIONS", regrid_iterations, & - "The number of regridding iterations to perform to generate "//& - "an initial grid that is consistent with the initial conditions.", & - default=1, do_not_log=just_read) - - call get_param(PF, mdl, "DT", dt, "Timestep", & - units="s", scale=US%s_to_T, fail_if_missing=.true.) - - if (new_sim .and. debug) & - call hchksum(h, "Pre-ALE_regrid: h ", G%HI, haloshift=1, scale=GV%H_to_m) - call ALE_regrid_accelerated(ALE_CSp, G, GV, h, tv, regrid_iterations, u, v, OBC, tracer_Reg, & - dt=dt, initial=.true.) - endif - endif + ! This is the end of the block of code that might have initialized fields + ! internally at the start of a new run. ! Initialized assimilative incremental update (oda_incupd) structure and ! register restart. @@ -517,9 +520,6 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & call restart_registry_lock(restart_CS) endif - ! This is the end of the block of code that might have initialized fields - ! internally at the start of a new run. - if (.not.new_sim) then ! This block restores the state from a restart file. ! This line calls a subroutine that reads the initial conditions ! from a previously generated file. @@ -548,7 +548,7 @@ subroutine MOM_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, & call pass_var(h, G%Domain) if (debug) then - call hchksum(h, "MOM_initialize_state: h ", G%HI, haloshift=1, scale=GV%H_to_m) + call hchksum(h, "MOM_initialize_state: h ", G%HI, haloshift=1, scale=GV%H_to_MKS) if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", G%HI, haloshift=1, scale=US%C_to_degC) if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", G%HI, haloshift=1, scale=US%S_to_ppt) if ( use_temperature .and. debug_layers) then ; do k=1,nz @@ -667,12 +667,14 @@ end subroutine MOM_initialize_state !> Reads the layer thicknesses or interface heights from a file. subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, file_has_thickness, & - just_read) + just_read, mass_file) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized, in height + !! or thickness units, depending on the value of + !! mass_file [Z ~> m] or [H ~> m or kg m-2]. real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -682,6 +684,8 @@ subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, f !! interface heights. logical, intent(in) :: just_read !< If true, this call will only read !! parameters without changing h. + logical, intent(in) :: mass_file !< If true, this file contains layer thicknesses in + !! units of mass per unit area. ! Local variables real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Interface heights, in depth units [Z ~> m]. @@ -723,12 +727,17 @@ subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, f "The variable name for layer thickness initial conditions.", & default="h", do_not_log=just_read) call get_param(param_file, mdl, "THICKNESS_IC_RESCALE", h_rescale, & - "A factor by which to rescale the initial thicknesses in the input "//& - "file to convert them to units of m.", & + 'A factor by which to rescale the initial thicknesses in the input file to '//& + 'convert them to units of kg/m2 (if THICKNESS_CONFIG="mass_file") or m.', & default=1.0, units="various", do_not_log=just_read) if (just_read) return ! All run-time parameters have been read, so return. - call MOM_read_data(filename, h_var, h(:,:,:), G%Domain, scale=h_rescale*GV%m_to_H) + if (mass_file) then + h_rescale = h_rescale*GV%kg_m2_to_H + else + h_rescale = h_rescale*US%m_to_Z + endif + call MOM_read_data(filename, h_var, h(:,:,:), G%Domain, scale=h_rescale) else call get_param(param_file, mdl, "ADJUST_THICKNESS", correct_thickness, & "If true, all mass below the bottom removed if the "//& @@ -763,9 +772,9 @@ subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, f do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) then eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + h(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (eta(i,j,K) - eta(i,j,K+1)) + h(i,j,k) = eta(i,j,K) - eta(i,j,K+1) endif enddo ; enddo ; enddo @@ -798,7 +807,7 @@ subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h, ht, dZ_ref_eta) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: eta !< Interface heights [Z ~> m]. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [Z ~> m] real, intent(in) :: ht !< Tolerance to exceed adjustment !! criteria [Z ~> m] real, optional, intent(in) :: dZ_ref_eta !< The difference between the @@ -857,10 +866,6 @@ subroutine adjustEtaToFitBathymetry(G, GV, US, eta, h, ht, dZ_ref_eta) endif enddo ; enddo - ! Now convert thicknesses to units of H. - do k=1,nz ; do j=js,je ; do i=is,ie - h(i,j,k) = h(i,j,k)*GV%Z_to_H - enddo ; enddo ; enddo call sum_across_PEs(dilations) if ((dilations > 0) .and. (is_root_pe())) then @@ -876,7 +881,7 @@ subroutine initialize_thickness_uniform(h, depth_tot, G, GV, param_file, just_re type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -915,9 +920,9 @@ subroutine initialize_thickness_uniform(h, depth_tot, G, GV, param_file, just_re eta1D(K) = e0(K) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then eta1D(K) = eta1D(K+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + h(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (eta1D(K) - eta1D(K+1)) + h(i,j,k) = eta1D(K) - eta1D(K+1) endif enddo enddo ; enddo @@ -929,9 +934,9 @@ end subroutine initialize_thickness_uniform subroutine initialize_thickness_list(h, depth_tot, G, GV, US, param_file, just_read) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -990,9 +995,9 @@ subroutine initialize_thickness_list(h, depth_tot, G, GV, US, param_file, just_r eta1D(K) = e0(K) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then eta1D(K) = eta1D(K+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + h(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (eta1D(K) - eta1D(K+1)) + h(i,j,k) = eta1D(K) - eta1D(K+1) endif enddo enddo ; enddo @@ -1005,14 +1010,17 @@ subroutine initialize_thickness_search call MOM_error(FATAL," MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED") end subroutine initialize_thickness_search -!> Converts thickness from geometric to pressure units -subroutine convert_thickness(h, G, GV, US, tv) +!> Converts thickness from geometric height units to thickness units +subroutine convert_thickness(dz, h, G, GV, US, tv) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(inout) :: h !< Input geometric layer thicknesses being converted - !! to layer pressure [H ~> m or kg m-2]. + intent(in) :: dz !< Input geometric layer thicknesses [Z ~> m]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & + intent(inout) :: h !< Output thicknesses in thickness units [H ~> m or kg m-2]. + !! This is essentially intent out, but declared as intent + !! inout to preserve any initalized values in halo points. type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables ! Local variables @@ -1022,8 +1030,6 @@ subroutine convert_thickness(h, G, GV, US, tv) real :: rho(SZI_(G)) ! The in situ density [R ~> kg m-3] real :: I_gEarth ! Unit conversion factors divided by the gravitational acceleration ! [H T2 R-1 L-2 ~> s2 m2 kg-1 or s2 m-1] - real :: HR_to_pres ! A conversion factor from the input geometric thicknesses times the layer - ! densities into pressure units [L2 T-2 H-1 ~> m s-2 or m4 kg-1 s-2]. 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 integer :: itt, max_itt @@ -1033,10 +1039,11 @@ subroutine convert_thickness(h, G, GV, US, tv) max_itt = 10 if (GV%Boussinesq) then - call MOM_error(FATAL,"Not yet converting thickness with Boussinesq approx.") + do k=1,nz ; do j=js,je ; do i=is,ie + h(i,j,k) = GV%Z_to_H * dz(i,j,k) + enddo ; enddo ; enddo else I_gEarth = GV%RZ_to_H / GV%g_Earth - HR_to_pres = GV%g_Earth * GV%H_to_Z if (associated(tv%eqn_of_state)) then do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 @@ -1049,7 +1056,8 @@ subroutine convert_thickness(h, G, GV, US, tv) call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_top(:,j), rho, & tv%eqn_of_state, EOSdom) do i=is,ie - p_bot(i,j) = p_top(i,j) + HR_to_pres * (h(i,j,k) * rho(i)) + ! This could be simplified, but it would change answers at roundoff. + p_bot(i,j) = p_top(i,j) + (GV%g_Earth*GV%H_to_Z) * ((GV%Z_to_H*dz(i,j,k)) * rho(i)) enddo enddo @@ -1062,7 +1070,7 @@ subroutine convert_thickness(h, G, GV, US, tv) ! Use Newton's method to correct the bottom value. ! The hydrostatic equation is sufficiently linear that no bounds-checking is needed. do i=is,ie - p_bot(i,j) = p_bot(i,j) + rho(i) * (HR_to_pres*h(i,j,k) - dz_geo(i,j)) + p_bot(i,j) = p_bot(i,j) + rho(i) * ((GV%g_Earth*GV%H_to_Z)*(GV%Z_to_H*dz(i,j,k)) - dz_geo(i,j)) enddo enddo ; endif enddo @@ -1073,7 +1081,7 @@ subroutine convert_thickness(h, G, GV, US, tv) enddo else do k=1,nz ; do j=js,je ; do i=is,ie - h(i,j,k) = h(i,j,k) * (GV%Rlay(k) / GV%Rho0) + h(i,j,k) = (GV%Z_to_H*dz(i,j,k)) * (GV%Rlay(k) / GV%Rho0) enddo ; enddo ; enddo endif endif @@ -2503,7 +2511,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just real, dimension(:,:,:), allocatable, target :: salt_z ! Input salinities [S ~> ppt] real, dimension(:,:,:), allocatable, target :: mask_z ! 1 for valid data points [nondim] real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zi ! Interface heights [Z ~> m]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zi ! Interface heights [Z ~> m] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Layer thicknesses in height units [Z ~> m] real, dimension(SZI_(G),SZJ_(G)) :: Z_bottom ! The (usually negative) height of the seafloor ! relative to the surface [Z ~> m]. integer, dimension(SZI_(G),SZJ_(G)) :: nlevs ! The number of levels in each column with valid data @@ -2514,7 +2523,8 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just real, dimension(:,:,:), allocatable, target :: tmpT1dIn ! Input temperatures on a model-sized grid [C ~> degC] real, dimension(:,:,:), allocatable, target :: tmpS1dIn ! Input salinities on a model-sized grid [S ~> ppt] real, dimension(:,:,:), allocatable :: tmp_mask_in ! The valid data mask on a model-sized grid [nondim] - real, dimension(:,:,:), allocatable :: h1 ! Thicknesses [H ~> m or kg m-2]. + real, dimension(:,:,:), allocatable :: dz1 ! Input grid thicknesses in depth units [Z ~> m] + real, dimension(:,:,:), allocatable :: h1 ! Thicknesses on the input grid [H ~> m or kg m-2]. real, dimension(:,:,:), allocatable :: dz_interface ! Change in position of interface due to ! regridding [H ~> m or kg m-2] real :: zTopOfCell, zBottomOfCell ! Heights in Z units [Z ~> m]. @@ -2721,7 +2731,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just if ((.not.useALEremapping) .and. adjust_temperature) & ! This call is just here to read and log the determine_temperature parameters call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), eos, tv%P_Ref, 0, & - h, 0, G, GV, US, PF, just_read=.true.) + 0, G, GV, US, PF, just_read=.true.) call cpu_clock_end(id_clock_routine) return ! All run-time parameters have been read, so return. endif @@ -2773,6 +2783,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! Build the source grid and copy data onto model-shaped arrays with vanished layers allocate( tmp_mask_in(isd:ied,jsd:jed,nkd), source=0.0 ) + allocate( dz1(isd:ied,jsd:jed,nkd), source=0.0 ) allocate( h1(isd:ied,jsd:jed,nkd), source=0.0 ) allocate( tmpT1dIn(isd:ied,jsd:jed,nkd), source=0.0 ) allocate( tmpS1dIn(isd:ied,jsd:jed,nkd), source=0.0 ) @@ -2793,10 +2804,10 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just tmpT1dIn(i,j,k) = temp_land_fill tmpS1dIn(i,j,k) = salt_land_fill endif - h1(i,j,k) = GV%Z_to_H * (zTopOfCell - zBottomOfCell) + dz1(i,j,k) = (zTopOfCell - zBottomOfCell) zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k enddo - h1(i,j,kd) = h1(i,j,kd) + GV%Z_to_H * max(0., zTopOfCell - Z_bottom(i,j) ) + dz1(i,j,kd) = dz1(i,j,kd) + max(0., zTopOfCell - Z_bottom(i,j) ) ! The max here is in case the data data is shallower than model endif ! mask2dT enddo ; enddo @@ -2811,20 +2822,27 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just allocate( hTarget(nz) ) hTarget = getCoordinateResolution( regridCS ) do j = js, je ; do i = is, ie - h(i,j,:) = 0. + dz(i,j,:) = 0. if (G%mask2dT(i,j) > 0.) then ! Build the target grid combining hTarget and topography zTopOfCell = 0. ; zBottomOfCell = 0. do k = 1, nz zBottomOfCell = max( zTopOfCell - hTarget(k), Z_bottom(i,j)) - h(i,j,k) = GV%Z_to_H * (zTopOfCell - zBottomOfCell) + dz(i,j,k) = zTopOfCell - zBottomOfCell zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k enddo else - h(i,j,:) = 0. + dz(i,j,:) = 0. endif ! mask2dT enddo ; enddo deallocate( hTarget ) + + do k=1,nkd ; do j=js,je ; do i=is,ie + h1(i,j,k) = GV%Z_to_H*dz1(i,j,k) + enddo ; enddo ; enddo + do k=1,nz ; do j=js,je ; do i=is,ie + h(i,j,k) = GV%Z_to_H*dz(i,j,k) + enddo ; enddo ; enddo endif ! Now remap from source grid to target grid, first setting reconstruction parameters @@ -2838,6 +2856,9 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just GV_loc%ke = nkd allocate( dz_interface(isd:ied,jsd:jed,nkd+1) ) ! Need for argument to regridding_main() but is not used + ! Convert thicknesses to units of H. + call convert_thickness(dz1, h1, G, GV_loc, US, tv_loc) + call regridding_preadjust_reqs(regridCS, do_conv_adj, ignore) if (do_conv_adj) call convective_adjustment(G, GV_loc, h1, tv_loc) call regridding_main( remapCS, regridCS, G, GV_loc, h1, tv_loc, h, dz_interface, & @@ -2850,6 +2871,7 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpS1dIn, h, tv%S, all_cells=remap_full_column, & old_remap=remap_old_alg, answer_date=remap_answer_date ) + deallocate( dz1 ) deallocate( h1 ) deallocate( tmpT1dIn ) deallocate( tmpS1dIn ) @@ -2886,15 +2908,16 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just deallocate(rho_z) + dz(:,:,:) = 0.0 if (correct_thickness) then - call adjustEtaToFitBathymetry(G, GV, US, zi, h, h_tolerance, dZ_ref_eta=G%Z_ref) + call adjustEtaToFitBathymetry(G, GV, US, zi, dz, h_tolerance, dZ_ref_eta=G%Z_ref) else do k=nz,1,-1 ; do j=js,je ; do i=is,ie if (zi(i,j,K) < (zi(i,j,K+1) + GV%Angstrom_Z)) then zi(i,j,K) = zi(i,j,K+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + dz(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (zi(i,j,K) - zi(i,j,K+1)) + dz(i,j,k) = zi(i,j,K) - zi(i,j,K+1) endif enddo ; enddo ; enddo inconsistent = 0 @@ -2926,9 +2949,12 @@ subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just ! Finally adjust to target density ks = 1 ; if (separate_mixed_layer) ks = GV%nk_rho_varies + 1 call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), eos, tv%P_Ref, niter, & - h, ks, G, GV, US, PF, just_read) + ks, G, GV, US, PF, just_read) endif + ! Now convert thicknesses to units of H. + call convert_thickness(dz, h, G, GV, US, tv) + endif ! useALEremapping deallocate(z_in, z_edges_in, temp_z, salt_z, mask_z) diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index c089181c16..fab7da3917 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/src/tracer/MOM_tracer_Z_init.F90 @@ -556,8 +556,8 @@ end function find_limited_slope !> This subroutine determines the potential temperature and salinity that !! is consistent with the target density using provided initial guess -subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_start, G, GV, US, & - PF, just_read, h_massless) +subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, k_start, G, GV, US, PF, & + just_read) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & @@ -565,20 +565,15 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(inout) :: salt !< salinity [S ~> ppt] real, dimension(SZK_(GV)), intent(in) :: R_tgt !< desired potential density [R ~> kg m-3]. - type(EOS_type), intent(in) :: EOS !< seawater equation of state control structure + type(EOS_type), intent(in) :: EOS !< seawater equation of state control structure real, intent(in) :: p_ref !< reference pressure [R L2 T-2 ~> Pa]. integer, intent(in) :: niter !< maximum number of iterations integer, intent(in) :: k_start !< starting index (i.e. below the buffer layer) - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< layer thickness, used only to avoid working on - !! massless layers [H ~> m or kg m-2] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: PF !< A structure indicating the open file !! to parse for model parameter values. logical, intent(in) :: just_read !< If true, this call will only read !! parameters without changing T or S. - real, optional, intent(in) :: h_massless !< A threshold below which a layer is - !! determined to be massless [H ~> m or kg m-2] ! Local variables (All of which need documentation!) real, dimension(SZI_(G),SZK_(GV)) :: & @@ -587,7 +582,6 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star dT, & ! An estimated change in temperature before bounding [C ~> degC] dS, & ! An estimated change in salinity before bounding [S ~> ppt] rho, & ! Layer densities with the current estimate of temperature and salinity [R ~> kg m-3] - hin, & ! A 2D copy of the layer thicknesses [H ~> m or kg m-2] drho_dT, & ! Partial derivative of density with temperature [R C-1 ~> kg m-3 degC-1] drho_dS ! Partial derivative of density with salinity [R S-1 ~> kg m-3 ppt-1] real, dimension(SZI_(G)) :: press ! Reference pressures [R L2 T-2 ~> Pa] @@ -675,7 +669,6 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star dS(:,:) = 0. ! Needs to be zero everywhere since there is a maxval(abs(dS)) later... T(:,:) = temp(:,j,:) S(:,:) = salt(:,j,:) - hin(:,:) = h(:,j,:) dT(:,:) = 0.0 adjust_salt = .true. iter_loop: do itt = 1,niter @@ -685,7 +678,7 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star EOS, EOSdom ) enddo do k=k_start,nz ; do i=is,ie -! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln) then +! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln) then if (abs(rho(i,k)-R_tgt(k))>tol_rho) then if (.not.fit_together) then dT(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dT(i,k), max_t_adj), -max_t_adj) @@ -713,7 +706,7 @@ subroutine determine_temperature(temp, salt, R_tgt, EOS, p_ref, niter, h, k_star EOS, EOSdom ) enddo do k=k_start,nz ; do i=is,ie -! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. hin(i,k)>h_massless .and. abs(T(i,k)-land_fill) < epsln ) then +! if (abs(rho(i,k)-R_tgt(k))>tol_rho .and. abs(T(i,k)-land_fill) < epsln ) then if (abs(rho(i,k)-R_tgt(k)) > tol_rho) then dS(i,k) = max(min((R_tgt(k)-rho(i,k)) / drho_dS(i,k), max_s_adj), -max_s_adj) S(i,k) = max(min(S(i,k)+dS(i,k), S_max), S_min) diff --git a/src/user/DOME2d_initialization.F90 b/src/user/DOME2d_initialization.F90 index 1382fe8e34..5cc63e734f 100644 --- a/src/user/DOME2d_initialization.F90 +++ b/src/user/DOME2d_initialization.F90 @@ -98,7 +98,7 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -158,16 +158,16 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then eta1D(k) = eta1D(k+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + h(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo x = ( G%geoLonT(i,j) - G%west_lon ) / G%len_lon if ( x <= dome2d_width_bay ) then - h(i,j,1:nz-1) = GV%Angstrom_H - h(i,j,nz) = GV%Z_to_H * dome2d_depth_bay * G%max_depth - (nz-1) * GV%Angstrom_H + h(i,j,1:nz-1) = GV%Angstrom_Z + h(i,j,nz) = dome2d_depth_bay * G%max_depth - (nz-1) * GV%Angstrom_Z endif enddo ; enddo @@ -180,16 +180,16 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju ! eta1D(k) = e0(k) ! if (eta1D(k) < (eta1D(k+1) + min_thickness)) then ! eta1D(k) = eta1D(k+1) + min_thickness - ! h(i,j,k) = GV%Z_to_H * min_thickness + ! h(i,j,k) = min_thickness ! else - ! h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + ! h(i,j,k) = eta1D(k) - eta1D(k+1) ! endif ! enddo ! ! x = G%geoLonT(i,j) / G%len_lon ! if ( x <= dome2d_width_bay ) then - ! h(i,j,1:nz-1) = GV%Z_to_H * min_thickness - ! h(i,j,nz) = GV%Z_to_H * (dome2d_depth_bay * G%max_depth - (nz-1) * min_thickness) + ! h(i,j,1:nz-1) = min_thickness + ! h(i,j,nz) = dome2d_depth_bay * G%max_depth - (nz-1) * min_thickness ! endif ! ! enddo ; enddo @@ -202,16 +202,16 @@ subroutine DOME2d_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ju eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + min_thickness)) then eta1D(k) = eta1D(k+1) + min_thickness - h(i,j,k) = GV%Z_to_H * min_thickness + h(i,j,k) = min_thickness else - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo enddo ; enddo case ( REGRIDDING_SIGMA ) do j=js,je ; do i=is,ie - h(i,j,:) = GV%Z_to_H*depth_tot(i,j) / nz + h(i,j,:) = depth_tot(i,j) / nz enddo ; enddo case default @@ -225,11 +225,11 @@ end subroutine DOME2d_initialize_thickness !> Initialize temperature and salinity in the 2d DOME configuration subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_file, just_read) - type(ocean_grid_type), intent(in) :: G !< Ocean grid structure + type(ocean_grid_type), intent(in) :: G !< Ocean grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> ppt] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> ppt] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Parameter file structure logical, intent(in) :: just_read !< If true, this call will @@ -287,7 +287,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi do j=js,je ; do i=is,ie xi0 = 0.0 do k = 1,nz - xi1 = xi0 + (GV%H_to_Z * h(i,j,k)) / G%max_depth + xi1 = xi0 + h(i,j,k) / G%max_depth S(i,j,k) = S_surf + 0.5 * S_range * (xi0 + xi1) xi0 = xi1 enddo @@ -298,7 +298,7 @@ subroutine DOME2d_initialize_temperature_salinity ( T, S, h, G, GV, US, param_fi do j=js,je ; do i=is,ie xi0 = 0.0 do k = 1,nz - xi1 = xi0 + (GV%H_to_Z * h(i,j,k)) / G%max_depth + xi1 = xi0 + h(i,j,k) / G%max_depth S(i,j,k) = S_surf + 0.5 * S_range * (xi0 + xi1) xi0 = xi1 enddo diff --git a/src/user/DOME_initialization.F90 b/src/user/DOME_initialization.F90 index 7f939ffef6..4a12387d9d 100644 --- a/src/user/DOME_initialization.F90 +++ b/src/user/DOME_initialization.F90 @@ -105,7 +105,7 @@ subroutine DOME_initialize_thickness(h, depth_tot, G, GV, param_file, just_read) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -141,9 +141,9 @@ subroutine DOME_initialize_thickness(h, depth_tot, G, GV, param_file, just_read) eta1D(K) = e0(K) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then eta1D(K) = eta1D(K+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + h(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (eta1D(K) - eta1D(K+1)) + h(i,j,k) = eta1D(K) - eta1D(K+1) endif enddo enddo ; enddo diff --git a/src/user/ISOMIP_initialization.F90 b/src/user/ISOMIP_initialization.F90 index bba357f490..7e3299b372 100644 --- a/src/user/ISOMIP_initialization.F90 +++ b/src/user/ISOMIP_initialization.F90 @@ -143,7 +143,7 @@ subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -170,7 +170,7 @@ subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke if (.not.just_read) & - call MOM_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness") + call MOM_mesg("ISOMIP_initialization.F90, ISOMIP_initialize_thickness: setting thickness") call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, & 'Minimum layer thickness', units='m', default=1.e-3, do_not_log=just_read, scale=US%m_to_Z) @@ -225,9 +225,9 @@ subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then eta1D(k) = eta1D(k+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + h(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo enddo ; enddo @@ -240,9 +240,9 @@ subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv eta1D(k) = -G%max_depth * real(k-1) / real(nz) if (eta1D(k) < (eta1D(k+1) + min_thickness)) then eta1D(k) = eta1D(k+1) + min_thickness - h(i,j,k) = GV%Z_to_H * min_thickness + h(i,j,k) = min_thickness else - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo enddo ; enddo @@ -250,7 +250,7 @@ subroutine ISOMIP_initialize_thickness ( h, depth_tot, G, GV, US, param_file, tv case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie - h(i,j,:) = GV%Z_to_H * depth_tot(i,j) / real(nz) + h(i,j,:) = depth_tot(i,j) / real(nz) enddo ; enddo case default @@ -269,7 +269,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, U type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> ppt] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), intent(in) :: depth_tot !< The nominal total bottom-to-top !! depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< Parameter file structure @@ -334,10 +334,10 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, U do j=js,je ; do i=is,ie xi0 = -depth_tot(i,j) do k = nz,1,-1 - xi0 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z ! Depth in middle of layer + xi0 = xi0 + 0.5 * h(i,j,k) ! Depth in middle of layer S(i,j,k) = S_sur + dS_dz * xi0 T(i,j,k) = T_sur + dT_dz * xi0 - xi0 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z ! Depth at top of layer + xi0 = xi0 + 0.5 * h(i,j,k) ! Depth at top of layer enddo enddo ; enddo @@ -372,10 +372,10 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, U xi0 = 0.0 do k = 1,nz !T0(k) = T_Ref; S0(k) = S_Ref - xi1 = xi0 + 0.5 * h(i,j,k) * GV%H_to_Z + xi1 = xi0 + 0.5 * h(i,j,k) S0(k) = S_sur - dS_dz * xi1 T0(k) = T_sur - dT_dz * xi1 - xi0 = xi0 + h(i,j,k) * GV%H_to_Z + xi0 = xi0 + h(i,j,k) ! write(mesg,*) 'S,T,xi0,xi1,k',S0(k),T0(k),xi0,xi1,k ! call MOM_mesg(mesg,5) enddo @@ -430,7 +430,7 @@ subroutine ISOMIP_initialize_temperature_salinity ( T, S, h, depth_tot, G, GV, U !i=G%iec; j=G%jec !do k = 1,nz ! call calculate_density(T(i,j,k), S(i,j,k),0.0,rho_tmp,eqn_of_state, scale=US%kg_m3_to_R) - ! write(mesg,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),US%C_to_degC*T(i,j,k),US%S_to_ppt*S(i,j,k),rho_tmp,GV%Rlay(k) + ! write(mesg,*) 'k,h,T,S,rho,Rlay',k,US%Z_to_m*h(i,j,k),US%C_to_degC*T(i,j,k),US%S_to_ppt*S(i,j,k),rho_tmp,GV%Rlay(k) ! call MOM_mesg(mesg,5) !enddo diff --git a/src/user/Neverworld_initialization.F90 b/src/user/Neverworld_initialization.F90 index fcd40cf8da..05de663d46 100644 --- a/src/user/Neverworld_initialization.F90 +++ b/src/user/Neverworld_initialization.F90 @@ -243,7 +243,7 @@ subroutine Neverworld_initialize_thickness(h, depth_tot, G, GV, US, param_file, type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< The thickness that is being - !! initialized [H ~> m or kg m-2]. + !! initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open @@ -288,12 +288,12 @@ subroutine Neverworld_initialize_thickness(h, depth_tot, G, GV, US, param_file, do j=js,je ; do i=is,ie e_interface = -depth_tot(i,j) do k=nz,2,-1 - h(i,j,k) = GV%Z_to_H * (e0(k) - e_interface) ! Nominal thickness + h(i,j,k) = e0(k) - e_interface ! Nominal thickness x = (G%geoLonT(i,j)-G%west_lon)/G%len_lon y = (G%geoLatT(i,j)-G%south_lat)/G%len_lat r1 = sqrt((x-0.7)**2+(y-0.2)**2) r2 = sqrt((x-0.3)**2+(y-0.25)**2) - h(i,j,k) = h(i,j,k) + pert_amp * (e0(k) - e0(nz+1)) * GV%Z_to_H * & + h(i,j,k) = h(i,j,k) + pert_amp * (e0(k) - e0(nz+1)) * & (spike(r1,0.15)-spike(r2,0.15)) ! Prescribed perturbation if (h_noise /= 0.) then rns = initializeRandomNumberStream( int( 4096*(x + (y+1.)) ) ) @@ -301,11 +301,11 @@ subroutine Neverworld_initialize_thickness(h, depth_tot, G, GV, US, param_file, noise = h_noise * 2. * ( noise - 0.5 ) ! range -h_noise to h_noise h(i,j,k) = ( 1. + noise ) * h(i,j,k) endif - h(i,j,k) = max( GV%Angstrom_H, h(i,j,k) ) ! Limit to non-negative - e_interface = e_interface + GV%H_to_Z * h(i,j,k) ! Actual position of upper interface + h(i,j,k) = max( GV%Angstrom_Z, h(i,j,k) ) ! Limit to non-negative + e_interface = e_interface + h(i,j,k) ! Actual position of upper interface enddo - h(i,j,1) = GV%Z_to_H * (e0(1) - e_interface) ! Nominal thickness - h(i,j,1) = max( GV%Angstrom_H, h(i,j,1) ) ! Limit to non-negative + h(i,j,1) = e0(1) - e_interface ! Nominal thickness + h(i,j,1) = max( GV%Angstrom_Z, h(i,j,1) ) ! Limit to non-negative enddo ; enddo end subroutine Neverworld_initialize_thickness diff --git a/src/user/Phillips_initialization.F90 b/src/user/Phillips_initialization.F90 index 62b55bb0a1..e0d2cafeae 100644 --- a/src/user/Phillips_initialization.F90 +++ b/src/user/Phillips_initialization.F90 @@ -39,7 +39,7 @@ subroutine Phillips_initialize_thickness(h, depth_tot, G, GV, US, param_file, ju type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2] + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -116,9 +116,9 @@ subroutine Phillips_initialize_thickness(h, depth_tot, G, GV, US, param_file, ju eta1D(K) = eta_im(j,K) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then eta1D(K) = eta1D(K+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + h(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (eta1D(K) - eta1D(K+1)) + h(i,j,k) = eta1D(K) - eta1D(K+1) endif enddo enddo ; enddo diff --git a/src/user/Rossby_front_2d_initialization.F90 b/src/user/Rossby_front_2d_initialization.F90 index 9ff99b583f..4f213d86d9 100644 --- a/src/user/Rossby_front_2d_initialization.F90 +++ b/src/user/Rossby_front_2d_initialization.F90 @@ -40,7 +40,7 @@ subroutine Rossby_front_initialize_thickness(h, G, GV, US, param_file, just_read type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2] + intent(out) :: h !< The thickness that is being initialized [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, intent(in) :: just_read !< If true, this call will only read @@ -83,7 +83,7 @@ subroutine Rossby_front_initialize_thickness(h, G, GV, US, param_file, just_read stretch = ( ( G%max_depth + eta ) / G%max_depth ) h0 = ( G%max_depth / real(nz) ) * stretch do k = 1, nz - h(i,j,k) = h0 * GV%Z_to_H + h(i,j,k) = h0 enddo enddo ; enddo @@ -94,7 +94,7 @@ subroutine Rossby_front_initialize_thickness(h, G, GV, US, param_file, just_read stretch = ( ( G%max_depth + eta ) / G%max_depth ) h0 = ( G%max_depth / real(nz) ) * stretch do k = 1, nz - h(i,j,k) = h0 * GV%Z_to_H + h(i,j,k) = h0 enddo enddo ; enddo @@ -114,7 +114,7 @@ subroutine Rossby_front_initialize_temperature_salinity(T, S, h, G, GV, US, & type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> ppt] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Thickness [Z ~> m] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Parameter file handle logical, intent(in) :: just_read !< If true, this call will @@ -125,7 +125,7 @@ subroutine Rossby_front_initialize_temperature_salinity(T, S, h, G, GV, US, & real :: S_ref ! Reference salinity within the surface layer [S ~> ppt] real :: T_range ! Range of temperatures over the vertical [C ~> degC] real :: zc ! Position of the middle of the cell [Z ~> m] - real :: zi ! Bottom interface position relative to the sea surface [H ~> m or kg m-2] + real :: zi ! Bottom interface position relative to the sea surface [Z ~> m] real :: dTdz ! Vertical temperature gradient [C Z-1 ~> degC m-1] character(len=40) :: verticalCoordinate @@ -149,8 +149,8 @@ subroutine Rossby_front_initialize_temperature_salinity(T, S, h, G, GV, US, & do j = G%jsc,G%jec ; do i = G%isc,G%iec zi = 0. do k = 1, nz - zi = zi - h(i,j,k) ! Bottom interface position - zc = GV%H_to_Z * (zi - 0.5*h(i,j,k)) ! Position of middle of cell + zi = zi - h(i,j,k) ! Bottom interface position + zc = zi - 0.5*h(i,j,k) ! Position of middle of cell zc = min( zc, -Hml(G, G%geoLatT(i,j)) ) ! Bound by depth of mixed layer T(i,j,k) = T_ref + dTdz * zc ! Linear temperature profile enddo diff --git a/src/user/SCM_CVMix_tests.F90 b/src/user/SCM_CVMix_tests.F90 index 8df8f90e3d..7b1b4b3946 100644 --- a/src/user/SCM_CVMix_tests.F90 +++ b/src/user/SCM_CVMix_tests.F90 @@ -57,7 +57,7 @@ subroutine SCM_CVMix_tests_TS_init(T, S, h, G, GV, US, param_file, just_read) type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> ppt] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Input parameter structure logical, intent(in) :: just_read !< If present and true, this call @@ -108,7 +108,7 @@ subroutine SCM_CVMix_tests_TS_init(T, S, h, G, GV, US, param_file, just_read) top = 0. ! Reference to surface bottom = 0. do k=1,nz - bottom = bottom - h(i,j,k)*GV%H_to_Z ! Interface below layer [Z ~> m] + bottom = bottom - h(i,j,k) ! Interface below layer [Z ~> m] zC = 0.5*( top + bottom ) ! Z of middle of layer [Z ~> m] DZ = min(0., zC + UpperLayerTempMLD) T(i,j,k) = max(LowerLayerMinTemp,LowerLayerTemp + LowerLayerdTdZ * DZ) diff --git a/src/user/adjustment_initialization.F90 b/src/user/adjustment_initialization.F90 index a958ebdebb..58389b7b5c 100644 --- a/src/user/adjustment_initialization.F90 +++ b/src/user/adjustment_initialization.F90 @@ -36,7 +36,7 @@ subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, intent(in) :: just_read !< If true, this call will only read @@ -71,7 +71,7 @@ subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke if (.not.just_read) & - call MOM_mesg("initialize_thickness_uniform: setting thickness") + call MOM_mesg("adjustment_initialize_thickness: setting thickness") ! Parameters used by main model initialization if (.not.just_read) call log_version(param_file, mdl, version, "") @@ -170,12 +170,12 @@ subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read do k=nz,1,-1 if (eta1D(k) > 0.) then eta1D(k) = max( eta1D(k+1) + min_thickness, 0. ) - h(i,j,k) = GV%Z_to_H * max( eta1D(k) - eta1D(k+1), min_thickness ) + h(i,j,k) = max( eta1D(k) - eta1D(k+1), min_thickness ) elseif (eta1D(k) <= (eta1D(k+1) + min_thickness)) then eta1D(k) = eta1D(k+1) + min_thickness - h(i,j,k) = GV%Z_to_H * min_thickness + h(i,j,k) = min_thickness else - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo enddo ; enddo @@ -187,7 +187,7 @@ subroutine adjustment_initialize_thickness ( h, G, GV, US, param_file, just_read enddo do j=js,je ; do i=is,ie do k=nz,1,-1 - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) enddo enddo ; enddo @@ -209,7 +209,7 @@ subroutine adjustment_initialize_temperature_salinity(T, S, h, depth_tot, G, GV, real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(out) :: S !< The salinity that is being initialized [S ~> ppt] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< The model thicknesses [H ~> m or kg m-2]. + intent(in) :: h !< The model thicknesses [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to @@ -275,7 +275,7 @@ subroutine adjustment_initialize_temperature_salinity(T, S, h, depth_tot, G, GV, do j=js,je ; do i=is,ie eta1d(nz+1) = -depth_tot(i,j) do k=nz,1,-1 - eta1d(k) = eta1d(k+1) + h(i,j,k)*GV%H_to_Z + eta1d(k) = eta1d(k+1) + h(i,j,k) enddo if (front_wave_length /= 0.) then y = ( 0.125 + G%geoLatT(i,j) / front_wave_length ) * ( 4. * acos(0.) ) @@ -296,7 +296,7 @@ subroutine adjustment_initialize_temperature_salinity(T, S, h, depth_tot, G, GV, x = 1. - min(1., x) T(i,j,k) = T_range * x enddo - ! x = GV%H_to_Z*sum(T(i,j,:)*h(i,j,:)) + ! x = sum(T(i,j,:)*h(i,j,:)) ! T(i,j,:) = (T(i,j,:) / x) * (G%max_depth*1.5/real(nz)) enddo ; enddo diff --git a/src/user/baroclinic_zone_initialization.F90 b/src/user/baroclinic_zone_initialization.F90 index 2ff4e1ec80..e2c6182231 100644 --- a/src/user/baroclinic_zone_initialization.F90 +++ b/src/user/baroclinic_zone_initialization.F90 @@ -86,7 +86,7 @@ subroutine baroclinic_zone_init_temperature_salinity(T, S, h, depth_tot, G, GV, real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & intent(out) :: S !< Salinity [S ~> ppt] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(in) :: h !< The model thicknesses [H ~> m or kg m-2] + intent(in) :: h !< The model thicknesses [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -135,8 +135,8 @@ subroutine baroclinic_zone_init_temperature_salinity(T, S, h, depth_tot, G, GV, fn = xs endif do k = nz, 1, -1 - zc = zi + 0.5*h(i,j,k)*GV%H_to_Z ! Position of middle of cell - zi = zi + h(i,j,k)*GV%H_to_Z ! Top interface position + zc = zi + 0.5*h(i,j,k) ! Position of middle of cell + zi = zi + h(i,j,k) ! Top interface position T(i,j,k) = T_ref + dTdz * zc & ! Linear temperature stratification + dTdx * x & ! Linear gradient + delta_T * fn ! Smooth fn of width L_zone diff --git a/src/user/benchmark_initialization.F90 b/src/user/benchmark_initialization.F90 index 3920b52729..333f53895e 100644 --- a/src/user/benchmark_initialization.F90 +++ b/src/user/benchmark_initialization.F90 @@ -84,7 +84,7 @@ subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, e type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -184,9 +184,9 @@ subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, e do k=1,nz ; e_pert(K) = 0.0 ; enddo - ! This sets the initial thickness (in [H ~> m or kg m-2]) of the layers. The thicknesses + ! This sets the initial thickness (in [Z ~> m]) of the layers. The thicknesses ! are set to insure that: - ! 1. each layer is at least GV%Angstrom_H thick, and + ! 1. each layer is at least GV%Angstrom_Z thick, and ! 2. the interfaces are where they should be based on the resting depths and ! interface height perturbations, as long at this doesn't interfere with 1. eta1D(nz+1) = -depth_tot(i,j) @@ -211,9 +211,9 @@ subroutine benchmark_initialize_thickness(h, depth_tot, G, GV, US, param_file, e if (eta1D(K) < eta1D(K+1) + GV%Angstrom_Z) & eta1D(K) = eta1D(K+1) + GV%Angstrom_Z - h(i,j,k) = max(GV%Z_to_H * (eta1D(K) - eta1D(K+1)), GV%Angstrom_H) + h(i,j,k) = max(eta1D(K) - eta1D(K+1), GV%Angstrom_Z) enddo - h(i,j,1) = max(GV%Z_to_H * (0.0 - eta1D(2)), GV%Angstrom_H) + h(i,j,1) = max(0.0 - eta1D(2), GV%Angstrom_Z) enddo ; enddo diff --git a/src/user/circle_obcs_initialization.F90 b/src/user/circle_obcs_initialization.F90 index 63c5c8a0d4..ab9ab385de 100644 --- a/src/user/circle_obcs_initialization.F90 +++ b/src/user/circle_obcs_initialization.F90 @@ -10,6 +10,7 @@ module circle_obcs_initialization use MOM_get_input, only : directories use MOM_grid, only : ocean_grid_type use MOM_tracer_registry, only : tracer_registry_type +use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : thermo_var_ptrs use MOM_verticalGrid, only : verticalGrid_type @@ -27,11 +28,12 @@ module circle_obcs_initialization contains !> This subroutine initializes layer thicknesses for the circle_obcs experiment. -subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, param_file, just_read) - type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. - type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. +subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, US, param_file, just_read) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -43,7 +45,7 @@ subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, param_file, jus ! negative because it is positive upward. real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface ! positive upward, in depth units [Z ~> m]. - real :: IC_amp ! The amplitude of the initial height displacement [H ~> m or kg m-2]. + real :: IC_amp ! The amplitude of the initial height displacement [Z ~> m]. real :: diskrad ! Radius of the elevated disk [km] or [degrees] or [m] real :: rad ! Distance from the center of the elevated disk [km] or [degrees] or [m] real :: lonC ! The x-position of a point [km] or [degrees] or [m] @@ -73,7 +75,7 @@ subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, param_file, jus call get_param(param_file, mdl, "DISK_IC_AMPLITUDE", IC_amp, & "Initial amplitude of interface height displacements "//& "in the circle_obcs test case.", & - units='m', default=5.0, scale=GV%m_to_H, do_not_log=just_read) + units='m', default=5.0, scale=US%m_to_Z, do_not_log=just_read) if (just_read) return ! All run-time parameters have been read, so return. @@ -88,9 +90,9 @@ subroutine circle_obcs_initialize_thickness(h, depth_tot, G, GV, param_file, jus eta1D(K) = e0(K) if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then eta1D(K) = eta1D(K+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + h(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (eta1D(K) - eta1D(K+1)) + h(i,j,k) = eta1D(K) - eta1D(K+1) endif enddo enddo ; enddo diff --git a/src/user/dense_water_initialization.F90 b/src/user/dense_water_initialization.F90 index 81aa4c2b3b..6feb2bdda6 100644 --- a/src/user/dense_water_initialization.F90 +++ b/src/user/dense_water_initialization.F90 @@ -105,7 +105,7 @@ subroutine dense_water_initialize_TS(G, GV, US, param_file, T, S, h, just_read) type(param_file_type), intent(in) :: param_file !< Parameter file structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Output temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Output salinity [S ~> ppt] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [Z ~> m] logical, intent(in) :: just_read !< If true, this call will !! only read parameters without changing T & S. ! Local variables @@ -137,7 +137,7 @@ subroutine dense_water_initialize_TS(G, GV, US, param_file, T, S, h, just_read) zi = 0. do k = 1,nz ! nondimensional middle of layer - zmid = zi + 0.5 * h(i,j,k) / (GV%Z_to_H * G%max_depth) + zmid = zi + 0.5 * h(i,j,k) / G%max_depth if (zmid < mld) then ! use reference salinity in the mixed layer @@ -147,7 +147,7 @@ subroutine dense_water_initialize_TS(G, GV, US, param_file, T, S, h, just_read) S(i,j,k) = S_ref + S_range * (zmid - mld) / (1.0 - mld) endif - zi = zi + h(i,j,k) / (GV%Z_to_H * G%max_depth) + zi = zi + h(i,j,k) / G%max_depth enddo enddo enddo diff --git a/src/user/dumbbell_initialization.F90 b/src/user/dumbbell_initialization.F90 index 0b65883eca..abd4f4f37e 100644 --- a/src/user/dumbbell_initialization.F90 +++ b/src/user/dumbbell_initialization.F90 @@ -96,7 +96,7 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -126,7 +126,7 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke if (.not.just_read) & - call MOM_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness") + call MOM_mesg("dumbbell_initialization.F90, dumbbell_initialize_thickness: setting thickness") if (.not.just_read) call log_version(param_file, mdl, version, "") call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, & @@ -174,7 +174,7 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, enddo endif do k=1,nz - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) enddo enddo enddo @@ -217,9 +217,9 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then eta1D(k) = eta1D(k+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + h(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo enddo ; enddo @@ -232,9 +232,9 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, eta1D(k) = -G%max_depth * real(k-1) / real(nz) if (eta1D(k) < (eta1D(k+1) + min_thickness)) then eta1D(k) = eta1D(k+1) + min_thickness - h(i,j,k) = GV%Z_to_H * min_thickness + h(i,j,k) = min_thickness else - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo enddo ; enddo @@ -242,7 +242,7 @@ subroutine dumbbell_initialize_thickness ( h, depth_tot, G, GV, US, param_file, case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie - h(i,j,:) = GV%Z_to_H * depth_tot(i,j) / real(nz) + h(i,j,:) = depth_tot(i,j) / real(nz) enddo ; enddo end select @@ -255,7 +255,7 @@ subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, US, param_ type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> ppt] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Parameter file structure logical, intent(in) :: just_read !< If true, this call will diff --git a/src/user/external_gwave_initialization.F90 b/src/user/external_gwave_initialization.F90 index 63cc89342a..437edc49b2 100644 --- a/src/user/external_gwave_initialization.F90 +++ b/src/user/external_gwave_initialization.F90 @@ -30,7 +30,7 @@ subroutine external_gwave_initialize_thickness(h, G, GV, US, param_file, just_re type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, intent(in) :: just_read !< If true, this call will only read @@ -73,7 +73,7 @@ subroutine external_gwave_initialize_thickness(h, G, GV, US, param_file, just_re enddo eta1D(nz+1) = -G%max_depth ! Force bottom interface to bottom do k=1,nz - h(i,j,k) = GV%Z_to_H * (eta1D(K) - eta1D(K+1)) + h(i,j,k) = eta1D(K) - eta1D(K+1) enddo enddo ; enddo diff --git a/src/user/lock_exchange_initialization.F90 b/src/user/lock_exchange_initialization.F90 index 3b41237c36..ab08d4068d 100644 --- a/src/user/lock_exchange_initialization.F90 +++ b/src/user/lock_exchange_initialization.F90 @@ -28,7 +28,7 @@ subroutine lock_exchange_initialize_thickness(h, G, GV, US, param_file, just_rea type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file !! to parse for model parameter values. logical, intent(in) :: just_read !< If true, this call will only read @@ -80,7 +80,7 @@ subroutine lock_exchange_initialize_thickness(h, G, GV, US, param_file, just_rea eta1D(K) = min( eta1D(K), eta1D(K-1) - GV%Angstrom_Z ) enddo do k=nz,1,-1 - h(i,j,k) = GV%Z_to_H * (eta1D(K) - eta1D(K+1)) + h(i,j,k) = eta1D(K) - eta1D(K+1) enddo enddo ; enddo diff --git a/src/user/seamount_initialization.F90 b/src/user/seamount_initialization.F90 index a1f978a784..d1971f25f9 100644 --- a/src/user/seamount_initialization.F90 +++ b/src/user/seamount_initialization.F90 @@ -84,7 +84,7 @@ subroutine seamount_initialize_thickness (h, depth_tot, G, GV, US, param_file, j type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open file @@ -105,7 +105,7 @@ subroutine seamount_initialize_thickness (h, depth_tot, G, GV, US, param_file, j is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke if (.not.just_read) & - call MOM_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness") + call MOM_mesg("seamount_initialization.F90, seamount_initialize_thickness: setting thickness") call get_param(param_file, mdl,"MIN_THICKNESS",min_thickness, & 'Minimum thickness for layer', & @@ -164,9 +164,9 @@ subroutine seamount_initialize_thickness (h, depth_tot, G, GV, US, param_file, j eta1D(k) = e0(k) if (eta1D(k) < (eta1D(k+1) + GV%Angstrom_Z)) then eta1D(k) = eta1D(k+1) + GV%Angstrom_Z - h(i,j,k) = GV%Angstrom_H + h(i,j,k) = GV%Angstrom_Z else - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo enddo ; enddo @@ -179,9 +179,9 @@ subroutine seamount_initialize_thickness (h, depth_tot, G, GV, US, param_file, j eta1D(k) = -G%max_depth * real(k-1) / real(nz) if (eta1D(k) < (eta1D(k+1) + min_thickness)) then eta1D(k) = eta1D(k+1) + min_thickness - h(i,j,k) = GV%Z_to_H * min_thickness + h(i,j,k) = min_thickness else - h(i,j,k) = GV%Z_to_H * (eta1D(k) - eta1D(k+1)) + h(i,j,k) = eta1D(k) - eta1D(k+1) endif enddo enddo ; enddo @@ -189,7 +189,7 @@ subroutine seamount_initialize_thickness (h, depth_tot, G, GV, US, param_file, j case ( REGRIDDING_SIGMA ) ! Initial thicknesses for sigma coordinates if (just_read) return ! All run-time parameters have been read, so return. do j=js,je ; do i=is,ie - h(i,j,:) = GV%Z_to_H * depth_tot(i,j) / real(nz) + h(i,j,:) = depth_tot(i,j) / real(nz) enddo ; enddo end select @@ -202,7 +202,7 @@ subroutine seamount_initialize_temperature_salinity(T, S, h, G, GV, US, param_fi type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> ppt] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< Parameter file structure logical, intent(in) :: just_read !< If true, this call will @@ -282,7 +282,7 @@ subroutine seamount_initialize_temperature_salinity(T, S, h, G, GV, US, param_fi do j=js,je ; do i=is,ie xi0 = 0.0 do k = 1,nz - xi1 = xi0 + GV%H_to_Z * h(i,j,k) / G%max_depth + xi1 = xi0 + h(i,j,k) / G%max_depth select case ( trim(density_profile) ) case ('linear') !S(i,j,k) = S_surf + S_range * 0.5 * (xi0 + xi1) diff --git a/src/user/sloshing_initialization.F90 b/src/user/sloshing_initialization.F90 index 357f247896..75e5889092 100644 --- a/src/user/sloshing_initialization.F90 +++ b/src/user/sloshing_initialization.F90 @@ -57,7 +57,7 @@ subroutine sloshing_initialize_thickness ( h, depth_tot, G, GV, US, param_file, type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure to parse for model parameter values. @@ -160,7 +160,7 @@ subroutine sloshing_initialize_thickness ( h, depth_tot, G, GV, US, param_file, ! 4. Define layers do k = 1,nz - h(i,j,k) = GV%Z_to_H * (z_inter(k) - z_inter(k+1)) + h(i,j,k) = z_inter(k) - z_inter(k+1) enddo enddo ; enddo @@ -179,7 +179,7 @@ subroutine sloshing_initialize_temperature_salinity ( T, S, h, G, GV, US, param_ type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< Potential temperature [C ~> degC]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< Salinity [S ~> ppt]. - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [Z ~> m]. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< A structure to parse !! for model parameter values. diff --git a/src/user/soliton_initialization.F90 b/src/user/soliton_initialization.F90 index b3b45da997..06a781ec94 100644 --- a/src/user/soliton_initialization.F90 +++ b/src/user/soliton_initialization.F90 @@ -32,7 +32,7 @@ subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US) type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thickness that is being initialized [Z ~> m] real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m] @@ -55,7 +55,7 @@ subroutine soliton_initialize_thickness(h, depth_tot, G, GV, US) y = G%geoLatT(i,j)-y0 val3 = exp(-val1*x) val4 = val2 * ( 2.0*val3 / (1.0 + (val3*val3)) )**2 - h(i,j,k) = GV%Z_to_H * (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + depth_tot(i,j)) + h(i,j,k) = (0.25*val4*(6.0*y*y + 3.0) * exp(-0.5*y*y) + depth_tot(i,j)) enddo enddo ; enddo @@ -63,12 +63,11 @@ end subroutine soliton_initialize_thickness !> Initialization of u and v in the equatorial Rossby soliton test -subroutine soliton_initialize_velocity(u, v, h, G, GV, US) +subroutine soliton_initialize_velocity(u, v, G, GV, US) type(ocean_grid_type), intent(in) :: G !< Grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(out) :: u !< i-component of velocity [L T-1 ~> m s-1] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(out) :: v !< j-component of velocity [L T-1 ~> m s-1] - real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Thickness [H ~> m or kg m-2] type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type ! Local variables diff --git a/src/user/user_initialization.F90 b/src/user/user_initialization.F90 index b9d16e548a..207f009c9c 100644 --- a/src/user/user_initialization.F90 +++ b/src/user/user_initialization.F90 @@ -76,12 +76,12 @@ subroutine USER_initialize_topography(D, G, param_file, max_depth, US) end subroutine USER_initialize_topography -!> initialize thicknesses. +!> Initialize thicknesses in depth units. These will be converted to thickness units later. subroutine USER_initialize_thickness(h, G, GV, param_file, just_read) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), & - intent(out) :: h !< The thicknesses being initialized [H ~> m or kg m-2]. + intent(out) :: h !< The thicknesses being initialized [Z ~> m] type(param_file_type), intent(in) :: param_file !< A structure indicating the open !! file to parse for model parameter values. logical, intent(in) :: just_read !< If true, this call will @@ -93,7 +93,8 @@ subroutine USER_initialize_thickness(h, G, GV, param_file, just_read) if (just_read) return ! All run-time parameters have been read, so return. - h(:,:,1) = 0.0 ! h should be set [H ~> m or kg m-2]. + h(:,:,1:GV%ke) = 0.0 ! h should be set in [Z ~> m]. It will be converted to thickness units + ! [H ~> m or kg m-2] once the temperatures and salinities are known. if (first_call) call write_user_log(param_file)