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

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

domzgr.F90 in branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 9316

Last change on this file since 9316 was 9316, checked in by mathiot, 6 years ago
update domzgr and namelist (rewrite of zgr_isf based on discussion with P. Holland, R. Smith, A. Siaahan and J. Gregory)
File size: 124.9 KB

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

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

Download in other formats: