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 @ 2380

Last change on this file since 2380 was 2380, checked in by acc, 13 years ago
nemo_v3_3beta. ORCA_R1 settings (step 2, see ticket #758). Introduces key_orca_r1 (46 level default, 75 level if key_orca_r1=75)
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	! ! =====================
	414	IF( n_cla == 0 ) THEN
	415	!
	416	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	417	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
	418	DO ji = mi0(ii0), mi1(ii1)
	419	DO jj = mj0(ij0), mj1(ij1)
	420	mbathy(ji,jj) = 15
	421	END DO
	422	END DO
	423	IF(lwp) WRITE(numout,*)
	424	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
	425	!
	426	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	427	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
	428	DO ji = mi0(ii0), mi1(ii1)
	429	DO jj = mj0(ij0), mj1(ij1)
	430	mbathy(ji,jj) = 12
	431	END DO
	432	END DO
	433	IF(lwp) WRITE(numout,*)
	434	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
	435	!
	436	ENDIF
	437	!
	438	ENDIF
	439	!
[454]	440	ENDIF
[1099]	441	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
	442	CALL iom_open( 'bathy_meter.nc', inum )
[473]	443	CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
	444	CALL iom_close (inum)
[2380]	445	! ! =====================
	446	IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration
	447	ii0 = 142 ; ii1 = 142 ! Close Halmera Strait
	448	ij0 = 51 ; ij1 = 53 ! =====================
	449	DO ji = mi0(ii0), mi1(ii1)
	450	DO jj = mj0(ij0), mj1(ij1)
	451	bathy(ji,jj) = 0.0
	452	END DO
	453	END DO
	454	IF(lwp) WRITE(numout,*)
	455	IF(lwp) WRITE(numout,*) ' orca_r1: Halmera strait closed at i=',ii0,' j=',ij0,'->',ij1
	456	ENDIF
[1348]	457	! ! =====================
[2380]	458	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
	459	! ! =====================
[1273]	460	IF( n_cla == 0 ) THEN
	461	!
	462	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
	463	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
	464	DO ji = mi0(ii0), mi1(ii1)
	465	DO jj = mj0(ij0), mj1(ij1)
	466	bathy(ji,jj) = 284.e0
	467	END DO
	468	END DO
	469	IF(lwp) WRITE(numout,*)
	470	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
	471	!
	472	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
	473	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
	474	DO ji = mi0(ii0), mi1(ii1)
	475	DO jj = mj0(ij0), mj1(ij1)
	476	bathy(ji,jj) = 137.e0
	477	END DO
	478	END DO
	479	IF(lwp) WRITE(numout,*)
	480	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
	481	!
	482	ENDIF
	483	!
	484	ENDIF
[1348]	485	!
	486	ENDIF
[3]	487	! ! =============== !
	488	ELSE ! error !
	489	! ! =============== !
[1099]	490	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
[473]	491	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
[3]	492	ENDIF
[1099]	493	!
	494	! ! =========================== !
	495	IF( nclosea == 0 ) THEN ! NO closed seas or lakes !
	496	DO jl = 1, jpncs ! =========================== !
[3]	497	DO jj = ncsj1(jl), ncsj2(jl)
	498	DO ji = ncsi1(jl), ncsi2(jl)
[1099]	499	mbathy(ji,jj) = 0 ! suppress closed seas and lakes from bathymetry
	500	bathy (ji,jj) = 0.e0
[3]	501	END DO
	502	END DO
	503	END DO
	504	ENDIF
[473]	505	#if defined key_orca_lev10
[1099]	506	! ! ================= !
	507	mbathy(:,:) = 10 * mbathy(:,:) ! key_orca_lev10 !
	508	! ! ================= !
	509	IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' CAUTION: 300 levels only with level bathymetry' )
[473]	510	#endif
[1099]	511	! ! =============== !
	512	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain !
	513	! ! =============== !
[3]	514
[2236]	515	#if ! defined key_c1d
	516	! ! =================== !
	517	CALL zgr_bat_ctl ! Bathymetry check !
	518	! ! =================== !
	519	#endif
[3]	520	END SUBROUTINE zgr_bat
	521
	522
	523	SUBROUTINE zgr_bat_zoom
	524	!!----------------------------------------------------------------------
	525	!! * ROUTINE zgr_bat_zoom *
	526	!!
	527	!! ** Purpose : - Close zoom domain boundary if necessary
	528	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
	529	!!
	530	!! ** Method :
	531	!!
	532	!! ** Action : - update mbathy: level bathymetry (in level index)
	533	!!----------------------------------------------------------------------
	534	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
	535	!!----------------------------------------------------------------------
[1099]	536	!
[3]	537	IF(lwp) WRITE(numout,*)
	538	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
	539	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
[1099]	540	!
[3]	541	! Zoom domain
	542	! ===========
[1099]	543	!
[3]	544	! Forced closed boundary if required
[1099]	545	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
	546	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
	547	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
	548	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
	549	!
[3]	550	! Configuration specific domain modifications
	551	! (here, ORCA arctic configuration: suppress Med Sea)
	552	IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
	553	SELECT CASE ( jp_cfg )
	554	! ! =======================
	555	CASE ( 2 ) ! ORCA_R2 configuration
	556	! ! =======================
	557	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
	558	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
	559	ij0 = 98 ; ij1 = 110
	560	! ! =======================
	561	CASE ( 05 ) ! ORCA_R05 configuration
	562	! ! =======================
	563	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
	564	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
	565	ij0 = 314 ; ij1 = 370
	566	END SELECT
	567	!
	568	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
	569	!
	570	ENDIF
[1099]	571	!
[3]	572	END SUBROUTINE zgr_bat_zoom
	573
	574
	575	SUBROUTINE zgr_bat_ctl
	576	!!----------------------------------------------------------------------
	577	!! * ROUTINE zgr_bat_ctl *
	578	!!
	579	!! ** Purpose : check the bathymetry in levels
	580	!!
	581	!! ** Method : The array mbathy is checked to verified its consistency
	582	!! with the model options. in particular:
	583	!! mbathy must have at least 1 land grid-points (mbathy<=0)
	584	!! along closed boundary.
	585	!! mbathy must be cyclic IF jperio=1.
	586	!! mbathy must be lower or equal to jpk-1.
	587	!! isolated ocean grid points are suppressed from mbathy
	588	!! since they are only connected to remaining
	589	!! ocean through vertical diffusion.
	590	!! C A U T I O N : mbathy will be modified during the initializa-
	591	!! tion phase to become the number of non-zero w-levels of a water
	592	!! column, with a minimum value of 1.
	593	!!
	594	!! ** Action : - update mbathy: level bathymetry (in level index)
	595	!! - update bathy : meter bathymetry (in meters)
	596	!!----------------------------------------------------------------------
[1099]	597	INTEGER :: ji, jj, jl ! dummy loop indices
	598	INTEGER :: icompt, ibtest, ikmax ! temporary integers
	599	REAL(wp), DIMENSION(jpi,jpj) :: zbathy ! temporary workspace
[3]	600	!!----------------------------------------------------------------------
	601
	602	IF(lwp) WRITE(numout,*)
	603	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
	604	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
	605
[1099]	606	! ! Suppress isolated ocean grid points
	607	IF(lwp) WRITE(numout,*)
	608	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
	609	IF(lwp) WRITE(numout,*)' -----------------------------------'
	610	icompt = 0
	611	DO jl = 1, 2
	612	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
	613	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
	614	mbathy(jpi,:) = mbathy( 2 ,:)
	615	ENDIF
	616	DO jj = 2, jpjm1
	617	DO ji = 2, jpim1
	618	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
	619	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
	620	IF( ibtest < mbathy(ji,jj) ) THEN
	621	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
	622	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
	623	mbathy(ji,jj) = ibtest
	624	icompt = icompt + 1
	625	ENDIF
	626	END DO
	627	END DO
	628	END DO
	629	IF( icompt == 0 ) THEN
	630	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
	631	ELSE
	632	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
	633	ENDIF
	634	IF( lk_mpp ) THEN
	635	zbathy(:,:) = FLOAT( mbathy(:,:) )
	636	CALL lbc_lnk( zbathy, 'T', 1. )
	637	mbathy(:,:) = INT( zbathy(:,:) )
	638	ENDIF
[3]	639
[1099]	640	! ! East-west cyclic boundary conditions
	641	IF( nperio == 0 ) THEN
	642	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
	643	IF( lk_mpp ) THEN
	644	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
	645	IF( jperio /= 1 ) mbathy(1,:) = 0
[411]	646	ENDIF
[1099]	647	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
	648	IF( jperio /= 1 ) mbathy(nlci,:) = 0
	649	ENDIF
[411]	650	ELSE
[1099]	651	IF( ln_zco .OR. ln_zps ) THEN
	652	mbathy( 1 ,:) = 0
	653	mbathy(jpi,:) = 0
	654	ELSE
	655	mbathy( 1 ,:) = jpkm1
	656	mbathy(jpi,:) = jpkm1
	657	ENDIF
[411]	658	ENDIF
[1099]	659	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
	660	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
	661	mbathy( 1 ,:) = mbathy(jpim1,:)
	662	mbathy(jpi,:) = mbathy( 2 ,:)
	663	ELSEIF( nperio == 2 ) THEN
	664	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
	665	ELSE
	666	IF(lwp) WRITE(numout,*) ' e r r o r'
	667	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
	668	! STOP 'dom_mba'
	669	ENDIF
	670
[1528]	671	! Boundary condition on mbathy
	672	IF( .NOT.lk_mpp ) THEN
	673	!!gm !!bug ??? think about it !
	674	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
	675	zbathy(:,:) = FLOAT( mbathy(:,:) )
	676	CALL lbc_lnk( zbathy, 'T', 1. )
	677	mbathy(:,:) = INT( zbathy(:,:) )
[3]	678	ENDIF
	679
	680	! Number of ocean level inferior or equal to jpkm1
	681	ikmax = 0
	682	DO jj = 1, jpj
	683	DO ji = 1, jpi
	684	ikmax = MAX( ikmax, mbathy(ji,jj) )
	685	END DO
	686	END DO
[1099]	687	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
[3]	688	IF( ikmax > jpkm1 ) THEN
	689	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
	690	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
	691	ELSE IF( ikmax < jpkm1 ) THEN
	692	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
	693	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
	694	ENDIF
	695
[1566]	696	IF( lwp .AND. nprint == 1 ) THEN ! control print
[3]	697	WRITE(numout,*)
[1099]	698	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
[3]	699	WRITE(numout,*) ' ------------------'
[1099]	700	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
[3]	701	WRITE(numout,*)
	702	ENDIF
[1099]	703	!
[3]	704	END SUBROUTINE zgr_bat_ctl
	705
	706
[454]	707	SUBROUTINE zgr_zco
	708	!!----------------------------------------------------------------------
	709	!! * ROUTINE zgr_zco *
	710	!!
	711	!! ** Purpose : define the z-coordinate system
	712	!!
[2240]	713	!! ** Method : set 3D coord. arrays to reference 1D array
[454]	714	!!----------------------------------------------------------------------
	715	INTEGER :: jk
	716	!!----------------------------------------------------------------------
[1099]	717	!
[2240]	718	DO jk = 1, jpk
	719	fsdept(:,:,jk) = gdept_0(jk)
	720	fsdepw(:,:,jk) = gdepw_0(jk)
	721	fsde3w(:,:,jk) = gdepw_0(jk)
	722	fse3t (:,:,jk) = e3t_0(jk)
	723	fse3u (:,:,jk) = e3t_0(jk)
	724	fse3v (:,:,jk) = e3t_0(jk)
	725	fse3f (:,:,jk) = e3t_0(jk)
	726	fse3w (:,:,jk) = e3w_0(jk)
	727	fse3uw(:,:,jk) = e3w_0(jk)
	728	fse3vw(:,:,jk) = e3w_0(jk)
	729	END DO
[1099]	730	!
[454]	731	END SUBROUTINE zgr_zco
	732
[1083]	733	!!----------------------------------------------------------------------
	734	!! Default option : zco, zps and/or sco available (gedp & e3 are 3D arrays)
	735	!!----------------------------------------------------------------------
[454]	736
[1083]	737	SUBROUTINE zgr_zps
	738	!!----------------------------------------------------------------------
	739	!! * ROUTINE zgr_zps *
	740	!!
	741	!! ** Purpose : the depth and vertical scale factor in partial step
	742	!! z-coordinate case
	743	!!
	744	!! ** Method : Partial steps : computes the 3D vertical scale factors
	745	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
	746	!! a partial step representation of bottom topography.
	747	!!
	748	!! The reference depth of model levels is defined from an analytical
	749	!! function the derivative of which gives the reference vertical
	750	!! scale factors.
	751	!! From depth and scale factors reference, we compute there new value
	752	!! with partial steps on 3d arrays ( i, j, k ).
	753	!!
	754	!! w-level: gdepw(i,j,k) = fsdep(k)
	755	!! e3w(i,j,k) = dk(fsdep)(k) = fse3(i,j,k)
	756	!! t-level: gdept(i,j,k) = fsdep(k+0.5)
	757	!! e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
	758	!!
	759	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
	760	!! we find the mbathy index of the depth at each grid point.
	761	!! This leads us to three cases:
	762	!!
	763	!! - bathy = 0 => mbathy = 0
	764	!! - 1 < mbathy < jpkm1
	765	!! - bathy > gdepw(jpk) => mbathy = jpkm1
	766	!!
	767	!! Then, for each case, we find the new depth at t- and w- levels
	768	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
	769	!! and f-points.
	770	!!
	771	!! This routine is given as an example, it must be modified
	772	!! following the user s desiderata. nevertheless, the output as
	773	!! well as the way to compute the model levels and scale factors
	774	!! must be respected in order to insure second order accuracy
	775	!! schemes.
	776	!!
	777	!! c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
	778	!! - - - - - - - gdept, gdepw and e3. are positives
	779	!!
[1099]	780	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
[1083]	781	!!----------------------------------------------------------------------
[1099]	782	INTEGER :: ji, jj, jk ! dummy loop indices
	783	INTEGER :: ik, it ! temporary integers
[2379]	784	INTEGER :: nlbelow ! temporary integer
[1099]	785	LOGICAL :: ll_print ! Allow control print for debugging
	786	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
	787	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
	788	REAL(wp) :: zmax, zmin ! Maximum and minimum depth
	789	REAL(wp) :: zdiff ! temporary scalar
[2379]	790	REAL(wp) :: zrefdep ! temporary scalar
[1099]	791	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zprt ! 3D workspace
	792	!!---------------------------------------------------------------------
[3]	793
[1099]	794	IF(lwp) WRITE(numout,*)
	795	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
	796	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
	797	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
[3]	798
[1099]	799	ll_print=.FALSE. ! Local variable for debugging
	800	!! ll_print=.TRUE.
[1083]	801
[1099]	802	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
[1083]	803	WRITE(numout,*)
	804	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
	805	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
	806	ENDIF
	807
	808
	809	! bathymetry in level (from bathy_meter)
	810	! ===================
[1099]	811	zmax = gdepw_0(jpk) + e3t_0(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_0(jpkm1) )
[1083]	812
[2379]	813	zrefdep = 25.
	814	nlbelow = MINLOC( gdepw_0, mask = gdepw_0 > zrefdep, dim = 1 ) ! shallowest W level Below zrefdep
	815	IF(lwp) write(numout,*) 'Minimum depth level selected: ', nlbelow
	816	zmin = gdepw_0(nlbelow) ! minimum depth = nlbelow-1 levels
	817
[1099]	818	mbathy(:,:) = jpkm1 ! initialize mbathy to the maximum ocean level available
[1083]	819
[1099]	820	! ! storage of land and island's number (zera and negative values) in mbathy
	821	WHERE( bathy(:,:) <= 0. ) mbathy(:,:) = INT( bathy(:,:) )
[1083]	822
	823	! bounded value of bathy
[1099]	824	!!gm bathy(:,:) = MIN( zmax, MAX( bathy(:,:), zmin ) ) !!gm : simpler as zmin is > 0
[1083]	825	DO jj = 1, jpj
	826	DO ji= 1, jpi
[1099]	827	IF( bathy(ji,jj) <= 0. ) THEN ; bathy(ji,jj) = 0.e0
	828	ELSE ; bathy(ji,jj) = MIN( zmax, MAX( bathy(ji,jj), zmin ) )
[1083]	829	ENDIF
	830	END DO
	831	END DO
	832
	833	! Compute mbathy for ocean points (i.e. the number of ocean levels)
	834	! find the number of ocean levels such that the last level thickness
	835	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
	836	! e3t_0 is the reference level thickness
	837	DO jk = jpkm1, 1, -1
	838	zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
	839	DO jj = 1, jpj
	840	DO ji = 1, jpi
	841	IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth ) mbathy(ji,jj) = jk-1
	842	END DO
	843	END DO
	844	END DO
	845
[1099]	846	! Scale factors and depth at T- and W-points
	847	DO jk = 1, jpk ! intitialization to the reference z-coordinate
	848	gdept(:,:,jk) = gdept_0(jk)
	849	gdepw(:,:,jk) = gdepw_0(jk)
	850	e3t (:,:,jk) = e3t_0 (jk)
	851	e3w (:,:,jk) = e3w_0 (jk)
[1083]	852	END DO
[1099]	853	hdept(:,:) = gdept(:,:,2 )
	854	hdepw(:,:) = gdepw(:,:,3 )
	855	!
	856	DO jj = 1, jpj
	857	DO ji = 1, jpi
	858	ik = mbathy(ji,jj)
	859	IF( ik > 0 ) THEN ! ocean point only
	860	! max ocean level case
	861	IF( ik == jpkm1 ) THEN
	862	zdepwp = bathy(ji,jj)
	863	ze3tp = bathy(ji,jj) - gdepw_0(ik)
	864	ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
	865	e3t(ji,jj,ik ) = ze3tp
	866	e3t(ji,jj,ik+1) = ze3tp
	867	e3w(ji,jj,ik ) = ze3wp
	868	e3w(ji,jj,ik+1) = ze3tp
	869	gdepw(ji,jj,ik+1) = zdepwp
	870	gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp
	871	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
	872	!
	873	ELSE ! standard case
	874	IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN ; gdepw(ji,jj,ik+1) = bathy(ji,jj)
	875	ELSE ; gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
	876	ENDIF
	877	!gm Bug? check the gdepw_0
	878	! ... on ik
	879	gdept(ji,jj,ik) = gdepw_0(ik) + ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
	880	& * ((gdept_0( ik ) - gdepw_0(ik) ) &
	881	& / ( gdepw_0( ik+1) - gdepw_0(ik) ))
	882	e3t (ji,jj,ik) = e3t_0 (ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
	883	& / ( gdepw_0( ik+1) - gdepw_0(ik) )
	884	e3w (ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) &
	885	& * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
	886	! ... on ik+1
	887	e3w (ji,jj,ik+1) = e3t (ji,jj,ik)
	888	e3t (ji,jj,ik+1) = e3t (ji,jj,ik)
	889	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
	890	ENDIF
	891	ENDIF
	892	END DO
	893	END DO
	894	!
	895	it = 0
	896	DO jj = 1, jpj
	897	DO ji = 1, jpi
	898	ik = mbathy(ji,jj)
	899	IF( ik > 0 ) THEN ! ocean point only
	900	hdept(ji,jj) = gdept(ji,jj,ik )
	901	hdepw(ji,jj) = gdepw(ji,jj,ik+1)
	902	e3tp (ji,jj) = e3t(ji,jj,ik )
	903	e3wp (ji,jj) = e3w(ji,jj,ik )
	904	! test
	905	zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik )
	906	IF( zdiff <= 0. .AND. lwp ) THEN
	907	it = it + 1
	908	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
	909	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
	910	WRITE(numout,*) ' gdept = ', gdept(ji,jj,ik), ' gdepw = ', gdepw(ji,jj,ik+1), ' zdiff = ', zdiff
	911	WRITE(numout,*) ' e3tp = ', e3t (ji,jj,ik), ' e3wp = ', e3w (ji,jj,ik )
	912	ENDIF
	913	ENDIF
	914	END DO
	915	END DO
[1083]	916
	917	! Scale factors and depth at U-, V-, UW and VW-points
[1099]	918	DO jk = 1, jpk ! initialisation to z-scale factors
[1083]	919	e3u (:,:,jk) = e3t_0(jk)
	920	e3v (:,:,jk) = e3t_0(jk)
	921	e3uw(:,:,jk) = e3w_0(jk)
	922	e3vw(:,:,jk) = e3w_0(jk)
	923	END DO
[1099]	924	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
	925	DO jj = 1, jpjm1
	926	DO ji = 1, fs_jpim1 ! vector opt.
	927	e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
	928	e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
	929	e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
	930	e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
	931	END DO
	932	END DO
	933	END DO
	934	! ! Boundary conditions
	935	CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. )
	936	CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. )
	937	!
	938	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
	939	WHERE( e3u (:,:,jk) == 0.e0 ) e3u (:,:,jk) = e3t_0(jk)
	940	WHERE( e3v (:,:,jk) == 0.e0 ) e3v (:,:,jk) = e3t_0(jk)
	941	WHERE( e3uw(:,:,jk) == 0.e0 ) e3uw(:,:,jk) = e3w_0(jk)
	942	WHERE( e3vw(:,:,jk) == 0.e0 ) e3vw(:,:,jk) = e3w_0(jk)
	943	END DO
	944
	945	! Scale factor at F-point
	946	DO jk = 1, jpk ! initialisation to z-scale factors
	947	e3f (:,:,jk) = e3t_0(jk)
	948	END DO
	949	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
	950	DO jj = 1, jpjm1
	951	DO ji = 1, fs_jpim1 ! vector opt.
	952	e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
	953	END DO
	954	END DO
	955	END DO
	956	CALL lbc_lnk( e3f, 'F', 1. ) ! Boundary conditions
	957	!
	958	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
	959	WHERE( e3f(:,:,jk) == 0.e0 ) e3f(:,:,jk) = e3t_0(jk)
	960	END DO
	961	!!gm bug ? : must be a do loop with mj0,mj1
	962	!
	963	e3t(:,mj0(1),:) = e3t(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
	964	e3w(:,mj0(1),:) = e3w(:,mj0(2),:)
	965	e3u(:,mj0(1),:) = e3u(:,mj0(2),:)
	966	e3v(:,mj0(1),:) = e3v(:,mj0(2),:)
	967	e3f(:,mj0(1),:) = e3f(:,mj0(2),:)
[1083]	968
[1099]	969	! Control of the sign
	970	IF( MINVAL( e3t (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3t <= 0' )
	971	IF( MINVAL( e3w (:,:,:) ) <= 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r e3w <= 0' )
	972	IF( MINVAL( gdept(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
	973	IF( MINVAL( gdepw(:,:,:) ) < 0.e0 ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw < 0' )
[1083]	974
[1099]	975	! Compute gdep3w (vertical sum of e3w)
	976	gdep3w(:,:,1) = 0.5 * e3w(:,:,1)
	977	DO jk = 2, jpk
	978	gdep3w(:,:,jk) = gdep3w(:,:,jk-1) + e3w(:,:,jk)
	979	END DO
[1083]	980
[1099]	981	! ! ================= !
	982	IF(lwp .AND. ll_print) THEN ! Control print !
	983	! ! ================= !
[1083]	984	DO jj = 1,jpj
	985	DO ji = 1, jpi
[1099]	986	ik = MAX( mbathy(ji,jj), 1 )
	987	zprt(ji,jj,1) = e3t (ji,jj,ik)
	988	zprt(ji,jj,2) = e3w (ji,jj,ik)
	989	zprt(ji,jj,3) = e3u (ji,jj,ik)
	990	zprt(ji,jj,4) = e3v (ji,jj,ik)
	991	zprt(ji,jj,5) = e3f (ji,jj,ik)
	992	zprt(ji,jj,6) = gdep3w(ji,jj,ik)
[1083]	993	END DO
	994	END DO
	995	WRITE(numout,*)
[1099]	996	WRITE(numout,*) 'domzgr e3t(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 e3w(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 e3u(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 e3v(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1003	WRITE(numout,*)
[1099]	1004	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1005	WRITE(numout,*)
[1099]	1006	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
[1083]	1007	ENDIF
[1099]	1008
	1009	! ! =============== !
	1010	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain !
	1011	! ! =============== !
[1083]	1012
[2236]	1013	#if ! defined key_c1d
	1014	! ! =================== !
	1015	CALL zgr_bat_ctl ! Bathymetry check !
	1016	! ! =================== !
	1017	#endif
[1083]	1018	END SUBROUTINE zgr_zps
	1019
	1020
	1021	FUNCTION fssig( pk ) RESULT( pf )
	1022	!!----------------------------------------------------------------------
	1023	!! * ROUTINE eos_init *
	1024	!!
	1025	!! ** Purpose : provide the analytical function in s-coordinate
	1026	!!
	1027	!! ** Method : the function provide the non-dimensional position of
	1028	!! T and W (i.e. between 0 and 1)
	1029	!! T-points at integer values (between 1 and jpk)
	1030	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1031	!!
	1032	!! Reference : ???
	1033	!!----------------------------------------------------------------------
	1034	REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate
	1035	REAL(wp) :: pf ! sigma value
	1036	!!----------------------------------------------------------------------
	1037	!
[1601]	1038	pf = ( TANH( rn_theta * ( -(pk-0.5) / REAL(jpkm1) + rn_thetb ) ) &
	1039	& - TANH( rn_thetb * rn_theta ) ) &
	1040	& * ( COSH( rn_theta ) &
	1041	& + COSH( rn_theta * ( 2.e0 * rn_thetb - 1.e0 ) ) ) &
	1042	& / ( 2.e0 * SINH( rn_theta ) )
[1083]	1043	!
	1044	END FUNCTION fssig
	1045
	1046
[1601]	1047	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
[1348]	1048	!!----------------------------------------------------------------------
	1049	!! * ROUTINE eos_init *
	1050	!!
	1051	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
	1052	!!
	1053	!! ** Method : the function provides the non-dimensional position of
	1054	!! T and W (i.e. between 0 and 1)
	1055	!! T-points at integer values (between 1 and jpk)
	1056	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
	1057	!!
	1058	!! Reference : ???
	1059	!!----------------------------------------------------------------------
	1060	REAL(wp), INTENT(in ) :: pk1 ! continuous "k" coordinate
[1601]	1061	REAL(wp), INTENT(in ) :: pbb ! Stretching coefficient
[1348]	1062	REAL(wp) :: pf1 ! sigma value
	1063	!!----------------------------------------------------------------------
	1064	!
[1601]	1065	IF ( rn_theta == 0 ) then ! uniform sigma
[1566]	1066	pf1 = -(pk1-0.5) / REAL( jpkm1 )
	1067	ELSE ! stretched sigma
[1601]	1068	pf1 = (1.0-pbb) * (sinh( rn_theta*(-(pk1-0.5)/REAL(jpkm1)) ) ) / sinh(rn_theta) + &
	1069	& pbb * ( (tanh( rn_theta( (-(pk1-0.5)/REAL(jpkm1)) + 0.5) ) - tanh(0.5rn_theta) ) / &
	1070	& (2tanh(0.5rn_theta) ) )
[1348]	1071	ENDIF
[1566]	1072	!
[1348]	1073	END FUNCTION fssig1
	1074
	1075
[454]	1076	SUBROUTINE zgr_sco
	1077	!!----------------------------------------------------------------------
	1078	!! * ROUTINE zgr_sco *
	1079	!!
	1080	!! ** Purpose : define the s-coordinate system
	1081	!!
	1082	!! ** Method : s-coordinate
	1083	!! The depth of model levels is defined as the product of an
	1084	!! analytical function by the local bathymetry, while the vertical
	1085	!! scale factors are defined as the product of the first derivative
	1086	!! of the analytical function by the bathymetry.
	1087	!! (this solution save memory as depth and scale factors are not
	1088	!! 3d fields)
	1089	!! - Read bathymetry (in meters) at t-point and compute the
	1090	!! bathymetry at u-, v-, and f-points.
	1091	!! hbatu = mi( hbatt )
	1092	!! hbatv = mj( hbatt )
	1093	!! hbatf = mi( mj( hbatt ) )
	1094	!! - Compute gsigt, gsigw, esigt, esigw from an analytical
[1083]	1095	!! function and its derivative given as function.
	1096	!! gsigt(k) = fssig (k )
	1097	!! gsigw(k) = fssig (k-0.5)
	1098	!! esigt(k) = fsdsig(k )
	1099	!! esigw(k) = fsdsig(k-0.5)
[454]	1100	!! This routine is given as an example, it must be modified
	1101	!! following the user s desiderata. nevertheless, the output as
	1102	!! well as the way to compute the model levels and scale factors
	1103	!! must be respected in order to insure second order a!!uracy
	1104	!! schemes.
	1105	!!
[1099]	1106	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
	1107	!!----------------------------------------------------------------------
	1108	INTEGER :: ji, jj, jk, jl ! dummy loop argument
	1109	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
	1110	REAL(wp) :: zcoeft, zcoefw, zrmax, ztaper ! temporary scalars
	1111	REAL(wp), DIMENSION(jpi,jpj) :: zenv, ztmp, zmsk ! 2D workspace
	1112	REAL(wp), DIMENSION(jpi,jpj) :: zri , zrj , zhbat ! - -
[1566]	1113	!!
[1629]	1114	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigw3
	1115	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsigt3
	1116	REAL(wp), DIMENSION(jpi,jpj,jpk) :: gsi3w3
	1117	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigt3
	1118	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigw3
	1119	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtu3
	1120	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtv3
	1121	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigtf3
	1122	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigwu3
	1123	REAL(wp), DIMENSION(jpi,jpj,jpk) :: esigwv3
[454]	1124	!!
[1601]	1125	NAMELIST/namzgr_sco/ rn_sbot_max, rn_sbot_min, rn_theta, rn_thetb, rn_rmax, ln_s_sigma, rn_bb, rn_hc
[454]	1126	!!----------------------------------------------------------------------
	1127
[1601]	1128	REWIND( numnam ) ! Read Namelist namzgr_sco : sigma-stretching parameters
	1129	READ ( numnam, namzgr_sco )
[454]	1130
[1099]	1131	IF(lwp) THEN ! control print
[454]	1132	WRITE(numout,*)
	1133	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
	1134	WRITE(numout,*) '~~~~~~~~~~~'
[1601]	1135	WRITE(numout,*) ' Namelist namzgr_sco'
	1136	WRITE(numout,*) ' sigma-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,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ', rn_theta
	1140	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ', rn_thetb
	1141	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ', rn_rmax
	1142	WRITE(numout,*) ' Hybrid s-sigma-coordinate ln_s_sigma = ', ln_s_sigma
	1143	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ', rn_bb
	1144	WRITE(numout,*) ' Critical depth rn_hc = ', rn_hc
[454]	1145	ENDIF
	1146
[1629]	1147	gsigw3 = 0.0d0 ; gsigt3 = 0.0d0 ; gsi3w3 = 0.0d0
	1148	esigt3 = 0.0d0 ; esigw3 = 0.0d0
	1149	esigtu3 = 0.0d0 ; esigtv3 = 0.0d0 ; esigtf3 = 0.0d0
	1150	esigwu3 = 0.0d0 ; esigwv3 = 0.0d0
	1151
[1601]	1152	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
	1153	hifu(:,:) = rn_sbot_min
	1154	hifv(:,:) = rn_sbot_min
	1155	hiff(:,:) = rn_sbot_min
[1348]	1156
	1157	! ! set maximum ocean depth
[1601]	1158	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
[454]	1159
[1461]	1160	DO jj = 1, jpj
	1161	DO ji = 1, jpi
	1162	IF (bathy(ji,jj) .gt. 0.0) THEN
[1601]	1163	bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
[1461]	1164	ENDIF
	1165	END DO
	1166	END DO
[1099]	1167	! ! =============================
	1168	! ! Define the envelop bathymetry (hbatt)
	1169	! ! =============================
[454]	1170	! use r-value to create hybrid coordinates
	1171	DO jj = 1, jpj
	1172	DO ji = 1, jpi
[1601]	1173	zenv(ji,jj) = MAX( bathy(ji,jj), rn_sbot_min )
[454]	1174	END DO
	1175	END DO
[1639]	1176	!
	1177	! Smooth the bathymetry (if required)
	1178	scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea)
	1179	scobot(:,:) = bathy(:,:) ! ocean bottom depth
	1180	!
[454]	1181	jl = 0
	1182	zrmax = 1.e0
	1183	! ! ================ !
[1601]	1184	DO WHILE ( jl <= 10000 .AND. zrmax > rn_rmax ) ! Iterative loop !
[454]	1185	! ! ================ !
	1186	jl = jl + 1
	1187	zrmax = 0.e0
	1188	zmsk(:,:) = 0.e0
	1189	DO jj = 1, nlcj
	1190	DO ji = 1, nlci
	1191	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
	1192	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
	1193	zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
	1194	zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
	1195	zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
[1601]	1196	IF( zri(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1.0
	1197	IF( zri(ji,jj) > rn_rmax ) zmsk(iip1,jj ) = 1.0
	1198	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,jj ) = 1.0
	1199	IF( zrj(ji,jj) > rn_rmax ) zmsk(ji ,ijp1) = 1.0
[454]	1200	END DO
	1201	END DO
	1202	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
	1203	! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
[1099]	1204	ztmp(:,:) = zmsk(:,:) ; CALL lbc_lnk( zmsk, 'T', 1. )
[454]	1205	DO jj = 1, nlcj
	1206	DO ji = 1, nlci
[1348]	1207	zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
[454]	1208	END DO
	1209	END DO
[1348]	1210	!
[1099]	1211	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
	1212	!
[454]	1213	DO jj = 1, nlcj
	1214	DO ji = 1, nlci
	1215	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
	1216	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
	1217	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
	1218	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
	1219	IF( zmsk(ji,jj) == 1.0 ) THEN
	1220	ztmp(ji,jj) = ( &
	1221	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
	1222	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2.e0 + zenv(iip1,jj )*zmsk(iip1,jj ) &
	1223	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
	1224	& ) / ( &
	1225	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
	1226	& + zmsk(iim1,jj ) + 2.e0 + zmsk(iip1,jj ) &
	1227	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
	1228	& )
	1229	ENDIF
	1230	END DO
	1231	END DO
[1348]	1232	!
[454]	1233	DO jj = 1, nlcj
	1234	DO ji = 1, nlci
	1235	IF( zmsk(ji,jj) == 1.0 ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
	1236	END DO
	1237	END DO
[1099]	1238	!
[454]	1239	! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero
[1099]	1240	ztmp(:,:) = zenv(:,:) ; CALL lbc_lnk( zenv, 'T', 1. )
[454]	1241	DO jj = 1, nlcj
	1242	DO ji = 1, nlci
	1243	IF( zenv(ji,jj) == 0.e0 ) zenv(ji,jj) = ztmp(ji,jj)
	1244	END DO
	1245	END DO
	1246	! ! ================ !
	1247	END DO ! End loop !
	1248	! ! ================ !
[1099]	1249	!
	1250	! ! envelop bathymetry saved in hbatt
[454]	1251	hbatt(:,:) = zenv(:,:)
[1099]	1252	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0.e0 ) THEN
	1253	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
	1254	DO jj = 1, jpj
	1255	DO ji = 1, jpi
	1256	ztaper = EXP( -(gphit(ji,jj)/8)**2 )
[1601]	1257	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
[1099]	1258	END DO
	1259	END DO
[516]	1260	ENDIF
[1099]	1261	!
	1262	IF(lwp) THEN ! Control print
[454]	1263	WRITE(numout,*)
	1264	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
	1265	WRITE(numout,*)
	1266	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
[1099]	1267	IF( nprint == 1 ) THEN
	1268	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
	1269	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
	1270	ENDIF
[454]	1271	ENDIF
	1272
[1099]	1273	! ! ==============================
	1274	! ! hbatu, hbatv, hbatf fields
	1275	! ! ==============================
[454]	1276	IF(lwp) THEN
	1277	WRITE(numout,*)
[1601]	1278	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
[454]	1279	ENDIF
[1601]	1280	hbatu(:,:) = rn_sbot_min
	1281	hbatv(:,:) = rn_sbot_min
	1282	hbatf(:,:) = rn_sbot_min
[454]	1283	DO jj = 1, jpjm1
[1694]	1284	DO ji = 1, jpim1 ! NO vector opt.
[1099]	1285	hbatu(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
	1286	hbatv(ji,jj) = 0.5 * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
	1287	hbatf(ji,jj) = 0.25* ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
	1288	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
[454]	1289	END DO
	1290	END DO
[1099]	1291	!
[454]	1292	! Apply lateral boundary condition
[1099]	1293	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
	1294	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk( hbatu, 'U', 1. )
[454]	1295	DO jj = 1, jpj
	1296	DO ji = 1, jpi
	1297	IF( hbatu(ji,jj) == 0.e0 ) THEN
[1601]	1298	IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = rn_sbot_min
[454]	1299	IF( zhbat(ji,jj) /= 0.e0 ) hbatu(ji,jj) = zhbat(ji,jj)
	1300	ENDIF
	1301	END DO
	1302	END DO
[1099]	1303	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk( hbatv, 'V', 1. )
[454]	1304	DO jj = 1, jpj
	1305	DO ji = 1, jpi
	1306	IF( hbatv(ji,jj) == 0.e0 ) THEN
[1601]	1307	IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = rn_sbot_min
[454]	1308	IF( zhbat(ji,jj) /= 0.e0 ) hbatv(ji,jj) = zhbat(ji,jj)
	1309	ENDIF
	1310	END DO
	1311	END DO
[1099]	1312	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk( hbatf, 'F', 1. )
[454]	1313	DO jj = 1, jpj
	1314	DO ji = 1, jpi
	1315	IF( hbatf(ji,jj) == 0.e0 ) THEN
[1601]	1316	IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = rn_sbot_min
[454]	1317	IF( zhbat(ji,jj) /= 0.e0 ) hbatf(ji,jj) = zhbat(ji,jj)
	1318	ENDIF
	1319	END DO
	1320	END DO
	1321
	1322	!!bug: key_helsinki a verifer
	1323	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
	1324	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
	1325	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
	1326	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
	1327
[516]	1328	IF( nprint == 1 .AND. lwp ) THEN
[1099]	1329	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
	1330	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
	1331	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
	1332	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
[516]	1333	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
	1334	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
	1335	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
	1336	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
	1337	ENDIF
[454]	1338	!! helsinki
	1339
[1099]	1340	! ! =======================
	1341	! ! s-ccordinate fields (gdep., e3.)
	1342	! ! =======================
	1343	!
	1344	! non-dimensional "sigma" for model level depth at w- and t-levels
[1348]	1345
[1601]	1346	IF( ln_s_sigma ) THEN ! Song and Haidvogel style stretched sigma for depths
	1347	! ! below rn_hc, with uniform sigma in shallower waters
	1348	DO ji = 1, jpi
	1349	DO jj = 1, jpj
[1348]	1350
[1601]	1351	IF (hbatt(ji,jj).GT.rn_hc) THEN !deep water, stretched sigma
[1348]	1352	DO jk = 1, jpk
[1601]	1353	gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
	1354	gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
[1348]	1355	END DO
	1356	ELSE ! shallow water, uniform sigma
	1357	DO jk = 1, jpk
	1358	gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
	1359	gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
	1360	END DO
	1361	ENDIF
	1362	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw3 1 jpk ', gsigw3(ji,jj,1), gsigw3(ji,jj,jpk)
	1363
	1364
	1365	DO jk = 1, jpkm1
	1366	esigt3(ji,jj,jk) = gsigw3(ji,jj,jk+1) - gsigw3(ji,jj,jk)
	1367	esigw3(ji,jj,jk+1) = gsigt3(ji,jj,jk+1) - gsigt3(ji,jj,jk)
	1368	END DO
[1639]	1369	esigw3(ji,jj,1 ) = 2.0_wp * (gsigt3(ji,jj,1 ) - gsigw3(ji,jj,1 ))
	1370	esigt3(ji,jj,jpk) = 2.0_wp * (gsigt3(ji,jj,jpk) - gsigw3(ji,jj,jpk))
[1348]	1371
	1372	! Coefficients for vertical depth as the sum of e3w scale factors
	1373	gsi3w3(ji,jj,1) = 0.5 * esigw3(ji,jj,1)
	1374	DO jk = 2, jpk
	1375	gsi3w3(ji,jj,jk) = gsi3w3(ji,jj,jk-1) + esigw3(ji,jj,jk)
	1376	END DO
	1377
	1378	DO jk = 1, jpk
	1379	zcoeft = ( REAL(jk,wp) - 0.5 ) / REAL(jpkm1,wp)
	1380	zcoefw = ( REAL(jk,wp) - 1.0 ) / REAL(jpkm1,wp)
[1601]	1381	gdept (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsigt3(ji,jj,jk)+rn_hczcoeft)
	1382	gdepw (ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsigw3(ji,jj,jk)+rn_hczcoefw)
[1639]	1383	gdep3w(ji,jj,jk) = (scosrf(ji,jj)+(hbatt(ji,jj)-rn_hc)gsi3w3(ji,jj,jk)+rn_hczcoeft)
[1348]	1384	END DO
	1385
	1386	ENDDO ! for all jj's
	1387	ENDDO ! for all ji's
	1388
	1389
	1390	DO ji=1,jpi
	1391	DO jj=1,jpj
	1392
	1393	DO jk = 1, jpk
	1394	esigtu3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji+1,jj)esigt3(ji+1,jj,jk) ) / &
	1395	( hbatt(ji,jj)+hbatt(ji+1,jj) )
	1396	esigtv3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji,jj+1)esigt3(ji,jj+1,jk) ) / &
	1397	( hbatt(ji,jj)+hbatt(ji,jj+1) )
	1398	esigtf3(ji,jj,jk) = ( hbatt(ji,jj)esigt3(ji,jj,jk)+hbatt(ji+1,jj)esigt3(ji+1,jj,jk) + &
	1399	hbatt(ji,jj+1)esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)esigt3(ji+1,jj+1,jk) ) / &
	1400	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
	1401	esigwu3(ji,jj,jk) = ( hbatt(ji,jj)esigw3(ji,jj,jk)+hbatt(ji+1,jj)esigw3(ji+1,jj,jk) ) / &
	1402	( hbatt(ji,jj)+hbatt(ji+1,jj) )
	1403	esigwv3(ji,jj,jk) = ( hbatt(ji,jj)esigw3(ji,jj,jk)+hbatt(ji,jj+1)esigw3(ji,jj+1,jk) ) / &
	1404	( hbatt(ji,jj)+hbatt(ji,jj+1) )
	1405
[1601]	1406	e3t(ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigt3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
	1407	e3u(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigtu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
	1408	e3v(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigtv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
	1409	e3f(ji,jj,jk)=((hbatf(ji,jj)-rn_hc)*esigtf3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
[1348]	1410	!
[1601]	1411	e3w (ji,jj,jk)=((hbatt(ji,jj)-rn_hc)*esigw3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
	1412	e3uw(ji,jj,jk)=((hbatu(ji,jj)-rn_hc)*esigwu3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
	1413	e3vw(ji,jj,jk)=((hbatv(ji,jj)-rn_hc)*esigwv3(ji,jj,jk) + rn_hc/FLOAT(jpkm1))
[1348]	1414	END DO
	1415
	1416	ENDDO
	1417	ENDDO
	1418
	1419	ELSE ! not ln_s_sigma
	1420
	1421	DO jk = 1, jpk
	1422	gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
	1423	gsigt(jk) = -fssig( REAL(jk,wp) )
	1424	END DO
	1425	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)
	1426	!
	1427	! Coefficients for vertical scale factors at w-, t- levels
[1099]	1428	!!gm bug : define it from analytical function, not like juste bellow....
	1429	!!gm or betteroffer the 2 possibilities....
[1348]	1430	DO jk = 1, jpkm1
	1431	esigt(jk ) = gsigw(jk+1) - gsigw(jk)
	1432	esigw(jk+1) = gsigt(jk+1) - gsigt(jk)
	1433	END DO
[1639]	1434	esigw( 1 ) = 2.0_wp * (gsigt(1) - gsigw(1))
	1435	esigt(jpk) = 2.0_wp * (gsigt(jpk) - gsigw(jpk))
[454]	1436
[1099]	1437	!!gm original form
	1438	!!org DO jk = 1, jpk
	1439	!!org esigt(jk)=fsdsig( FLOAT(jk) )
	1440	!!org esigw(jk)=fsdsig( FLOAT(jk)-0.5 )
	1441	!!org END DO
	1442	!!gm
[1348]	1443	!
	1444	! Coefficients for vertical depth as the sum of e3w scale factors
	1445	gsi3w(1) = 0.5 * esigw(1)
	1446	DO jk = 2, jpk
	1447	gsi3w(jk) = gsi3w(jk-1) + esigw(jk)
	1448	END DO
[1099]	1449	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
[1348]	1450	DO jk = 1, jpk
	1451	zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
	1452	zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
	1453	gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigt(jk)+hift(:,:)zcoeft)
	1454	gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigw(jk)+hift(:,:)zcoefw)
[1639]	1455	gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsi3w(jk)+hift(:,:)zcoeft)
[1348]	1456	END DO
[1099]	1457	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
[1348]	1458	DO jj = 1, jpj
	1459	DO ji = 1, jpi
	1460	DO jk = 1, jpk
	1461	e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
	1462	e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
	1463	e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
	1464	e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
	1465	!
	1466	e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
	1467	e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
	1468	e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
	1469	END DO
[454]	1470	END DO
	1471	END DO
[1348]	1472
	1473	ENDIF ! ln_s_sigma
	1474
	1475
[1099]	1476	!
[1461]	1477	!! H. Liu, POL. April 2009. Added for passing the scale check for the new released vvl code.
	1478
	1479	fsdept(:,:,:) = gdept (:,:,:)
	1480	fsdepw(:,:,:) = gdepw (:,:,:)
	1481	fsde3w(:,:,:) = gdep3w(:,:,:)
	1482	fse3t (:,:,:) = e3t (:,:,:)
	1483	fse3u (:,:,:) = e3u (:,:,:)
	1484	fse3v (:,:,:) = e3v (:,:,:)
	1485	fse3f (:,:,:) = e3f (:,:,:)
	1486	fse3w (:,:,:) = e3w (:,:,:)
	1487	fse3uw(:,:,:) = e3uw (:,:,:)
	1488	fse3vw(:,:,:) = e3vw (:,:,:)
	1489	!!
[1099]	1490	! HYBRID :
[454]	1491	DO jj = 1, jpj
	1492	DO ji = 1, jpi
	1493	DO jk = 1, jpkm1
	1494	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
	1495	IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0
	1496	END DO
	1497	END DO
	1498	END DO
[1099]	1499	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
	1500	& ' MAX ', MAXVAL( mbathy(:,:) )
[454]	1501
	1502
[1632]	1503	! ! ===========
	1504	IF( lzoom ) CALL zgr_bat_zoom ! Zoom domain
	1505	! ! ===========
[454]	1506
[2236]	1507	#if ! defined key_c1d
	1508	! ! =================== !
	1509	CALL zgr_bat_ctl ! Bathymetry check !
	1510	! ! =================== !
	1511	#endif
[1632]	1512
	1513	! ! =============
	1514	IF(lwp) THEN ! Control print
	1515	! ! =============
[454]	1516	WRITE(numout,*)
	1517	WRITE(numout,*) ' domzgr: vertical coefficients for model level'
[1099]	1518	WRITE(numout, "(9x,' level gsigt gsigw esigt esigw gsi3w')" )
	1519	WRITE(numout, "(10x,i4,5f11.4)" ) ( jk, gsigt(jk), gsigw(jk), esigt(jk), esigw(jk), gsi3w(jk), jk=1,jpk )
[454]	1520	ENDIF
[1099]	1521	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
	1522	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
	1523	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
	1524	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ' , MINVAL( fsde3w(:,:,:) )
	1525	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t (:,:,:) ), ' f ' , MINVAL( fse3f (:,:,:) ), &
	1526	& ' u ', MINVAL( fse3u (:,:,:) ), ' u ' , MINVAL( fse3v (:,:,:) ), &
	1527	& ' uw', MINVAL( fse3uw(:,:,:) ), ' vw' , MINVAL( fse3vw(:,:,:) ), &
	1528	& ' w ', MINVAL( fse3w (:,:,:) )
[454]	1529
[1099]	1530	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
	1531	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ' , MAXVAL( fsde3w(:,:,:) )
	1532	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t (:,:,:) ), ' f ' , MAXVAL( fse3f (:,:,:) ), &
	1533	& ' u ', MAXVAL( fse3u (:,:,:) ), ' u ' , MAXVAL( fse3v (:,:,:) ), &
	1534	& ' uw', MAXVAL( fse3uw(:,:,:) ), ' vw' , MAXVAL( fse3vw(:,:,:) ), &
	1535	& ' w ', MAXVAL( fse3w (:,:,:) )
	1536	ENDIF
	1537	!
	1538	IF(lwp) THEN ! selected vertical profiles
[454]	1539	WRITE(numout,*)
	1540	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
	1541	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1542	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1543	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(1,1,jk), fsdepw(1,1,jk), &
	1544	& fse3t (1,1,jk), fse3w (1,1,jk), jk=1,jpk )
[473]	1545	DO jj = mj0(20), mj1(20)
	1546	DO ji = mi0(20), mi1(20)
	1547	WRITE(numout,*)
	1548	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	1549	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1550	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1551	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
	1552	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
[473]	1553	END DO
	1554	END DO
	1555	DO jj = mj0(74), mj1(74)
	1556	DO ji = mi0(100), mi1(100)
	1557	WRITE(numout,*)
	1558	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
	1559	WRITE(numout,*) ' ~~~~~~ --------------------'
[1099]	1560	WRITE(numout,"(9x,' level gdept gdepw gde3w e3t e3w ')")
	1561	WRITE(numout,"(10x,i4,4f9.2)") ( jk, fsdept(ji,jj,jk), fsdepw(ji,jj,jk), &
	1562	& fse3t (ji,jj,jk), fse3w (ji,jj,jk), jk=1,jpk )
[473]	1563	END DO
	1564	END DO
[454]	1565	ENDIF
	1566
[1099]	1567	!!gm bug? no more necessary? if ! defined key_helsinki
[454]	1568	DO jk = 1, jpk
	1569	DO jj = 1, jpj
	1570	DO ji = 1, jpi
	1571	IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
[1099]	1572	WRITE(ctmp1,*) 'zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
	1573	CALL ctl_stop( ctmp1 )
[454]	1574	ENDIF
	1575	IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
[1099]	1576	WRITE(ctmp1,*) 'zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
	1577	CALL ctl_stop( ctmp1 )
[454]	1578	ENDIF
	1579	END DO
	1580	END DO
	1581	END DO
[1099]	1582	!!gm bug #endif
	1583	!
[454]	1584	END SUBROUTINE zgr_sco
	1585
	1586
[3]	1587	!!======================================================================
	1588	END MODULE domzgr

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

Download in other formats: