portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 1 | ! This file is part of MOM6, the Modular Ocean Model version 6. | |
| 2 | ! See the LICENSE file for licensing information. | |
| 3 | ! SPDX-License-Identifier: Apache-2.0 | |
| 4 | ||
| 5 | !> Initializes fixed aspects of the related to its vertical coordinate. | |
| 6 | module MOM_coord_initialization | |
| 7 | ||
| 8 | use MOM_debugging, only : chksum | |
| 9 | use MOM_EOS, only : calculate_density, EOS_type | |
| 10 | use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe | |
| 11 | use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint | |
| 12 | use MOM_file_parser, only : get_param, read_param, log_param, param_file_type, log_version | |
| 13 | use MOM_io, only : create_MOM_file, file_exists | |
| 14 | use MOM_io, only : MOM_netCDF_file, MOM_field | |
| 15 | use MOM_io, only : MOM_read_data, MOM_write_field, vardesc, var_desc, SINGLE_FILE | |
| 16 | use MOM_string_functions, only : slasher, uppercase | |
| 17 | use MOM_unit_scaling, only : unit_scale_type | |
| 18 | use MOM_variables, only : thermo_var_ptrs | |
| 19 | use MOM_verticalGrid, only : verticalGrid_type, setVerticalGridAxes | |
| 20 | use user_initialization, only : user_set_coord | |
| 21 | use BFB_initialization, only : BFB_set_coord | |
| 22 | ||
| 23 | implicit none ; private | |
| 24 | ||
| 25 | public MOM_initialize_coord, write_vertgrid_file | |
| 26 | ||
| 27 | ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional | |
| 28 | ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with | |
| 29 | ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units | |
| 30 | ! vary with the Boussinesq approximation, the Boussinesq variant is given first. | |
| 31 | ||
| 32 | character(len=40) :: mdl = "MOM_coord_initialization" !< This module's name. | |
| 33 | ||
| 34 | contains | |
| 35 | ||
| 36 | !> MOM_initialize_coord sets up time-invariant quantities related to MOM6's | |
| 37 | !! vertical coordinate. | |
| 38 | 1 | subroutine MOM_initialize_coord(GV, US, PF, tv, max_depth) |
| 39 | type(verticalGrid_type), intent(inout) :: GV !< Ocean vertical grid structure. | |
| 40 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 41 | type(param_file_type), intent(in) :: PF !< A structure indicating the open file | |
| 42 | !! to parse for model parameter values. | |
| 43 | type(thermo_var_ptrs), intent(inout) :: tv !< The thermodynamic variable structure. | |
| 44 | real, intent(in) :: max_depth !< The ocean's maximum depth [Z ~> m]. | |
| 45 | ! Local | |
| 46 | character(len=200) :: config | |
| 47 | logical :: debug | |
| 48 | ! This include declares and sets the variable "version". | |
| 49 | # include "version_variable.h" | |
| 50 | integer :: nz | |
| 51 | ||
| 52 | 1 | nz = GV%ke |
| 53 | ||
| 54 | 1 | call callTree_enter("MOM_initialize_coord(), MOM_coord_initialization.F90") |
| 55 | 1 | call log_version(PF, mdl, version, "") |
| 56 | 1 | call get_param(PF, mdl, "DEBUG", debug, default=.false.) |
| 57 | ||
| 58 | ! Set-up the layer densities, GV%Rlay, and reduced gravities, GV%g_prime. | |
| 59 | call get_param(PF, mdl, "COORD_CONFIG", config, & | |
| 60 | "This specifies how layers are to be defined: \n"//& | |
| 61 | " \t ALE or none - used to avoid defining layers in ALE mode \n"//& | |
| 62 | " \t file - read coordinate information from the file \n"//& | |
| 63 | " \t\t specified by (COORD_FILE).\n"//& | |
| 64 | " \t BFB - Custom coords for buoyancy-forced basin case \n"//& | |
| 65 | " \t\t based on SST_S, T_BOT and DRHO_DT.\n"//& | |
| 66 | " \t linear - linear based on interfaces not layers \n"//& | |
| 67 | " \t layer_ref - linear based on layer densities \n"//& | |
| 68 | " \t ts_ref - use reference temperature and salinity \n"//& | |
| 69 | " \t ts_range - use range of temperature and salinity \n"//& | |
| 70 | " \t\t (T_REF and S_REF) to determine surface density \n"//& | |
| 71 | " \t\t and GINT calculate internal densities. \n"//& | |
| 72 | " \t gprime - use reference density (RHO_0) for surface \n"//& | |
| 73 | " \t\t density and GINT calculate internal densities. \n"//& | |
| 74 | " \t ts_profile - use temperature and salinity profiles \n"//& | |
| 75 | " \t\t (read from COORD_FILE) to set layer densities. \n"//& | |
| 76 | " \t USER - call a user modified routine.", & | |
| 77 | 1 | default="none") |
| 78 | 2 | select case ( trim(config) ) |
| 79 | case ("gprime") | |
| 80 | 0 | call set_coord_from_gprime(GV%Rlay, GV%g_prime, GV, US, PF) |
| 81 | case ("layer_ref") | |
| 82 | 0 | call set_coord_from_layer_density(GV%Rlay, GV%g_prime, GV, US, PF) |
| 83 | case ("linear") | |
| 84 | 0 | call set_coord_linear(GV%Rlay, GV%g_prime, GV, US, PF) |
| 85 | case ("ts_ref") | |
| 86 | 0 | call set_coord_from_TS_ref(GV%Rlay, GV%g_prime, GV, US, PF, tv%eqn_of_state, tv%P_Ref) |
| 87 | case ("ts_profile") | |
| 88 | 0 | call set_coord_from_TS_profile(GV%Rlay, GV%g_prime, GV, US, PF, tv%eqn_of_state, tv%P_Ref) |
| 89 | case ("ts_range") | |
| 90 | 0 | call set_coord_from_TS_range(GV%Rlay, GV%g_prime, GV, US, PF, tv%eqn_of_state, tv%P_Ref) |
| 91 | case ("file") | |
| 92 | 0 | call set_coord_from_file(GV%Rlay, GV%g_prime, GV, US, PF) |
| 93 | case ("USER") | |
| 94 | 0 | call user_set_coord(GV%Rlay, GV%g_prime, GV, US, PF) |
| 95 | case ("BFB") | |
| 96 | 0 | call BFB_set_coord(GV%Rlay, GV%g_prime, GV, US, PF) |
| 97 | case ("none", "ALE") | |
| 98 | 1 | call set_coord_to_none(GV%Rlay, GV%g_prime, GV, US, PF) |
| 99 | case default ; call MOM_error(FATAL,"MOM_initialize_coord: "// & | |
| 100 | 2 | "Unrecognized coordinate setup"//trim(config)) |
| 101 | end select | |
| 102 | ! There are nz+1 values of g_prime because it is an interface field, but the value at the bottom | |
| 103 | ! should not matter. This is here just to avoid having an uninitialized value in some output. | |
| 104 | 1 | GV%g_prime(nz+1) = 10.0*GV%g_Earth |
| 105 | ||
| 106 | 1 | if (debug) call chksum(US%R_to_kg_m3*GV%Rlay(:), "MOM_initialize_coord: Rlay ", 1, nz) |
| 107 | 1 | if (debug) call chksum(US%m_to_Z*US%L_to_m**2*US%s_to_T**2*GV%g_prime(:), "MOM_initialize_coord: g_prime ", 1, nz) |
| 108 | 1 | call setVerticalGridAxes( GV%Rlay, GV, scale=US%R_to_kg_m3 ) |
| 109 | ||
| 110 | ! Copy the maximum depth across from the input argument | |
| 111 | 1 | GV%max_depth = max_depth |
| 112 | ||
| 113 | 1 | call callTree_leave('MOM_initialize_coord()') |
| 114 | ||
| 115 | 1 | end subroutine MOM_initialize_coord |
| 116 | ||
| 117 | ! The set_coord routines deal with initializing aspects of the vertical grid. | |
| 118 | ||
| 119 | !> Sets the layer densities (Rlay) and the interface reduced gravities (g). | |
| 120 | 0 | subroutine set_coord_from_gprime(Rlay, g_prime, GV, US, param_file) |
| 121 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 122 | real, dimension(GV%ke), intent(out) :: Rlay !< The layers' target coordinate values | |
| 123 | !! (potential density) [R ~> kg m-3]. | |
| 124 | real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity across the interfaces | |
| 125 | !! [L2 Z-1 T-2 ~> m s-2]. | |
| 126 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 127 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 128 | ! Local variables | |
| 129 | real :: g_int ! Reduced gravities across the internal interfaces [L2 Z-1 T-2 ~> m s-2]. | |
| 130 | real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. | |
| 131 | real :: Rlay_Ref ! The target density of the surface layer [R ~> kg m-3]. | |
| 132 | character(len=40) :: mdl = "set_coord_from_gprime" ! This subroutine's name. | |
| 133 | integer :: k, nz | |
| 134 | 0 | nz = GV%ke |
| 135 | ||
| 136 | 0 | call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90") |
| 137 | ||
| 138 | call get_param(param_file, mdl, "GFS" , g_fs, & | |
| 139 | "The reduced gravity at the free surface.", units="m s-2", & | |
| 140 | 0 | default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m) |
| 141 | call get_param(param_file, mdl, "GINT", g_int, & | |
| 142 | "The reduced gravity across internal interfaces.", & | |
| 143 | 0 | units="m s-2", fail_if_missing=.true., scale=US%m_s_to_L_T**2*US%Z_to_m) |
| 144 | call get_param(param_file, mdl, "LIGHTEST_DENSITY", Rlay_Ref, & | |
| 145 | "The reference potential density used for layer 1.", & | |
| 146 | 0 | units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R) |
| 147 | ||
| 148 | 0 | g_prime(1) = g_fs |
| 149 | 0 | do k=2,nz ; g_prime(k) = g_int ; enddo |
| 150 | 0 | Rlay(1) = Rlay_Ref |
| 151 | 0 | if (GV%Boussinesq .or. GV%semi_Boussinesq) then |
| 152 | 0 | do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo |
| 153 | else | |
| 154 | 0 | do k=2,nz |
| 155 | 0 | Rlay(k) = Rlay(k-1) * ((GV%g_Earth + 0.5*g_prime(k)) / (GV%g_Earth - 0.5*g_prime(k))) |
| 156 | enddo | |
| 157 | endif | |
| 158 | ||
| 159 | 0 | call callTree_leave(trim(mdl)//'()') |
| 160 | ||
| 161 | 0 | end subroutine set_coord_from_gprime |
| 162 | ||
| 163 | !> Sets the layer densities (Rlay) and the interface reduced gravities (g). | |
| 164 | 0 | subroutine set_coord_from_layer_density(Rlay, g_prime, GV, US, param_file) |
| 165 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 166 | real, dimension(GV%ke), intent(out) :: Rlay !< The layers' target coordinate values | |
| 167 | !! (potential density) [R ~> kg m-3]. | |
| 168 | real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity across the interfaces | |
| 169 | !! [L2 Z-1 T-2 ~> m s-2]. | |
| 170 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 171 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 172 | ||
| 173 | ! Local variables | |
| 174 | real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. | |
| 175 | real :: Rlay_Ref! The surface layer's target density [R ~> kg m-3]. | |
| 176 | real :: RLay_range ! The range of densities [R ~> kg m-3]. | |
| 177 | character(len=40) :: mdl = "set_coord_from_layer_density" ! This subroutine's name. | |
| 178 | integer :: k, nz | |
| 179 | 0 | nz = GV%ke |
| 180 | ||
| 181 | 0 | call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90") |
| 182 | ||
| 183 | call get_param(param_file, mdl, "GFS", g_fs, & | |
| 184 | "The reduced gravity at the free surface.", units="m s-2", & | |
| 185 | 0 | default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m) |
| 186 | call get_param(param_file, mdl, "LIGHTEST_DENSITY", Rlay_Ref, & | |
| 187 | "The reference potential density used for layer 1.", & | |
| 188 | 0 | units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R) |
| 189 | call get_param(param_file, mdl, "DENSITY_RANGE", Rlay_range, & | |
| 190 | "The range of reference potential densities in the layers.", & | |
| 191 | 0 | units="kg m-3", default=2.0, scale=US%kg_m3_to_R) |
| 192 | ||
| 193 | 0 | Rlay(1) = Rlay_Ref |
| 194 | 0 | do k=2,nz |
| 195 | 0 | Rlay(k) = Rlay(k-1) + RLay_range/(real(nz-1)) |
| 196 | enddo | |
| 197 | ! These statements set the interface reduced gravities. ! | |
| 198 | 0 | g_prime(1) = g_fs |
| 199 | 0 | if (GV%Boussinesq .or. GV%semi_Boussinesq) then |
| 200 | 0 | do k=2,nz |
| 201 | 0 | g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) |
| 202 | enddo | |
| 203 | else | |
| 204 | 0 | do k=2,nz |
| 205 | 0 | g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1)) |
| 206 | enddo | |
| 207 | endif | |
| 208 | ||
| 209 | 0 | call callTree_leave(trim(mdl)//'()') |
| 210 | 0 | end subroutine set_coord_from_layer_density |
| 211 | ||
| 212 | !> Sets the layer densities (Rlay) and the interface reduced gravities (g) from a profile of g'. | |
| 213 | 0 | subroutine set_coord_from_TS_ref(Rlay, g_prime, GV, US, param_file, eqn_of_state, P_Ref) |
| 214 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. | |
| 215 | real, dimension(GV%ke), intent(out) :: Rlay !< The layers' target coordinate values | |
| 216 | !! (potential density) [R ~> kg m-3]. | |
| 217 | real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity across the interfaces | |
| 218 | !! [L2 Z-1 T-2 ~> m s-2]. | |
| 219 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 220 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 221 | type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure | |
| 222 | real, intent(in) :: P_Ref !< The coordinate-density reference pressure | |
| 223 | !! [R L2 T-2 ~> Pa]. | |
| 224 | ||
| 225 | ! Local variables | |
| 226 | real :: T_ref ! Reference temperature [C ~> degC] | |
| 227 | real :: S_ref ! Reference salinity [S ~> ppt] | |
| 228 | real :: g_int ! Reduced gravities across the internal interfaces [L2 Z-1 T-2 ~> m s-2]. | |
| 229 | real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. | |
| 230 | character(len=40) :: mdl = "set_coord_from_TS_ref" ! This subroutine's name. | |
| 231 | integer :: k, nz | |
| 232 | 0 | nz = GV%ke |
| 233 | ||
| 234 | 0 | call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90") |
| 235 | ||
| 236 | call get_param(param_file, mdl, "T_REF", T_ref, & | |
| 237 | "The initial temperature of the lightest layer.", & | |
| 238 | 0 | units="degC", scale=US%degC_to_C, fail_if_missing=.true.) |
| 239 | call get_param(param_file, mdl, "S_REF", S_ref, & | |
| 240 | 0 | "The initial salinities.", units="ppt", default=35.0, scale=US%ppt_to_S) |
| 241 | call get_param(param_file, mdl, "GFS", g_fs, & | |
| 242 | "The reduced gravity at the free surface.", units="m s-2", & | |
| 243 | 0 | default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m) |
| 244 | call get_param(param_file, mdl, "GINT", g_int, & | |
| 245 | "The reduced gravity across internal interfaces.", & | |
| 246 | 0 | units="m s-2", fail_if_missing=.true., scale=US%m_s_to_L_T**2*US%Z_to_m) |
| 247 | ||
| 248 | ! These statements set the interface reduced gravities. ! | |
| 249 | 0 | g_prime(1) = g_fs |
| 250 | 0 | do k=2,nz ; g_prime(k) = g_int ; enddo |
| 251 | ||
| 252 | ! The uppermost layer's density is set here. Subsequent layers' ! | |
| 253 | ! densities are determined from this value and the g values. ! | |
| 254 | ! T0 = 28.228 ; S0 = 34.5848 ; Pref = P_Ref | |
| 255 | 0 | call calculate_density(T_ref, S_ref, P_ref, Rlay(1), eqn_of_state) |
| 256 | ||
| 257 | ! These statements set the layer densities. ! | |
| 258 | 0 | if (GV%Boussinesq .or. GV%semi_Boussinesq) then |
| 259 | 0 | do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo |
| 260 | else | |
| 261 | 0 | do k=2,nz |
| 262 | 0 | Rlay(k) = Rlay(k-1) * ((GV%g_Earth + 0.5*g_prime(k)) / (GV%g_Earth - 0.5*g_prime(k))) |
| 263 | enddo | |
| 264 | endif | |
| 265 | ||
| 266 | 0 | call callTree_leave(trim(mdl)//'()') |
| 267 | 0 | end subroutine set_coord_from_TS_ref |
| 268 | ||
| 269 | !> Sets the layer densities (Rlay) and the interface reduced gravities (g) from a T-S profile. | |
| 270 | 0 | subroutine set_coord_from_TS_profile(Rlay, g_prime, GV, US, param_file, eqn_of_state, P_Ref) |
| 271 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 272 | real, dimension(GV%ke), intent(out) :: Rlay !< Layer potential density [R ~> kg m-3]. | |
| 273 | real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each | |
| 274 | !! interface [L2 Z-1 T-2 ~> m s-2]. | |
| 275 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 276 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 277 | type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure | |
| 278 | real, intent(in) :: P_Ref !< The coordinate-density reference pressure | |
| 279 | !! [R L2 T-2 ~> Pa]. | |
| 280 | ||
| 281 | ! Local variables | |
| 282 | 0 | real, dimension(GV%ke) :: T0 ! A profile of temperatures [C ~> degC] |
| 283 | 0 | real, dimension(GV%ke) :: S0 ! A profile of salinities [S ~> ppt] |
| 284 | 0 | real, dimension(GV%ke) :: Pref ! A array of reference pressures [R L2 T-2 ~> Pa] |
| 285 | real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. | |
| 286 | integer :: k, nz | |
| 287 | character(len=40) :: mdl = "set_coord_from_TS_profile" ! This subroutine's name. | |
| 288 | character(len=200) :: filename, coord_file, inputdir ! Strings for file/path | |
| 289 | character(len=64) :: temp_var, salt_var ! Temperature and salinity names in files | |
| 290 | ||
| 291 | 0 | nz = GV%ke |
| 292 | ||
| 293 | 0 | call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90") |
| 294 | ||
| 295 | call get_param(param_file, mdl, "GFS", g_fs, & | |
| 296 | "The reduced gravity at the free surface.", units="m s-2", & | |
| 297 | 0 | default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m) |
| 298 | call get_param(param_file, mdl, "COORD_FILE", coord_file, & | |
| 299 | "The file from which the coordinate temperatures and salinities are read.", & | |
| 300 | 0 | fail_if_missing=.true.) |
| 301 | call get_param(param_file, mdl, "TEMP_COORD_VAR", temp_var, & | |
| 302 | "The coordinate reference profile variable name for potential temperature.", & | |
| 303 | 0 | default="PTEMP") |
| 304 | call get_param(param_file, mdl, "SALT_COORD_VAR", salt_var, & | |
| 305 | "The coordinate reference profile variable name for salinity.", & | |
| 306 | 0 | default="SALT") |
| 307 | ||
| 308 | 0 | call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") |
| 309 | 0 | filename = trim(slasher(inputdir))//trim(coord_file) |
| 310 | 0 | call log_param(param_file, mdl, "INPUTDIR/COORD_FILE", filename) |
| 311 | ||
| 312 | 0 | call MOM_read_data(filename, temp_var, T0(:), scale=US%degC_to_C) |
| 313 | 0 | call MOM_read_data(filename, salt_var, S0(:), scale=US%ppt_to_S) |
| 314 | ||
| 315 | 0 | if (.not.file_exists(filename)) call MOM_error(FATAL, & |
| 316 | 0 | " set_coord_from_TS_profile: Unable to open " //trim(filename)) |
| 317 | ! These statements set the interface reduced gravities. ! | |
| 318 | 0 | g_prime(1) = g_fs |
| 319 | 0 | do k=1,nz ; Pref(k) = P_Ref ; enddo |
| 320 | 0 | call calculate_density(T0, S0, Pref, Rlay, eqn_of_state, (/1,nz/) ) |
| 321 | 0 | if (GV%Boussinesq .or. GV%semi_Boussinesq) then |
| 322 | 0 | do k=2,nz |
| 323 | 0 | g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) |
| 324 | enddo | |
| 325 | else | |
| 326 | 0 | do k=2,nz |
| 327 | 0 | g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1)) |
| 328 | enddo | |
| 329 | endif | |
| 330 | ||
| 331 | 0 | call callTree_leave(trim(mdl)//'()') |
| 332 | 0 | end subroutine set_coord_from_TS_profile |
| 333 | ||
| 334 | !> Sets the layer densities (Rlay) and the interface reduced gravities (g) from a linear T-S profile. | |
| 335 | 0 | subroutine set_coord_from_TS_range(Rlay, g_prime, GV, US, param_file, eqn_of_state, P_Ref) |
| 336 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 337 | real, dimension(GV%ke), intent(out) :: Rlay !< Layer potential density [R ~> kg m-3]. | |
| 338 | real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each | |
| 339 | !! interface [L2 Z-1 T-2 ~> m s-2]. | |
| 340 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 341 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 342 | type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure | |
| 343 | real, intent(in) :: P_Ref !< The coordinate-density reference pressure | |
| 344 | !! [R L2 T-2 ~> Pa]. | |
| 345 | ||
| 346 | ! Local variables | |
| 347 | 0 | real, dimension(GV%ke) :: T0 ! A profile of temperatures [C ~> degC] |
| 348 | 0 | real, dimension(GV%ke) :: S0 ! A profile of salinities [S ~> ppt] |
| 349 | 0 | real, dimension(GV%ke) :: Pref ! A array of reference pressures [R L2 T-2 ~> Pa] |
| 350 | real :: S_Ref ! Default salinity range parameters [S ~> ppt]. | |
| 351 | real :: T_Ref ! Default temperature range parameters [C ~> degC]. | |
| 352 | real :: S_Light, S_Dense ! Salinity range parameters [S ~> ppt]. | |
| 353 | real :: T_Light, T_Dense ! Temperature range parameters [C ~> degC]. | |
| 354 | real :: res_rat ! The ratio of density space resolution in the denser part | |
| 355 | ! of the range to that in the lighter part of the range. | |
| 356 | ! Setting this greater than 1 increases the resolution for | |
| 357 | ! the denser water [nondim]. | |
| 358 | real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. | |
| 359 | real :: a1, frac_dense, k_frac ! Nondimensional temporary variables [nondim] | |
| 360 | integer :: k, nz, k_light | |
| 361 | character(len=40) :: mdl = "set_coord_from_TS_range" ! This subroutine's name. | |
| 362 | ||
| 363 | 0 | nz = GV%ke |
| 364 | ||
| 365 | 0 | call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90") |
| 366 | ||
| 367 | call get_param(param_file, mdl, "T_REF", T_Ref, & | |
| 368 | "The default initial temperatures.", & | |
| 369 | 0 | units="degC", default=10.0, scale=US%degC_to_C) |
| 370 | call get_param(param_file, mdl, "TS_RANGE_T_LIGHT", T_Light, & | |
| 371 | "The initial temperature of the lightest layer when "//& | |
| 372 | "COORD_CONFIG is set to ts_range.", & | |
| 373 | 0 | units="degC", default=US%C_to_degC*T_Ref, scale=US%degC_to_C) |
| 374 | call get_param(param_file, mdl, "TS_RANGE_T_DENSE", T_Dense, & | |
| 375 | "The initial temperature of the densest layer when "//& | |
| 376 | "COORD_CONFIG is set to ts_range.", & | |
| 377 | 0 | units="degC", default=US%C_to_degC*T_Ref, scale=US%degC_to_C) |
| 378 | ||
| 379 | call get_param(param_file, mdl, "S_REF", S_Ref, & | |
| 380 | "The default initial salinities.", & | |
| 381 | 0 | units="ppt", default=35.0, scale=US%ppt_to_S) |
| 382 | call get_param(param_file, mdl, "TS_RANGE_S_LIGHT", S_Light, & | |
| 383 | "The initial lightest salinities when COORD_CONFIG is set to ts_range.", & | |
| 384 | 0 | units="ppt", default=US%S_to_ppt*S_Ref, scale=US%ppt_to_S) |
| 385 | call get_param(param_file, mdl, "TS_RANGE_S_DENSE", S_Dense, & | |
| 386 | "The initial densest salinities when COORD_CONFIG is set to ts_range.", & | |
| 387 | 0 | units="ppt", default=US%S_to_ppt*S_Ref, scale=US%ppt_to_S) |
| 388 | ||
| 389 | call get_param(param_file, mdl, "TS_RANGE_RESOLN_RATIO", res_rat, & | |
| 390 | "The ratio of density space resolution in the densest "//& | |
| 391 | "part of the range to that in the lightest part of the "//& | |
| 392 | "range when COORD_CONFIG is set to ts_range. Values "//& | |
| 393 | "greater than 1 increase the resolution of the denser water.",& | |
| 394 | 0 | default=1.0, units="nondim") |
| 395 | ||
| 396 | call get_param(param_file, mdl, "GFS", g_fs, & | |
| 397 | "The reduced gravity at the free surface.", units="m s-2", & | |
| 398 | 0 | default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m) |
| 399 | ||
| 400 | 0 | if ((GV%nk_rho_varies > 0) .and. (nz < GV%nk_rho_varies+2)) & |
| 401 | 0 | call MOM_error(FATAL, "set_coord_from_TS_range requires that NZ >= NKML+NKBL+2.") |
| 402 | ||
| 403 | 0 | k_light = GV%nk_rho_varies + 1 |
| 404 | ||
| 405 | ! Set T0(k) to range from T_LIGHT to T_DENSE, and similarly for S0(k). | |
| 406 | 0 | T0(k_light) = T_Light ; S0(k_light) = S_Light |
| 407 | 0 | a1 = 2.0 * res_rat / (1.0 + res_rat) |
| 408 | 0 | do k=k_light+1,nz |
| 409 | 0 | k_frac = real(k-k_light)/real(nz-k_light) |
| 410 | 0 | frac_dense = a1 * k_frac + (1.0 - a1) * k_frac**2 |
| 411 | 0 | T0(k) = frac_dense * (T_Dense - T_Light) + T_Light |
| 412 | 0 | S0(k) = frac_dense * (S_Dense - S_Light) + S_Light |
| 413 | enddo | |
| 414 | ||
| 415 | 0 | g_prime(1) = g_fs |
| 416 | 0 | do k=1,nz ; Pref(k) = P_Ref ; enddo |
| 417 | 0 | call calculate_density(T0, S0, Pref, Rlay, eqn_of_state, (/k_light,nz/) ) |
| 418 | ! Extrapolate target densities for the variable density mixed and buffer layers. | |
| 419 | 0 | do k=k_light-1,1,-1 |
| 420 | 0 | Rlay(k) = 2.0*Rlay(k+1) - Rlay(k+2) |
| 421 | enddo | |
| 422 | 0 | if (GV%Boussinesq .or. GV%semi_Boussinesq) then |
| 423 | 0 | do k=2,nz |
| 424 | 0 | g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) |
| 425 | enddo | |
| 426 | else | |
| 427 | 0 | do k=2,nz |
| 428 | 0 | g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1)) |
| 429 | enddo | |
| 430 | endif | |
| 431 | ||
| 432 | 0 | call callTree_leave(trim(mdl)//'()') |
| 433 | 0 | end subroutine set_coord_from_TS_range |
| 434 | ||
| 435 | ! Sets the layer densities (Rlay) and the interface reduced gravities (g) from data in file. | |
| 436 | 0 | subroutine set_coord_from_file(Rlay, g_prime, GV, US, param_file) |
| 437 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 438 | real, dimension(GV%ke), intent(out) :: Rlay !< Layer potential density [R ~> kg m-3]. | |
| 439 | real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each | |
| 440 | !! interface [L2 Z-1 T-2 ~> m s-2]. | |
| 441 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 442 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 443 | ||
| 444 | ! Local variables | |
| 445 | real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. | |
| 446 | integer :: k, nz | |
| 447 | character(len=40) :: mdl = "set_coord_from_file" ! This subroutine's name. | |
| 448 | character(len=40) :: coord_var | |
| 449 | character(len=200) :: filename,coord_file,inputdir ! Strings for file/path | |
| 450 | 0 | nz = GV%ke |
| 451 | ||
| 452 | 0 | call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90") |
| 453 | ||
| 454 | call get_param(param_file, mdl, "GFS", g_fs, & | |
| 455 | "The reduced gravity at the free surface.", units="m s-2", & | |
| 456 | 0 | default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m) |
| 457 | 0 | call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".") |
| 458 | 0 | inputdir = slasher(inputdir) |
| 459 | call get_param(param_file, mdl, "COORD_FILE", coord_file, & | |
| 460 | "The file from which the coordinate densities are read.", & | |
| 461 | 0 | fail_if_missing=.true.) |
| 462 | call get_param(param_file, mdl, "COORD_VAR", coord_var, & | |
| 463 | "The variable in COORD_FILE that is to be used for the "//& | |
| 464 | 0 | "coordinate densities.", default="Layer") |
| 465 | 0 | filename = trim(inputdir)//trim(coord_file) |
| 466 | 0 | call log_param(param_file, mdl, "INPUTDIR/COORD_FILE", filename) |
| 467 | 0 | if (.not.file_exists(filename)) call MOM_error(FATAL, & |
| 468 | 0 | " set_coord_from_file: Unable to open "//trim(filename)) |
| 469 | ||
| 470 | 0 | call MOM_read_data(filename, coord_var, Rlay, scale=US%kg_m3_to_R) |
| 471 | 0 | g_prime(1) = g_fs |
| 472 | 0 | if (GV%Boussinesq .or. GV%semi_Boussinesq) then |
| 473 | 0 | do k=2,nz |
| 474 | 0 | g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) |
| 475 | enddo | |
| 476 | else | |
| 477 | 0 | do k=2,nz |
| 478 | 0 | g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1)) |
| 479 | enddo | |
| 480 | endif | |
| 481 | 0 | do k=1,nz ; if (g_prime(k) <= 0.0) then |
| 482 | call MOM_error(FATAL, "MOM_initialization set_coord_from_file: "//& | |
| 483 | "Zero or negative g_primes read from variable "//"Layer"//" in file "//& | |
| 484 | 0 | trim(filename)) |
| 485 | endif ; enddo | |
| 486 | ||
| 487 | 0 | call callTree_leave(trim(mdl)//'()') |
| 488 | 0 | end subroutine set_coord_from_file |
| 489 | ||
| 490 | !> Sets the layer densities (Rlay) and the interface | |
| 491 | !! reduced gravities (g) according to a linear profile starting at a | |
| 492 | !! reference surface layer density and spanning a range of densities | |
| 493 | !! to the bottom defined by the parameter RLAY_RANGE | |
| 494 | !! (defaulting to 2.0 if not defined) | |
| 495 | 0 | subroutine set_coord_linear(Rlay, g_prime, GV, US, param_file) |
| 496 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 497 | real, dimension(GV%ke), intent(out) :: Rlay !< Layer potential density [R ~> kg m-3]. | |
| 498 | real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each | |
| 499 | !! interface [L2 Z-1 T-2 ~> m s-2]. | |
| 500 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 501 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 502 | ||
| 503 | ! Local variables | |
| 504 | character(len=40) :: mdl = "set_coord_linear" ! This subroutine | |
| 505 | real :: Rlay_ref, Rlay_range ! A reference density and its range [R ~> kg m-3] | |
| 506 | real :: g_fs ! The reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2] | |
| 507 | integer :: k, nz | |
| 508 | 0 | nz = GV%ke |
| 509 | ||
| 510 | 0 | call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90") |
| 511 | ||
| 512 | call get_param(param_file, mdl, "LIGHTEST_DENSITY", Rlay_Ref, & | |
| 513 | "The reference potential density used for the surface interface.", & | |
| 514 | 0 | units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R) |
| 515 | call get_param(param_file, mdl, "DENSITY_RANGE", Rlay_range, & | |
| 516 | "The range of reference potential densities across all interfaces.", & | |
| 517 | 0 | units="kg m-3", default=2.0, scale=US%kg_m3_to_R) |
| 518 | call get_param(param_file, mdl, "GFS", g_fs, & | |
| 519 | "The reduced gravity at the free surface.", units="m s-2", & | |
| 520 | 0 | default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m) |
| 521 | ||
| 522 | ! This following sets the target layer densities such that the | |
| 523 | ! surface interface has density Rlay_ref and the bottom | |
| 524 | ! is Rlay_range larger | |
| 525 | 0 | do k=1,nz |
| 526 | 0 | Rlay(k) = Rlay_Ref + RLay_range*((real(k)-0.5)/real(nz)) |
| 527 | enddo | |
| 528 | ! These statements set the interface reduced gravities. | |
| 529 | 0 | g_prime(1) = g_fs |
| 530 | 0 | if (GV%Boussinesq .or. GV%semi_Boussinesq) then |
| 531 | 0 | do k=2,nz |
| 532 | 0 | g_prime(k) = (GV%g_Earth/GV%Rho0) * (Rlay(k) - Rlay(k-1)) |
| 533 | enddo | |
| 534 | else | |
| 535 | 0 | do k=2,nz |
| 536 | 0 | g_prime(k) = 2.0*GV%g_Earth * (Rlay(k) - Rlay(k-1)) / (Rlay(k) + Rlay(k-1)) |
| 537 | enddo | |
| 538 | endif | |
| 539 | ||
| 540 | 0 | call callTree_leave(trim(mdl)//'()') |
| 541 | 0 | end subroutine set_coord_linear |
| 542 | ||
| 543 | !> Sets Rlay to Rho0 and g_prime to zero except for the free surface. | |
| 544 | !! This is for use only in ALE mode where Rlay should not be used and g_prime(1) alone | |
| 545 | !! might be used. | |
| 546 | 1 | subroutine set_coord_to_none(Rlay, g_prime, GV, US, param_file) |
| 547 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 548 | real, dimension(GV%ke), intent(out) :: Rlay !< Layer potential density [R ~> kg m-3]. | |
| 549 | real, dimension(GV%ke+1), intent(out) :: g_prime !< The reduced gravity at each | |
| 550 | !! interface [L2 Z-1 T-2 ~> m s-2]. | |
| 551 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 552 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 553 | ! Local variables | |
| 554 | real :: g_fs ! Reduced gravity across the free surface [L2 Z-1 T-2 ~> m s-2]. | |
| 555 | real :: Rlay_Ref ! The target density of the surface layer [R ~> kg m-3]. | |
| 556 | character(len=40) :: mdl = "set_coord_to_none" ! This subroutine's name. | |
| 557 | integer :: k, nz | |
| 558 | 1 | nz = GV%ke |
| 559 | ||
| 560 | 1 | call callTree_enter(trim(mdl)//"(), MOM_coord_initialization.F90") |
| 561 | ||
| 562 | call get_param(param_file, mdl, "GFS" , g_fs, & | |
| 563 | "The reduced gravity at the free surface.", units="m s-2", & | |
| 564 | 1 | default=GV%g_Earth*US%L_T_to_m_s**2*US%m_to_Z, scale=US%m_s_to_L_T**2*US%Z_to_m) |
| 565 | call get_param(param_file, mdl, "LIGHTEST_DENSITY", Rlay_Ref, & | |
| 566 | "The reference potential density used for layer 1.", & | |
| 567 | 1 | units="kg m-3", default=US%R_to_kg_m3*GV%Rho0, scale=US%kg_m3_to_R) |
| 568 | ||
| 569 | 1 | g_prime(1) = g_fs |
| 570 | 75 | do k=2,nz ; g_prime(k) = 0. ; enddo |
| 571 | 1 | Rlay(1) = Rlay_Ref |
| 572 | 1 | if (GV%Boussinesq .or. GV%semi_Boussinesq) then |
| 573 | 75 | do k=2,nz ; Rlay(k) = Rlay(k-1) + g_prime(k)*(GV%Rho0/GV%g_Earth) ; enddo |
| 574 | else | |
| 575 | 0 | do k=2,nz |
| 576 | 0 | Rlay(k) = Rlay(k-1) * ((GV%g_Earth + 0.5*g_prime(k)) / (GV%g_Earth - 0.5*g_prime(k))) |
| 577 | enddo | |
| 578 | endif | |
| 579 | ||
| 580 | 1 | call callTree_leave(trim(mdl)//'()') |
| 581 | ||
| 582 | 1 | end subroutine set_coord_to_none |
| 583 | ||
| 584 | !> Writes out a file containing any available data related | |
| 585 | !! to the vertical grid used by the MOM ocean model. | |
| 586 | 0 | subroutine write_vertgrid_file(GV, US, param_file, directory) |
| 587 | type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure | |
| 588 | type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type | |
| 589 | type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters | |
| 590 | character(len=*), intent(in) :: directory !< The directory into which to place the file. | |
| 591 | ! Local variables | |
| 592 | character(len=240) :: filepath | |
| 593 | 0 | type(vardesc) :: vars(2) |
| 594 | 0 | type(MOM_field) :: fields(2) |
| 595 | 0 | type(MOM_netCDF_file) :: IO_handle ! The I/O handle of the fileset |
| 596 | ||
| 597 | 0 | filepath = trim(directory) // trim("Vertical_coordinate.nc") |
| 598 | ||
| 599 | 0 | vars(1) = var_desc("R","kilogram meter-3","Target Potential Density",'1','L','1') |
| 600 | 0 | vars(2) = var_desc("g","meter second-2","Reduced gravity",'1','i','1') |
| 601 | ||
| 602 | call create_MOM_file(IO_handle, trim(filepath), vars, 2, fields, & | |
| 603 | 0 | SINGLE_FILE, GV=GV) |
| 604 | ||
| 605 | 0 | call MOM_write_field(IO_handle, fields(1), GV%Rlay, unscale=US%R_to_kg_m3) |
| 606 | 0 | call MOM_write_field(IO_handle, fields(2), GV%g_prime, unscale=US%L_T_to_m_s**2*US%m_to_Z) |
| 607 | ||
| 608 | 0 | call IO_handle%close() |
| 609 | ||
| 610 | 0 | end subroutine write_vertgrid_file |
| 611 | ||
| 612 | end module MOM_coord_initialization |