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 utils/tools/DOMAINcfg/src – NEMO

Context Navigation

source: utils/tools/DOMAINcfg/src/domzgr.F90 @ 12414

Last change on this file since 12414 was 12414, checked in by smueller, 4 years ago
Reintegration of 2019 development branch ^{/utils/tools_MERGE_2019 into the tools directory (}/utils/tools)
File size: 98.1 KB

Rev	Line
[6951]	1	MODULE domzgr
	2	!!==============================================================================
	3	!! * MODULE domzgr *
	4	!! Ocean domain : definition of the vertical coordinate system
	5	!!==============================================================================
	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)
	9	!! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module
	10	!! - ! 2002-09 (A. de Miranda) rigid-lid + islands
	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
	16	!! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level
	17	!! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function
	18	!! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case
[12101]	19	!! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capability
[6951]	20	!! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying
	21	!!----------------------------------------------------------------------
	22
	23	!!----------------------------------------------------------------------
	24	!! dom_zgr : defined the ocean vertical coordinate system
	25	!! zgr_bat : bathymetry fields (levels and meters)
	26	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
	27	!! zgr_bat_ctl : check the bathymetry files
	28	!! zgr_bot_level: deepest ocean level for t-, u, and v-points
	29	!! zgr_z : reference z-coordinate
	30	!! zgr_zco : z-coordinate
	31	!! zgr_zps : z-coordinate with partial steps
	32	!! zgr_sco : s-coordinate
	33	!! fssig : tanh stretch function
	34	!! fssig1 : Song and Haidvogel 1994 stretch function
	35	!! fgamma : Siddorn and Furner 2012 stretching function
	36	!!---------------------------------------------------------------------
	37	USE dom_oce ! ocean domain
	38	!
	39	USE in_out_manager ! I/O manager
	40	USE iom ! I/O library
	41	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	42	USE lib_mpp ! distributed memory computing library
[12101]	43	USE lib_fortran
[10727]	44	USE dombat
[12101]	45	USE domisf
[6951]	46
	47	IMPLICIT NONE
	48	PRIVATE
	49
	50	PUBLIC dom_zgr ! called by dom_init.F90
[12101]	51	!
[6951]	52	! !!* Namelist namzgr_sco *
	53	LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
	54	LOGICAL :: ln_s_sf12 ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
	55	!
	56	REAL(wp) :: rn_sbot_min ! minimum depth of s-bottom surface (>0) (m)
	57	REAL(wp) :: rn_sbot_max ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
	58	REAL(wp) :: rn_rmax ! maximum cut-off r-value allowed (0<rn_rmax<1)
	59	REAL(wp) :: rn_hc ! Critical depth for transition from sigma to stretched coordinates
[10727]	60
[6951]	61	! Song and Haidvogel 1994 stretching parameters
	62	REAL(wp) :: rn_theta ! surface control parameter (0<=rn_theta<=20)
	63	REAL(wp) :: rn_thetb ! bottom control parameter (0<=rn_thetb<= 1)
	64	REAL(wp) :: rn_bb ! stretching parameter
	65	! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
	66	! Siddorn and Furner stretching parameters
	67	LOGICAL :: ln_sigcrit ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
	68	REAL(wp) :: rn_alpha ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
	69	REAL(wp) :: rn_efold ! efold length scale for transition to stretched coord
	70	REAL(wp) :: rn_zs ! depth of surface grid box
	71	! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
	72	REAL(wp) :: rn_zb_a ! bathymetry scaling factor for calculating Zb
	73	REAL(wp) :: rn_zb_b ! offset for calculating Zb
	74
	75	!! * Substitutions
	76	!!----------------------------------------------------------------------
	77	!! * vectopt_loop_substitute *
	78	!!----------------------------------------------------------------------
	79	!! ** purpose : substitute the inner loop start/end indices with CPP macro
	80	!! allow unrolling of do-loop (useful with vector processors)
	81	!!----------------------------------------------------------------------
	82	!!----------------------------------------------------------------------
	83	!! NEMO/OPA 3.7 , NEMO Consortium (2014)
	84	!! $Id: vectopt_loop_substitute.h90 4990 2014-12-15 16:42:49Z timgraham $
[9598]	85	!! Software governed by the CeCILL licence (./LICENSE)
[6951]	86	!!----------------------------------------------------------------------
	87	!!----------------------------------------------------------------------
	88	!! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
[7107]	89	!! $Id: domzgr.F90 6827 2016-08-01 13:37:15Z flavoni $
[9598]	90	!! Software governed by the CeCILL licence (./LICENSE)
[6951]	91	!!----------------------------------------------------------------------
	92	CONTAINS
	93
	94	SUBROUTINE dom_zgr
	95	!!----------------------------------------------------------------------
	96	!! * ROUTINE dom_zgr *
	97	!!
	98	!! ** Purpose : set the depth of model levels and the resulting
	99	!! vertical scale factors.
	100	!!
	101	!! ** Method : - reference 1D vertical coordinate (gdep._1d, e3._1d)
	102	!! - read/set ocean depth and ocean levels (bathy, mbathy)
	103	!! - vertical coordinate (gdep., e3.) depending on the
	104	!! coordinate chosen :
	105	!! ln_zco=T z-coordinate
	106	!! ln_zps=T z-coordinate with partial steps
	107	!! ln_zco=T s-coordinate
	108	!!
	109	!! ** Action : define gdep., e3., mbathy and bathy
	110	!!----------------------------------------------------------------------
	111	INTEGER :: ioptio, ibat ! local integer
	112	INTEGER :: ios
	113	!
	114	NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh
	115	!!----------------------------------------------------------------------
	116	!
	117	!
	118	REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate
	119	READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
	120	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
	121
	122	REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate
	123	READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
	124	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
	125	IF(lwm) WRITE ( numond, namzgr )
	126
	127	IF(lwp) THEN ! Control print
	128	WRITE(numout,*)
	129	WRITE(numout,*) 'dom_zgr : vertical coordinate'
	130	WRITE(numout,*) '~~~~~~~'
	131	WRITE(numout,*) ' Namelist namzgr : set vertical coordinate'
	132	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
	133	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
	134	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
	135	WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav
	136	WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh
	137	ENDIF
	138
	139	IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time'
	140
	141	ioptio = 0 ! Check Vertical coordinate options
	142	IF( ln_zco ) ioptio = ioptio + 1
	143	IF( ln_zps ) ioptio = ioptio + 1
	144	IF( ln_sco ) ioptio = ioptio + 1
	145	IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
	146	!
[12101]	147	IF ( ln_isfcav ) CALL zgr_isf_nam
[6951]	148	ioptio = 0
	149	IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1
	150	IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1
	151	IF( ioptio > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' )
	152	!
	153	! Build the vertical coordinate system
	154	! ------------------------------------
	155	CALL zgr_z ! Reference z-coordinate system (always called)
	156	CALL zgr_bat ! Bathymetry fields (levels and meters)
	157	IF( ln_zco ) CALL zgr_zco ! z-coordinate
	158	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
	159	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
[12101]	160	CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points
[6951]	161	!
	162	! final adjustment of mbathy & check
	163	! -----------------------------------
	164	CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
	165	CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points
	166	!
	167	IF( nprint == 1 .AND. lwp ) THEN
	168	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	169	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
[12101]	170	& ' w ', MINVAL( gdepw_0(:,:,:) )
[6951]	171	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), &
	172	& ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), &
	173	& ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), &
	174	& ' w ', MINVAL( e3w_0(:,:,:) )
	175
	176	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
[12101]	177	& ' w ', MAXVAL( gdepw_0(:,:,:) )
[6951]	178	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), &
	179	& ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), &
	180	& ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), &
	181	& ' w ', MAXVAL( e3w_0(:,:,:) )
	182	ENDIF
	183	!
	184	END SUBROUTINE dom_zgr
	185
	186
	187	SUBROUTINE zgr_z
	188	!!----------------------------------------------------------------------
	189	!! * ROUTINE zgr_z *
	190	!!
	191	!! ** Purpose : set the depth of model levels and the resulting
	192	!! vertical scale factors.
	193	!!
	194	!! ** Method : z-coordinate system (use in all type of coordinate)
	195	!! The depth of model levels is defined from an analytical
	196	!! function the derivative of which gives the scale factors.
	197	!! both depth and scale factors only depend on k (1d arrays).
	198	!! w-level: gdepw_1d = gdep(k)
	199	!! e3w_1d(k) = dk(gdep)(k) = e3(k)
	200	!! t-level: gdept_1d = gdep(k+0.5)
	201	!! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
	202	!!
	203	!! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
	204	!! - e3t_1d , e3w_1d : scale factors at T- and W-levels (m)
	205	!!
	206	!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
	207	!!----------------------------------------------------------------------
	208	INTEGER :: jk ! dummy loop indices
	209	REAL(wp) :: zt, zw ! temporary scalars
	210	REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
	211	REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
	212	REAL(wp) :: zrefdep ! depth of the reference level (~10m)
	213	REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
	214	!!----------------------------------------------------------------------
	215	!
	216	! Set variables from parameters
	217	! ------------------------------
	218	zkth = ppkth ; zacr = ppacr
	219	zdzmin = ppdzmin ; zhmax = pphmax
	220	zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters
	221
	222	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
	223	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
	224	IF( ppa1 == pp_to_be_computed .AND. &
	225	& ppa0 == pp_to_be_computed .AND. &
	226	& ppsur == pp_to_be_computed ) THEN
	227	!
	228	za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
	229	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
	230	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
	231	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
	232	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
	233	ELSE
	234	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
	235	za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter
	236	ENDIF
	237
	238	IF(lwp) THEN ! Parameter print
	239	WRITE(numout,*)
	240	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
	241	WRITE(numout,*) ' ~~~~~~~'
	242	IF( ppkth == 0._wp ) THEN
	243	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
	244	WRITE(numout,*) ' Total depth :', zhmax
	245	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
	246	ELSE
	247	IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
	248	WRITE(numout,*) ' zsur, za0, za1 computed from '
	249	WRITE(numout,*) ' zdzmin = ', zdzmin
	250	WRITE(numout,*) ' zhmax = ', zhmax
	251	ENDIF
	252	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
	253	WRITE(numout,*) ' zsur = ', zsur
	254	WRITE(numout,*) ' za0 = ', za0
	255	WRITE(numout,*) ' za1 = ', za1
	256	WRITE(numout,*) ' zkth = ', zkth
	257	WRITE(numout,*) ' zacr = ', zacr
	258	IF( ldbletanh ) THEN
	259	WRITE(numout,*) ' (Double tanh za2 = ', za2
	260	WRITE(numout,*) ' parameters) zkth2= ', zkth2
	261	WRITE(numout,*) ' zacr2= ', zacr2
	262	ENDIF
	263	ENDIF
	264	ENDIF
	265
	266
	267	! Reference z-coordinate (depth - scale factor at T- and W-points)
	268	! ======================
	269	IF( ppkth == 0._wp ) THEN ! uniform vertical grid
	270
	271
	272
	273	za1 = zhmax / FLOAT(jpk-1)
	274
	275	DO jk = 1, jpk
	276	zw = FLOAT( jk )
	277	zt = FLOAT( jk ) + 0.5_wp
	278	gdepw_1d(jk) = ( zw - 1 ) * za1
	279	gdept_1d(jk) = ( zt - 1 ) * za1
	280	e3w_1d (jk) = za1
	281	e3t_1d (jk) = za1
	282	END DO
	283	ELSE ! Madec & Imbard 1996 function
	284	IF( .NOT. ldbletanh ) THEN
	285	DO jk = 1, jpk
	286	zw = REAL( jk , wp )
	287	zt = REAL( jk , wp ) + 0.5_wp
	288	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
	289	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
	290	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
	291	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
	292	END DO
	293	ELSE
	294	DO jk = 1, jpk
	295	zw = FLOAT( jk )
	296	zt = FLOAT( jk ) + 0.5_wp
	297	! Double tanh function
	298	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) &
	299	& + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) )
	300	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) &
	301	& + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) )
	302	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) &
	303	& + za2 * TANH( (zw-zkth2) / zacr2 )
	304	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) &
	305	& + za2 * TANH( (zt-zkth2) / zacr2 )
	306	END DO
	307	ENDIF
	308	gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero
	309	ENDIF
	310
[7200]	311	IF ( ln_isfcav .OR. ln_e3_dep ) THEN ! e3. = dk[gdep]
[7189]	312	!
[6951]	313	DO jk = 1, jpkm1
	314	e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)
	315	END DO
	316	e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO
	317
	318	DO jk = 2, jpk
	319	e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)
	320	END DO
	321	e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))
	322	END IF
	323
	324	!!gm BUG in s-coordinate this does not work!
	325	! deepest/shallowest W level Above/Below ~10m
	326	zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
	327	nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
	328	nla10 = nlb10 - 1 ! deepest W level Above ~10m
	329	!!gm end bug
	330
	331	IF(lwp) THEN ! control print
	332	WRITE(numout,*)
	333	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
	334	WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
	335	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
	336	ENDIF
	337	DO jk = 1, jpk ! control positivity
	338	IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 ' )
	339	IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
	340	END DO
	341	!
	342	END SUBROUTINE zgr_z
	343
	344
	345	SUBROUTINE zgr_bat
	346	!!----------------------------------------------------------------------
	347	!! * ROUTINE zgr_bat *
	348	!!
	349	!! ** Purpose : set bathymetry both in levels and meters
	350	!!
	351	!! ** Method : read or define mbathy and bathy arrays
	352	!! * level bathymetry:
	353	!! The ocean basin geometry is given by a two-dimensional array,
	354	!! mbathy, which is defined as follow :
	355	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
	356	!! at t-point (ji,jj).
	357	!! = 0 over the continental t-point.
	358	!! The array mbathy is checked to verified its consistency with
	359	!! model option. in particular:
	360	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	361	!! along closed boundary.
	362	!! mbathy must be cyclic IF jperio=1.
	363	!! mbathy must be lower or equal to jpk-1.
	364	!! isolated ocean grid points are suppressed from mbathy
	365	!! since they are only connected to remaining
	366	!! ocean through vertical diffusion.
	367	!! ntopo=-1 : rectangular channel or bassin with a bump
	368	!! ntopo= 0 : flat rectangular channel or basin
	369	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
	370	!! bathy is read in 'bathy_meter.nc' NetCDF file
	371	!!
	372	!! ** Action : - mbathy: level bathymetry (in level index)
	373	!! - bathy : meter bathymetry (in meters)
	374	!!----------------------------------------------------------------------
	375	INTEGER :: ji, jj, jk ! dummy loop indices
	376	INTEGER :: inum ! temporary logical unit
	377	INTEGER :: ierror ! error flag
	378	INTEGER :: ii_bump, ij_bump, ih ! bump center position
	379	INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices
	380	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
	381	REAL(wp) :: zi, zj, zh, zhmin ! local scalars
	382	INTEGER , ALLOCATABLE, DIMENSION(:,:) :: idta ! global domain integer data
	383	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta ! global domain scalar data
	384	!!----------------------------------------------------------------------
	385	!
	386	IF(lwp) WRITE(numout,*)
	387	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
	388	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
	389	! ! ================== !
	390	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
	391	! ! ================== !
	392	! ! global domain level and meter bathymetry (idta,zdta)
	393	!
[12101]	394	ALLOCATE( idta(jpiglo,jpjglo), STAT=ierror )
[6951]	395	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate idta array' )
[12101]	396	ALLOCATE( zdta(jpiglo,jpjglo), STAT=ierror )
[6951]	397	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zdta array' )
	398	!
	399	IF( ntopo == 0 ) THEN ! flat basin
	400	IF(lwp) WRITE(numout,*)
	401	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
	402	IF( rn_bathy > 0.01 ) THEN
	403	IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist'
	404	zdta(:,:) = rn_bathy
	405	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
	406	idta(:,:) = jpkm1
	407	ELSE ! z-coordinate (zco or zps): step-like topography
	408	idta(:,:) = jpkm1
	409	DO jk = 1, jpkm1
	410	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
	411	END DO
	412	ENDIF
	413	ELSE
	414	IF(lwp) WRITE(numout,*) ' Depth = depthw(jpkm1)'
	415	idta(:,:) = jpkm1 ! before last level
	416	zdta(:,:) = gdepw_1d(jpk) ! last w-point depth
	417	h_oce = gdepw_1d(jpk)
	418	ENDIF
	419	ELSE ! bump centered in the basin
	420	IF(lwp) WRITE(numout,*)
	421	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
[12101]	422	ii_bump = jpiglo / 2 ! i-index of the bump center
	423	ij_bump = jpjglo / 2 ! j-index of the bump center
[6951]	424	r_bump = 50000._wp ! bump radius (meters)
	425	h_bump = 2700._wp ! bump height (meters)
	426	h_oce = gdepw_1d(jpk) ! background ocean depth (meters)
	427	IF(lwp) WRITE(numout,*) ' bump characteristics: '
	428	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
	429	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
	430	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
	431	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
	432	!
[12101]	433	DO jj = 1, jpjglo ! zdta :
	434	DO ji = 1, jpiglo
[6951]	435	zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
	436	zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
	437	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
	438	END DO
	439	END DO
	440	! ! idta :
	441	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
	442	idta(:,:) = jpkm1
	443	ELSE ! z-coordinate (zco or zps): step-like topography
	444	idta(:,:) = jpkm1
	445	DO jk = 1, jpkm1
	446	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
	447	END DO
	448	ENDIF
	449	ENDIF
[12101]	450	!
[6951]	451	! ! set GLOBAL boundary conditions
	452	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
	453	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp
[12101]	454	idta( : ,jpjglo) = 0 ; zdta( : ,jpjglo) = 0._wp
[6951]	455	ELSEIF( jperio == 2 ) THEN
	456	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
[12101]	457	idta( : ,jpjglo) = 0 ; zdta( : ,jpjglo) = 0._wp
[6951]	458	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp
[12101]	459	idta(jpiglo, : ) = 0 ; zdta(jpiglo, : ) = 0._wp
[6951]	460	ELSE
	461	ih = 0 ; zh = 0._wp
	462	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
	463	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
[12101]	464	idta( : ,jpjglo) = ih ; zdta( : ,jpjglo) = zh
[6951]	465	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
[12101]	466	idta(jpiglo, : ) = ih ; zdta(jpiglo, : ) = zh
[6951]	467	ENDIF
	468
	469	! ! local domain level and meter bathymetries (mbathy,bathy)
	470	mbathy(:,:) = 0 ! set to zero extra halo points
	471	bathy (:,:) = 0._wp ! (require for mpp case)
	472	DO jj = 1, nlcj ! interior values
	473	DO ji = 1, nlci
	474	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
	475	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
	476	END DO
	477	END DO
	478	risfdep(:,:)=0.e0
	479	misfdep(:,:)=1
	480	!
	481	DEALLOCATE( idta, zdta )
	482	!
	483	! ! ================ !
[10727]	484	ELSEIF( ntopo == 1 .OR. ntopo ==2 ) THEN ! read in file ! (over the local domain)
[6951]	485	! ! ================ !
	486	!
	487	IF( ln_zco ) THEN ! zco : read level bathymetry
	488	CALL iom_open ( 'bathy_level.nc', inum )
	489	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
	490	CALL iom_close( inum )
	491	mbathy(:,:) = INT( bathy(:,:) )
	492	! initialisation isf variables
	493	risfdep(:,:)=0._wp ; misfdep(:,:)=1
	494	! ! =====================
	495	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	496	! ! =====================
	497	!
	498	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	499	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
	500	DO ji = mi0(ii0), mi1(ii1)
	501	DO jj = mj0(ij0), mj1(ij1)
	502	mbathy(ji,jj) = 15
	503	END DO
	504	END DO
	505	IF(lwp) WRITE(numout,*)
	506	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
	507	!
	508	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	509	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
	510	DO ji = mi0(ii0), mi1(ii1)
	511	DO jj = mj0(ij0), mj1(ij1)
	512	mbathy(ji,jj) = 12
	513	END DO
	514	END DO
	515	IF(lwp) WRITE(numout,*)
	516	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
	517	!
	518	ENDIF
	519	!
	520	ENDIF
	521	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
[10727]	522	#if defined key_agrif
	523	IF (agrif_root()) THEN
	524	#endif
	525	IF( ntopo == 1) THEN
	526	CALL iom_open ( 'bathy_meter.nc', inum )
	527	IF ( ln_isfcav ) THEN
	528	CALL iom_get ( inum, jpdom_data, 'Bathymetry_isf', bathy, lrowattr=.false. )
	529	ELSE
	530	CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr )
	531	END IF
	532	CALL iom_close( inum )
[6951]	533	ELSE
[10727]	534	CALL dom_bat
	535	ENDIF
	536	#if defined key_agrif
	537	ELSE
	538	IF( ntopo == 1) THEN
	539	CALL agrif_create_bathy_meter()
	540	ELSE
	541	CALL dom_bat
	542	ENDIF
	543	ENDIF
	544	#endif
[6951]	545	!
	546	! initialisation isf variables
	547	risfdep(:,:)=0._wp ; misfdep(:,:)=1
	548	!
	549	IF ( ln_isfcav ) THEN
	550	CALL iom_open ( 'isf_draft_meter.nc', inum )
	551	CALL iom_get ( inum, jpdom_data, 'isf_draft', risfdep )
	552	CALL iom_close( inum )
	553	END IF
	554	!
	555	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	556	!
	557	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	558	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
	559	DO ji = mi0(ii0), mi1(ii1)
	560	DO jj = mj0(ij0), mj1(ij1)
	561	bathy(ji,jj) = 284._wp
	562	END DO
	563	END DO
	564	IF(lwp) WRITE(numout,*)
	565	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
	566	!
	567	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	568	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
	569	DO ji = mi0(ii0), mi1(ii1)
	570	DO jj = mj0(ij0), mj1(ij1)
	571	bathy(ji,jj) = 137._wp
	572	END DO
	573	END DO
	574	IF(lwp) WRITE(numout,*)
	575	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
	576	!
	577	ENDIF
	578	!
	579	ENDIF
	580	! ! =============== !
	581	ELSE ! error !
	582	! ! =============== !
	583	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
	584	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
	585	ENDIF
	586	!
	587	IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
	588	IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
	589	ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth
	590	ENDIF
[12101]	591	zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels
[6951]	592	WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
[12101]	593	ELSE WHERE ( risfdep == 0._wp ); bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
[6951]	594	END WHERE
	595	IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
	596	ENDIF
	597	!
	598	END SUBROUTINE zgr_bat
	599
	600	SUBROUTINE zgr_bat_ctl
	601	!!----------------------------------------------------------------------
	602	!! * ROUTINE zgr_bat_ctl *
	603	!!
	604	!! ** Purpose : check the bathymetry in levels
	605	!!
	606	!! ** Method : The array mbathy is checked to verified its consistency
	607	!! with the model options. in particular:
	608	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	609	!! along closed boundary.
	610	!! mbathy must be cyclic IF jperio=1.
	611	!! mbathy must be lower or equal to jpk-1.
	612	!! isolated ocean grid points are suppressed from mbathy
	613	!! since they are only connected to remaining
	614	!! ocean through vertical diffusion.
	615	!! C A U T I O N : mbathy will be modified during the initializa-
	616	!! tion phase to become the number of non-zero w-levels of a water
	617	!! column, with a minimum value of 1.
	618	!!
	619	!! ** Action : - update mbathy: level bathymetry (in level index)
	620	!! - update bathy : meter bathymetry (in meters)
	621	!!----------------------------------------------------------------------
	622	INTEGER :: ji, jj, jl ! dummy loop indices
	623	INTEGER :: icompt, ibtest, ikmax ! temporary integers
[12101]	624	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zbathy
[6951]	625	!!----------------------------------------------------------------------
	626	!
[12101]	627	ALLOCATE(zbathy(jpi,jpj))
[6951]	628	!
	629	IF(lwp) WRITE(numout,*)
	630	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
	631	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
	632	! ! Suppress isolated ocean grid points
	633	IF(lwp) WRITE(numout,*)
	634	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
	635	IF(lwp) WRITE(numout,*)' -----------------------------------'
	636	icompt = 0
	637	DO jl = 1, 2
[12101]	638	IF( l_Iperio ) THEN
[6951]	639	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
	640	mbathy(jpi,:) = mbathy( 2 ,:)
	641	ENDIF
[12101]	642	zbathy(:,:) = FLOAT( mbathy(:,:) )
	643	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
	644	mbathy(:,:) = INT( zbathy(:,:) )
	645
[6951]	646	DO jj = 2, jpjm1
	647	DO ji = 2, jpim1
	648	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
	649	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
	650	IF( ibtest < mbathy(ji,jj) ) THEN
	651	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
	652	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
	653	mbathy(ji,jj) = ibtest
	654	icompt = icompt + 1
	655	ENDIF
	656	END DO
	657	END DO
	658	END DO
[12101]	659
	660	IF( lk_mpp ) CALL mpp_sum( 'domzgr', icompt )
[6951]	661	IF( icompt == 0 ) THEN
	662	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
	663	ELSE
	664	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
	665	ENDIF
[12101]	666
	667	zbathy(:,:) = FLOAT( mbathy(:,:) )
	668	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
	669	mbathy(:,:) = INT( zbathy(:,:) )
	670
[6951]	671	! ! East-west cyclic boundary conditions
[12101]	672	IF( jperio == 0 ) THEN
	673	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: jperio = ', jperio
[6951]	674	IF( lk_mpp ) THEN
	675	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
	676	IF( jperio /= 1 ) mbathy(1,:) = 0
	677	ENDIF
	678	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
	679	IF( jperio /= 1 ) mbathy(nlci,:) = 0
	680	ENDIF
	681	ELSE
	682	IF( ln_zco .OR. ln_zps ) THEN
	683	mbathy( 1 ,:) = 0
	684	mbathy(jpi,:) = 0
	685	ELSE
	686	mbathy( 1 ,:) = jpkm1
	687	mbathy(jpi,:) = jpkm1
	688	ENDIF
	689	ENDIF
[12101]	690	ELSEIF( l_Iperio ) THEN
	691	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio
[6951]	692	mbathy( 1 ,:) = mbathy(jpim1,:)
	693	mbathy(jpi,:) = mbathy( 2 ,:)
[12101]	694	ELSEIF( jperio == 2 ) THEN
	695	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: jperio = ', jperio
[6951]	696	ELSE
	697	IF(lwp) WRITE(numout,*) ' e r r o r'
[12101]	698	IF(lwp) WRITE(numout,*) ' parameter , jperio = ', jperio
[6951]	699	! STOP 'dom_mba'
	700	ENDIF
[12101]	701
[6951]	702	! Boundary condition on mbathy
	703	IF( .NOT.lk_mpp ) THEN
	704	!!gm !!bug ??? think about it !
	705	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
	706	zbathy(:,:) = FLOAT( mbathy(:,:) )
[12101]	707	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
[6951]	708	mbathy(:,:) = INT( zbathy(:,:) )
	709	ENDIF
[12101]	710
[6951]	711	! Number of ocean level inferior or equal to jpkm1
[12101]	712	zbathy(:,:) = FLOAT( mbathy(:,:) )
	713	ikmax = glob_max( 'domzgr', zbathy(:,:) )
	714
[6951]	715	IF( ikmax > jpkm1 ) THEN
	716	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
	717	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
	718	ELSE IF( ikmax < jpkm1 ) THEN
	719	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
	720	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
	721	ENDIF
	722	!
[12101]	723	DEALLOCATE( zbathy )
[6951]	724	!
	725	END SUBROUTINE zgr_bat_ctl
	726
	727
	728	SUBROUTINE zgr_bot_level
	729	!!----------------------------------------------------------------------
	730	!! * ROUTINE zgr_bot_level *
	731	!!
	732	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
	733	!!
	734	!! ** Method : computes from mbathy with a minimum value of 1 over land
	735	!!
	736	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
	737	!! ocean level at t-, u- & v-points
	738	!! (min value = 1 over land)
	739	!!----------------------------------------------------------------------
	740	INTEGER :: ji, jj ! dummy loop indices
[12101]	741	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmbk
[6951]	742	!!----------------------------------------------------------------------
	743	!
[12101]	744	ALLOCATE( zmbk(jpi,jpj) )
[6951]	745	!
	746	IF(lwp) WRITE(numout,*)
	747	IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
	748	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
	749	!
	750	mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
	751
	752	! ! bottom k-index of W-level = mbkt+1
	753	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
	754	DO ji = 1, jpim1
	755	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
	756	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
	757	END DO
	758	END DO
	759	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
[12101]	760	zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
	761	zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
[6951]	762	!
[12101]	763	DEALLOCATE( zmbk )
[6951]	764	!
	765	END SUBROUTINE zgr_bot_level
	766
	767
	768	SUBROUTINE zgr_top_level
	769	!!----------------------------------------------------------------------
	770	!! * ROUTINE zgr_top_level *
	771	!!
	772	!! ** Purpose : defines the vertical index of ocean top (mik. arrays)
	773	!!
	774	!! ** Method : computes from misfdep with a minimum value of 1
	775	!!
	776	!! ** Action : mikt, miku, mikv : vertical indices of the shallowest
	777	!! ocean level at t-, u- & v-points
	778	!! (min value = 1)
	779	!!----------------------------------------------------------------------
	780	INTEGER :: ji, jj ! dummy loop indices
[12101]	781	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmik
[6951]	782	!!----------------------------------------------------------------------
	783	!
[12101]	784	ALLOCATE( zmik(jpi,jpj) )
[6951]	785	!
	786	IF(lwp) WRITE(numout,*)
	787	IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
	788	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
	789	!
	790	mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1)
	791	! ! top k-index of W-level (=mikt)
	792	DO jj = 1, jpjm1 ! top k-index of U- (U-) level
	793	DO ji = 1, jpim1
	794	miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) )
	795	mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) )
	796	mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) )
	797	END DO
	798	END DO
	799
	800	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
[12101]	801	zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 )
	802	zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 )
	803	zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 )
[6951]	804	!
[12101]	805	DEALLOCATE( zmik )
[6951]	806	!
	807	END SUBROUTINE zgr_top_level
	808
	809
	810	SUBROUTINE zgr_zco
	811	!!----------------------------------------------------------------------
	812	!! * ROUTINE zgr_zco *
	813	!!
	814	!! ** Purpose : define the reference z-coordinate system
	815	!!
	816	!! ** Method : set 3D coord. arrays to reference 1D array
	817	!!----------------------------------------------------------------------
	818	INTEGER :: jk
	819	!!----------------------------------------------------------------------
	820	!
	821	DO jk = 1, jpk
	822	gdept_0(:,:,jk) = gdept_1d(jk)
	823	gdepw_0(:,:,jk) = gdepw_1d(jk)
	824	e3t_0 (:,:,jk) = e3t_1d (jk)
	825	e3u_0 (:,:,jk) = e3t_1d (jk)
	826	e3v_0 (:,:,jk) = e3t_1d (jk)
	827	e3f_0 (:,:,jk) = e3t_1d (jk)
	828	e3w_0 (:,:,jk) = e3w_1d (jk)
	829	e3uw_0 (:,:,jk) = e3w_1d (jk)
	830	e3vw_0 (:,:,jk) = e3w_1d (jk)
	831	END DO
	832	!
	833	END SUBROUTINE zgr_zco
	834
	835
	836	SUBROUTINE zgr_zps
	837	!!----------------------------------------------------------------------
	838	!! * ROUTINE zgr_zps *
	839	!!
	840	!! ** Purpose : the depth and vertical scale factor in partial step
	841	!! reference z-coordinate case
	842	!!
	843	!! ** Method : Partial steps : computes the 3D vertical scale factors
	844	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
	845	!! a partial step representation of bottom topography.
	846	!!
	847	!! The reference depth of model levels is defined from an analytical
	848	!! function the derivative of which gives the reference vertical
	849	!! scale factors.
	850	!! From depth and scale factors reference, we compute there new value
	851	!! with partial steps on 3d arrays ( i, j, k ).
	852	!!
	853	!! w-level: gdepw_0(i,j,k) = gdep(k)
	854	!! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k)
	855	!! t-level: gdept_0(i,j,k) = gdep(k+0.5)
	856	!! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
	857	!!
	858	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
	859	!! we find the mbathy index of the depth at each grid point.
	860	!! This leads us to three cases:
	861	!!
	862	!! - bathy = 0 => mbathy = 0
	863	!! - 1 < mbathy < jpkm1
	864	!! - bathy > gdepw_0(jpk) => mbathy = jpkm1
	865	!!
	866	!! Then, for each case, we find the new depth at t- and w- levels
	867	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
	868	!! and f-points.
	869	!!
	870	!! This routine is given as an example, it must be modified
	871	!! following the user s desiderata. nevertheless, the output as
	872	!! well as the way to compute the model levels and scale factors
	873	!! must be respected in order to insure second order accuracy
	874	!! schemes.
	875	!!
	876	!! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
	877	!! - - - - - - - gdept_0, gdepw_0 and e3. are positives
	878	!!
	879	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
	880	!!----------------------------------------------------------------------
	881	INTEGER :: ji, jj, jk ! dummy loop indices
	882	INTEGER :: ik, it, ikb, ikt ! temporary integers
	883	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
	884	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
	885	REAL(wp) :: zdiff ! temporary scalar
	886	REAL(wp) :: zmax ! temporary scalar
[12101]	887	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zprt
[6951]	888	!!---------------------------------------------------------------------
	889	!
[12101]	890	ALLOCATE( zprt(jpi,jpj,jpk) )
[6951]	891	!
	892	IF(lwp) WRITE(numout,*)
	893	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
	894	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
	895	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
	896
[12101]	897	! compute position of the ice shelf grounding line
	898	! set bathy and isfdraft to 0 where grounded
	899	IF ( ln_isfcav ) CALL zgr_isf_zspace
	900
[6951]	901	! bathymetry in level (from bathy_meter)
	902	! ===================
	903	zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
	904	bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
	905	WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
	906	ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
	907	END WHERE
	908
	909	! Compute mbathy for ocean points (i.e. the number of ocean levels)
	910	! find the number of ocean levels such that the last level thickness
	911	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
	912	! e3t_1d is the reference level thickness
	913	DO jk = jpkm1, 1, -1
	914	zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
	915	WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
	916	END DO
	917
[12101]	918	! Check compatibility between bathy and iceshelf draft
	919	! insure at least 2 wet level on the vertical under an ice shelf
	920	! compute misfdep and adjust isf draft if needed
	921	IF ( ln_isfcav ) CALL zgr_isf_kspace
	922
[6951]	923	! Scale factors and depth at T- and W-points
	924	DO jk = 1, jpk ! intitialization to the reference z-coordinate
	925	gdept_0(:,:,jk) = gdept_1d(jk)
	926	gdepw_0(:,:,jk) = gdepw_1d(jk)
	927	e3t_0 (:,:,jk) = e3t_1d (jk)
	928	e3w_0 (:,:,jk) = e3w_1d (jk)
	929	END DO
	930
	931	! Scale factors and depth at T- and W-points
[12101]	932	DO jj = 1, jpj
	933	DO ji = 1, jpi
	934	ik = mbathy(ji,jj)
	935	IF( ik > 0 ) THEN ! ocean point only
	936	! max ocean level case
	937	IF( ik == jpkm1 ) THEN
	938	zdepwp = bathy(ji,jj)
	939	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
	940	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
	941	e3t_0(ji,jj,ik ) = ze3tp
	942	e3t_0(ji,jj,ik+1) = ze3tp
	943	e3w_0(ji,jj,ik ) = ze3wp
	944	e3w_0(ji,jj,ik+1) = ze3tp
	945	gdepw_0(ji,jj,ik+1) = zdepwp
	946	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
	947	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
	948	!
	949	ELSE ! standard case
	950	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
	951	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
[6951]	952	ENDIF
[12101]	953	!gm Bug? check the gdepw_1d
	954	! ... on ik
	955	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
	956	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
	957	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
	958	e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
	959	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
	960	e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
	961	& * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
	962	! ... on ik+1
	963	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
	964	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
	965	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
[6951]	966	ENDIF
[12101]	967	ENDIF
[6951]	968	END DO
[12101]	969	END DO
	970	!
	971	it = 0
	972	DO jj = 1, jpj
	973	DO ji = 1, jpi
	974	ik = mbathy(ji,jj)
	975	IF( ik > 0 ) THEN ! ocean point only
	976	e3tp (ji,jj) = e3t_0(ji,jj,ik)
	977	e3wp (ji,jj) = e3w_0(ji,jj,ik)
	978	! test
	979	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
	980	IF( zdiff <= 0._wp .AND. lwp ) THEN
	981	it = it + 1
	982	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
	983	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
	984	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
	985	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
[6951]	986	ENDIF
[12101]	987	ENDIF
[6951]	988	END DO
[12101]	989	END DO
[6951]	990	!
[12101]	991	! compute top scale factor if ice shelf
	992	IF (ln_isfcav) CALL zps_isf
	993	!
[6951]	994	! Scale factors and depth at U-, V-, UW and VW-points
	995	DO jk = 1, jpk ! initialisation to z-scale factors
	996	e3u_0 (:,:,jk) = e3t_1d(jk)
	997	e3v_0 (:,:,jk) = e3t_1d(jk)
	998	e3uw_0(:,:,jk) = e3w_1d(jk)
	999	e3vw_0(:,:,jk) = e3w_1d(jk)
	1000	END DO
	1001
	1002	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
	1003	DO jj = 1, jpjm1
	1004	DO ji = 1, jpim1 ! vector opt.
	1005	e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
	1006	e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
	1007	e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
	1008	e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
	1009	END DO
	1010	END DO
	1011	END DO
	1012
[12101]	1013	! update e3uw in case only 2 cells in the water column
	1014	IF ( ln_isfcav ) CALL zps_isf_e3uv_w
[6951]	1015	!
[12101]	1016	CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp ) ! lateral boundary conditions
	1017	CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
	1018	!
[6951]	1019	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
	1020	WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)
	1021	WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)
	1022	WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)
	1023	WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)
	1024	END DO
	1025
	1026	! Scale factor at F-point
	1027	DO jk = 1, jpk ! initialisation to z-scale factors
	1028	e3f_0(:,:,jk) = e3t_1d(jk)
	1029	END DO
	1030	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
	1031	DO jj = 1, jpjm1
	1032	DO ji = 1, jpim1 ! vector opt.
	1033	e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
	1034	END DO
	1035	END DO
	1036	END DO
[12101]	1037	CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp ) ! Lateral boundary conditions
[6951]	1038	!
	1039	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
	1040	WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)
	1041	END DO
	1042	!!gm bug ? : must be a do loop with mj0,mj1
	1043	!
	1044	e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
	1045	e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)
	1046	e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)
	1047	e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)
	1048	e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)
	1049
	1050	! Control of the sign
	1051	IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )
	1052	IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )
	1053	IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )
	1054	IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )
	1055	!
[12101]	1056	! if in the future gde3w_0 need to be compute, use the function defined in NEMO
	1057	! for now gde3w_0 computation is removed as not an output of domcfg
	1058
	1059	DEALLOCATE( zprt )
[6951]	1060	!
	1061	END SUBROUTINE zgr_zps
	1062
	1063	SUBROUTINE zgr_sco
	1064	!!----------------------------------------------------------------------
	1065	!! * ROUTINE zgr_sco *
	1066	!!
	1067	!! ** Purpose : define the s-coordinate system
	1068	!!
	1069	!! ** Method : s-coordinate
	1070	!! The depth of model levels is defined as the product of an
	1071	!! analytical function by the local bathymetry, while the vertical
	1072	!! scale factors are defined as the product of the first derivative
	1073	!! of the analytical function by the bathymetry.
	1074	!! (this solution save memory as depth and scale factors are not
	1075	!! 3d fields)
	1076	!! - Read bathymetry (in meters) at t-point and compute the
	1077	!! bathymetry at u-, v-, and f-points.
	1078	!! hbatu = mi( hbatt )
	1079	!! hbatv = mj( hbatt )
	1080	!! hbatf = mi( mj( hbatt ) )
	1081	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
	1082	!! function and its derivative given as function.
	1083	!! z_gsigt(k) = fssig (k )
	1084	!! z_gsigw(k) = fssig (k-0.5)
	1085	!! z_esigt(k) = fsdsig(k )
	1086	!! z_esigw(k) = fsdsig(k-0.5)
	1087	!! Three options for stretching are give, and they can be modified
	1088	!! following the users requirements. Nevertheless, the output as
	1089	!! well as the way to compute the model levels and scale factors
	1090	!! must be respected in order to insure second order accuracy
	1091	!! schemes.
	1092	!!
	1093	!! The three methods for stretching available are:
	1094	!!
	1095	!! s_sh94 (Song and Haidvogel 1994)
	1096	!! a sinh/tanh function that allows sigma and stretched sigma
	1097	!!
	1098	!! s_sf12 (Siddorn and Furner 2012?)
	1099	!! allows the maintenance of fixed surface and or
	1100	!! bottom cell resolutions (cf. geopotential coordinates)
	1101	!! within an analytically derived stretched S-coordinate framework.
	1102	!!
	1103	!! s_tanh (Madec et al 1996)
	1104	!! a cosh/tanh function that gives stretched coordinates
	1105	!!
	1106	!!----------------------------------------------------------------------
	1107	INTEGER :: ji, jj, jk, jl ! dummy loop argument
	1108	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
	1109	INTEGER :: ios ! Local integer output status for namelist read
	1110	REAL(wp) :: zrmax, ztaper ! temporary scalars
	1111	REAL(wp) :: zrfact
	1112	!
[12101]	1113	REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
	1114	REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
[6951]	1115	!!
	1116	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
	1117	& rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
	1118	!!----------------------------------------------------------------------
	1119	!
[12101]	1120	ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj) )
[6951]	1121	!
	1122	REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
	1123	READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
	1124	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
	1125
	1126	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
	1127	READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
	1128	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
	1129	IF(lwm) WRITE ( numond, namzgr_sco )
	1130
	1131	IF(lwp) THEN ! control print
	1132	WRITE(numout,*)
	1133	WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
	1134	WRITE(numout,*) '~~~~~~~~~~~'
	1135	WRITE(numout,*) ' Namelist namzgr_sco'
	1136	WRITE(numout,*) ' stretching coeffs '
	1137	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
	1138	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
	1139	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
	1140	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
	1141	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
	1142	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
	1143	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
	1144	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
	1145	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
	1146	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
	1147	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
	1148	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
	1149	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
	1150	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
	1151	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
	1152	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
	1153	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
	1154	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
	1155	ENDIF
	1156
	1157	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
	1158	hifu(:,:) = rn_sbot_min
	1159	hifv(:,:) = rn_sbot_min
	1160	hiff(:,:) = rn_sbot_min
	1161
	1162	! ! set maximum ocean depth
	1163	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
	1164
	1165	DO jj = 1, jpj
	1166	DO ji = 1, jpi
	1167	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
	1168	END DO
	1169	END DO
	1170	! ! =============================
	1171	! ! Define the envelop bathymetry (hbatt)
	1172	! ! =============================
	1173	! use r-value to create hybrid coordinates
	1174	zenv(:,:) = bathy(:,:)
	1175	!
	1176	! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
	1177	DO jj = 1, jpj
	1178	DO ji = 1, jpi
	1179	IF( bathy(ji,jj) == 0._wp ) THEN
	1180	iip1 = MIN( ji+1, jpi )
	1181	ijp1 = MIN( jj+1, jpj )
	1182	iim1 = MAX( ji-1, 1 )
	1183	ijm1 = MAX( jj-1, 1 )
	1184	!!gm BUG fix see ticket #1617
	1185	IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) &
	1186	& + bathy(iim1,jj ) + bathy(iip1,jj ) &
	1187	& + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) &
	1188	& zenv(ji,jj) = rn_sbot_min
	1189	!!gm
	1190	!!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + &
	1191	!!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
	1192	!!gm zenv(ji,jj) = rn_sbot_min
	1193	!!gm ENDIF
	1194	!!gm end
	1195	ENDIF
	1196	END DO
	1197	END DO
	1198
	1199	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
[12101]	1200	CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' )
[6951]	1201	!
	1202	! smooth the bathymetry (if required)
	1203	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
	1204	scobot(:,:) = bathy(:,:) ! ocean bottom depth
	1205	!
	1206	jl = 0
	1207	zrmax = 1._wp
	1208	!
	1209	!
	1210	! set scaling factor used in reducing vertical gradients
	1211	zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
	1212	!
	1213	! initialise temporary evelope depth arrays
	1214	ztmpi1(:,:) = zenv(:,:)
	1215	ztmpi2(:,:) = zenv(:,:)
	1216	ztmpj1(:,:) = zenv(:,:)
	1217	ztmpj2(:,:) = zenv(:,:)
	1218	!
	1219	! initialise temporary r-value arrays
	1220	zri(:,:) = 1._wp
	1221	zrj(:,:) = 1._wp
	1222	! ! ================ !
	1223	DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop !
	1224	! ! ================ !
	1225	jl = jl + 1
	1226	zrmax = 0._wp
	1227	! we set zrmax from previous r-values (zri and zrj) first
	1228	! if set after current r-value calculation (as previously)
	1229	! we could exit DO WHILE prematurely before checking r-value
	1230	! of current zenv
	1231	DO jj = 1, nlcj
	1232	DO ji = 1, nlci
	1233	zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
	1234	END DO
	1235	END DO
	1236	zri(:,:) = 0._wp
	1237	zrj(:,:) = 0._wp
	1238	DO jj = 1, nlcj
	1239	DO ji = 1, nlci
	1240	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
	1241	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
	1242	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
	1243	zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
	1244	END IF
	1245	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
	1246	zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
	1247	END IF
	1248	IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact
	1249	IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact
	1250	IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact
	1251	IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact
	1252	END DO
	1253	END DO
[10727]	1254	! IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
[6951]	1255	!
	1256	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax
	1257	!
	1258	DO jj = 1, nlcj
	1259	DO ji = 1, nlci
	1260	zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
	1261	END DO
	1262	END DO
	1263	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
[12101]	1264	CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' )
[6951]	1265	! ! ================ !
	1266	END DO ! End loop !
	1267	! ! ================ !
	1268	DO jj = 1, jpj
	1269	DO ji = 1, jpi
	1270	zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
	1271	END DO
	1272	END DO
	1273	!
	1274	! Envelope bathymetry saved in hbatt
	1275	hbatt(:,:) = zenv(:,:)
	1276	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
	1277	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
	1278	DO jj = 1, jpj
	1279	DO ji = 1, jpi
	1280	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
	1281	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
	1282	END DO
	1283	END DO
	1284	ENDIF
	1285	!
	1286	! ! ==============================
	1287	! ! hbatu, hbatv, hbatf fields
	1288	! ! ==============================
	1289	IF(lwp) THEN
	1290	WRITE(numout,*)
	1291	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
	1292	ENDIF
	1293	hbatu(:,:) = rn_sbot_min
	1294	hbatv(:,:) = rn_sbot_min
	1295	hbatf(:,:) = rn_sbot_min
	1296	DO jj = 1, jpjm1
	1297	DO ji = 1, jpim1 ! NO vector opt.
	1298	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
	1299	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
	1300	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
	1301	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
	1302	END DO
	1303	END DO
	1304
	1305	!
	1306	! Apply lateral boundary condition
	1307	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
[12101]	1308	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp )
[6951]	1309	DO jj = 1, jpj
	1310	DO ji = 1, jpi
	1311	IF( hbatu(ji,jj) == 0._wp ) THEN
	1312	!No worries about the following line when ln_wd == .true.
	1313	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
	1314	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
	1315	ENDIF
	1316	END DO
	1317	END DO
[12101]	1318	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp )
[6951]	1319	DO jj = 1, jpj
	1320	DO ji = 1, jpi
	1321	IF( hbatv(ji,jj) == 0._wp ) THEN
	1322	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
	1323	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
	1324	ENDIF
	1325	END DO
	1326	END DO
[12101]	1327	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp )
[6951]	1328	DO jj = 1, jpj
	1329	DO ji = 1, jpi
	1330	IF( hbatf(ji,jj) == 0._wp ) THEN
	1331	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
	1332	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
	1333	ENDIF
	1334	END DO
	1335	END DO
	1336
	1337	!!bug: key_helsinki a verifer
	1338	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
	1339	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
	1340	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
	1341	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
	1342
	1343	IF( nprint == 1 .AND. lwp ) THEN
	1344	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
	1345	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
	1346	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
	1347	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
	1348	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
	1349	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
	1350	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
	1351	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
	1352	ENDIF
	1353	!! helsinki
	1354
	1355	! ! =======================
	1356	! ! s-ccordinate fields (gdep., e3.)
	1357	! ! =======================
	1358	!
	1359	! non-dimensional "sigma" for model level depth at w- and t-levels
	1360
	1361
	1362	!========================================================================
	1363	! Song and Haidvogel 1994 (ln_s_sh94=T)
	1364	! Siddorn and Furner 2012 (ln_sf12=T)
	1365	! or tanh function (both false)
	1366	!========================================================================
	1367	IF ( ln_s_sh94 ) THEN
	1368	CALL s_sh94()
	1369	ELSE IF ( ln_s_sf12 ) THEN
	1370	CALL s_sf12()
	1371	ELSE
	1372	CALL s_tanh()
	1373	ENDIF
	1374
[12101]	1375	CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp )
	1376	CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp )
	1377	CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp )
	1378	CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp )
	1379	CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp )
	1380	CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp )
	1381	CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
[6951]	1382	!
	1383	WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp
	1384	WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp
	1385	WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp
	1386	WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp
	1387	WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp
	1388	WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp
	1389	WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp
	1390	!!
	1391	! HYBRID :
	1392	DO jj = 1, jpj
	1393	DO ji = 1, jpi
	1394	DO jk = 1, jpkm1
[12101]	1395	IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
[6951]	1396	END DO
	1397	END DO
	1398	END DO
	1399	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
	1400	& ' MAX ', MAXVAL( mbathy(:,:) )
	1401
	1402	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
	1403	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	1404	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
[12101]	1405	& ' w ', MINVAL( gdepw_0(:,:,:) )
[6951]	1406	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
	1407	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
	1408	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
	1409	& ' w ', MINVAL( e3w_0 (:,:,:) )
	1410
	1411	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
[12101]	1412	& ' w ', MAXVAL( gdepw_0(:,:,:) )
[6951]	1413	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
	1414	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
	1415	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
	1416	& ' w ', MAXVAL( e3w_0 (:,:,:) )
	1417	ENDIF
	1418	! END DO
	1419	IF(lwp) THEN ! selected vertical profiles
	1420	WRITE(numout,*)
	1421	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
	1422	WRITE(numout,*) ' ~~~~~~ --------------------'
	1423	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
	1424	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
	1425	& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
	1426	DO jj = mj0(20), mj1(20)
	1427	DO ji = mi0(20), mi1(20)
	1428	WRITE(numout,*)
	1429	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	1430	WRITE(numout,*) ' ~~~~~~ --------------------'
	1431	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
	1432	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
	1433	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
	1434	END DO
	1435	END DO
	1436	DO jj = mj0(74), mj1(74)
	1437	DO ji = mi0(100), mi1(100)
	1438	WRITE(numout,*)
	1439	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	1440	WRITE(numout,*) ' ~~~~~~ --------------------'
	1441	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
	1442	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
	1443	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
	1444	END DO
	1445	END DO
	1446	ENDIF
	1447	!
	1448	!================================================================================
	1449	! check the coordinate makes sense
	1450	!================================================================================
	1451	DO ji = 1, jpi
	1452	DO jj = 1, jpj
	1453	!
	1454	IF( hbatt(ji,jj) > 0._wp) THEN
	1455	DO jk = 1, mbathy(ji,jj)
	1456	! check coordinate is monotonically increasing
[12101]	1457	IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
[6951]	1458	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
	1459	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
[12101]	1460	WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
	1461	WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
[6951]	1462	CALL ctl_stop( ctmp1 )
	1463	ENDIF
	1464	! and check it has never gone negative
[12101]	1465	IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
[6951]	1466	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
	1467	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
[12101]	1468	WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
	1469	WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
[6951]	1470	CALL ctl_stop( ctmp1 )
	1471	ENDIF
	1472	! and check it never exceeds the total depth
[12101]	1473	IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
[6951]	1474	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
	1475	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
[12101]	1476	WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
[6951]	1477	CALL ctl_stop( ctmp1 )
	1478	ENDIF
	1479	END DO
	1480	!
	1481	DO jk = 1, mbathy(ji,jj)-1
	1482	! and check it never exceeds the total depth
[12101]	1483	IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
[6951]	1484	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
	1485	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
[12101]	1486	WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
[6951]	1487	CALL ctl_stop( ctmp1 )
	1488	ENDIF
	1489	END DO
	1490	ENDIF
	1491	END DO
	1492	END DO
	1493	!
[12101]	1494	DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
[6951]	1495	!
	1496	END SUBROUTINE zgr_sco
	1497
	1498
	1499	SUBROUTINE s_sh94()
	1500	!!----------------------------------------------------------------------
	1501	!! * ROUTINE s_sh94 *
	1502	!!
	1503	!! ** Purpose : stretch the s-coordinate system
	1504	!!
	1505	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
	1506	!! mixed S/sigma coordinate
	1507	!!
	1508	!! Reference : Song and Haidvogel 1994.
	1509	!!----------------------------------------------------------------------
	1510	INTEGER :: ji, jj, jk ! dummy loop argument
	1511	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
	1512	REAL(wp) :: ztmpu, ztmpv, ztmpf
	1513	REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
	1514	!
[12101]	1515	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
	1516	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
[6951]	1517	!!----------------------------------------------------------------------
	1518
[12101]	1519	ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
	1520	ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) )
	1521	ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
[6951]	1522
[12101]	1523	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp
[6951]	1524	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
	1525	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
	1526	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
	1527	!
	1528	DO ji = 1, jpi
	1529	DO jj = 1, jpj
	1530	!
	1531	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
	1532	DO jk = 1, jpk
	1533	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
	1534	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
	1535	END DO
	1536	ELSE ! shallow water, uniform sigma
	1537	DO jk = 1, jpk
	1538	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
	1539	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
	1540	END DO
	1541	ENDIF
	1542	!
	1543	DO jk = 1, jpkm1
	1544	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
	1545	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
	1546	END DO
	1547	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
	1548	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
	1549	!
	1550	DO jk = 1, jpk
	1551	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
	1552	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
	1553	gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
	1554	gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
	1555	END DO
	1556	!
	1557	END DO ! for all jj's
	1558	END DO ! for all ji's
	1559
	1560	DO ji = 1, jpim1
	1561	DO jj = 1, jpjm1
	1562	! extended for Wetting/Drying case
	1563	ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
	1564	ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
	1565	ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
	1566	ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
	1567	ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
	1568	ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
	1569	& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
	1570	DO jk = 1, jpk
	1571	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
	1572	& / ztmpu
	1573	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
	1574	& / ztmpu
	1575
	1576	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
	1577	& / ztmpv
	1578	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
	1579	& / ztmpv
	1580
	1581	z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
	1582	& + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
	1583	& + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
	1584	& + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
	1585
	1586	!
	1587	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1588	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1589	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1590	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1591	!
	1592	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1593	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1594	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
	1595	END DO
	1596	END DO
	1597	END DO
	1598	!
[12101]	1599	DEALLOCATE( z_gsigw3, z_gsigt3 )
	1600	DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
[6951]	1601	!
	1602	END SUBROUTINE s_sh94
	1603
	1604
	1605	SUBROUTINE s_sf12
	1606	!!----------------------------------------------------------------------
	1607	!! * ROUTINE s_sf12 *
	1608	!!
	1609	!! ** Purpose : stretch the s-coordinate system
	1610	!!
	1611	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
	1612	!! mixed S/sigma/Z coordinate
	1613	!!
	1614	!! This method allows the maintenance of fixed surface and or
	1615	!! bottom cell resolutions (cf. geopotential coordinates)
	1616	!! within an analytically derived stretched S-coordinate framework.
	1617	!!
	1618	!!
	1619	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
	1620	!!----------------------------------------------------------------------
	1621	INTEGER :: ji, jj, jk ! dummy loop argument
	1622	REAL(wp) :: zsmth ! smoothing around critical depth
	1623	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
	1624	REAL(wp) :: ztmpu, ztmpv, ztmpf
	1625	REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
	1626	!
[12101]	1627	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
	1628	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
[6951]	1629	!!----------------------------------------------------------------------
	1630	!
[12101]	1631	ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
	1632	ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk))
	1633	ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
[6951]	1634
[12101]	1635	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp
[6951]	1636	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
	1637	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
	1638	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
	1639
	1640	DO ji = 1, jpi
	1641	DO jj = 1, jpj
	1642
	1643	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
	1644
	1645	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
	1646	! could be changed by users but care must be taken to do so carefully
	1647	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
	1648
	1649	zzs = rn_zs / hbatt(ji,jj)
	1650
	1651	IF (rn_efold /= 0.0_wp) THEN
	1652	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
	1653	ELSE
	1654	zsmth = 1.0_wp
	1655	ENDIF
	1656
	1657	DO jk = 1, jpk
	1658	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
	1659	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
	1660	ENDDO
	1661	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
	1662	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
	1663
	1664	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
	1665
	1666	DO jk = 1, jpk
	1667	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
	1668	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
	1669	END DO
	1670
	1671	ELSE ! shallow water, z coordinates
	1672
	1673	DO jk = 1, jpk
	1674	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
	1675	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
	1676	END DO
	1677
	1678	ENDIF
	1679
	1680	DO jk = 1, jpkm1
	1681	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
	1682	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
	1683	END DO
	1684	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
	1685	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
	1686
	1687	DO jk = 1, jpk
	1688	gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
	1689	gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
	1690	END DO
	1691
	1692	ENDDO ! for all jj's
	1693	ENDDO ! for all ji's
	1694
	1695	DO ji=1,jpi-1
	1696	DO jj=1,jpj-1
	1697
	1698	! extend to suit for Wetting/Drying case
	1699	ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
	1700	ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
	1701	ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
	1702	ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
	1703	ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
	1704	ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
	1705	& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
	1706	DO jk = 1, jpk
	1707	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
	1708	& / ztmpu
	1709	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
	1710	& / ztmpu
	1711	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
	1712	& / ztmpv
	1713	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
	1714	& / ztmpv
	1715
	1716	z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
	1717	& + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
	1718	& + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
	1719	& + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
	1720
	1721	! Code prior to wetting and drying option (for reference)
	1722	!z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
	1723	! /( hbatt(ji,jj)+hbatt(ji+1,jj) )
	1724	!
	1725	!z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
	1726	! /( hbatt(ji,jj)+hbatt(ji+1,jj) )
	1727	!
	1728	!z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
	1729	! /( hbatt(ji,jj)+hbatt(ji,jj+1) )
	1730	!
	1731	!z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
	1732	! /( hbatt(ji,jj)+hbatt(ji,jj+1) )
	1733	!
	1734	!z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
	1735	! & +hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
	1736	! +hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
	1737	! & +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) &
	1738	! /( hbatt(ji ,jj )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
	1739
	1740	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
	1741	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
	1742	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
	1743	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
	1744	!
	1745	e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
	1746	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
	1747	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
	1748	END DO
	1749
	1750	ENDDO
	1751	ENDDO
	1752	!
[12101]	1753	CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.)
	1754	CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.)
	1755	CALL lbc_lnk('domzgr',e3w_0 ,'T',1.)
	1756	CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.)
[6951]	1757	!
[12101]	1758	DEALLOCATE( z_gsigw3, z_gsigt3 )
	1759	DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
[6951]	1760	!
	1761	END SUBROUTINE s_sf12
	1762
	1763
	1764	SUBROUTINE s_tanh()
	1765	!!----------------------------------------------------------------------
	1766	!! * ROUTINE s_tanh*
	1767	!!
	1768	!! ** Purpose : stretch the s-coordinate system
	1769	!!
	1770	!! ** Method : s-coordinate stretch
	1771	!!
	1772	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
	1773	!!----------------------------------------------------------------------
	1774	INTEGER :: ji, jj, jk ! dummy loop argument
	1775	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
[12101]	1776	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt
	1777	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
[6951]	1778	!!----------------------------------------------------------------------
	1779
[12101]	1780	ALLOCATE( z_gsigw(jpk), z_gsigt(jpk) )
	1781	ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
[6951]	1782
[12101]	1783	z_gsigw = 0._wp ; z_gsigt = 0._wp
[6951]	1784	z_esigt = 0._wp ; z_esigw = 0._wp
	1785
	1786	DO jk = 1, jpk
	1787	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
	1788	z_gsigt(jk) = -fssig( REAL(jk,wp) )
	1789	END DO
	1790	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
	1791	!
	1792	! Coefficients for vertical scale factors at w-, t- levels
	1793	!!gm bug : define it from analytical function, not like juste bellow....
	1794	!!gm or betteroffer the 2 possibilities....
	1795	DO jk = 1, jpkm1
	1796	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
	1797	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
	1798	END DO
	1799	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
	1800	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
	1801	!
	1802	DO jk = 1, jpk
	1803	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
	1804	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
	1805	gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
	1806	gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
	1807	END DO
[12101]	1808
[6951]	1809	DO jj = 1, jpj
	1810	DO ji = 1, jpi
	1811	DO jk = 1, jpk
	1812	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
	1813	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
	1814	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
	1815	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
	1816	!
	1817	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
	1818	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
	1819	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
	1820	END DO
	1821	END DO
	1822	END DO
	1823	!
[12101]	1824	DEALLOCATE( z_gsigw, z_gsigt )
	1825	DEALLOCATE( z_esigt, z_esigw )
[6951]	1826	!
	1827	END SUBROUTINE s_tanh
	1828
	1829
	1830	FUNCTION fssig( pk ) RESULT( pf )
	1831	!!----------------------------------------------------------------------
	1832	!! * ROUTINE fssig *
	1833	!!
	1834	!! ** Purpose : provide the analytical function in s-coordinate
	1835	!!
	1836	!! ** Method : the function provide the non-dimensional position of
	1837	!! T and W (i.e. between 0 and 1)
	1838	!! T-points at integer values (between 1 and jpk)
	1839	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1840	!!----------------------------------------------------------------------
	1841	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
	1842	REAL(wp) :: pf ! sigma value
	1843	!!----------------------------------------------------------------------
	1844	!
	1845	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
	1846	& - TANH( rn_thetb * rn_theta ) ) &
	1847	& * ( COSH( rn_theta ) &
	1848	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
	1849	& / ( 2._wp * SINH( rn_theta ) )
	1850	!
	1851	END FUNCTION fssig
	1852
	1853
	1854	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
	1855	!!----------------------------------------------------------------------
	1856	!! * ROUTINE fssig1 *
	1857	!!
	1858	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
	1859	!!
	1860	!! ** Method : the function provides the non-dimensional position of
	1861	!! T and W (i.e. between 0 and 1)
	1862	!! T-points at integer values (between 1 and jpk)
	1863	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1864	!!----------------------------------------------------------------------
	1865	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
	1866	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
	1867	REAL(wp) :: pf1 ! sigma value
	1868	!!----------------------------------------------------------------------
	1869	!
	1870	IF ( rn_theta == 0 ) then ! uniform sigma
	1871	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
	1872	ELSE ! stretched sigma
	1873	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
	1874	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
	1875	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
	1876	ENDIF
	1877	!
	1878	END FUNCTION fssig1
	1879
	1880
	1881	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
	1882	!!----------------------------------------------------------------------
	1883	!! * ROUTINE fgamma *
	1884	!!
	1885	!! ** Purpose : provide analytical function for the s-coordinate
	1886	!!
	1887	!! ** Method : the function provides the non-dimensional position of
	1888	!! T and W (i.e. between 0 and 1)
	1889	!! T-points at integer values (between 1 and jpk)
	1890	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1891	!!
	1892	!! This method allows the maintenance of fixed surface and or
	1893	!! bottom cell resolutions (cf. geopotential coordinates)
	1894	!! within an analytically derived stretched S-coordinate framework.
	1895	!!
	1896	!! Reference : Siddorn and Furner, in prep
	1897	!!----------------------------------------------------------------------
	1898	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
	1899	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
	1900	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
	1901	REAL(wp), INTENT(in ) :: pzs ! surface box depth
	1902	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
	1903	!
	1904	INTEGER :: jk ! dummy loop index
	1905	REAL(wp) :: za1,za2,za3 ! local scalar
	1906	REAL(wp) :: zn1,zn2 ! - -
	1907	REAL(wp) :: za,zb,zx ! - -
	1908	!!----------------------------------------------------------------------
	1909	!
	1910	zn1 = 1._wp / REAL( jpkm1, wp )
	1911	zn2 = 1._wp - zn1
	1912	!
	1913	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
	1914	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
	1915	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
	1916	!
	1917	za = pzb - za3*(pzs-za1)-za2
	1918	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
	1919	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
	1920	zx = 1.0_wp-za/2.0_wp-zb
	1921	!
	1922	DO jk = 1, jpk
	1923	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
	1924	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
	1925	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
	1926	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
	1927	END DO
	1928	!
	1929	END FUNCTION fgamma
	1930
	1931	!!======================================================================
	1932	END MODULE domzgr

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

Download in other formats: