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 NEMO/branches/2019/ENHANCE-03_domcfg/src – NEMO

Context Navigation

source: NEMO/branches/2019/ENHANCE-03_domcfg/src/domzgr.F90 @ 11129

Last change on this file since 11129 was 11129, checked in by mathiot, 5 years ago
simplification of domcfg (rm all var_n and var_b as it is not needed) (ticket #2143)
File size: 136.4 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
[10727]	19	!! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye
[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
[10727]	38	! USE closea ! closed seas
[6951]	39	!
	40	USE in_out_manager ! I/O manager
	41	USE iom ! I/O library
	42	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	43	USE lib_mpp ! distributed memory computing library
[11129]	44	USE lib_fortran
[10727]	45	USE dombat
[6951]	46
	47	IMPLICIT NONE
	48	PRIVATE
	49
	50	PUBLIC dom_zgr ! called by dom_init.F90
	51
	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	INTEGER , PUBLIC :: ntopo !: = 0/1 ,compute/read the bathymetry file
	61	REAL(wp), PUBLIC :: e3zps_min !: miminum thickness for partial steps (meters)
	62	REAL(wp), PUBLIC :: e3zps_rat !: minimum thickness ration for partial steps
[11129]	63	INTEGER, PUBLIC :: nperio !: type of lateral boundary condition
[10727]	64
[6951]	65	! Song and Haidvogel 1994 stretching parameters
	66	REAL(wp) :: rn_theta ! surface control parameter (0<=rn_theta<=20)
	67	REAL(wp) :: rn_thetb ! bottom control parameter (0<=rn_thetb<= 1)
	68	REAL(wp) :: rn_bb ! stretching parameter
	69	! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
	70	! Siddorn and Furner stretching parameters
	71	LOGICAL :: ln_sigcrit ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
	72	REAL(wp) :: rn_alpha ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
	73	REAL(wp) :: rn_efold ! efold length scale for transition to stretched coord
	74	REAL(wp) :: rn_zs ! depth of surface grid box
	75	! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
	76	REAL(wp) :: rn_zb_a ! bathymetry scaling factor for calculating Zb
	77	REAL(wp) :: rn_zb_b ! offset for calculating Zb
	78
	79	!! * Substitutions
	80	!!----------------------------------------------------------------------
	81	!! * vectopt_loop_substitute *
	82	!!----------------------------------------------------------------------
	83	!! ** purpose : substitute the inner loop start/end indices with CPP macro
	84	!! allow unrolling of do-loop (useful with vector processors)
	85	!!----------------------------------------------------------------------
	86	!!----------------------------------------------------------------------
	87	!! NEMO/OPA 3.7 , NEMO Consortium (2014)
	88	!! $Id: vectopt_loop_substitute.h90 4990 2014-12-15 16:42:49Z timgraham $
[9598]	89	!! Software governed by the CeCILL licence (./LICENSE)
[6951]	90	!!----------------------------------------------------------------------
	91	!!----------------------------------------------------------------------
	92	!! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
[7107]	93	!! $Id: domzgr.F90 6827 2016-08-01 13:37:15Z flavoni $
[9598]	94	!! Software governed by the CeCILL licence (./LICENSE)
[6951]	95	!!----------------------------------------------------------------------
	96	CONTAINS
	97
	98	SUBROUTINE dom_zgr
	99	!!----------------------------------------------------------------------
	100	!! * ROUTINE dom_zgr *
	101	!!
	102	!! ** Purpose : set the depth of model levels and the resulting
	103	!! vertical scale factors.
	104	!!
	105	!! ** Method : - reference 1D vertical coordinate (gdep._1d, e3._1d)
	106	!! - read/set ocean depth and ocean levels (bathy, mbathy)
	107	!! - vertical coordinate (gdep., e3.) depending on the
	108	!! coordinate chosen :
	109	!! ln_zco=T z-coordinate
	110	!! ln_zps=T z-coordinate with partial steps
	111	!! ln_zco=T s-coordinate
	112	!!
	113	!! ** Action : define gdep., e3., mbathy and bathy
	114	!!----------------------------------------------------------------------
	115	INTEGER :: ioptio, ibat ! local integer
	116	INTEGER :: ios
	117	!
	118	NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh
	119	!!----------------------------------------------------------------------
	120	!
	121	!
	122	REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate
	123	READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
	124	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
	125
	126	REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate
	127	READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
	128	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
	129	IF(lwm) WRITE ( numond, namzgr )
	130
	131	IF(lwp) THEN ! Control print
	132	WRITE(numout,*)
	133	WRITE(numout,*) 'dom_zgr : vertical coordinate'
	134	WRITE(numout,*) '~~~~~~~'
	135	WRITE(numout,*) ' Namelist namzgr : set vertical coordinate'
	136	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
	137	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
	138	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
	139	WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav
	140	WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh
	141	ENDIF
	142
	143	IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time'
	144
	145	ioptio = 0 ! Check Vertical coordinate options
	146	IF( ln_zco ) ioptio = ioptio + 1
	147	IF( ln_zps ) ioptio = ioptio + 1
	148	IF( ln_sco ) ioptio = ioptio + 1
	149	IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
	150	!
	151	ioptio = 0
	152	IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1
	153	IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1
	154	IF( ioptio > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' )
	155	!
	156	! Build the vertical coordinate system
	157	! ------------------------------------
	158	CALL zgr_z ! Reference z-coordinate system (always called)
	159	CALL zgr_bat ! Bathymetry fields (levels and meters)
	160	IF( ln_zco ) CALL zgr_zco ! z-coordinate
	161	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
	162	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
	163	!
	164	! final adjustment of mbathy & check
	165	! -----------------------------------
	166	IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain
	167	CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points
	168	CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
	169	CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points
	170	!
	171	IF( nprint == 1 .AND. lwp ) THEN
	172	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	173	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
	174	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) )
	175	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), &
	176	& ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), &
	177	& ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), &
	178	& ' w ', MINVAL( e3w_0(:,:,:) )
	179
	180	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
	181	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) )
	182	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), &
	183	& ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), &
	184	& ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), &
	185	& ' w ', MAXVAL( e3w_0(:,:,:) )
	186	ENDIF
	187	!
	188	END SUBROUTINE dom_zgr
	189
	190
	191	SUBROUTINE zgr_z
	192	!!----------------------------------------------------------------------
	193	!! * ROUTINE zgr_z *
	194	!!
	195	!! ** Purpose : set the depth of model levels and the resulting
	196	!! vertical scale factors.
	197	!!
	198	!! ** Method : z-coordinate system (use in all type of coordinate)
	199	!! The depth of model levels is defined from an analytical
	200	!! function the derivative of which gives the scale factors.
	201	!! both depth and scale factors only depend on k (1d arrays).
	202	!! w-level: gdepw_1d = gdep(k)
	203	!! e3w_1d(k) = dk(gdep)(k) = e3(k)
	204	!! t-level: gdept_1d = gdep(k+0.5)
	205	!! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
	206	!!
	207	!! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
	208	!! - e3t_1d , e3w_1d : scale factors at T- and W-levels (m)
	209	!!
	210	!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
	211	!!----------------------------------------------------------------------
	212	INTEGER :: jk ! dummy loop indices
	213	REAL(wp) :: zt, zw ! temporary scalars
	214	REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
	215	REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
	216	REAL(wp) :: zrefdep ! depth of the reference level (~10m)
	217	REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
	218	!!----------------------------------------------------------------------
	219	!
	220	! Set variables from parameters
	221	! ------------------------------
	222	zkth = ppkth ; zacr = ppacr
	223	zdzmin = ppdzmin ; zhmax = pphmax
	224	zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters
	225
	226	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
	227	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
	228	IF( ppa1 == pp_to_be_computed .AND. &
	229	& ppa0 == pp_to_be_computed .AND. &
	230	& ppsur == pp_to_be_computed ) THEN
	231	!
	232	za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
	233	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
	234	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
	235	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
	236	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
	237	ELSE
	238	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
	239	za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter
	240	ENDIF
	241
	242	IF(lwp) THEN ! Parameter print
	243	WRITE(numout,*)
	244	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
	245	WRITE(numout,*) ' ~~~~~~~'
	246	IF( ppkth == 0._wp ) THEN
	247	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
	248	WRITE(numout,*) ' Total depth :', zhmax
	249	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
	250	ELSE
	251	IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
	252	WRITE(numout,*) ' zsur, za0, za1 computed from '
	253	WRITE(numout,*) ' zdzmin = ', zdzmin
	254	WRITE(numout,*) ' zhmax = ', zhmax
	255	ENDIF
	256	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
	257	WRITE(numout,*) ' zsur = ', zsur
	258	WRITE(numout,*) ' za0 = ', za0
	259	WRITE(numout,*) ' za1 = ', za1
	260	WRITE(numout,*) ' zkth = ', zkth
	261	WRITE(numout,*) ' zacr = ', zacr
	262	IF( ldbletanh ) THEN
	263	WRITE(numout,*) ' (Double tanh za2 = ', za2
	264	WRITE(numout,*) ' parameters) zkth2= ', zkth2
	265	WRITE(numout,*) ' zacr2= ', zacr2
	266	ENDIF
	267	ENDIF
	268	ENDIF
	269
	270
	271	! Reference z-coordinate (depth - scale factor at T- and W-points)
	272	! ======================
	273	IF( ppkth == 0._wp ) THEN ! uniform vertical grid
	274
	275
	276
	277	za1 = zhmax / FLOAT(jpk-1)
	278
	279	DO jk = 1, jpk
	280	zw = FLOAT( jk )
	281	zt = FLOAT( jk ) + 0.5_wp
	282	gdepw_1d(jk) = ( zw - 1 ) * za1
	283	gdept_1d(jk) = ( zt - 1 ) * za1
	284	e3w_1d (jk) = za1
	285	e3t_1d (jk) = za1
	286	END DO
	287	ELSE ! Madec & Imbard 1996 function
	288	IF( .NOT. ldbletanh ) THEN
	289	DO jk = 1, jpk
	290	zw = REAL( jk , wp )
	291	zt = REAL( jk , wp ) + 0.5_wp
	292	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
	293	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
	294	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
	295	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
	296	END DO
	297	ELSE
	298	DO jk = 1, jpk
	299	zw = FLOAT( jk )
	300	zt = FLOAT( jk ) + 0.5_wp
	301	! Double tanh function
	302	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) &
	303	& + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) )
	304	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) &
	305	& + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) )
	306	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) &
	307	& + za2 * TANH( (zw-zkth2) / zacr2 )
	308	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) &
	309	& + za2 * TANH( (zt-zkth2) / zacr2 )
	310	END DO
	311	ENDIF
	312	gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero
	313	ENDIF
	314
[7200]	315	IF ( ln_isfcav .OR. ln_e3_dep ) THEN ! e3. = dk[gdep]
[7189]	316	!
	317	!==>>> need to be like this to compute the pressure gradient with ISF.
	318	! If not, level beneath the ISF are not aligned (sum(e3t) /= depth)
	319	! define e3t_0 and e3w_0 as the differences between gdept and gdepw respectively
	320	!
[6951]	321	DO jk = 1, jpkm1
	322	e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)
	323	END DO
	324	e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO
	325
	326	DO jk = 2, jpk
	327	e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)
	328	END DO
	329	e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))
	330	END IF
	331
	332	!!gm BUG in s-coordinate this does not work!
	333	! deepest/shallowest W level Above/Below ~10m
	334	zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
	335	nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
	336	nla10 = nlb10 - 1 ! deepest W level Above ~10m
	337	!!gm end bug
	338
	339	IF(lwp) THEN ! control print
	340	WRITE(numout,*)
	341	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
	342	WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
	343	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
	344	ENDIF
	345	DO jk = 1, jpk ! control positivity
	346	IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 ' )
	347	IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
	348	END DO
	349	!
	350	END SUBROUTINE zgr_z
	351
	352
	353	SUBROUTINE zgr_bat
	354	!!----------------------------------------------------------------------
	355	!! * ROUTINE zgr_bat *
	356	!!
	357	!! ** Purpose : set bathymetry both in levels and meters
	358	!!
	359	!! ** Method : read or define mbathy and bathy arrays
	360	!! * level bathymetry:
	361	!! The ocean basin geometry is given by a two-dimensional array,
	362	!! mbathy, which is defined as follow :
	363	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
	364	!! at t-point (ji,jj).
	365	!! = 0 over the continental t-point.
	366	!! The array mbathy is checked to verified its consistency with
	367	!! model option. in particular:
	368	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	369	!! along closed boundary.
	370	!! mbathy must be cyclic IF jperio=1.
	371	!! mbathy must be lower or equal to jpk-1.
	372	!! isolated ocean grid points are suppressed from mbathy
	373	!! since they are only connected to remaining
	374	!! ocean through vertical diffusion.
	375	!! ntopo=-1 : rectangular channel or bassin with a bump
	376	!! ntopo= 0 : flat rectangular channel or basin
	377	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
	378	!! bathy is read in 'bathy_meter.nc' NetCDF file
	379	!!
	380	!! ** Action : - mbathy: level bathymetry (in level index)
	381	!! - bathy : meter bathymetry (in meters)
	382	!!----------------------------------------------------------------------
	383	INTEGER :: ji, jj, jk ! dummy loop indices
	384	INTEGER :: inum ! temporary logical unit
	385	INTEGER :: ierror ! error flag
	386	INTEGER :: ii_bump, ij_bump, ih ! bump center position
	387	INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices
	388	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
	389	REAL(wp) :: zi, zj, zh, zhmin ! local scalars
	390	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: idta ! global domain integer data
	391	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! global domain scalar data
	392	!!----------------------------------------------------------------------
	393	!
	394	IF(lwp) WRITE(numout,*)
	395	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
	396	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
	397	! ! ================== !
	398	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
	399	! ! ================== !
	400	! ! global domain level and meter bathymetry (idta,zdta)
	401	!
[11129]	402	ALLOCATE( idta(jpiglo,jpjglo), STAT=ierror )
[6951]	403	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' )
[11129]	404	ALLOCATE( zdta(jpiglo,jpjglo), STAT=ierror )
[6951]	405	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' )
	406	!
	407	IF( ntopo == 0 ) THEN ! flat basin
	408	IF(lwp) WRITE(numout,*)
	409	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
	410	IF( rn_bathy > 0.01 ) THEN
	411	IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist'
	412	zdta(:,:) = rn_bathy
	413	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
	414	idta(:,:) = jpkm1
	415	ELSE ! z-coordinate (zco or zps): step-like topography
	416	idta(:,:) = jpkm1
	417	DO jk = 1, jpkm1
	418	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
	419	END DO
	420	ENDIF
	421	ELSE
	422	IF(lwp) WRITE(numout,*) ' Depth = depthw(jpkm1)'
	423	idta(:,:) = jpkm1 ! before last level
	424	zdta(:,:) = gdepw_1d(jpk) ! last w-point depth
	425	h_oce = gdepw_1d(jpk)
	426	ENDIF
	427	ELSE ! bump centered in the basin
	428	IF(lwp) WRITE(numout,*)
	429	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
[11129]	430	ii_bump = jpiglo / 2 ! i-index of the bump center
	431	ij_bump = jpjglo / 2 ! j-index of the bump center
[6951]	432	r_bump = 50000._wp ! bump radius (meters)
	433	h_bump = 2700._wp ! bump height (meters)
	434	h_oce = gdepw_1d(jpk) ! background ocean depth (meters)
	435	IF(lwp) WRITE(numout,*) ' bump characteristics: '
	436	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
	437	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
	438	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
	439	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
	440	!
[11129]	441	DO jj = 1, jpjglo ! zdta :
	442	DO ji = 1, jpiglo
[6951]	443	zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
	444	zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
	445	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
	446	END DO
	447	END DO
	448	! ! idta :
	449	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
	450	idta(:,:) = jpkm1
	451	ELSE ! z-coordinate (zco or zps): step-like topography
	452	idta(:,:) = jpkm1
	453	DO jk = 1, jpkm1
	454	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
	455	END DO
	456	ENDIF
	457	ENDIF
[11129]	458	!
[6951]	459	! ! set GLOBAL boundary conditions
	460	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
	461	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp
[11129]	462	idta( : ,jpjglo) = 0 ; zdta( : ,jpjglo) = 0._wp
[6951]	463	ELSEIF( jperio == 2 ) THEN
	464	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
[11129]	465	idta( : ,jpjglo) = 0 ; zdta( : ,jpjglo) = 0._wp
[6951]	466	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp
[11129]	467	idta(jpiglo, : ) = 0 ; zdta(jpiglo, : ) = 0._wp
[6951]	468	ELSE
	469	ih = 0 ; zh = 0._wp
	470	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
	471	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
[11129]	472	idta( : ,jpjglo) = ih ; zdta( : ,jpjglo) = zh
[6951]	473	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
[11129]	474	idta(jpiglo, : ) = ih ; zdta(jpiglo, : ) = zh
[6951]	475	ENDIF
	476
	477	! ! local domain level and meter bathymetries (mbathy,bathy)
	478	mbathy(:,:) = 0 ! set to zero extra halo points
	479	bathy (:,:) = 0._wp ! (require for mpp case)
	480	DO jj = 1, nlcj ! interior values
	481	DO ji = 1, nlci
	482	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
	483	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
	484	END DO
	485	END DO
	486	risfdep(:,:)=0.e0
	487	misfdep(:,:)=1
	488	!
[7107]	489	! (ISF) TODO build ice draft netcdf file for isomip and build the corresponding part of code
	490	IF( cp_cfg == "isomip" .AND. ln_isfcav ) THEN
	491	risfdep(:,:)=200.e0
	492	misfdep(:,:)=1
	493	ij0 = 1 ; ij1 = 40
	494	DO jj = mj0(ij0), mj1(ij1)
	495	risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
	496	END DO
	497	WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp
	498	!
	499	ELSEIF ( cp_cfg == "isomip2" .AND. ln_isfcav ) THEN
	500	!
	501	risfdep(:,:)=0.e0
	502	misfdep(:,:)=1
	503	ij0 = 1 ; ij1 = 40
	504	DO jj = mj0(ij0), mj1(ij1)
	505	risfdep(:,jj)=700.0_wp-(gphit(:,jj)+80.0_wp)*125.0_wp
	506	END DO
	507	WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp
	508	END IF
	509	!
[6951]	510	DEALLOCATE( idta, zdta )
	511	!
	512	! ! ================ !
[10727]	513	ELSEIF( ntopo == 1 .OR. ntopo ==2 ) THEN ! read in file ! (over the local domain)
[6951]	514	! ! ================ !
	515	!
	516	IF( ln_zco ) THEN ! zco : read level bathymetry
	517	CALL iom_open ( 'bathy_level.nc', inum )
	518	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
	519	CALL iom_close( inum )
	520	mbathy(:,:) = INT( bathy(:,:) )
	521	! initialisation isf variables
	522	risfdep(:,:)=0._wp ; misfdep(:,:)=1
	523	! ! =====================
	524	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	525	! ! =====================
	526	!
	527	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	528	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
	529	DO ji = mi0(ii0), mi1(ii1)
	530	DO jj = mj0(ij0), mj1(ij1)
	531	mbathy(ji,jj) = 15
	532	END DO
	533	END DO
	534	IF(lwp) WRITE(numout,*)
	535	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
	536	!
	537	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	538	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
	539	DO ji = mi0(ii0), mi1(ii1)
	540	DO jj = mj0(ij0), mj1(ij1)
	541	mbathy(ji,jj) = 12
	542	END DO
	543	END DO
	544	IF(lwp) WRITE(numout,*)
	545	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
	546	!
	547	ENDIF
	548	!
	549	ENDIF
	550	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
[10727]	551	#if defined key_agrif
	552	IF (agrif_root()) THEN
	553	#endif
	554	IF( ntopo == 1) THEN
	555	CALL iom_open ( 'bathy_meter.nc', inum )
	556	IF ( ln_isfcav ) THEN
	557	CALL iom_get ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
	558	ELSE
	559	CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr )
	560	END IF
	561	CALL iom_close( inum )
[6951]	562	ELSE
[10727]	563	CALL dom_bat
	564	ENDIF
	565	#if defined key_agrif
	566	ELSE
	567	IF( ntopo == 1) THEN
	568	CALL agrif_create_bathy_meter()
	569	ELSE
	570	CALL dom_bat
	571	ENDIF
	572	ENDIF
	573	#endif
[6951]	574	!
	575	! initialisation isf variables
	576	risfdep(:,:)=0._wp ; misfdep(:,:)=1
	577	!
	578	IF ( ln_isfcav ) THEN
	579	CALL iom_open ( 'isf_draft_meter.nc', inum )
	580	CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep )
	581	CALL iom_close( inum )
	582	WHERE( bathy(:,:) <= 0._wp ) risfdep(:,:) = 0._wp
	583
	584	! set grounded point to 0
	585	! (a treshold could be set here if needed, or set it offline based on the grounded fraction)
	586	WHERE ( bathy(:,:) <= risfdep(:,:) + rn_isfhmin )
	587	misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
	588	mbathy (:,:) = 0 ; bathy (:,:) = 0._wp
	589	END WHERE
	590	END IF
	591	!
	592	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	593	!
	594	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	595	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
	596	DO ji = mi0(ii0), mi1(ii1)
	597	DO jj = mj0(ij0), mj1(ij1)
	598	bathy(ji,jj) = 284._wp
	599	END DO
	600	END DO
	601	IF(lwp) WRITE(numout,*)
	602	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
	603	!
	604	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	605	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
	606	DO ji = mi0(ii0), mi1(ii1)
	607	DO jj = mj0(ij0), mj1(ij1)
	608	bathy(ji,jj) = 137._wp
	609	END DO
	610	END DO
	611	IF(lwp) WRITE(numout,*)
	612	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
	613	!
	614	ENDIF
	615	!
	616	ENDIF
	617	! ! =============== !
	618	ELSE ! error !
	619	! ! =============== !
	620	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
	621	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
	622	ENDIF
	623	!
[10727]	624	! IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==!
[6951]	625	!
	626	IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
	627	IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
	628	ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth
	629	ENDIF
	630	zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels
	631	WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
	632	ELSE WHERE ; bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
	633	END WHERE
	634	IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
	635	ENDIF
	636	!
	637	END SUBROUTINE zgr_bat
	638
	639
	640	SUBROUTINE zgr_bat_zoom
	641	!!----------------------------------------------------------------------
	642	!! * ROUTINE zgr_bat_zoom *
	643	!!
	644	!! ** Purpose : - Close zoom domain boundary if necessary
	645	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
	646	!!
	647	!! ** Method :
	648	!!
	649	!! ** Action : - update mbathy: level bathymetry (in level index)
	650	!!----------------------------------------------------------------------
	651	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
	652	!!----------------------------------------------------------------------
	653	!
	654	IF(lwp) WRITE(numout,*)
	655	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
	656	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
	657	!
	658	! Zoom domain
	659	! ===========
	660	!
	661	! Forced closed boundary if required
	662	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
	663	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
	664	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
	665	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
	666	!
	667	! Configuration specific domain modifications
	668	! (here, ORCA arctic configuration: suppress Med Sea)
	669	IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
	670	SELECT CASE ( jp_cfg )
	671	! ! =======================
	672	CASE ( 2 ) ! ORCA_R2 configuration
	673	! ! =======================
	674	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
	675	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
	676	ij0 = 98 ; ij1 = 110
	677	! ! =======================
	678	CASE ( 05 ) ! ORCA_R05 configuration
	679	! ! =======================
	680	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
	681	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
	682	ij0 = 314 ; ij1 = 370
	683	END SELECT
	684	!
	685	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
	686	!
	687	ENDIF
	688	!
	689	END SUBROUTINE zgr_bat_zoom
	690
	691
	692	SUBROUTINE zgr_bat_ctl
	693	!!----------------------------------------------------------------------
	694	!! * ROUTINE zgr_bat_ctl *
	695	!!
	696	!! ** Purpose : check the bathymetry in levels
	697	!!
	698	!! ** Method : The array mbathy is checked to verified its consistency
	699	!! with the model options. in particular:
	700	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	701	!! along closed boundary.
	702	!! mbathy must be cyclic IF jperio=1.
	703	!! mbathy must be lower or equal to jpk-1.
	704	!! isolated ocean grid points are suppressed from mbathy
	705	!! since they are only connected to remaining
	706	!! ocean through vertical diffusion.
	707	!! C A U T I O N : mbathy will be modified during the initializa-
	708	!! tion phase to become the number of non-zero w-levels of a water
	709	!! column, with a minimum value of 1.
	710	!!
	711	!! ** Action : - update mbathy: level bathymetry (in level index)
	712	!! - update bathy : meter bathymetry (in meters)
	713	!!----------------------------------------------------------------------
	714	INTEGER :: ji, jj, jl ! dummy loop indices
	715	INTEGER :: icompt, ibtest, ikmax ! temporary integers
[11129]	716	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zbathy
[6951]	717	!!----------------------------------------------------------------------
	718	!
[11129]	719	ALLOCATE(zbathy(jpi,jpj))
[6951]	720	!
	721	IF(lwp) WRITE(numout,*)
	722	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
	723	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
	724	! ! Suppress isolated ocean grid points
	725	IF(lwp) WRITE(numout,*)
	726	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
	727	IF(lwp) WRITE(numout,*)' -----------------------------------'
	728	icompt = 0
	729	DO jl = 1, 2
[11129]	730	IF( l_Iperio ) THEN
[6951]	731	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
	732	mbathy(jpi,:) = mbathy( 2 ,:)
	733	ENDIF
[11129]	734	zbathy(:,:) = FLOAT( mbathy(:,:) )
	735	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
	736	mbathy(:,:) = INT( zbathy(:,:) )
	737
[6951]	738	DO jj = 2, jpjm1
	739	DO ji = 2, jpim1
	740	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
	741	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
	742	IF( ibtest < mbathy(ji,jj) ) THEN
	743	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
	744	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
	745	mbathy(ji,jj) = ibtest
	746	icompt = icompt + 1
	747	ENDIF
	748	END DO
	749	END DO
	750	END DO
[11129]	751
	752	IF( lk_mpp ) CALL mpp_sum( 'domzgr', icompt )
[6951]	753	IF( icompt == 0 ) THEN
	754	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
	755	ELSE
	756	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
	757	ENDIF
[11129]	758
	759	zbathy(:,:) = FLOAT( mbathy(:,:) )
	760	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
	761	mbathy(:,:) = INT( zbathy(:,:) )
	762
[6951]	763	! ! East-west cyclic boundary conditions
[11129]	764	IF( jperio == 0 ) THEN
	765	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: jperio = ', jperio
[6951]	766	IF( lk_mpp ) THEN
	767	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
	768	IF( jperio /= 1 ) mbathy(1,:) = 0
	769	ENDIF
	770	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
	771	IF( jperio /= 1 ) mbathy(nlci,:) = 0
	772	ENDIF
	773	ELSE
	774	IF( ln_zco .OR. ln_zps ) THEN
	775	mbathy( 1 ,:) = 0
	776	mbathy(jpi,:) = 0
	777	ELSE
	778	mbathy( 1 ,:) = jpkm1
	779	mbathy(jpi,:) = jpkm1
	780	ENDIF
	781	ENDIF
[11129]	782	ELSEIF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
	783	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio
[6951]	784	mbathy( 1 ,:) = mbathy(jpim1,:)
	785	mbathy(jpi,:) = mbathy( 2 ,:)
[11129]	786	ELSEIF( jperio == 2 ) THEN
	787	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: jperio = ', jperio
[6951]	788	ELSE
	789	IF(lwp) WRITE(numout,*) ' e r r o r'
[11129]	790	IF(lwp) WRITE(numout,*) ' parameter , jperio = ', jperio
[6951]	791	! STOP 'dom_mba'
	792	ENDIF
[11129]	793
[6951]	794	! Boundary condition on mbathy
	795	IF( .NOT.lk_mpp ) THEN
	796	!!gm !!bug ??? think about it !
	797	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
	798	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	799	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
[6951]	800	mbathy(:,:) = INT( zbathy(:,:) )
	801	ENDIF
[11129]	802
[6951]	803	! Number of ocean level inferior or equal to jpkm1
[11129]	804	zbathy(:,:) = FLOAT( mbathy(:,:) )
	805	ikmax = glob_max( 'domzgr', zbathy(:,:) )
	806
[6951]	807	IF( ikmax > jpkm1 ) THEN
	808	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
	809	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
	810	ELSE IF( ikmax < jpkm1 ) THEN
	811	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
	812	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
	813	ENDIF
	814	!
[11129]	815	DEALLOCATE( zbathy )
[6951]	816	!
	817	END SUBROUTINE zgr_bat_ctl
	818
	819
	820	SUBROUTINE zgr_bot_level
	821	!!----------------------------------------------------------------------
	822	!! * ROUTINE zgr_bot_level *
	823	!!
	824	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
	825	!!
	826	!! ** Method : computes from mbathy with a minimum value of 1 over land
	827	!!
	828	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
	829	!! ocean level at t-, u- & v-points
	830	!! (min value = 1 over land)
	831	!!----------------------------------------------------------------------
	832	INTEGER :: ji, jj ! dummy loop indices
[11129]	833	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmbk
[6951]	834	!!----------------------------------------------------------------------
	835	!
[11129]	836	ALLOCATE( zmbk(jpi,jpj) )
[6951]	837	!
	838	IF(lwp) WRITE(numout,*)
	839	IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
	840	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
	841	!
	842	mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
	843
	844	! ! bottom k-index of W-level = mbkt+1
	845	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
	846	DO ji = 1, jpim1
	847	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
	848	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
	849	END DO
	850	END DO
	851	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
[11129]	852	zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
	853	zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
[6951]	854	!
[11129]	855	DEALLOCATE( zmbk )
[6951]	856	!
	857	END SUBROUTINE zgr_bot_level
	858
	859
	860	SUBROUTINE zgr_top_level
	861	!!----------------------------------------------------------------------
	862	!! * ROUTINE zgr_top_level *
	863	!!
	864	!! ** Purpose : defines the vertical index of ocean top (mik. arrays)
	865	!!
	866	!! ** Method : computes from misfdep with a minimum value of 1
	867	!!
	868	!! ** Action : mikt, miku, mikv : vertical indices of the shallowest
	869	!! ocean level at t-, u- & v-points
	870	!! (min value = 1)
	871	!!----------------------------------------------------------------------
	872	INTEGER :: ji, jj ! dummy loop indices
[11129]	873	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmik
[6951]	874	!!----------------------------------------------------------------------
	875	!
[11129]	876	ALLOCATE( zmik(jpi,jpj) )
[6951]	877	!
	878	IF(lwp) WRITE(numout,*)
	879	IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
	880	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
	881	!
	882	mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1)
	883	! ! top k-index of W-level (=mikt)
	884	DO jj = 1, jpjm1 ! top k-index of U- (U-) level
	885	DO ji = 1, jpim1
	886	miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) )
	887	mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) )
	888	mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) )
	889	END DO
	890	END DO
	891
	892	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
[11129]	893	zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 )
	894	zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 )
	895	zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 )
[6951]	896	!
[11129]	897	DEALLOCATE( zmik )
[6951]	898	!
	899	END SUBROUTINE zgr_top_level
	900
	901
	902	SUBROUTINE zgr_zco
	903	!!----------------------------------------------------------------------
	904	!! * ROUTINE zgr_zco *
	905	!!
	906	!! ** Purpose : define the reference z-coordinate system
	907	!!
	908	!! ** Method : set 3D coord. arrays to reference 1D array
	909	!!----------------------------------------------------------------------
	910	INTEGER :: jk
	911	!!----------------------------------------------------------------------
	912	!
	913	DO jk = 1, jpk
	914	gdept_0(:,:,jk) = gdept_1d(jk)
	915	gdepw_0(:,:,jk) = gdepw_1d(jk)
	916	gde3w_0(:,:,jk) = gdepw_1d(jk)
	917	e3t_0 (:,:,jk) = e3t_1d (jk)
	918	e3u_0 (:,:,jk) = e3t_1d (jk)
	919	e3v_0 (:,:,jk) = e3t_1d (jk)
	920	e3f_0 (:,:,jk) = e3t_1d (jk)
	921	e3w_0 (:,:,jk) = e3w_1d (jk)
	922	e3uw_0 (:,:,jk) = e3w_1d (jk)
	923	e3vw_0 (:,:,jk) = e3w_1d (jk)
	924	END DO
	925	!
	926	END SUBROUTINE zgr_zco
	927
	928
	929	SUBROUTINE zgr_zps
	930	!!----------------------------------------------------------------------
	931	!! * ROUTINE zgr_zps *
	932	!!
	933	!! ** Purpose : the depth and vertical scale factor in partial step
	934	!! reference z-coordinate case
	935	!!
	936	!! ** Method : Partial steps : computes the 3D vertical scale factors
	937	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
	938	!! a partial step representation of bottom topography.
	939	!!
	940	!! The reference depth of model levels is defined from an analytical
	941	!! function the derivative of which gives the reference vertical
	942	!! scale factors.
	943	!! From depth and scale factors reference, we compute there new value
	944	!! with partial steps on 3d arrays ( i, j, k ).
	945	!!
	946	!! w-level: gdepw_0(i,j,k) = gdep(k)
	947	!! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k)
	948	!! t-level: gdept_0(i,j,k) = gdep(k+0.5)
	949	!! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
	950	!!
	951	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
	952	!! we find the mbathy index of the depth at each grid point.
	953	!! This leads us to three cases:
	954	!!
	955	!! - bathy = 0 => mbathy = 0
	956	!! - 1 < mbathy < jpkm1
	957	!! - bathy > gdepw_0(jpk) => mbathy = jpkm1
	958	!!
	959	!! Then, for each case, we find the new depth at t- and w- levels
	960	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
	961	!! and f-points.
	962	!!
	963	!! This routine is given as an example, it must be modified
	964	!! following the user s desiderata. nevertheless, the output as
	965	!! well as the way to compute the model levels and scale factors
	966	!! must be respected in order to insure second order accuracy
	967	!! schemes.
	968	!!
	969	!! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
	970	!! - - - - - - - gdept_0, gdepw_0 and e3. are positives
	971	!!
	972	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
	973	!!----------------------------------------------------------------------
	974	INTEGER :: ji, jj, jk ! dummy loop indices
	975	INTEGER :: ik, it, ikb, ikt ! temporary integers
	976	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
	977	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
	978	REAL(wp) :: zdiff ! temporary scalar
	979	REAL(wp) :: zmax ! temporary scalar
[11129]	980	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zprt
[6951]	981	!!---------------------------------------------------------------------
	982	!
[11129]	983	ALLOCATE( zprt(jpi,jpj,jpk) )
[6951]	984	!
	985	IF(lwp) WRITE(numout,*)
	986	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
	987	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
	988	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
	989
	990	! bathymetry in level (from bathy_meter)
	991	! ===================
	992	zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
	993	bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
	994	WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
	995	ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
	996	END WHERE
	997
	998	! Compute mbathy for ocean points (i.e. the number of ocean levels)
	999	! find the number of ocean levels such that the last level thickness
	1000	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
	1001	! e3t_1d is the reference level thickness
	1002	DO jk = jpkm1, 1, -1
	1003	zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
	1004	WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
	1005	END DO
	1006
	1007	! Scale factors and depth at T- and W-points
	1008	DO jk = 1, jpk ! intitialization to the reference z-coordinate
	1009	gdept_0(:,:,jk) = gdept_1d(jk)
	1010	gdepw_0(:,:,jk) = gdepw_1d(jk)
	1011	e3t_0 (:,:,jk) = e3t_1d (jk)
	1012	e3w_0 (:,:,jk) = e3w_1d (jk)
	1013	END DO
	1014
	1015	! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf
	1016	IF ( ln_isfcav ) CALL zgr_isf
	1017
	1018	! Scale factors and depth at T- and W-points
	1019	IF ( .NOT. ln_isfcav ) THEN
	1020	DO jj = 1, jpj
	1021	DO ji = 1, jpi
	1022	ik = mbathy(ji,jj)
	1023	IF( ik > 0 ) THEN ! ocean point only
	1024	! max ocean level case
	1025	IF( ik == jpkm1 ) THEN
	1026	zdepwp = bathy(ji,jj)
	1027	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
	1028	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
	1029	e3t_0(ji,jj,ik ) = ze3tp
	1030	e3t_0(ji,jj,ik+1) = ze3tp
	1031	e3w_0(ji,jj,ik ) = ze3wp
	1032	e3w_0(ji,jj,ik+1) = ze3tp
	1033	gdepw_0(ji,jj,ik+1) = zdepwp
	1034	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
	1035	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
	1036	!
	1037	ELSE ! standard case
	1038	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
	1039	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
	1040	ENDIF
	1041	!gm Bug? check the gdepw_1d
	1042	! ... on ik
	1043	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
	1044	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
	1045	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
	1046	e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
	1047	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
	1048	e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
	1049	& * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
	1050	! ... on ik+1
	1051	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
	1052	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
	1053	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
	1054	ENDIF
	1055	ENDIF
	1056	END DO
	1057	END DO
	1058	!
	1059	it = 0
	1060	DO jj = 1, jpj
	1061	DO ji = 1, jpi
	1062	ik = mbathy(ji,jj)
	1063	IF( ik > 0 ) THEN ! ocean point only
	1064	e3tp (ji,jj) = e3t_0(ji,jj,ik)
	1065	e3wp (ji,jj) = e3w_0(ji,jj,ik)
	1066	! test
	1067	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
	1068	IF( zdiff <= 0._wp .AND. lwp ) THEN
	1069	it = it + 1
	1070	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
	1071	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
	1072	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
	1073	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
	1074	ENDIF
	1075	ENDIF
	1076	END DO
	1077	END DO
	1078	END IF
	1079	!
	1080	! Scale factors and depth at U-, V-, UW and VW-points
	1081	DO jk = 1, jpk ! initialisation to z-scale factors
	1082	e3u_0 (:,:,jk) = e3t_1d(jk)
	1083	e3v_0 (:,:,jk) = e3t_1d(jk)
	1084	e3uw_0(:,:,jk) = e3w_1d(jk)
	1085	e3vw_0(:,:,jk) = e3w_1d(jk)
	1086	END DO
	1087
	1088	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
	1089	DO jj = 1, jpjm1
	1090	DO ji = 1, jpim1 ! vector opt.
	1091	e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
	1092	e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
	1093	e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
	1094	e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
	1095	END DO
	1096	END DO
	1097	END DO
	1098	IF ( ln_isfcav ) THEN
	1099	! (ISF) define e3uw (adapted for 2 cells in the water column)
	1100	DO jj = 2, jpjm1
	1101	DO ji = 2, jpim1 ! vector opt.
	1102	ikb = MAX(mbathy (ji,jj),mbathy (ji+1,jj))
	1103	ikt = MAX(misfdep(ji,jj),misfdep(ji+1,jj))
	1104	IF (ikb == ikt+1) e3uw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji+1,jj ,ikb ) ) &
	1105	& - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji+1,jj ,ikb-1) )
	1106	ikb = MAX(mbathy (ji,jj),mbathy (ji,jj+1))
	1107	ikt = MAX(misfdep(ji,jj),misfdep(ji,jj+1))
	1108	IF (ikb == ikt+1) e3vw_0(ji,jj,ikb) = MIN( gdept_0(ji,jj,ikb ), gdept_0(ji ,jj+1,ikb ) ) &
	1109	& - MAX( gdept_0(ji,jj,ikb-1), gdept_0(ji ,jj+1,ikb-1) )
	1110	END DO
	1111	END DO
	1112	END IF
	1113
[11129]	1114	CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp ) ! lateral boundary conditions
	1115	CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
[6951]	1116	!
	1117
	1118	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
	1119	WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)
	1120	WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)
	1121	WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)
	1122	WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)
	1123	END DO
	1124
	1125	! Scale factor at F-point
	1126	DO jk = 1, jpk ! initialisation to z-scale factors
	1127	e3f_0(:,:,jk) = e3t_1d(jk)
	1128	END DO
	1129	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
	1130	DO jj = 1, jpjm1
	1131	DO ji = 1, jpim1 ! vector opt.
	1132	e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
	1133	END DO
	1134	END DO
	1135	END DO
[11129]	1136	CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp ) ! Lateral boundary conditions
[6951]	1137	!
	1138	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
	1139	WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)
	1140	END DO
	1141	!!gm bug ? : must be a do loop with mj0,mj1
	1142	!
	1143	e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
	1144	e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)
	1145	e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)
	1146	e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)
	1147	e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)
	1148
	1149	! Control of the sign
	1150	IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )
	1151	IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )
	1152	IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )
	1153	IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )
	1154
	1155	! Compute gde3w_0 (vertical sum of e3w)
	1156	IF ( ln_isfcav ) THEN ! if cavity
	1157	WHERE( misfdep == 0 ) misfdep = 1
	1158	DO jj = 1,jpj
	1159	DO ji = 1,jpi
	1160	gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)
	1161	DO jk = 2, misfdep(ji,jj)
	1162	gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
	1163	END DO
	1164	IF( misfdep(ji,jj) >= 2 ) gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))
	1165	DO jk = misfdep(ji,jj) + 1, jpk
	1166	gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)
	1167	END DO
	1168	END DO
	1169	END DO
	1170	ELSE ! no cavity
	1171	gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
	1172	DO jk = 2, jpk
	1173	gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
	1174	END DO
	1175	END IF
	1176	!
[11129]	1177	DEALLOCATE( zprt )
[6951]	1178	!
	1179	END SUBROUTINE zgr_zps
	1180
	1181
	1182	SUBROUTINE zgr_isf
	1183	!!----------------------------------------------------------------------
	1184	!! * ROUTINE zgr_isf *
	1185	!!
	1186	!! ** Purpose : check the bathymetry in levels
	1187	!!
	1188	!! ** Method : THe water column have to contained at least 2 cells
	1189	!! Bathymetry and isfdraft are modified (dig/close) to respect
	1190	!! this criterion.
	1191	!!
	1192	!! ** Action : - test compatibility between isfdraft and bathy
	1193	!! - bathy and isfdraft are modified
	1194	!!----------------------------------------------------------------------
	1195	INTEGER :: ji, jj, jl, jk ! dummy loop indices
	1196	INTEGER :: ik, it ! temporary integers
	1197	INTEGER :: icompt, ibtest ! (ISF)
	1198	INTEGER :: ibtestim1, ibtestip1 ! (ISF)
	1199	INTEGER :: ibtestjm1, ibtestjp1 ! (ISF)
	1200	REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t
	1201	REAL(wp) :: zmax ! Maximum and minimum depth
	1202	REAL(wp) :: zbathydiff ! isf temporary scalar
	1203	REAL(wp) :: zrisfdepdiff ! isf temporary scalar
	1204	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
	1205	REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t
	1206	REAL(wp) :: zdiff ! temporary scalar
[11129]	1207	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH)
	1208	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH)
[6951]	1209	!!---------------------------------------------------------------------
	1210	!
[11129]	1211	ALLOCATE( zbathy(jpi,jpj), zmask(jpi,jpj), zrisfdep(jpi,jpj) )
	1212	ALLOCATE( zmisfdep(jpi,jpj), zmbathy(jpi,jpj) )
[6951]	1213	!
	1214	! (ISF) compute misfdep
	1215	WHERE( risfdep(:,:) == 0._wp .AND. bathy(:,:) /= 0 ) ; misfdep(:,:) = 1 ! open water : set misfdep to 1
	1216	ELSEWHERE ; misfdep(:,:) = 2 ! iceshelf : initialize misfdep to second level
	1217	END WHERE
	1218
	1219	! Compute misfdep for ocean points (i.e. first wet level)
	1220	! find the first ocean level such that the first level thickness
	1221	! is larger than the bot_level of e3zps_min and e3zps_rat * e3t_0 (where
	1222	! e3t_0 is the reference level thickness
	1223	DO jk = 2, jpkm1
	1224	zdepth = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
	1225	WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1
	1226	END DO
	1227	WHERE ( 0._wp < risfdep(:,:) .AND. risfdep(:,:) <= e3t_1d(1) )
	1228	risfdep(:,:) = 0. ; misfdep(:,:) = 1
	1229	END WHERE
	1230
	1231	! remove very shallow ice shelf (less than ~ 10m if 75L)
	1232	WHERE (risfdep(:,:) <= 10._wp .AND. misfdep(:,:) > 1)
	1233	misfdep = 0; risfdep = 0.0_wp;
	1234	mbathy = 0; bathy = 0.0_wp;
	1235	END WHERE
	1236	WHERE (bathy(:,:) <= 30.0_wp .AND. gphit < -60._wp)
	1237	misfdep = 0; risfdep = 0.0_wp;
	1238	mbathy = 0; bathy = 0.0_wp;
	1239	END WHERE
	1240
	1241	! basic check for the compatibility of bathy and risfdep. I think it should be offline because it is not perfect and cannot solved all the situation
	1242	icompt = 0
	1243	! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together
	1244	DO jl = 1, 10
	1245	! check at each iteration if isf is grounded or not (1cm treshold have to be update after first coupling experiments)
	1246	WHERE (bathy(:,:) <= risfdep(:,:) + rn_isfhmin)
	1247	misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp
	1248	mbathy (:,:) = 0 ; bathy (:,:) = 0._wp
	1249	END WHERE
	1250	WHERE (mbathy(:,:) <= 0)
	1251	misfdep(:,:) = 0; risfdep(:,:) = 0._wp
	1252	mbathy (:,:) = 0; bathy (:,:) = 0._wp
	1253	END WHERE
	1254	IF( lk_mpp ) THEN
	1255	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1256	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1257	misfdep(:,:) = INT( zbathy(:,:) )
	1258
[11129]	1259	CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
	1260	CALL lbc_lnk( 'domzgr',bathy, 'T', 1. )
[6951]	1261
	1262	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1263	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1264	mbathy(:,:) = INT( zbathy(:,:) )
	1265	ENDIF
	1266	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
	1267	misfdep( 1 ,:) = misfdep(jpim1,:) ! local domain is cyclic east-west
	1268	misfdep(jpi,:) = misfdep( 2 ,:)
	1269	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
	1270	mbathy(jpi,:) = mbathy( 2 ,:)
	1271	ENDIF
	1272
	1273	! split last cell if possible (only where water column is 2 cell or less)
	1274	! if coupled to ice sheet, we do not modify the bathymetry (can be discuss).
	1275	IF ( .NOT. ln_iscpl) THEN
	1276	DO jk = jpkm1, 1, -1
	1277	zmax = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
	1278	WHERE( gdepw_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= zmax .AND. misfdep + 1 >= mbathy)
	1279	mbathy(:,:) = jk
	1280	bathy(:,:) = zmax
	1281	END WHERE
	1282	END DO
	1283	END IF
	1284
	1285	! split top cell if possible (only where water column is 2 cell or less)
	1286	DO jk = 2, jpkm1
	1287	zmax = gdepw_1d(jk+1) - MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
	1288	WHERE( gdepw_1d(jk+1) > risfdep(:,:) .AND. risfdep(:,:) >= zmax .AND. misfdep + 1 >= mbathy)
	1289	misfdep(:,:) = jk
	1290	risfdep(:,:) = zmax
	1291	END WHERE
	1292	END DO
	1293
	1294
	1295	! Case where bathy and risfdep compatible but not the level variable mbathy/misfdep because of partial cell condition
	1296	DO jj = 1, jpj
	1297	DO ji = 1, jpi
	1298	! find the minimum change option:
	1299	! test bathy
	1300	IF (risfdep(ji,jj) > 1) THEN
	1301	IF ( .NOT. ln_iscpl ) THEN
	1302	zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) &
	1303	& + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
	1304	zrisfdepdiff=ABS(risfdep(ji,jj) - (gdepw_1d(misfdep(ji,jj) ) &
	1305	& - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
	1306	IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN
	1307	IF (zbathydiff <= zrisfdepdiff) THEN
	1308	bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat )
	1309	mbathy(ji,jj)= mbathy(ji,jj) + 1
	1310	ELSE
	1311	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
	1312	misfdep(ji,jj) = misfdep(ji,jj) - 1
	1313	END IF
	1314	ENDIF
	1315	ELSE
	1316	IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN
	1317	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )
	1318	misfdep(ji,jj) = misfdep(ji,jj) - 1
	1319	END IF
	1320	END IF
	1321	END IF
	1322	END DO
	1323	END DO
	1324
	1325	! At least 2 levels for water thickness at T, U, and V point.
	1326	DO jj = 1, jpj
	1327	DO ji = 1, jpi
	1328	! find the minimum change option:
	1329	! test bathy
	1330	IF( misfdep(ji,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
	1331	IF ( .NOT. ln_iscpl ) THEN
	1332	zbathydiff =ABS(bathy(ji,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &
	1333	& + MIN( e3zps_min,e3t_1d(mbathy (ji,jj)+1)*e3zps_rat )))
	1334	zrisfdepdiff=ABS(risfdep(ji,jj) - ( gdepw_1d(misfdep(ji,jj) ) &
	1335	& - MIN( e3zps_min,e3t_1d(misfdep(ji,jj)-1)*e3zps_rat )))
	1336	IF (zbathydiff <= zrisfdepdiff) THEN
	1337	mbathy(ji,jj) = mbathy(ji,jj) + 1
	1338	bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
	1339	ELSE
	1340	misfdep(ji,jj)= misfdep(ji,jj) - 1
	1341	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
	1342	END IF
	1343	ELSE
	1344	misfdep(ji,jj)= misfdep(ji,jj) - 1
	1345	risfdep(ji,jj)= gdepw_1d(misfdep(ji,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj))*e3zps_rat )
	1346	END IF
	1347	ENDIF
	1348	END DO
	1349	END DO
	1350
	1351	! point V mbathy(ji,jj) == misfdep(ji,jj+1)
	1352	DO jj = 1, jpjm1
	1353	DO ji = 1, jpim1
	1354	IF( misfdep(ji,jj+1) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
	1355	IF ( .NOT. ln_iscpl ) THEN
	1356	zbathydiff =ABS(bathy(ji,jj ) - ( gdepw_1d(mbathy (ji,jj)+1) &
	1357	& + MIN( e3zps_min, e3t_1d(mbathy (ji,jj )+1)*e3zps_rat )))
	1358	zrisfdepdiff=ABS(risfdep(ji,jj+1) - ( gdepw_1d(misfdep(ji,jj+1)) &
	1359	& - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1)-1)*e3zps_rat )))
	1360	IF (zbathydiff <= zrisfdepdiff) THEN
	1361	mbathy(ji,jj) = mbathy(ji,jj) + 1
	1362	bathy(ji,jj) = gdepw_1d(mbathy (ji,jj )) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj )+1)*e3zps_rat )
	1363	ELSE
	1364	misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1
	1365	risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
	1366	END IF
	1367	ELSE
	1368	misfdep(ji,jj+1) = misfdep(ji,jj+1) - 1
	1369	risfdep (ji,jj+1) = gdepw_1d(misfdep(ji,jj+1)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj+1))*e3zps_rat )
	1370	END IF
	1371	ENDIF
	1372	END DO
	1373	END DO
	1374
	1375	IF( lk_mpp ) THEN
	1376	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1377	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1378	misfdep(:,:) = INT( zbathy(:,:) )
	1379
[11129]	1380	CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
	1381	CALL lbc_lnk( 'domzgr',bathy, 'T', 1. )
[6951]	1382
	1383	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1384	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1385	mbathy(:,:) = INT( zbathy(:,:) )
	1386	ENDIF
	1387	! point V misdep(ji,jj) == mbathy(ji,jj+1)
	1388	DO jj = 1, jpjm1
	1389	DO ji = 1, jpim1
	1390	IF( misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) > 1) THEN
	1391	IF ( .NOT. ln_iscpl ) THEN
	1392	zbathydiff =ABS( bathy(ji,jj+1) - ( gdepw_1d(mbathy (ji,jj+1)+1) &
	1393	& + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )))
	1394	zrisfdepdiff=ABS(risfdep(ji,jj ) - ( gdepw_1d(misfdep(ji,jj ) ) &
	1395	& - MIN( e3zps_min, e3t_1d(misfdep(ji,jj )-1)*e3zps_rat )))
	1396	IF (zbathydiff <= zrisfdepdiff) THEN
	1397	mbathy (ji,jj+1) = mbathy(ji,jj+1) + 1
	1398	bathy (ji,jj+1) = gdepw_1d(mbathy (ji,jj+1) ) + MIN( e3zps_min, e3t_1d(mbathy (ji,jj+1)+1)*e3zps_rat )
	1399	ELSE
	1400	misfdep(ji,jj) = misfdep(ji,jj) - 1
	1401	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )
	1402	END IF
	1403	ELSE
	1404	misfdep(ji,jj) = misfdep(ji,jj) - 1
	1405	risfdep(ji,jj) = gdepw_1d(misfdep(ji,jj )+1) - MIN( e3zps_min, e3t_1d(misfdep(ji,jj ) )*e3zps_rat )
	1406	END IF
	1407	ENDIF
	1408	END DO
	1409	END DO
	1410
	1411
	1412	IF( lk_mpp ) THEN
	1413	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1414	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1415	misfdep(:,:) = INT( zbathy(:,:) )
	1416
[11129]	1417	CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
	1418	CALL lbc_lnk( 'domzgr',bathy, 'T', 1. )
[6951]	1419
	1420	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1421	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1422	mbathy(:,:) = INT( zbathy(:,:) )
	1423	ENDIF
	1424
	1425	! point U mbathy(ji,jj) == misfdep(ji,jj+1)
	1426	DO jj = 1, jpjm1
	1427	DO ji = 1, jpim1
	1428	IF( misfdep(ji+1,jj) == mbathy(ji,jj) .AND. mbathy(ji,jj) > 1) THEN
	1429	IF ( .NOT. ln_iscpl ) THEN
	1430	zbathydiff =ABS( bathy(ji ,jj) - ( gdepw_1d(mbathy (ji,jj)+1) &
	1431	& + MIN( e3zps_min, e3t_1d(mbathy (ji ,jj)+1)*e3zps_rat )))
	1432	zrisfdepdiff=ABS(risfdep(ji+1,jj) - ( gdepw_1d(misfdep(ji+1,jj)) &
	1433	& - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj)-1)*e3zps_rat )))
	1434	IF (zbathydiff <= zrisfdepdiff) THEN
	1435	mbathy(ji,jj) = mbathy(ji,jj) + 1
	1436	bathy(ji,jj) = gdepw_1d(mbathy (ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj) +1)*e3zps_rat )
	1437	ELSE
	1438	misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
	1439	risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
	1440	END IF
	1441	ELSE
	1442	misfdep(ji+1,jj)= misfdep(ji+1,jj) - 1
	1443	risfdep(ji+1,jj) = gdepw_1d(misfdep(ji+1,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji+1,jj))*e3zps_rat )
	1444	ENDIF
	1445	ENDIF
	1446	ENDDO
	1447	ENDDO
	1448
	1449	IF( lk_mpp ) THEN
	1450	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1451	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1452	misfdep(:,:) = INT( zbathy(:,:) )
	1453
[11129]	1454	CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
	1455	CALL lbc_lnk( 'domzgr',bathy, 'T', 1. )
[6951]	1456
	1457	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1458	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1459	mbathy(:,:) = INT( zbathy(:,:) )
	1460	ENDIF
	1461
	1462	! point U misfdep(ji,jj) == bathy(ji,jj+1)
	1463	DO jj = 1, jpjm1
	1464	DO ji = 1, jpim1
	1465	IF( misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) > 1) THEN
	1466	IF ( .NOT. ln_iscpl ) THEN
	1467	zbathydiff =ABS( bathy(ji+1,jj) - ( gdepw_1d(mbathy (ji+1,jj)+1) &
	1468	& + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj)+1)*e3zps_rat )))
	1469	zrisfdepdiff=ABS(risfdep(ji ,jj) - ( gdepw_1d(misfdep(ji ,jj) ) &
	1470	& - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj)-1)*e3zps_rat )))
	1471	IF (zbathydiff <= zrisfdepdiff) THEN
	1472	mbathy(ji+1,jj) = mbathy (ji+1,jj) + 1
	1473	bathy (ji+1,jj) = gdepw_1d(mbathy (ji+1,jj) ) + MIN( e3zps_min, e3t_1d(mbathy (ji+1,jj) +1)*e3zps_rat )
	1474	ELSE
	1475	misfdep(ji,jj) = misfdep(ji ,jj) - 1
	1476	risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )
	1477	END IF
	1478	ELSE
	1479	misfdep(ji,jj) = misfdep(ji ,jj) - 1
	1480	risfdep(ji,jj) = gdepw_1d(misfdep(ji ,jj)+1) - MIN( e3zps_min, e3t_1d(misfdep(ji ,jj) )*e3zps_rat )
	1481	ENDIF
	1482	ENDIF
	1483	ENDDO
	1484	ENDDO
	1485
	1486	IF( lk_mpp ) THEN
	1487	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1488	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1489	misfdep(:,:) = INT( zbathy(:,:) )
	1490
[11129]	1491	CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
	1492	CALL lbc_lnk( 'domzgr',bathy, 'T', 1. )
[6951]	1493
	1494	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1495	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1496	mbathy(:,:) = INT( zbathy(:,:) )
	1497	ENDIF
	1498	END DO
	1499	! end dig bathy/ice shelf to be compatible
	1500	! now fill single point in "coastline" of ice shelf, bathy, hole, and test again one cell tickness
	1501	DO jl = 1,20
	1502
	1503	! remove single point "bay" on isf coast line in the ice shelf draft'
	1504	DO jk = 2, jpk
	1505	WHERE (misfdep==0) misfdep=jpk
	1506	zmask=0._wp
	1507	WHERE (misfdep <= jk) zmask=1
	1508	DO jj = 2, jpjm1
	1509	DO ji = 2, jpim1
	1510	IF (misfdep(ji,jj) == jk) THEN
	1511	ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
	1512	IF (ibtest <= 1) THEN
	1513	risfdep(ji,jj)=gdepw_1d(jk+1) ; misfdep(ji,jj)=jk+1
	1514	IF (misfdep(ji,jj) > mbathy(ji,jj)) misfdep(ji,jj) = jpk
	1515	END IF
	1516	END IF
	1517	END DO
	1518	END DO
	1519	END DO
	1520	WHERE (misfdep==jpk)
	1521	misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
	1522	END WHERE
	1523	IF( lk_mpp ) THEN
	1524	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1525	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1526	misfdep(:,:) = INT( zbathy(:,:) )
	1527
[11129]	1528	CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
	1529	CALL lbc_lnk('domzgr', bathy, 'T', 1. )
[6951]	1530
	1531	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1532	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1533	mbathy(:,:) = INT( zbathy(:,:) )
	1534	ENDIF
	1535
	1536	! remove single point "bay" on bathy coast line beneath an ice shelf'
	1537	DO jk = jpk,1,-1
	1538	zmask=0._wp
	1539	WHERE (mbathy >= jk ) zmask=1
	1540	DO jj = 2, jpjm1
	1541	DO ji = 2, jpim1
	1542	IF (mbathy(ji,jj) == jk .AND. misfdep(ji,jj) >= 2) THEN
	1543	ibtest = zmask(ji-1,jj) + zmask(ji+1,jj) + zmask(ji,jj-1) + zmask(ji,jj+1)
	1544	IF (ibtest <= 1) THEN
	1545	bathy(ji,jj)=gdepw_1d(jk) ; mbathy(ji,jj)=jk-1
	1546	IF (misfdep(ji,jj) > mbathy(ji,jj)) mbathy(ji,jj) = 0
	1547	END IF
	1548	END IF
	1549	END DO
	1550	END DO
	1551	END DO
	1552	WHERE (mbathy==0)
	1553	misfdep=0 ; risfdep=0._wp ; mbathy=0 ; bathy=0._wp
	1554	END WHERE
	1555	IF( lk_mpp ) THEN
	1556	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1557	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1558	misfdep(:,:) = INT( zbathy(:,:) )
	1559
[11129]	1560	CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
	1561	CALL lbc_lnk( 'domzgr',bathy, 'T', 1. )
[6951]	1562
	1563	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1564	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1565	mbathy(:,:) = INT( zbathy(:,:) )
	1566	ENDIF
	1567
	1568	! fill hole in ice shelf
	1569	zmisfdep = misfdep
	1570	zrisfdep = risfdep
	1571	WHERE (zmisfdep <= 1._wp) zmisfdep=jpk
	1572	DO jj = 2, jpjm1
	1573	DO ji = 2, jpim1
	1574	ibtestim1 = zmisfdep(ji-1,jj ) ; ibtestip1 = zmisfdep(ji+1,jj )
	1575	ibtestjm1 = zmisfdep(ji ,jj-1) ; ibtestjp1 = zmisfdep(ji ,jj+1)
	1576	IF( zmisfdep(ji,jj) >= mbathy(ji-1,jj ) ) ibtestim1 = jpk
	1577	IF( zmisfdep(ji,jj) >= mbathy(ji+1,jj ) ) ibtestip1 = jpk
	1578	IF( zmisfdep(ji,jj) >= mbathy(ji ,jj-1) ) ibtestjm1 = jpk
	1579	IF( zmisfdep(ji,jj) >= mbathy(ji ,jj+1) ) ibtestjp1 = jpk
	1580	ibtest=MIN(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
	1581	IF( ibtest == jpk .AND. misfdep(ji,jj) >= 2) THEN
	1582	mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp
	1583	END IF
	1584	IF( zmisfdep(ji,jj) < ibtest .AND. misfdep(ji,jj) >= 2) THEN
	1585	misfdep(ji,jj) = ibtest
	1586	risfdep(ji,jj) = gdepw_1d(ibtest)
	1587	ENDIF
	1588	ENDDO
	1589	ENDDO
	1590
	1591	IF( lk_mpp ) THEN
	1592	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1593	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1594	misfdep(:,:) = INT( zbathy(:,:) )
	1595
[11129]	1596	CALL lbc_lnk( 'domzgr',risfdep, 'T', 1. )
	1597	CALL lbc_lnk( 'domzgr',bathy, 'T', 1. )
[6951]	1598
	1599	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1600	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1601	mbathy(:,:) = INT( zbathy(:,:) )
	1602	ENDIF
	1603	!
	1604	!! fill hole in bathymetry
	1605	zmbathy (:,:)=mbathy (:,:)
	1606	DO jj = 2, jpjm1
	1607	DO ji = 2, jpim1
	1608	ibtestim1 = zmbathy(ji-1,jj ) ; ibtestip1 = zmbathy(ji+1,jj )
	1609	ibtestjm1 = zmbathy(ji ,jj-1) ; ibtestjp1 = zmbathy(ji ,jj+1)
	1610	IF( zmbathy(ji,jj) < misfdep(ji-1,jj ) ) ibtestim1 = 0
	1611	IF( zmbathy(ji,jj) < misfdep(ji+1,jj ) ) ibtestip1 = 0
	1612	IF( zmbathy(ji,jj) < misfdep(ji ,jj-1) ) ibtestjm1 = 0
	1613	IF( zmbathy(ji,jj) < misfdep(ji ,jj+1) ) ibtestjp1 = 0
	1614	ibtest=MAX(ibtestim1, ibtestip1, ibtestjm1, ibtestjp1)
	1615	IF( ibtest == 0 .AND. misfdep(ji,jj) >= 2) THEN
	1616	mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0.0_wp ; misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0.0_wp ;
	1617	END IF
	1618	IF( ibtest < zmbathy(ji,jj) .AND. misfdep(ji,jj) >= 2) THEN
	1619	mbathy(ji,jj) = ibtest
	1620	bathy(ji,jj) = gdepw_1d(ibtest+1)
	1621	ENDIF
	1622	END DO
	1623	END DO
	1624	IF( lk_mpp ) THEN
	1625	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1626	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1627	misfdep(:,:) = INT( zbathy(:,:) )
	1628
[11129]	1629	CALL lbc_lnk( 'domzgr',risfdep, 'T', 1. )
	1630	CALL lbc_lnk( 'domzgr',bathy, 'T', 1. )
[6951]	1631
	1632	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1633	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1634	mbathy(:,:) = INT( zbathy(:,:) )
	1635	ENDIF
	1636	! if not compatible after all check (ie U point water column less than 2 cells), mask U
	1637	DO jj = 1, jpjm1
	1638	DO ji = 1, jpim1
	1639	IF (mbathy(ji,jj) == misfdep(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN
	1640	mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;
	1641	END IF
	1642	END DO
	1643	END DO
	1644	IF( lk_mpp ) THEN
	1645	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1646	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1647	misfdep(:,:) = INT( zbathy(:,:) )
	1648
[11129]	1649	CALL lbc_lnk('domzgr', risfdep, 'T', 1. )
	1650	CALL lbc_lnk('domzgr', bathy, 'T', 1. )
[6951]	1651
	1652	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1653	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1654	mbathy(:,:) = INT( zbathy(:,:) )
	1655	ENDIF
	1656	! if not compatible after all check (ie U point water column less than 2 cells), mask U
	1657	DO jj = 1, jpjm1
	1658	DO ji = 1, jpim1
	1659	IF (misfdep(ji,jj) == mbathy(ji+1,jj) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji+1,jj) >= 1) THEN
	1660	mbathy(ji+1,jj) = mbathy(ji+1,jj) - 1; bathy(ji+1,jj) = gdepw_1d(mbathy(ji+1,jj)+1) ;
	1661	END IF
	1662	END DO
	1663	END DO
	1664	IF( lk_mpp ) THEN
	1665	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1666	CALL lbc_lnk('domzgr', zbathy, 'T', 1. )
[6951]	1667	misfdep(:,:) = INT( zbathy(:,:) )
	1668
[11129]	1669	CALL lbc_lnk('domzgr', risfdep,'T', 1. )
	1670	CALL lbc_lnk( 'domzgr',bathy, 'T', 1. )
[6951]	1671
	1672	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1673	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1674	mbathy(:,:) = INT( zbathy(:,:) )
	1675	ENDIF
	1676	! if not compatible after all check (ie V point water column less than 2 cells), mask V
	1677	DO jj = 1, jpjm1
	1678	DO ji = 1, jpi
	1679	IF (mbathy(ji,jj) == misfdep(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN
	1680	mbathy(ji,jj) = mbathy(ji,jj) - 1 ; bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)+1) ;
	1681	END IF
	1682	END DO
	1683	END DO
	1684	IF( lk_mpp ) THEN
	1685	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1686	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1687	misfdep(:,:) = INT( zbathy(:,:) )
	1688
[11129]	1689	CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
	1690	CALL lbc_lnk('domzgr', bathy, 'T', 1. )
[6951]	1691
	1692	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1693	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1694	mbathy(:,:) = INT( zbathy(:,:) )
	1695	ENDIF
	1696	! if not compatible after all check (ie V point water column less than 2 cells), mask V
	1697	DO jj = 1, jpjm1
	1698	DO ji = 1, jpi
	1699	IF (misfdep(ji,jj) == mbathy(ji,jj+1) .AND. mbathy(ji,jj) >= 1 .AND. mbathy(ji,jj+1) >= 1) THEN
	1700	mbathy(ji,jj+1) = mbathy(ji,jj+1) - 1 ; bathy(ji,jj+1) = gdepw_1d(mbathy(ji,jj+1)+1) ;
	1701	END IF
	1702	END DO
	1703	END DO
	1704	IF( lk_mpp ) THEN
	1705	zbathy(:,:) = FLOAT( misfdep(:,:) )
[11129]	1706	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1707	misfdep(:,:) = INT( zbathy(:,:) )
	1708
[11129]	1709	CALL lbc_lnk( 'domzgr',risfdep,'T', 1. )
	1710	CALL lbc_lnk( 'domzgr',bathy, 'T', 1. )
[6951]	1711
	1712	zbathy(:,:) = FLOAT( mbathy(:,:) )
[11129]	1713	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1. )
[6951]	1714	mbathy(:,:) = INT( zbathy(:,:) )
	1715	ENDIF
	1716	! if not compatible after all check, mask T
	1717	DO jj = 1, jpj
	1718	DO ji = 1, jpi
	1719	IF (mbathy(ji,jj) <= misfdep(ji,jj)) THEN
	1720	misfdep(ji,jj) = 0 ; risfdep(ji,jj) = 0._wp ; mbathy(ji,jj) = 0 ; bathy(ji,jj) = 0._wp ;
	1721	END IF
	1722	END DO
	1723	END DO
	1724
	1725	WHERE (mbathy(:,:) == 1)
	1726	mbathy = 0; bathy = 0.0_wp ; misfdep = 0 ; risfdep = 0.0_wp
	1727	END WHERE
	1728	END DO
	1729	! end check compatibility ice shelf/bathy
	1730	! remove very shallow ice shelf (less than ~ 10m if 75L)
	1731	WHERE (risfdep(:,:) <= 10._wp)
	1732	misfdep = 1; risfdep = 0.0_wp;
	1733	END WHERE
	1734
	1735	IF( icompt == 0 ) THEN
	1736	IF(lwp) WRITE(numout,*)' no points with ice shelf too close to bathymetry'
	1737	ELSE
	1738	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points with ice shelf thickness reduced to avoid bathymetry'
	1739	ENDIF
	1740
	1741	! compute scale factor and depth at T- and W- points
	1742	DO jj = 1, jpj
	1743	DO ji = 1, jpi
	1744	ik = mbathy(ji,jj)
	1745	IF( ik > 0 ) THEN ! ocean point only
	1746	! max ocean level case
	1747	IF( ik == jpkm1 ) THEN
	1748	zdepwp = bathy(ji,jj)
	1749	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
	1750	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
	1751	e3t_0(ji,jj,ik ) = ze3tp
	1752	e3t_0(ji,jj,ik+1) = ze3tp
	1753	e3w_0(ji,jj,ik ) = ze3wp
	1754	e3w_0(ji,jj,ik+1) = ze3tp
	1755	gdepw_0(ji,jj,ik+1) = zdepwp
	1756	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
	1757	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
	1758	!
	1759	ELSE ! standard case
	1760	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
	1761	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
	1762	ENDIF
	1763	! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
	1764	!gm Bug? check the gdepw_1d
	1765	! ... on ik
	1766	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
	1767	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
	1768	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
	1769	e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik )
	1770	e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1)
	1771	! ... on ik+1
	1772	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
	1773	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
	1774	ENDIF
	1775	ENDIF
	1776	END DO
	1777	END DO
	1778	!
	1779	it = 0
	1780	DO jj = 1, jpj
	1781	DO ji = 1, jpi
	1782	ik = mbathy(ji,jj)
	1783	IF( ik > 0 ) THEN ! ocean point only
	1784	e3tp (ji,jj) = e3t_0(ji,jj,ik)
	1785	e3wp (ji,jj) = e3w_0(ji,jj,ik)
	1786	! test
	1787	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
	1788	IF( zdiff <= 0._wp .AND. lwp ) THEN
	1789	it = it + 1
	1790	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
	1791	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
	1792	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
	1793	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
	1794	ENDIF
	1795	ENDIF
	1796	END DO
	1797	END DO
	1798	!
	1799	! (ISF) Definition of e3t, u, v, w for ISF case
	1800	DO jj = 1, jpj
	1801	DO ji = 1, jpi
	1802	ik = misfdep(ji,jj)
	1803	IF( ik > 1 ) THEN ! ice shelf point only
	1804	IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik)
	1805	gdepw_0(ji,jj,ik) = risfdep(ji,jj)
	1806	!gm Bug? check the gdepw_0
	1807	! ... on ik
	1808	gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) &
	1809	& * ( gdepw_1d(ik+1) - gdept_1d(ik) ) &
	1810	& / ( gdepw_1d(ik+1) - gdepw_1d(ik) )
	1811	e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik)
	1812	e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik)
	1813
	1814	IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column)
	1815	e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik)
	1816	ENDIF
	1817	! ... on ik / ik-1
	1818	e3w_0 (ji,jj,ik ) = e3t_0 (ji,jj,ik) !2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik))
[10727]	1819	gdept_0(ji,jj,ik-1) = gdept_0(ji,jj,ik) - e3w_0(ji,jj,ik)
[6951]	1820	e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1)
[10727]	1821	e3w_0 (ji,jj,ik-1) = gdept_0(ji,jj,ik-1) - gdept_1d(ik-2)
	1822	gdepw_0(ji,jj,ik-1) = gdepw_0(ji,jj,ik) - e3t_0(ji,jj,ik-1)
[6951]	1823	ENDIF
	1824	END DO
	1825	END DO
	1826
	1827	it = 0
	1828	DO jj = 1, jpj
	1829	DO ji = 1, jpi
	1830	ik = misfdep(ji,jj)
	1831	IF( ik > 1 ) THEN ! ice shelf point only
	1832	e3tp (ji,jj) = e3t_0(ji,jj,ik )
	1833	e3wp (ji,jj) = e3w_0(ji,jj,ik+1 )
	1834	! test
	1835	zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik )
	1836	IF( zdiff <= 0. .AND. lwp ) THEN
	1837	it = it + 1
	1838	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
	1839	WRITE(numout,*) ' risfdep = ', risfdep(ji,jj)
	1840	WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
	1841	WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj)
	1842	ENDIF
	1843	ENDIF
	1844	END DO
	1845	END DO
	1846
[11129]	1847	DEALLOCATE( zbathy, zmask, zrisfdep )
	1848	DEALLOCATE( zmisfdep, zmbathy )
[6951]	1849	!
	1850	END SUBROUTINE zgr_isf
	1851
	1852
	1853	SUBROUTINE zgr_sco
	1854	!!----------------------------------------------------------------------
	1855	!! * ROUTINE zgr_sco *
	1856	!!
	1857	!! ** Purpose : define the s-coordinate system
	1858	!!
	1859	!! ** Method : s-coordinate
	1860	!! The depth of model levels is defined as the product of an
	1861	!! analytical function by the local bathymetry, while the vertical
	1862	!! scale factors are defined as the product of the first derivative
	1863	!! of the analytical function by the bathymetry.
	1864	!! (this solution save memory as depth and scale factors are not
	1865	!! 3d fields)
	1866	!! - Read bathymetry (in meters) at t-point and compute the
	1867	!! bathymetry at u-, v-, and f-points.
	1868	!! hbatu = mi( hbatt )
	1869	!! hbatv = mj( hbatt )
	1870	!! hbatf = mi( mj( hbatt ) )
	1871	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
	1872	!! function and its derivative given as function.
	1873	!! z_gsigt(k) = fssig (k )
	1874	!! z_gsigw(k) = fssig (k-0.5)
	1875	!! z_esigt(k) = fsdsig(k )
	1876	!! z_esigw(k) = fsdsig(k-0.5)
	1877	!! Three options for stretching are give, and they can be modified
	1878	!! following the users requirements. Nevertheless, the output as
	1879	!! well as the way to compute the model levels and scale factors
	1880	!! must be respected in order to insure second order accuracy
	1881	!! schemes.
	1882	!!
	1883	!! The three methods for stretching available are:
	1884	!!
	1885	!! s_sh94 (Song and Haidvogel 1994)
	1886	!! a sinh/tanh function that allows sigma and stretched sigma
	1887	!!
	1888	!! s_sf12 (Siddorn and Furner 2012?)
	1889	!! allows the maintenance of fixed surface and or
	1890	!! bottom cell resolutions (cf. geopotential coordinates)
	1891	!! within an analytically derived stretched S-coordinate framework.
	1892	!!
	1893	!! s_tanh (Madec et al 1996)
	1894	!! a cosh/tanh function that gives stretched coordinates
	1895	!!
	1896	!!----------------------------------------------------------------------
	1897	INTEGER :: ji, jj, jk, jl ! dummy loop argument
	1898	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
	1899	INTEGER :: ios ! Local integer output status for namelist read
	1900	REAL(wp) :: zrmax, ztaper ! temporary scalars
	1901	REAL(wp) :: zrfact
	1902	!
[11129]	1903	REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
	1904	REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
[6951]	1905	!!
	1906	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
	1907	& rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
	1908	!!----------------------------------------------------------------------
	1909	!
[11129]	1910	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]	1911	!
	1912	REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
	1913	READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
	1914	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
	1915
	1916	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
	1917	READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
	1918	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
	1919	IF(lwm) WRITE ( numond, namzgr_sco )
	1920
	1921	IF(lwp) THEN ! control print
	1922	WRITE(numout,*)
	1923	WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
	1924	WRITE(numout,*) '~~~~~~~~~~~'
	1925	WRITE(numout,*) ' Namelist namzgr_sco'
	1926	WRITE(numout,*) ' stretching coeffs '
	1927	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
	1928	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
	1929	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
	1930	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
	1931	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
	1932	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
	1933	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
	1934	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
	1935	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
	1936	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
	1937	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
	1938	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
	1939	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
	1940	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
	1941	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
	1942	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
	1943	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
	1944	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
	1945	ENDIF
	1946
	1947	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
	1948	hifu(:,:) = rn_sbot_min
	1949	hifv(:,:) = rn_sbot_min
	1950	hiff(:,:) = rn_sbot_min
	1951
	1952	! ! set maximum ocean depth
	1953	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
	1954
	1955	DO jj = 1, jpj
	1956	DO ji = 1, jpi
	1957	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
	1958	END DO
	1959	END DO
	1960	! ! =============================
	1961	! ! Define the envelop bathymetry (hbatt)
	1962	! ! =============================
	1963	! use r-value to create hybrid coordinates
	1964	zenv(:,:) = bathy(:,:)
	1965	!
	1966	! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
	1967	DO jj = 1, jpj
	1968	DO ji = 1, jpi
	1969	IF( bathy(ji,jj) == 0._wp ) THEN
	1970	iip1 = MIN( ji+1, jpi )
	1971	ijp1 = MIN( jj+1, jpj )
	1972	iim1 = MAX( ji-1, 1 )
	1973	ijm1 = MAX( jj-1, 1 )
	1974	!!gm BUG fix see ticket #1617
	1975	IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) &
	1976	& + bathy(iim1,jj ) + bathy(iip1,jj ) &
	1977	& + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) &
	1978	& zenv(ji,jj) = rn_sbot_min
	1979	!!gm
	1980	!!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + &
	1981	!!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
	1982	!!gm zenv(ji,jj) = rn_sbot_min
	1983	!!gm ENDIF
	1984	!!gm end
	1985	ENDIF
	1986	END DO
	1987	END DO
	1988
	1989	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
[11129]	1990	CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' )
[6951]	1991	!
	1992	! smooth the bathymetry (if required)
	1993	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
	1994	scobot(:,:) = bathy(:,:) ! ocean bottom depth
	1995	!
	1996	jl = 0
	1997	zrmax = 1._wp
	1998	!
	1999	!
	2000	! set scaling factor used in reducing vertical gradients
	2001	zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
	2002	!
	2003	! initialise temporary evelope depth arrays
	2004	ztmpi1(:,:) = zenv(:,:)
	2005	ztmpi2(:,:) = zenv(:,:)
	2006	ztmpj1(:,:) = zenv(:,:)
	2007	ztmpj2(:,:) = zenv(:,:)
	2008	!
	2009	! initialise temporary r-value arrays
	2010	zri(:,:) = 1._wp
	2011	zrj(:,:) = 1._wp
	2012	! ! ================ !
	2013	DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop !
	2014	! ! ================ !
	2015	jl = jl + 1
	2016	zrmax = 0._wp
	2017	! we set zrmax from previous r-values (zri and zrj) first
	2018	! if set after current r-value calculation (as previously)
	2019	! we could exit DO WHILE prematurely before checking r-value
	2020	! of current zenv
	2021	DO jj = 1, nlcj
	2022	DO ji = 1, nlci
	2023	zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
	2024	END DO
	2025	END DO
	2026	zri(:,:) = 0._wp
	2027	zrj(:,:) = 0._wp
	2028	DO jj = 1, nlcj
	2029	DO ji = 1, nlci
	2030	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
	2031	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
	2032	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
	2033	zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
	2034	END IF
	2035	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
	2036	zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
	2037	END IF
	2038	IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact
	2039	IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact
	2040	IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact
	2041	IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact
	2042	END DO
	2043	END DO
[10727]	2044	! IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
[6951]	2045	!
	2046	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax
	2047	!
	2048	DO jj = 1, nlcj
	2049	DO ji = 1, nlci
	2050	zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
	2051	END DO
	2052	END DO
	2053	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
[11129]	2054	CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' )
[6951]	2055	! ! ================ !
	2056	END DO ! End loop !
	2057	! ! ================ !
	2058	DO jj = 1, jpj
	2059	DO ji = 1, jpi
	2060	zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
	2061	END DO
	2062	END DO
	2063	!
	2064	! Envelope bathymetry saved in hbatt
	2065	hbatt(:,:) = zenv(:,:)
	2066	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
	2067	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
	2068	DO jj = 1, jpj
	2069	DO ji = 1, jpi
	2070	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
	2071	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
	2072	END DO
	2073	END DO
	2074	ENDIF
	2075	!
	2076	! ! ==============================
	2077	! ! hbatu, hbatv, hbatf fields
	2078	! ! ==============================
	2079	IF(lwp) THEN
	2080	WRITE(numout,*)
	2081	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
	2082	ENDIF
	2083	hbatu(:,:) = rn_sbot_min
	2084	hbatv(:,:) = rn_sbot_min
	2085	hbatf(:,:) = rn_sbot_min
	2086	DO jj = 1, jpjm1
	2087	DO ji = 1, jpim1 ! NO vector opt.
	2088	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
	2089	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
	2090	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
	2091	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
	2092	END DO
	2093	END DO
	2094
	2095	!
	2096	! Apply lateral boundary condition
	2097	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
[11129]	2098	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp )
[6951]	2099	DO jj = 1, jpj
	2100	DO ji = 1, jpi
	2101	IF( hbatu(ji,jj) == 0._wp ) THEN
	2102	!No worries about the following line when ln_wd == .true.
	2103	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
	2104	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
	2105	ENDIF
	2106	END DO
	2107	END DO
[11129]	2108	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp )
[6951]	2109	DO jj = 1, jpj
	2110	DO ji = 1, jpi
	2111	IF( hbatv(ji,jj) == 0._wp ) THEN
	2112	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
	2113	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
	2114	ENDIF
	2115	END DO
	2116	END DO
[11129]	2117	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp )
[6951]	2118	DO jj = 1, jpj
	2119	DO ji = 1, jpi
	2120	IF( hbatf(ji,jj) == 0._wp ) THEN
	2121	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
	2122	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
	2123	ENDIF
	2124	END DO
	2125	END DO
	2126
	2127	!!bug: key_helsinki a verifer
	2128	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
	2129	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
	2130	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
	2131	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
	2132
	2133	IF( nprint == 1 .AND. lwp ) THEN
	2134	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
	2135	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
	2136	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
	2137	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
	2138	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
	2139	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
	2140	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
	2141	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
	2142	ENDIF
	2143	!! helsinki
	2144
	2145	! ! =======================
	2146	! ! s-ccordinate fields (gdep., e3.)
	2147	! ! =======================
	2148	!
	2149	! non-dimensional "sigma" for model level depth at w- and t-levels
	2150
	2151
	2152	!========================================================================
	2153	! Song and Haidvogel 1994 (ln_s_sh94=T)
	2154	! Siddorn and Furner 2012 (ln_sf12=T)
	2155	! or tanh function (both false)
	2156	!========================================================================
	2157	IF ( ln_s_sh94 ) THEN
	2158	CALL s_sh94()
	2159	ELSE IF ( ln_s_sf12 ) THEN
	2160	CALL s_sf12()
	2161	ELSE
	2162	CALL s_tanh()
	2163	ENDIF
	2164
[11129]	2165	CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp )
	2166	CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp )
	2167	CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp )
	2168	CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp )
	2169	CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp )
	2170	CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp )
	2171	CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
[6951]	2172	!
	2173	WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp
	2174	WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp
	2175	WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp
	2176	WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp
	2177	WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp
	2178	WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp
	2179	WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp
	2180	!!
	2181	! HYBRID :
	2182	DO jj = 1, jpj
	2183	DO ji = 1, jpi
	2184	DO jk = 1, jpkm1
[11129]	2185	IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
[6951]	2186	END DO
	2187	END DO
	2188	END DO
	2189	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
	2190	& ' MAX ', MAXVAL( mbathy(:,:) )
	2191
	2192	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
	2193	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	2194	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
	2195	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) )
	2196	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
	2197	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
	2198	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
	2199	& ' w ', MINVAL( e3w_0 (:,:,:) )
	2200
	2201	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
	2202	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) )
	2203	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
	2204	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
	2205	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
	2206	& ' w ', MAXVAL( e3w_0 (:,:,:) )
	2207	ENDIF
	2208	! END DO
	2209	IF(lwp) THEN ! selected vertical profiles
	2210	WRITE(numout,*)
	2211	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
	2212	WRITE(numout,*) ' ~~~~~~ --------------------'
	2213	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
	2214	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
	2215	& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
	2216	DO jj = mj0(20), mj1(20)
	2217	DO ji = mi0(20), mi1(20)
	2218	WRITE(numout,*)
	2219	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	2220	WRITE(numout,*) ' ~~~~~~ --------------------'
	2221	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
	2222	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
	2223	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
	2224	END DO
	2225	END DO
	2226	DO jj = mj0(74), mj1(74)
	2227	DO ji = mi0(100), mi1(100)
	2228	WRITE(numout,*)
	2229	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	2230	WRITE(numout,*) ' ~~~~~~ --------------------'
	2231	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
	2232	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
	2233	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
	2234	END DO
	2235	END DO
	2236	ENDIF
	2237	!
	2238	!================================================================================
	2239	! check the coordinate makes sense
	2240	!================================================================================
	2241	DO ji = 1, jpi
	2242	DO jj = 1, jpj
	2243	!
	2244	IF( hbatt(ji,jj) > 0._wp) THEN
	2245	DO jk = 1, mbathy(ji,jj)
	2246	! check coordinate is monotonically increasing
[11129]	2247	IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
[6951]	2248	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
	2249	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
[11129]	2250	WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
	2251	WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
[6951]	2252	CALL ctl_stop( ctmp1 )
	2253	ENDIF
	2254	! and check it has never gone negative
[11129]	2255	IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
[6951]	2256	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
	2257	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
[11129]	2258	WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
	2259	WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
[6951]	2260	CALL ctl_stop( ctmp1 )
	2261	ENDIF
	2262	! and check it never exceeds the total depth
[11129]	2263	IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
[6951]	2264	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
	2265	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
[11129]	2266	WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
[6951]	2267	CALL ctl_stop( ctmp1 )
	2268	ENDIF
	2269	END DO
	2270	!
	2271	DO jk = 1, mbathy(ji,jj)-1
	2272	! and check it never exceeds the total depth
[11129]	2273	IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
[6951]	2274	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
	2275	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
[11129]	2276	WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
[6951]	2277	CALL ctl_stop( ctmp1 )
	2278	ENDIF
	2279	END DO
	2280	ENDIF
	2281	END DO
	2282	END DO
	2283	!
[11129]	2284	DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
[6951]	2285	!
	2286	END SUBROUTINE zgr_sco
	2287
	2288
	2289	SUBROUTINE s_sh94()
	2290	!!----------------------------------------------------------------------
	2291	!! * ROUTINE s_sh94 *
	2292	!!
	2293	!! ** Purpose : stretch the s-coordinate system
	2294	!!
	2295	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
	2296	!! mixed S/sigma coordinate
	2297	!!
	2298	!! Reference : Song and Haidvogel 1994.
	2299	!!----------------------------------------------------------------------
	2300	INTEGER :: ji, jj, jk ! dummy loop argument
	2301	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
	2302	REAL(wp) :: ztmpu, ztmpv, ztmpf
	2303	REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
	2304	!
[11129]	2305	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
	2306	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
[6951]	2307	!!----------------------------------------------------------------------
	2308
[11129]	2309	ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk), z_gsi3w3 (jpi,jpj,jpk) )
	2310	ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) )
	2311	ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
[6951]	2312
	2313	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
	2314	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
	2315	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
	2316	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
	2317	!
	2318	DO ji = 1, jpi
	2319	DO jj = 1, jpj
	2320	!
	2321	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
	2322	DO jk = 1, jpk
	2323	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
	2324	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
	2325	END DO
	2326	ELSE ! shallow water, uniform sigma
	2327	DO jk = 1, jpk
	2328	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
	2329	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
	2330	END DO
	2331	ENDIF
	2332	!
	2333	DO jk = 1, jpkm1
	2334	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
	2335	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
	2336	END DO
	2337	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
	2338	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
	2339	!
	2340	! Coefficients for vertical depth as the sum of e3w scale factors
	2341	z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
	2342	DO jk = 2, jpk
	2343	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
	2344	END DO
	2345	!
	2346	DO jk = 1, jpk
	2347	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
	2348	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
	2349	gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
	2350	gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
	2351	gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsi3w3(ji,jj,jk)+rn_hczcoeft )
	2352	END DO
	2353	!
	2354	END DO ! for all jj's
	2355	END DO ! for all ji's
	2356
	2357	DO ji = 1, jpim1
	2358	DO jj = 1, jpjm1
	2359	! extended for Wetting/Drying case
	2360	ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
	2361	ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
	2362	ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
	2363	ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
	2364	ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
	2365	ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
	2366	& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
	2367	DO jk = 1, jpk
	2368	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
	2369	& / ztmpu
	2370	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
	2371	& / ztmpu
	2372
	2373	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
	2374	& / ztmpv
	2375	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
	2376	& / ztmpv
	2377
	2378	z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
	2379	& + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
	2380	& + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
	2381	& + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
	2382
	2383	!
	2384	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	2385	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	2386	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	2387	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	2388	!
	2389	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	2390	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	2391	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	2392	END DO
	2393	END DO
	2394	END DO
	2395	!
[11129]	2396	DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3 )
	2397	DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
[6951]	2398	!
	2399	END SUBROUTINE s_sh94
	2400
	2401
	2402	SUBROUTINE s_sf12
	2403	!!----------------------------------------------------------------------
	2404	!! * ROUTINE s_sf12 *
	2405	!!
	2406	!! ** Purpose : stretch the s-coordinate system
	2407	!!
	2408	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
	2409	!! mixed S/sigma/Z coordinate
	2410	!!
	2411	!! This method allows the maintenance of fixed surface and or
	2412	!! bottom cell resolutions (cf. geopotential coordinates)
	2413	!! within an analytically derived stretched S-coordinate framework.
	2414	!!
	2415	!!
	2416	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
	2417	!!----------------------------------------------------------------------
	2418	INTEGER :: ji, jj, jk ! dummy loop argument
	2419	REAL(wp) :: zsmth ! smoothing around critical depth
	2420	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
	2421	REAL(wp) :: ztmpu, ztmpv, ztmpf
	2422	REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
	2423	!
[11129]	2424	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
	2425	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
[6951]	2426	!!----------------------------------------------------------------------
	2427	!
[11129]	2428	ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk), z_gsi3w3 (jpi,jpj,jpk) )
	2429	ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk))
	2430	ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
[6951]	2431
	2432	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
	2433	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
	2434	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
	2435	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
	2436
	2437	DO ji = 1, jpi
	2438	DO jj = 1, jpj
	2439
	2440	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
	2441
	2442	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
	2443	! could be changed by users but care must be taken to do so carefully
	2444	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
	2445
	2446	zzs = rn_zs / hbatt(ji,jj)
	2447
	2448	IF (rn_efold /= 0.0_wp) THEN
	2449	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
	2450	ELSE
	2451	zsmth = 1.0_wp
	2452	ENDIF
	2453
	2454	DO jk = 1, jpk
	2455	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
	2456	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
	2457	ENDDO
	2458	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
	2459	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
	2460
	2461	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
	2462
	2463	DO jk = 1, jpk
	2464	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
	2465	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
	2466	END DO
	2467
	2468	ELSE ! shallow water, z coordinates
	2469
	2470	DO jk = 1, jpk
	2471	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
	2472	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
	2473	END DO
	2474
	2475	ENDIF
	2476
	2477	DO jk = 1, jpkm1
	2478	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
	2479	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
	2480	END DO
	2481	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
	2482	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
	2483
	2484	! Coefficients for vertical depth as the sum of e3w scale factors
	2485	z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
	2486	DO jk = 2, jpk
	2487	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
	2488	END DO
	2489
	2490	DO jk = 1, jpk
	2491	gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
	2492	gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
	2493	gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
	2494	END DO
	2495
	2496	ENDDO ! for all jj's
	2497	ENDDO ! for all ji's
	2498
	2499	DO ji=1,jpi-1
	2500	DO jj=1,jpj-1
	2501
	2502	! extend to suit for Wetting/Drying case
	2503	ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
	2504	ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
	2505	ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
	2506	ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
	2507	ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
	2508	ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
	2509	& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
	2510	DO jk = 1, jpk
	2511	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
	2512	& / ztmpu
	2513	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
	2514	& / ztmpu
	2515	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
	2516	& / ztmpv
	2517	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
	2518	& / ztmpv
	2519
	2520	z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
	2521	& + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
	2522	& + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
	2523	& + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
	2524
	2525	! Code prior to wetting and drying option (for reference)
	2526	!z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
	2527	! /( hbatt(ji,jj)+hbatt(ji+1,jj) )
	2528	!
	2529	!z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
	2530	! /( hbatt(ji,jj)+hbatt(ji+1,jj) )
	2531	!
	2532	!z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
	2533	! /( hbatt(ji,jj)+hbatt(ji,jj+1) )
	2534	!
	2535	!z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
	2536	! /( hbatt(ji,jj)+hbatt(ji,jj+1) )
	2537	!
	2538	!z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
	2539	! & +hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
	2540	! +hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
	2541	! & +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) &
	2542	! /( hbatt(ji ,jj )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
	2543
	2544	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
	2545	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
	2546	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
	2547	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
	2548	!
	2549	e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
	2550	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
	2551	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
	2552	END DO
	2553
	2554	ENDDO
	2555	ENDDO
	2556	!
[11129]	2557	CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.)
	2558	CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.)
	2559	CALL lbc_lnk('domzgr',e3w_0 ,'T',1.)
	2560	CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.)
[6951]	2561	!
[11129]	2562	DEALLOCATE( z_gsigw3, z_gsigt3, z_gsi3w3 )
	2563	DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
[6951]	2564	!
	2565	END SUBROUTINE s_sf12
	2566
	2567
	2568	SUBROUTINE s_tanh()
	2569	!!----------------------------------------------------------------------
	2570	!! * ROUTINE s_tanh*
	2571	!!
	2572	!! ** Purpose : stretch the s-coordinate system
	2573	!!
	2574	!! ** Method : s-coordinate stretch
	2575	!!
	2576	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
	2577	!!----------------------------------------------------------------------
	2578	INTEGER :: ji, jj, jk ! dummy loop argument
	2579	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
[11129]	2580	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
	2581	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
[6951]	2582	!!----------------------------------------------------------------------
	2583
[11129]	2584	ALLOCATE( z_gsigw(jpk), z_gsigt(jpk), z_gsi3w(jpk) )
	2585	ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
[6951]	2586
	2587	z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp
	2588	z_esigt = 0._wp ; z_esigw = 0._wp
	2589
	2590	DO jk = 1, jpk
	2591	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
	2592	z_gsigt(jk) = -fssig( REAL(jk,wp) )
	2593	END DO
	2594	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
	2595	!
	2596	! Coefficients for vertical scale factors at w-, t- levels
	2597	!!gm bug : define it from analytical function, not like juste bellow....
	2598	!!gm or betteroffer the 2 possibilities....
	2599	DO jk = 1, jpkm1
	2600	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
	2601	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
	2602	END DO
	2603	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
	2604	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
	2605	!
	2606	! Coefficients for vertical depth as the sum of e3w scale factors
	2607	z_gsi3w(1) = 0.5_wp * z_esigw(1)
	2608	DO jk = 2, jpk
	2609	z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
	2610	END DO
	2611	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
	2612	DO jk = 1, jpk
	2613	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
	2614	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
	2615	gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
	2616	gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
	2617	gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsi3w(jk) + hift(:,:)zcoeft )
	2618	END DO
	2619	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
	2620	DO jj = 1, jpj
	2621	DO ji = 1, jpi
	2622	DO jk = 1, jpk
	2623	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
	2624	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
	2625	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
	2626	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
	2627	!
	2628	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
	2629	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
	2630	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
	2631	END DO
	2632	END DO
	2633	END DO
	2634	!
[11129]	2635	DEALLOCATE( z_gsigw, z_gsigt, z_gsi3w )
	2636	DEALLOCATE( z_esigt, z_esigw )
[6951]	2637	!
	2638	END SUBROUTINE s_tanh
	2639
	2640
	2641	FUNCTION fssig( pk ) RESULT( pf )
	2642	!!----------------------------------------------------------------------
	2643	!! * ROUTINE fssig *
	2644	!!
	2645	!! ** Purpose : provide the analytical function in s-coordinate
	2646	!!
	2647	!! ** Method : the function provide the non-dimensional position of
	2648	!! T and W (i.e. between 0 and 1)
	2649	!! T-points at integer values (between 1 and jpk)
	2650	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	2651	!!----------------------------------------------------------------------
	2652	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
	2653	REAL(wp) :: pf ! sigma value
	2654	!!----------------------------------------------------------------------
	2655	!
	2656	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
	2657	& - TANH( rn_thetb * rn_theta ) ) &
	2658	& * ( COSH( rn_theta ) &
	2659	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
	2660	& / ( 2._wp * SINH( rn_theta ) )
	2661	!
	2662	END FUNCTION fssig
	2663
	2664
	2665	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
	2666	!!----------------------------------------------------------------------
	2667	!! * ROUTINE fssig1 *
	2668	!!
	2669	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
	2670	!!
	2671	!! ** Method : the function provides the non-dimensional position of
	2672	!! T and W (i.e. between 0 and 1)
	2673	!! T-points at integer values (between 1 and jpk)
	2674	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	2675	!!----------------------------------------------------------------------
	2676	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
	2677	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
	2678	REAL(wp) :: pf1 ! sigma value
	2679	!!----------------------------------------------------------------------
	2680	!
	2681	IF ( rn_theta == 0 ) then ! uniform sigma
	2682	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
	2683	ELSE ! stretched sigma
	2684	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
	2685	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
	2686	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
	2687	ENDIF
	2688	!
	2689	END FUNCTION fssig1
	2690
	2691
	2692	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
	2693	!!----------------------------------------------------------------------
	2694	!! * ROUTINE fgamma *
	2695	!!
	2696	!! ** Purpose : provide analytical function for the s-coordinate
	2697	!!
	2698	!! ** Method : the function provides the non-dimensional position of
	2699	!! T and W (i.e. between 0 and 1)
	2700	!! T-points at integer values (between 1 and jpk)
	2701	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	2702	!!
	2703	!! This method allows the maintenance of fixed surface and or
	2704	!! bottom cell resolutions (cf. geopotential coordinates)
	2705	!! within an analytically derived stretched S-coordinate framework.
	2706	!!
	2707	!! Reference : Siddorn and Furner, in prep
	2708	!!----------------------------------------------------------------------
	2709	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
	2710	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
	2711	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
	2712	REAL(wp), INTENT(in ) :: pzs ! surface box depth
	2713	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
	2714	!
	2715	INTEGER :: jk ! dummy loop index
	2716	REAL(wp) :: za1,za2,za3 ! local scalar
	2717	REAL(wp) :: zn1,zn2 ! - -
	2718	REAL(wp) :: za,zb,zx ! - -
	2719	!!----------------------------------------------------------------------
	2720	!
	2721	zn1 = 1._wp / REAL( jpkm1, wp )
	2722	zn2 = 1._wp - zn1
	2723	!
	2724	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
	2725	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
	2726	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
	2727	!
	2728	za = pzb - za3*(pzs-za1)-za2
	2729	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
	2730	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
	2731	zx = 1.0_wp-za/2.0_wp-zb
	2732	!
	2733	DO jk = 1, jpk
	2734	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
	2735	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
	2736	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
	2737	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
	2738	END DO
	2739	!
	2740	END FUNCTION fgamma
	2741
	2742	!!======================================================================
	2743	END MODULE domzgr

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

Download in other formats: