MOM6
Functions/Subroutines | Variables
isomip_initialization Module Reference

Detailed Description

The module configures the ISOMIP test case.

Functions/Subroutines

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

Variables

character(len=40) mdl = "ISOMIP_initialization"
 

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(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]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 434 of file ISOMIP_initialization.F90.

References mom_ale_sponge::initialize_ale_sponge(), mom_sponge::initialize_sponge(), mom_error_handler::mom_error(), regrid_consts::regridding_sigma_shelf_zstar, mom_ale_sponge::set_up_ale_sponge_field(), and mom_sponge::set_up_sponge_field().

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

Definition at line 271 of file ISOMIP_initialization.F90.

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

Referenced by mom_state_initialization::mom_initialize_state().

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

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

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

Initialization of topography.

Parameters
[in]gThe dynamic horizontal grid type
[out]dOcean bottom depth in m
[in]param_fileParameter file structure
[in]max_depthMaximum depth of model in m

Definition at line 64 of file ISOMIP_initialization.F90.

References mom_error_handler::mom_mesg().

Referenced by mom_fixed_initialization::mom_initialize_topography().

64  type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
65  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
66  intent(out) :: d !< Ocean bottom depth in m
67  type(param_file_type), intent(in) :: param_file !< Parameter file structure
68  real, intent(in) :: max_depth !< Maximum depth of model in m
69 
70 ! This subroutine sets up the ISOMIP topography
71  real :: min_depth ! The minimum and maximum depths in m.
72 
73 ! The following variables are used to set up the bathymetry in the ISOMIP example.
74 ! check this paper: http://www.geosci-model-dev-discuss.net/8/9859/2015/gmdd-8-9859-2015.pdf
75 
76  real :: bmax ! max depth of bedrock topography
77  real :: b0,b2,b4,b6 ! first, second, third and fourth bedrock topography coeff
78  real :: xbar ! characteristic along-flow lenght scale of the bedrock
79  real :: dc ! depth of the trough compared with side walls
80  real :: fc ! characteristic width of the side walls of the channel
81  real :: wc ! half-width of the trough
82  real :: ly ! domain width (across ice flow)
83  real :: bx, by, xtil ! dummy vatiables
84  logical :: is_2d ! If true, use 2D setup
85 
86 ! G%ieg and G%jeg are the last indices in the global domain
87 
88 ! This include declares and sets the variable "version".
89 #include "version_variable.h"
90  character(len=40) :: mdl = "ISOMIP_initialize_topography" ! This subroutine's name.
91  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
92  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
93  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
94 
95 
96  call mom_mesg(" ISOMIP_initialization.F90, ISOMIP_initialize_topography: setting topography", 5)
97 
98  call log_version(param_file, mdl, version, "")
99  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
100  "The minimum depth of the ocean.", units="m", default=0.0)
101  call get_param(param_file, mdl, "ISOMIP_2D",is_2d,'If true, use a 2D setup.', default=.false.)
102 
103 ! The following variables should be transformed into runtime parameters?
104  bmax=720.0; b0=-150.0; b2=-728.8; b4=343.91; b6=-50.57
105  xbar=300.0e3; dc=500.0; fc=4.0e3; wc=24.0e3; ly=80.0e3
106  bx = 0.0; by = 0.0; xtil = 0.0
107 
108 
109  if (is_2d) then
110  do j=js,je ; do i=is,ie
111  ! 2D setup
112  xtil = g%geoLonT(i,j)*1.0e3/xbar
113  !xtil = 450*1.0e3/xbar
114  bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6
115  !by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + &
116  ! (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc)))
117 
118  ! slice at y = 40 km
119  by = (dc/(1.+exp(-2.*(40.0*1.0e3- ly/2. - wc)/fc))) + &
120  (dc/(1.+exp(2.*(40.0*1.0e3- ly/2. + wc)/fc)))
121 
122  d(i,j) = -max(bx+by,-bmax)
123  if (d(i,j) > max_depth) d(i,j) = max_depth
124  if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
125  enddo ; enddo
126 
127  else
128  do j=js,je ; do i=is,ie
129  ! 3D setup
130  ! #### TEST #######
131  !if (G%geoLonT(i,j)<500.) then
132  ! xtil = 500.*1.0e3/xbar
133  !else
134  ! xtil = G%geoLonT(i,j)*1.0e3/xbar
135  !endif
136  ! ##### TEST #####
137 
138  xtil = g%geoLonT(i,j)*1.0e3/xbar
139 
140  bx = b0+b2*xtil**2 + b4*xtil**4 + b6*xtil**6
141  by = (dc/(1.+exp(-2.*(g%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + &
142  (dc/(1.+exp(2.*(g%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc)))
143 
144  d(i,j) = -max(bx+by,-bmax)
145  if (d(i,j) > max_depth) d(i,j) = max_depth
146  if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
147  enddo ; enddo
148  endif
149 
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

Definition at line 47 of file ISOMIP_initialization.F90.

Referenced by isomip_initialize_temperature_salinity(), and isomip_initialize_thickness().

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