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

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

domzgr.F90 in branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 7412

Last change on this file since 7412 was 7412, checked in by lovato, 8 years ago
Merge dev_NOC_CMCC_merge_2016 into branch
Property svn:keywords set to `Id`
File size: 147.5 KB

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

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

Download in other formats: