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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 2392

Last change on this file since 2392 was 2392, checked in by gm, 13 years ago
v3.3beta: Cross Land Advection (ticket #127) full rewriting + MPP bug corrections
Property svn:keywords set to `Id`
File size: 78.9 KB

Rev	Line
[3]	1	MODULE domzgr
	2	!!==============================================================================
	3	!! * MODULE domzgr *
	4	!! Ocean initialization : domain initialization
	5	!!==============================================================================
[1566]	6	!! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate
	7	!! ! 1997-07 (G. Madec) lbc_lnk call
	8	!! ! 1997-04 (J.-O. Beismann)
	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
[1099]	16	!!----------------------------------------------------------------------
[3]	17
	18	!!----------------------------------------------------------------------
[1099]	19	!! dom_zgr : defined the ocean vertical coordinate system
[3]	20	!! zgr_bat : bathymetry fields (levels and meters)
	21	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
	22	!! zgr_bat_ctl : check the bathymetry files
	23	!! zgr_z : reference z-coordinate
[454]	24	!! zgr_zco : z-coordinate
[3]	25	!! zgr_zps : z-coordinate with partial steps
[454]	26	!! zgr_sco : s-coordinate
[1083]	27	!! fssig : sigma coordinate non-dimensional function
[1099]	28	!! dfssig : derivative of the sigma coordinate function !!gm (currently missing!)
[3]	29	!!---------------------------------------------------------------------
	30	USE oce ! ocean dynamics and tracers
	31	USE dom_oce ! ocean space and time domain
	32	USE in_out_manager ! I/O manager
	33	USE lib_mpp ! distributed memory computing library
	34	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[1099]	35	USE closea ! closed seas
[3]	36
	37	IMPLICIT NONE
	38	PRIVATE
	39
[1099]	40	PUBLIC dom_zgr ! called by dom_init.F90
[3]	41
[1099]	42	!!gm DOCTOR name for the namelist parameter...
[1601]	43	! !!! Namelist namzgr_sco
	44	REAL(wp) :: rn_sbot_min = 300. ! minimum depth of s-bottom surface (>0) (m)
	45	REAL(wp) :: rn_sbot_max = 5250. ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
	46	REAL(wp) :: rn_theta = 6.0 ! surface control parameter (0<=rn_theta<=20)
	47	REAL(wp) :: rn_thetb = 0.75 ! bottom control parameter (0<=rn_thetb<= 1)
	48	REAL(wp) :: rn_rmax = 0.15 ! maximum cut-off r-value allowed (0<rn_rmax<1)
	49	LOGICAL :: ln_s_sigma = .false. ! use hybrid s-sigma -coordinate & stretching function fssig1 (ln_sco=T)
	50	REAL(wp) :: rn_bb = 0.8 ! stretching parameter for song and haidvogel stretching
	51	! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
	52	REAL(wp) :: rn_hc = 150. ! Critical depth for s-sigma coordinates
[454]	53
[3]	54	!! * Substitutions
	55	# include "domzgr_substitute.h90"
	56	# include "vectopt_loop_substitute.h90"
	57	!!----------------------------------------------------------------------
[2287]	58	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1146]	59	!! $Id$
[2287]	60	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	61	!!----------------------------------------------------------------------
	62
	63	CONTAINS
	64
	65	SUBROUTINE dom_zgr
	66	!!----------------------------------------------------------------------
	67	!! * ROUTINE dom_zgr *
	68	!!
	69	!! ** Purpose : set the depth of model levels and the resulting
	70	!! vertical scale factors.
	71	!!
[1099]	72	!! ** Method : - reference 1D vertical coordinate (gdep._0, e3._0)
	73	!! - read/set ocean depth and ocean levels (bathy, mbathy)
	74	!! - vertical coordinate (gdep., e3.) depending on the
	75	!! coordinate chosen :
[2240]	76	!! ln_zco=T z-coordinate
[1099]	77	!! ln_zps=T z-coordinate with partial steps
	78	!! ln_zco=T s-coordinate
[3]	79	!!
[1099]	80	!! ** Action : define gdep., e3., mbathy and bathy
	81	!!----------------------------------------------------------------------
	82	INTEGER :: ioptio = 0 ! temporary integer
[3]	83	!!
[1601]	84	NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
[3]	85	!!----------------------------------------------------------------------
	86
[1601]	87	REWIND ( numnam ) ! Read Namelist namzgr : vertical coordinate'
	88	READ ( numnam, namzgr )
[454]	89
[1099]	90	IF(lwp) THEN ! Control print
[454]	91	WRITE(numout,*)
	92	WRITE(numout,*) 'dom_zgr : vertical coordinate'
	93	WRITE(numout,*) '~~~~~~~'
[1601]	94	WRITE(numout,*) ' Namelist namzgr : set vertical coordinate'
[454]	95	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
	96	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
	97	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
	98	ENDIF
	99
[1099]	100	ioptio = 0 ! Check Vertical coordinate options
[454]	101	IF( ln_zco ) ioptio = ioptio + 1
	102	IF( ln_zps ) ioptio = ioptio + 1
	103	IF( ln_sco ) ioptio = ioptio + 1
[1099]	104	IF ( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
[454]	105
[3]	106	! Build the vertical coordinate system
	107	! ------------------------------------
[454]	108	CALL zgr_z ! Reference z-coordinate system (always called)
	109	CALL zgr_bat ! Bathymetry fields (levels and meters)
	110	IF( ln_zco ) CALL zgr_zco ! z-coordinate
	111	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
	112	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
[1348]	113
	114	!!bug gm control print:
	115	IF( nprint == 1 .AND. lwp ) THEN
	116	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	117	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
	118	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
	119	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ), &
	120	& ' u ', MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ), &
	121	& ' uw', MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)), &
	122	& ' w ', MINVAL( fse3w(:,:,:) )
	123
	124	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
	125	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
	126	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ), &
	127	& ' u ', MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ), &
	128	& ' uw', MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)), &
	129	& ' w ', MAXVAL( fse3w(:,:,:) )
	130	ENDIF
	131	!!!bug gm
	132
[3]	133	END SUBROUTINE dom_zgr
	134
	135
	136	SUBROUTINE zgr_z
	137	!!----------------------------------------------------------------------
	138	!! * ROUTINE zgr_z *
	139	!!
	140	!! ** Purpose : set the depth of model levels and the resulting
	141	!! vertical scale factors.
	142	!!
	143	!! ** Method : z-coordinate system (use in all type of coordinate)
	144	!! The depth of model levels is defined from an analytical
	145	!! function the derivative of which gives the scale factors.
	146	!! both depth and scale factors only depend on k (1d arrays).
[454]	147	!! w-level: gdepw_0 = fsdep(k)
	148	!! e3w_0(k) = dk(fsdep)(k) = fse3(k)
	149	!! t-level: gdept_0 = fsdep(k+0.5)
	150	!! e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
[3]	151	!!
[454]	152	!! ** Action : - gdept_0, gdepw_0 : depth of T- and W-point (m)
[1099]	153	!! - e3t_0 , e3w_0 : scale factors at T- and W-levels (m)
[3]	154	!!
[1099]	155	!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
[3]	156	!!----------------------------------------------------------------------
	157	INTEGER :: jk ! dummy loop indices
	158	REAL(wp) :: zt, zw ! temporary scalars
[1099]	159	REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
	160	REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
[1577]	161	REAL(wp) :: zrefdep ! depth of the reference level (~10m)
[2379]	162	REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
[3]	163	!!----------------------------------------------------------------------
	164
	165	! Set variables from parameters
	166	! ------------------------------
	167	zkth = ppkth ; zacr = ppacr
	168	zdzmin = ppdzmin ; zhmax = pphmax
[2379]	169	zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters
[3]	170
	171	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
	172	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
	173	!
[1099]	174	IF( ppa1 == pp_to_be_computed .AND. &
[3]	175	& ppa0 == pp_to_be_computed .AND. &
	176	& ppsur == pp_to_be_computed ) THEN
[1099]	177	!
	178	za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
	179	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
	180	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
	181	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
	182	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
	183	ELSE
[3]	184	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
[2379]	185	za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter
[1099]	186	ENDIF
[3]	187
[1099]	188	IF(lwp) THEN ! Parameter print
[3]	189	WRITE(numout,*)
	190	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
	191	WRITE(numout,*) ' ~~~~~~~'
[454]	192	IF( ppkth == 0. ) THEN
[250]	193	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
	194	WRITE(numout,*) ' Total depth :', zhmax
	195	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
	196	ELSE
[454]	197	IF( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0. ) THEN
[250]	198	WRITE(numout,*) ' zsur, za0, za1 computed from '
	199	WRITE(numout,*) ' zdzmin = ', zdzmin
	200	WRITE(numout,*) ' zhmax = ', zhmax
	201	ENDIF
	202	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
	203	WRITE(numout,*) ' zsur = ', zsur
	204	WRITE(numout,*) ' za0 = ', za0
	205	WRITE(numout,*) ' za1 = ', za1
	206	WRITE(numout,*) ' zkth = ', zkth
	207	WRITE(numout,*) ' zacr = ', zacr
[2379]	208	IF( ldbletanh ) THEN
	209	WRITE(numout,*) ' (Double tanh za2 = ', za2
	210	WRITE(numout,*) ' parameters) zkth2= ', zkth2
	211	WRITE(numout,*) ' zacr2= ', zacr2
	212	ENDIF
[3]	213	ENDIF
	214	ENDIF
	215
	216
	217	! Reference z-coordinate (depth - scale factor at T- and W-points)
	218	! ======================
[454]	219	IF( ppkth == 0.e0 ) THEN ! uniform vertical grid
	220	za1 = zhmax / FLOAT(jpk-1)
[250]	221	DO jk = 1, jpk
	222	zw = FLOAT( jk )
	223	zt = FLOAT( jk ) + 0.5
[454]	224	gdepw_0(jk) = ( zw - 1 ) * za1
	225	gdept_0(jk) = ( zt - 1 ) * za1
	226	e3w_0 (jk) = za1
	227	e3t_0 (jk) = za1
[250]	228	END DO
[1099]	229	ELSE ! Madec & Imbard 1996 function
[2379]	230	IF( .NOT. ldbletanh ) THEN
	231	DO jk = 1, jpk
	232	zw = FLOAT( jk )
	233	zt = FLOAT( jk ) + 0.5
	234	gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
	235	gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
	236	e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
	237	e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
	238	END DO
	239	ELSE
	240	DO jk = 1, jpk
	241	zw = FLOAT( jk )
	242	zt = FLOAT( jk ) + 0.5
	243	! Double tanh function
	244	gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) &
	245	& + za2 * zacr2* LOG ( COSH( (zw-zkth2)/zacr2 ) ) )
	246	gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) &
	247	& + za2 * zacr2* LOG ( COSH( (zt-zkth2)/zacr2 ) ) )
	248	e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth) / zacr ) &
	249	& + za2 * TANH( (zw-zkth2)/zacr2 )
	250	e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth) / zacr ) &
	251	& + za2 * TANH( (zt-zkth2)/zacr2 )
	252	END DO
	253	ENDIF
[1099]	254	gdepw_0(1) = 0.e0 ! force first w-level to be exactly at zero
[250]	255	ENDIF
	256
[1601]	257	!!gm BUG in s-coordinate this does not work!
[2379]	258	! deepest/shallowest W level Above/Below ~10m
[1577]	259	zrefdep = 10. - ( 0.1*MINVAL(e3w_0) ) ! ref. depth with tolerance (10% of minimum layer thickness)
[2379]	260	nlb10 = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
	261	nla10 = nlb10 - 1 ! deepest W level Above ~10m
[1601]	262	!!gm end bug
[1577]	263
[1099]	264	IF(lwp) THEN ! control print
[3]	265	WRITE(numout,*)
	266	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
	267	WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" )
[454]	268	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
[3]	269	ENDIF
[1099]	270	DO jk = 1, jpk ! control positivity
	271	IF( e3w_0 (jk) <= 0.e0 .OR. e3t_0 (jk) <= 0.e0 ) CALL ctl_stop( ' e3w or e3t =< 0 ' )
	272	IF( gdepw_0(jk) < 0.e0 .OR. gdept_0(jk) < 0.e0 ) CALL ctl_stop( ' gdepw or gdept < 0 ' )
[3]	273	END DO
[1099]	274	!
[3]	275	END SUBROUTINE zgr_z
	276
	277
	278	SUBROUTINE zgr_bat
	279	!!----------------------------------------------------------------------
	280	!! * ROUTINE zgr_bat *
	281	!!
	282	!! ** Purpose : set bathymetry both in levels and meters
	283	!!
	284	!! ** Method : read or define mbathy and bathy arrays
	285	!! * level bathymetry:
	286	!! The ocean basin geometry is given by a two-dimensional array,
	287	!! mbathy, which is defined as follow :
	288	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
	289	!! at t-point (ji,jj).
	290	!! = 0 over the continental t-point.
	291	!! The array mbathy is checked to verified its consistency with
	292	!! model option. in particular:
	293	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	294	!! along closed boundary.
	295	!! mbathy must be cyclic IF jperio=1.
	296	!! mbathy must be lower or equal to jpk-1.
	297	!! isolated ocean grid points are suppressed from mbathy
	298	!! since they are only connected to remaining
	299	!! ocean through vertical diffusion.
	300	!! ntopo=-1 : rectangular channel or bassin with a bump
	301	!! ntopo= 0 : flat rectangular channel or basin
[128]	302	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
[3]	303	!! bathy is read in 'bathy_meter.nc' NetCDF file
	304	!! C A U T I O N : mbathy will be modified during the initializa-
	305	!! tion phase to become the number of non-zero w-levels of a water
	306	!! column, with a minimum value of 1.
	307	!!
	308	!! ** Action : - mbathy: level bathymetry (in level index)
	309	!! - bathy : meter bathymetry (in meters)
	310	!!----------------------------------------------------------------------
[473]	311	USE iom
[1099]	312	!!
	313	INTEGER :: ji, jj, jl, jk ! dummy loop indices
	314	INTEGER :: inum ! temporary logical unit
[1348]	315	INTEGER :: ii_bump, ij_bump, ih ! bump center position
[1273]	316	INTEGER :: ii0, ii1, ij0, ij1 ! indices
[1099]	317	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
	318	REAL(wp) :: zi , zj , zh ! temporary scalars
	319	INTEGER , DIMENSION(jpidta,jpjdta) :: idta ! global domain integer data
	320	REAL(wp), DIMENSION(jpidta,jpjdta) :: zdta ! global domain scalar data
[3]	321	!!----------------------------------------------------------------------
	322
	323	IF(lwp) WRITE(numout,*)
	324	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
	325	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
	326
[1099]	327	! ! ================== !
	328	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
	329	! ! ================== !
	330	! ! global domain level and meter bathymetry (idta,zdta)
	331	!
[3]	332	IF( ntopo == 0 ) THEN ! flat basin
	333	IF(lwp) WRITE(numout,*)
	334	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
[1099]	335	idta(:,:) = jpkm1 ! before last level
	336	zdta(:,:) = gdepw_0(jpk) ! last w-point depth
[592]	337	h_oce = gdepw_0(jpk)
[1099]	338	ELSE ! bump centered in the basin
[3]	339	IF(lwp) WRITE(numout,*)
	340	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
[1099]	341	ii_bump = jpidta / 2 ! i-index of the bump center
	342	ij_bump = jpjdta / 2 ! j-index of the bump center
	343	r_bump = 50000.e0 ! bump radius (meters)
	344	h_bump = 2700.e0 ! bump height (meters)
	345	h_oce = gdepw_0(jpk) ! background ocean depth (meters)
[3]	346	IF(lwp) WRITE(numout,*) ' bump characteristics: '
	347	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
	348	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
	349	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
	350	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
[1099]	351	!
	352	DO jj = 1, jpjdta ! zdta :
[3]	353	DO ji = 1, jpidta
[592]	354	zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
	355	zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
[3]	356	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
	357	END DO
	358	END DO
[1099]	359	! ! idta :
	360	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
[454]	361	idta(:,:) = jpkm1
[1099]	362	ELSE ! z-coordinate (zco or zps): step-like topography
[454]	363	idta(:,:) = jpkm1
	364	DO jk = 1, jpkm1
	365	DO jj = 1, jpjdta
	366	DO ji = 1, jpidta
	367	IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) ) idta(ji,jj) = jk
	368	END DO
[3]	369	END DO
	370	END DO
[454]	371	ENDIF
[3]	372	ENDIF
[1099]	373	! ! set GLOBAL boundary conditions
	374	! ! Caution : idta on the global domain: use of jperio, not nperio
[3]	375	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
[30]	376	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1.e0
	377	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0
[3]	378	ELSEIF( jperio == 2 ) THEN
[30]	379	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
	380	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0
	381	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0.e0
	382	idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0.e0
[3]	383	ELSE
[454]	384	ih = 0 ; zh = 0.e0
[525]	385	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
[454]	386	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
	387	idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh
	388	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
	389	idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh
[3]	390	ENDIF
	391
[1099]	392	! ! local domain level and meter bathymetries (mbathy,bathy)
	393	mbathy(:,:) = 0 ! set to zero extra halo points
	394	bathy (:,:) = 0.e0 ! (require for mpp case)
	395	DO jj = 1, nlcj ! interior values
[473]	396	DO ji = 1, nlci
	397	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
	398	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
	399	END DO
	400	END DO
[1099]	401	!
	402	! ! ================ !
	403	ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain)
	404	! ! ================ !
	405	!
	406	IF( ln_zco ) THEN ! zco : read level bathymetry
	407	CALL iom_open( 'bathy_level.nc', inum )
[473]	408	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
	409	CALL iom_close (inum)
	410	mbathy(:,:) = INT( bathy(:,:) )
[1273]	411	! ! =====================
	412	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	413	! ! =====================
[2392]	414	IF( nn_cla == 0 ) THEN
[1273]	415	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	416	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
	417	DO ji = mi0(ii0), mi1(ii1)
	418	DO jj = mj0(ij0), mj1(ij1)
	419	mbathy(ji,jj) = 15
	420	END DO
	421	END DO
	422	IF(lwp) WRITE(numout,*)
[2392]	423	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
[1273]	424	!
	425	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	426	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
	427	DO ji = mi0(ii0), mi1(ii1)
	428	DO jj = mj0(ij0), mj1(ij1)
	429	mbathy(ji,jj) = 12
	430	END DO
	431	END DO
	432	IF(lwp) WRITE(numout,*)
[2392]	433	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
[1273]	434	ENDIF
	435	!
	436	ENDIF
	437	!
[454]	438	ENDIF
[1099]	439	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
	440	CALL iom_open( 'bathy_meter.nc', inum )
[473]	441	CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
	442	CALL iom_close (inum)
[2392]	443	! ! =====================
[2380]	444	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
[2392]	445	ii0 = 142 ; ii1 = 142 ! =====================
	446	ij0 = 51 ; ij1 = 53
	447	DO ji = mi0(ii0), mi1(ii1) ! Close Halmera Strait
[2380]	448	DO jj = mj0(ij0), mj1(ij1)
	449	bathy(ji,jj) = 0.0
	450	END DO
	451	END DO
	452	IF(lwp) WRITE(numout,*)
[2392]	453	IF(lwp) WRITE(numout,*) ' orca_r1: Halmera strait closed at i=',ii0,' j=',ij0,'->',ij1
[2380]	454	ENDIF
[1348]	455	! ! =====================
[2380]	456	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	457	! ! =====================
[2392]	458	IF( nn_cla == 0 ) THEN
	459	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	460	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
[1273]	461	DO ji = mi0(ii0), mi1(ii1)
	462	DO jj = mj0(ij0), mj1(ij1)
	463	bathy(ji,jj) = 284.e0
	464	END DO
	465	END DO
	466	IF(lwp) WRITE(numout,*)
[2392]	467	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
[1273]	468	!
[2392]	469	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	470	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
[1273]	471	DO ji = mi0(ii0), mi1(ii1)
	472	DO jj = mj0(ij0), mj1(ij1)
	473	bathy(ji,jj) = 137.e0
	474	END DO
	475	END DO
	476	IF(lwp) WRITE(numout,*)
	477	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
	478	ENDIF
	479	!
	480	ENDIF
[1348]	481	!
	482	ENDIF
[3]	483	! ! =============== !
	484	ELSE ! error !
	485	! ! =============== !
[1099]	486	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
[473]	487	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
[3]	488	ENDIF
[1099]	489	!
	490	! ! =========================== !
	491	IF( nclosea == 0 ) THEN ! NO closed seas or lakes !
	492	DO jl = 1, jpncs ! =========================== !
[3]	493	DO jj = ncsj1(jl), ncsj2(jl)
	494	DO ji = ncsi1(jl), ncsi2(jl)
[1099]	495	mbathy(ji,jj) = 0 ! suppress closed seas and lakes from bathymetry
	496	bathy (ji,jj) = 0.e0
[3]	497	END DO
	498	END DO
	499	END DO
	500	ENDIF
[473]	501	#if defined key_orca_lev10
[1099]	502	! ! ================= !
	503	mbathy(:,:) = 10 * mbathy(:,:) ! key_orca_lev10 !
	504	! ! ================= !
	505	IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' )
[473]	506	#endif
[1099]	507	! ! =============== !
	508	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain !
	509	! ! =============== !
[3]	510
[2236]	511	#if ! defined key_c1d
	512	! ! =================== !
	513	CALL zgr_bat_ctl ! Bathymetry check !
	514	! ! =================== !
	515	#endif
[3]	516	END SUBROUTINE zgr_bat
	517
	518
	519	SUBROUTINE zgr_bat_zoom
	520	!!----------------------------------------------------------------------
	521	!! * ROUTINE zgr_bat_zoom *
	522	!!
	523	!! ** Purpose : - Close zoom domain boundary if necessary
	524	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
	525	!!
	526	!! ** Method :
	527	!!
	528	!! ** Action : - update mbathy: level bathymetry (in level index)
	529	!!----------------------------------------------------------------------
	530	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
	531	!!----------------------------------------------------------------------
[1099]	532	!
[3]	533	IF(lwp) WRITE(numout,*)
	534	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
	535	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
[1099]	536	!
[3]	537	! Zoom domain
	538	! ===========
[1099]	539	!
[3]	540	! Forced closed boundary if required
[1099]	541	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
	542	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
	543	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
	544	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
	545	!
[3]	546	! Configuration specific domain modifications
	547	! (here, ORCA arctic configuration: suppress Med Sea)
	548	IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
	549	SELECT CASE ( jp_cfg )
	550	! ! =======================
	551	CASE ( 2 ) ! ORCA_R2 configuration
	552	! ! =======================
	553	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
	554	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
	555	ij0 = 98 ; ij1 = 110
	556	! ! =======================
	557	CASE ( 05 ) ! ORCA_R05 configuration
	558	! ! =======================
	559	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
	560	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
	561	ij0 = 314 ; ij1 = 370
	562	END SELECT
	563	!
	564	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
	565	!
	566	ENDIF
[1099]	567	!
[3]	568	END SUBROUTINE zgr_bat_zoom
	569
	570
	571	SUBROUTINE zgr_bat_ctl
	572	!!----------------------------------------------------------------------
	573	!! * ROUTINE zgr_bat_ctl *
	574	!!
	575	!! ** Purpose : check the bathymetry in levels
	576	!!
	577	!! ** Method : The array mbathy is checked to verified its consistency
	578	!! with the model options. in particular:
	579	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	580	!! along closed boundary.
	581	!! mbathy must be cyclic IF jperio=1.
	582	!! mbathy must be lower or equal to jpk-1.
	583	!! isolated ocean grid points are suppressed from mbathy
	584	!! since they are only connected to remaining
	585	!! ocean through vertical diffusion.
	586	!! C A U T I O N : mbathy will be modified during the initializa-
	587	!! tion phase to become the number of non-zero w-levels of a water
	588	!! column, with a minimum value of 1.
	589	!!
	590	!! ** Action : - update mbathy: level bathymetry (in level index)
	591	!! - update bathy : meter bathymetry (in meters)
	592	!!----------------------------------------------------------------------
[1099]	593	INTEGER :: ji, jj, jl ! dummy loop indices
	594	INTEGER :: icompt, ibtest, ikmax ! temporary integers
	595	REAL(wp), DIMENSION(jpi,jpj) :: zbathy ! temporary workspace
[3]	596	!!----------------------------------------------------------------------
	597
	598	IF(lwp) WRITE(numout,*)
	599	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
	600	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
	601
[1099]	602	! ! Suppress isolated ocean grid points
	603	IF(lwp) WRITE(numout,*)
	604	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
	605	IF(lwp) WRITE(numout,*)' -----------------------------------'
	606	icompt = 0
	607	DO jl = 1, 2
	608	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
	609	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
	610	mbathy(jpi,:) = mbathy( 2 ,:)
	611	ENDIF
	612	DO jj = 2, jpjm1
	613	DO ji = 2, jpim1
	614	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
	615	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
	616	IF( ibtest < mbathy(ji,jj) ) THEN
	617	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
	618	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
	619	mbathy(ji,jj) = ibtest
	620	icompt = icompt + 1
	621	ENDIF
	622	END DO
	623	END DO
	624	END DO
	625	IF( icompt == 0 ) THEN
	626	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
	627	ELSE
	628	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
	629	ENDIF
	630	IF( lk_mpp ) THEN
	631	zbathy(:,:) = FLOAT( mbathy(:,:) )
	632	CALL lbc_lnk( zbathy, 'T', 1. )
	633	mbathy(:,:) = INT( zbathy(:,:) )
	634	ENDIF
[3]	635
[1099]	636	! ! East-west cyclic boundary conditions
	637	IF( nperio == 0 ) THEN
	638	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
	639	IF( lk_mpp ) THEN
	640	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
	641	IF( jperio /= 1 ) mbathy(1,:) = 0
[411]	642	ENDIF
[1099]	643	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
	644	IF( jperio /= 1 ) mbathy(nlci,:) = 0
	645	ENDIF
[411]	646	ELSE
[1099]	647	IF( ln_zco .OR. ln_zps ) THEN
	648	mbathy( 1 ,:) = 0
	649	mbathy(jpi,:) = 0
	650	ELSE
	651	mbathy( 1 ,:) = jpkm1
	652	mbathy(jpi,:) = jpkm1
	653	ENDIF
[411]	654	ENDIF
[1099]	655	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
	656	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
	657	mbathy( 1 ,:) = mbathy(jpim1,:)
	658	mbathy(jpi,:) = mbathy( 2 ,:)
	659	ELSEIF( nperio == 2 ) THEN
	660	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
	661	ELSE
	662	IF(lwp) WRITE(numout,*) ' e r r o r'
	663	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
	664	! STOP 'dom_mba'
	665	ENDIF
	666
[1528]	667	! Boundary condition on mbathy
	668	IF( .NOT.lk_mpp ) THEN
	669	!!gm !!bug ??? think about it !
	670	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
	671	zbathy(:,:) = FLOAT( mbathy(:,:) )
	672	CALL lbc_lnk( zbathy, 'T', 1. )
	673	mbathy(:,:) = INT( zbathy(:,:) )
[3]	674	ENDIF
	675
	676	! Number of ocean level inferior or equal to jpkm1
	677	ikmax = 0
	678	DO jj = 1, jpj
	679	DO ji = 1, jpi
	680	ikmax = MAX( ikmax, mbathy(ji,jj) )
	681	END DO
	682	END DO
[1099]	683	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
[3]	684	IF( ikmax > jpkm1 ) THEN
	685	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
	686	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
	687	ELSE IF( ikmax < jpkm1 ) THEN
	688	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
	689	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
	690	ENDIF
	691
[1566]	692	IF( lwp .AND. nprint == 1 ) THEN ! control print
[3]	693	WRITE(numout,*)
[1099]	694	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
[3]	695	WRITE(numout,*) ' ------------------'
[1099]	696	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
[3]	697	WRITE(numout,*)
	698	ENDIF
[1099]	699	!
[3]	700	END SUBROUTINE zgr_bat_ctl
	701
	702
[454]	703	SUBROUTINE zgr_zco
	704	!!----------------------------------------------------------------------
	705	!! * ROUTINE zgr_zco *
	706	!!
	707	!! ** Purpose : define the z-coordinate system
	708	!!
[2240]	709	!! ** Method : set 3D coord. arrays to reference 1D array
[454]	710	!!----------------------------------------------------------------------
	711	INTEGER :: jk
	712	!!----------------------------------------------------------------------
[1099]	713	!
[2240]	714	DO jk = 1, jpk
	715	fsdept(:,:,jk) = gdept_0(jk)
	716	fsdepw(:,:,jk) = gdepw_0(jk)
	717	fsde3w(:,:,jk) = gdepw_0(jk)
	718	fse3t (:,:,jk) = e3t_0(jk)
	719	fse3u (:,:,jk) = e3t_0(jk)
	720	fse3v (:,:,jk) = e3t_0(jk)
	721	fse3f (:,:,jk) = e3t_0(jk)
	722	fse3w (:,:,jk) = e3w_0(jk)
	723	fse3uw(:,:,jk) = e3w_0(jk)
	724	fse3vw(:,:,jk) = e3w_0(jk)
	725	END DO
[1099]	726	!
[454]	727	END SUBROUTINE zgr_zco
	728
[1083]	729	!!----------------------------------------------------------------------
	730	!! Default option : zco, zps and/or sco available (gedp & e3 are 3D arrays)
	731	!!----------------------------------------------------------------------
[454]	732
[1083]	733	SUBROUTINE zgr_zps
	734	!!----------------------------------------------------------------------
	735	!! * ROUTINE zgr_zps *
	736	!!
	737	!! ** Purpose : the depth and vertical scale factor in partial step
	738	!! z-coordinate case
	739	!!
	740	!! ** Method : Partial steps : computes the 3D vertical scale factors
	741	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
	742	!! a partial step representation of bottom topography.
	743	!!
	744	!! The reference depth of model levels is defined from an analytical
	745	!! function the derivative of which gives the reference vertical
	746	!! scale factors.
	747	!! From depth and scale factors reference, we compute there new value
	748	!! with partial steps on 3d arrays ( i, j, k ).
	749	!!
	750	!! w-level: gdepw(i,j,k) = fsdep(k)
	751	!! e3w(i,j,k) = dk(fsdep)(k) = fse3(i,j,k)
	752	!! t-level: gdept(i,j,k) = fsdep(k+0.5)
	753	!! e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
	754	!!
	755	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
	756	!! we find the mbathy index of the depth at each grid point.
	757	!! This leads us to three cases:
	758	!!
	759	!! - bathy = 0 => mbathy = 0
	760	!! - 1 < mbathy < jpkm1
	761	!! - bathy > gdepw(jpk) => mbathy = jpkm1
	762	!!
	763	!! Then, for each case, we find the new depth at t- and w- levels
	764	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
	765	!! and f-points.
	766	!!
	767	!! This routine is given as an example, it must be modified
	768	!! following the user s desiderata. nevertheless, the output as
	769	!! well as the way to compute the model levels and scale factors
	770	!! must be respected in order to insure second order accuracy
	771	!! schemes.
	772	!!
	773	!! c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
	774	!! - - - - - - - gdept, gdepw and e3. are positives
	775	!!
[1099]	776	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
[1083]	777	!!----------------------------------------------------------------------
[1099]	778	INTEGER :: ji, jj, jk ! dummy loop indices
	779	INTEGER :: ik, it ! temporary integers
[2379]	780	INTEGER :: nlbelow ! temporary integer
[1099]	781	LOGICAL :: ll_print ! Allow control print for debugging
	782	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
	783	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
	784	REAL(wp) :: zmax, zmin ! Maximum and minimum depth
	785	REAL(wp) :: zdiff ! temporary scalar
[2379]	786	REAL(wp) :: zrefdep ! temporary scalar
[1099]	787	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprt ! 3D workspace
	788	!!---------------------------------------------------------------------
[3]	789
[1099]	790	IF(lwp) WRITE(numout,*)
	791	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
	792	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
	793	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
[3]	794
[1099]	795	ll_print=.FALSE. ! Local variable for debugging
	796	!! ll_print=.TRUE.
[1083]	797
[1099]	798	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
[1083]	799	WRITE(numout,*)
	800	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
	801	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
	802	ENDIF
	803
	804
	805	! bathymetry in level (from bathy_meter)
	806	! ===================
[1099]	807	zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
[1083]	808
[2379]	809	zrefdep = 25.
	810	nlbelow = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Below zrefdep
	811	IF(lwp) write(numout,*) 'Minimum depth level selected: ', nlbelow
	812	zmin = gdepw_0(nlbelow) ! minimum depth = nlbelow-1 levels
	813
[1099]	814	mbathy(:,:) = jpkm1 ! initialize mbathy to the maximum ocean level available
[1083]	815
[1099]	816	! ! storage of land and island's number (zera and negative values) in mbathy
	817	WHERE( bathy(:,:) <= 0. ) mbathy(:,:) = INT( bathy(:,:) )
[1083]	818
	819	! bounded value of bathy
[1099]	820	!!gm bathy(:,:) = MIN( zmax, MAX( bathy(:,:), zmin ) ) !!gm : simpler as zmin is > 0
[1083]	821	DO jj = 1, jpj
	822	DO ji= 1, jpi
[1099]	823	IF( bathy(ji,jj) <= 0. ) THEN ; bathy(ji,jj) = 0.e0
	824	ELSE ; bathy(ji,jj) = MIN( zmax, MAX( bathy(ji,jj), zmin ) )
[1083]	825	ENDIF
	826	END DO
	827	END DO
	828
	829	! Compute mbathy for ocean points (i.e. the number of ocean levels)
	830	! find the number of ocean levels such that the last level thickness
	831	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
	832	! e3t_0 is the reference level thickness
	833	DO jk = jpkm1, 1, -1
	834	zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
	835	DO jj = 1, jpj
	836	DO ji = 1, jpi
	837	IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth ) mbathy(ji,jj) = jk-1
	838	END DO
	839	END DO
	840	END DO
	841
[1099]	842	! Scale factors and depth at T- and W-points
	843	DO jk = 1, jpk ! intitialization to the reference z-coordinate
	844	gdept(:,:,jk) = gdept_0(jk)
	845	gdepw(:,:,jk) = gdepw_0(jk)
	846	e3t (:,:,jk) = e3t_0 (jk)
	847	e3w (:,:,jk) = e3w_0 (jk)
[1083]	848	END DO
[1099]	849	hdept(:,:) = gdept(:,:,2 )
	850	hdepw(:,:) = gdepw(:,:,3 )
	851	!
	852	DO jj = 1, jpj
	853	DO ji = 1, jpi
	854	ik = mbathy(ji,jj)
	855	IF( ik > 0 ) THEN ! ocean point only
	856	! max ocean level case
	857	IF( ik == jpkm1 ) THEN
	858	zdepwp = bathy(ji,jj)
	859	ze3tp = bathy(ji,jj) - gdepw_0(ik)
	860	ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
	861	e3t(ji,jj,ik ) = ze3tp
	862	e3t(ji,jj,ik+1) = ze3tp
	863	e3w(ji,jj,ik ) = ze3wp
	864	e3w(ji,jj,ik+1) = ze3tp
	865	gdepw(ji,jj,ik+1) = zdepwp
	866	gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp
	867	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
	868	!
	869	ELSE ! standard case
	870	IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj)
	871	ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
	872	ENDIF
	873	!gm Bug? check the gdepw_0
	874	! ... on ik
	875	gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
	876	& * ((gdept_0( ik ) - gdepw_0(ik) ) &
	877	& / ( gdepw_0( ik+1) - gdepw_0(ik) ))
	878	e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
	879	& / ( gdepw_0( ik+1) - gdepw_0(ik) )
	880	e3w (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) &
	881	& * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
	882	! ... on ik+1
	883	e3w (ji,jj,ik+1) = e3t (ji,jj,ik)
	884	e3t (ji,jj,ik+1) = e3t (ji,jj,ik)
	885	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
	886	ENDIF
	887	ENDIF
	888	END DO
	889	END DO
	890	!
	891	it = 0
	892	DO jj = 1, jpj
	893	DO ji = 1, jpi
	894	ik = mbathy(ji,jj)
	895	IF( ik > 0 ) THEN ! ocean point only
	896	hdept(ji,jj) = gdept(ji,jj,ik )
	897	hdepw(ji,jj) = gdepw(ji,jj,ik+1)
	898	e3tp (ji,jj) = e3t(ji,jj,ik )
	899	e3wp (ji,jj) = e3w(ji,jj,ik )
	900	! test
	901	zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik )
	902	IF( zdiff <= 0. .AND. lwp ) THEN
	903	it = it + 1
	904	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
	905	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
	906	WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
	907	WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik )
	908	ENDIF
	909	ENDIF
	910	END DO
	911	END DO
[1083]	912
	913	! Scale factors and depth at U-, V-, UW and VW-points
[1099]	914	DO jk = 1, jpk ! initialisation to z-scale factors
[1083]	915	e3u (:,:,jk) = e3t_0(jk)
	916	e3v (:,:,jk) = e3t_0(jk)
	917	e3uw(:,:,jk) = e3w_0(jk)
	918	e3vw(:,:,jk) = e3w_0(jk)
	919	END DO
[1099]	920	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
	921	DO jj = 1, jpjm1
	922	DO ji = 1, fs_jpim1 ! vector opt.
	923	e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
	924	e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
	925	e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
	926	e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
	927	END DO
	928	END DO
	929	END DO
	930	! ! Boundary conditions
	931	CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. )
	932	CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. )
	933	!
	934	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
	935	WHERE( e3u (:,:,jk) == 0.e0 ) e3u (:,:,jk) = e3t_0(jk)
	936	WHERE( e3v (:,:,jk) == 0.e0 ) e3v (:,:,jk) = e3t_0(jk)
	937	WHERE( e3uw(:,:,jk) == 0.e0 ) e3uw(:,:,jk) = e3w_0(jk)
	938	WHERE( e3vw(:,:,jk) == 0.e0 ) e3vw(:,:,jk) = e3w_0(jk)
	939	END DO
	940
	941	! Scale factor at F-point
	942	DO jk = 1, jpk ! initialisation to z-scale factors
	943	e3f (:,:,jk) = e3t_0(jk)
	944	END DO
	945	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
	946	DO jj = 1, jpjm1
	947	DO ji = 1, fs_jpim1 ! vector opt.
	948	e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
	949	END DO
	950	END DO
	951	END DO
	952	CALL lbc_lnk( e3f, 'F', 1. ) ! Boundary conditions
	953	!
	954	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
	955	WHERE( e3f(:,:,jk) == 0.e0 ) e3f(:,:,jk) = e3t_0(jk)
	956	END DO
	957	!!gm bug ? : must be a do loop with mj0,mj1
	958	!
	959	e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
	960	e3w(:,mj0(1),:) = e3w(:,mj0(2),:)
	961	e3u(:,mj0(1),:) = e3u(:,mj0(2),:)
	962	e3v(:,mj0(1),:) = e3v(:,mj0(2),:)
	963	e3f(:,mj0(1),:) = e3f(:,mj0(2),:)
[1083]	964
[1099]	965	! Control of the sign
	966	IF( MINVAL( e3t (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' )
	967	IF( MINVAL( e3w (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' )
	968	IF( MINVAL( gdept(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
	969	IF( MINVAL( gdepw(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
[1083]	970
[1099]	971	! Compute gdep3w (vertical sum of e3w)
	972	gdep3w(:,:,1) = 0.5 * e3w(:,:,1)
	973	DO jk = 2, jpk
	974	gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk)
	975	END DO
[1083]	976
[1099]	977	! ! ================= !
	978	IF(lwp .AND. ll_print) THEN ! Control print !
	979	! ! ================= !
[1083]	980	DO jj = 1,jpj
	981	DO ji = 1, jpi
[1099]	982	ik = MAX( mbathy(ji,jj), 1 )
	983	zprt(ji,jj,1) = e3t (ji,jj,ik)
	984	zprt(ji,jj,2) = e3w (ji,jj,ik)
	985	zprt(ji,jj,3) = e3u (ji,jj,ik)
	986	zprt(ji,jj,4) = e3v (ji,jj,ik)
	987	zprt(ji,jj,5) = e3f (ji,jj,ik)
	988	zprt(ji,jj,6) = gdep3w(ji,jj,ik)
[1083]	989	END DO
	990	END DO
	991	WRITE(numout,*)
[1099]	992	WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	993	WRITE(numout,*)
[1099]	994	WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	995	WRITE(numout,*)
[1099]	996	WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	997	WRITE(numout,*)
[1099]	998	WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	999	WRITE(numout,*)
[1099]	1000	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1001	WRITE(numout,*)
[1099]	1002	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1003	ENDIF
[1099]	1004
	1005	! ! =============== !
	1006	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain !
	1007	! ! =============== !
[1083]	1008
[2236]	1009	#if ! defined key_c1d
	1010	! ! =================== !
	1011	CALL zgr_bat_ctl ! Bathymetry check !
	1012	! ! =================== !
	1013	#endif
[1083]	1014	END SUBROUTINE zgr_zps
	1015
	1016
	1017	FUNCTION fssig( pk ) RESULT( pf )
	1018	!!----------------------------------------------------------------------
	1019	!! * ROUTINE eos_init *
	1020	!!
	1021	!! ** Purpose : provide the analytical function in s-coordinate
	1022	!!
	1023	!! ** Method : the function provide the non-dimensional position of
	1024	!! T and W (i.e. between 0 and 1)
	1025	!! T-points at integer values (between 1 and jpk)
	1026	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1027	!!
	1028	!! Reference : ???
	1029	!!----------------------------------------------------------------------
	1030	REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate
	1031	REAL(wp) :: pf ! sigma value
	1032	!!----------------------------------------------------------------------
	1033	!
[1601]	1034	pf = ( TANH( rn_theta * ( -(pk-0.5) / REAL(jpkm1) + rn_thetb ) ) &
	1035	& - TANH( rn_thetb * rn_theta ) ) &
	1036	& * ( COSH( rn_theta ) &
	1037	& + COSH( rn_theta * ( 2.e0 * rn_thetb - 1.e0 ) ) ) &
	1038	& / ( 2.e0 * SINH( rn_theta ) )
[1083]	1039	!
	1040	END FUNCTION fssig
	1041
	1042
[1601]	1043	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
[1348]	1044	!!----------------------------------------------------------------------
	1045	!! * ROUTINE eos_init *
	1046	!!
	1047	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
	1048	!!
	1049	!! ** Method : the function provides the non-dimensional position of
	1050	!! T and W (i.e. between 0 and 1)
	1051	!! T-points at integer values (between 1 and jpk)
	1052	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1053	!!
	1054	!! Reference : ???
	1055	!!----------------------------------------------------------------------
	1056	REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate
[1601]	1057	REAL(wp), INTENT(in ) :: pbb ! Stretching coefficient
[1348]	1058	REAL(wp) :: pf1 ! sigma value
	1059	!!----------------------------------------------------------------------
	1060	!
[1601]	1061	IF ( rn_theta == 0 ) then ! uniform sigma
[1566]	1062	pf1 = -(pk1-0.5) / REAL( jpkm1 )
	1063	ELSE ! stretched sigma
[1601]	1064	pf1 = (1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(rn_theta) + &
	1065	& pbb * ( (tanh( rn_theta( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5rn_theta) ) / &
	1066	& (2tanh(0.5rn_theta) ) )
[1348]	1067	ENDIF
[1566]	1068	!
[1348]	1069	END FUNCTION fssig1
	1070
	1071
[454]	1072	SUBROUTINE zgr_sco
	1073	!!----------------------------------------------------------------------
	1074	!! * ROUTINE zgr_sco *
	1075	!!
	1076	!! ** Purpose : define the s-coordinate system
	1077	!!
	1078	!! ** Method : s-coordinate
	1079	!! The depth of model levels is defined as the product of an
	1080	!! analytical function by the local bathymetry, while the vertical
	1081	!! scale factors are defined as the product of the first derivative
	1082	!! of the analytical function by the bathymetry.
	1083	!! (this solution save memory as depth and scale factors are not
	1084	!! 3d fields)
	1085	!! - Read bathymetry (in meters) at t-point and compute the
	1086	!! bathymetry at u-, v-, and f-points.
	1087	!! hbatu = mi( hbatt )
	1088	!! hbatv = mj( hbatt )
	1089	!! hbatf = mi( mj( hbatt ) )
	1090	!! - Compute gsigt, gsigw, esigt, esigw from an analytical
[1083]	1091	!! function and its derivative given as function.
	1092	!! gsigt(k) = fssig (k )
	1093	!! gsigw(k) = fssig (k-0.5)
	1094	!! esigt(k) = fsdsig(k )
	1095	!! esigw(k) = fsdsig(k-0.5)
[454]	1096	!! This routine is given as an example, it must be modified
	1097	!! following the user s desiderata. nevertheless, the output as
	1098	!! well as the way to compute the model levels and scale factors
	1099	!! must be respected in order to insure second order a!!uracy
	1100	!! schemes.
	1101	!!
[1099]	1102	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
	1103	!!----------------------------------------------------------------------
	1104	INTEGER :: ji, jj, jk, jl ! dummy loop argument
	1105	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
	1106	REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper ! temporary scalars
	1107	REAL(wp), DIMENSION(jpi,jpj) :: zenv, ztmp, zmsk ! 2D workspace
	1108	REAL(wp), DIMENSION(jpi,jpj) :: zri , zrj , zhbat ! - -
[1566]	1109	!!
[1629]	1110	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigw3
	1111	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigt3
	1112	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsi3w3
	1113	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigt3
	1114	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigw3
	1115	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtu3
	1116	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtv3
	1117	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtf3
	1118	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigwu3
	1119	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigwv3
[454]	1120	!!
[1601]	1121	NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
[454]	1122	!!----------------------------------------------------------------------
	1123
[1601]	1124	REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters
	1125	READ ( numnam, namzgr_sco )
[454]	1126
[1099]	1127	IF(lwp) THEN ! control print
[454]	1128	WRITE(numout,*)
	1129	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
	1130	WRITE(numout,*) '~~~~~~~~~~~'
[1601]	1131	WRITE(numout,*) ' Namelist namzgr_sco'
	1132	WRITE(numout,*) ' sigma-stretching coeffs '
	1133	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ' ,rn_sbot_max
	1134	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ' ,rn_sbot_min
	1135	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta
	1136	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb
	1137	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax
	1138	WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma
	1139	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb
	1140	WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc
[454]	1141	ENDIF
	1142
[1629]	1143	gsigw3 = 0.0d0 ; gsigt3 = 0.0d0 ; gsi3w3 = 0.0d0
	1144	esigt3 = 0.0d0 ; esigw3 = 0.0d0
	1145	esigtu3 = 0.0d0 ; esigtv3 = 0.0d0 ; esigtf3 = 0.0d0
	1146	esigwu3 = 0.0d0 ; esigwv3 = 0.0d0
	1147
[1601]	1148	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
	1149	hifu(:,:) = rn_sbot_min
	1150	hifv(:,:) = rn_sbot_min
	1151	hiff(:,:) = rn_sbot_min
[1348]	1152
	1153	! ! set maximum ocean depth
[1601]	1154	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
[454]	1155
[1461]	1156	DO jj = 1, jpj
	1157	DO ji = 1, jpi
	1158	IF (bathy(ji,jj) .gt. 0.0) THEN
[1601]	1159	bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
[1461]	1160	ENDIF
	1161	END DO
	1162	END DO
[1099]	1163	! ! =============================
	1164	! ! Define the envelop bathymetry (hbatt)
	1165	! ! =============================
[454]	1166	! use r-value to create hybrid coordinates
	1167	DO jj = 1, jpj
	1168	DO ji = 1, jpi
[1601]	1169	zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
[454]	1170	END DO
	1171	END DO
[1639]	1172	!
	1173	! Smooth the bathymetry (if required)
	1174	scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea)
	1175	scobot(:,:) = bathy(:,:) ! ocean bottom depth
	1176	!
[454]	1177	jl = 0
	1178	zrmax = 1.e0
	1179	! ! ================ !
[1601]	1180	DO WHILE ( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop !
[454]	1181	! ! ================ !
	1182	jl = jl + 1
	1183	zrmax = 0.e0
	1184	zmsk(:,:) = 0.e0
	1185	DO jj = 1, nlcj
	1186	DO ji = 1, nlci
	1187	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
	1188	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
	1189	zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
	1190	zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
	1191	zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
[1601]	1192	IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1.0
	1193	IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1.0
	1194	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1.0
	1195	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1.0
[454]	1196	END DO
	1197	END DO
	1198	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
	1199	! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
[1099]	1200	ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1. )
[454]	1201	DO jj = 1, nlcj
	1202	DO ji = 1, nlci
[1348]	1203	zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
[454]	1204	END DO
	1205	END DO
[1348]	1206	!
[1099]	1207	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
	1208	!
[454]	1209	DO jj = 1, nlcj
	1210	DO ji = 1, nlci
	1211	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
	1212	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
	1213	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
	1214	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
	1215	IF( zmsk(ji,jj) == 1.0 ) THEN
	1216	ztmp(ji,jj) = ( &
	1217	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
	1218	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2.e0 + zenv(iip1,jj )*zmsk(iip1,jj ) &
	1219	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
	1220	& ) / ( &
	1221	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
	1222	& + zmsk(iim1,jj ) + 2.e0 + zmsk(iip1,jj ) &
	1223	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
	1224	& )
	1225	ENDIF
	1226	END DO
	1227	END DO
[1348]	1228	!
[454]	1229	DO jj = 1, nlcj
	1230	DO ji = 1, nlci
	1231	IF( zmsk(ji,jj) == 1.0 ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
	1232	END DO
	1233	END DO
[1099]	1234	!
[454]	1235	! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero
[1099]	1236	ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1. )
[454]	1237	DO jj = 1, nlcj
	1238	DO ji = 1, nlci
	1239	IF( zenv(ji,jj) == 0.e0 ) zenv(ji,jj) = ztmp(ji,jj)
	1240	END DO
	1241	END DO
	1242	! ! ================ !
	1243	END DO ! End loop !
	1244	! ! ================ !
[1099]	1245	!
	1246	! ! envelop bathymetry saved in hbatt
[454]	1247	hbatt(:,:) = zenv(:,:)
[1099]	1248	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) THEN
	1249	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
	1250	DO jj = 1, jpj
	1251	DO ji = 1, jpi
	1252	ztaper = EXP( -(gphit(ji,jj)/8)**2 )
[1601]	1253	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
[1099]	1254	END DO
	1255	END DO
[516]	1256	ENDIF
[1099]	1257	!
	1258	IF(lwp) THEN ! Control print
[454]	1259	WRITE(numout,*)
	1260	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
	1261	WRITE(numout,*)
	1262	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
[1099]	1263	IF( nprint == 1 ) THEN
	1264	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
	1265	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
	1266	ENDIF
[454]	1267	ENDIF
	1268
[1099]	1269	! ! ==============================
	1270	! ! hbatu, hbatv, hbatf fields
	1271	! ! ==============================
[454]	1272	IF(lwp) THEN
	1273	WRITE(numout,*)
[1601]	1274	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
[454]	1275	ENDIF
[1601]	1276	hbatu(:,:) = rn_sbot_min
	1277	hbatv(:,:) = rn_sbot_min
	1278	hbatf(:,:) = rn_sbot_min
[454]	1279	DO jj = 1, jpjm1
[1694]	1280	DO ji = 1, jpim1 ! NO vector opt.
[1099]	1281	hbatu(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
	1282	hbatv(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
	1283	hbatf(ji,jj) = 0.25* ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
	1284	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
[454]	1285	END DO
	1286	END DO
[1099]	1287	!
[454]	1288	! Apply lateral boundary condition
[1099]	1289	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
	1290	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1. )
[454]	1291	DO jj = 1, jpj
	1292	DO ji = 1, jpi
	1293	IF( hbatu(ji,jj) == 0.e0 ) THEN
[1601]	1294	IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = rn_sbot_min
[454]	1295	IF( zhbat(ji,jj) /= 0.e0 ) hbatu(ji,jj) = zhbat(ji,jj)
	1296	ENDIF
	1297	END DO
	1298	END DO
[1099]	1299	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1. )
[454]	1300	DO jj = 1, jpj
	1301	DO ji = 1, jpi
	1302	IF( hbatv(ji,jj) == 0.e0 ) THEN
[1601]	1303	IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = rn_sbot_min
[454]	1304	IF( zhbat(ji,jj) /= 0.e0 ) hbatv(ji,jj) = zhbat(ji,jj)
	1305	ENDIF
	1306	END DO
	1307	END DO
[1099]	1308	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1. )
[454]	1309	DO jj = 1, jpj
	1310	DO ji = 1, jpi
	1311	IF( hbatf(ji,jj) == 0.e0 ) THEN
[1601]	1312	IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = rn_sbot_min
[454]	1313	IF( zhbat(ji,jj) /= 0.e0 ) hbatf(ji,jj) = zhbat(ji,jj)
	1314	ENDIF
	1315	END DO
	1316	END DO
	1317
	1318	!!bug: key_helsinki a verifer
	1319	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
	1320	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
	1321	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
	1322	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
	1323
[516]	1324	IF( nprint == 1 .AND. lwp ) THEN
[1099]	1325	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
	1326	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
	1327	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
	1328	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
[516]	1329	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
	1330	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
	1331	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
	1332	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
	1333	ENDIF
[454]	1334	!! helsinki
	1335
[1099]	1336	! ! =======================
	1337	! ! s-ccordinate fields (gdep., e3.)
	1338	! ! =======================
	1339	!
	1340	! non-dimensional "sigma" for model level depth at w- and t-levels
[1348]	1341
[1601]	1342	IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths
	1343	! ! below rn_hc, with uniform sigma in shallower waters
	1344	DO ji = 1, jpi
	1345	DO jj = 1, jpj
[1348]	1346
[1601]	1347	IF (hbatt(ji,jj).GT.rn_hc) THEN !deep water, stretched sigma
[1348]	1348	DO jk = 1, jpk
[1601]	1349	gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
	1350	gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
[1348]	1351	END DO
	1352	ELSE ! shallow water, uniform sigma
	1353	DO jk = 1, jpk
	1354	gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
	1355	gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
	1356	END DO
	1357	ENDIF
	1358	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
	1359
	1360
	1361	DO jk = 1, jpkm1
	1362	esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
	1363	esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
	1364	END DO
[1639]	1365	esigw3(ji,jj,1 ) = 2.0_wp * (gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ))
	1366	esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk))
[1348]	1367
	1368	! Coefficients for vertical depth as the sum of e3w scale factors
	1369	gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1)
	1370	DO jk = 2, jpk
	1371	gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
	1372	END DO
	1373
	1374	DO jk = 1, jpk
	1375	zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp)
	1376	zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp)
[1601]	1377	gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsigt3(ji,jj,jk)+rn_hczcoeft)
	1378	gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsigw3(ji,jj,jk)+rn_hczcoefw)
[1639]	1379	gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsi3w3(ji,jj,jk)+rn_hczcoeft)
[1348]	1380	END DO
	1381
	1382	ENDDO ! for all jj's
	1383	ENDDO ! for all ji's
	1384
	1385
	1386	DO ji=1,jpi
	1387	DO jj=1,jpj
	1388
	1389	DO jk = 1, jpk
	1390	esigtu3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji+1,jj)esigt3(ji+1,jj,jk) ) / &
	1391	( hbatt(ji,jj)+hbatt(ji+1,jj) )
	1392	esigtv3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji,jj+1)esigt3(ji,jj+1,jk) ) / &
	1393	( hbatt(ji,jj)+hbatt(ji,jj+1) )
	1394	esigtf3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji+1,jj)esigt3(ji+1,jj,jk) + &
	1395	hbatt(ji,jj+1)esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)esigt3(ji+1,jj+1,jk) ) / &
	1396	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
	1397	esigwu3(ji,jj,jk) = ( hbatt(ji,jj)esigw3(ji,jj,jk)+hbatt(ji+1,jj)esigw3(ji+1,jj,jk) ) / &
	1398	( hbatt(ji,jj)+hbatt(ji+1,jj) )
	1399	esigwv3(ji,jj,jk) = ( hbatt(ji,jj)esigw3(ji,jj,jk)+hbatt(ji,jj+1)esigw3(ji,jj+1,jk) ) / &
	1400	( hbatt(ji,jj)+hbatt(ji,jj+1) )
	1401
[1601]	1402	e3t(ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigt3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
	1403	e3u(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
	1404	e3v(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
	1405	e3f(ji,jj,jk)=((hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
[1348]	1406	!
[1601]	1407	e3w (ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigw3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
	1408	e3uw(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
	1409	e3vw(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
[1348]	1410	END DO
	1411
	1412	ENDDO
	1413	ENDDO
	1414
	1415	ELSE ! not ln_s_sigma
	1416
	1417	DO jk = 1, jpk
	1418	gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
	1419	gsigt(jk) = -fssig( REAL(jk,wp) )
	1420	END DO
	1421	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)
	1422	!
	1423	! Coefficients for vertical scale factors at w-, t- levels
[1099]	1424	!!gm bug : define it from analytical function, not like juste bellow....
	1425	!!gm or betteroffer the 2 possibilities....
[1348]	1426	DO jk = 1, jpkm1
	1427	esigt(jk ) = gsigw(jk+1) - gsigw(jk)
	1428	esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
	1429	END DO
[1639]	1430	esigw( 1 ) = 2.0_wp * (gsigt(1) - gsigw(1))
	1431	esigt(jpk) = 2.0_wp * (gsigt(jpk) - gsigw(jpk))
[454]	1432
[1099]	1433	!!gm original form
	1434	!!org DO jk = 1, jpk
	1435	!!org esigt(jk)=fsdsig( FLOAT(jk) )
	1436	!!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
	1437	!!org END DO
	1438	!!gm
[1348]	1439	!
	1440	! Coefficients for vertical depth as the sum of e3w scale factors
	1441	gsi3w(1) = 0.5 * esigw(1)
	1442	DO jk = 2, jpk
	1443	gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
	1444	END DO
[1099]	1445	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
[1348]	1446	DO jk = 1, jpk
	1447	zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
	1448	zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
	1449	gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigt(jk)+hift(:,:)zcoeft)
	1450	gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigw(jk)+hift(:,:)zcoefw)
[1639]	1451	gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsi3w(jk)+hift(:,:)zcoeft)
[1348]	1452	END DO
[1099]	1453	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
[1348]	1454	DO jj = 1, jpj
	1455	DO ji = 1, jpi
	1456	DO jk = 1, jpk
	1457	e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
	1458	e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
	1459	e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
	1460	e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
	1461	!
	1462	e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
	1463	e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
	1464	e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
	1465	END DO
[454]	1466	END DO
	1467	END DO
[1348]	1468
	1469	ENDIF ! ln_s_sigma
	1470
	1471
[1099]	1472	!
[1461]	1473	!! H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code.
	1474
	1475	fsdept(:,:,:) = gdept (:,:,:)
	1476	fsdepw(:,:,:) = gdepw (:,:,:)
	1477	fsde3w(:,:,:) = gdep3w(:,:,:)
	1478	fse3t (:,:,:) = e3t (:,:,:)
	1479	fse3u (:,:,:) = e3u (:,:,:)
	1480	fse3v (:,:,:) = e3v (:,:,:)
	1481	fse3f (:,:,:) = e3f (:,:,:)
	1482	fse3w (:,:,:) = e3w (:,:,:)
	1483	fse3uw(:,:,:) = e3uw (:,:,:)
	1484	fse3vw(:,:,:) = e3vw (:,:,:)
	1485	!!
[1099]	1486	! HYBRID :
[454]	1487	DO jj = 1, jpj
	1488	DO ji = 1, jpi
	1489	DO jk = 1, jpkm1
	1490	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
	1491	IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0
	1492	END DO
	1493	END DO
	1494	END DO
[1099]	1495	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
	1496	& ' MAX ', MAXVAL( mbathy(:,:) )
[454]	1497
	1498
[1632]	1499	! ! ===========
	1500	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain
	1501	! ! ===========
[454]	1502
[2236]	1503	#if ! defined key_c1d
	1504	! ! =================== !
	1505	CALL zgr_bat_ctl ! Bathymetry check !
	1506	! ! =================== !
	1507	#endif
[1632]	1508
	1509	! ! =============
	1510	IF(lwp) THEN ! Control print
	1511	! ! =============
[454]	1512	WRITE(numout,*)
	1513	WRITE(numout,*) ' domzgr: vertical coefficients for model level'
[1099]	1514	WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" )
	1515	WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
[454]	1516	ENDIF
[1099]	1517	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
	1518	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	1519	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
	1520	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) )
	1521	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), &
	1522	& ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), &
	1523	& ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), &
	1524	& ' w ', MINVAL( fse3w (:,:,:) )
[454]	1525
[1099]	1526	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
	1527	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) )
	1528	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), &
	1529	& ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), &
	1530	& ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), &
	1531	& ' w ', MAXVAL( fse3w (:,:,:) )
	1532	ENDIF
	1533	!
	1534	IF(lwp) THEN ! selected vertical profiles
[454]	1535	WRITE(numout,*)
	1536	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
	1537	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1538	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1539	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), &
	1540	& fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
[473]	1541	DO jj = mj0(20), mj1(20)
	1542	DO ji = mi0(20), mi1(20)
	1543	WRITE(numout,*)
	1544	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	1545	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1546	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1547	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
	1548	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
[473]	1549	END DO
	1550	END DO
	1551	DO jj = mj0(74), mj1(74)
	1552	DO ji = mi0(100), mi1(100)
	1553	WRITE(numout,*)
	1554	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	1555	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1556	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1557	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
	1558	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
[473]	1559	END DO
	1560	END DO
[454]	1561	ENDIF
	1562
[1099]	1563	!!gm bug? no more necessary? if ! defined key_helsinki
[454]	1564	DO jk = 1, jpk
	1565	DO jj = 1, jpj
	1566	DO ji = 1, jpi
	1567	IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
[1099]	1568	WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
	1569	CALL ctl_stop( ctmp1 )
[454]	1570	ENDIF
	1571	IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
[1099]	1572	WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
	1573	CALL ctl_stop( ctmp1 )
[454]	1574	ENDIF
	1575	END DO
	1576	END DO
	1577	END DO
[1099]	1578	!!gm bug #endif
	1579	!
[454]	1580	END SUBROUTINE zgr_sco
	1581
	1582
[3]	1583	!!======================================================================
	1584	END MODULE domzgr

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

Download in other formats: