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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6485

Last change on this file since 6485 was 6152, checked in by acc, 8 years ago
Add wetting and drying option from dev_r5803_NOC_WAD branch. Logically isolated code changes in domvvl.F90, domzgr.F90, dynhpg.F90, dynspg_ts.F90, sshwzv.F90 and nemogcm.F90. New module wet_dry.F90 in DYN. Fully SETTE tested with code deactivated (ln_wad=.false.). No test case yet available to justify activating option (still under development)
Property svn:keywords set to `Id`
File size: 139.6 KB

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

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

Download in other formats: