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

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

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

Download in other formats: