MOM6
isomip_initialization Module Reference

Detailed Description

Configures the ISOMIP test case.

See this paper for details: http://www.geosci-model-dev-discuss.net/8/9859/2015/gmdd-8-9859-2015.pdf

Functions/Subroutines

subroutine, public isomip_initialize_topography (D, G, param_file, max_depth, US)
 Initialization of topography for the ISOMIP configuration. More...
 
subroutine, public isomip_initialize_thickness (h, G, GV, US, param_file, tv, just_read_params)
 Initialization of thicknesses. More...
 
subroutine, public isomip_initialize_temperature_salinity (T, S, h, G, GV, US, param_file, eqn_of_state, just_read_params)
 Initial values for temperature and salinity. More...
 
subroutine, public isomip_initialize_sponges (G, GV, US, tv, PF, use_ALE, CSp, ACSp)
 Sets up the the inverse restoration time (Idamp), and. More...
 

Variables

character(len=40) mdl = "ISOMIP_initialization"
 This module's name. More...
 

Function/Subroutine Documentation

◆ isomip_initialize_sponges()

subroutine, public isomip_initialization::isomip_initialize_sponges ( type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(thermo_var_ptrs), intent(in)  tv,
type(param_file_type), intent(in)  PF,
logical, intent(in)  use_ALE,
type(sponge_cs), pointer  CSp,
type(ale_sponge_cs), pointer  ACSp 
)

Sets up the the inverse restoration time (Idamp), and.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]tvA structure containing pointers to any available thermodynamic fields, potential temperature and salinity or mixed layer density. Absent fields have NULL ptrs.
[in]pfA structure indicating the open file to parse for model parameter values.
[in]use_aleIf true, indicates model is in ALE mode
cspLayer-mode sponge structure
acspALE-mode sponge structure

Definition at line 425 of file ISOMIP_initialization.F90.

425  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
426  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
427  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
428  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers
429  !! to any available thermodynamic
430  !! fields, potential temperature and
431  !! salinity or mixed layer density.
432  !! Absent fields have NULL ptrs.
433  type(param_file_type), intent(in) :: PF !< A structure indicating the
434  !! open file to parse for model
435  !! parameter values.
436  logical, intent(in) :: use_ALE !< If true, indicates model is in ALE mode
437  type(sponge_CS), pointer :: CSp !< Layer-mode sponge structure
438  type(ALE_sponge_CS), pointer :: ACSp !< ALE-mode sponge structure
439  ! Local variables
440  real :: T(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for temp
441  real :: S(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for salt
442  real :: RHO(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for RHO
443  real :: h(SZI_(G),SZJ_(G),SZK_(G)) ! A temporary array for thickness
444  real :: Idamp(SZI_(G),SZJ_(G)) ! The inverse damping rate [s-1].
445  real :: TNUDG ! Nudging time scale, days
446  real :: S_sur, T_sur ! Surface salinity and temerature in sponge
447  real :: S_bot, T_bot ! Bottom salinity and temerature in sponge
448  real :: t_ref, s_ref ! reference T and S
449  real :: rho_sur, rho_bot ! Surface and bottom densities [R ~> kg m-3]
450  real :: rho_range ! The range of densities [R ~> kg m-3]
451  real :: dT_dz, dS_dz ! Gradients of T and S in degC/Z and PPT/Z.
452 
453  real :: e0(SZK_(G)+1) ! The resting interface heights [Z ~> m], usually
454  ! negative because it is positive upward.
455  real :: eta1D(SZK_(G)+1) ! Interface height relative to the sea surface, positive upward [Z ~> m].
456  real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1) ! A temporary array for eta [Z ~> m].
457  real :: min_depth, dummy1, z
458  real :: damp, rho_dummy, min_thickness, rho_tmp, xi0
459  character(len=40) :: verticalCoordinate, filename, state_file
460  character(len=40) :: temp_var, salt_var, eta_var, inputdir
461 
462  character(len=40) :: mdl = "ISOMIP_initialize_sponges" ! This subroutine's name.
463  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
464 
465  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
466  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
467 
468  call get_param(pf, mdl, "MIN_THICKNESS", min_thickness, "Minimum layer thickness", &
469  units="m", default=1.e-3, scale=us%m_to_Z)
470 
471  call get_param(pf, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
472  default=default_coordinate_mode)
473 
474  call get_param(pf, mdl, "ISOMIP_TNUDG", tnudg, "Nudging time scale for sponge layers (days)", default=0.0)
475 
476  call get_param(pf, mdl, "T_REF", t_ref, "Reference temperature", default=10.0,&
477  do_not_log=.true.)
478 
479  call get_param(pf, mdl, "S_REF", s_ref, "Reference salinity", default=35.0,&
480  do_not_log=.true.)
481 
482  call get_param(pf, mdl, "ISOMIP_S_SUR_SPONGE", s_sur, &
483  'Surface salinity in sponge layer.', default=s_ref) ! units="ppt")
484 
485  call get_param(pf, mdl, "ISOMIP_S_BOT_SPONGE", s_bot, &
486  'Bottom salinity in sponge layer.', default=s_ref) ! units="ppt")
487 
488  call get_param(pf, mdl, "ISOMIP_T_SUR_SPONGE", t_sur, &
489  'Surface temperature in sponge layer.', default=t_ref) ! units="degC")
490 
491  call get_param(pf, mdl, "ISOMIP_T_BOT_SPONGE", t_bot, &
492  'Bottom temperature in sponge layer.', default=t_ref) ! units="degC")
493 
494  t(:,:,:) = 0.0 ; s(:,:,:) = 0.0 ; idamp(:,:) = 0.0; rho(:,:,:) = 0.0
495 
496 ! Set up sponges for ISOMIP configuration
497  call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, &
498  "The minimum depth of the ocean.", units="m", default=0.0, scale=us%m_to_Z)
499 
500  if (associated(csp)) call mom_error(fatal, &
501  "ISOMIP_initialize_sponges called with an associated control structure.")
502  if (associated(acsp)) call mom_error(fatal, &
503  "ISOMIP_initialize_sponges called with an associated ALE-sponge control structure.")
504 
505  ! Here the inverse damping time [s-1], is set. Set Idamp to 0 !
506  ! wherever there is no sponge, and the subroutines that are called !
507  ! will automatically set up the sponges only where Idamp is positive!
508  ! and mask2dT is 1.
509 
510  do i=is,ie; do j=js,je
511  if (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0) then
512 
513  ! 1 / day
514  dummy1=(g%geoLonT(i,j)-790.0)/(800.0-790.0)
515  damp = 1.0/tnudg * max(0.0,dummy1)
516 
517  else ; damp=0.0
518  endif
519 
520  ! convert to 1 / seconds
521  if (g%bathyT(i,j) > min_depth) then
522  idamp(i,j) = damp/86400.0
523  else ; idamp(i,j) = 0.0 ; endif
524 
525  enddo ; enddo
526 
527  ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT
528  call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state, scale=us%kg_m3_to_R)
529  !write (mesg,*) 'Surface density in sponge:', rho_sur
530  ! call MOM_mesg(mesg,5)
531  call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state, scale=us%kg_m3_to_R)
532  !write (mesg,*) 'Bottom density in sponge:', rho_bot
533  ! call MOM_mesg(mesg,5)
534  rho_range = rho_bot - rho_sur
535  !write (mesg,*) 'Density range in sponge:', rho_range
536  ! call MOM_mesg(mesg,5)
537 
538  if (use_ale) then
539 
540  select case ( coordinatemode(verticalcoordinate) )
541 
542  case ( regridding_rho )
543  ! Construct notional interface positions
544  e0(1) = 0.
545  do k=2,nz
546  e0(k) = -g%max_depth * ( 0.5 * ( gv%Rlay(k-1) + gv%Rlay(k) ) - rho_sur ) / rho_range
547  e0(k) = min( 0., e0(k) ) ! Bound by surface
548  e0(k) = max( -g%max_depth, e0(k) ) ! Bound by possible deepest point in model
549  ! write(mesg,*) 'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',&
550  ! G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)
551  ! call MOM_mesg(mesg,5)
552  enddo
553  e0(nz+1) = -g%max_depth
554 
555  ! Calculate thicknesses
556  do j=js,je ; do i=is,ie
557  eta1d(nz+1) = -g%bathyT(i,j)
558  do k=nz,1,-1
559  eta1d(k) = e0(k)
560  if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
561  eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
562  h(i,j,k) = gv%Angstrom_H
563  else
564  h(i,j,k) = gv%Z_to_H*(eta1d(k) - eta1d(k+1))
565  endif
566  enddo
567  enddo ; enddo
568 
569  case ( regridding_zstar, regridding_sigma_shelf_zstar ) ! Initial thicknesses for z coordinates
570  do j=js,je ; do i=is,ie
571  eta1d(nz+1) = -g%bathyT(i,j)
572  do k=nz,1,-1
573  eta1d(k) = -g%max_depth * real(k-1) / real(nz)
574  if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
575  eta1d(k) = eta1d(k+1) + min_thickness
576  h(i,j,k) = min_thickness * gv%Z_to_H
577  else
578  h(i,j,k) = gv%Z_to_H*(eta1d(k) - eta1d(k+1))
579  endif
580  enddo
581  enddo ; enddo
582 
583  case ( regridding_sigma ) ! Initial thicknesses for sigma coordinates
584  do j=js,je ; do i=is,ie
585  h(i,j,:) = gv%Z_to_H * (g%bathyT(i,j) / dfloat(nz))
586  enddo ; enddo
587 
588  case default
589  call mom_error(fatal,"ISOMIP_initialize_sponges: "// &
590  "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
591 
592  end select
593 
594  ! This call sets up the damping rates and interface heights.
595  ! This sets the inverse damping timescale fields in the sponges.
596  call initialize_ale_sponge(idamp, g, pf, acsp, h, nz)
597 
598  ds_dz = (s_sur - s_bot) / g%max_depth
599  dt_dz = (t_sur - t_bot) / g%max_depth
600  do j=js,je ; do i=is,ie
601  xi0 = -g%bathyT(i,j)
602  do k = nz,1,-1
603  xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z ! Depth in middle of layer
604  s(i,j,k) = s_sur + ds_dz * xi0
605  t(i,j,k) = t_sur + dt_dz * xi0
606  xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z ! Depth at top of layer
607  enddo
608  enddo ; enddo
609  ! for debugging
610  !i=G%iec; j=G%jec
611  !do k = 1,nz
612  ! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state, scale=US%kg_m3_to_R)
613  ! write(mesg,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k)
614  ! call MOM_mesg(mesg,5)
615  !enddo
616 
617  ! Now register all of the fields which are damped in the sponge. !
618  ! By default, momentum is advected vertically within the sponge, but !
619  ! momentum is typically not damped within the sponge. !
620 
621  ! The remaining calls to set_up_sponge_field can be in any order. !
622  if ( associated(tv%T) ) then
623  call set_up_ale_sponge_field(t, g, tv%T, acsp)
624  endif
625  if ( associated(tv%S) ) then
626  call set_up_ale_sponge_field(s, g, tv%S, acsp)
627  endif
628 
629  else ! layer mode
630  ! 1) Read eta, salt and temp from IC file
631  call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
632  inputdir = slasher(inputdir)
633  ! GM: get two different files, one with temp and one with salt values
634  ! this is work around to avoid having wrong values near the surface
635  ! because of the FIT_SALINITY option. To get salt values right in the
636  ! sponge, FIT_SALINITY=False. The oposite is true for temp. One can
637  ! combined the *correct* temp and salt values in one file instead.
638  call get_param(pf, mdl, "ISOMIP_SPONGE_FILE", state_file, &
639  "The name of the file with temps., salts. and interfaces to "//&
640  "damp toward.", fail_if_missing=.true.)
641  call get_param(pf, mdl, "SPONGE_PTEMP_VAR", temp_var, &
642  "The name of the potential temperature variable in "//&
643  "SPONGE_STATE_FILE.", default="Temp")
644  call get_param(pf, mdl, "SPONGE_SALT_VAR", salt_var, &
645  "The name of the salinity variable in "//&
646  "SPONGE_STATE_FILE.", default="Salt")
647  call get_param(pf, mdl, "SPONGE_ETA_VAR", eta_var, &
648  "The name of the interface height variable in "//&
649  "SPONGE_STATE_FILE.", default="eta")
650 
651  !read temp and eta
652  filename = trim(inputdir)//trim(state_file)
653  if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
654  "ISOMIP_initialize_sponges: Unable to open "//trim(filename))
655  call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
656  call mom_read_data(filename, temp_var, t(:,:,:), g%Domain)
657  call mom_read_data(filename, salt_var, s(:,:,:), g%Domain)
658 
659  ! for debugging
660  !i=G%iec; j=G%jec
661  !do k = 1,nz
662  ! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state, scale=US%kg_m3_to_R)
663  ! write(mesg,*) 'Sponge - k,eta,T,S,rho,Rlay',k,eta(i,j,k),T(i,j,k),&
664  ! S(i,j,k),rho_tmp,GV%Rlay(k)
665  ! call MOM_mesg(mesg,5)
666  !enddo
667 
668  ! Set the inverse damping rates so that the model will know where to
669  ! apply the sponges, along with the interface heights.
670  call initialize_sponge(idamp, eta, g, pf, csp, gv)
671  ! Apply sponge in tracer fields
672  call set_up_sponge_field(t, tv%T, g, nz, csp)
673  call set_up_sponge_field(s, tv%S, g, nz, csp)
674 
675  endif
676 

References mom_sponge::initialize_sponge(), mdl, mom_error_handler::mom_error(), regrid_consts::regridding_sigma_shelf_zstar, and mom_sponge::set_up_sponge_field().

Here is the call graph for this function:

◆ isomip_initialize_temperature_salinity()

subroutine, public isomip_initialization::isomip_initialize_temperature_salinity ( real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  T,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(out)  S,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(eos_type), pointer  eqn_of_state,
logical, intent(in), optional  just_read_params 
)

Initial values for temperature and salinity.

Parameters
[in]gOcean grid structure
[in]gvVertical grid structure
[in]usA dimensional unit scaling type
[out]tPotential temperature [degC]
[out]sSalinity [ppt]
[in]hLayer thickness [H ~> m or kg m-2]
[in]param_fileParameter file structure
eqn_of_stateEquation of state structure
[in]just_read_paramsIf present and true, this call will only read parameters without changing T & S.

Definition at line 255 of file ISOMIP_initialization.F90.

255  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
256  type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
257  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
258  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: T !< Potential temperature [degC]
259  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity [ppt]
260  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
261  type(param_file_type), intent(in) :: param_file !< Parameter file structure
262  type(EOS_type), pointer :: eqn_of_state !< Equation of state structure
263  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
264  !! only read parameters without changing T & S.
265  ! Local variables
266  integer :: i, j, k, is, ie, js, je, nz, itt
267  real :: x, ds, dt
268  real :: rho_sur, rho_bot ! Surface and bottom densities [R ~> kg m-3]
269  real :: xi0, xi1 ! Heights in depth units [Z ~> m].
270  real :: S_sur, S_bot ! Salinity at the surface and bottom [ppt]
271  real :: T_sur, T_bot ! Temperature at the bottom [degC]
272  real :: dT_dz ! Vertical gradient of temperature [degC Z-1 ~> degC m-1].
273  real :: dS_dz ! Vertical gradient of salinity [ppt Z-1 ~> ppt m-1].
274  real :: z ! vertical position in z space [Z ~> m]
275  character(len=256) :: mesg ! The text of an error message
276  character(len=40) :: verticalCoordinate, density_profile
277  real :: rho_tmp
278  logical :: just_read ! If true, just read parameters but set nothing.
279  logical :: fit_salin ! If true, accept the prescribed temperature and fit the salinity.
280  real :: T0(SZK_(G)), S0(SZK_(G))
281  real :: drho_dT(SZK_(G)) ! Derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
282  real :: drho_dS(SZK_(G)) ! Derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
283  real :: rho_guess(SZK_(G)) ! Potential density at T0 & S0 [R ~> kg m-3].
284  real :: pres(SZK_(G)) ! An array of the reference pressure [Pa]. (zero here)
285  real :: drho_dT1 ! A prescribed derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
286  real :: drho_dS1 ! A prescribed derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
287  real :: T_Ref, S_Ref
288  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
289  pres(:) = 0.0
290 
291  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
292 
293  call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
294  default=default_coordinate_mode, do_not_log=just_read)
295  call get_param(param_file, mdl, "ISOMIP_T_SUR",t_sur, &
296  'Temperature at the surface (interface)', default=-1.9, do_not_log=just_read)
297  call get_param(param_file, mdl, "ISOMIP_S_SUR", s_sur, &
298  'Salinity at the surface (interface)', default=33.8, do_not_log=just_read)
299  call get_param(param_file, mdl, "ISOMIP_T_BOT", t_bot, &
300  'Temperature at the bottom (interface)', default=-1.9, do_not_log=just_read)
301  call get_param(param_file, mdl, "ISOMIP_S_BOT", s_bot, &
302  'Salinity at the bottom (interface)', default=34.55, do_not_log=just_read)
303 
304  call calculate_density(t_sur,s_sur,0.0,rho_sur,eqn_of_state, scale=us%kg_m3_to_R)
305  ! write(mesg,*) 'Density in the surface layer:', rho_sur
306  ! call MOM_mesg(mesg,5)
307  call calculate_density(t_bot,s_bot,0.0,rho_bot,eqn_of_state, scale=us%kg_m3_to_R)
308  ! write(mesg,*) 'Density in the bottom layer::', rho_bot
309  ! call MOM_mesg(mesg,5)
310 
311  select case ( coordinatemode(verticalcoordinate) )
312 
313  case ( regridding_rho, regridding_zstar, regridding_sigma_shelf_zstar, regridding_sigma )
314  if (just_read) return ! All run-time parameters have been read, so return.
315 
316  ds_dz = (s_sur - s_bot) / g%max_depth
317  dt_dz = (t_sur - t_bot) / g%max_depth
318  do j=js,je ; do i=is,ie
319  xi0 = -g%bathyT(i,j)
320  do k = nz,1,-1
321  xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z ! Depth in middle of layer
322  s(i,j,k) = s_sur + ds_dz * xi0
323  t(i,j,k) = t_sur + dt_dz * xi0
324  xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z ! Depth at top of layer
325  enddo
326  enddo ; enddo
327 
328  case ( regridding_layer )
329  call get_param(param_file, mdl, "FIT_SALINITY", fit_salin, &
330  "If true, accept the prescribed temperature and fit the "//&
331  "salinity; otherwise take salinity and fit temperature.", &
332  default=.false., do_not_log=just_read)
333  call get_param(param_file, mdl, "DRHO_DS", drho_ds1, &
334  "Partial derivative of density with salinity.", &
335  units="kg m-3 PSU-1", scale=us%kg_m3_to_R, fail_if_missing=.not.just_read, do_not_log=just_read)
336  call get_param(param_file, mdl, "DRHO_DT", drho_dt1, &
337  "Partial derivative of density with temperature.", &
338  units="kg m-3 K-1", scale=us%kg_m3_to_R, fail_if_missing=.not.just_read, do_not_log=just_read)
339  call get_param(param_file, mdl, "T_REF", t_ref, &
340  "A reference temperature used in initialization.", &
341  units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
342  call get_param(param_file, mdl, "S_REF", s_ref, &
343  "A reference salinity used in initialization.", units="PSU", &
344  default=35.0, do_not_log=just_read)
345  if (just_read) return ! All run-time parameters have been read, so return.
346 
347  ! write(mesg,*) 'read drho_dS, drho_dT', drho_dS1, drho_dT1
348  ! call MOM_mesg(mesg,5)
349 
350  ds_dz = (s_sur - s_bot) / g%max_depth
351  dt_dz = (t_sur - t_bot) / g%max_depth
352 
353  do j=js,je ; do i=is,ie
354  xi0 = 0.0
355  do k = 1,nz
356  !T0(k) = T_Ref; S0(k) = S_Ref
357  xi1 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
358  s0(k) = s_sur - ds_dz * xi1
359  t0(k) = t_sur - dt_dz * xi1
360  xi0 = xi0 + h(i,j,k) * gv%H_to_Z
361  ! write(mesg,*) 'S,T,xi0,xi1,k',S0(k),T0(k),xi0,xi1,k
362  ! call MOM_mesg(mesg,5)
363  enddo
364 
365  call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,1,eqn_of_state, scale=us%kg_m3_to_R)
366  ! write(mesg,*) 'computed drho_dS, drho_dT', drho_dS(1), drho_dT(1)
367  ! call MOM_mesg(mesg,5)
368  call calculate_density(t0(1),s0(1),0.,rho_guess(1),eqn_of_state, scale=us%kg_m3_to_R)
369 
370  if (fit_salin) then
371  ! A first guess of the layers' salinity.
372  do k=nz,1,-1
373  s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds1)
374  enddo
375  ! Refine the guesses for each layer.
376  do itt=1,6
377  call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
378  call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
379  do k=1,nz
380  s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds1)
381  enddo
382  enddo
383 
384  else
385  ! A first guess of the layers' temperatures.
386  do k=nz,1,-1
387  t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt1
388  enddo
389 
390  do itt=1,6
391  call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
392  call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
393  do k=1,nz
394  t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
395  enddo
396  enddo
397  endif
398 
399  do k=1,nz
400  t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
401  enddo
402 
403  enddo ; enddo
404 
405  case default
406  call mom_error(fatal,"isomip_initialize: "// &
407  "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
408 
409  end select
410 
411  ! for debugging
412  !i=G%iec; j=G%jec
413  !do k = 1,nz
414  ! call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,eqn_of_state, scale=US%kg_m3_to_R)
415  ! write(mesg,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k)
416  ! call MOM_mesg(mesg,5)
417  !enddo
418 

References mdl, mom_error_handler::mom_error(), and regrid_consts::regridding_sigma_shelf_zstar.

Referenced by mom_state_initialization::mom_initialize_state().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ isomip_initialize_thickness()

subroutine, public isomip_initialization::isomip_initialize_thickness ( real, dimension( g %isd: g %ied, g %jsd: g %jed, gv %ke), intent(out)  h,
type(ocean_grid_type), intent(in)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(thermo_var_ptrs), intent(in)  tv,
logical, intent(in), optional  just_read_params 
)

Initialization of thicknesses.

Parameters
[in]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[out]hThe thickness that is being initialized [H ~> m or kg m-2].
[in]param_fileA structure indicating the open file to parse for model parameter values.
[in]tvA structure containing pointers to any available thermodynamic fields, including the eqn. of state.
[in]just_read_paramsIf present and true, this call will only read parameters without changing h.

Definition at line 134 of file ISOMIP_initialization.F90.

134  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
135  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
136  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
137  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
138  intent(out) :: h !< The thickness that is being initialized [H ~> m or kg m-2].
139  type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
140  !! to parse for model parameter values.
141  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
142  !! available thermodynamic fields, including
143  !! the eqn. of state.
144  logical, optional, intent(in) :: just_read_params !< If present and true, this call will
145  !! only read parameters without changing h.
146  ! Local variables
147  real :: e0(SZK_(G)+1) ! The resting interface heights, in depth units [Z ~> m],
148  ! usually negative because it is positive upward.
149  real :: eta1D(SZK_(G)+1)! Interface height relative to the sea surface
150  ! positive upward, in depth units [Z ~> m].
151  integer :: i, j, k, is, ie, js, je, nz, tmp1
152  real :: x
153  real :: min_thickness, s_sur, s_bot, t_sur, t_bot
154  real :: rho_sur, rho_bot ! Surface and bottom densities [R ~> kg m-3]
155  real :: rho_range ! The range of densities [R ~> kg m-3]
156  logical :: just_read ! If true, just read parameters but set nothing.
157  character(len=256) :: mesg ! The text of an error message
158  character(len=40) :: verticalCoordinate
159 
160  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
161 
162  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
163 
164  if (.not.just_read) &
165  call mom_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness")
166 
167  call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, &
168  'Minimum layer thickness', units='m', default=1.e-3, do_not_log=just_read, scale=us%m_to_Z)
169  call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
170  default=default_coordinate_mode, do_not_log=just_read)
171 
172  select case ( coordinatemode(verticalcoordinate) )
173 
174  case ( regridding_layer, regridding_rho ) ! Initial thicknesses for isopycnal coordinates
175  call get_param(param_file, mdl, "ISOMIP_T_SUR",t_sur, &
176  'Temperature at the surface (interface)', default=-1.9, do_not_log=just_read)
177  call get_param(param_file, mdl, "ISOMIP_S_SUR", s_sur, &
178  'Salinity at the surface (interface)', default=33.8, do_not_log=just_read)
179  call get_param(param_file, mdl, "ISOMIP_T_BOT", t_bot, &
180  'Temperature at the bottom (interface)', default=-1.9, do_not_log=just_read)
181  call get_param(param_file, mdl, "ISOMIP_S_BOT", s_bot,&
182  'Salinity at the bottom (interface)', default=34.55, do_not_log=just_read)
183 
184  if (just_read) return ! All run-time parameters have been read, so return.
185 
186  ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT
187  call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state, scale=us%kg_m3_to_R)
188  ! write(mesg,*) 'Surface density is:', rho_sur
189  ! call MOM_mesg(mesg,5)
190  call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state, scale=us%kg_m3_to_R)
191  ! write(mesg,*) 'Bottom density is:', rho_bot
192  ! call MOM_mesg(mesg,5)
193  rho_range = rho_bot - rho_sur
194  ! write(mesg,*) 'Density range is:', rho_range
195  ! call MOM_mesg(mesg,5)
196 
197  ! Construct notional interface positions
198  e0(1) = 0.
199  do k=2,nz
200  e0(k) = -g%max_depth * ( 0.5 * ( gv%Rlay(k-1) + gv%Rlay(k) ) - rho_sur ) / rho_range
201  e0(k) = min( 0., e0(k) ) ! Bound by surface
202  e0(k) = max( -g%max_depth, e0(k) ) ! Bound by possible deepest point in model
203  ! write(mesg,*) 'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)', &
204  ! G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)
205  ! call MOM_mesg(mesg,5)
206  enddo
207  e0(nz+1) = -g%max_depth
208 
209  ! Calculate thicknesses
210  do j=js,je ; do i=is,ie
211  eta1d(nz+1) = -g%bathyT(i,j)
212  do k=nz,1,-1
213  eta1d(k) = e0(k)
214  if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
215  eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
216  h(i,j,k) = gv%Angstrom_H
217  else
218  h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
219  endif
220  enddo
221  enddo ; enddo
222 
223  case ( regridding_zstar, regridding_sigma_shelf_zstar ) ! Initial thicknesses for z coordinates
224  if (just_read) return ! All run-time parameters have been read, so return.
225  do j=js,je ; do i=is,ie
226  eta1d(nz+1) = -g%bathyT(i,j)
227  do k=nz,1,-1
228  eta1d(k) = -g%max_depth * real(k-1) / real(nz)
229  if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
230  eta1d(k) = eta1d(k+1) + min_thickness
231  h(i,j,k) = gv%Z_to_H * min_thickness
232  else
233  h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
234  endif
235  enddo
236  enddo ; enddo
237 
238  case ( regridding_sigma ) ! Initial thicknesses for sigma coordinates
239  if (just_read) return ! All run-time parameters have been read, so return.
240  do j=js,je ; do i=is,ie
241  h(i,j,:) = gv%Z_to_H * g%bathyT(i,j) / dfloat(nz)
242  enddo ; enddo
243 
244  case default
245  call mom_error(fatal,"isomip_initialize: "// &
246  "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")
247 
248  end select
249 

References mdl, mom_error_handler::mom_error(), mom_error_handler::mom_mesg(), and regrid_consts::regridding_sigma_shelf_zstar.

Here is the call graph for this function:

◆ isomip_initialize_topography()

subroutine, public isomip_initialization::isomip_initialize_topography ( real, dimension(g%isd:g%ied,g%jsd:g%jed), intent(out)  D,
type(dyn_horgrid_type), intent(in)  G,
type(param_file_type), intent(in)  param_file,
real, intent(in)  max_depth,
type(unit_scale_type), intent(in), optional  US 
)

Initialization of topography for the ISOMIP configuration.

Parameters
[in]gThe dynamic horizontal grid type
[out]dOcean bottom depth in m or Z if US is present
[in]param_fileParameter file structure
[in]max_depthMaximum model depth in the units of D
[in]usA dimensional unit scaling type

Definition at line 45 of file ISOMIP_initialization.F90.

45  type(dyn_horgrid_type), intent(in) :: G !< The dynamic horizontal grid type
46  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
47  intent(out) :: D !< Ocean bottom depth in m or Z if US is present
48  type(param_file_type), intent(in) :: param_file !< Parameter file structure
49  real, intent(in) :: max_depth !< Maximum model depth in the units of D
50  type(unit_scale_type), optional, intent(in) :: US !< A dimensional unit scaling type
51 
52  ! Local variables
53  real :: min_depth ! The minimum and maximum depths [Z ~> m].
54  real :: m_to_Z ! A dimensional rescaling factor.
55  ! The following variables are used to set up the bathymetry in the ISOMIP example.
56  real :: bmax ! max depth of bedrock topography
57  real :: b0,b2,b4,b6 ! first, second, third and fourth bedrock topography coeff
58  real :: xbar ! characteristic along-flow lenght scale of the bedrock
59  real :: dc ! depth of the trough compared with side walls [Z ~> m].
60  real :: fc ! characteristic width of the side walls of the channel
61  real :: wc ! half-width of the trough
62  real :: ly ! domain width (across ice flow)
63  real :: bx, by ! dummy vatiables [Z ~> m].
64  real :: xtil ! dummy vatiable
65  logical :: is_2D ! If true, use 2D setup
66 ! This include declares and sets the variable "version".
67 #include "version_variable.h"
68  character(len=40) :: mdl = "ISOMIP_initialize_topography" ! This subroutine's name.
69  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
70  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
71  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
72 
73  call mom_mesg(" ISOMIP_initialization.F90, ISOMIP_initialize_topography: setting topography", 5)
74 
75  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
76 
77  call log_version(param_file, mdl, version, "")
78  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
79  "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_z)
80  call get_param(param_file, mdl, "ISOMIP_2D",is_2d,'If true, use a 2D setup.', default=.false.)
81 
82  ! The following variables should be transformed into runtime parameters?
83  bmax = 720.0*m_to_z ; dc = 500.0*m_to_z
84  b0 = -150.0*m_to_z ; b2 = -728.8*m_to_z ; b4 = 343.91*m_to_z ; b6 = -50.57*m_to_z
85  xbar = 300.0e3 ; fc = 4.0e3 ; wc = 24.0e3 ; ly = 80.0e3
86  bx = 0.0 ; by = 0.0 ; xtil = 0.0
87 
88 
89  if (is_2d) then
90  do j=js,je ; do i=is,ie
91  ! 2D setup
92  xtil = g%geoLonT(i,j)*1.0e3/xbar
93  !xtil = 450*1.0e3/xbar
94  bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6
95  !by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + &
96  ! (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc)))
97 
98  ! slice at y = 40 km
99  by = (dc / (1.+exp(-2.*(40.0*1.0e3- ly/2. - wc)/fc))) + &
100  (dc / (1.+exp(2.*(40.0*1.0e3- ly/2. + wc)/fc)))
101 
102  d(i,j) = -max(bx+by, -bmax)
103  if (d(i,j) > max_depth) d(i,j) = max_depth
104  if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
105  enddo ; enddo
106 
107  else
108  do j=js,je ; do i=is,ie
109  ! 3D setup
110  ! ===== TEST =====
111  !if (G%geoLonT(i,j)<500.) then
112  ! xtil = 500.*1.0e3/xbar
113  !else
114  ! xtil = G%geoLonT(i,j)*1.0e3/xbar
115  !endif
116  ! ===== TEST =====
117 
118  xtil = g%geoLonT(i,j)*1.0e3/xbar
119 
120  bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6
121  by = (dc / (1.+exp(-2.*(g%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + &
122  (dc / (1.+exp(2.*(g%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc)))
123 
124  d(i,j) = -max(bx+by, -bmax)
125  if (d(i,j) > max_depth) d(i,j) = max_depth
126  if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
127  enddo ; enddo
128  endif
129 

References mdl, and mom_error_handler::mom_mesg().

Referenced by mom_fixed_initialization::mom_initialize_topography().

Here is the call graph for this function:
Here is the caller graph for this function:

Variable Documentation

◆ mdl

character(len=40) isomip_initialization::mdl = "ISOMIP_initialization"
private

This module's name.

Definition at line 28 of file ISOMIP_initialization.F90.

28 character(len=40) :: mdl = "ISOMIP_initialization" !< This module's name.

Referenced by isomip_initialize_sponges(), isomip_initialize_temperature_salinity(), isomip_initialize_thickness(), and isomip_initialize_topography().