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

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

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

Download in other formats: