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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90 @ 2950

Last change on this file since 2950 was 2950, checked in by acc, 13 years ago
Fixed incorrect ORCA1 settings in domzgr. See ticket #875
Property svn:keywords set to `Id`
File size: 81.3 KB

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

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

Download in other formats: