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/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/2015/dev_r5803_UKMO_AGRIF_Vert_interp/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 7836

Last change on this file since 7836 was 7825, checked in by timgraham, 7 years ago
Extra temporary change in domzgr for my GYRE_AGRIF config Commented out Jerome's old hack for different jpk in nemogcm.F90
Property svn:keywords set to `Id`
File size: 141.0 KB

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

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

Download in other formats: