Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

domzgr.F90 @ 12871

Last change on this file since 12871 was 12871, checked in by mathiot, 4 years ago
ticket #2457: add fix suggested in ticket.
File size: 98.1 KB

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

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

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

Context Navigation

source: utils/tools_ticket2457/DOMAINcfg/src/domzgr.F90 @ 12871

Download in other formats: