New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

domzgr.F90 in branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2012/dev_UKMO_2012/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 3618

Last change on this file since 3618 was 3618, checked in by johnsiddorn, 11 years ago
updates to reflect reviewers comments on coding standards and documentation
Property svn:keywords set to `Id`
File size: 94.9 KB

Rev	Line
[3]	1	MODULE domzgr
	2	!!==============================================================================
	3	!! * MODULE domzgr *
	4	!! Ocean initialization : domain initialization
	5	!!==============================================================================
[1566]	6	!! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate
	7	!! ! 1997-07 (G. Madec) lbc_lnk call
	8	!! ! 1997-04 (J.-O. Beismann)
[2528]	9	!! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module
	10	!! - ! 2002-09 (A. de Miranda) rigid-lid + islands
[1566]	11	!! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module
	12	!! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function
	13	!! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco
	14	!! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style
	15	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
[2528]	16	!! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level
[3618]	17	!! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function
[1099]	18	!!----------------------------------------------------------------------
[3]	19
	20	!!----------------------------------------------------------------------
[1099]	21	!! dom_zgr : defined the ocean vertical coordinate system
[3]	22	!! zgr_bat : bathymetry fields (levels and meters)
	23	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
	24	!! zgr_bat_ctl : check the bathymetry files
[2528]	25	!! zgr_bot_level: deepest ocean level for t-, u, and v-points
[3]	26	!! zgr_z : reference z-coordinate
[454]	27	!! zgr_zco : z-coordinate
[3]	28	!! zgr_zps : z-coordinate with partial steps
[454]	29	!! zgr_sco : s-coordinate
[3618]	30	!! fssig : tanh stretch function
	31	!! fssig1 : Song and Haidvogel 1994 stretch function
	32	!! fgamma : Siddorn and Furner 2012 stretching function
[3]	33	!!---------------------------------------------------------------------
[2528]	34	USE oce ! ocean variables
	35	USE dom_oce ! ocean domain
	36	USE closea ! closed seas
	37	USE c1d ! 1D vertical configuration
	38	USE in_out_manager ! I/O manager
	39	USE iom ! I/O library
	40	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	41	USE lib_mpp ! distributed memory computing library
[3294]	42	USE wrk_nemo ! Memory allocation
	43	USE timing ! Timing
[3]	44
	45	IMPLICIT NONE
	46	PRIVATE
	47
[2715]	48	PUBLIC dom_zgr ! called by dom_init.F90
[3]	49
[2528]	50	! !!* Namelist namzgr_sco *
[3600]	51	LOGICAL :: ln_s_sh94 = .false. ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
	52	LOGICAL :: ln_s_sf12 = .true. ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
	53	!
[2528]	54	REAL(wp) :: rn_sbot_min = 300._wp ! minimum depth of s-bottom surface (>0) (m)
	55	REAL(wp) :: rn_sbot_max = 5250._wp ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
[3600]	56	REAL(wp) :: rn_rmax = 0.15_wp ! maximum cut-off r-value allowed (0<rn_rmax<1)
	57	REAL(wp) :: rn_hc = 150._wp ! Critical depth for transition from sigma to stretched coordinates
	58	! Song and Haidvogel 1994 stretching parameters
[2528]	59	REAL(wp) :: rn_theta = 6.00_wp ! surface control parameter (0<=rn_theta<=20)
	60	REAL(wp) :: rn_thetb = 0.75_wp ! bottom control parameter (0<=rn_thetb<= 1)
[3600]	61	REAL(wp) :: rn_bb = 0.80_wp ! stretching parameter
[2528]	62	! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
[3600]	63	! Siddorn and Furner stretching parameters
[3618]	64	LOGICAL :: ln_sigcrit = .false. ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
[3600]	65	REAL(wp) :: rn_alpha = 4.4_wp ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
	66	REAL(wp) :: rn_efold = 0.0_wp ! efold length scale for transition to stretched coord
	67	REAL(wp) :: rn_zs = 1.0_wp ! depth of surface grid box
	68	! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
	69	REAL(wp) :: rn_zb_a = 0.024_wp ! bathymetry scaling factor for calculating Zb
	70	REAL(wp) :: rn_zb_b = -0.2_wp ! offset for calculating Zb
[2715]	71
	72	!! * Substitutions
[3]	73	# include "domzgr_substitute.h90"
	74	# include "vectopt_loop_substitute.h90"
	75	!!----------------------------------------------------------------------
[2715]	76	!! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
[1146]	77	!! $Id$
[2528]	78	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	79	!!----------------------------------------------------------------------
	80	CONTAINS
	81
	82	SUBROUTINE dom_zgr
	83	!!----------------------------------------------------------------------
	84	!! * ROUTINE dom_zgr *
	85	!!
	86	!! ** Purpose : set the depth of model levels and the resulting
	87	!! vertical scale factors.
	88	!!
[1099]	89	!! ** Method : - reference 1D vertical coordinate (gdep._0, e3._0)
	90	!! - read/set ocean depth and ocean levels (bathy, mbathy)
	91	!! - vertical coordinate (gdep., e3.) depending on the
	92	!! coordinate chosen :
[2528]	93	!! ln_zco=T z-coordinate
[1099]	94	!! ln_zps=T z-coordinate with partial steps
	95	!! ln_zco=T s-coordinate
[3]	96	!!
[1099]	97	!! ** Action : define gdep., e3., mbathy and bathy
	98	!!----------------------------------------------------------------------
	99	INTEGER :: ioptio = 0 ! temporary integer
[2528]	100	!
[1601]	101	NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
[3]	102	!!----------------------------------------------------------------------
[3294]	103	!
	104	IF( nn_timing == 1 ) CALL timing_start('dom_zgr')
	105	!
[2528]	106	REWIND( numnam ) ! Read Namelist namzgr : vertical coordinate'
	107	READ ( numnam, namzgr )
[454]	108
[1099]	109	IF(lwp) THEN ! Control print
[454]	110	WRITE(numout,*)
	111	WRITE(numout,*) 'dom_zgr : vertical coordinate'
	112	WRITE(numout,*) '~~~~~~~'
[1601]	113	WRITE(numout,*) ' Namelist namzgr : set vertical coordinate'
[454]	114	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
	115	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
	116	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
	117	ENDIF
	118
[1099]	119	ioptio = 0 ! Check Vertical coordinate options
[454]	120	IF( ln_zco ) ioptio = ioptio + 1
	121	IF( ln_zps ) ioptio = ioptio + 1
	122	IF( ln_sco ) ioptio = ioptio + 1
[2528]	123	IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
	124	!
[3]	125	! Build the vertical coordinate system
	126	! ------------------------------------
[2528]	127	CALL zgr_z ! Reference z-coordinate system (always called)
	128	CALL zgr_bat ! Bathymetry fields (levels and meters)
	129	IF( ln_zco ) CALL zgr_zco ! z-coordinate
	130	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
	131	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
[2465]	132	!
[2528]	133	! final adjustment of mbathy & check
	134	! -----------------------------------
	135	IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain
	136	IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isoated ocean points
	137	CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
	138	!
	139	!
[2715]	140
[1348]	141	IF( nprint == 1 .AND. lwp ) THEN
	142	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	143	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
	144	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
	145	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ), &
	146	& ' u ', MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ), &
	147	& ' uw', MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)), &
	148	& ' w ', MINVAL( fse3w(:,:,:) )
	149
	150	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
	151	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
	152	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ), &
	153	& ' u ', MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ), &
	154	& ' uw', MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)), &
	155	& ' w ', MAXVAL( fse3w(:,:,:) )
	156	ENDIF
[2528]	157	!
[3294]	158	IF( nn_timing == 1 ) CALL timing_stop('dom_zgr')
	159	!
[3]	160	END SUBROUTINE dom_zgr
	161
	162
	163	SUBROUTINE zgr_z
	164	!!----------------------------------------------------------------------
	165	!! * ROUTINE zgr_z *
	166	!!
	167	!! ** Purpose : set the depth of model levels and the resulting
	168	!! vertical scale factors.
	169	!!
	170	!! ** Method : z-coordinate system (use in all type of coordinate)
	171	!! The depth of model levels is defined from an analytical
	172	!! function the derivative of which gives the scale factors.
	173	!! both depth and scale factors only depend on k (1d arrays).
[454]	174	!! w-level: gdepw_0 = fsdep(k)
	175	!! e3w_0(k) = dk(fsdep)(k) = fse3(k)
	176	!! t-level: gdept_0 = fsdep(k+0.5)
	177	!! e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
[3]	178	!!
[454]	179	!! ** Action : - gdept_0, gdepw_0 : depth of T- and W-point (m)
[1099]	180	!! - e3t_0 , e3w_0 : scale factors at T- and W-levels (m)
[3]	181	!!
[1099]	182	!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
[3]	183	!!----------------------------------------------------------------------
	184	INTEGER :: jk ! dummy loop indices
	185	REAL(wp) :: zt, zw ! temporary scalars
[1099]	186	REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
	187	REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
[1577]	188	REAL(wp) :: zrefdep ! depth of the reference level (~10m)
[2528]	189	REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
[3]	190	!!----------------------------------------------------------------------
[3294]	191	!
	192	IF( nn_timing == 1 ) CALL timing_start('zgr_z')
	193	!
[3]	194	! Set variables from parameters
	195	! ------------------------------
	196	zkth = ppkth ; zacr = ppacr
	197	zdzmin = ppdzmin ; zhmax = pphmax
[2528]	198	zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters
[3]	199
	200	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
	201	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
[1099]	202	IF( ppa1 == pp_to_be_computed .AND. &
[3]	203	& ppa0 == pp_to_be_computed .AND. &
	204	& ppsur == pp_to_be_computed ) THEN
[1099]	205	!
	206	za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
	207	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
	208	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
	209	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
	210	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
	211	ELSE
[3]	212	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
[2528]	213	za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter
[1099]	214	ENDIF
[3]	215
[1099]	216	IF(lwp) THEN ! Parameter print
[3]	217	WRITE(numout,*)
	218	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
	219	WRITE(numout,*) ' ~~~~~~~'
[2528]	220	IF( ppkth == 0._wp ) THEN
[250]	221	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
	222	WRITE(numout,*) ' Total depth :', zhmax
	223	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
	224	ELSE
[2528]	225	IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
[250]	226	WRITE(numout,*) ' zsur, za0, za1 computed from '
	227	WRITE(numout,*) ' zdzmin = ', zdzmin
	228	WRITE(numout,*) ' zhmax = ', zhmax
	229	ENDIF
	230	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
	231	WRITE(numout,*) ' zsur = ', zsur
	232	WRITE(numout,*) ' za0 = ', za0
	233	WRITE(numout,*) ' za1 = ', za1
	234	WRITE(numout,*) ' zkth = ', zkth
	235	WRITE(numout,*) ' zacr = ', zacr
[2528]	236	IF( ldbletanh ) THEN
	237	WRITE(numout,*) ' (Double tanh za2 = ', za2
	238	WRITE(numout,*) ' parameters) zkth2= ', zkth2
	239	WRITE(numout,*) ' zacr2= ', zacr2
	240	ENDIF
[3]	241	ENDIF
	242	ENDIF
	243
	244
	245	! Reference z-coordinate (depth - scale factor at T- and W-points)
	246	! ======================
[2528]	247	IF( ppkth == 0._wp ) THEN ! uniform vertical grid
[454]	248	za1 = zhmax / FLOAT(jpk-1)
[250]	249	DO jk = 1, jpk
	250	zw = FLOAT( jk )
[2528]	251	zt = FLOAT( jk ) + 0.5_wp
[454]	252	gdepw_0(jk) = ( zw - 1 ) * za1
	253	gdept_0(jk) = ( zt - 1 ) * za1
	254	e3w_0 (jk) = za1
	255	e3t_0 (jk) = za1
[250]	256	END DO
[1099]	257	ELSE ! Madec & Imbard 1996 function
[2528]	258	IF( .NOT. ldbletanh ) THEN
	259	DO jk = 1, jpk
	260	zw = REAL( jk , wp )
	261	zt = REAL( jk , wp ) + 0.5_wp
	262	gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
	263	gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
	264	e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
	265	e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
	266	END DO
	267	ELSE
	268	DO jk = 1, jpk
	269	zw = FLOAT( jk )
	270	zt = FLOAT( jk ) + 0.5_wp
	271	! Double tanh function
	272	gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) &
	273	& + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) )
	274	gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) &
	275	& + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) )
	276	e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) &
	277	& + za2 * TANH( (zw-zkth2) / zacr2 )
	278	e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) &
	279	& + za2 * TANH( (zt-zkth2) / zacr2 )
	280	END DO
	281	ENDIF
	282	gdepw_0(1) = 0._wp ! force first w-level to be exactly at zero
[250]	283	ENDIF
	284
[1601]	285	!!gm BUG in s-coordinate this does not work!
[2528]	286	! deepest/shallowest W level Above/Below ~10m
	287	zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_0 ) ! ref. depth with tolerance (10% of minimum layer thickness)
	288	nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
	289	nla10 = nlb10 - 1 ! deepest W level Above ~10m
[1601]	290	!!gm end bug
[1577]	291
[1099]	292	IF(lwp) THEN ! control print
[3]	293	WRITE(numout,*)
	294	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
	295	WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" )
[454]	296	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
[3]	297	ENDIF
[1099]	298	DO jk = 1, jpk ! control positivity
[2528]	299	IF( e3w_0 (jk) <= 0._wp .OR. e3t_0 (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w or e3t =< 0 ' )
	300	IF( gdepw_0(jk) < 0._wp .OR. gdept_0(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw or gdept < 0 ' )
[3]	301	END DO
[1099]	302	!
[3294]	303	IF( nn_timing == 1 ) CALL timing_stop('zgr_z')
	304	!
[3]	305	END SUBROUTINE zgr_z
	306
	307
	308	SUBROUTINE zgr_bat
	309	!!----------------------------------------------------------------------
	310	!! * ROUTINE zgr_bat *
	311	!!
	312	!! ** Purpose : set bathymetry both in levels and meters
	313	!!
	314	!! ** Method : read or define mbathy and bathy arrays
	315	!! * level bathymetry:
	316	!! The ocean basin geometry is given by a two-dimensional array,
	317	!! mbathy, which is defined as follow :
	318	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
	319	!! at t-point (ji,jj).
	320	!! = 0 over the continental t-point.
	321	!! The array mbathy is checked to verified its consistency with
	322	!! model option. in particular:
	323	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	324	!! along closed boundary.
	325	!! mbathy must be cyclic IF jperio=1.
	326	!! mbathy must be lower or equal to jpk-1.
	327	!! isolated ocean grid points are suppressed from mbathy
	328	!! since they are only connected to remaining
	329	!! ocean through vertical diffusion.
	330	!! ntopo=-1 : rectangular channel or bassin with a bump
	331	!! ntopo= 0 : flat rectangular channel or basin
[128]	332	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
[3]	333	!! bathy is read in 'bathy_meter.nc' NetCDF file
	334	!!
	335	!! ** Action : - mbathy: level bathymetry (in level index)
	336	!! - bathy : meter bathymetry (in meters)
	337	!!----------------------------------------------------------------------
[1099]	338	INTEGER :: ji, jj, jl, jk ! dummy loop indices
	339	INTEGER :: inum ! temporary logical unit
[1348]	340	INTEGER :: ii_bump, ij_bump, ih ! bump center position
[2528]	341	INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices
[1099]	342	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
[2528]	343	REAL(wp) :: zi, zj, zh, zhmin ! local scalars
[3294]	344	INTEGER , POINTER, DIMENSION(:,:) :: idta ! global domain integer data
	345	REAL(wp), POINTER, DIMENSION(:,:) :: zdta ! global domain scalar data
[3]	346	!!----------------------------------------------------------------------
[3294]	347	!
	348	IF( nn_timing == 1 ) CALL timing_start('zgr_bat')
	349	!
	350	CALL wrk_alloc( jpidta, jpjdta, idta )
	351	CALL wrk_alloc( jpidta, jpjdta, zdta )
	352	!
[3]	353	IF(lwp) WRITE(numout,*)
	354	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
	355	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
	356
[1099]	357	! ! ================== !
	358	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
	359	! ! ================== !
	360	! ! global domain level and meter bathymetry (idta,zdta)
	361	!
[3]	362	IF( ntopo == 0 ) THEN ! flat basin
	363	IF(lwp) WRITE(numout,*)
	364	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
[1099]	365	idta(:,:) = jpkm1 ! before last level
	366	zdta(:,:) = gdepw_0(jpk) ! last w-point depth
[592]	367	h_oce = gdepw_0(jpk)
[1099]	368	ELSE ! bump centered in the basin
[3]	369	IF(lwp) WRITE(numout,*)
	370	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
[1099]	371	ii_bump = jpidta / 2 ! i-index of the bump center
	372	ij_bump = jpjdta / 2 ! j-index of the bump center
[2528]	373	r_bump = 50000._wp ! bump radius (meters)
	374	h_bump = 2700._wp ! bump height (meters)
[1099]	375	h_oce = gdepw_0(jpk) ! background ocean depth (meters)
[3]	376	IF(lwp) WRITE(numout,*) ' bump characteristics: '
	377	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
	378	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
	379	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
	380	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
[1099]	381	!
	382	DO jj = 1, jpjdta ! zdta :
[3]	383	DO ji = 1, jpidta
[592]	384	zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
	385	zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
[3]	386	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
	387	END DO
	388	END DO
[1099]	389	! ! idta :
	390	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
[454]	391	idta(:,:) = jpkm1
[1099]	392	ELSE ! z-coordinate (zco or zps): step-like topography
[454]	393	idta(:,:) = jpkm1
	394	DO jk = 1, jpkm1
[2528]	395	WHERE( gdept_0(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_0(jk+1) ) idta(:,:) = jk
[3]	396	END DO
[454]	397	ENDIF
[3]	398	ENDIF
[1099]	399	! ! set GLOBAL boundary conditions
	400	! ! Caution : idta on the global domain: use of jperio, not nperio
[3]	401	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
[2528]	402	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp
	403	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp
[3]	404	ELSEIF( jperio == 2 ) THEN
[30]	405	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
[2528]	406	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp
	407	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp
	408	idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0._wp
[3]	409	ELSE
[2528]	410	ih = 0 ; zh = 0._wp
[525]	411	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
[454]	412	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
	413	idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh
	414	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
	415	idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh
[3]	416	ENDIF
	417
[1099]	418	! ! local domain level and meter bathymetries (mbathy,bathy)
	419	mbathy(:,:) = 0 ! set to zero extra halo points
[2528]	420	bathy (:,:) = 0._wp ! (require for mpp case)
[1099]	421	DO jj = 1, nlcj ! interior values
[473]	422	DO ji = 1, nlci
	423	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
	424	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
	425	END DO
	426	END DO
[1099]	427	!
	428	! ! ================ !
	429	ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain)
	430	! ! ================ !
	431	!
	432	IF( ln_zco ) THEN ! zco : read level bathymetry
[2528]	433	CALL iom_open ( 'bathy_level.nc', inum )
	434	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
	435	CALL iom_close( inum )
[473]	436	mbathy(:,:) = INT( bathy(:,:) )
[3421]	437	!
[1273]	438	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
[3421]	439	!
[2528]	440	IF( nn_cla == 0 ) THEN
[1273]	441	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	442	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
	443	DO ji = mi0(ii0), mi1(ii1)
	444	DO jj = mj0(ij0), mj1(ij1)
	445	mbathy(ji,jj) = 15
	446	END DO
	447	END DO
	448	IF(lwp) WRITE(numout,*)
[2528]	449	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
[1273]	450	!
	451	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	452	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
	453	DO ji = mi0(ii0), mi1(ii1)
	454	DO jj = mj0(ij0), mj1(ij1)
	455	mbathy(ji,jj) = 12
	456	END DO
	457	END DO
	458	IF(lwp) WRITE(numout,*)
[2528]	459	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
[1273]	460	ENDIF
	461	!
	462	ENDIF
	463	!
[454]	464	ENDIF
[1099]	465	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
[2528]	466	CALL iom_open ( 'bathy_meter.nc', inum )
	467	CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
	468	CALL iom_close( inum )
[3421]	469	!
[2528]	470	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
[3421]	471	!
[2528]	472	IF( nn_cla == 0 ) THEN
	473	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	474	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
[1273]	475	DO ji = mi0(ii0), mi1(ii1)
	476	DO jj = mj0(ij0), mj1(ij1)
[2528]	477	bathy(ji,jj) = 284._wp
[1273]	478	END DO
	479	END DO
	480	IF(lwp) WRITE(numout,*)
[2528]	481	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
[1273]	482	!
[2528]	483	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	484	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
[1273]	485	DO ji = mi0(ii0), mi1(ii1)
	486	DO jj = mj0(ij0), mj1(ij1)
[2528]	487	bathy(ji,jj) = 137._wp
[1273]	488	END DO
	489	END DO
	490	IF(lwp) WRITE(numout,*)
	491	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
	492	ENDIF
	493	!
	494	ENDIF
[1348]	495	!
	496	ENDIF
[3]	497	! ! =============== !
	498	ELSE ! error !
	499	! ! =============== !
[1099]	500	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
[473]	501	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
[3]	502	ENDIF
[1099]	503	!
[3421]	504	IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==!
	505	!
	506	IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
[2712]	507	IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
	508	ELSE ; ik = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth
	509	ENDIF
	510	zhmin = gdepw_0(ik+1) ! minimum depth = ik+1 w-levels
	511	WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
	512	ELSE WHERE ; bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
	513	END WHERE
	514	IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
[2528]	515	ENDIF
	516	!
[3294]	517	CALL wrk_dealloc( jpidta, jpjdta, idta )
	518	CALL wrk_dealloc( jpidta, jpjdta, zdta )
	519	!
	520	IF( nn_timing == 1 ) CALL timing_stop('zgr_bat')
	521	!
[3]	522	END SUBROUTINE zgr_bat
	523
	524
	525	SUBROUTINE zgr_bat_zoom
	526	!!----------------------------------------------------------------------
	527	!! * ROUTINE zgr_bat_zoom *
	528	!!
	529	!! ** Purpose : - Close zoom domain boundary if necessary
	530	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
	531	!!
	532	!! ** Method :
	533	!!
	534	!! ** Action : - update mbathy: level bathymetry (in level index)
	535	!!----------------------------------------------------------------------
	536	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
	537	!!----------------------------------------------------------------------
[1099]	538	!
[3]	539	IF(lwp) WRITE(numout,*)
	540	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
	541	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
[1099]	542	!
[3]	543	! Zoom domain
	544	! ===========
[1099]	545	!
[3]	546	! Forced closed boundary if required
[1099]	547	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
	548	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
	549	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
	550	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
	551	!
[3]	552	! Configuration specific domain modifications
	553	! (here, ORCA arctic configuration: suppress Med Sea)
	554	IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
	555	SELECT CASE ( jp_cfg )
	556	! ! =======================
	557	CASE ( 2 ) ! ORCA_R2 configuration
	558	! ! =======================
	559	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
	560	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
	561	ij0 = 98 ; ij1 = 110
	562	! ! =======================
	563	CASE ( 05 ) ! ORCA_R05 configuration
	564	! ! =======================
	565	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
	566	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
	567	ij0 = 314 ; ij1 = 370
	568	END SELECT
	569	!
	570	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
	571	!
	572	ENDIF
[1099]	573	!
[3]	574	END SUBROUTINE zgr_bat_zoom
	575
	576
	577	SUBROUTINE zgr_bat_ctl
	578	!!----------------------------------------------------------------------
	579	!! * ROUTINE zgr_bat_ctl *
	580	!!
	581	!! ** Purpose : check the bathymetry in levels
	582	!!
	583	!! ** Method : The array mbathy is checked to verified its consistency
	584	!! with the model options. in particular:
	585	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	586	!! along closed boundary.
	587	!! mbathy must be cyclic IF jperio=1.
	588	!! mbathy must be lower or equal to jpk-1.
	589	!! isolated ocean grid points are suppressed from mbathy
	590	!! since they are only connected to remaining
	591	!! ocean through vertical diffusion.
	592	!! C A U T I O N : mbathy will be modified during the initializa-
	593	!! tion phase to become the number of non-zero w-levels of a water
	594	!! column, with a minimum value of 1.
	595	!!
	596	!! ** Action : - update mbathy: level bathymetry (in level index)
	597	!! - update bathy : meter bathymetry (in meters)
	598	!!----------------------------------------------------------------------
[2715]	599	!!
[1099]	600	INTEGER :: ji, jj, jl ! dummy loop indices
	601	INTEGER :: icompt, ibtest, ikmax ! temporary integers
[3294]	602	REAL(wp), POINTER, DIMENSION(:,:) :: zbathy
[3]	603	!!----------------------------------------------------------------------
[3294]	604	!
	605	IF( nn_timing == 1 ) CALL timing_start('zgr_bat_ctl')
	606	!
	607	CALL wrk_alloc( jpi, jpj, zbathy )
	608	!
[3]	609	IF(lwp) WRITE(numout,*)
	610	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
	611	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
	612
[1099]	613	! ! Suppress isolated ocean grid points
	614	IF(lwp) WRITE(numout,*)
	615	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
	616	IF(lwp) WRITE(numout,*)' -----------------------------------'
	617	icompt = 0
	618	DO jl = 1, 2
	619	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
	620	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
	621	mbathy(jpi,:) = mbathy( 2 ,:)
	622	ENDIF
	623	DO jj = 2, jpjm1
	624	DO ji = 2, jpim1
	625	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
	626	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
	627	IF( ibtest < mbathy(ji,jj) ) THEN
	628	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
	629	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
	630	mbathy(ji,jj) = ibtest
	631	icompt = icompt + 1
	632	ENDIF
	633	END DO
	634	END DO
	635	END DO
	636	IF( icompt == 0 ) THEN
	637	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
	638	ELSE
	639	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
	640	ENDIF
	641	IF( lk_mpp ) THEN
	642	zbathy(:,:) = FLOAT( mbathy(:,:) )
[2528]	643	CALL lbc_lnk( zbathy, 'T', 1._wp )
[1099]	644	mbathy(:,:) = INT( zbathy(:,:) )
	645	ENDIF
[3]	646
[1099]	647	! ! East-west cyclic boundary conditions
	648	IF( nperio == 0 ) THEN
	649	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
	650	IF( lk_mpp ) THEN
	651	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
	652	IF( jperio /= 1 ) mbathy(1,:) = 0
[411]	653	ENDIF
[1099]	654	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
	655	IF( jperio /= 1 ) mbathy(nlci,:) = 0
	656	ENDIF
[411]	657	ELSE
[1099]	658	IF( ln_zco .OR. ln_zps ) THEN
	659	mbathy( 1 ,:) = 0
	660	mbathy(jpi,:) = 0
	661	ELSE
	662	mbathy( 1 ,:) = jpkm1
	663	mbathy(jpi,:) = jpkm1
	664	ENDIF
[411]	665	ENDIF
[1099]	666	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
	667	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
	668	mbathy( 1 ,:) = mbathy(jpim1,:)
	669	mbathy(jpi,:) = mbathy( 2 ,:)
	670	ELSEIF( nperio == 2 ) THEN
	671	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
	672	ELSE
	673	IF(lwp) WRITE(numout,*) ' e r r o r'
	674	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
	675	! STOP 'dom_mba'
	676	ENDIF
	677
[1528]	678	! Boundary condition on mbathy
	679	IF( .NOT.lk_mpp ) THEN
	680	!!gm !!bug ??? think about it !
	681	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
	682	zbathy(:,:) = FLOAT( mbathy(:,:) )
[2528]	683	CALL lbc_lnk( zbathy, 'T', 1._wp )
[1528]	684	mbathy(:,:) = INT( zbathy(:,:) )
[3]	685	ENDIF
	686
	687	! Number of ocean level inferior or equal to jpkm1
	688	ikmax = 0
	689	DO jj = 1, jpj
	690	DO ji = 1, jpi
	691	ikmax = MAX( ikmax, mbathy(ji,jj) )
	692	END DO
	693	END DO
[1099]	694	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
[3]	695	IF( ikmax > jpkm1 ) THEN
	696	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
	697	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
	698	ELSE IF( ikmax < jpkm1 ) THEN
	699	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
	700	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
	701	ENDIF
	702
[1566]	703	IF( lwp .AND. nprint == 1 ) THEN ! control print
[3]	704	WRITE(numout,*)
[1099]	705	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
[3]	706	WRITE(numout,*) ' ------------------'
[1099]	707	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
[3]	708	WRITE(numout,*)
	709	ENDIF
[1099]	710	!
[3294]	711	CALL wrk_dealloc( jpi, jpj, zbathy )
[2715]	712	!
[3294]	713	IF( nn_timing == 1 ) CALL timing_stop('zgr_bat_ctl')
	714	!
[3]	715	END SUBROUTINE zgr_bat_ctl
	716
	717
[2528]	718	SUBROUTINE zgr_bot_level
	719	!!----------------------------------------------------------------------
	720	!! * ROUTINE zgr_bot_level *
	721	!!
	722	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
	723	!!
	724	!! ** Method : computes from mbathy with a minimum value of 1 over land
	725	!!
	726	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
	727	!! ocean level at t-, u- & v-points
	728	!! (min value = 1 over land)
	729	!!----------------------------------------------------------------------
[2715]	730	!!
[2528]	731	INTEGER :: ji, jj ! dummy loop indices
[3294]	732	REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
[2528]	733	!!----------------------------------------------------------------------
	734	!
[3294]	735	IF( nn_timing == 1 ) CALL timing_start('zgr_bot_level')
[2715]	736	!
[3294]	737	CALL wrk_alloc( jpi, jpj, zmbk )
	738	!
[2528]	739	IF(lwp) WRITE(numout,*)
	740	IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
	741	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
	742	!
	743	mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
	744	! ! bottom k-index of W-level = mbkt+1
	745	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
	746	DO ji = 1, jpim1
	747	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
	748	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
	749	END DO
	750	END DO
	751	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
	752	zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
	753	zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
	754	!
[3294]	755	CALL wrk_dealloc( jpi, jpj, zmbk )
[2715]	756	!
[3294]	757	IF( nn_timing == 1 ) CALL timing_stop('zgr_bot_level')
	758	!
[2528]	759	END SUBROUTINE zgr_bot_level
	760
	761
[454]	762	SUBROUTINE zgr_zco
	763	!!----------------------------------------------------------------------
	764	!! * ROUTINE zgr_zco *
	765	!!
	766	!! ** Purpose : define the z-coordinate system
	767	!!
[2528]	768	!! ** Method : set 3D coord. arrays to reference 1D array
[454]	769	!!----------------------------------------------------------------------
	770	INTEGER :: jk
	771	!!----------------------------------------------------------------------
[1099]	772	!
[3294]	773	IF( nn_timing == 1 ) CALL timing_start('zgr_zco')
	774	!
[2528]	775	DO jk = 1, jpk
[3294]	776	gdept(:,:,jk) = gdept_0(jk)
	777	gdepw(:,:,jk) = gdepw_0(jk)
	778	gdep3w(:,:,jk) = gdepw_0(jk)
	779	e3t (:,:,jk) = e3t_0(jk)
	780	e3u (:,:,jk) = e3t_0(jk)
	781	e3v (:,:,jk) = e3t_0(jk)
	782	e3f (:,:,jk) = e3t_0(jk)
	783	e3w (:,:,jk) = e3w_0(jk)
	784	e3uw(:,:,jk) = e3w_0(jk)
	785	e3vw(:,:,jk) = e3w_0(jk)
[2528]	786	END DO
[1099]	787	!
[3294]	788	IF( nn_timing == 1 ) CALL timing_stop('zgr_zco')
	789	!
[454]	790	END SUBROUTINE zgr_zco
	791
	792
[1083]	793	SUBROUTINE zgr_zps
	794	!!----------------------------------------------------------------------
	795	!! * ROUTINE zgr_zps *
	796	!!
	797	!! ** Purpose : the depth and vertical scale factor in partial step
	798	!! z-coordinate case
	799	!!
	800	!! ** Method : Partial steps : computes the 3D vertical scale factors
	801	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
	802	!! a partial step representation of bottom topography.
	803	!!
	804	!! The reference depth of model levels is defined from an analytical
	805	!! function the derivative of which gives the reference vertical
	806	!! scale factors.
	807	!! From depth and scale factors reference, we compute there new value
	808	!! with partial steps on 3d arrays ( i, j, k ).
	809	!!
	810	!! w-level: gdepw(i,j,k) = fsdep(k)
	811	!! e3w(i,j,k) = dk(fsdep)(k) = fse3(i,j,k)
	812	!! t-level: gdept(i,j,k) = fsdep(k+0.5)
	813	!! e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
	814	!!
	815	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
	816	!! we find the mbathy index of the depth at each grid point.
	817	!! This leads us to three cases:
	818	!!
	819	!! - bathy = 0 => mbathy = 0
	820	!! - 1 < mbathy < jpkm1
	821	!! - bathy > gdepw(jpk) => mbathy = jpkm1
	822	!!
	823	!! Then, for each case, we find the new depth at t- and w- levels
	824	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
	825	!! and f-points.
	826	!!
	827	!! This routine is given as an example, it must be modified
	828	!! following the user s desiderata. nevertheless, the output as
	829	!! well as the way to compute the model levels and scale factors
	830	!! must be respected in order to insure second order accuracy
	831	!! schemes.
	832	!!
	833	!! c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
	834	!! - - - - - - - gdept, gdepw and e3. are positives
	835	!!
[1099]	836	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
[1083]	837	!!----------------------------------------------------------------------
[2715]	838	!!
[1099]	839	INTEGER :: ji, jj, jk ! dummy loop indices
	840	INTEGER :: ik, it ! temporary integers
	841	LOGICAL :: ll_print ! Allow control print for debugging
	842	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
	843	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
[2528]	844	REAL(wp) :: zmax ! Maximum depth
[1099]	845	REAL(wp) :: zdiff ! temporary scalar
[2528]	846	REAL(wp) :: zrefdep ! temporary scalar
[3294]	847	REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt
[1099]	848	!!---------------------------------------------------------------------
[3294]	849	!
	850	IF( nn_timing == 1 ) CALL timing_start('zgr_zps')
	851	!
	852	CALL wrk_alloc( jpi, jpj, jpk, zprt )
	853	!
[1099]	854	IF(lwp) WRITE(numout,*)
	855	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
	856	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
	857	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
[3]	858
[2528]	859	ll_print = .FALSE. ! Local variable for debugging
[1083]	860
[1099]	861	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
[1083]	862	WRITE(numout,*)
	863	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
	864	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
	865	ENDIF
	866
	867
	868	! bathymetry in level (from bathy_meter)
	869	! ===================
[2528]	870	zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
	871	bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
	872	WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
	873	ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
	874	END WHERE
[1083]	875
	876	! Compute mbathy for ocean points (i.e. the number of ocean levels)
	877	! find the number of ocean levels such that the last level thickness
	878	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
	879	! e3t_0 is the reference level thickness
	880	DO jk = jpkm1, 1, -1
	881	zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
[2528]	882	WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
[1083]	883	END DO
	884
[1099]	885	! Scale factors and depth at T- and W-points
	886	DO jk = 1, jpk ! intitialization to the reference z-coordinate
	887	gdept(:,:,jk) = gdept_0(jk)
	888	gdepw(:,:,jk) = gdepw_0(jk)
	889	e3t (:,:,jk) = e3t_0 (jk)
	890	e3w (:,:,jk) = e3w_0 (jk)
[1083]	891	END DO
[1099]	892	!
	893	DO jj = 1, jpj
	894	DO ji = 1, jpi
	895	ik = mbathy(ji,jj)
	896	IF( ik > 0 ) THEN ! ocean point only
	897	! max ocean level case
	898	IF( ik == jpkm1 ) THEN
	899	zdepwp = bathy(ji,jj)
	900	ze3tp = bathy(ji,jj) - gdepw_0(ik)
[2528]	901	ze3wp = 0.5_wp * e3w_0(ik) * ( 1._wp + ( ze3tp/e3t_0(ik) ) )
[1099]	902	e3t(ji,jj,ik ) = ze3tp
	903	e3t(ji,jj,ik+1) = ze3tp
	904	e3w(ji,jj,ik ) = ze3wp
	905	e3w(ji,jj,ik+1) = ze3tp
	906	gdepw(ji,jj,ik+1) = zdepwp
	907	gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp
	908	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
	909	!
	910	ELSE ! standard case
	911	IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj)
	912	ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
	913	ENDIF
	914	!gm Bug? check the gdepw_0
	915	! ... on ik
	916	gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
	917	& * ((gdept_0( ik ) - gdepw_0(ik) ) &
	918	& / ( gdepw_0( ik+1) - gdepw_0(ik) ))
	919	e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
	920	& / ( gdepw_0( ik+1) - gdepw_0(ik) )
[2528]	921	e3w (ji,jj,ik) = 0.5_wp * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2._wp * gdepw_0(ik) ) &
	922	& * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
[1099]	923	! ... on ik+1
	924	e3w (ji,jj,ik+1) = e3t (ji,jj,ik)
	925	e3t (ji,jj,ik+1) = e3t (ji,jj,ik)
	926	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
	927	ENDIF
	928	ENDIF
	929	END DO
	930	END DO
	931	!
	932	it = 0
	933	DO jj = 1, jpj
	934	DO ji = 1, jpi
	935	ik = mbathy(ji,jj)
	936	IF( ik > 0 ) THEN ! ocean point only
	937	e3tp (ji,jj) = e3t(ji,jj,ik )
	938	e3wp (ji,jj) = e3w(ji,jj,ik )
	939	! test
	940	zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik )
[2528]	941	IF( zdiff <= 0._wp .AND. lwp ) THEN
[1099]	942	it = it + 1
	943	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
	944	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
	945	WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
	946	WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik )
	947	ENDIF
	948	ENDIF
	949	END DO
	950	END DO
[1083]	951
	952	! Scale factors and depth at U-, V-, UW and VW-points
[1099]	953	DO jk = 1, jpk ! initialisation to z-scale factors
[2528]	954	e3u (:,:,jk) = e3t_0(jk)
	955	e3v (:,:,jk) = e3t_0(jk)
	956	e3uw(:,:,jk) = e3w_0(jk)
	957	e3vw(:,:,jk) = e3w_0(jk)
[1083]	958	END DO
[1099]	959	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
	960	DO jj = 1, jpjm1
	961	DO ji = 1, fs_jpim1 ! vector opt.
[2528]	962	e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk) )
	963	e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk) )
[1099]	964	e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
	965	e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
	966	END DO
	967	END DO
	968	END DO
[2528]	969	CALL lbc_lnk( e3u , 'U', 1._wp ) ; CALL lbc_lnk( e3uw, 'U', 1._wp ) ! lateral boundary conditions
	970	CALL lbc_lnk( e3v , 'V', 1._wp ) ; CALL lbc_lnk( e3vw, 'V', 1._wp )
[1099]	971	!
	972	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
[2528]	973	WHERE( e3u (:,:,jk) == 0._wp ) e3u (:,:,jk) = e3t_0(jk)
	974	WHERE( e3v (:,:,jk) == 0._wp ) e3v (:,:,jk) = e3t_0(jk)
	975	WHERE( e3uw(:,:,jk) == 0._wp ) e3uw(:,:,jk) = e3w_0(jk)
	976	WHERE( e3vw(:,:,jk) == 0._wp ) e3vw(:,:,jk) = e3w_0(jk)
[1099]	977	END DO
	978
	979	! Scale factor at F-point
	980	DO jk = 1, jpk ! initialisation to z-scale factors
[2528]	981	e3f(:,:,jk) = e3t_0(jk)
[1099]	982	END DO
	983	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
	984	DO jj = 1, jpjm1
	985	DO ji = 1, fs_jpim1 ! vector opt.
	986	e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
	987	END DO
	988	END DO
	989	END DO
[2528]	990	CALL lbc_lnk( e3f, 'F', 1._wp ) ! Lateral boundary conditions
[1099]	991	!
	992	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
[2528]	993	WHERE( e3f(:,:,jk) == 0._wp ) e3f(:,:,jk) = e3t_0(jk)
[1099]	994	END DO
	995	!!gm bug ? : must be a do loop with mj0,mj1
	996	!
	997	e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
	998	e3w(:,mj0(1),:) = e3w(:,mj0(2),:)
	999	e3u(:,mj0(1),:) = e3u(:,mj0(2),:)
	1000	e3v(:,mj0(1),:) = e3v(:,mj0(2),:)
	1001	e3f(:,mj0(1),:) = e3f(:,mj0(2),:)
[1083]	1002
[1099]	1003	! Control of the sign
[2528]	1004	IF( MINVAL( e3t (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' )
	1005	IF( MINVAL( e3w (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' )
	1006	IF( MINVAL( gdept(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
	1007	IF( MINVAL( gdepw(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
[1083]	1008
[1099]	1009	! Compute gdep3w (vertical sum of e3w)
[2528]	1010	gdep3w(:,:,1) = 0.5_wp * e3w(:,:,1)
[1099]	1011	DO jk = 2, jpk
	1012	gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk)
	1013	END DO
[2528]	1014
[1099]	1015	! ! ================= !
	1016	IF(lwp .AND. ll_print) THEN ! Control print !
	1017	! ! ================= !
[1083]	1018	DO jj = 1,jpj
	1019	DO ji = 1, jpi
[1099]	1020	ik = MAX( mbathy(ji,jj), 1 )
	1021	zprt(ji,jj,1) = e3t (ji,jj,ik)
	1022	zprt(ji,jj,2) = e3w (ji,jj,ik)
	1023	zprt(ji,jj,3) = e3u (ji,jj,ik)
	1024	zprt(ji,jj,4) = e3v (ji,jj,ik)
	1025	zprt(ji,jj,5) = e3f (ji,jj,ik)
	1026	zprt(ji,jj,6) = gdep3w(ji,jj,ik)
[1083]	1027	END DO
	1028	END DO
	1029	WRITE(numout,*)
[1099]	1030	WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1031	WRITE(numout,*)
[3294]	1032	WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1033	WRITE(numout,*)
[3294]	1034	WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1035	WRITE(numout,*)
[3294]	1036	WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1037	WRITE(numout,*)
[3294]	1038	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1039	WRITE(numout,*)
[3294]	1040	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1041	ENDIF
[2528]	1042	!
[3294]	1043	CALL wrk_dealloc( jpi, jpj, jpk, zprt )
[2715]	1044	!
[3294]	1045	IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')
	1046	!
[1083]	1047	END SUBROUTINE zgr_zps
	1048
[454]	1049	SUBROUTINE zgr_sco
	1050	!!----------------------------------------------------------------------
	1051	!! * ROUTINE zgr_sco *
	1052	!!
	1053	!! ** Purpose : define the s-coordinate system
	1054	!!
	1055	!! ** Method : s-coordinate
	1056	!! The depth of model levels is defined as the product of an
	1057	!! analytical function by the local bathymetry, while the vertical
	1058	!! scale factors are defined as the product of the first derivative
	1059	!! of the analytical function by the bathymetry.
	1060	!! (this solution save memory as depth and scale factors are not
	1061	!! 3d fields)
	1062	!! - Read bathymetry (in meters) at t-point and compute the
	1063	!! bathymetry at u-, v-, and f-points.
	1064	!! hbatu = mi( hbatt )
	1065	!! hbatv = mj( hbatt )
	1066	!! hbatf = mi( mj( hbatt ) )
[3618]	1067	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
[1083]	1068	!! function and its derivative given as function.
[3618]	1069	!! z_gsigt(k) = fssig (k )
	1070	!! z_gsigw(k) = fssig (k-0.5)
	1071	!! z_esigt(k) = fsdsig(k )
	1072	!! z_esigw(k) = fsdsig(k-0.5)
[3600]	1073	!! Three options for stretching are give, and they can be modified
	1074	!! following the users requirements. Nevertheless, the output as
[454]	1075	!! well as the way to compute the model levels and scale factors
[3600]	1076	!! must be respected in order to insure second order accuracy
[454]	1077	!! schemes.
	1078	!!
[3600]	1079	!! The three methods for stretching available are:
	1080	!!
	1081	!! s_sh94 (Song and Haidvogel 1994)
	1082	!! a sinh/tanh function that allows sigma and stretched sigma
	1083	!!
	1084	!! s_sf12 (Siddorn and Furner 2012?)
	1085	!! allows the maintenance of fixed surface and or
	1086	!! bottom cell resolutions (cf. geopotential coordinates)
	1087	!! within an analytically derived stretched S-coordinate framework.
	1088	!!
	1089	!! s_tanh (Madec et al 1996)
	1090	!! a cosh/tanh function that gives stretched coordinates
	1091	!!
[1099]	1092	!!----------------------------------------------------------------------
[2715]	1093	!
[1099]	1094	INTEGER :: ji, jj, jk, jl ! dummy loop argument
	1095	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
[3600]	1096	REAL(wp) :: zrmax, ztaper ! temporary scalars
[2715]	1097	!
[3294]	1098	REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
[2715]	1099
[3600]	1100	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
	1101	rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
	1102	!!----------------------------------------------------------------------
[3294]	1103	!
	1104	IF( nn_timing == 1 ) CALL timing_start('zgr_sco')
	1105	!
	1106	CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat )
	1107	!
[2715]	1108	REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters
[1601]	1109	READ ( numnam, namzgr_sco )
[454]	1110
[2715]	1111	IF(lwp) THEN ! control print
[454]	1112	WRITE(numout,*)
	1113	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
	1114	WRITE(numout,*) '~~~~~~~~~~~'
[1601]	1115	WRITE(numout,*) ' Namelist namzgr_sco'
[3600]	1116	WRITE(numout,*) ' stretching coeffs '
	1117	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
	1118	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
	1119	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
	1120	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
	1121	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
	1122	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
	1123	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
	1124	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
	1125	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
	1126	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
[3618]	1127	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
[3600]	1128	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
	1129	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
	1130	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
	1131	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
	1132	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
	1133	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
	1134	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
[454]	1135	ENDIF
	1136
[1601]	1137	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
	1138	hifu(:,:) = rn_sbot_min
	1139	hifv(:,:) = rn_sbot_min
	1140	hiff(:,:) = rn_sbot_min
[1348]	1141
	1142	! ! set maximum ocean depth
[1601]	1143	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
[454]	1144
[1461]	1145	DO jj = 1, jpj
	1146	DO ji = 1, jpi
[2715]	1147	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
[1461]	1148	END DO
	1149	END DO
[1099]	1150	! ! =============================
	1151	! ! Define the envelop bathymetry (hbatt)
	1152	! ! =============================
[454]	1153	! use r-value to create hybrid coordinates
	1154	DO jj = 1, jpj
	1155	DO ji = 1, jpi
[1601]	1156	zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
[454]	1157	END DO
	1158	END DO
[1639]	1159	!
	1160	! Smooth the bathymetry (if required)
[2528]	1161	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
[1639]	1162	scobot(:,:) = bathy(:,:) ! ocean bottom depth
	1163	!
[454]	1164	jl = 0
[2528]	1165	zrmax = 1._wp
[454]	1166	! ! ================ !
[2528]	1167	DO WHILE( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop !
[454]	1168	! ! ================ !
	1169	jl = jl + 1
[2528]	1170	zrmax = 0._wp
	1171	zmsk(:,:) = 0._wp
[454]	1172	DO jj = 1, nlcj
	1173	DO ji = 1, nlci
	1174	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
	1175	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
	1176	zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
	1177	zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
	1178	zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
[2528]	1179	IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp
	1180	IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1._wp
	1181	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1._wp
	1182	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1._wp
[454]	1183	END DO
	1184	END DO
	1185	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
	1186	! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
[2528]	1187	ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1._wp )
[454]	1188	DO jj = 1, nlcj
	1189	DO ji = 1, nlci
[1348]	1190	zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
[454]	1191	END DO
	1192	END DO
[1348]	1193	!
[1099]	1194	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
	1195	!
[454]	1196	DO jj = 1, nlcj
	1197	DO ji = 1, nlci
	1198	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
	1199	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
	1200	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
	1201	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
[2528]	1202	IF( zmsk(ji,jj) == 1._wp ) THEN
[454]	1203	ztmp(ji,jj) = ( &
	1204	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
[2528]	1205	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2._wp + zenv(iip1,jj )*zmsk(iip1,jj ) &
[454]	1206	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
	1207	& ) / ( &
	1208	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
[2528]	1209	& + zmsk(iim1,jj ) + 2._wp + zmsk(iip1,jj ) &
[454]	1210	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
	1211	& )
	1212	ENDIF
	1213	END DO
	1214	END DO
[1348]	1215	!
[454]	1216	DO jj = 1, nlcj
	1217	DO ji = 1, nlci
[2528]	1218	IF( zmsk(ji,jj) == 1._wp ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
[454]	1219	END DO
	1220	END DO
[1099]	1221	!
[454]	1222	! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero
[2528]	1223	ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1._wp )
[454]	1224	DO jj = 1, nlcj
	1225	DO ji = 1, nlci
[2528]	1226	IF( zenv(ji,jj) == 0._wp ) zenv(ji,jj) = ztmp(ji,jj)
[454]	1227	END DO
	1228	END DO
	1229	! ! ================ !
	1230	END DO ! End loop !
	1231	! ! ================ !
[1099]	1232	!
	1233	! ! envelop bathymetry saved in hbatt
[454]	1234	hbatt(:,:) = zenv(:,:)
[2528]	1235	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
[1099]	1236	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
	1237	DO jj = 1, jpj
	1238	DO ji = 1, jpi
[2528]	1239	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2 )
	1240	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
[1099]	1241	END DO
	1242	END DO
[516]	1243	ENDIF
[1099]	1244	!
	1245	IF(lwp) THEN ! Control print
[454]	1246	WRITE(numout,*)
	1247	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
	1248	WRITE(numout,*)
[2528]	1249	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
[1099]	1250	IF( nprint == 1 ) THEN
	1251	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
	1252	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
	1253	ENDIF
[454]	1254	ENDIF
	1255
[1099]	1256	! ! ==============================
	1257	! ! hbatu, hbatv, hbatf fields
	1258	! ! ==============================
[454]	1259	IF(lwp) THEN
	1260	WRITE(numout,*)
[1601]	1261	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
[454]	1262	ENDIF
[1601]	1263	hbatu(:,:) = rn_sbot_min
	1264	hbatv(:,:) = rn_sbot_min
	1265	hbatf(:,:) = rn_sbot_min
[454]	1266	DO jj = 1, jpjm1
[1694]	1267	DO ji = 1, jpim1 ! NO vector opt.
[2528]	1268	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
	1269	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
	1270	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
	1271	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
[454]	1272	END DO
	1273	END DO
[1099]	1274	!
[454]	1275	! Apply lateral boundary condition
[1099]	1276	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
[2528]	1277	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1._wp )
[454]	1278	DO jj = 1, jpj
	1279	DO ji = 1, jpi
[2528]	1280	IF( hbatu(ji,jj) == 0._wp ) THEN
	1281	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
	1282	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
[454]	1283	ENDIF
	1284	END DO
	1285	END DO
[2528]	1286	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1._wp )
[454]	1287	DO jj = 1, jpj
	1288	DO ji = 1, jpi
[2528]	1289	IF( hbatv(ji,jj) == 0._wp ) THEN
	1290	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
	1291	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
[454]	1292	ENDIF
	1293	END DO
	1294	END DO
[2528]	1295	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1._wp )
[454]	1296	DO jj = 1, jpj
	1297	DO ji = 1, jpi
[2528]	1298	IF( hbatf(ji,jj) == 0._wp ) THEN
	1299	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
	1300	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
[454]	1301	ENDIF
	1302	END DO
	1303	END DO
	1304
	1305	!!bug: key_helsinki a verifer
	1306	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
	1307	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
	1308	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
	1309	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
	1310
[516]	1311	IF( nprint == 1 .AND. lwp ) THEN
[1099]	1312	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
	1313	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
	1314	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
	1315	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
[516]	1316	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
	1317	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
	1318	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
	1319	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
	1320	ENDIF
[454]	1321	!! helsinki
	1322
[1099]	1323	! ! =======================
	1324	! ! s-ccordinate fields (gdep., e3.)
	1325	! ! =======================
	1326	!
	1327	! non-dimensional "sigma" for model level depth at w- and t-levels
[1348]	1328
	1329
[3600]	1330	!========================================================================
	1331	! Song and Haidvogel 1994 (ln_s_sh94=T)
	1332	! Siddorn and Furner 2012 (ln_sf12=T)
	1333	! or tanh function (both false)
	1334	!========================================================================
	1335	IF ( ln_s_sh94 ) THEN
	1336	CALL s_sh94()
	1337	ELSE IF ( ln_s_sf12 ) THEN
	1338	CALL s_sf12()
	1339	ELSE
	1340	CALL s_tanh()
	1341	ENDIF
[2528]	1342
[3600]	1343	CALL lbc_lnk( e3t , 'T', 1._wp )
	1344	CALL lbc_lnk( e3u , 'U', 1._wp )
	1345	CALL lbc_lnk( e3v , 'V', 1._wp )
	1346	CALL lbc_lnk( e3f , 'F', 1._wp )
	1347	CALL lbc_lnk( e3w , 'W', 1._wp )
	1348	CALL lbc_lnk( e3uw, 'U', 1._wp )
	1349	CALL lbc_lnk( e3vw, 'V', 1._wp )
[2715]	1350
[3600]	1351	fsdepw(:,:,:) = gdepw (:,:,:)
	1352	fsde3w(:,:,:) = gdep3w(:,:,:)
[1099]	1353	!
[3294]	1354	where (e3t (:,:,:).eq.0.0) e3t(:,:,:) = 1.0
	1355	where (e3u (:,:,:).eq.0.0) e3u(:,:,:) = 1.0
	1356	where (e3v (:,:,:).eq.0.0) e3v(:,:,:) = 1.0
	1357	where (e3f (:,:,:).eq.0.0) e3f(:,:,:) = 1.0
	1358	where (e3w (:,:,:).eq.0.0) e3w(:,:,:) = 1.0
	1359	where (e3uw (:,:,:).eq.0.0) e3uw(:,:,:) = 1.0
	1360	where (e3vw (:,:,:).eq.0.0) e3vw(:,:,:) = 1.0
[1461]	1361
[3294]	1362
[1461]	1363	fsdept(:,:,:) = gdept (:,:,:)
	1364	fsdepw(:,:,:) = gdepw (:,:,:)
	1365	fsde3w(:,:,:) = gdep3w(:,:,:)
	1366	fse3t (:,:,:) = e3t (:,:,:)
	1367	fse3u (:,:,:) = e3u (:,:,:)
	1368	fse3v (:,:,:) = e3v (:,:,:)
	1369	fse3f (:,:,:) = e3f (:,:,:)
	1370	fse3w (:,:,:) = e3w (:,:,:)
	1371	fse3uw(:,:,:) = e3uw (:,:,:)
	1372	fse3vw(:,:,:) = e3vw (:,:,:)
	1373	!!
[1099]	1374	! HYBRID :
[454]	1375	DO jj = 1, jpj
	1376	DO ji = 1, jpi
	1377	DO jk = 1, jpkm1
	1378	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
[2528]	1379	IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0
[454]	1380	END DO
	1381	END DO
	1382	END DO
[1099]	1383	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
	1384	& ' MAX ', MAXVAL( mbathy(:,:) )
[454]	1385
[1099]	1386	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
	1387	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	1388	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
	1389	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) )
	1390	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), &
	1391	& ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), &
	1392	& ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), &
	1393	& ' w ', MINVAL( fse3w (:,:,:) )
[454]	1394
[1099]	1395	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
	1396	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) )
	1397	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), &
	1398	& ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), &
	1399	& ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), &
	1400	& ' w ', MAXVAL( fse3w (:,:,:) )
	1401	ENDIF
[3600]	1402	! END DO
[1099]	1403	IF(lwp) THEN ! selected vertical profiles
[454]	1404	WRITE(numout,*)
	1405	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
	1406	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1407	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1408	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), &
	1409	& fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
[473]	1410	DO jj = mj0(20), mj1(20)
	1411	DO ji = mi0(20), mi1(20)
	1412	WRITE(numout,*)
	1413	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	1414	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1415	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1416	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
	1417	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
[473]	1418	END DO
	1419	END DO
	1420	DO jj = mj0(74), mj1(74)
	1421	DO ji = mi0(100), mi1(100)
	1422	WRITE(numout,*)
	1423	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	1424	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1425	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1426	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
	1427	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
[473]	1428	END DO
	1429	END DO
[454]	1430	ENDIF
	1431
[3600]	1432	!================================================================================
	1433	! check the coordinate makes sense
	1434	!================================================================================
	1435	DO ji = 1, jpi
[454]	1436	DO jj = 1, jpj
[3600]	1437
	1438	IF( hbatt(ji,jj) > 0._wp) THEN
	1439	DO jk = 1, mbathy(ji,jj)
	1440	! check coordinate is monotonically increasing
	1441	IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
	1442	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
	1443	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
	1444	WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
	1445	WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
	1446	CALL ctl_stop( ctmp1 )
	1447	ENDIF
	1448	! and check it has never gone negative
	1449	IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
	1450	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
	1451	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
	1452	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
	1453	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
	1454	CALL ctl_stop( ctmp1 )
	1455	ENDIF
	1456	! and check it never exceeds the total depth
	1457	IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
	1458	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
	1459	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
	1460	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
	1461	CALL ctl_stop( ctmp1 )
	1462	ENDIF
	1463	END DO
	1464
	1465	DO jk = 1, mbathy(ji,jj)-1
	1466	! and check it never exceeds the total depth
	1467	IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
	1468	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
	1469	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
	1470	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
	1471	CALL ctl_stop( ctmp1 )
	1472	ENDIF
	1473	END DO
	1474
	1475	ENDIF
	1476
[454]	1477	END DO
	1478	END DO
[1099]	1479	!
[3294]	1480	CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat )
[2715]	1481	!
[3294]	1482	IF( nn_timing == 1 ) CALL timing_stop('zgr_sco')
	1483	!
[454]	1484	END SUBROUTINE zgr_sco
	1485
[3600]	1486	!!======================================================================
	1487	SUBROUTINE s_sh94()
	1488
	1489	!!----------------------------------------------------------------------
	1490	!! * ROUTINE s_sh94 *
	1491	!!
	1492	!! ** Purpose : stretch the s-coordinate system
	1493	!!
	1494	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
	1495	!! mixed S/sigma coordinate
	1496	!!
	1497	!! Reference : Song and Haidvogel 1994.
	1498	!!----------------------------------------------------------------------
	1499	!
	1500	INTEGER :: ji, jj, jk ! dummy loop argument
	1501	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
	1502	!
[3618]	1503	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
	1504	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
[3600]	1505
[3618]	1506	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
	1507	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
[3600]	1508
[3618]	1509	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
	1510	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
	1511	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
	1512	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
[3600]	1513
	1514	DO ji = 1, jpi
	1515	DO jj = 1, jpj
	1516
	1517	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
	1518	DO jk = 1, jpk
[3618]	1519	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
	1520	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
[3600]	1521	END DO
	1522	ELSE ! shallow water, uniform sigma
	1523	DO jk = 1, jpk
[3618]	1524	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
	1525	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
[3600]	1526	END DO
	1527	ENDIF
	1528	!
	1529	DO jk = 1, jpkm1
[3618]	1530	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
	1531	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
[3600]	1532	END DO
[3618]	1533	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
	1534	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
[3600]	1535	!
	1536	! Coefficients for vertical depth as the sum of e3w scale factors
[3618]	1537	z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
[3600]	1538	DO jk = 2, jpk
[3618]	1539	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
[3600]	1540	END DO
	1541	!
	1542	DO jk = 1, jpk
	1543	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
	1544	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
[3618]	1545	gdept (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
	1546	gdepw (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
	1547	gdep3w(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsi3w3(ji,jj,jk)+rn_hczcoeft )
[3600]	1548	END DO
	1549	!
	1550	END DO ! for all jj's
	1551	END DO ! for all ji's
	1552
	1553	DO ji = 1, jpim1
	1554	DO jj = 1, jpjm1
	1555	DO jk = 1, jpk
[3618]	1556	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
[3600]	1557	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
[3618]	1558	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
[3600]	1559	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
[3618]	1560	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
	1561	& + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
[3600]	1562	& / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
[3618]	1563	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
[3600]	1564	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
[3618]	1565	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
[3600]	1566	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
	1567	!
[3618]	1568	e3t(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1569	e3u(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1570	e3v(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1571	e3f(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
[3600]	1572	!
[3618]	1573	e3w (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1574	e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1575	e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
[3600]	1576	END DO
	1577	END DO
	1578	END DO
	1579
[3618]	1580	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
	1581	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
[3600]	1582
	1583	END SUBROUTINE s_sh94
	1584
	1585	SUBROUTINE s_sf12
	1586
	1587	!!----------------------------------------------------------------------
	1588	!! * ROUTINE s_sf12 *
	1589	!!
	1590	!! ** Purpose : stretch the s-coordinate system
	1591	!!
	1592	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
	1593	!! mixed S/sigma/Z coordinate
	1594	!!
	1595	!! This method allows the maintenance of fixed surface and or
	1596	!! bottom cell resolutions (cf. geopotential coordinates)
	1597	!! within an analytically derived stretched S-coordinate framework.
	1598	!!
	1599	!!
	1600	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
	1601	!!----------------------------------------------------------------------
	1602	!
	1603	INTEGER :: ji, jj, jk ! dummy loop argument
[3618]	1604	REAL(wp) :: zsmth ! smoothing around critical depth
	1605	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
[3600]	1606	!
[3618]	1607	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
	1608	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
[3600]	1609
	1610	!
[3618]	1611	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
	1612	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
[3600]	1613
[3618]	1614	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
	1615	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
	1616	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
	1617	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
[3600]	1618
	1619	DO ji = 1, jpi
	1620	DO jj = 1, jpj
	1621
	1622	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
	1623
[3618]	1624	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
[3600]	1625	! could be changed by users but care must be taken to do so carefully
[3618]	1626	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
[3600]	1627
[3618]	1628	zzs = rn_zs / hbatt(ji,jj)
[3600]	1629
	1630	IF (rn_efold /= 0.0_wp) THEN
[3618]	1631	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
[3600]	1632	ELSE
[3618]	1633	zsmth = 1.0_wp
[3600]	1634	ENDIF
	1635
	1636	DO jk = 1, jpk
[3618]	1637	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
	1638	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
[3600]	1639	ENDDO
[3618]	1640	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
	1641	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
[3600]	1642
	1643	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
	1644
	1645	DO jk = 1, jpk
[3618]	1646	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
	1647	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
[3600]	1648	END DO
	1649
	1650	ELSE ! shallow water, z coordinates
	1651
	1652	DO jk = 1, jpk
[3618]	1653	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
	1654	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
[3600]	1655	END DO
	1656
	1657	ENDIF
	1658
	1659	DO jk = 1, jpkm1
[3618]	1660	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
	1661	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
[3600]	1662	END DO
[3618]	1663	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
	1664	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
[3600]	1665
	1666	! Coefficients for vertical depth as the sum of e3w scale factors
[3618]	1667	z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
[3600]	1668	DO jk = 2, jpk
[3618]	1669	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
[3600]	1670	END DO
	1671
	1672	DO jk = 1, jpk
[3618]	1673	gdept (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
	1674	gdepw (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
	1675	gdep3w(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
[3600]	1676	END DO
	1677
	1678	ENDDO ! for all jj's
	1679	ENDDO ! for all ji's
	1680
	1681	DO ji=1,jpi
	1682	DO jj=1,jpj
	1683
	1684	DO jk = 1, jpk
[3618]	1685	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) / &
[3600]	1686	( hbatt(ji,jj)+hbatt(ji+1,jj) )
[3618]	1687	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) / &
[3600]	1688	( hbatt(ji,jj)+hbatt(ji,jj+1) )
[3618]	1689	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) + &
	1690	hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) / &
[3600]	1691	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
[3618]	1692	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) / &
[3600]	1693	( hbatt(ji,jj)+hbatt(ji+1,jj) )
[3618]	1694	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) / &
[3600]	1695	( hbatt(ji,jj)+hbatt(ji,jj+1) )
	1696
[3618]	1697	e3t(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
	1698	e3u(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
	1699	e3v(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
	1700	e3f(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
[3600]	1701	!
[3618]	1702	e3w(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
	1703	e3uw(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
	1704	e3vw(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
[3600]	1705	END DO
	1706
	1707	ENDDO
	1708	ENDDO
[3618]	1709	! ! =============
[3600]	1710
[3618]	1711	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
	1712	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
[3600]	1713
	1714	END SUBROUTINE s_sf12
	1715
	1716	SUBROUTINE s_tanh()
	1717
	1718	!!----------------------------------------------------------------------
	1719	!! * ROUTINE s_tanh*
	1720	!!
	1721	!! ** Purpose : stretch the s-coordinate system
	1722	!!
	1723	!! ** Method : s-coordinate stretch
	1724	!!
	1725	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
	1726	!!----------------------------------------------------------------------
	1727
	1728	INTEGER :: ji, jj, jk ! dummy loop argument
	1729	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
	1730
[3618]	1731	REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
	1732	REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
	1733
	1734	CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
	1735	CALL wrk_alloc( jpk, z_esigt, z_esigw )
	1736
	1737	z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp
	1738	z_esigt = 0._wp ; z_esigw = 0._wp
	1739
[3600]	1740	DO jk = 1, jpk
[3618]	1741	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
	1742	z_gsigt(jk) = -fssig( REAL(jk,wp) )
[3600]	1743	END DO
[3618]	1744	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
[3600]	1745	!
	1746	! Coefficients for vertical scale factors at w-, t- levels
	1747	!!gm bug : define it from analytical function, not like juste bellow....
	1748	!!gm or betteroffer the 2 possibilities....
	1749	DO jk = 1, jpkm1
[3618]	1750	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
	1751	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
[3600]	1752	END DO
[3618]	1753	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
	1754	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
[3600]	1755	!
	1756	! Coefficients for vertical depth as the sum of e3w scale factors
[3618]	1757	z_gsi3w(1) = 0.5_wp * z_esigw(1)
[3600]	1758	DO jk = 2, jpk
[3618]	1759	z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
[3600]	1760	END DO
	1761	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
	1762	DO jk = 1, jpk
	1763	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
	1764	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
[3618]	1765	gdept (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
	1766	gdepw (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
	1767	gdep3w(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsi3w(jk) + hift(:,:)zcoeft )
[3600]	1768	END DO
	1769	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
	1770	DO jj = 1, jpj
	1771	DO ji = 1, jpi
	1772	DO jk = 1, jpk
[3618]	1773	e3t(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
	1774	e3u(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
	1775	e3v(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
	1776	e3f(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
[3600]	1777	!
[3618]	1778	e3w (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
	1779	e3uw(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
	1780	e3vw(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
[3600]	1781	END DO
	1782	END DO
	1783	END DO
[3618]	1784
	1785	CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
	1786	CALL wrk_dealloc( jpk, z_esigt, z_esigw )
	1787
[3600]	1788	END SUBROUTINE s_tanh
	1789
	1790	FUNCTION fssig( pk ) RESULT( pf )
	1791	!!----------------------------------------------------------------------
	1792	!! * ROUTINE fssig *
	1793	!!
	1794	!! ** Purpose : provide the analytical function in s-coordinate
	1795	!!
	1796	!! ** Method : the function provide the non-dimensional position of
	1797	!! T and W (i.e. between 0 and 1)
	1798	!! T-points at integer values (between 1 and jpk)
	1799	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1800	!!----------------------------------------------------------------------
	1801	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
	1802	REAL(wp) :: pf ! sigma value
	1803	!!----------------------------------------------------------------------
	1804	!
	1805	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
	1806	& - TANH( rn_thetb * rn_theta ) ) &
	1807	& * ( COSH( rn_theta ) &
	1808	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
	1809	& / ( 2._wp * SINH( rn_theta ) )
	1810	!
	1811	END FUNCTION fssig
	1812
	1813
	1814	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
	1815	!!----------------------------------------------------------------------
	1816	!! * ROUTINE fssig1 *
	1817	!!
	1818	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
	1819	!!
	1820	!! ** Method : the function provides the non-dimensional position of
	1821	!! T and W (i.e. between 0 and 1)
	1822	!! T-points at integer values (between 1 and jpk)
	1823	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1824	!!----------------------------------------------------------------------
	1825	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
	1826	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
	1827	REAL(wp) :: pf1 ! sigma value
	1828	!!----------------------------------------------------------------------
	1829	!
	1830	IF ( rn_theta == 0 ) then ! uniform sigma
	1831	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
	1832	ELSE ! stretched sigma
	1833	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
	1834	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
	1835	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
	1836	ENDIF
	1837	!
	1838	END FUNCTION fssig1
	1839
	1840
[3618]	1841	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
[3600]	1842	!!----------------------------------------------------------------------
	1843	!! * ROUTINE fgamma *
	1844	!!
	1845	!! ** Purpose : provide analytical function for the s-coordinate
	1846	!!
	1847	!! ** Method : the function provides the non-dimensional position of
	1848	!! T and W (i.e. between 0 and 1)
	1849	!! T-points at integer values (between 1 and jpk)
	1850	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1851	!!
	1852	!! This method allows the maintenance of fixed surface and or
	1853	!! bottom cell resolutions (cf. geopotential coordinates)
	1854	!! within an analytically derived stretched S-coordinate framework.
	1855	!!
	1856	!! Reference : Siddorn and Furner, in prep
	1857	!!----------------------------------------------------------------------
[3618]	1858	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
	1859	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
	1860	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
	1861	REAL(wp), INTENT(in ) :: pzs ! surface box depth
	1862	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
	1863	REAL(wp) :: za1,za2,za3 ! local variables
	1864	REAL(wp) :: zn1,zn2 ! local variables
	1865	REAL(wp) :: za,zb,zx ! local variables
[3600]	1866	integer :: jk
	1867	!!----------------------------------------------------------------------
	1868	!
	1869
[3618]	1870	zn1 = 1./(jpk-1.)
	1871	zn2 = 1. - zn1
[3600]	1872
[3618]	1873	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
	1874	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
	1875	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
[3600]	1876
[3618]	1877	za = pzb - za3*(pzs-za1)-za2
	1878	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
	1879	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
	1880	zx = 1.0_wp-za/2.0_wp-zb
[3600]	1881
	1882	DO jk = 1, jpk
[3618]	1883	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(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) )
	1884	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
[3600]	1885	ENDDO
	1886
	1887	!
	1888	END FUNCTION fgamma
	1889
[3]	1890	!!======================================================================
	1891	END MODULE domzgr

Note: See TracBrowser for help on using the repository browser.

Download in other formats: