!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! !!
!! GNU General Public License !!
!! !!
!! This file is part of the Flexible Modeling System (FMS). !!
!! !!
!! FMS is free software; you can redistribute it and/or modify !!
!! it and are expected to follow the terms of the GNU General Public !!
!! License as published by the Free Software Foundation. !!
!! !!
!! FMS is distributed in the hope that it will be useful, !!
!! but WITHOUT ANY WARRANTY; without even the implied warranty of !!
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the !!
!! GNU General Public License for more details. !!
!! !!
!! You should have received a copy of the GNU General Public License !!
!! along with FMS; if not, write to: !!
!! Free Software Foundation, Inc. !!
!! 59 Temple Place, Suite 330 !!
!! Boston, MA 02111-1307 USA !!
!! or see: !!
!! http://www.gnu.org/licenses/gpl.txt !!
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module ocean_thickness_mod
!
! S.M. Griffies
!
!
!
! Determine thickness of grid cells.
!
!
!
! This module determines the thickness of grid cells.
! Thicknesses are generally time dependent and functions
! of the vertical coordinate.
!
!
!
!
!
! S.M. Griffies, Elements of mom4p1 (2006)
!
!
!
!
!
!
! For debugging.
!
!
! For debugging pressure coordinate models. Lots of grid information
! printed.
!
!
! Set true to write a restart. False setting only for rare
! cases where wish to benchmark model without measuring the cost
! of writing restarts and associated chksums.
! Default is write_a_restart=.true.
!
!
!
! For debugging, set the thickness of top cell in
! geopotential model to time independent values.
! This option is needed if use the kappa_sort diagnostic.
! Default linear_free_surface=.false.
!
!
!
! To determine whether use energetic method or finite volume method
! to compute the thickness of a grid cell. Options are
! thickness_method=energetic or thickness_method=finitevolume.
! There is little overall difference in results for pbot and eta.
! However, it has been found that for realistic bottom topography
! simulations, the vertical velocity component is very noisy with
! the finitevolume approach. So this approach is considered
! experimental. The default is thickness_method='energetic'.
!
!
! For case where with to only have the dzt be determined by the full step
! bottom topography. This nml option is provided only for backwards
! compatibility with older mom experiments using the full step topog.
!
!
! For cases where wish to run model even with negative thickness.
! Default enforce_positive_dzt=.false.
!
!
! For determining how strict we are to check for the thickness of
! a column when initializing pressure based vertical coordinate models.
!
!
!
! Minimum dzt set when enforce_positive_dzt set true.
! Default thickness_dzt_min=1.0.
!
!
!
! For pressure-based models, we can (with some work) initialize the
! model to have a zero surface height. The recomended approach is to
! allow the surface height to be whatever it wants to be, and let adjustments
! smooth it over time. Default initialize_zero_eta=.false.
!
!
! For determining a modified bottom depth array
! that is required to ensure pressure model, based on initial
! in-situ density, retains a nontrivial bottom cell thickness
! in the case when initialize_zero_eta=0.0
! Default thickness_dzt_min_init=5.0.
!
!
! Expedient to allow for the computation of ht_mod.
! in the case when initialize_zero_eta=.true. Here,
! we run the pressure based model with a rescaled mass that
! is sufficient to maintain non-negative dzt, at least for
! a short period. This allows for one to run a day integration
! to produce ht_mod. rescale_mass_to_get_ht_mod=.true. will
! produce spurious results in general due to problems with the
! pressure gradient computation. So it is not recommended for
! more than initial day or so. Default rescale_mass_to_get_ht_mod=.false.
!
!
!
! For reading in a basin mask of use to re-define rho0 in
! isolated regions such as the Black Sea. This is used for
! modifying the definition of the pressure or pstar levels
! during the initialization of the thicknesses dst. This
! approach is appropriate in general, but has only been tested
! when modify the pressure levels within a fully enclosed basin.
! Default read_rescale_rho0_mask=.false.
!
!
! For specifying the rescale_rho0_mask based on reading in
! the GFDL regional mask. Default rescale_rho0_mask_gfdl=.false.
!
!
! For rescaling rho0 in a basin with a number rescale_rho0_basin_label.
! For the Black Sea using GFDL basin masks in OM3,
! rescale_rho0_basin_label=7.0. Default rescale_rho0_basin_label=-1.0
!
!
! Fractional value for rescaling rho0 in the a region.
! Default rescale_rho0_value=1.0.
!
!
!
! For sigma coordinates, have minimum depth so that have layers defined
! globally. Masks will zero out results over land, but for numerics
! it is useful to compute everywhere. Default depth_min_for_sigma=0.01.
!
!
!
! To read in an initial rho0(z) profile to assist in defining the
! initial settings for the pressure increments dst, for use in
! setting the pressure-based vertical coordinate grids. Ideally,
! this profile is determined by the level averaged density in
! the initial conditions. Note that it is essential to have
! rho0_profile have a sensible value at all depths even if there
! is no water there, since there are places where we divide by
! rho0_profile in rock. Also, be mindful that with denser water
! at depth, the pressure levels will be coarser at depth than if
! using the trivial density profile rho0(k)=rho0.
! This option is experimental, so it is recommended that user
! maintain the default read_rho0_profile=.false.
!
!
!
! For testing purposes, have this option compute pbot0=g*rho0*ht
! with rho0= constant. Default pbot0_simple=.false.
!
!
!
! A bug in certain versions of MOM4p1 was present, whereby
! Thickness%dzwu(i,j,k=0) was never updated, except for GEOPOTENTIAL
! vertical coordinates. This logical, whose default is
! update_dzwu_k0=.true., is provided for legacy purposes.
! To test the older results, have update_dzwu_k0=.false.
!
!
!
! Maximum bad grid cells printout for identifying problematic simulations.
! Default max_num_bad_print=25.
!
!
!
use constants_mod, only: epsln, grav, c2dbars
use diag_manager_mod, only: register_diag_field, register_static_field, send_data
use fms_mod, only: write_version_number, error_mesg, FATAL, WARNING
use fms_mod, only: read_data
use fms_mod, only: open_namelist_file, close_file, check_nml_error, file_exist
use fms_io_mod, only: register_restart_field, save_restart, restore_state
use fms_io_mod, only: restart_file_type, reset_field_pointer
use mpp_domains_mod, only: mpp_update_domains, mpp_global_field, mpp_global_min, mpp_global_max, domain2d
use mpp_mod, only: stdout, stdlog, mpp_error, mpp_chksum, mpp_max, mpp_min, mpp_pe
use ocean_domains_mod, only: get_local_indices
use ocean_grids_mod, only: update_boundaries
use ocean_parameters_mod, only: GEOPOTENTIAL, ZSTAR, ZSIGMA, PRESSURE, PSTAR, PSIGMA
use ocean_parameters_mod, only: PRESSURE_BASED, DEPTH_BASED, TERRAIN_FOLLOWING
use ocean_parameters_mod, only: TWO_LEVEL, THREE_LEVEL
use ocean_parameters_mod, only: ENERGETIC, FINITEVOLUME
use ocean_parameters_mod, only: missing_value, rho0, rho0r
use ocean_tracer_util_mod,only: dzt_min_max
use ocean_types_mod, only: ocean_time_type, ocean_domain_type, ocean_external_mode_type
use ocean_types_mod, only: ocean_grid_type, ocean_thickness_type, ocean_density_type
use ocean_types_mod, only: ocean_time_steps_type
use ocean_util_mod, only: write_timestamp
use ocean_workspace_mod, only: wrk1, wrk2, wrk3, wrk4, wrk1_2d
implicit none
private
#include
! for vertical coordinate choice
integer :: vert_coordinate
integer :: vert_coordinate_class
integer :: vert_coordinate_type
! useful constant
real :: grav_r
real :: grav_rho0
! for remap operator
integer :: halo
! to write some output to all processors
integer :: unit=6
! for diagnostics
integer :: id_dst = -1
integer :: id_dzt = -1
integer :: id_dzt_pbc = -1
integer :: id_dztlo = -1
integer :: id_dztup = -1
integer :: id_dstlo = -1
integer :: id_dstup = -1
integer :: id_dzwt = -1
integer :: id_dzt_dst = -1
integer :: id_dzu = -1
integer :: id_dzwu = -1
integer :: id_rho_dzt = -1
integer :: id_rho_dzu = -1
integer :: id_rho_dzt_tendency = -1
integer :: id_pbot0 = -1
integer :: id_depth_st = -1
integer :: id_depth_swt = -1
integer :: id_depth_zt = -1
integer :: id_geodepth_zt = -1
integer :: id_depth_zwt = -1
integer :: id_depth_zu = -1
integer :: id_depth_zwu = -1
integer :: id_mass_u = -1
integer :: id_thicku = -1
integer :: id_rescale_mass = -1
integer :: id_rescale_rho0_mask= -1
integer :: id_ht_mod = -1
logical :: used
!for restart file
integer :: id_restart_rho_dzt = 0
type(restart_file_type), save :: Thk_restart
! for discretization of time tendency ("threelevel" or "twolevel")
integer :: tendency
! conversion for prinout
real :: convert_factor=1.0
! for setting dst with pressure-based vertical coordinate models
real, allocatable, dimension(:) :: rho0_profile
! for the basin mask used to rescale rho0
real, allocatable, dimension(:,:) :: rescale_rho0_mask
real, allocatable, dimension(:,:) :: data
character(len=128) :: version=&
'$Id: ocean_thickness.F90,v 1.1.2.2.2.53.6.5.6.1.18.5.4.1 2009/10/10 00:42:01 nnz Exp $'
character (len=128) :: tagname = &
'$Name: mom4p1_pubrel_dec2009_nnz $'
type(ocean_domain_type), pointer :: Dom =>NULL()
public ocean_thickness_init
public ocean_thickness_init_adjust
public ocean_thickness_end
public update_tcell_thickness
public update_ucell_thickness
public rho_dzt_tendency
public ocean_thickness_restart
private thickness_restart
private thickness_initialize
private dst_land_adjust
private thickness_chksum
private thickness_details
private REMAP_ZT_TO_ZU
logical :: module_is_initialized = .false.
logical :: rescale_mass_to_get_ht_mod = .false.
character(len=32) :: thickness_method='energetic'
logical :: linear_free_surface = .false.
logical :: full_step_topography = .false.
logical :: enforce_positive_dzt = .false.
logical :: debug_this_module = .false.
logical :: debug_this_module_detail = .false.
logical :: write_a_restart = .true.
logical :: read_rho0_profile = .false.
logical :: pbot0_simple = .false.
logical :: initialize_zero_eta = .false.
logical :: read_rescale_rho0_mask = .false.
logical :: rescale_rho0_mask_gfdl = .false.
logical :: update_dzwu_k0 = .true.
integer :: max_num_bad_print = 25
real :: rescale_rho0_value = 1.0
real :: rescale_rho0_basin_label = -1.0
real :: thickness_dzt_min_init = 5.0 ! metre
real :: thickness_dzt_min = 1.0 ! metre
real :: depth_min_for_sigma = 0.01 ! metre
real :: epsilon_init_thickness = 1.e-5
namelist /ocean_thickness_nml/ debug_this_module, debug_this_module_detail, write_a_restart, &
full_step_topography, initialize_zero_eta, enforce_positive_dzt, &
depth_min_for_sigma, thickness_method, read_rho0_profile, &
thickness_dzt_min, thickness_dzt_min_init, &
rescale_mass_to_get_ht_mod, pbot0_simple, epsilon_init_thickness, &
read_rescale_rho0_mask, rescale_rho0_mask_gfdl, &
rescale_rho0_basin_label, rescale_rho0_value, linear_free_surface,&
max_num_bad_print, update_dzwu_k0
contains
!#######################################################################
!
!
!
! Initialize the thickness type.
!
! For pressure-based vertical coordinates, this initialization here
! is preliminary.
!
!
subroutine ocean_thickness_init (Time, Time_steps, Domain, Grid, Ext_mode, Thickness, &
ver_coordinate, ver_coordinate_class, ver_coordinate_type, &
debug)
type(ocean_time_type), intent(in) :: Time
type(ocean_time_steps_type), intent(in) :: Time_steps
type(ocean_domain_type), intent(in), target :: Domain
type(ocean_grid_type), intent(inout) :: Grid
type(ocean_external_mode_type), intent(in) :: Ext_mode
type(ocean_thickness_type), intent(inout) :: Thickness
integer, intent(in) :: ver_coordinate
integer, intent(in) :: ver_coordinate_class
integer, intent(in) :: ver_coordinate_type
logical, intent(in), optional :: debug
integer :: ioun, io_status, ierr
integer :: i,j,k
logical :: error_flag=.false.
integer :: stdoutunit,stdlogunit
stdoutunit=stdout();stdlogunit=stdlog()
if ( module_is_initialized ) then
call mpp_error(FATAL, &
'==>Error ocean_thickness_init: module has been initialized')
endif
module_is_initialized = .TRUE.
call write_version_number(version, tagname)
! provide for namelist over-ride of defaults
ioun = open_namelist_file()
read (ioun,ocean_thickness_nml,IOSTAT=io_status)
write (stdoutunit,'(/)')
write (stdoutunit,ocean_thickness_nml)
write (stdlogunit,ocean_thickness_nml)
ierr = check_nml_error(io_status, 'ocean_thickness_nml')
call close_file(ioun)
if (PRESENT(debug) .and. .not. debug_this_module) then
debug_this_module = debug
endif
if(debug_this_module) then
write(stdoutunit,'(a)') &
'==>Note: running ocean_thickness_mod with debug_this_module=.true.'
endif
if(debug_this_module_detail) then
write(stdoutunit,'(a)') &
'==>Note: running ocean_thickness_mod with debug_this_module_detail=.true.'
endif
if(.not. write_a_restart) then
write(stdoutunit,'(a)') '==>Note: running ocean_thickness with write_a_restart=.false.'
write(stdoutunit,'(a)') ' Will NOT write restart file, and so cannot restart the run.'
endif
grav_r = 1.0/grav
grav_rho0 = grav*rho0
tendency = Time_steps%tendency
vert_coordinate = ver_coordinate
vert_coordinate_class = ver_coordinate_class
vert_coordinate_type = ver_coordinate_type
if(vert_coordinate_class==DEPTH_BASED) then
convert_factor=1.0
elseif(vert_coordinate_class==PRESSURE_BASED) then
convert_factor=c2dbars
endif
Dom => Domain
halo = Dom%xhalo ! (xhalo=yhalo assumed in mom4)
if(enforce_positive_dzt) then
write(stdoutunit,'(/a)') '==>Warning: running with enforce_positive_dzt=.true. '
write(stdoutunit,'(a)') ' This option artifically truncates cell size to enforce dzt > 0.'
write(stdoutunit,'(a,f10.4/)') ' Minimum dzt is set to ', thickness_dzt_min
endif
if(vert_coordinate == ZSIGMA .or. vert_coordinate == PSIGMA) then
write(stdoutunit,'(a,f12.4,a)') &
'==>To remove division by zero with sigma-coord, add fictitious water of depth(m) ' &
,depth_min_for_sigma, ' over land. This water does not contribute to budgets.'
if(Grid%tripolar) then
call mpp_error(WARNING, &
'==>Warning ocean_thickness_mod: zsigma and psigma may have problems with tripolar. Testing incomplete.')
endif
endif
if(vert_coordinate == GEOPOTENTIAL .and. linear_free_surface) then
write(stdoutunit,'(a)') &
'==>Running GEOPOTENTIAL model with linear_free_surface, so cell thicknesses are constant in time.'
call mpp_error(WARNING, &
'==>Running GEOPOTENTIAL model with linear_free_surface, so cell thicknesses are constant in time.')
endif
if(thickness_method=='energetic') then
Thickness%method = ENERGETIC
write(stdoutunit,'(a)') &
'==>Note: running ocean_thickness with thickness_method=energetic.'
elseif(thickness_method=='finitevolume') then
Thickness%method = FINITEVOLUME
write(stdoutunit,'(a)') &
'==>Warning: ocean_thickness w/ thickness_method=finitevolume is experimental. Testing incomplete.'
call mpp_error(WARNING, &
'==>Warning ocean_thickness w/ thickness_method=finitevolume is experimental. Testing incomplete.')
endif
#ifndef MOM4_STATIC_ARRAYS
call get_local_indices(Domain, isd, ied, jsd, jed, isc, iec, jsc, jec)
nk = Grid%nk
allocate (Thickness%rho_dzt(isd:ied,jsd:jed,nk,3) )
allocate (Thickness%rho_dztr(isd:ied,jsd:jed,nk) )
allocate (Thickness%rho_dzu(isd:ied,jsd:jed,nk,3) )
allocate (Thickness%rho_dzur(isd:ied,jsd:jed,nk) )
allocate (Thickness%rho_dzt_tendency(isd:ied,jsd:jed,nk) )
allocate (Thickness%sea_lev(isd:ied,jsd:jed) )
allocate (Thickness%dzt(isd:ied,jsd:jed,nk) )
allocate (Thickness%dzu(isd:ied,jsd:jed,nk) )
allocate (Thickness%dzwt(isd:ied,jsd:jed,0:nk) )
allocate (Thickness%dzwu(isd:ied,jsd:jed,0:nk) )
allocate (Thickness%dztlo(isd:ied,jsd:jed,nk) )
allocate (Thickness%dztup(isd:ied,jsd:jed,nk) )
allocate (Thickness%geodepth_zt(isd:ied,jsd:jed,nk) )
allocate (Thickness%depth_zt(isd:ied,jsd:jed,nk) )
allocate (Thickness%depth_zwt(isd:ied,jsd:jed,nk) )
allocate (Thickness%depth_zu(isd:ied,jsd:jed,nk) )
allocate (Thickness%depth_zwu(isd:ied,jsd:jed,nk) )
allocate (Thickness%depth_st(isd:ied,jsd:jed,nk) )
allocate (Thickness%depth_swt(isd:ied,jsd:jed,nk) )
allocate (Thickness%dst(isd:ied,jsd:jed,nk) )
allocate (Thickness%dswt(isd:ied,jsd:jed,0:nk) )
allocate (Thickness%dzt_dst(isd:ied,jsd:jed,nk) )
allocate (Thickness%dstlo(isd:ied,jsd:jed,nk) )
allocate (Thickness%dstup(isd:ied,jsd:jed,nk) )
allocate (Thickness%pbot0(isd:ied,jsd:jed))
allocate (Thickness%pbot0r(isd:ied,jsd:jed))
allocate (Thickness%mass_source(isd:ied,jsd:jed,nk) )
allocate (Thickness%mass_u(isd:ied,jsd:jed,3) )
allocate (Thickness%thicku(isd:ied,jsd:jed,3) )
#endif
! set rescale_rho0_mask=rescale_rho0_value in those regions where we modify
! the rho0 value by a fraction. Otherwise rho0 is unscaled.
allocate (rescale_rho0_mask(isd:ied,jsd:jed) )
rescale_rho0_mask(:,:) = 1.0
if(read_rescale_rho0_mask .and. vert_coordinate_class==PRESSURE_BASED) then
allocate (data(isd:ied,jsd:jed) )
data(:,:) = 0.0
call read_data('INPUT/basin_mask','basin_mask',data,Domain%domain2d)
! the following is specific to the mask used at GFDL
if(rescale_rho0_mask_gfdl) then
write(stdoutunit,'(a,f12.6)') &
'==>Note: running ocean_thickness with rho0 in selected basin modified to rho0(kg/m3) = ',&
rho0*rescale_rho0_value
do j=jsc,jec
do i=isc,iec
if(data(i,j)==rescale_rho0_basin_label) then
rescale_rho0_mask(i,j)=rescale_rho0_value
else
rescale_rho0_mask(i,j)=1.0
endif
enddo
enddo
endif
call mpp_update_domains(rescale_rho0_mask(:,:), Dom%domain2d)
endif
id_rescale_rho0_mask= register_static_field ('ocean_model', 'rescale_rho0_mask', &
Grid%tracer_axes(1:2),'fraction of rho0 used for specifying pressure and pstar levels',&
'dimensionles',missing_value=missing_value, range=(/-1.e1,1.e8/))
if (id_rescale_rho0_mask > 0) then
used = send_data (id_rescale_rho0_mask, rescale_rho0_mask(isc:iec,jsc:jec), &
Time%model_time, rmask=Grid%tmask(isc:iec,jsc:jec,1))
endif
! rho0_profile for setting dst with pressure based vertical coordinate models
allocate(rho0_profile(nk))
rho0_profile(:) = rho0
if(read_rho0_profile) then
write(stdoutunit,'(a)') &
'==>Warning: ocean_thickness_mod: read rho0 profile to set dst w/ pressure models.'
write(stdoutunit,'(a)') &
' This option is experimental, and NOT generally supported.'
if(vert_coordinate_class==DEPTH_BASED) then
write(stdoutunit,'(a)') &
' Since using DEPTH_BASED vertical coordinates, rho0_profile = rho0.'
endif
if(vert_coordinate_class==PRESSURE_BASED) then
call read_data('INPUT/rho0_profile.nc','rho0_profile', rho0_profile)
write(stdoutunit,'(/a/)') 'rho0_profile is used to define pressure grid and pressure gradients.'
do k=1,nk
write(stdoutunit,'(a,i4,a,e22.12)') 'rho0_profile(',k,') = ',rho0_profile(k)
if(rho0_profile(k) <= 0.0) then
call mpp_error(FATAL, &
'==>ocean_thickness_mod: rho0_profile must be > 0.0, even in rock, since we divide by rho0_profile.')
endif
enddo
endif
endif
if(pbot0_simple) then
write(stdoutunit,'(/a)') '==>Warning from ocean_thickness_mod: pbot0_simple=.true.'
endif
! initialize vertical grid arrays
call thickness_initialize(Grid, Thickness)
! modify time dependent vertical grid arrays based on restart values
call thickness_restart(Time, Grid, Ext_mode, Thickness)
if(debug_this_module) then
write(stdoutunit,*) ' '
write(stdoutunit,*) ' From ocean_thickness_mod: initial thickness checksum'
call write_timestamp(Time%model_time)
call thickness_chksum(Time, Grid, Thickness)
endif
! check that ht is deep enough to maintain chosen minimum size cell.
error_flag=.false.
do j=jsc,jec
do i=isc,iec
if(Grid%ht(i,j) > 0.0 .and. Grid%ht(i,j) < thickness_dzt_min) then
write(unit,'(/a,i4,a,i4,e22.12,a,e22.12)') &
'==>Error: ocean_thickness_init: ht(',i+Dom%ioff,',',j+Dom%joff,') = ',Grid%ht(i,j), &
'is less than the chosen setting for thickness_dzt_min = ',thickness_dzt_min
error_flag=.true.
endif
enddo
enddo
if(error_flag) then
call mpp_error(FATAL,&
'==>ocean_thickness_init: ocean is not deep enough to support setting for thickness_dzt_min')
endif
! check that ht is deep enough to maintain chosen minimum size initial bottom cell.
if(vert_coordinate_class==PRESSURE_BASED) then
error_flag=.false.
do j=jsc,jec
do i=isc,iec
if(Grid%ht(i,j) > 0.0 .and. Grid%ht(i,j) < thickness_dzt_min_init) then
write(unit,'(/a,i4,a,i4,e22.12,a,e22.12)') &
'==>Error: ocean_thickness_init: ht(',i+Dom%ioff,',',j+Dom%joff,') = ',Grid%ht(i,j), &
'is less than the chosen setting for thickness_dzt_min_init = ',thickness_dzt_min_init
error_flag=.true.
endif
enddo
enddo
if(error_flag) then
call mpp_error(FATAL,&
'==>ocean_thickness_init: ocean is not deep enough to support setting for thickness_dzt_min_init')
endif
endif
! ids for diagnostic manager
id_dzu = register_diag_field ('ocean_model', 'dzu', Grid%vel_axes_uv(1:3), Time%model_time, &
'ocean u-cell thickness', 'm', &
missing_value=missing_value, range=(/-1e1,1e5/))
id_dzwu = register_diag_field ('ocean_model', 'dzwu', Grid%vel_axes_uv(1:3), Time%model_time, &
'thickness between velocity points', 'm', &
missing_value=missing_value, range=(/-1e1,1e5/))
id_dzt = register_diag_field ('ocean_model', 'dzt', Grid%tracer_axes(1:3), Time%model_time, &
'ocean t-cell thickness', 'm', &
missing_value=missing_value, range=(/-1e1,1e5/), &
standard_name='cell_thickness')
id_dzt_pbc = register_diag_field ('ocean_model', 'dzt_pbc', Grid%tracer_axes(1:2), Time%model_time, &
'ocean bottom t-cell thickness', 'm', &
missing_value=missing_value, range=(/-1e1,1e5/))
id_dztlo = register_diag_field ('ocean_model', 'dztlo', Grid%tracer_axes(1:3), Time%model_time, &
'ocean distance from t-cell center to t-cell bottom', 'm', &
missing_value=missing_value, range=(/-1e1,1e5/))
id_dztup = register_diag_field ('ocean_model', 'dztup', Grid%tracer_axes(1:3), Time%model_time, &
'ocean distance from t-cell center to t-cell top', 'm', &
missing_value=missing_value, range=(/-1e1,1e5/))
id_dstlo = register_diag_field ('ocean_model', 'dstlo', Grid%tracer_axes(1:3), Time%model_time, &
's-distance from t-cell center to t-cell bottom', 's-units', &
missing_value=missing_value, range=(/-1e5,1e5/))
id_dstup = register_diag_field ('ocean_model', 'dstup', Grid%tracer_axes(1:3), Time%model_time, &
's-distance from t-cell center to t-cell top', 's-units', &
missing_value=missing_value, range=(/-1e5,1e5/))
id_dzwt = register_diag_field ('ocean_model', 'dzwt', Grid%tracer_axes(1:3), Time%model_time, &
'thickness between tracer points', 'm', &
missing_value=missing_value, range=(/-1e1,1e5/))
id_dst = register_diag_field ('ocean_model', 'dst', Grid%tracer_axes(1:3), Time%model_time, &
'ocean s-cell thickness', 's-units', &
missing_value=missing_value, range=(/-1e5,1e5/))
id_geodepth_zt = register_diag_field ('ocean_model', 'geodepth_zt', Grid%tracer_axes(1:3), Time%model_time, &
'ocean t-cell depth relative to z=0', 'm', &
missing_value=missing_value, range=(/-1e1,1.e8/))
if(vert_coordinate == GEOPOTENTIAL) then
id_depth_zt = register_diag_field ('ocean_model', 'depth_zt', Grid%tracer_axes(1:3), Time%model_time, &
'ocean t-cell depth relative to z=0', 'm', &
missing_value=missing_value, range=(/-1e1,1.e8/))
id_depth_st = register_diag_field ('ocean_model', 'depth_st', Grid%tracer_axes(1:3), Time%model_time, &
'ocean t-cell depth relative to z=0', 'm', &
missing_value=missing_value, range=(/-1e1,1.e10/))
id_depth_swt = register_diag_field ('ocean_model', 'depth_swt', Grid%tracer_axes(1:3), Time%model_time, &
'ocean t-cell bottom relative to z=0', 'm', &
missing_value=missing_value, range=(/-1e1,1.e10/))
id_depth_zwt = register_diag_field ('ocean_model', 'depth_zwt', Grid%tracer_axes(1:3), Time%model_time, &
'depth to t-cell bottom relative to z=0', 'm', &
missing_value=missing_value, range=(/-1e1,1.e10/))
id_depth_zu = register_diag_field ('ocean_model', 'depth_zu', Grid%vel_axes_uv(1:3), Time%model_time, &
'ocean u-cell depth relative to z=0', 'm', &
missing_value=missing_value, range=(/-1e1,1.e10/))
id_depth_zwu = register_diag_field ('ocean_model', 'depth_zwu', Grid%vel_axes_uv(1:3), Time%model_time, &
'depth to ocean u-cell bottom reative to z=0', 'm', &
missing_value=missing_value, range=(/-1e1,1.e10/))
else
id_depth_zt = register_diag_field ('ocean_model', 'depth_zt', Grid%tracer_axes(1:3), Time%model_time, &
'ocean t-cell depth relative to column top', 'm', &
missing_value=missing_value, range=(/-1e1,1.e10/))
id_depth_st = register_diag_field ('ocean_model', 'depth_st', Grid%tracer_axes(1:3), Time%model_time, &
'ocean t-cell s-depth', 's-units', &
missing_value=missing_value, range=(/-1e8,1.e8/))
id_depth_swt = register_diag_field ('ocean_model', 'depth_swt', Grid%tracer_axes(1:3), Time%model_time, &
'ocean t-cell bottom s-depth', 's-units', &
missing_value=missing_value, range=(/-1e8,1.e8/))
id_depth_zwt = register_diag_field ('ocean_model', 'depth_zwt', Grid%tracer_axes(1:3), Time%model_time, &
'depth to t-cell bottom relative to column top', 'm', &
missing_value=missing_value, range=(/-1e1,1.e10/))
id_depth_zu = register_diag_field ('ocean_model', 'depth_zu', Grid%vel_axes_uv(1:3), Time%model_time, &
'ocean u-cell depth relative to column top', 'm', &
missing_value=missing_value, range=(/-1e1,1.e10/))
id_depth_zwu = register_diag_field ('ocean_model', 'depth_zwu', Grid%vel_axes_uv(1:3), Time%model_time, &
'depth to ocean u-cell bottom reative to column top', 'm', &
missing_value=missing_value, range=(/-1e1,1.e10/))
endif
id_dzt_dst = register_diag_field ('ocean_model', 'dzt_dst', Grid%tracer_axes(1:3), Time%model_time, &
'ocean t-cell specific thickness', 'm/s-coordinate', &
missing_value=missing_value, range=(/-1.e8,1.e8/))
id_rho_dzu = register_diag_field ('ocean_model', 'rho_dzu', Grid%vel_axes_uv(1:3), Time%model_time, &
'ocean u-cell rho*thickness', '(kg/m^3)*m', &
missing_value=missing_value, range=(/-1.e8,1.e8/))
id_rho_dzt = register_diag_field ('ocean_model', 'rho_dzt', Grid%tracer_axes(1:3), Time%model_time, &
'ocean t-cell rho*thickness', '(kg/m^3)*m', &
missing_value=missing_value, range=(/-1.e8,1.e8/), &
standard_name='sea_water_mass_per_unit_area')
id_mass_u = register_diag_field ('ocean_model', 'mass_u', Grid%vel_axes_uv(1:2), Time%model_time, &
'mass per area on C-cell column', 'kg/m^2', &
missing_value=missing_value, range=(/-1e1,1e8/))
id_thicku = register_diag_field ('ocean_model', 'thicku', Grid%vel_axes_uv(1:2), Time%model_time, &
'thickness from z=eta_u to z=-H of U-cell column', 'm', &
missing_value=missing_value, range=(/-1e1,1e8/))
id_rho_dzt_tendency = register_diag_field ('ocean_model', 'rho_dzt_tendency', Grid%tracer_axes(1:3), &
Time%model_time, 'rho_dzt_tendency', '(kg/m^3)*(m/s)', &
missing_value=missing_value, range=(/-1e8,1e8/))
id_pbot0 = register_static_field ('ocean_model', 'pbot0', Grid%tracer_axes(1:2), &
'reference bottom pressure t-cells', 'dbar', &
missing_value=missing_value, range=(/-1.e8,1.e8/))
end subroutine ocean_thickness_init
! NAME="ocean_thickness_init"
!#######################################################################
!
!
!
! Initialize vertical thicknesses of grid cells.
! For Boussinesq models, this code assumes the
! surface heights eta_t and eta_u are zero.
!
! The values here are relevant for the time=0 initialization
! of the model. Some time independent arrays are also set here,
! but they are over-written if there is a restart file.
!
! For pressure based vertical coordinates, the results here
! assume density = rho0_profile(k). The values are readjusted in
! ocean_thickness_init_adjust after we have determined the initial
! in situ density. This readjustment may involve enforcing an
! initial eta field of zero value, or the default which allows
! for eta to initially be nonzero.
!
!
!
subroutine thickness_initialize(Grid, Thickness)
type(ocean_grid_type), intent(inout) :: Grid
type(ocean_thickness_type), intent(inout) :: Thickness
real :: sigma_sign, density_tmp
integer :: i,j,k,n,kb,kmin
! initialization to zero
Thickness%sea_lev(:,:) = 0.0
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzt_tendency(i,j,k) = 0.0
Thickness%mass_source(i,j,k) = 0.0
enddo
enddo
enddo
! vertical increments for non-terrain following coordinates
if(vert_coordinate_type /= TERRAIN_FOLLOWING) then
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,0) = Grid%dzw(0)
Thickness%dzwu(i,j,0) = Grid%dzw(0)
enddo
enddo
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%dzt(i,j,k) = Grid%dzt(k)
Thickness%depth_zt(i,j,k) = Grid%zt(k)
Thickness%dzwt(i,j,k) = Grid%dzw(k)
Thickness%depth_zwt(i,j,k) = Grid%zw(k)
enddo
enddo
enddo
! modifications for partial step representation of bottom topography
if(full_step_topography) then
do j=jsd,jed
do i=isd,ied
kb = Grid%kmt(i,j)
if (kb > 1) then
Grid%ht(i,j) = Thickness%depth_zwt(i,j,kb)
endif
enddo
enddo
else
do j=jsd,jed
do i=isd,ied
kb = Grid%kmt(i,j)
if (kb > 1) then
Thickness%depth_zwt(i,j,kb)= Grid%ht(i,j)
Thickness%dzt(i,j,kb) = &
Thickness%depth_zwt(i,j,kb) - Thickness%depth_zwt(i,j,kb-1)
Thickness%depth_zt(i,j,kb) = &
Thickness%depth_zwt(i,j,kb-1) + Grid%fracdz(kb,0)*Thickness%dzt(i,j,kb)
Thickness%dzwt(i,j,kb-1) = &
Thickness%depth_zt(i,j,kb) - Thickness%depth_zt(i,j,kb-1)
Thickness%dzwt(i,j,kb) = &
Thickness%depth_zwt(i,j,kb) - Thickness%depth_zt(i,j,kb)
endif
enddo
enddo
endif
! vertical increments for terrain following coordinates
else
! -1 <= ZSIGMA <= 0, so dst and dsw are positive
! 0 <= PSIGMA <= 1, so dst and dsw are negative
! Recall s-coordinate is dimensionless when use ZSIGMA or PSIGMA.
if(vert_coordinate==ZSIGMA) then
sigma_sign =-1.0
elseif(vert_coordinate==PSIGMA) then
sigma_sign = 1.0
endif
! set depth_st and depth_swt according to Grid arrays
k=0
do j=jsd,jed
do i=isd,ied
Thickness%dswt(i,j,k) = Grid%dsw(k)
enddo
enddo
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_st(i,j,k) = Grid%st(k)
Thickness%depth_swt(i,j,k) = Grid%sw(k)
Thickness%dst(i,j,k) = Grid%dst(k)
Thickness%dswt(i,j,k) = Grid%dsw(k)
enddo
enddo
enddo
! now get the z-increments
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%dzt_dst(i,j,k) = -sigma_sign*max(Grid%ht(i,j),depth_min_for_sigma)
Thickness%dzt(i,j,k) = Thickness%dzt_dst(i,j,k)*Thickness%dst(i,j,k)
enddo
enddo
enddo
wrk1(:,:,:) = epsln
do k=1,nk
do j=jsd,jed
do i=isd,ied
if(Grid%tmask(i,j,k) > 0.0) then
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
endif
enddo
enddo
enddo
do k=1,nk-1
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,k) = 2.0*Thickness%dswt(i,j,k)/(wrk1(i,j,k)+wrk1(i,j,k+1))
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,0) = Thickness%dswt(i,j,0) *Thickness%dzt_dst(i,j,1)
Thickness%dzwt(i,j,nk) = Thickness%dswt(i,j,nk)*Thickness%dzt_dst(i,j,nk)
enddo
enddo
do j=jsd,jed
do i=isd,ied
Thickness%dzwu(i,j,0) = Thickness%dzwt(i,j,0)
enddo
enddo
! vertically sum the z-increments to get z-depths
k=1
do j=jsd,jed
do i=isd,ied
Thickness%depth_zt(i,j,k) = Thickness%dzwt(i,j,k-1)
Thickness%depth_zwt(i,j,k) = Thickness%dzt(i,j,k)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_zt(i,j,k) = Thickness%depth_zt(i,j,k-1) + Thickness%dzwt(i,j,k-1)
Thickness%depth_zwt(i,j,k) = Thickness%depth_zwt(i,j,k-1) + Thickness%dzt(i,j,k)
enddo
enddo
enddo
endif ! endif for ZSIGMA and PSIGMA
! half-thicknesses of T-cells
k=1
do j=jsd,jed
do i=isd,ied
Thickness%dztlo(i,j,k) = Thickness%depth_zwt(i,j,k)-Thickness%depth_zt(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dzwt(i,j,k-1)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%dztlo(i,j,k) = Thickness%depth_zwt(i,j,k)-Thickness%depth_zt(i,j,k)
Thickness%dztup(i,j,k) = Thickness%depth_zt(i,j,k) -Thickness%depth_zwt(i,j,k-1)
enddo
enddo
enddo
! U-cell thicknesses as minimum of surrounding T-cell thicknesses.
! This approach follows that used in all versions of MOM using
! thicknesses that are functions of (i,j,k).
Thickness%dzu = Thickness%dzt
Thickness%dzwu = Thickness%dzwt
do k=1,nk
do j=jsd,jec
do i=isd,iec
Thickness%dzu(i,j,k) = min(Thickness%dzt(i,j,k), &
Thickness%dzt(i+1,j,k), &
Thickness%dzt(i,j+1,k), &
Thickness%dzt(i+1,j+1,k))
Thickness%dzwu(i,j,k) = min(Thickness%dzwt(i,j,k), &
Thickness%dzwt(i+1,j,k), &
Thickness%dzwt(i,j+1,k), &
Thickness%dzwt(i+1,j+1,k))
enddo
enddo
enddo
if(update_dzwu_k0) then
k=0
do j=jsd,jec
do i=isd,iec
Thickness%dzwu(i,j,k) = min(Thickness%dzwt(i,j,k), &
Thickness%dzwt(i+1,j,k), &
Thickness%dzwt(i,j+1,k), &
Thickness%dzwt(i+1,j+1,k))
enddo
enddo
endif
! update across domain boundaries
do k=1,nk
call update_boundaries (Dom, Grid, Thickness%dzu(:,:,k), -1, -1,0,0)
enddo
do k=0,nk
call update_boundaries (Dom, Grid, Thickness%dzwu(:,:,k), -1, -1,0,0)
enddo
! U-cell depths based on U-cell thicknesses
Thickness%depth_zu = 0.0
Thickness%depth_zwu = 0.0
do j=jsd,jed
do i=isd,ied
Thickness%depth_zu(i,j,1) = Thickness%dzwu(i,j,0)
Thickness%depth_zwu(i,j,1) = Thickness%dzu(i,j,1)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_zu(i,j,k) = Thickness%depth_zu(i,j,k-1) + Thickness%dzwu(i,j,k-1)
Thickness%depth_zwu(i,j,k) = Thickness%depth_zwu(i,j,k-1) + Thickness%dzu(i,j,k)
enddo
enddo
enddo
! thickness of U-cell column
do n=1,3
do j=jsd,jed
do i=isd,ied
k=Grid%kmu(i,j)
if(k > 1) then
Thickness%thicku(i,j,n) = Thickness%depth_zwu(i,j,k)
else
Thickness%thicku(i,j,n) = 0.0
endif
enddo
enddo
enddo
! compute rho*dz values with rho0_profile*rescale_rho0_mask for initialization
Thickness%mass_u(:,:,:) = 0.0
do n=1,3
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzt(i,j,k,n) = rho0_profile(k)*rescale_rho0_mask(i,j)*Thickness%dzt(i,j,k)
Thickness%rho_dzu(i,j,k,n) = rho0_profile(k)*rescale_rho0_mask(i,j)*Thickness%dzu(i,j,k)
Thickness%mass_u(i,j,n) = Thickness%mass_u(i,j,n) &
+ Grid%umask(i,j,k)*Thickness%rho_dzu(i,j,k,n)
enddo
enddo
enddo
enddo
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,1)+epsln)
Thickness%rho_dzur(i,j,k) = 1.0/(Thickness%rho_dzu(i,j,k,1)+epsln)
enddo
enddo
enddo
! arrays for vertical s-coordinates;
! ZSIGMA and PSIGMA already computed these arrays
if(vert_coordinate==GEOPOTENTIAL .or. vert_coordinate==ZSTAR) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%dzt_dst(i,j,k) = 1.0
Thickness%dst(i,j,k) = Thickness%dzt(i,j,k)
Thickness%depth_st(i,j,k) = Thickness%depth_zt(i,j,k)
Thickness%depth_swt(i,j,k) = Thickness%depth_zwt(i,j,k)
enddo
enddo
enddo
do k=0,nk
do j=jsd,jed
do i=isd,ied
Thickness%dswt(i,j,k) = Thickness%dzwt(i,j,k)
enddo
enddo
enddo
elseif(vert_coordinate==PRESSURE .or. vert_coordinate==PSTAR) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
density_tmp = rho0_profile(k)*rescale_rho0_mask(i,j)
Thickness%dzt_dst(i,j,k) = -grav_r/(density_tmp+epsln)
Thickness%dst(i,j,k) = -grav*density_tmp*Thickness%dzt(i,j,k)
enddo
enddo
enddo
do k=0,nk
kmin = max(1,k)
do j=jsd,jed
do i=isd,ied
Thickness%dswt(i,j,k) = -grav*rho0_profile(kmin)*rescale_rho0_mask(i,j)*Thickness%dzwt(i,j,k)
enddo
enddo
enddo
! vertically sum the s-increments to get s-depths
k=1
do j=jsd,jed
do i=isd,ied
Thickness%depth_st(i,j,k) = -Thickness%dswt(i,j,k-1)
Thickness%depth_swt(i,j,k) = -Thickness%dst(i,j,k)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_st(i,j,k) = Thickness%depth_st(i,j,k-1) - Thickness%dswt(i,j,k-1)
Thickness%depth_swt(i,j,k) = Thickness%depth_swt(i,j,k-1) - Thickness%dst(i,j,k)
enddo
enddo
enddo
endif
! half-thicknesses of T-cells in s-space
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%dstlo(i,j,k) = Thickness%dztlo(i,j,k)/Thickness%dzt_dst(i,j,k)
Thickness%dstup(i,j,k) = Thickness%dztup(i,j,k)/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
! initialize depth from z=0 to t-point.
! this depth is used for computing the geopotential
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%geodepth_zt(i,j,k) = Thickness%depth_zt(i,j,k)
enddo
enddo
enddo
! tentative setting for reference bottom pressure (Pa)
! on T-cells, as well as its inverse. Will recompute
! this field in ocean_thickness_init_adjust
Thickness%pbot0(:,:) = 0.0
Thickness%pbot0r(:,:) = 0.0
if(pbot0_simple) then
do j=jsd,jed
do i=isd,ied
Thickness%pbot0(i,j) = Grid%tmask(i,j,1)*grav*rho0*Grid%ht(i,j)
enddo
enddo
else
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%pbot0(i,j) = Thickness%pbot0(i,j) + Grid%tmask(i,j,k)*grav*Thickness%rho_dzt(i,j,k,1)
enddo
enddo
enddo
endif
do j=jsd,jed
do i=isd,ied
if(Grid%kmt(i,j) > 1) then
Thickness%pbot0r(i,j) = 1.0/Thickness%pbot0(i,j)
endif
enddo
enddo
end subroutine thickness_initialize
! NAME="thickness_initialize"
!#######################################################################
!
!
!
!
! When initializing the pressure model, we can choose to initialize
! with a zero surface height, which requires some adjustment of
! bottom depths, or allow the surface height to be nonzero.
!
! The default is a nonzero surface height via
! initialize_zero_eta=.false.
!
!===============================================================
! The following comments refer to the possible issues arising
! from enforcing eta=0 on initialiation via
! initialize_zero_eta=.true.
!
! If we wish to enforce eta=0 on initialization, then we determine
! the depth of ocean needed to ensure a minimum bottom cell
! thickness as set by nml paramter thickness_dzt_min_init.
! This bottom depth is a function of the initial in-situ density
! as well as the nml thickness_dzt_min_init.
!
! If Grid%ht is sufficient according to the specifications of
! thickness_dzt_min_init and the initial in-situ density, then
! model is likely to remain stable, in that there is very small
! chance of losing bottom cell.
! (assuming thickness_dzt_min_init is > 5.0 or so).
!
! If Grid%ht is too shallow, again given the initial rho and
! pressure increments dst, then it is recommended that the user
! modify either the initial density (not easy) or the bottom
! topography (easy). The modifications are often quite trivial,
! depending on thickness_dzt_min_init and how light certain regions
! of the model are. The array "ht_mod" can be saved in the
! diag_table and compared to Grid%ht to see what modifications
! are required. If modificatations are to be made to the original
! grid_spec.nc file, simply use the NCO commands as follows:
!
! 1/ remove the original depth_t array from the original grid spec
! file, and write out a new grid spec file absent the depth_t array.
! ncea -x -v depth_t grid_spec_old.nc grid_spec_new.nc
!
! 2/ append to the new grid spec file the modified bottom topography
! array ht_mod, which was generated in an earlier run of the model.
! ncea -A -v ht_mod history_file_with_ht_mod.nc grid_spec_new.nc
!
! 3/ rename ht_mod to depth_t to conform to standard mom4p1 grid_spec
! name convention.
! ncvarrename new_grid.nc ht_mod depth_t
!
! If there are no modifications required to the depth_t array,
! then it is unlikely that the pressure model will evolve to the
! point of loosing the bottom cell and thus the model will remain
! stable. This conclusion is certainly based on the degree to which
! the model is run with changes in water masses.
!
!===============================================================
!
! NOTE: We generally compute values for the grid increments over
! land. The land values that are in the global halos are not
! going to be available from the restart files. Hence, they
! will need to be recomputed when restarting the model. This
! recalculation is done in dst_land_adjust.
!
! The land values are needed for use in computing U-cell grid
! factors using remapping operators. If we naively place zeros
! in the land, then, for example, rho_dzu next to land will be
! wrong, as it will be based on averaging a non-zero interior
! rho_dzt value with a rho_dzt incorretly set to zero, and so
! lead to spuriously small rho_dzu values next to land. Note that
! the placement of nonzero values of the grid values over land
! is something that is trivially done for the z-models. More care
! is needed with pressure models.
!
!
!
subroutine ocean_thickness_init_adjust(Grid, Time, Dens, Ext_mode, Thickness)
type(ocean_grid_type), intent(inout) :: Grid
type(ocean_time_type), intent(in) :: Time
type(ocean_density_type), intent(in) :: Dens
type(ocean_external_mode_type), intent(inout) :: Ext_mode
type(ocean_thickness_type), intent(inout) :: Thickness
logical :: warning_flag=.false.
integer :: i, j, k, kb, m, n
integer :: tau, taup1
real :: thick_tmp, density_tmp, delta_tmp
real :: max_fraction_differ, rescale_mass_min
real, dimension(isd:ied,jsd:jed) :: fraction_differ
real, dimension(isd:ied,jsd:jed) :: rescale_mass
real, dimension(isd:ied,jsd:jed) :: ht_mod
integer :: stdoutunit
stdoutunit=stdout()
if (file_exist('INPUT/ocean_thickness.res.nc')) then
if (id_pbot0 > 0) then
used = send_data (id_pbot0, c2dbars*Thickness%pbot0(isc:iec,jsc:jec), &
Time%model_time, rmask=Grid%tmask(isc:iec,jsc:jec,1))
endif
return
endif
tau = Time%tau
taup1 = Time%taup1
! for diagnostic purposes, compute initial bottom pressure
if(vert_coordinate_class==DEPTH_BASED .or. vert_coordinate_type==TERRAIN_FOLLOWING) then
if(.not. pbot0_simple) then
Thickness%pbot0(:,:) = 0.0
Thickness%pbot0r(:,:) = 0.0
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%pbot0(i,j) = Thickness%pbot0(i,j) &
+ Grid%tmask(i,j,k)*grav*Dens%rho(i,j,k,tau)*Thickness%dzt(i,j,k)
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
if(Grid%kmt(i,j) > 1) then
Thickness%pbot0r(i,j) = 1.0/Thickness%pbot0(i,j)
endif
enddo
enddo
endif
if (id_pbot0 > 0) then
used = send_data (id_pbot0, c2dbars*Thickness%pbot0(isc:iec,jsc:jec), &
Time%model_time, rmask=Grid%tmask(isc:iec,jsc:jec,1))
endif
return
endif
write(stdoutunit,'(/a/)') &
'==>Note: ocean_thickness_init_adjust adjusting time=0 vgrid using in situ density.'
wrk1(:,:,:) = 0.0
wrk2(:,:,:) = 0.0
wrk3(:,:,:) = 0.0
wrk4(:,:,:) = 0.0
! save original dzt computed from thickness_initialize
do k=1,nk
do j=jsd,jed
do i=isd,ied
wrk2(i,j,k) = Thickness%dzt(i,j,k)
wrk3(i,j,k) = Thickness%dztlo(i,j,k)
wrk4(i,j,k) = Thickness%dztup(i,j,k)
enddo
enddo
enddo
! fill arrays everywhere
do k=1,nk
do j=jsd,jed
do i=isd,ied
if(Grid%tmask(i,j,k) > 0.0) then
density_tmp = Dens%rho(i,j,k,tau)
else
density_tmp = rho0_profile(k)*rescale_rho0_mask(i,j)
endif
Thickness%dzt_dst(i,j,k) = -grav_r/(density_tmp+epsln)
Thickness%dzt(i,j,k) = Thickness%dst(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%rho_dzt(i,j,k,:) = density_tmp*Thickness%dzt(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,tau)+epsln)
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
rescale_mass(:,:) = 1.0
fraction_differ(:,:) = 0.0
ht_mod(:,:) = Grid%ht(:,:)
! code for enforcing eta=0 at initialization
if(initialize_zero_eta) then
! check to see if column is too thick to add
! a bottom cell of thickness thickness_dzt_min_init,
! assuming the surface height is initially eta_t=0.
do j=jsd,jed
do i=isd,ied
kb=Grid%kmt(i,j)
if(kb>1 .and. Grid%ht(i,j) > thickness_dzt_min_init) then
thick_tmp = thickness_dzt_min_init
do k=1,kb-1
thick_tmp = thick_tmp + Thickness%dzt(i,j,k)
enddo
if(thick_tmp >= (1.0+epsilon_init_thickness)*Grid%ht(i,j)) then
write(*,*)' '
write(*,'(a,i4,a,i4,a,f16.8,a,f16.8)') &
'==>ocean_thickness_init_adjust: thick_tmp(' &
,i+Dom%ioff,',',j+Dom%joff,') = ',thick_tmp, ' is >= Grid%ht = ',Grid%ht(i,j)
ht_mod(i,j) = thick_tmp
rescale_mass(i,j) = Grid%ht(i,j)/thick_tmp
fraction_differ(i,j) = (thick_tmp-Grid%ht(i,j))/Grid%ht(i,j)
endif
endif
enddo
enddo
max_fraction_differ = mpp_global_max(Dom%domain2d,fraction_differ)
write(stdoutunit,'(/a)') &
'==>ocean_thickness_init_adjust: WARNING ABOUT INITIAL TOPOGRAPHY AND DENSITY AND VERTICAL LEVELS'
write(stdoutunit,'(a)') &
'During initialization, sum(dzt[k=1,kmt-1]) resulted in columns that are too thick to fit a bottom cell with'
write(stdoutunit,'(a,e24.12)') &
'thickness equal to or greater than the namelist thickness_dzt_min_init = ',thickness_dzt_min_init
write(stdoutunit,'(a)') &
'This situation may be of no consequence to the stability of the pressure-based model simulation, if '
write(stdoutunit,'(a,e24.12,a)') &
'the maximum fractional difference ',max_fraction_differ,' is small, and thickness_dzt_min_init is greater than roughly 5m.'
write(stdoutunit,'(a)') &
'If these conditions are NOT met, there is the distinct possibility of losing bottom cells should density'
write(stdoutunit,'(a)') &
'evolve away from the initial conditions. A means for ensuring added model stability is to use the modified '
write(stdoutunit,'(a)') &
'bottom depths found by adding "ht_mod" to diag_table. These depths should replace "depth_t" in the original'
write(stdoutunit,'(a)') &
'grid specification file. Use the NCO command "ncea" to extract the original depth_t and replace with ht_mod.'
if(Time%init .and. rescale_mass_to_get_ht_mod) then
rescale_mass_min = mpp_global_min(Dom%domain2d,rescale_mass)
if(rescale_mass_min < 1.0) then
write(stdoutunit,'(a)') &
'==>Warning: rescale_mass_to_get_ht_mod=.true. Pressure gradients are wrong. Set false to run long.'
write(stdoutunit,'(a,e22.12)') &
'==>Warning: rescale_mass_min = ',rescale_mass_min
call mpp_error(WARNING, &
'==>ocean_thickness_init_adjust: rescale_mass_to_get_ht_mod=.true. Pressure gradients are wrong.')
endif
else
rescale_mass_min = 1.0
endif
do k=1,nk
rho0_profile(k) = rho0_profile(k)*rescale_mass_min
enddo
! endif for initialize_zero_eta
endif
! refill grid arrays over full column. reconsider bottom cell k=kmt later.
! NOTE: this step just reproduces what was done above the initialize_zero_eta
! if-block, in the case when rescale_mass_min=1.0
do k=1,nk
do j=jsd,jed
do i=isd,ied
density_tmp = rho0_profile(k)*rescale_rho0_mask(i,j)
Thickness%dst(i,j,k) = -grav*density_tmp*wrk2(i,j,k)
Thickness%dstlo(i,j,k) = -grav*density_tmp*wrk3(i,j,k)
Thickness%dstup(i,j,k) = -grav*density_tmp*wrk4(i,j,k)
if(Grid%tmask(i,j,k) > 0.0) then
density_tmp = Dens%rho(i,j,k,tau)
endif
Thickness%dzt_dst(i,j,k) = -grav_r/(density_tmp+epsln)
Thickness%dzt(i,j,k) = Thickness%dst(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%rho_dzt(i,j,k,:) = density_tmp*Thickness%dzt(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,tau)+epsln)
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
! Recompute depth_st and depth_swt by vertically summing s-increments.
! Note that this readjustment still maintains constant depth_st and constant
! depth_swt when k=constant, except at the bottom cell, since wrk2, wrk3, and
! wrk4 are only functions of k-level except at the bottom. We revisit the
! bottom cell issue below.
k=1
do j=jsd,jed
do i=isd,ied
Thickness%depth_st(i,j,k) = -Thickness%dstup(i,j,k)
Thickness%depth_swt(i,j,k) = -Thickness%dst(i,j,k)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_st(i,j,k) = Thickness%depth_swt(i,j,k-1) - Thickness%dstup(i,j,k)
Thickness%depth_swt(i,j,k) = Thickness%depth_st(i,j,k) - Thickness%dstlo(i,j,k)
enddo
enddo
enddo
! modifications for partial step representation of bottom topography
do j=jsd,jed
do i=isd,ied
kb=Grid%kmt(i,j)
if(kb>1) then
thick_tmp = 0.0
do k=1,kb-1
thick_tmp = thick_tmp + Thickness%dzt(i,j,k)
enddo
if(thick_tmp > Grid%ht(i,j) .and. initialize_zero_eta) then
warning_flag=.true.
write(unit,'(/a,i4,a,i4,a,i4, a, e22.12,a,i4,a,i4,a,e22.12)') &
'==>Warning: ocean_thickness_init_adjust. with kb = ', kb, ' thick_tmp(',i+Dom%ioff,',',&
j+Dom%joff,')= ',thick_tmp, ' is greater than ht(' ,i+Dom%ioff,',', &
j+Dom%joff,')= ',Grid%ht(i,j)
if(debug_this_module_detail) then
call thickness_details(Grid, Time, Thickness, Ext_mode, Dens, i, j, kb, &
'From ocean_thickness_init_adjust', unit)
endif
endif
k=kb
delta_tmp = Grid%ht(i,j) - thick_tmp
if(delta_tmp >= thickness_dzt_min_init) then
Thickness%dzt(i,j,k) = delta_tmp
else
Thickness%dzt(i,j,k) = thickness_dzt_min_init
endif
Thickness%dst(i,j,k) = -grav*Dens%rho(i,j,k,tau)*Thickness%dzt(i,j,k)
Thickness%rho_dzt(i,j,k,:) = Dens%rho(i,j,k,tau)*Thickness%dzt(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,tau)+epsln)
endif
enddo
enddo
if(debug_this_module_detail) then
call dzt_min_max(Time, Thickness, &
'From ocean_thickness_init_adjust, dzt_min_max information')
endif
if(warning_flag) then
call mpp_error(WARNING, &
'==>ocean_thickness_init_adjust: there potentially remains a problem with dzt.')
endif
! compute new reference bottom pressure
if(.not. pbot0_simple) then
Thickness%pbot0(:,:) = 0.0
Thickness%pbot0r(:,:) = 0.0
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%pbot0(i,j) = Thickness%pbot0(i,j) + Grid%tmask(i,j,k)*grav*Thickness%rho_dzt(i,j,k,1)
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
if(Grid%kmt(i,j) > 1) then
Thickness%pbot0r(i,j) = 1.0/Thickness%pbot0(i,j)
endif
enddo
enddo
endif
! adjust bottom cell s-dimensions (recall dst<0 and dswt<0 for pressure coordinates).
! do not alter the intermediate cell s-thicknesses, as these remain independent of i,j.
do j=jsd,jed
do i=isd,ied
kb=Grid%kmt(i,j)
if(kb>1) then
Thickness%depth_swt(i,j,kb) = Thickness%pbot0(i,j)
Thickness%depth_st(i,j,kb) = Thickness%depth_swt(i,j,kb-1) - Grid%fracdz(kb,0)*Thickness%dst(i,j,kb)
Thickness%dswt(i,j,kb-1) = -(Thickness%depth_st(i,j,kb) - Thickness%depth_st(i,j,kb-1))
Thickness%dswt(i,j,kb) = -(Thickness%depth_swt(i,j,kb) - Thickness%depth_st(i,j,kb))
Thickness%dstup(i,j,kb) = -(Thickness%depth_st(i,j,kb) -Thickness%depth_swt(i,j,kb-1))
Thickness%dstlo(i,j,kb) = -(Thickness%depth_swt(i,j,kb)-Thickness%depth_st(i,j,kb))
Thickness%dztup(i,j,kb) = Thickness%dstup(i,j,kb)*Thickness%dzt_dst(i,j,kb)
Thickness%dztlo(i,j,kb) = Thickness%dstlo(i,j,kb)*Thickness%dzt_dst(i,j,kb)
endif
enddo
enddo
! distances between tracer points
if(Thickness%method==ENERGETIC) then
do k=1,nk-1
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,k) = 2.0*Thickness%dswt(i,j,k)/(wrk1(i,j,k)+wrk1(i,j,k+1))
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,0) = Thickness%dswt(i,j,0)*Thickness%dzt_dst(i,j,1)
kb=Grid%kmt(i,j)
if(kb > 0) then
Thickness%dzwt(i,j,kb) = Thickness%dswt(i,j,kb)*Thickness%dzt_dst(i,j,kb)
endif
enddo
enddo
elseif(Thickness%method==FINITEVOLUME) then
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,0) = Thickness%dztup(i,j,1)
enddo
enddo
do k=1,nk-1
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k+1)
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
kb=Grid%kmt(i,j)
if(kb > 0) then
Thickness%dzwt(i,j,kb) = Thickness%dztlo(i,j,kb)
endif
enddo
enddo
endif
! adjust remaining z-dimensions throughout
! column through vertical sum of z-increments.
do j=jsd,jed
do i=isd,ied
Thickness%depth_zt(i,j,1) = Thickness%dzwt(i,j,0)
Thickness%depth_zwt(i,j,1) = Thickness%dzt(i,j,1)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_zt(i,j,k) = Thickness%depth_zt(i,j,k-1) + Thickness%dzwt(i,j,k-1)
Thickness%depth_zwt(i,j,k) = Thickness%depth_zwt(i,j,k-1) + Thickness%dzt(i,j,k)
enddo
enddo
enddo
! U-cell thicknesses computed as in update_ucell_thickness
do k=1,nk
do j=jsd,jec
do i=isd,iec
Thickness%dzu(i,j,k) = min(Thickness%dzt(i,j,k), &
Thickness%dzt(i+1,j,k), &
Thickness%dzt(i,j+1,k), &
Thickness%dzt(i+1,j+1,k))
Thickness%dzwu(i,j,k) = min(Thickness%dzwt(i,j,k), &
Thickness%dzwt(i+1,j,k), &
Thickness%dzwt(i,j+1,k), &
Thickness%dzwt(i+1,j+1,k))
enddo
enddo
enddo
if(update_dzwu_k0) then
k=0
do j=jsd,jec
do i=isd,iec
Thickness%dzwu(i,j,k) = min(Thickness%dzwt(i,j,k), &
Thickness%dzwt(i+1,j,k), &
Thickness%dzwt(i,j+1,k), &
Thickness%dzwt(i+1,j+1,k))
enddo
enddo
endif
do m=1,3
do k=1,nk
do j=jsd,jec
do i=isd,iec
Thickness%rho_dzu(i,j,k,m) = min(Thickness%rho_dzt(i,j,k,m), &
Thickness%rho_dzt(i+1,j,k,m), &
Thickness%rho_dzt(i,j+1,k,m), &
Thickness%rho_dzt(i+1,j+1,k,m))
enddo
enddo
enddo
enddo
call mpp_update_domains(Thickness%dzu(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%dzwu(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%rho_dzu(:,:,:,:), Dom%domain2d)
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzur(i,j,k) = Grid%umask(i,j,k)/(Thickness%rho_dzu(i,j,k,taup1)+epsln)
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
Thickness%depth_zu(i,j,1) = Thickness%dzwu(i,j,0)
Thickness%depth_zwu(i,j,1) = Thickness%dzu(i,j,1)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_zu(i,j,k) = Thickness%depth_zu(i,j,k-1) + Thickness%dzwu(i,j,k-1)
Thickness%depth_zwu(i,j,k) = Thickness%depth_zwu(i,j,k-1) + Thickness%dzu(i,j,k)
enddo
enddo
enddo
! thickness of U-cell column
do n=1,3
do j=jsd,jed
do i=isd,ied
k=Grid%kmu(i,j)
if(k > 1) then
Thickness%thicku(i,j,n) = Thickness%depth_zwu(i,j,k)
else
Thickness%thicku(i,j,n) = 0.0
endif
enddo
enddo
enddo
! mass per area of U-cell column
Thickness%mass_u(:,:,:) = 0.0
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%mass_u(i,j,1) = Thickness%mass_u(i,j,1) + Grid%umask(i,j,k)*Thickness%rho_dzu(i,j,k,tau)
Thickness%mass_u(i,j,2) = Thickness%mass_u(i,j,1)
Thickness%mass_u(i,j,3) = Thickness%mass_u(i,j,1)
enddo
enddo
enddo
! reinitialize the bottom pressure used for external mode calculation
do m=1,3
Ext_mode%pbot_t(:,:,m) = Thickness%pbot0(:,:)
Ext_mode%anompb(:,:,m) = Thickness%pbot0(:,:) - rho0*grav*Grid%ht(:,:)
Ext_mode%anompb_bar(:,:,m) = Thickness%pbot0(:,:) - rho0*grav*Grid%ht(:,:)
Ext_mode%pbot_u(:,:,m) = REMAP_ZT_TO_ZU(Thickness%pbot0(:,:),Grid)
enddo
id_rescale_mass= register_static_field ('ocean_model', 'rescale_mass', Grid%tracer_axes(1:2), &
'rescale_mass factor as function of Grid%ht and initial rho',&
'dimensionless', missing_value=missing_value, range=(/-1.e8,1.e8/))
if (id_rescale_mass > 0) used = send_data (id_rescale_mass, rescale_mass(isc:iec,jsc:jec), &
Time%model_time, rmask=Grid%tmask(isc:iec,jsc:jec,1))
id_ht_mod = register_static_field ('ocean_model', 'ht_mod', Grid%tracer_axes(1:2), &
'modified bottom depth for pressure model', 'm', &
missing_value=missing_value, range=(/-1.e1,1.e8/))
if (id_ht_mod > 0) used = send_data (id_ht_mod, ht_mod(isc:iec,jsc:jec), &
Time%model_time, rmask=Grid%tmask(isc:iec,jsc:jec,1))
if (id_pbot0 > 0) then
used = send_data (id_pbot0, c2dbars*Thickness%pbot0(isc:iec,jsc:jec), &
Time%model_time, rmask=Grid%tmask(isc:iec,jsc:jec,1))
endif
if(debug_this_module) then
write(stdoutunit,*) ' '
write(stdoutunit,*) ' From ocean_thickness_init_adjust: readjusted initial thickness checksum'
call write_timestamp(Time%model_time)
call thickness_chksum(Time, Grid, Thickness)
endif
if(debug_this_module_detail) then
call thickness_details(Grid, Time, Thickness, Ext_mode, Dens, isc, jsc, nk, &
'From ocean_thickness_init_adjust', stdoutunit)
call dzt_min_max(Time, Thickness,'From ocean_thickness_init_adjust, dzt_min_max information')
endif
if(debug_this_module_detail) then
write(unit,'(a)')''
write(unit,'(a)')'---------Bottom values for grid thicknesses---------'
do j=jsd,jed
do i=isd,ied
kb=Grid%kmt(i,j)
if(kb>1) then
if(Thickness%depth_swt(i,j,kb) > 0.0) then
write(unit,'(a)')''
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'depth_swt(kb-1)(',i+Dom%ioff,',',j+Dom%joff,',',kb-1,') = ',Thickness%depth_swt(i,j,kb-1)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'depth_swt(kb) (',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',Thickness%depth_swt(i,j,kb)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'depth_st(kb) (',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',Thickness%depth_st(i,j,kb)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dswt(kb) (',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',Thickness%dswt(i,j,kb)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dst(kb) (',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',Thickness%dst(i,j,kb)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dstup(kb) (',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',Thickness%dstup(i,j,kb)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dstlo(kb) (',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',Thickness%dstlo(i,j,kb)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dstlo+dstup(kb)(',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',&
Thickness%dstlo(i,j,kb)+Thickness%dstup(i,j,kb)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dzt(kb) (',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',Thickness%dzt(i,j,kb)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dztup(kb) (',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',Thickness%dztup(i,j,kb)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dztlo(kb) (',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',Thickness%dztlo(i,j,kb)
write(unit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dztlo+dztup(kb)(',i+Dom%ioff,',',j+Dom%joff,',',kb,') = ',&
Thickness%dztlo(i,j,kb)+Thickness%dztup(i,j,kb)
endif
endif
enddo
enddo
endif
end subroutine ocean_thickness_init_adjust
! NAME="ocean_thickness_init_adjust"
!#######################################################################
!
!
!
! Read basic elements of thickness derived type from restart,
! then set remaining elements of the type.
!
! Note that some array elements may be time independent.
! Their values have already been set in thickness_initialize
! or ocean_thickness_init_adjust, and they should not be
! over-written here (in particular, they should not be
! reset here to zero).
!
! To ensure reproducibility across restarts, this routine
! uses the same logic as update_tcell_thickness and
! update_ucell_thickness. This is particularly important
! for the terrain following calculations.
!
! Note that ocean_thickness_init is called prior to
! ocean_operators_init. Hence, we cannot use any
! operators in this subroutine.
!
!
!
subroutine thickness_restart (Time, Grid, Ext_mode, Thickness)
type(ocean_time_type), intent(in) :: Time
type(ocean_grid_type), intent(in) :: Grid
type(ocean_external_mode_type), intent(in) :: Ext_mode
type(ocean_thickness_type), intent(inout) :: Thickness
character*128 file_name
integer :: tau, taum1, taup1, id_restart
integer :: i,j,k,n,kb
real :: density_tmp
integer :: stdoutunit
stdoutunit=stdout()
tau = Time%tau
taum1 = Time%taum1
taup1 = Time%taup1
file_name = 'ocean_thickness.res.nc'
if(tendency==THREE_LEVEL) then
id_restart_rho_dzt = register_restart_field(Thk_restart, file_name, 'rho_dzt', Thickness%rho_dzt(:,:,:,tau), &
Thickness%rho_dzt(:,:,:,taup1), Dom%domain2d)
elseif(tendency==TWO_LEVEL) then
id_restart_rho_dzt = register_restart_field(Thk_restart, file_name, 'rho_dzt', &
Thickness%rho_dzt(:,:,:,taup1), Dom%domain2d)
end if
id_restart = register_restart_field(Thk_restart, file_name, 'dzt_dst', &
Thickness%dzt_dst(:,:,:), Dom%domain2d)
id_restart = register_restart_field(Thk_restart, file_name, 'dstlo', &
Thickness%dstlo(:,:,:), Dom%domain2d)
id_restart = register_restart_field(Thk_restart, file_name, 'dstup', &
Thickness%dstup(:,:,:), Dom%domain2d)
id_restart = register_restart_field(Thk_restart, file_name, 'dswt', &
Thickness%dswt(:,:,:), Dom%domain2d)
id_restart = register_restart_field(Thk_restart, file_name, 'dst', &
Thickness%dst(:,:,:), Dom%domain2d)
id_restart = register_restart_field(Thk_restart, file_name, 'pbot0', &
Thickness%pbot0(:,:), Dom%domain2d)
if(vert_coordinate_class==PRESSURE_BASED) then
id_restart = register_restart_field(Thk_restart, file_name, 'depth_st', &
Thickness%depth_st(:,:,:), Dom%domain2d)
id_restart = register_restart_field(Thk_restart, file_name, 'depth_swt', &
Thickness%depth_swt(:,:,:), Dom%domain2d)
end if
if (.NOT. file_exist('INPUT/ocean_thickness.res.nc')) then
if (.NOT.Time%init) then
call mpp_error(FATAL,&
'Expecting file INPUT/ocean_thickness.res.nc to exist.&
&This file was not found and Time%init=.false.')
endif
return
endif
call restore_state(Thk_restart)
if(tendency==THREE_LEVEL) then
write (stdoutunit,'(a)') ' Reading THREE_LEVEL restart for rho_dzt from INPUT/ocean_thickness.res.nc'
write (stdoutunit,'(a)') ' If not an initial condition, then expect two time records for rho_dzt.'
call mpp_update_domains(Thickness%rho_dzt(:,:,:,tau), Dom%domain2d)
call mpp_update_domains(Thickness%rho_dzt(:,:,:,taup1), Dom%domain2d)
elseif(tendency==TWO_LEVEL) then
write (stdoutunit,'(a)') ' Reading TWO_LEVEL restart for rho_dzt from INPUT/ocean_thickness.res.nc'
write (stdoutunit,'(a)') ' Expect only one time record for rho_dzt.'
! initialize to nonzero to eliminate NaNs when mask
! out processors and then change layout upon a restart
where (Thickness%rho_dzt(:,:,:,taup1)==0.) Thickness%rho_dzt(:,:,:,taup1) = rho0*thickness_dzt_min_init
call mpp_update_domains(Thickness%rho_dzt(:,:,:,taup1), Dom%domain2d)
Thickness%rho_dzt(:,:,:,tau) = Thickness%rho_dzt(:,:,:,taup1)
endif
! initialize to nonzero to eliminate NaNs when mask
! out processors and then change layout upon a restart
where (Thickness%dzt_dst(:,:,:) ==0.) Thickness%dzt_dst(:,:,:) = thickness_dzt_min_init
where (Thickness%dstlo(:,:,:) ==0.) Thickness%dstlo(:,:,:) = thickness_dzt_min_init
where (Thickness%dstup(:,:,:) ==0.) Thickness%dstup(:,:,:) = thickness_dzt_min_init
where (Thickness%dst(:,:,:) ==0.) Thickness%dst(:,:,:) = thickness_dzt_min_init
where (Thickness%dswt(:,:,:) ==0.) Thickness%dswt(:,:,:) = thickness_dzt_min_init
call mpp_update_domains(Thickness%dzt_dst(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%dstlo(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%dstup(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%dst(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%dswt(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%pbot0(:,:), Dom%domain2d)
if(vert_coordinate_class==PRESSURE_BASED) then
call mpp_update_domains(Thickness%depth_st(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%depth_swt(:,:,:), Dom%domain2d)
endif
! for adjusting values of pressure-coordinate intervals over land
call dst_land_adjust(Grid, Thickness)
! inverse pbot0 values
Thickness%pbot0r = 0.0
do j=jsd,jed
do i=isd,ied
if(Grid%kmt(i,j) > 1) then
Thickness%pbot0r(i,j) = 1.0/Thickness%pbot0(i,j)
endif
enddo
enddo
! T-cell dz values
if(vert_coordinate==GEOPOTENTIAL) then
do j=jsd,jed
do i=isd,ied
Thickness%dztlo(i,j,1) = Thickness%dstlo(i,j,1)
Thickness%dztup(i,j,1) = Thickness%dstup(i,j,1)
Thickness%dzt(i,j,1) = Thickness%dztlo(i,j,1) + Thickness%dztup(i,j,1)
Thickness%rho_dztr(i,j,1) = 1.0/(Thickness%rho_dzt(i,j,1,taup1)+epsln)
wrk1(i,j,:) = 1.0
enddo
enddo
elseif(vert_coordinate==ZSIGMA) then
! note that it is important to do the "max" step here,
! as the restart will not get the halos correct, so we must do
! this calculation here to get dzt_dst in the halos.
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%dzt_dst(i,j,k) = max(depth_min_for_sigma,Grid%ht(i,j)+Ext_mode%eta_t(i,j,taup1))
Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dzt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,taup1)+epsln)
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
elseif(vert_coordinate==PRESSURE .or. vert_coordinate==PSTAR) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
if(Grid%tmask(i,j,k)==0.0) then
density_tmp = rho0_profile(k)*rescale_rho0_mask(i,j)
Thickness%dzt_dst(i,j,k) = -grav_r/(density_tmp+epsln)
endif
Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dzt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,taup1)+epsln)
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
elseif(vert_coordinate==PSIGMA) then
! note that it is important to do the Grid%tmask(i,j,k)==0.0 steps,
! as the restart will not get the halos correct, so we must do
! this calculation here to get dzt_dst in the halos.
do k=1,nk
do j=jsd,jed
do i=isd,ied
if(Grid%tmask(i,j,k) == 0.0) then
Thickness%dzt_dst(i,j,k) = -depth_min_for_sigma
endif
Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dzt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k)
if(Grid%tmask(i,j,k) == 0.0) then
Thickness%rho_dzt(i,j,k,taup1) = rho0_profile(k)*Thickness%dzt(i,j,k)
endif
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,taup1)+epsln)
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
else
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dzt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,taup1)+epsln)
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
endif
! vertical distance (m) between t-cell points
if(Thickness%method==ENERGETIC) then
do k=1,nk-1
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,k) = 2.0*Thickness%dswt(i,j,k)/(wrk1(i,j,k)+wrk1(i,j,k+1))
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,0) = Thickness%dswt(i,j,0)*Thickness%dzt_dst(i,j,1)
kb=Grid%kmt(i,j)
if(kb > 0) then
Thickness%dzwt(i,j,kb) = Thickness%dswt(i,j,kb)*Thickness%dzt_dst(i,j,kb)
endif
enddo
enddo
elseif(Thickness%method==FINITEVOLUME) then
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,0) = Thickness%dztup(i,j,1)
enddo
enddo
do k=1,nk-1
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k+1)
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
kb=Grid%kmt(i,j)
if(kb > 0) then
Thickness%dzwt(i,j,kb) = Thickness%dztlo(i,j,kb)
endif
enddo
enddo
endif
! vertical depth (m) of t and w points
if(vert_coordinate /= GEOPOTENTIAL) then
! depth from z=eta to t-point and t-cell bottom
do j=jsd,jed
do i=isd,ied
Thickness%depth_zt(i,j,1) = Thickness%dzwt(i,j,0)
Thickness%depth_zwt(i,j,1) = Thickness%dzt(i,j,1)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_zt(i,j,k) = Thickness%depth_zt(i,j,k-1) + Thickness%dzwt(i,j,k-1)
Thickness%depth_zwt(i,j,k) = Thickness%depth_zwt(i,j,k-1) + Thickness%dzt(i,j,k)
enddo
enddo
enddo
! depth from z=0 to t point
do j=jsd,jed
do i=isd,ied
Thickness%geodepth_zt(i,j,1) = Thickness%depth_zt(i,j,1) - Ext_mode%eta_t(i,j,taup1)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%geodepth_zt(i,j,k) = Thickness%geodepth_zt(i,j,k-1) + Thickness%dzwt(i,j,k-1)
enddo
enddo
enddo
endif
! U-cell thicknesses computed as in update_ucell_thickness
if(vert_coordinate==GEOPOTENTIAL) then
do j=jsd,jed
do i=isd,ied
Thickness%dzwu(i,j,0) = Thickness%depth_zt(i,j,1) + Ext_mode%eta_u(i,j,taup1)
Thickness%dzu(i,j,1) = Thickness%depth_zwt(i,j,1) + Ext_mode%eta_u(i,j,taup1)
Thickness%rho_dzu(i,j,1,taup1) = rho0*Thickness%dzu(i,j,1)
Thickness%rho_dzu(i,j,1,tau) = rho0*(Thickness%depth_zwt(i,j,1) + Ext_mode%eta_u(i,j,tau))
Thickness%rho_dzur(i,j,1) = 1.0/(Thickness%rho_dzu(i,j,1,taup1)+epsln)
enddo
enddo
if(linear_free_surface) then
do j=jsd,jed
do i=isd,ied
Thickness%dzwu(i,j,0) = Thickness%depth_zt(i,j,1)
Thickness%dzu(i,j,1) = Thickness%depth_zwt(i,j,1)
Thickness%rho_dzu(i,j,1,taup1) = rho0*Thickness%dzu(i,j,1)
Thickness%rho_dzu(i,j,1,tau) = rho0*Thickness%depth_zwt(i,j,1)
Thickness%rho_dzur(i,j,1) = 1.0/(Thickness%rho_dzu(i,j,1,taup1)+epsln)
enddo
enddo
endif
else
if(vert_coordinate==ZSIGMA .or. vert_coordinate==PSIGMA) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%dzu(i,j,k) = Thickness%dzt(i,j,k)
Thickness%dzwu(i,j,k) = Thickness%dzwt(i,j,k)
Thickness%rho_dzu(i,j,k,taup1) = Thickness%rho_dzt(i,j,k,taup1)
enddo
enddo
enddo
k=0
do j=jsd,jed
do i=isd,ied
Thickness%dzwu(i,j,k) = Thickness%dzwt(i,j,k)
enddo
enddo
else
do k=1,nk
do j=jsd,jec
do i=isd,iec
Thickness%dzu(i,j,k) = min(Thickness%dzt(i,j,k), &
Thickness%dzt(i+1,j,k), &
Thickness%dzt(i,j+1,k), &
Thickness%dzt(i+1,j+1,k))
Thickness%dzwu(i,j,k) = min(Thickness%dzwt(i,j,k), &
Thickness%dzwt(i+1,j,k), &
Thickness%dzwt(i,j+1,k), &
Thickness%dzwt(i+1,j+1,k))
Thickness%rho_dzu(i,j,k,taup1) = min(Thickness%rho_dzt(i,j,k,taup1), &
Thickness%rho_dzt(i+1,j,k,taup1), &
Thickness%rho_dzt(i,j+1,k,taup1), &
Thickness%rho_dzt(i+1,j+1,k,taup1))
enddo
enddo
enddo
if(update_dzwu_k0) then
k=0
do j=jsd,jec
do i=isd,iec
Thickness%dzwu(i,j,k) = min(Thickness%dzwt(i,j,k), &
Thickness%dzwt(i+1,j,k), &
Thickness%dzwt(i,j+1,k), &
Thickness%dzwt(i+1,j+1,k))
enddo
enddo
endif
call mpp_update_domains(Thickness%dzu(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%dzwu(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%rho_dzu(:,:,:,taup1), Dom%domain2d)
endif
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzur(i,j,k) = 1.0/(Thickness%rho_dzu(i,j,k,taup1)+epsln)
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
Thickness%depth_zu(i,j,1) = Thickness%dzwu(i,j,0)
Thickness%depth_zwu(i,j,1) = Thickness%dzu(i,j,1)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_zu(i,j,k) = Thickness%depth_zu(i,j,k-1) + Thickness%dzwu(i,j,k-1)
Thickness%depth_zwu(i,j,k) = Thickness%depth_zwu(i,j,k-1) + Thickness%dzu(i,j,k)
enddo
enddo
enddo
endif ! endif for choice of vertical coordinate
! thickness of U-cell column
do n=1,3
do j=jsd,jed
do i=isd,ied
k=Grid%kmu(i,j)
if(k > 1) then
Thickness%thicku(i,j,n) = Thickness%depth_zwu(i,j,k)
else
Thickness%thicku(i,j,n) = 0.0
endif
enddo
enddo
enddo
! mass per area of U-cell column
Thickness%mass_u(:,:,:) = 0.0
do n=1,3
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%mass_u(i,j,n) = Thickness%mass_u(i,j,n) + Grid%umask(i,j,k)*Thickness%rho_dzu(i,j,k,n)
enddo
enddo
enddo
enddo
end subroutine thickness_restart
! NAME="thickness_restart"
!#######################################################################
!
!
!
! For computing the grid values inside land in global halos.
! These values are not zero, and so they are not available from
! restart file. They are needed for pressure-based vertical
! coordinate models in order to get the U-cell values at
! global land boundaries via the remap operator.
!
! This step is necessitated by the modifications to dst
! and other arrays made in ocean_thickness_init_adjust.
!
!
subroutine dst_land_adjust(Grid, Thickness)
type(ocean_grid_type), intent(in) :: Grid
type(ocean_thickness_type), intent(inout) :: Thickness
integer :: i, j, k
real :: density_tmp
if(vert_coordinate_class==DEPTH_BASED) return
if(vert_coordinate_type==TERRAIN_FOLLOWING) return
do k=1,nk
do j=jsd,jed
do i=isd,ied
density_tmp = rho0_profile(k)*rescale_rho0_mask(i,j)
if(Grid%tmask(i,j,k) == 0.0) then
Thickness%dst(i,j,k) = -grav*density_tmp*Thickness%dzt(i,j,k)
Thickness%dzt_dst(i,j,k) = -grav_r/(density_tmp+epsln)
Thickness%dzt(i,j,k) = Thickness%dst(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%rho_dzt(i,j,k,:) = density_tmp*Thickness%dzt(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,1)+epsln)
endif
enddo
enddo
enddo
end subroutine dst_land_adjust
! NAME="dst_land_adjust"
!#######################################################################
!
!
!
! Update time dependent thickness of T cells.
!
! Notes
!
! 1. For pressure-based coordinates, must use rho(tau) since
! rho(taup1) has not yet been computed. This is a limitation of
! the "z-like" algorithm approach used in mom4p1. It is minor,
! however, since rho(tau) is very close to rho(taup1). Also, this
! approach comes with the advantage of rendering vertical
! physical processes (e.g., diffusion) linear for the general
! coordinates in mom4p1, just as with geopotential models.
!
! 2. For all coordinates, need to place something reasonable over
! land for dz increments to preclude division by zero.
! For GEOPOTENTIAL and ZSTAR, we just use the time=zero value for dzt.
! For ZSIGMA we set a minimum depth for a fictitious layer over land.
! For PRESSURE and PSTAR, we set rho=rho0_profile and use the time=0 dzt values.
! For PSIGMA, we rho=rho0_profile and use the fictitious layer over land.
! The treatment of these land-points only contributes to the treatment
! of the U-cell grid quantities in the halos and on land, as they are
! computed via the REMAP_ZT_TO_ZU operator.
!
! 3. We need to update s-coordinate increments for GEOPOTENTIAL and
! PRESSURE. It is only for these two coordinates that we modify the
! endpoint values for the s-grid increments.
!
!
!
subroutine update_tcell_thickness (Time, Grid, Ext_mode, Dens, Thickness)
type(ocean_time_type), intent(in) :: Time
type(ocean_grid_type), intent(in) :: Grid
type(ocean_external_mode_type), intent(in) :: Ext_mode
type(ocean_density_type), intent(in) :: Dens
type(ocean_thickness_type), intent(inout) :: Thickness
integer :: i, j, k, kb
integer :: tau, taup1
integer :: num_bad=0
logical :: huge_undulations
logical :: error_flag=.false.
real :: density_tmp
character(len=128) :: errorstring
integer :: stdoutunit
stdoutunit=stdout()
tau = Time%tau
taup1 = Time%taup1
huge_undulations=.false.
! write tau values of the vertical grid increments before update to taup1
if (id_dst > 0) used = send_data (id_dst, Thickness%dst(:,:,:), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_dzt > 0) used = send_data (id_dzt, Thickness%dzt(:,:,:), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_dzt_pbc > 0) then
wrk1_2d(:,:) = 0.0
do j=jsd,jed
do i=isd,ied
kb = Grid%kmt(i,j)
if (kb > 1) then
wrk1_2d(i,j) = Thickness%dzt(i,j,kb)
endif
enddo
enddo
used = send_data (id_dzt_pbc, wrk1_2d(:,:), &
Time%model_time, rmask=Grid%tmask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
endif
if (id_dztlo > 0) used = send_data (id_dztlo, Thickness%dztlo(:,:,:), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_dztup > 0) used = send_data (id_dztup, Thickness%dztup(:,:,:), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_dzwt > 0) used = send_data (id_dzwt, Thickness%dzwt(:,:,1:nk), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_rho_dzt > 0) used = send_data (id_rho_dzt, Thickness%rho_dzt(:,:,:,tau), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_dzt_dst > 0) used = send_data (id_dzt_dst, Thickness%dzt_dst(:,:,:), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_depth_zt > 0) used = send_data (id_depth_zt, Thickness%depth_zt(:,:,:), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_geodepth_zt > 0) used = send_data (id_geodepth_zt, Thickness%geodepth_zt(:,:,:), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_depth_st > 0) used = send_data (id_depth_st, Thickness%depth_st(:,:,:), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_depth_swt > 0) used = send_data (id_depth_swt, Thickness%depth_swt(:,:,:),&
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_depth_zwt > 0) used= send_data (id_depth_zwt, Thickness%depth_zwt(:,:,:), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
! update coordinate increments for GEOPOTENTIAL and PRESSURE.
! only for these two coordinates do we need to modify the
! endpoint values for the s-grid increments. Other mom4p1
! coordinates have dst and dswt constant in time.
if(vert_coordinate==GEOPOTENTIAL) then
do j=jsd,jed
do i=isd,ied
Thickness%dswt(i,j,0) = Thickness%depth_zt(i,j,1) + Ext_mode%eta_t(i,j,taup1)
Thickness%dst(i,j,1) = Thickness%depth_zwt(i,j,1) + Ext_mode%eta_t(i,j,taup1)
Thickness%dstup(i,j,1) = Thickness%dswt(i,j,0)
Thickness%dstlo(i,j,1) = Thickness%dst(i,j,1) - Thickness%dstup(i,j,1)
enddo
enddo
if(linear_free_surface) then
do j=jsd,jed
do i=isd,ied
Thickness%dswt(i,j,0) = Thickness%depth_zt(i,j,1)
Thickness%dst(i,j,1) = Thickness%depth_zwt(i,j,1)
Thickness%dstup(i,j,1) = Thickness%dswt(i,j,0)
Thickness%dstlo(i,j,1) = Thickness%dst(i,j,1) - Thickness%dstup(i,j,1)
enddo
enddo
endif
elseif(vert_coordinate==PRESSURE) then
do j=jsd,jed
do i=isd,ied
Thickness%dswt(i,j,0) = -Thickness%depth_st(i,j,1) + Ext_mode%patm_t(i,j,taup1)
Thickness%dst(i,j,1) = -Thickness%depth_swt(i,j,1) + Ext_mode%patm_t(i,j,taup1)
Thickness%dstup(i,j,1) = Thickness%dswt(i,j,0)
Thickness%dstlo(i,j,1) = Thickness%dst(i,j,1) - Thickness%dstup(i,j,1)
kb = Grid%kmt(i,j)
if(kb > 1) then
Thickness%dswt(i,j,kb) = Thickness%depth_st(i,j,kb) - Ext_mode%pbot_t(i,j,taup1)
Thickness%dst(i,j,kb) = Thickness%depth_swt(i,j,kb-1) - Ext_mode%pbot_t(i,j,taup1)
Thickness%dstlo(i,j,kb) = Thickness%dswt(i,j,kb)
Thickness%dstup(i,j,kb) = Thickness%dst(i,j,kb) - Thickness%dstlo(i,j,kb)
endif
enddo
enddo
endif
! update specific thicknesses and vertical increments
if(vert_coordinate==GEOPOTENTIAL) then
do j=jsd,jed
do i=isd,ied
Thickness%dztlo(i,j,1) = Thickness%dstlo(i,j,1)
Thickness%dztup(i,j,1) = Thickness%dstup(i,j,1)
Thickness%dzt(i,j,1) = Thickness%dztlo(i,j,1) + Thickness%dztup(i,j,1)
Thickness%rho_dzt(i,j,1,taup1) = rho0*Thickness%dzt(i,j,1)
Thickness%rho_dztr(i,j,1) = 1.0/(Thickness%rho_dzt(i,j,1,taup1)+epsln)
wrk1(i,j,:) = 1.0
enddo
enddo
elseif(vert_coordinate==ZSTAR) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%dzt_dst(i,j,k) = 1.0 + Ext_mode%eta_t(i,j,taup1)*Grid%htr(i,j)
Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dzt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k)
Thickness%rho_dzt(i,j,k,taup1) = rho0*Thickness%dzt(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,taup1)+epsln)
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
elseif(vert_coordinate==ZSIGMA) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%dzt_dst(i,j,k) = max(depth_min_for_sigma,Grid%ht(i,j)+Ext_mode%eta_t(i,j,taup1))
Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dzt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k)
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
Thickness%rho_dzt(i,j,k,taup1) = rho0*Thickness%dzt(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,taup1)+epsln)
enddo
enddo
enddo
elseif(vert_coordinate==PRESSURE) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
! leave values over land untouched
if(Grid%tmask(i,j,k) > 0.0) then
Thickness%dzt_dst(i,j,k) = -grav_r/(Dens%rho(i,j,k,tau)+epsln)
Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dzt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k)
Thickness%rho_dzt(i,j,k,taup1) = Dens%rho(i,j,k,tau)*Thickness%dzt(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,taup1)+epsln)
endif
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
elseif(vert_coordinate==PSTAR) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
! leave values over land untouched
if(Grid%tmask(i,j,k) > 0.0) then
Thickness%dzt_dst(i,j,k) = -grav_r/(Dens%rho(i,j,k,tau)+epsln) &
*(Ext_mode%pbot_t(i,j,taup1)-Ext_mode%patm_t(i,j,taup1)) &
*Thickness%pbot0r(i,j)
Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dzt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k)
Thickness%rho_dzt(i,j,k,taup1) = Dens%rho(i,j,k,tau)*Thickness%dzt(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,taup1)+epsln)
endif
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
elseif(vert_coordinate==PSIGMA) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
if(Grid%tmask(i,j,k) > 0.0) then
Thickness%dzt_dst(i,j,k) = -grav_r/(Dens%rho(i,j,k,tau)+epsln) &
*(Ext_mode%pbot_t(i,j,taup1)-Ext_mode%patm_t(i,j,taup1))
density_tmp = Dens%rho(i,j,k,tau)
else
Thickness%dzt_dst(i,j,k) = -depth_min_for_sigma
density_tmp = rho0_profile(k)
endif
Thickness%dztlo(i,j,k) = Thickness%dstlo(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dztup(i,j,k) = Thickness%dstup(i,j,k)*Thickness%dzt_dst(i,j,k)
Thickness%dzt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k)
Thickness%rho_dzt(i,j,k,taup1) = density_tmp*Thickness%dzt(i,j,k)
Thickness%rho_dztr(i,j,k) = 1.0/(Thickness%rho_dzt(i,j,k,taup1)+epsln)
wrk1(i,j,k) = 1.0/Thickness%dzt_dst(i,j,k)
enddo
enddo
enddo
endif
! vertical distance (m) between t-cell points
if(Thickness%method==ENERGETIC) then
do k=1,nk-1
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,k) = 2.0*Thickness%dswt(i,j,k)/(wrk1(i,j,k)+wrk1(i,j,k+1))
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,0) = Thickness%dswt(i,j,0)*Thickness%dzt_dst(i,j,1)
kb=Grid%kmt(i,j)
if(kb > 0) then
Thickness%dzwt(i,j,kb) = Thickness%dswt(i,j,kb)*Thickness%dzt_dst(i,j,kb)
endif
enddo
enddo
elseif(Thickness%method==FINITEVOLUME) then
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,0) = Thickness%dztup(i,j,1)
enddo
enddo
do k=1,nk-1
do j=jsd,jed
do i=isd,ied
Thickness%dzwt(i,j,k) = Thickness%dztlo(i,j,k) + Thickness%dztup(i,j,k+1)
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
kb=Grid%kmt(i,j)
if(kb > 0) then
Thickness%dzwt(i,j,kb) = Thickness%dztlo(i,j,kb)
endif
enddo
enddo
endif
! vertical depth (m) of t and w points as
! computed via vertical sum of z-increments
if(vert_coordinate /= GEOPOTENTIAL) then
! compute depth from z=eta to t-point and to bottom of t-cell
do j=jsd,jed
do i=isd,ied
Thickness%depth_zt(i,j,1) = Thickness%dzwt(i,j,0)
Thickness%depth_zwt(i,j,1) = Thickness%dzt(i,j,1)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_zt(i,j,k) = Thickness%depth_zt(i,j,k-1) + Thickness%dzwt(i,j,k-1)
Thickness%depth_zwt(i,j,k) = Thickness%depth_zwt(i,j,k-1) + Thickness%dzt(i,j,k)
enddo
enddo
enddo
endif
! check for negative specific thickness, which occurs when eta_t is too negative,
! in which case the free surface is deeper than the ocean bottom topography.
do j=jsc,jec
do i=isc,iec
if(Ext_mode%eta_t(i,j,taup1) < -Grid%ht(i,j)) then
write(stdoutunit,*) &
'==>Error from ocean_thickness: Surface undulations too negative; model unstable'
write(stdoutunit,*) &
'eta_t(',i+Dom%ioff,',',j+Dom%joff,') = ',Ext_mode%eta_t(i,j,taup1)
write(stdoutunit,*) &
'depth(',i+Dom%ioff,',',j+Dom%joff,') = ',Grid%ht(i,j)
huge_undulations=.true.
endif
enddo
enddo
if(huge_undulations) then
call mpp_error(FATAL, &
'==>Error from ocean_thickness_mod: Free surface penetrating rock! Model unstable.')
endif
! check for negative thickness at top and bottom
error_flag=.false.
num_bad=0
do j=jsc,jec
do i=isc,iec
kb = Grid%kmt(i,j)
if(kb>0 .and. num_bad > max_num_bad_print) then
if(Thickness%dzt(i,j,1) < 0.0) then
num_bad = num_bad+1
write(stdoutunit,'(/a,i3,a,i3,a,i3,a,e10.4)') &
'==>Error: negative thickness at top: dzt(',&
i+Dom%ioff,',',j+Dom%joff,',',k,') =', Thickness%dzt(i,j,1)
if(enforce_positive_dzt) then
write(stdoutunit,'(a)') 'Resetting Thickness%dzt to ', thickness_dzt_min
Thickness%dzt(i,j,1) = thickness_dzt_min
else
write(errorstring,'(a,i3,a,i3,a,e10.4)') &
'(',i+Dom%ioff,',',j+Dom%joff,',1) =', Thickness%dzt(i,j,1)
call thickness_details(Grid, Time, Thickness, Ext_mode, Dens, i, j, kb, &
'From update_tcell_thickness with dzt<0', unit)
error_flag=.true.
endif
elseif(Thickness%dzt(i,j,kb) < 0.0) then
num_bad = num_bad+1
write(stdoutunit,'(/a,i3,a,i3,a,i3,a,e10.4)') &
'==>Error: negative thickness at bottom: dzt(',&
i+Dom%ioff,',',j+Dom%joff,',',kb,') =', Thickness%dzt(i,j,kb)
if(enforce_positive_dzt) then
write(stdoutunit,'(a)') 'Resetting Thickness%dzt to ', thickness_dzt_min
Thickness%dzt(i,j,kb) = thickness_dzt_min
else
write(errorstring,'(a,i3,a,i3,a,i3,a,e10.4)') &
'(',i+Dom%ioff,',',j+Dom%joff,',',kb,') =', Thickness%dzt(i,j,kb)
call thickness_details(Grid, Time, Thickness, Ext_mode, Dens, i, j, kb, &
'From update_tcell_thickness with dzt<0', unit)
error_flag=.true.
endif
endif
endif
enddo
enddo
if(error_flag) then
call mpp_error(FATAL, &
'==>Error from ocean_thickness_mod: dzt'//trim(errorstring)//'<0. grep for "dzt" to find location')
endif
if(debug_this_module_detail) then
call thickness_details(Grid, Time, Thickness, Ext_mode, Dens, isc, jsc, nk, &
'From update_tcell_thickness', stdoutunit)
call dzt_min_max(Time, Thickness, &
'From update_tcell_thickness with dzt<0, dzt_min_max information')
endif
end subroutine update_tcell_thickness
! NAME="update_tcell_thickness"
!#######################################################################
!
!
!
! Update time dependent thickness of U cells.
!
! Notes
!
! The computation of the depth arrays is a bit ad hoc. Here are
! the methods and their rationale.
!
! 1. For terrain following coordinates, the remap operator will entrain
! the very shallow layer T-points over land into the four-point averaging.
! This will cause the resulting U-cell depths next to land to be far shallower
! than what they should be. To avoid this problem, we set the U-cell depths the
! same as the T-cell depths. This specification causes no problems for budgets
! or energetic balance since Grid%umask array is already defined according to
! the usual B-grid specification.
!
! 2. For non-terrain and non-geopotential coordinates, compute U-cell thicknesses
! as min of surrounding T-cell thicknesses. This method follows that
! used with earlier MOM versions for the bottom partial step topography.
! Experiments with ZSTAR reveal that we must use the min function for dzu
! computation, rather than the alternative of a remap operator. If using
! the remap operator and nontrivial topography, then the velocity field can
! develop nontrivial noise. We have found that can compute dzwu using the
! remap operator without introducing noise. But we choose to use the
! min function to maintain compatibility with traditional approach in
! geopotential vertical coordinate versions of MOM.
!
!
subroutine update_ucell_thickness (Time, Grid, Ext_mode, Thickness)
type(ocean_time_type), intent(in) :: Time
type(ocean_grid_type), intent(in) :: Grid
type(ocean_external_mode_type), intent(in) :: Ext_mode
type(ocean_thickness_type), intent(inout) :: Thickness
integer :: i, j, k
integer :: tau, taup1
integer :: stdoutunit
stdoutunit=stdout()
tau = Time%tau
taup1 = Time%taup1
! write tau values of grid increments before update to taup1
if (id_dzu > 0) used = send_data (id_dzu, Thickness%dzu(:,:,:), &
Time%model_time, rmask=Grid%umask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_dzwu > 0) used = send_data (id_dzwu, Thickness%dzwu(:,:,1:nk), &
Time%model_time, rmask=Grid%umask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_rho_dzu > 0) used = send_data (id_rho_dzu, Thickness%rho_dzu(:,:,:,tau), &
Time%model_time, rmask=Grid%umask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_depth_zu > 0) used = send_data (id_depth_zu, Thickness%depth_zu(:,:,:), &
Time%model_time, rmask=Grid%umask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_depth_zwu > 0) used = send_data (id_depth_zwu, Thickness%depth_zwu(:,:,:), &
Time%model_time, rmask=Grid%umask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
if (id_mass_u > 0) used = send_data (id_mass_u, Thickness%mass_u(:,:,tau), &
Time%model_time, rmask=Grid%umask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
if (id_thicku > 0) used = send_data (id_thicku, Thickness%thicku(:,:,tau), &
Time%model_time, rmask=Grid%umask(:,:,1), &
is_in=isc, js_in=jsc, ie_in=iec, je_in=jec)
if(vert_coordinate==GEOPOTENTIAL) then
do j=jsd,jed
do i=isd,ied
Thickness%dzwu(i,j,0) = Thickness%depth_zt(i,j,1) + Ext_mode%eta_u(i,j,taup1)
Thickness%dzu(i,j,1) = Thickness%depth_zwt(i,j,1) + Ext_mode%eta_u(i,j,taup1)
Thickness%rho_dzu(i,j,1,taup1) = rho0*Thickness%dzu(i,j,1)
Thickness%rho_dzur(i,j,1) = 1.0/(Thickness%rho_dzu(i,j,1,taup1)+epsln)
enddo
enddo
if(linear_free_surface) then
do j=jsd,jed
do i=isd,ied
Thickness%dzwu(i,j,0) = Thickness%depth_zt(i,j,1)
Thickness%dzu(i,j,1) = Thickness%depth_zwt(i,j,1)
Thickness%rho_dzu(i,j,1,taup1) = rho0*Thickness%dzu(i,j,1)
Thickness%rho_dzur(i,j,1) = 1.0/(Thickness%rho_dzu(i,j,1,taup1)+epsln)
enddo
enddo
endif
else
if(vert_coordinate==ZSIGMA .or. vert_coordinate==PSIGMA) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%dzu(i,j,k) = Thickness%dzt(i,j,k)
Thickness%dzwu(i,j,k) = Thickness%dzwt(i,j,k)
Thickness%rho_dzu(i,j,k,taup1) = Thickness%rho_dzt(i,j,k,taup1)
enddo
enddo
enddo
k=0
do j=jsd,jed
do i=isd,ied
Thickness%dzwu(i,j,k) = Thickness%dzwt(i,j,k)
enddo
enddo
else
! U-cell thicknesses as minimum of surrounding T-cell thicknesses
do k=1,nk
do j=jsd,jec
do i=isd,iec
Thickness%dzu(i,j,k) = min(Thickness%dzt(i,j,k), &
Thickness%dzt(i+1,j,k), &
Thickness%dzt(i,j+1,k), &
Thickness%dzt(i+1,j+1,k))
Thickness%dzwu(i,j,k) = min(Thickness%dzwt(i,j,k), &
Thickness%dzwt(i+1,j,k), &
Thickness%dzwt(i,j+1,k), &
Thickness%dzwt(i+1,j+1,k))
Thickness%rho_dzu(i,j,k,taup1) = min(Thickness%rho_dzt(i,j,k,taup1), &
Thickness%rho_dzt(i+1,j,k,taup1), &
Thickness%rho_dzt(i,j+1,k,taup1), &
Thickness%rho_dzt(i+1,j+1,k,taup1))
enddo
enddo
enddo
if(update_dzwu_k0) then
k=0
do j=jsd,jec
do i=isd,iec
Thickness%dzwu(i,j,k) = min(Thickness%dzwt(i,j,k), &
Thickness%dzwt(i+1,j,k), &
Thickness%dzwt(i,j+1,k), &
Thickness%dzwt(i+1,j+1,k))
enddo
enddo
endif
call mpp_update_domains(Thickness%dzu(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%dzwu(:,:,:), Dom%domain2d)
call mpp_update_domains(Thickness%rho_dzu(:,:,:,taup1), Dom%domain2d)
endif
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzur(i,j,k) = Grid%umask(i,j,k)/(Thickness%rho_dzu(i,j,k,taup1)+epsln)
enddo
enddo
enddo
do j=jsd,jed
do i=isd,ied
Thickness%depth_zu(i,j,1) = Thickness%dzwu(i,j,0)
Thickness%depth_zwu(i,j,1) = Thickness%dzu(i,j,1)
enddo
enddo
do k=2,nk
do j=jsd,jed
do i=isd,ied
Thickness%depth_zu(i,j,k) = Thickness%depth_zu(i,j,k-1) + Thickness%dzwu(i,j,k-1)
Thickness%depth_zwu(i,j,k) = Thickness%depth_zwu(i,j,k-1) + Thickness%dzu(i,j,k)
enddo
enddo
enddo
endif
! thickness of U-cell column
do j=jsd,jed
do i=isd,ied
k=Grid%kmu(i,j)
if(k > 1) then
Thickness%thicku(i,j,taup1) = Thickness%depth_zwu(i,j,k)
else
Thickness%thicku(i,j,taup1) = 0.0
endif
enddo
enddo
! mass per area of U-cell column
Thickness%mass_u(:,:,taup1) = 0.0
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%mass_u(i,j,taup1) = Thickness%mass_u(i,j,taup1) &
+Grid%umask(i,j,k)*Thickness%rho_dzu(i,j,k,taup1)
enddo
enddo
enddo
if(debug_this_module) then
write(stdoutunit,*) ' '
write(stdoutunit,*) ' From ocean_thickness_mod: thickness checksum'
call write_timestamp(Time%model_time)
call thickness_chksum(Time, Grid, Thickness)
endif
end subroutine update_ucell_thickness
! NAME="update_ucell_thickness"
!#######################################################################
!
!
!
! Compute tendency for rho_dzt. This tendency is a function of the
! vertical coordinate.
!
!
subroutine rho_dzt_tendency(Time, Grid, Ext_mode, Thickness)
type(ocean_time_type), intent(in) :: Time
type(ocean_grid_type), intent(in) :: Grid
type(ocean_external_mode_type), intent(in) :: Ext_mode
type(ocean_thickness_type), intent(inout) :: Thickness
integer :: i, j, k
integer :: tau, taup1
tau = Time%tau
taup1 = Time%taup1
if(vert_coordinate==GEOPOTENTIAL) then
k=1
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzt_tendency(i,j,k) = Grid%tmask(i,j,k)*rho0*Ext_mode%deta_dt(i,j)
enddo
enddo
elseif(vert_coordinate==ZSTAR) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzt_tendency(i,j,k) = Grid%tmask(i,j,k)*rho0*Ext_mode%deta_dt(i,j) &
*Grid%htr(i,j)*Thickness%dst(i,j,k)
enddo
enddo
enddo
elseif(vert_coordinate==ZSIGMA) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzt_tendency(i,j,k) = Grid%tmask(i,j,k)*rho0*Ext_mode%deta_dt(i,j) &
*Thickness%dst(i,j,k)
enddo
enddo
enddo
elseif(vert_coordinate==PRESSURE) then
k=1
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzt_tendency(i,j,k) = -grav_r*Grid%tmask(i,j,k)*Ext_mode%dpatm_dt(i,j)
enddo
enddo
do j=jsd,jed
do i=isd,ied
k=Grid%kmt(i,j)
if(k>0) then
Thickness%rho_dzt_tendency(i,j,k) = grav_r*Grid%tmask(i,j,k)*Ext_mode%dpbot_dt(i,j)
endif
enddo
enddo
elseif(vert_coordinate==PSTAR) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzt_tendency(i,j,k) = -grav_r*Grid%tmask(i,j,k)*Thickness%dst(i,j,k) &
*(Ext_mode%dpbot_dt(i,j)-Ext_mode%dpatm_dt(i,j))*Thickness%pbot0r(i,j)
enddo
enddo
enddo
elseif(vert_coordinate==PSIGMA) then
do k=1,nk
do j=jsd,jed
do i=isd,ied
Thickness%rho_dzt_tendency(i,j,k) = -grav_r*Grid%tmask(i,j,k)*Thickness%dst(i,j,k) &
*(Ext_mode%dpbot_dt(i,j)-Ext_mode%dpatm_dt(i,j))
enddo
enddo
enddo
endif
if (id_rho_dzt_tendency > 0) used = send_data (id_rho_dzt_tendency, Thickness%rho_dzt_tendency(:,:,:), &
Time%model_time, rmask=Grid%tmask(:,:,:), &
is_in=isc, js_in=jsc, ks_in=1, ie_in=iec, je_in=jec, ke_in=nk)
end subroutine rho_dzt_tendency
! NAME="rho_dzt_tendency"
!#######################################################################
!
!
!
!
! Compute checksum for thickness components .
! Only print checksums for fields that should agree across restarts.
!
!
subroutine thickness_chksum(Time, Grid, Thickness)
type(ocean_time_type), intent(in) :: Time
type(ocean_grid_type), intent(in) :: Grid
type(ocean_thickness_type), intent(in) :: Thickness
integer :: taum1, tau, taup1
integer :: stdoutunit
stdoutunit=stdout()
taum1 = Time%taum1
tau = Time%tau
taup1 = Time%taup1
wrk1 = 0.0
if(tendency == THREE_LEVEL) then
wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzt(isc:iec,jsc:jec,:,tau)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%rho_dzt(tau) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzt(isc:iec,jsc:jec,:,taup1)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%rho_dzt(taup1) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzu(isc:iec,jsc:jec,:,tau)*Grid%umask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%rho_dzu(tau) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzu(isc:iec,jsc:jec,:,taup1)*Grid%umask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%rho_dzu(taup1) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,1) = Thickness%mass_u(isc:iec,jsc:jec,tau)*Grid%umask(isc:iec,jsc:jec,1)
write(stdoutunit,*) 'chksum for Thickness%mass_u(tau) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,1))
wrk1(isc:iec,jsc:jec,1) = Thickness%mass_u(isc:iec,jsc:jec,taup1)*Grid%umask(isc:iec,jsc:jec,1)
write(stdoutunit,*) 'chksum for Thickness%mass_u(taup1) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,1))
else
wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzt(isc:iec,jsc:jec,:,taup1)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%rho_dzt(taup1) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzu(isc:iec,jsc:jec,:,taup1)*Grid%umask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%rho_dzu(taup1) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,1) = Thickness%mass_u(isc:iec,jsc:jec,taup1)*Grid%umask(isc:iec,jsc:jec,1)
write(stdoutunit,*) 'chksum for Thickness%mass_u(taup1) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,1))
endif
wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dztr(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%rho_dztr = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzur(isc:iec,jsc:jec,:)*Grid%umask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%rho_dzur = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%rho_dzt_tendency(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%rho_dzt_tendency = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%dzt(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%dzt = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%dztlo(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%dztlo = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%dztup(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%dztup = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%dzt_dst(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%dzt_dst = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,1) = Thickness%dzwt(isc:iec,jsc:jec,0)*Grid%tmask(isc:iec,jsc:jec,1)
write(stdoutunit,*) 'chksum for Thickness%dzwt(k=0) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,1))
wrk1(isc:iec,jsc:jec,:) = Thickness%dzwt(isc:iec,jsc:jec,1:nk)*Grid%tmask(isc:iec,jsc:jec,1:nk)
write(stdoutunit,*) 'chksum for Thickness%dzwt(k=1:nk) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%dzu(isc:iec,jsc:jec,:)*Grid%umask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%dzu = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,1) = Thickness%dzwu(isc:iec,jsc:jec,0)*Grid%umask(isc:iec,jsc:jec,1)
write(stdoutunit,*) 'chksum for Thickness%dzwu(k=0) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,1))
wrk1(isc:iec,jsc:jec,:) = Thickness%dzwu(isc:iec,jsc:jec,1:nk)*Grid%umask(isc:iec,jsc:jec,1:nk)
write(stdoutunit,*) 'chksum for Thickness%dzwu(k=1:nk) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%depth_zt(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%depth_zt = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%geodepth_zt(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%geodepth_zt = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%depth_zu(isc:iec,jsc:jec,:)*Grid%umask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%depth_zu = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%depth_zwt(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%depth_zwt = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%depth_zwu(isc:iec,jsc:jec,:)*Grid%umask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%depth_zwu = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%depth_st(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%depth_st = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%depth_swt(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%depth_swt = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%dst(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%dst = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%dstlo(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%dstlo = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,:) = Thickness%dstup(isc:iec,jsc:jec,:)*Grid%tmask(isc:iec,jsc:jec,:)
write(stdoutunit,*) 'chksum for Thickness%dstup = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,1) = Thickness%dswt(isc:iec,jsc:jec,0)*Grid%tmask(isc:iec,jsc:jec,1)
write(stdoutunit,*) 'chksum for Thickness%dswt(k=0) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,1))
wrk1(isc:iec,jsc:jec,:) = Thickness%dswt(isc:iec,jsc:jec,1:nk)*Grid%tmask(isc:iec,jsc:jec,1:nk)
write(stdoutunit,*) 'chksum for Thickness%dswt(k=1:nk) = ', mpp_chksum(wrk1(isc:iec,jsc:jec,:))
wrk1(isc:iec,jsc:jec,1) = Thickness%pbot0(isc:iec,jsc:jec)*Grid%tmask(isc:iec,jsc:jec,1)
write(stdoutunit,*) 'chksum for Thickness%pbot0 = ', mpp_chksum(wrk1(isc:iec,jsc:jec,1))
write(stdoutunit,*) ' '
return
end subroutine thickness_chksum
! NAME="thickness_chksum"
!#######################################################################
!
!
!
!
! For debugging, we print here some details of the grid at a particular
! (i,j) point.
!
!
subroutine thickness_details(Grid, Time, Thickness, Ext_mode, Dens, &
ipoint, jpoint, kb, filecaller, outunit)
type(ocean_grid_type), intent(in) :: Grid
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
type(ocean_external_mode_type), intent(in) :: Ext_mode
type(ocean_density_type), intent(in) :: Dens
integer, intent(in) :: ipoint
integer, intent(in) :: jpoint
integer, intent(in) :: kb
character(len=*), intent(in) :: filecaller
integer, intent(in) :: outunit
integer :: i,j,k
integer :: taup1, tau
taup1 = Time%taup1
tau = Time%tau
write(outunit,*) ' '
write(outunit,'(a)') trim(filecaller)
call write_timestamp(Time%model_time)
i=ipoint
j=jpoint
write(outunit,*) ' '
write(outunit,'(a,i4,a,i4,a,f24.8)') &
'ht(',i+Dom%ioff,',',j+Dom%joff,')(metres) = ',&
Grid%ht(i,j)
write(outunit,'(a,i4,a,i4,a,f24.8)') &
'eta_t(',i+Dom%ioff,',',j+Dom%joff,')(metres) = ',&
Ext_mode%eta_t(i,j,taup1)
write(outunit,'(a,i4,a,i4,a,f24.8)') &
'pbot_t(',i+Dom%ioff,',',j+Dom%joff,')(dbars) = ',&
Ext_mode%pbot_t(i,j,taup1)
write(outunit,'(a,i4,a,i4,a,f24.8)') &
'patm_t(',i+Dom%ioff,',',j+Dom%joff,')(dbars) = ',&
Ext_mode%patm_t(i,j,taup1)
write(outunit,'(a,i4,a,i4,a,f24.8)') &
'pbot0(',i+Dom%ioff,',',j+Dom%joff,')(dbars) = ',&
Thickness%pbot0(i,j)*c2dbars
write(outunit,'(a,i4,a,i4,a,f24.8)') &
'pbot0r(',i+Dom%ioff,',',j+Dom%joff,')(1/dbar)= ',&
Thickness%pbot0r(i,j)/c2dbars
write(outunit,'(a,i4,a,i4,a,f24.8)') &
'xt(',i+Dom%ioff,',',j+Dom%joff,')(longitude) = ',&
Grid%xt(i,j)
write(outunit,'(a,i4,a,i4,a,f24.8)') &
'yt(',i+Dom%ioff,',',j+Dom%joff,')(latitude) = ', &
Grid%yt(i,j)
write(outunit,*) ' '
do k=1,kb
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'depth_st(',i+Dom%ioff,',',j+Dom%joff,',',k,')(s-units) = ',&
Thickness%depth_st(i,j,k)*convert_factor
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'depth_swt(',i+Dom%ioff,',',j+Dom%joff,',',k,')(s-units) = ',&
Thickness%depth_swt(i,j,k)*convert_factor
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'rho_t(',i+Dom%ioff,',',j+Dom%joff,',',k,')(kg/m^3) = ',&
Dens%rho(i,j,k,tau)
enddo
write(outunit,*) ' '
do k=1,kb
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dst (',i+Dom%ioff,',',j+Dom%joff,',',k,')(s-units)= ',&
Thickness%dst(i,j,k)*convert_factor
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dswt (',i+Dom%ioff,',',j+Dom%joff,',',k,')(s-units)= ',&
Thickness%dswt(i,j,k)*convert_factor
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dstlo(',i+Dom%ioff,',',j+Dom%joff,',',k,')(s-units)= ',&
Thickness%dstlo(i,j,k)*convert_factor
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dstup(',i+Dom%ioff,',',j+Dom%joff,',',k,')(s-units)= ',&
Thickness%dstup(i,j,k)*convert_factor
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dstlo + dstup(',i+Dom%ioff,',',j+Dom%joff,',',k,') = ',&
(Thickness%dstlo(i,j,k)+Thickness%dstup(i,j,k))*convert_factor
enddo
write(outunit,*) ' '
do k=1,kb
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'depth_zt(',i+Dom%ioff,',',j+Dom%joff,',',k,')(metres) = ',&
Thickness%depth_zt(i,j,k)
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'depth_zu(',i+Dom%ioff,',',j+Dom%joff,',',k,')(metres) = ',&
Thickness%depth_zu(i,j,k)
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'depth_zwt(',i+Dom%ioff,',',j+Dom%joff,',',k,')(metres)= ',&
Thickness%depth_zwt(i,j,k)
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'depth_zwu(',i+Dom%ioff,',',j+Dom%joff,',',k,')(metres)= ',&
Thickness%depth_zwu(i,j,k)
enddo
write(outunit,*) ' '
do k=1,kb
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dzt(',i+Dom%ioff,',',j+Dom%joff,',',k,')(metres) = ',Thickness%dzt(i,j,k)
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dzu(',i+Dom%ioff,',',j+Dom%joff,',',k,')(metres) = ',Thickness%dzu(i,j,k)
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dzwt(',i+Dom%ioff,',',j+Dom%joff,',',k,')(metres) = ',Thickness%dzwt(i,j,k)
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dzwu(',i+Dom%ioff,',',j+Dom%joff,',',k,')(metres) = ',Thickness%dzwu(i,j,k)
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dztlo(',i+Dom%ioff,',',j+Dom%joff,',',k,')(metres)= ',Thickness%dztlo(i,j,k)
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dztup(',i+Dom%ioff,',',j+Dom%joff,',',k,')(metres)= ',Thickness%dztup(i,j,k)
write(outunit,'(a,i4,a,i4,a,i4,a,f24.8)') &
'dztlo+dztup(',i+Dom%ioff,',',j+Dom%joff,',',k,') = ',&
Thickness%dztlo(i,j,k)+Thickness%dztup(i,j,k)
enddo
end subroutine thickness_details
! NAME="thickness_details"
!#######################################################################
!
!
! Write out restart files registered through register_restart_file
!
subroutine ocean_thickness_restart(Time, Thickness, time_stamp)
type(ocean_time_type), intent(in) :: Time
type(ocean_thickness_type), intent(in) :: Thickness
character(len=*), intent(in), optional :: time_stamp
integer :: tau, taup1
tau = Time%tau
taup1 = Time%taup1
if(tendency==THREE_LEVEL) then
call reset_field_pointer(Thk_restart, id_restart_rho_dzt, Thickness%rho_dzt(:,:,:,tau), &
Thickness%rho_dzt(:,:,:,taup1) )
elseif(tendency==TWO_LEVEL) then
call reset_field_pointer(Thk_restart, id_restart_rho_dzt, Thickness%rho_dzt(:,:,:,taup1) )
end if
call save_restart(Thk_restart, time_stamp)
end subroutine ocean_thickness_restart
! NAME="ocean_thickness_restart"
!#######################################################################
!
!
!
! Write basic elements of thickness derived type to restart
!
!
subroutine ocean_thickness_end (Time, Grid, Thickness)
type(ocean_time_type), intent(in) :: Time
type(ocean_grid_type), intent(in) :: Grid
type(ocean_thickness_type), intent(inout) :: Thickness
integer :: stdoutunit
stdoutunit=stdout()
if (.not.module_is_initialized) then
call mpp_error(FATAL, &
'==>Error from ocean_thickness_mod (ocean_thickness_end): module must be initialized')
endif
if(.not. write_a_restart) then
write(stdoutunit,'(/a)') &
'==>Warning from ocean_thickness_mod (ocean_thickness_end): NO restart written.'
call mpp_error(WARNING, &
'==>Warning from ocean_thickness_mod (ocean_thickness_end): NO restart written.')
return
endif
call ocean_thickness_restart(Time, Thickness)
write(stdoutunit,*) ' '
write(stdoutunit,*) 'From ocean_thickness_mod: ending thickness checksums'
call write_timestamp(Time%model_time)
call thickness_chksum(Time, Grid, Thickness)
nullify(Dom)
end subroutine ocean_thickness_end
! NAME="ocean_thickness_end"
!#######################################################################
!
!
!
! REMAP_ZT_TO_ZU remaps a T-cell thickness or vertical depth on a
! T-cell to the corresponding place on U-cell.
!
! This is the same operator as REMAP_BT_TO_BU.
! It is needed for ocean_thickness_init since this routine is called
! prior to ocean_operators_init. This is a bit awkward, but
! ocean_operators needs thickness values, so it must be called
! after thickness is initialized.
!
!
function REMAP_ZT_TO_ZU(a,Grid)
real, dimension(isd:,jsd:), intent(in) :: a
type(ocean_grid_type), intent(in) :: Grid
real, dimension(isd:ied,jsd:jed) :: REMAP_ZT_TO_ZU
integer :: i, j
do j=jsc-halo,jec+halo-1
do i=isc-halo,iec+halo-1
REMAP_ZT_TO_ZU(i,j) = (a(i,j) *Grid%dte(i,j) *Grid%dus(i,j) &
+ a(i+1,j) *Grid%dtw(i+1,j) *Grid%dus(i,j) &
+ a(i,j+1) *Grid%dte(i,j+1) *Grid%dun(i,j) &
+ a(i+1,j+1)*Grid%dtw(i+1,j+1)*Grid%dun(i,j) &
)*Grid%daur(i,j)
enddo
enddo
REMAP_ZT_TO_ZU(iec+halo,:) = 0.0
REMAP_ZT_TO_ZU(:,jec+halo) = 0.0
end function REMAP_ZT_TO_ZU
! NAME="REMAP_ZT_TO_ZU"
end module ocean_thickness_mod