ns5001k-admin/Nemo-NAA10km_nemo_4.2.2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Nemo-NAA10km_nemo_4.2.2/tools/DOMAINcfg/src/domzgr.F90
Robinson Hordoir 4fa5daa6c3 First commit of Nemo-NAA10km_nemo_4.2.2
2025-09-05 09:42:21 +02:00

2216 lines
111 KiB
Fortran
Executable File
Raw Blame History

MODULE domzgr
!!==============================================================================
!! *** MODULE domzgr ***
!! Ocean domain : definition of the vertical coordinate system
!!==============================================================================
!! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate
!! ! 1997-07 (G. Madec) lbc_lnk call
!! ! 1997-04 (J.-O. Beismann)
!! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module
!! - ! 2002-09 (A. de Miranda) rigid-lid + islands
!! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module
!! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function
!! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco
!! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
!! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level
!! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function
!! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case
!! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye
!! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
!! dom_zgr : defined the ocean vertical coordinate system
!! zgr_bat : bathymetry fields (levels and meters)
!! zgr_bat_zoom : modify the bathymetry field if zoom domain
!! zgr_bat_ctl : check the bathymetry files
!! zgr_bot_level: deepest ocean level for t-, u, and v-points
!! zgr_z : reference z-coordinate
!! zgr_zco : z-coordinate
!! zgr_zps : z-coordinate with partial steps
!! zgr_sco : s-coordinate
!! fssig : tanh stretch function
!! fssig1 : Song and Haidvogel 1994 stretch function
!! fgamma : Siddorn and Furner 2012 stretching function
!!---------------------------------------------------------------------
USE dom_oce ! ocean domain
USE depth_e3 ! depth <=> e3
!
USE in_out_manager ! I/O manager
USE iom ! I/O library
USE lbclnk ! ocean lateral boundary conditions (or mpp link)
USE lib_mpp ! distributed memory computing library
USE lib_fortran
USE dombat
USE domisf
#if defined key_agrif
USE agrif_connect
USE agrif_domzgr
#endif
IMPLICIT NONE
PRIVATE
PUBLIC dom_zgr ! called by dom_init.F90
!
! !!* Namelist namzgr_sco *
LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
LOGICAL :: ln_s_sf12 ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
!
REAL(wp), PUBLIC :: rn_sbot_min ! minimum depth of s-bottom surface (>0) (m)
REAL(wp) :: rn_sbot_max ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
REAL(wp) :: rn_rmax ! maximum cut-off r-value allowed (0<rn_rmax<1)
REAL(wp) :: rn_hc ! Critical depth for transition from sigma to stretched coordinates
! Song and Haidvogel 1994 stretching parameters
REAL(wp) :: rn_theta ! surface control parameter (0<=rn_theta<=20)
REAL(wp) :: rn_thetb ! bottom control parameter (0<=rn_thetb<= 1)
REAL(wp) :: rn_bb ! stretching parameter
! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
! Siddorn and Furner stretching parameters
LOGICAL :: ln_sigcrit ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
REAL(wp) :: rn_alpha ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
REAL(wp) :: rn_efold ! efold length scale for transition to stretched coord
REAL(wp) :: rn_zs ! depth of surface grid box
! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
REAL(wp) :: rn_zb_a ! bathymetry scaling factor for calculating Zb
REAL(wp) :: rn_zb_b ! offset for calculating Zb
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dommsk.F90 13305 2020-07-14 17:12:25Z acc $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dom_zgr( k_top, k_bot )
!!----------------------------------------------------------------------
!! *** ROUTINE dom_zgr ***
!!
!! ** Purpose : set the depth of model levels and the resulting
!! vertical scale factors.
!!
!! ** Method : - reference 1D vertical coordinate (gdep._1d, e3._1d)
!! - read/set ocean depth and ocean levels (bathy, mbathy)
!! - vertical coordinate (gdep., e3.) depending on the
!! coordinate chosen :
!! ln_zco=T z-coordinate
!! ln_zps=T z-coordinate with partial steps
!! ln_zco=T s-coordinate
!!
!! ** Action : define gdep., e3., mbathy and bathy
!!----------------------------------------------------------------------
INTEGER, DIMENSION(:,:), INTENT(out) :: k_top, k_bot ! ocean first and last level indices
!
INTEGER :: ioptio, ibat ! local integer
INTEGER :: ios
!
INTEGER :: ji,jj,jk
REAL(wp) :: zrefdep ! depth of the reference level (~10m)
NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
!!----------------------------------------------------------------------
!
!
! REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate
READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist')
! REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate
READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist')
IF(lwm) WRITE ( numond, namzgr )
#if defined key_agrif
IF (.NOT.Agrif_root()) THEN
IF (.NOT.(ln_dept_mid.AND.ln_e3_dep).OR. &
& .NOT.(Agrif_Parent(ln_dept_mid).AND.Agrif_Parent(ln_e3_dep)) ) THEN
CALL ctl_stop( 'STOP', 'namcfg: with AGRIF zooms, set ln_dept_mid = ln_e3_dep = T' )
ENDIF
ENDIF
#endif
IF(ln_read_cfg) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>>> Read vertical mesh in ', TRIM( cn_domcfg ), ' file'
!
CALL zgr_read ( ln_zco , ln_zps , ln_sco, ln_isfcav, &
& gdept_1d, gdepw_1d, e3t_1d, e3w_1d , & ! 1D gridpoints depth
& gdept_0 , gdepw_0 , & ! gridpoints depth
& e3t_0 , e3u_0 , e3v_0 , e3f_0 , & ! vertical scale factors
& e3w_0 , e3uw_0 , e3vw_0 , & ! vertical scale factors
& k_top , k_bot ) ! 1st & last ocean level
!
!!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears
! ! Compute gde3w_0 (vertical sum of e3w)
! gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
! DO jk = 2, jpk
! gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
! END DO
!
! ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top)
CALL zgr_top_bot( k_top, k_bot ) ! with a minimum value set to 1
! ! deepest/shallowest W level Above/Below ~10m
!!gm BUG in s-coordinate this does not work!
zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
nla10 = nlb10 - 1 ! deepest W level Above ~10m
!!gm end bug
ENDIF
IF(lwp) THEN ! Control print
WRITE(numout,*)
WRITE(numout,*) 'dom_zgr : vertical coordinate'
WRITE(numout,*) '~~~~~~~'
WRITE(numout,*) ' Namelist namzgr : set vertical coordinate'
WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav
ENDIF
ioptio = 0 ! Check Vertical coordinate options
IF( ln_zco ) ioptio = ioptio + 1
IF( ln_zps ) ioptio = ioptio + 1
IF( ln_sco ) ioptio = ioptio + 1
IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
!
IF ( ln_isfcav ) CALL zgr_isf_nam
ioptio = 0
IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1
IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1
IF( ioptio > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' )
#if defined key_agrif
! IF ( (.NOT.Agrif_Root()).AND.((.NOT.ln_zps).OR.(.NOT.Agrif_parent(ln_zps)))) THEN
! CALL ctl_stop( 'STOP', 'AGRIF zooms require ln_zps=T for both Child and Parent')
! ENDIF
#endif
IF(.NOT.ln_read_cfg) THEN
!
! Build the vertical coordinate system
! ------------------------------------
CALL zgr_z ! Reference z-coordinate system (always called)
CALL zgr_bat ! Bathymetry fields (levels and meters)
IF( ln_zco ) CALL zgr_zco ! z-coordinate
IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points
!
! final adjustment of mbathy & check
! -----------------------------------
CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
k_bot = mbkt
CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points
k_top = mikt
WHERE( bathy(:,:) <= 0._wp ) k_top(:,:) = 0 ! set k_top to zero over land
ENDIF
!
#if defined key_agrif
IF ( Agrif_root() ) THEN
#endif
IF ( ln_read_cfg.AND.ln_e3_dep.AND.ln_dept_mid ) THEN
! Define depths at cell centers
DO jk = 1, jpkm1
gdept_1d(jk) = 0.5_wp * (gdepw_1d(jk)+gdepw_1d(jk+1))
END DO
DO jk = 2, jpk
e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)
END DO
e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
gdept_0(ji,jj,jk) = 0.5_wp * ( gdepw_0(ji,jj,jk) + gdepw_0(ji,jj,jk+1) )
END_3D
e3w_0(:,:,1) = e3t_0(:,:,1)
e3uw_0(:,:,1) = e3u_0(:,:,1)
e3vw_0(:,:,1) = e3v_0(:,:,1)
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 2, jpk )
e3w_0(ji,jj,jk) = 0.5_wp * (e3t_0(ji,jj,jk)+e3t_0(ji,jj,jk-1))
e3uw_0(ji,jj,jk) = 0.5_wp * (e3u_0(ji,jj,jk)+e3u_0(ji,jj,jk-1))
e3vw_0(ji,jj,jk) = 0.5_wp * (e3v_0(ji,jj,jk)+e3v_0(ji,jj,jk-1))
END_3D
ENDIF
#if defined key_agrif
ENDIF
#endif
!
IF( lwp ) THEN
WRITE(numout,*) ' MIN val k_top ', MINVAL( k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) )
WRITE(numout,*) ' MIN val k_bot ', MINVAL( k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) )
WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
& ' w ', MINVAL( gdepw_0(:,:,:) )
WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), &
& ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), &
& ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), &
& ' w ', MINVAL( e3w_0(:,:,:) )
WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
& ' w ', MAXVAL( gdepw_0(:,:,:) )
WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), &
& ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), &
& ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), &
& ' w ', MAXVAL( e3w_0(:,:,:) )
ENDIF
!
END SUBROUTINE dom_zgr
SUBROUTINE zgr_read( ld_zco , ld_zps , ld_sco , ld_isfcav, & ! type of vertical coordinate
& pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d , & ! 1D reference vertical coordinate
& pdept , pdepw , & ! 3D t & w-points depth
& pe3t , pe3u , pe3v , pe3f , & ! vertical scale factors
& pe3w , pe3uw , pe3vw , & ! - - -
& k_top , k_bot ) ! top & bottom ocean level
!!---------------------------------------------------------------------
!! *** ROUTINE zgr_read ***
!!
!! ** Purpose : Read the vertical information in the domain configuration file
!!
!!----------------------------------------------------------------------
LOGICAL , INTENT(out) :: ld_zco, ld_zps, ld_sco ! vertical coordinate flags
LOGICAL , INTENT(out) :: ld_isfcav ! under iceshelf cavity flag
REAL(wp), DIMENSION(:) , INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m]
REAL(wp), DIMENSION(:) , INTENT(out) :: pe3t_1d , pe3w_1d ! 1D vertical scale factors [m]
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pdept, pdepw ! grid-point depth [m]
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3t , pe3u , pe3v , pe3f ! vertical scale factors [m]
REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3w , pe3uw, pe3vw ! - - -
INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top , k_bot ! first & last ocean level
!
INTEGER :: jk ! dummy loop index
INTEGER :: inum ! local logical unit
REAL(WP) :: z_zco, z_zps, z_sco, z_cav
REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
!!----------------------------------------------------------------------
!
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) ' zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file'
WRITE(numout,*) ' ~~~~~~~~'
ENDIF
!
CALL iom_open( cn_domcfg, inum )
!
! !* type of vertical coordinate
CALL iom_get( inum, 'ln_zco' , z_zco )
CALL iom_get( inum, 'ln_zps' , z_zps )
CALL iom_get( inum, 'ln_sco' , z_sco )
IF( z_zco == 0._wp ) THEN ; ld_zco = .false. ; ELSE ; ld_zco = .true. ; ENDIF
IF( z_zps == 0._wp ) THEN ; ld_zps = .false. ; ELSE ; ld_zps = .true. ; ENDIF
IF( z_sco == 0._wp ) THEN ; ld_sco = .false. ; ELSE ; ld_sco = .true. ; ENDIF
!
! !* ocean cavities under iceshelves
CALL iom_get( inum, 'ln_isfcav', z_cav )
IF( z_cav == 0._wp ) THEN ; ld_isfcav = .false. ; ELSE ; ld_isfcav = .true. ; ENDIF
!
! !* vertical scale factors
CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , pe3t_1d ) ! 1D reference coordinate
CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , pe3w_1d )
!
CALL iom_get( inum, jpdom_global, 'e3t_0' , pe3t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy ) ! 3D coordinate
CALL iom_get( inum, jpdom_global, 'e3u_0' , pe3u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e3v_0' , pe3v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e3f_0' , pe3f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e3w_0' , pe3w , cd_type = 'W', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e3uw_0' , pe3uw, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global, 'e3vw_0' , pe3vw, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
!
! !* depths
! !- old depth definition (obsolescent feature)
IF( iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0 .AND. &
& iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0 .AND. &
& iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0 .AND. &
& iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0 ) THEN
CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', &
& ' depths at t- and w-points read in the domain configuration file')
CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )
CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d )
CALL iom_get( inum, jpdom_global , 'gdept_0' , pdept, kfill = jpfillcopy )
CALL iom_get( inum, jpdom_global , 'gdepw_0' , pdepw, kfill = jpfillcopy )
!
ELSE !- depths computed from e3. scale factors
CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d ) ! 1D reference depth
CALL e3_to_depth( pe3t , pe3w , pdept , pdepw ) ! 3D depths
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) ' Reference 1D z-coordinate depth and scale factors:'
WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
ENDIF
ENDIF
!
! !* ocean top and bottom level
CALL iom_get( inum, jpdom_global, 'top_level' , z2d ) ! 1st wet T-points (ISF)
k_top(:,:) = NINT( z2d(:,:) )
CALL iom_get( inum, jpdom_global, 'bottom_level' , z2d ) ! last wet T-points
k_bot(:,:) = NINT( z2d(:,:) )
!
! reference depth for negative bathy (wetting and drying only)
! IF( ll_wd ) CALL iom_get( inum, 'rn_wd_ref_depth' , ssh_ref )
!
CALL iom_close( inum )
!
END SUBROUTINE zgr_read
SUBROUTINE zgr_top_bot( k_top, k_bot )
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_top_bot ***
!!
!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
!!
!! ** Method : computes from k_top and k_bot with a minimum value of 1 over land
!!
!! ** Action : mikt, miku, mikv : vertical indices of the shallowest
!! ocean level at t-, u- & v-points
!! (min value = 1)
!! ** Action : mbkt, mbku, mbkv, mbkf : vertical indices of the deeptest
!! ocean level at t-, u-, v- & f-points
!! (min value = 1 over land)
!!----------------------------------------------------------------------
INTEGER , DIMENSION(:,:), INTENT(in) :: k_top, k_bot ! top & bottom ocean level indices
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp), DIMENSION(jpi,jpj) :: zk ! workspace
!!----------------------------------------------------------------------
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' zgr_top_bot : ocean top and bottom k-index of T-, U-, V- and W-levels '
IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
!
mikt(:,:) = MAX( k_top(:,:) , 1 ) ! top ocean k-index of T-level (=1 over land)
!
mbkt(:,:) = MAX( k_bot(:,:) , 1 ) ! bottom ocean k-index of T-level (=1 over land)
! ! N.B. top k-index of W-level = mikt
! ! bottom k-index of W-level = mbkt+1
DO jj = 1, jpjm1
DO ji = 1, jpim1
miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) )
mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) )
mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) )
!
mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
mbkf(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj), mbkt(ji+1,jj ), mbkt(ji+1,jj+1) )
END DO
END DO
! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
zk(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'U', 1. ) ; miku(:,:) = MAX( NINT( zk(:,:) ), 1 )
zk(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'V', 1. ) ; mikv(:,:) = MAX( NINT( zk(:,:) ), 1 )
zk(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'F', 1. ) ; mikf(:,:) = MAX( NINT( zk(:,:) ), 1 )
!
zk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'U', 1. ) ; mbku(:,:) = MAX( NINT( zk(:,:) ), 1 )
zk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'V', 1. ) ; mbkv(:,:) = MAX( NINT( zk(:,:) ), 1 )
zk(:,:) = REAL( mbkf(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'F', 1. ) ; mbkf(:,:) = MAX( NINT( zk(:,:) ), 1 )
!
END SUBROUTINE zgr_top_bot
SUBROUTINE zgr_z
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_z ***
!!
!! ** Purpose : set the depth of model levels and the resulting
!! vertical scale factors.
!!
!! ** Method : z-coordinate system (use in all type of coordinate)
!! The depth of model levels is defined from an analytical
!! function the derivative of which gives the scale factors.
!! both depth and scale factors only depend on k (1d arrays).
!! w-level: gdepw_1d = gdep(k)
!! e3w_1d(k) = dk(gdep)(k) = e3(k)
!! t-level: gdept_1d = gdep(k+0.5)
!! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
!!
!! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
!! - e3t_1d , e3w_1d : scale factors at T- and W-levels (m)
!!
!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
!!----------------------------------------------------------------------
INTEGER :: jk ! dummy loop indices
REAL(wp) :: zt, zw ! temporary scalars
REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
REAL(wp) :: zrefdep ! depth of the reference level (~10m)
REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
!!----------------------------------------------------------------------
!
! Set variables from parameters
! ------------------------------
zkth = ppkth ; zacr = ppacr
zdzmin = ppdzmin ; zhmax = pphmax
zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters
! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
IF( ppa1 == pp_to_be_computed .AND. &
& ppa0 == pp_to_be_computed .AND. &
& ppsur == pp_to_be_computed ) THEN
!
za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
ELSE
za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter
ENDIF
IF(lwp) THEN ! Parameter print
WRITE(numout,*)
WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
WRITE(numout,*) ' ~~~~~~~'
IF( ppkth == 0._wp ) THEN
WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
WRITE(numout,*) ' Total depth :', zhmax
WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
ELSE
IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
WRITE(numout,*) ' zsur, za0, za1 computed from '
WRITE(numout,*) ' zdzmin = ', zdzmin
WRITE(numout,*) ' zhmax = ', zhmax
ENDIF
WRITE(numout,*) ' Value of coefficients for vertical mesh:'
WRITE(numout,*) ' zsur = ', zsur
WRITE(numout,*) ' za0 = ', za0
WRITE(numout,*) ' za1 = ', za1
WRITE(numout,*) ' zkth = ', zkth
WRITE(numout,*) ' zacr = ', zacr
IF( ldbletanh ) THEN
WRITE(numout,*) ' (Double tanh za2 = ', za2
WRITE(numout,*) ' parameters) zkth2= ', zkth2
WRITE(numout,*) ' zacr2= ', zacr2
ENDIF
ENDIF
ENDIF
! Reference z-coordinate (depth - scale factor at T- and W-points)
! ======================
IF( ppkth == 0._wp ) THEN ! uniform vertical grid
za1 = zhmax / FLOAT(jpkdta-1)
DO jk = 1, jpk
zw = FLOAT( jk )
zt = FLOAT( jk ) + 0.5_wp
gdepw_1d(jk) = ( zw - 1 ) * za1
gdept_1d(jk) = ( zt - 1 ) * za1
e3w_1d (jk) = za1
e3t_1d (jk) = za1
END DO
ELSE ! Madec & Imbard 1996 function
IF( .NOT. ldbletanh ) THEN
DO jk = 1, jpk
zw = REAL( jk , wp )
zt = REAL( jk , wp ) + 0.5_wp
gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
END DO
ELSE
DO jk = 1, jpk
zw = FLOAT( jk )
zt = FLOAT( jk ) + 0.5_wp
! Double tanh function
gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) &
& + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) )
gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) &
& + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) )
e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) &
& + za2 * TANH( (zw-zkth2) / zacr2 )
e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) &
& + za2 * TANH( (zt-zkth2) / zacr2 )
END DO
ENDIF
gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero
ENDIF
IF ( ln_isfcav .OR. ln_e3_dep ) THEN ! e3. = dk[gdep]
!
DO jk = 1, jpkm1
e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)
END DO
e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO
IF ( ln_dept_mid ) THEN
DO jk = 1, jpkm1
gdept_1d(jk) = 0.5_wp * (gdepw_1d(jk)+gdepw_1d(jk+1))
END DO
ENDIF
DO jk = 2, jpk
e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)
END DO
e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))
END IF
!!gm BUG in s-coordinate this does not work!
! deepest/shallowest W level Above/Below ~10m
zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
nla10 = nlb10 - 1 ! deepest W level Above ~10m
!!gm end bug
IF(lwp) THEN ! control print
WRITE(numout,*)
WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
ENDIF
DO jk = 1, jpk ! control positivity
IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 ' )
IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
END DO
!
END SUBROUTINE zgr_z
SUBROUTINE zgr_bat
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_bat ***
!!
!! ** Purpose : set bathymetry both in levels and meters
!!
!! ** Method : read or define mbathy and bathy arrays
!! * level bathymetry:
!! The ocean basin geometry is given by a two-dimensional array,
!! mbathy, which is defined as follow :
!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
!! at t-point (ji,jj).
!! = 0 over the continental t-point.
!! The array mbathy is checked to verified its consistency with
!! model option. in particular:
!! mbathy must have at least 1 land grid-points (mbathy<=0)
!! along closed boundary.
!! mbathy must be cyclic IF jperio=1.
!! mbathy must be lower or equal to jpk-1.
!! isolated ocean grid points are suppressed from mbathy
!! since they are only connected to remaining
!! ocean through vertical diffusion.
!! ntopo=-1 : rectangular channel or bassin with a bump
!! ntopo= 0 : flat rectangular channel or basin
!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
!! bathy is read in 'bathy_meter.nc' NetCDF file
!!
!! ** Action : - mbathy: level bathymetry (in level index)
!! - bathy : meter bathymetry (in meters)
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: inum ! temporary logical unit
INTEGER :: ierror ! error flag
INTEGER :: ii_bump, ij_bump, ih ! bump center position
INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices
REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
REAL(wp) :: zi, zj, zh, zhmin ! local scalars
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrand, zbatv
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztemp
!!----------------------------------------------------------------------
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
IF(lwp) WRITE(numout,*) ' ~~~~~~~'
! ! ================== !
IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
! ! ================== !
!
!
IF( ntopo == 0 ) THEN ! flat basin
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
IF( rn_bathy > 0.01 ) THEN
IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist'
bathy(:,:) = rn_bathy
IF( ln_sco ) THEN ! s-coordinate (zsc ):
mbathy(:,:) = jpkm1
ELSE ! z-coordinate (zco or zps): step-like topography
mbathy(:,:) = jpkm1
DO jk = 1, jpkm1
WHERE( gdept_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= gdept_1d(jk+1) ) mbathy(:,:) = jk
END DO
ENDIF
ELSE
IF(lwp) WRITE(numout,*) ' Depth = depthw(jpkm1)'
mbathy(:,:) = jpkm1 ! before last level
bathy(:,:) = gdepw_1d(jpk) ! last w-point depth
h_oce = gdepw_1d(jpk)
ENDIF
ELSE ! flat basin with random noise
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with random noise'
ALLOCATE( zrand(jpiglo,jpjglo), STAT=ierror )
IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zrand array' )
h_oce = gdepw_1d(jpk) ! background ocean depth (meters)
!
CALL RANDOM_SEED()
CALL RANDOM_NUMBER(zrand)
DO_2D( 0, 0, 0, 0 )
bathy(ji,jj) = h_oce + 0.1_wp *h_oce * (zrand(mig(ji),mjg(jj))-1._wp)
END_2D
IF ( cp_cfg=='OVERFLOW' ) THEN
DO jj=1,jpj
bathy(:,jj) = +( 500. + 0.5 * 1500. * ( 1.0 + tanh( (glamt(:,3) - 40.) / 7. ) ))
END DO
ENDIF
IF ( cp_cfg=='DOME' ) THEN
ALLOCATE(zbatv(jpi,jpj))
zbatv(:,:) = MIN(3600._wp, MAX( 600._wp, 600._wp -gphiv(:,:)*1.e3*0.01 ))
DO jj =2,jpj
bathy(:,jj) = 0.5_wp*(zbatv(:,jj) + zbatv(:,jj-1))
END DO
CALL lbc_lnk( 'zgr_bat', bathy, 'T', 1._wp )
WHERE (gphit(:,:) >0._wp) bathy(:,:) = 0._wp
! Dig inlet:
WHERE ((gphit(:,:)>0._wp).AND.(glamt(:,:)>-50._wp).AND.(glamt(:,:)<50._wp)) bathy(:,:) = 600._wp
DEALLOCATE(zbatv)
ENDIF
! CALL lbc_lnk( 'zgr_bat', bathy, 'T', 1._wp )
!
DEALLOCATE(zrand)
!
IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
mbathy(:,:) = jpkm1
ELSE ! z-coordinate (zco or zps): step-like topography
mbathy(:,:) = jpkm1
DO jk = 1, jpkm1
WHERE( gdept_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= gdept_1d(jk+1) ) mbathy(:,:) = jk
END DO
ENDIF
ENDIF
!
risfdep(:,:)=0.e0
misfdep(:,:)=1
!
! ! ================ !
ELSEIF( ntopo > 0 ) THEN ! read in file ! (over the local domain)
! ! ================ !
!
IF( ln_zco ) THEN ! zco : read level bathymetry
CALL iom_open ( cn_topolvl, inum )
CALL iom_get ( inum, jpdom_auto, cn_bathlvl, bathy )
CALL iom_close( inum )
mbathy(:,:) = INT( bathy(:,:) )
! initialisation isf variables
risfdep(:,:)=0._wp ; misfdep(:,:)=1
! ! =====================
IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
! ! =====================
!
ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
DO ji = mi0(ii0), mi1(ii1)
DO jj = mj0(ij0), mj1(ij1)
mbathy(ji,jj) = 15
END DO
END DO
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
!
ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
DO ji = mi0(ii0), mi1(ii1)
DO jj = mj0(ij0), mj1(ij1)
mbathy(ji,jj) = 12
END DO
END DO
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
!
ENDIF
!
ENDIF
IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
#if defined key_agrif
IF (agrif_root()) THEN
#endif
IF( ntopo == 1) THEN
CALL iom_open ( cn_topo, inum )
CALL iom_get ( inum, jpdom_auto, cn_bath, bathy )
CALL iom_close( inum )
ELSE
CALL dom_bat
ENDIF
#if defined key_agrif
ELSE
IF( ntopo == 1) THEN
CALL iom_open ( cn_topo, inum )
CALL iom_get ( inum, jpdom_auto, cn_bath, bathy )
CALL iom_close( inum )
ELSE
IF( ntopo == 3) THEN
CALL agrif_create_bathy_meter()
ELSE
CALL dom_bat
ENDIF
ENDIF
ENDIF
#endif
!
! initialisation isf variables
risfdep(:,:)=0._wp ; misfdep(:,:)=1
!
IF ( ln_isfcav ) THEN
CALL iom_open ( cn_fisfd, inum )
CALL iom_get ( inum, jpdom_auto, cn_visfd, risfdep )
CALL iom_close( inum )
END IF
!
IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
!
ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
DO ji = mi0(ii0), mi1(ii1)
DO jj = mj0(ij0), mj1(ij1)
bathy(ji,jj) = 284._wp
END DO
END DO
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
!
ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
DO ji = mi0(ii0), mi1(ii1)
DO jj = mj0(ij0), mj1(ij1)
bathy(ji,jj) = 137._wp
END DO
END DO
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
!
ENDIF
!
ENDIF
! ! =============== !
ELSE ! error !
! ! =============== !
WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
ENDIF
!
IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth
ENDIF
zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels
WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
ELSE WHERE ( risfdep == 0._wp ); bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
END WHERE
IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
ENDIF
#if defined key_agrif
! IF (Agrif_Root()) THEN
! IF (Agrif_level()==1) THEN
! ALLOCATE(ztemp(jpi,jpj))
! CALL smooth_bat_negative(bathy, ztemp, 0.1_wp, 10._wp)
! WHERE (bathy/=0._wp) bathy(:,:) = ztemp(:,:)
! DEALLOCATE(ztemp)
! ENDIF
IF (( .NOT.ln_sco).AND.(.NOT.Agrif_Root() )) THEN
CALL agrif_bathymetry_connect
IF ( ln_remove_closedseas ) CALL remove_closedseas
ENDIF
#endif
!
END SUBROUTINE zgr_bat
SUBROUTINE zgr_bat_ctl
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_bat_ctl ***
!!
!! ** Purpose : check the bathymetry in levels
!!
!! ** Method : The array mbathy is checked to verified its consistency
!! with the model options. in particular:
!! mbathy must have at least 1 land grid-points (mbathy<=0)
!! along closed boundary.
!! mbathy must be cyclic IF jperio=1.
!! mbathy must be lower or equal to jpk-1.
!! isolated ocean grid points are suppressed from mbathy
!! since they are only connected to remaining
!! ocean through vertical diffusion.
!! C A U T I O N : mbathy will be modified during the initializa-
!! tion phase to become the number of non-zero w-levels of a water
!! column, with a minimum value of 1.
!!
!! ** Action : - update mbathy: level bathymetry (in level index)
!! - update bathy : meter bathymetry (in meters)
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jl ! dummy loop indices
INTEGER :: icompt, ibtest, ikmax ! temporary integers
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zbathy
!!----------------------------------------------------------------------
!
ALLOCATE(zbathy(jpi,jpj))
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
! ! Suppress isolated ocean grid points
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
IF(lwp) WRITE(numout,*)' -----------------------------------'
icompt = 0
DO jl = 1, 2
IF( l_Iperio ) THEN
mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
mbathy(jpi,:) = mbathy( 2 ,:)
ENDIF
zbathy(:,:) = FLOAT( mbathy(:,:) )
CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
mbathy(:,:) = INT( zbathy(:,:) )
DO jj = 2, jpjm1
DO ji = 2, jpim1
ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
& mbathy(ji,jj-1), mbathy(ji,jj+1) )
IF( ibtest < mbathy(ji,jj) ) THEN
IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
mbathy(ji,jj) = ibtest
icompt = icompt + 1
ENDIF
END DO
END DO
END DO
IF( lk_mpp ) CALL mpp_sum( 'domzgr', icompt )
IF( icompt == 0 ) THEN
IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
ELSE
IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
ENDIF
IF( lk_mpp ) THEN
zbathy(:,:) = FLOAT( mbathy(:,:) )
CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp )
mbathy(:,:) = INT( zbathy(:,:) )
ENDIF
! ! East-west cyclic boundary conditions
IF( jperio == 0 ) THEN
IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: jperio = ', jperio
! IF( ln_zco .OR. ln_zps ) THEN
mbathy( mi0( 1+nn_hls):mi1( 1+nn_hls),:) = 0
mbathy( mi0(jpiglo-nn_hls):mi1(jpiglo-nn_hls),:) = 0
mbathy(:, mj0( 1+nn_hls):mj1( 1+nn_hls)) = 0
mbathy(:, mj0(jpjglo-nn_hls):mj1(jpjglo-nn_hls)) = 0
! ELSE
! mbathy( mi0( 1+nn_hls):mi1( 1+nn_hls),:) = jpkm1
! mbathy( mi0(jpiglo-nn_hls):mi1(jpiglo-nn_hls),:) = jpkm1
! ENDIF
ELSEIF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio
! mbathy( 1 ,:) = mbathy(jpim1,:)
! mbathy(jpi,:) = mbathy( 2 ,:)
ELSEIF( jperio == 2 ) THEN
IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: jperio = ', jperio
ELSE
IF(lwp) WRITE(numout,*) ' e r r o r'
IF(lwp) WRITE(numout,*) ' parameter , jperio = ', jperio
! STOP 'dom_mba'
ENDIF
! Boundary condition on mbathy
IF( .NOT.lk_mpp ) THEN
!!gm !!bug ??? think about it !
! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
zbathy(:,:) = FLOAT( mbathy(:,:) )
CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
mbathy(:,:) = INT( zbathy(:,:) )
ENDIF
! Number of ocean level inferior or equal to jpkm1
zbathy(:,:) = FLOAT( mbathy(:,:) )
ikmax = MAXVAL(zbathy(:,:))
CALL mpp_max( 'domzgr',ikmax)
IF( ikmax > jpkm1 ) THEN
IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
ELSE IF( ikmax < jpkm1 ) THEN
IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
ENDIF
!
DEALLOCATE( zbathy )
!
END SUBROUTINE zgr_bat_ctl
SUBROUTINE zgr_bot_level
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_bot_level ***
!!
!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
!!
!! ** Method : computes from mbathy with a minimum value of 1 over land
!!
!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
!! ocean level at t-, u- & v-points
!! (min value = 1 over land)
!!----------------------------------------------------------------------
INTEGER :: ji, jj ! dummy loop indices
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmbk
!!----------------------------------------------------------------------
!
ALLOCATE( zmbk(jpi,jpj) )
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
!
mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
! ! bottom k-index of W-level = mbkt+1
DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
DO ji = 1, jpim1
mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
mbkf(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj), mbkt(ji+1,jj ), mbkt(ji+1,jj+1) )
END DO
END DO
! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
zmbk(:,:) = REAL( mbkf(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'F',1.) ; mbkf (:,:) = MAX( INT( zmbk(:,:) ), 1 )
!
DEALLOCATE( zmbk )
!
END SUBROUTINE zgr_bot_level
SUBROUTINE zgr_top_level
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_top_level ***
!!
!! ** Purpose : defines the vertical index of ocean top (mik. arrays)
!!
!! ** Method : computes from misfdep with a minimum value of 1
!!
!! ** Action : mikt, miku, mikv : vertical indices of the shallowest
!! ocean level at t-, u- & v-points
!! (min value = 1)
!!----------------------------------------------------------------------
INTEGER :: ji, jj ! dummy loop indices
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmik
!!----------------------------------------------------------------------
!
ALLOCATE( zmik(jpi,jpj) )
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
!
mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1)
! ! top k-index of W-level (=mikt)
DO jj = 1, jpjm1 ! top k-index of U- (U-) level
DO ji = 1, jpim1
miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) )
mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) )
mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) )
END DO
END DO
! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 )
zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 )
zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 )
!
DEALLOCATE( zmik )
!
END SUBROUTINE zgr_top_level
SUBROUTINE zgr_zco
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_zco ***
!!
!! ** Purpose : define the reference z-coordinate system
!!
!! ** Method : set 3D coord. arrays to reference 1D array
!!----------------------------------------------------------------------
INTEGER :: jk
!!----------------------------------------------------------------------
!
DO jk = 1, jpk
gdept_0(:,:,jk) = gdept_1d(jk)
gdepw_0(:,:,jk) = gdepw_1d(jk)
e3t_0 (:,:,jk) = e3t_1d (jk)
e3u_0 (:,:,jk) = e3t_1d (jk)
e3v_0 (:,:,jk) = e3t_1d (jk)
e3f_0 (:,:,jk) = e3t_1d (jk)
e3w_0 (:,:,jk) = e3w_1d (jk)
e3uw_0 (:,:,jk) = e3w_1d (jk)
e3vw_0 (:,:,jk) = e3w_1d (jk)
END DO
!
END SUBROUTINE zgr_zco
SUBROUTINE zgr_zps
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_zps ***
!!
!! ** Purpose : the depth and vertical scale factor in partial step
!! reference z-coordinate case
!!
!! ** Method : Partial steps : computes the 3D vertical scale factors
!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
!! a partial step representation of bottom topography.
!!
!! The reference depth of model levels is defined from an analytical
!! function the derivative of which gives the reference vertical
!! scale factors.
!! From depth and scale factors reference, we compute there new value
!! with partial steps on 3d arrays ( i, j, k ).
!!
!! w-level: gdepw_0(i,j,k) = gdep(k)
!! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k)
!! t-level: gdept_0(i,j,k) = gdep(k+0.5)
!! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
!!
!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
!! we find the mbathy index of the depth at each grid point.
!! This leads us to three cases:
!!
!! - bathy = 0 => mbathy = 0
!! - 1 < mbathy < jpkm1
!! - bathy > gdepw_0(jpk) => mbathy = jpkm1
!!
!! Then, for each case, we find the new depth at t- and w- levels
!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
!! and f-points.
!!
!! This routine is given as an example, it must be modified
!! following the user s desiderata. nevertheless, the output as
!! well as the way to compute the model levels and scale factors
!! must be respected in order to insure second order accuracy
!! schemes.
!!
!! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
!! - - - - - - - gdept_0, gdepw_0 and e3. are positives
!!
!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ik, it, ikb, ikt ! temporary integers
REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
REAL(wp) :: zdiff ! temporary scalar
REAL(wp) :: zmax ! temporary scalar
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zprt
!!---------------------------------------------------------------------
!
ALLOCATE( zprt(jpi,jpj,jpk) )
!
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
! compute position of the ice shelf grounding line
! set bathy and isfdraft to 0 where grounded
IF ( ln_isfcav ) CALL zgr_isf_zspace
! bathymetry in level (from bathy_meter)
! ===================
zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
END WHERE
! Compute mbathy for ocean points (i.e. the number of ocean levels)
! find the number of ocean levels such that the last level thickness
! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
! e3t_1d is the reference level thickness
DO jk = jpkm1, 1, -1
zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
END DO
! Check compatibility between bathy and iceshelf draft
! insure at least 2 wet level on the vertical under an ice shelf
! compute misfdep and adjust isf draft if needed
IF ( ln_isfcav ) CALL zgr_isf_kspace
! Scale factors and depth at T- and W-points
DO jk = 1, jpk ! intitialization to the reference z-coordinate
gdept_0(:,:,jk) = gdept_1d(jk)
gdepw_0(:,:,jk) = gdepw_1d(jk)
e3t_0 (:,:,jk) = e3t_1d (jk)
e3w_0 (:,:,jk) = e3w_1d (jk)
END DO
! Scale factors and depth at T- and W-points
DO jj = 1, jpj
DO ji = 1, jpi
ik = mbathy(ji,jj)
IF( ik > 0 ) THEN ! ocean point only
! max ocean level case
IF( ik == jpkm1 ) THEN
zdepwp = bathy(ji,jj)
ze3tp = bathy(ji,jj) - gdepw_1d(ik)
ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
e3t_0(ji,jj,ik ) = ze3tp
e3t_0(ji,jj,ik+1) = ze3tp
IF ( ln_e3_dep.AND.ln_dept_mid ) THEN
gdept_0(ji,jj,ik) = gdepw_1d(ik) + 0.5_wp * ze3tp
e3w_0(ji,jj,ik ) = gdept_0(ji,jj,ik) - gdept_0(ji,jj,ik-1)
ELSE
gdept_0(ji,jj,ik) = gdept_1d(ik-1) + ze3wp
e3w_0(ji,jj,ik ) = ze3wp
ENDIF
gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
e3w_0(ji,jj,ik+1) = ze3tp
gdepw_0(ji,jj,ik+1) = zdepwp
!
ELSE ! standard case
IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
ENDIF
!gm Bug? check the gdepw_1d
! ... on ik
e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
& / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
IF ( ln_e3_dep.AND.ln_dept_mid ) THEN
gdept_0(ji,jj,ik) = gdepw_1d(ik) + 0.5_wp * e3t_0(ji,jj,ik)
e3w_0(ji,jj,ik) = gdept_0(ji,jj,ik) - gdept_0(ji,jj,ik-1)
ELSE
gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
& * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
ENDIF
! ... on ik+1
e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
ENDIF
ENDIF
END DO
END DO
!
it = 0
DO jj = 1, jpj
DO ji = 1, jpi
ik = mbathy(ji,jj)
IF( ik > 0 ) THEN ! ocean point only
e3tp (ji,jj) = e3t_0(ji,jj,ik)
e3wp (ji,jj) = e3w_0(ji,jj,ik)
! test
zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
IF( zdiff <= 0._wp .AND. lwp ) THEN
it = it + 1
WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
WRITE(numout,*) ' bathy = ', bathy(ji,jj)
WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
ENDIF
ENDIF
END DO
END DO
!
! compute top scale factor if ice shelf
IF (ln_isfcav) CALL zps_isf
!
! Scale factors and depth at U-, V-, UW and VW-points
DO jk = 1, jpk ! initialisation to z-scale factors
e3u_0 (:,:,jk) = e3t_1d(jk)
e3v_0 (:,:,jk) = e3t_1d(jk)
e3uw_0(:,:,jk) = e3w_1d(jk)
e3vw_0(:,:,jk) = e3w_1d(jk)
END DO
DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
DO jj = 1, jpjm1
DO ji = 1, jpim1 ! vector opt.
e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
END DO
END DO
END DO
! update e3uw in case only 2 cells in the water column
IF ( ln_isfcav ) CALL zps_isf_e3uv_w
!
CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp ) ! lateral boundary conditions
CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
!
DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)
WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)
WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)
WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)
END DO
! Scale factor at F-point
DO jk = 1, jpk ! initialisation to z-scale factors
e3f_0(:,:,jk) = e3t_1d(jk)
END DO
DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
DO jj = 1, jpjm1
DO ji = 1, jpim1 ! vector opt.
e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
END DO
END DO
END DO
CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp ) ! Lateral boundary conditions
!
DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)
END DO
!!gm bug ? : must be a do loop with mj0,mj1
!
e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)
e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)
e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)
e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)
! Control of the sign
IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )
IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )
IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )
IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )
!
! if in the future gde3w_0 need to be compute, use the function defined in NEMO
! for now gde3w_0 computation is removed as not an output of domcfg
DEALLOCATE( zprt )
!
END SUBROUTINE zgr_zps
SUBROUTINE zgr_sco
!!----------------------------------------------------------------------
!! *** ROUTINE zgr_sco ***
!!
!! ** Purpose : define the s-coordinate system
!!
!! ** Method : s-coordinate
!! The depth of model levels is defined as the product of an
!! analytical function by the local bathymetry, while the vertical
!! scale factors are defined as the product of the first derivative
!! of the analytical function by the bathymetry.
!! (this solution save memory as depth and scale factors are not
!! 3d fields)
!! - Read bathymetry (in meters) at t-point and compute the
!! bathymetry at u-, v-, and f-points.
!! hbatu = mi( hbatt )
!! hbatv = mj( hbatt )
!! hbatf = mi( mj( hbatt ) )
!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
!! function and its derivative given as function.
!! z_gsigt(k) = fssig (k )
!! z_gsigw(k) = fssig (k-0.5)
!! z_esigt(k) = fsdsig(k )
!! z_esigw(k) = fsdsig(k-0.5)
!! Three options for stretching are give, and they can be modified
!! following the users requirements. Nevertheless, the output as
!! well as the way to compute the model levels and scale factors
!! must be respected in order to insure second order accuracy
!! schemes.
!!
!! The three methods for stretching available are:
!!
!! s_sh94 (Song and Haidvogel 1994)
!! a sinh/tanh function that allows sigma and stretched sigma
!!
!! s_sf12 (Siddorn and Furner 2012?)
!! allows the maintenance of fixed surface and or
!! bottom cell resolutions (cf. geopotential coordinates)
!! within an analytically derived stretched S-coordinate framework.
!!
!! s_tanh (Madec et al 1996)
!! a cosh/tanh function that gives stretched coordinates
!!
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk, jl ! dummy loop argument
INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
INTEGER :: ios ! Local integer output status for namelist read
REAL(wp) :: zrmax, ztaper ! temporary scalars
REAL(wp) :: zrfact
!
REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
!!
NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
& rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
!!----------------------------------------------------------------------
!
ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj) )
!
!REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist')
!REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist')
IF(lwm) WRITE ( numond, namzgr_sco )
IF(lwp) THEN ! control print
WRITE(numout,*)
WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
WRITE(numout,*) '~~~~~~~~~~~'
WRITE(numout,*) ' Namelist namzgr_sco'
WRITE(numout,*) ' stretching coeffs '
WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
WRITE(numout,*) ' Bottom cell (Zb) (m) = H*rn_zb_a + rn_zb_b'
ENDIF
hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
hifu(:,:) = rn_sbot_min
hifv(:,:) = rn_sbot_min
hiff(:,:) = rn_sbot_min
! ! set maximum ocean depth
#if defined key_agrif
! Set boundary conditions prior filtering from parent:
IF (.NOT.Agrif_Root()) THEN
CALL agrif_bathymetry_connect
IF ( ln_remove_closedseas ) CALL remove_closedseas
ENDIF
#endif
bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
DO jj = 1, jpj
DO ji = 1, jpi
IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
END DO
END DO
! ! =============================
! ! Define the envelop bathymetry (hbatt)
! ! =============================
!
! smooth the bathymetry (if required)
scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
scobot(:,:) = bathy(:,:) ! ocean bottom depth
!
! Envelope bathymetry saved in hbatt
CALL smooth_bat_negative(bathy, hbatt, rn_rmax, rn_sbot_min)
IF ((ntopo>0).AND.MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
DO jj = 1, jpj
DO ji = 1, jpi
ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
END DO
END DO
ENDIF
#if defined key_agrif
IF (.NOT.Agrif_Root()) THEN
WHERE (bathy(:,:)>0._wp ) bathy(:,:) = hbatt(:,:)
CALL agrif_bathymetry_connect
IF ( ln_remove_closedseas ) CALL remove_closedseas
! Discard the enveloppe in the agrif case
hbatt(:,:) = bathy(:,:)
WHERE ( bathy(:,:) == 0._wp ) hbatt(:,:) = rn_sbot_min
scobot(:,:) = bathy(:,:)
ENDIF
#endif
!
! ! ==============================
! ! hbatu, hbatv, hbatf fields
! ! ==============================
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
ENDIF
hbatu(:,:) = rn_sbot_min
hbatv(:,:) = rn_sbot_min
hbatf(:,:) = rn_sbot_min
DO jj = 1, jpjm1
DO ji = 1, jpim1 ! NO vector opt.
hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
END DO
END DO
!
! Apply lateral boundary condition
!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp )
DO jj = 1, jpj
DO ji = 1, jpi
IF( hbatu(ji,jj) == 0._wp ) THEN
!No worries about the following line when ln_wd == .true.
IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
ENDIF
END DO
END DO
zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp )
DO jj = 1, jpj
DO ji = 1, jpi
IF( hbatv(ji,jj) == 0._wp ) THEN
IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
ENDIF
END DO
END DO
zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp )
DO jj = 1, jpj
DO ji = 1, jpi
IF( hbatf(ji,jj) == 0._wp ) THEN
IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
ENDIF
END DO
END DO
!!bug: key_helsinki a verifer
hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
IF( lwp ) THEN
WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
ENDIF
!! helsinki
! ! =======================
! ! s-ccordinate fields (gdep., e3.)
! ! =======================
!
! non-dimensional "sigma" for model level depth at w- and t-levels
!========================================================================
! Song and Haidvogel 1994 (ln_s_sh94=T)
! Siddorn and Furner 2012 (ln_sf12=T)
! or tanh function (both false)
!========================================================================
IF ( ln_s_sh94 ) THEN
CALL s_sh94()
ELSE IF ( ln_s_sf12 ) THEN
CALL s_sf12()
ELSE
CALL s_tanh()
ENDIF
CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp )
CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp )
CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp )
CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp )
CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp )
CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp )
CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
!
WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp
WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp
WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp
WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp
WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp
WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp
WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp
!!
! HYBRID :
DO jj = 1, jpj
DO ji = 1, jpi
DO jk = 1, jpkm1
IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
END DO
END DO
END DO
WHERE (bathy(:,:)<=0) mbathy(:,:) = 0
IF(lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
& ' MAX ', MAXVAL( mbathy(:,:) )
IF( lwp ) THEN ! min max values over the local domain
WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
& ' w ', MINVAL( gdepw_0(:,:,:) )
WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
& ' w ', MINVAL( e3w_0 (:,:,:) )
WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
& ' w ', MAXVAL( gdepw_0(:,:,:) )
WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
& ' w ', MAXVAL( e3w_0 (:,:,:) )
ENDIF
! END DO
IF(lwp) THEN ! selected vertical profiles
WRITE(numout,*)
WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
WRITE(numout,*) ' ~~~~~~ --------------------'
WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
DO jj = mj0(20), mj1(20)
DO ji = mi0(20), mi1(20)
WRITE(numout,*)
WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
WRITE(numout,*) ' ~~~~~~ --------------------'
WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
END DO
END DO
! DO jj = mj0(74), mj1(74)
! DO ji = mi0(10), mi1(10)
! WRITE(numout,*)
! WRITE(numout,*) ' domzgr: vertical coordinates : point (10,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
! WRITE(numout,*) ' ~~~~~~ --------------------'
! WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
! WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
! & e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
! END DO
! END DO
ENDIF
!
!================================================================================
! check the coordinate makes sense
!================================================================================
DO ji = 1, jpi
DO jj = 1, jpj
!
IF( hbatt(ji,jj) > 0._wp) THEN
DO jk = 1, mbathy(ji,jj)
! check coordinate is monotonically increasing
IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
CALL ctl_stop( ctmp1 )
ENDIF
! and check it has never gone negative
IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
CALL ctl_stop( ctmp1 )
ENDIF
! and check it never exceeds the total depth
IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
CALL ctl_stop( ctmp1 )
ENDIF
END DO
!
DO jk = 1, mbathy(ji,jj)-1
! and check it never exceeds the total depth
IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
CALL ctl_stop( ctmp1 )
ENDIF
END DO
ENDIF
END DO
END DO
!
DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
!
END SUBROUTINE zgr_sco
SUBROUTINE s_sh94()
!!----------------------------------------------------------------------
!! *** ROUTINE s_sh94 ***
!!
!! ** Purpose : stretch the s-coordinate system
!!
!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
!! mixed S/sigma coordinate
!!
!! Reference : Song and Haidvogel 1994.
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop argument
REAL(wp) :: zcoeft, zcoefw ! temporary scalars
REAL(wp) :: ztmpu, ztmpv, ztmpf
REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
!
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
!!----------------------------------------------------------------------
ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) )
ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp
z_esigt3 = 0._wp ; z_esigw3 = 0._wp
z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
!
DO ji = 1, jpi
DO jj = 1, jpj
!
IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
DO jk = 1, jpk
z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
END DO
ELSE ! shallow water, uniform sigma
DO jk = 1, jpk
z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
END DO
ENDIF
IF ( ln_e3_dep.AND.ln_dept_mid ) THEN
DO jk = 1, jpkm1
z_gsigt3(ji,jj,jk) = 0.5_wp * (z_gsigw3(ji,jj,jk) +z_gsigw3(ji,jj,jk+1))
END DO
ENDIF
!
DO jk = 1, jpkm1
z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
END DO
z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
!
DO jk = 1, jpk
zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft )
gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw )
END DO
!
END DO ! for all jj's
END DO ! for all ji's
DO ji = 1, jpim1
DO jj = 1, jpjm1
! extended for Wetting/Drying case
ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
DO jk = 1, jpk
z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
& / ztmpu
z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
& / ztmpu
z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
& / ztmpv
z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
& / ztmpv
z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
& + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
& + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
& + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
!
e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
!
e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
END DO
END DO
END DO
!
DEALLOCATE( z_gsigw3, z_gsigt3 )
DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
!
END SUBROUTINE s_sh94
SUBROUTINE s_sf12
!!----------------------------------------------------------------------
!! *** ROUTINE s_sf12 ***
!!
!! ** Purpose : stretch the s-coordinate system
!!
!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
!! mixed S/sigma/Z coordinate
!!
!! This method allows the maintenance of fixed surface and or
!! bottom cell resolutions (cf. geopotential coordinates)
!! within an analytically derived stretched S-coordinate framework.
!!
!!
!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop argument
REAL(wp) :: zsmth ! smoothing around critical depth
REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
REAL(wp) :: ztmpu, ztmpv, ztmpf
REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
!
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
!!----------------------------------------------------------------------
!
ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk))
ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp
z_esigt3 = 0._wp ; z_esigw3 = 0._wp
z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
DO ji = 1, jpi
DO jj = 1, jpj
IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
! could be changed by users but care must be taken to do so carefully
zzb = 1.0_wp-(zzb/hbatt(ji,jj))
zzs = rn_zs / hbatt(ji,jj)
IF (rn_efold /= 0.0_wp) THEN
zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
ELSE
zsmth = 1.0_wp
ENDIF
DO jk = 1, jpk
z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
ENDDO
z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
DO jk = 1, jpk
z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
END DO
ELSE ! shallow water, z coordinates
DO jk = 1, jpk
z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
END DO
ENDIF
IF ( ln_e3_dep.AND.ln_dept_mid ) THEN
DO jk = 1, jpkm1
z_gsigt3(ji,jj,jk) = 0.5_wp * (z_gsigw3(ji,jj,jk) +z_gsigw3(ji,jj,jk+1))
END DO
ENDIF
DO jk = 1, jpkm1
z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
END DO
z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
DO jk = 1, jpk
gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
END DO
ENDDO ! for all jj's
ENDDO ! for all ji's
DO ji=1,jpi-1
DO jj=1,jpj-1
! extend to suit for Wetting/Drying case
ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
DO jk = 1, jpk
z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
& / ztmpu
z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
& / ztmpu
z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
& / ztmpv
z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
& / ztmpv
z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
& + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
& + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
& + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
! Code prior to wetting and drying option (for reference)
!z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigt3(ji+1,jj,jk) ) &
! /( hbatt(ji,jj)+hbatt(ji+1,jj) )
!
!z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)*z_esigw3(ji+1,jj,jk) ) &
! /( hbatt(ji,jj)+hbatt(ji+1,jj) )
!
!z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigt3(ji,jj+1,jk) ) &
! /( hbatt(ji,jj)+hbatt(ji,jj+1) )
!
!z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)*z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)*z_esigw3(ji,jj+1,jk) ) &
! /( hbatt(ji,jj)+hbatt(ji,jj+1) )
!
!z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
! & +hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
! +hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
! & +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) &
! /( hbatt(ji ,jj )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
!
e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
END DO
ENDDO
ENDDO
!
CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.)
CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.)
CALL lbc_lnk('domzgr',e3w_0 ,'T',1.)
CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.)
!
DEALLOCATE( z_gsigw3, z_gsigt3 )
DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
!
END SUBROUTINE s_sf12
SUBROUTINE s_tanh()
!!----------------------------------------------------------------------
!! *** ROUTINE s_tanh***
!!
!! ** Purpose : stretch the s-coordinate system
!!
!! ** Method : s-coordinate stretch
!!
!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop argument
REAL(wp) :: zcoeft, zcoefw ! temporary scalars
REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt
REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
!!----------------------------------------------------------------------
ALLOCATE( z_gsigw(jpk), z_gsigt(jpk) )
ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
z_gsigw = 0._wp ; z_gsigt = 0._wp
z_esigt = 0._wp ; z_esigw = 0._wp
DO jk = 1, jpk
z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
z_gsigt(jk) = -fssig( REAL(jk,wp) )
END DO
IF ( ln_e3_dep.AND.ln_dept_mid ) THEN
DO jk = 1, jpkm1
z_gsigt(jk) = 0.5_wp * (z_gsigw(jk)+z_gsigw(jk+1))
END DO
ENDIF
IF( lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
!
! Coefficients for vertical scale factors at w-, t- levels
!!gm bug : define it from analytical function, not like juste bellow....
!!gm or betteroffer the 2 possibilities....
DO jk = 1, jpkm1
z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
END DO
z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
!
DO jk = 1, jpk
zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft )
gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw )
END DO
DO jj = 1, jpj
DO ji = 1, jpi
DO jk = 1, jpk
e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
!
e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
END DO
END DO
END DO
!
DEALLOCATE( z_gsigw, z_gsigt )
DEALLOCATE( z_esigt, z_esigw )
!
END SUBROUTINE s_tanh
FUNCTION fssig( pk ) RESULT( pf )
!!----------------------------------------------------------------------
!! *** ROUTINE fssig ***
!!
!! ** Purpose : provide the analytical function in s-coordinate
!!
!! ** Method : the function provide the non-dimensional position of
!! T and W (i.e. between 0 and 1)
!! T-points at integer values (between 1 and jpk)
!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
!!----------------------------------------------------------------------
REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
REAL(wp) :: pf ! sigma value
!!----------------------------------------------------------------------
!
pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
& - TANH( rn_thetb * rn_theta ) ) &
& * ( COSH( rn_theta ) &
& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
& / ( 2._wp * SINH( rn_theta ) )
!
END FUNCTION fssig
FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
!!----------------------------------------------------------------------
!! *** ROUTINE fssig1 ***
!!
!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
!!
!! ** Method : the function provides the non-dimensional position of
!! T and W (i.e. between 0 and 1)
!! T-points at integer values (between 1 and jpk)
!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
!!----------------------------------------------------------------------
REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
REAL(wp) :: pf1 ! sigma value
!!----------------------------------------------------------------------
!
IF ( rn_theta == 0 ) then ! uniform sigma
pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
ELSE ! stretched sigma
pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
& + pbb * ( (TANH( rn_theta*( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp * rn_theta ) ) &
& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
ENDIF
!
END FUNCTION fssig1
FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
!!----------------------------------------------------------------------
!! *** ROUTINE fgamma ***
!!
!! ** Purpose : provide analytical function for the s-coordinate
!!
!! ** Method : the function provides the non-dimensional position of
!! T and W (i.e. between 0 and 1)
!! T-points at integer values (between 1 and jpk)
!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
!!
!! This method allows the maintenance of fixed surface and or
!! bottom cell resolutions (cf. geopotential coordinates)
!! within an analytically derived stretched S-coordinate framework.
!!
!! Reference : Siddorn and Furner, in prep
!!----------------------------------------------------------------------
REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
REAL(wp) :: p_gamma(jpk) ! stretched coordinate
REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
REAL(wp), INTENT(in ) :: pzs ! surface box depth
REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
!
INTEGER :: jk ! dummy loop index
REAL(wp) :: za1,za2,za3 ! local scalar
REAL(wp) :: zn1,zn2 ! - -
REAL(wp) :: za,zb,zx ! - -
!!----------------------------------------------------------------------
!
zn1 = 1._wp / REAL( jpkm1, wp )
zn2 = 1._wp - zn1
!
za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp)
za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp)
za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1)
!
za = pzb - za3*(pzs-za1)-za2
za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) )
zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1)
zx = 1.0_wp-za/2.0_wp-zb
!
DO jk = 1, jpk
p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + &
& zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &
& (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )
p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth)
END DO
!
END FUNCTION fgamma
SUBROUTINE smooth_bat_negative(batin, batout, rmax, batmin)
!!----------------------------------------------------------------------
!! *** ROUTINE smooth_bat_negative ***
!!
!! ** Purpose : Smooth bathymetry according to rmax factor
!! such that the result is greater than the input
!! bathymetry
!!
!! ** Method :
!!
!! Reference :
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: batin
REAL(wp), DIMEnSION(jpi,jpj), INTENT(out) :: batout
REAL(wp) , INTENT(in) :: rmax, batmin
INTEGER :: ji, jj, jl ! dummy loop argument
INTEGER :: iip1, ijp1, iim1, ijm1
REAL(wp) :: zrmax, zrfact
REAL(wp), DIMENSION(jpi,jpj) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2, zri, zrj
!!----------------------------------------------------------------------
!
batout(:,:) = batin(:,:)
!
! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
DO jj = 1, jpj
DO ji = 1, jpi
IF( batin(ji,jj) == 0._wp ) THEN
iip1 = MIN( ji+1, jpi )
ijp1 = MIN( jj+1, jpj )
iim1 = MAX( ji-1, 1 )
ijm1 = MAX( jj-1, 1 )
IF( ( + batin(iim1,ijm1) + batin(ji,ijp1) + batin(iip1,ijp1) &
& + batin(iim1,jj ) + batin(iip1,jj ) &
& + batin(iim1,ijm1) + batin(ji,ijm1) + batin(iip1,ijp1) ) > 0._wp ) &
& batout(ji,jj) = batmin
ENDIF
END DO
END DO
!
! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
CALL lbc_lnk( 'domzgr',batout, 'T', 1._wp, kfillmode=jpfillnothing )
!
jl = 0
zrmax = 1._wp
!
! set scaling factor used in reducing vertical gradients
zrfact = ( 1._wp - rmax ) / ( 1._wp + rmax )
!
! initialise temporary evelope depth arrays
ztmpi1(:,:) = batout(:,:)
ztmpi2(:,:) = batout(:,:)
ztmpj1(:,:) = batout(:,:)
ztmpj2(:,:) = batout(:,:)
!
! initialise temporary r-value arrays
zri(:,:) = 1._wp
zrj(:,:) = 1._wp
! ! ================ !
DO WHILE( jl <= 10000 .AND. ( zrmax - rmax ) > 1.e-8_wp ) ! Iterative loop !
! ! ================ !
jl = jl + 1
zrmax = 0._wp
! we set zrmax from previous r-values (zri and zrj) first
! if set after current r-value calculation (as previously)
! we could exit DO WHILE prematurely before checking r-value
! of current batout
DO_2D( 0, 0, 0, 0 )
zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
END_2D
zri(:,:) = 0._wp
zrj(:,:) = 0._wp
DO_2D( 0, 0, 0, 0 )
iip1 = MIN( ji+1, jpi ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
ijp1 = MIN( jj+1, jpj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
IF( (batout(ji,jj) > 0._wp) .AND. (batout(iip1,jj) > 0._wp)) THEN
zri(ji,jj) = ( batout(iip1,jj ) - batout(ji,jj) ) / ( batout(iip1,jj ) + batout(ji,jj) )
END IF
IF( (batout(ji,jj) > 0._wp) .AND. (batout(ji,ijp1) > 0._wp)) THEN
zrj(ji,jj) = ( batout(ji ,ijp1) - batout(ji,jj) ) / ( batout(ji ,ijp1) + batout(ji,jj) )
END IF
IF( zri(ji,jj) > rmax ) ztmpi1(ji ,jj ) = batout(iip1,jj ) * zrfact
IF( zri(ji,jj) < -rmax ) ztmpi2(iip1,jj ) = batout(ji ,jj ) * zrfact
IF( zrj(ji,jj) > rmax ) ztmpj1(ji ,jj ) = batout(ji ,ijp1) * zrfact
IF( zrj(ji,jj) < -rmax ) ztmpj2(ji ,ijp1) = batout(ji ,jj ) * zrfact
END_2D
!
IF(lwp)WRITE(numout,*) 'smooth_bat_negative: iter= ',jl, ' rmax= ', zrmax
!
DO_2D( 0, 0, 0, 0 )
batout(ji,jj) = MAX(batout(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
END_2D
! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
CALL lbc_lnk( 'smooth_bat_negative',batout, 'T', 1._wp, kfillmode=jpfillnothing)
! ! ================ !
END DO ! End loop !
! ! ================ !
DO jj = 1, jpj
DO ji = 1, jpi
batout(ji,jj) = MAX( batout(ji,jj), batmin ) ! set all points to avoid undefined scale value warnings
END DO
END DO
!
END SUBROUTINE smooth_bat_negative
!!======================================================================
END MODULE domzgr
Reference in New Issue View Git Blame Copy Permalink