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

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

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

Context Navigation

source: branches/2016/dev_r6409_SIMPLIF_2_usrdef/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 6624

Last change on this file since 6624 was 6624, checked in by gm, 8 years ago
#1692 - branch SIMPLIF_2_usrdef: add domain_cfg.nc file which includes jperio, and vert. coord. logicals. +code cleaning
Property svn:keywords set to `Id`
File size: 139.3 KB

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

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

Download in other formats: