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/2014/dev_r4627_COMODO_UPW2D/NEMOGCM/CONFIG/UPW2D_eo_NOadv_vvl_EXP/MY_SRC – NEMO

Context Navigation

source: branches/2014/dev_r4627_COMODO_UPW2D/NEMOGCM/CONFIG/UPW2D_eo_NOadv_vvl_EXP/MY_SRC/domzgr.F90 @ 4648

Last change on this file since 4648 was 4648, checked in by flavoni, 10 years ago
add new experience for cas test Upwelling, for WP item CNRS-7
File size: 94.5 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	!! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function
18	!! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case
19	!!----------------------------------------------------------------------
20
21	!!----------------------------------------------------------------------
22	!! dom_zgr : defined the ocean vertical coordinate system
23	!! zgr_bat : bathymetry fields (levels and meters)
24	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
25	!! zgr_bat_ctl : check the bathymetry files
26	!! zgr_bot_level: deepest ocean level for t-, u, and v-points
27	!! zgr_z : reference z-coordinate
28	!! zgr_zco : z-coordinate
29	!! zgr_zps : z-coordinate with partial steps
30	!! zgr_sco : s-coordinate
31	!! fssig : tanh stretch function
32	!! fssig1 : Song and Haidvogel 1994 stretch function
33	!! fgamma : Siddorn and Furner 2012 stretching function
34	!!---------------------------------------------------------------------
35	USE oce ! ocean variables
36	USE dom_oce ! ocean domain
37	USE closea ! closed seas
38	USE c1d ! 1D vertical configuration
39	USE in_out_manager ! I/O manager
40	USE iom ! I/O library
41	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
42	USE lib_mpp ! distributed memory computing library
43	USE wrk_nemo ! Memory allocation
44	USE timing ! Timing
45
46	IMPLICIT NONE
47	PRIVATE
48
49	PUBLIC dom_zgr ! called by dom_init.F90
50
51	! !!* Namelist namzgr_sco *
52	LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
53	LOGICAL :: ln_s_sf12 ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
54	!
55	REAL(wp) :: rn_sbot_min ! minimum depth of s-bottom surface (>0) (m)
56	REAL(wp) :: rn_sbot_max ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
57	REAL(wp) :: rn_rmax ! maximum cut-off r-value allowed (0<rn_rmax<1)
58	REAL(wp) :: rn_hc ! Critical depth for transition from sigma to stretched coordinates
59	! Song and Haidvogel 1994 stretching parameters
60	REAL(wp) :: rn_theta ! surface control parameter (0<=rn_theta<=20)
61	REAL(wp) :: rn_thetb ! bottom control parameter (0<=rn_thetb<= 1)
62	REAL(wp) :: rn_bb ! stretching parameter
63	! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
64	! Siddorn and Furner stretching parameters
65	LOGICAL :: ln_sigcrit ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
66	REAL(wp) :: rn_alpha ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
67	REAL(wp) :: rn_efold ! efold length scale for transition to stretched coord
68	REAL(wp) :: rn_zs ! depth of surface grid box
69	! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
70	REAL(wp) :: rn_zb_a ! bathymetry scaling factor for calculating Zb
71	REAL(wp) :: rn_zb_b ! offset for calculating Zb
72
73	!! * Substitutions
74	# include "domzgr_substitute.h90"
75	# include "vectopt_loop_substitute.h90"
76	!!----------------------------------------------------------------------
77	!! NEMO/OPA 3.3.1 , NEMO Consortium (2011)
78	!! $Id: domzgr.F90 4292 2013-11-20 16:28:04Z cetlod $
79	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
80	!!----------------------------------------------------------------------
81	CONTAINS
82
83	SUBROUTINE dom_zgr
84	!!----------------------------------------------------------------------
85	!! * ROUTINE dom_zgr *
86	!!
87	!! ** Purpose : set the depth of model levels and the resulting
88	!! vertical scale factors.
89	!!
90	!! ** Method : - reference 1D vertical coordinate (gdep._1d, e3._1d)
91	!! - read/set ocean depth and ocean levels (bathy, mbathy)
92	!! - vertical coordinate (gdep., e3.) depending on the
93	!! coordinate chosen :
94	!! ln_zco=T z-coordinate
95	!! ln_zps=T z-coordinate with partial steps
96	!! ln_zco=T s-coordinate
97	!!
98	!! ** Action : define gdep., e3., mbathy and bathy
99	!!----------------------------------------------------------------------
100	INTEGER :: ioptio, ibat ! local integer
101	INTEGER :: ios
102	!
103	NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco
104	!!----------------------------------------------------------------------
105	!
106	IF( nn_timing == 1 ) CALL timing_start('dom_zgr')
107	!
108	REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate
109	READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
110	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist', lwp )
111
112	REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate
113	READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
114	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist', lwp )
115	WRITE ( numond, namzgr )
116
117	IF(lwp) THEN ! Control print
118	WRITE(numout,*)
119	WRITE(numout,*) 'dom_zgr : vertical coordinate'
120	WRITE(numout,*) '~~~~~~~'
121	WRITE(numout,*) ' Namelist namzgr : set vertical coordinate'
122	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
123	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
124	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
125	ENDIF
126
127	ioptio = 0 ! Check Vertical coordinate options
128	IF( ln_zco ) ioptio = ioptio + 1
129	IF( ln_zps ) ioptio = ioptio + 1
130	IF( ln_sco ) ioptio = ioptio + 1
131	IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
132	!
133	! Build the vertical coordinate system
134	! ------------------------------------
135	CALL zgr_z ! Reference z-coordinate system (always called)
136	CALL zgr_bat ! Bathymetry fields (levels and meters)
137	IF( lk_c1d ) CALL lbc_lnk( bathy , 'T', 1._wp ) ! 1D config.: same bathy value over the 3x3 domain
138	IF( ln_zco ) CALL zgr_zco ! z-coordinate
139	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
140	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
141	!
142	! final adjustment of mbathy & check
143	! -----------------------------------
144	IF( lzoom ) CALL zgr_bat_zoom ! correct mbathy in case of zoom subdomain
145	IF( .NOT.lk_c1d ) CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points
146	CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
147	!
148	IF( lk_c1d ) THEN ! 1D config.: same mbathy value over the 3x3 domain
149	ibat = mbathy(2,2)
150	mbathy(:,:) = ibat
151	END IF
152	!
153	IF( nprint == 1 .AND. lwp ) THEN
154	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
155	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
156	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )
157	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), &
158	& ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), &
159	& ' uw', MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL( e3vw_0(:,:,:)), &
160	& ' w ', MINVAL( e3w_0(:,:,:) )
161
162	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
163	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )
164	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), &
165	& ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), &
166	& ' uw', MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)), &
167	& ' w ', MAXVAL( e3w_0(:,:,:) )
168	ENDIF
169	!
170	IF( nn_timing == 1 ) CALL timing_stop('dom_zgr')
171	!
172	END SUBROUTINE dom_zgr
173
174
175	SUBROUTINE zgr_z
176	!!----------------------------------------------------------------------
177	!! * ROUTINE zgr_z *
178	!!
179	!! ** Purpose : set the depth of model levels and the resulting
180	!! vertical scale factors.
181	!!
182	!! ** Method : z-coordinate system (use in all type of coordinate)
183	!! The depth of model levels is defined from an analytical
184	!! function the derivative of which gives the scale factors.
185	!! both depth and scale factors only depend on k (1d arrays).
186	!! w-level: gdepw_1d = gdep(k)
187	!! e3w_1d(k) = dk(gdep)(k) = e3(k)
188	!! t-level: gdept_1d = gdep(k+0.5)
189	!! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
190	!!
191	!! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
192	!! - e3t_1d , e3w_1d : scale factors at T- and W-levels (m)
193	!!
194	!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
195	!!----------------------------------------------------------------------
196	INTEGER :: jk ! dummy loop indices
197	REAL(wp) :: zt, zw ! temporary scalars
198	REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
199	REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
200	REAL(wp) :: zrefdep ! depth of the reference level (~10m)
201	REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
202	!!----------------------------------------------------------------------
203	!
204	IF( nn_timing == 1 ) CALL timing_start('zgr_z')
205	!
206	! Set variables from parameters
207	! ------------------------------
208	zkth = ppkth ; zacr = ppacr
209	zdzmin = ppdzmin ; zhmax = pphmax
210	zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters
211
212	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
213	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
214	IF( ppa1 == pp_to_be_computed .AND. &
215	& ppa0 == pp_to_be_computed .AND. &
216	& ppsur == pp_to_be_computed ) THEN
217	!
218	za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
219	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
220	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
221	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
222	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
223	ELSE
224	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
225	za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter
226	ENDIF
227
228	IF(lwp) THEN ! Parameter print
229	WRITE(numout,*)
230	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
231	WRITE(numout,*) ' ~~~~~~~'
232	IF( ppkth == 0._wp ) THEN
233	!SF added ocean at the comment for layers
234	WRITE(numout,*) ' Uniform grid with ',jpk-1,' ocean layers'
235	WRITE(numout,*) ' Total depth :', zhmax
236	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
237	ELSE
238	IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
239	WRITE(numout,*) ' zsur, za0, za1 computed from '
240	WRITE(numout,*) ' zdzmin = ', zdzmin
241	WRITE(numout,*) ' zhmax = ', zhmax
242	ENDIF
243	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
244	WRITE(numout,*) ' zsur = ', zsur
245	WRITE(numout,*) ' za0 = ', za0
246	WRITE(numout,*) ' za1 = ', za1
247	WRITE(numout,*) ' zkth = ', zkth
248	WRITE(numout,*) ' zacr = ', zacr
249	IF( ldbletanh ) THEN
250	WRITE(numout,*) ' (Double tanh za2 = ', za2
251	WRITE(numout,*) ' parameters) zkth2= ', zkth2
252	WRITE(numout,*) ' zacr2= ', zacr2
253	ENDIF
254	ENDIF
255	ENDIF
256
257
258	! Reference z-coordinate (depth - scale factor at T- and W-points)
259	! ======================
260	IF( ppkth == 0._wp ) THEN ! uniform vertical grid
261	za1 = zhmax / FLOAT(jpk-1)
262	DO jk = 1, jpk
263	zw = FLOAT( jk )
264	zt = FLOAT( jk ) + 0.5_wp
265	gdepw_1d(jk) = ( zw - 1 ) * za1
266	gdept_1d(jk) = ( zt - 1 ) * za1
267	e3w_1d (jk) = za1
268	e3t_1d (jk) = za1
269	END DO
270	ELSE ! Madec & Imbard 1996 function
271	IF( .NOT. ldbletanh ) THEN
272	DO jk = 1, jpk
273	zw = REAL( jk , wp )
274	zt = REAL( jk , wp ) + 0.5_wp
275	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
276	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
277	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
278	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
279	END DO
280	ELSE
281	DO jk = 1, jpk
282	zw = FLOAT( jk )
283	zt = FLOAT( jk ) + 0.5_wp
284	! Double tanh function
285	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) &
286	& + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) )
287	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) &
288	& + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) )
289	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) &
290	& + za2 * TANH( (zw-zkth2) / zacr2 )
291	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) &
292	& + za2 * TANH( (zt-zkth2) / zacr2 )
293	END DO
294	ENDIF
295	gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero
296	ENDIF
297
298	!!gm BUG in s-coordinate this does not work!
299	! deepest/shallowest W level Above/Below ~10m
300	zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
301	nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
302	nla10 = nlb10 - 1 ! deepest W level Above ~10m
303	!!gm end bug
304
305	IF(lwp) THEN ! control print
306	WRITE(numout,*)
307	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
308	WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
309	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
310	ENDIF
311	DO jk = 1, jpk ! control positivity
312	IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 ' )
313	IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
314	END DO
315	!
316	IF( nn_timing == 1 ) CALL timing_stop('zgr_z')
317	!
318	END SUBROUTINE zgr_z
319
320
321	SUBROUTINE zgr_bat
322	!!----------------------------------------------------------------------
323	!! * ROUTINE zgr_bat *
324	!!
325	!! ** Purpose : set bathymetry both in levels and meters
326	!!
327	!! ** Method : read or define mbathy and bathy arrays
328	!! * level bathymetry:
329	!! The ocean basin geometry is given by a two-dimensional array,
330	!! mbathy, which is defined as follow :
331	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
332	!! at t-point (ji,jj).
333	!! = 0 over the continental t-point.
334	!! The array mbathy is checked to verified its consistency with
335	!! model option. in particular:
336	!! mbathy must have at least 1 land grid-points (mbathy<=0)
337	!! along closed boundary.
338	!! mbathy must be cyclic IF jperio=1.
339	!! mbathy must be lower or equal to jpk-1.
340	!! isolated ocean grid points are suppressed from mbathy
341	!! since they are only connected to remaining
342	!! ocean through vertical diffusion.
343	!! ntopo=-1 : rectangular channel or bassin with a bump
344	!! ntopo= 0 : flat rectangular channel or basin
345	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
346	!! bathy is read in 'bathy_meter.nc' NetCDF file
347	!!
348	!! ** Action : - mbathy: level bathymetry (in level index)
349	!! - bathy : meter bathymetry (in meters)
350	!!----------------------------------------------------------------------
351	INTEGER :: ji, jj, jl, jk ! dummy loop indices
352	INTEGER :: inum ! temporary logical unit
353	INTEGER :: ii_bump, ij_bump, ih ! bump center position
354	INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices
355	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
356	REAL(wp) :: zi, zj, zh, zhmin ! local scalars
357	INTEGER , POINTER, DIMENSION(:,:) :: idta ! global domain integer data
358	REAL(wp), POINTER, DIMENSION(:,:) :: zdta ! global domain scalar data
359	!!----------------------------------------------------------------------
360	!
361	IF( nn_timing == 1 ) CALL timing_start('zgr_bat')
362	!
363	CALL wrk_alloc( jpidta, jpjdta, idta )
364	CALL wrk_alloc( jpidta, jpjdta, zdta )
365	!
366	IF(lwp) WRITE(numout,*)
367	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
368	IF(lwp) WRITE(numout,*) ' SF '
369	IF(lwp) WRITE(numout,*) ' value of ntopo: ', ntopo
370	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
371	! ! ================== !
372	IF( ntopo == 0 .OR. ntopo == -1 .OR. ntopo == -2) THEN ! defined by hand !
373	! ! ================== !
374	! ! global domain level and meter bathymetry (idta,zdta)
375	!
376	IF( ntopo == 0 ) THEN ! flat basin
377	IF(lwp) WRITE(numout,*)
378	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
379	IF( rn_bathy > 0.01 ) THEN
380	IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist'
381	zdta(:,:) = rn_bathy
382	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
383	idta(:,:) = jpkm1
384	ELSE ! z-coordinate (zco or zps): step-like topography
385	idta(:,:) = jpkm1
386	DO jk = 1, jpkm1
387	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
388	END DO
389	ENDIF
390	ELSE
391	IF(lwp) WRITE(numout,*) ' Depth = depthw(jpkm1)'
392	idta(:,:) = jpkm1 ! before last level
393	zdta(:,:) = gdepw_1d(jpk) ! last w-point depth
394	h_oce = gdepw_1d(jpk)
395	ENDIF
396	!SF ELSE ! bump centered in the basin
397	ELSEIF( ntopo == 1 ) THEN ! bump centered in the basin
398	IF(lwp) WRITE(numout,*)
399	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
400	ii_bump = jpidta / 2 ! i-index of the bump center
401	ij_bump = jpjdta / 2 ! j-index of the bump center
402	r_bump = 50000._wp ! bump radius (meters)
403	h_bump = 2700._wp ! bump height (meters)
404	h_oce = gdepw_1d(jpk) ! background ocean depth (meters)
405	IF(lwp) WRITE(numout,*) ' bump characteristics: '
406	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
407	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
408	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
409	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
410	!
411	DO jj = 1, jpjdta ! zdta :
412	DO ji = 1, jpidta
413	zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
414	zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
415	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
416	END DO
417	END DO
418	! ! idta :
419	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
420	idta(:,:) = jpkm1
421	ELSE ! z-coordinate (zco or zps): step-like topography
422	idta(:,:) = jpkm1
423	DO jk = 1, jpkm1
424	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
425	END DO
426	ENDIF
427	!ENDIF
428	ELSEIF( ntopo == -2 ) THEN ! slope
429	IF(lwp) WRITE(numout,*) ' sono in ntopo = -2'
430	IF(lwp) WRITE(numout,*)
431	IF(lwp) WRITE(numout,*) ' bathymetry field: slope'
432	h_oce = gdepw_1d(jpk) ! background ocean depth (meters)
433	IF(lwp) WRITE(numout,*) ' h_oce = ', h_oce , ' meters'
434	!
435	zdta(:,:) = 0.
436	DO jj = 1, jpjdta ! zdta :
437	DO ji = 1, jpidta
438	!SF old north-south zdta(ji,jj) = 1.e-3 * ( h_oce - REAL(ji, wp) + 1.5_wp) * ppe1_m
439	zdta(ji,jj) = 1.e-3 * ( h_oce - REAL(jj, wp) + 1.5_wp) * ppe1_m
440	END DO
441	END DO
442	!SF old north-south IF(lwp) WRITE(numout,*) ' zdta_1 = ', zdta(1,1) , ' meters'
443	!SF old north-south IF(lwp) WRITE(numout,*) ' zdta_201 = ', zdta(201,1) , ' meters'
444	!SF old north-south IF(lwp) WRITE(numout,*) ' zdta_202 = ', zdta(202,1) , ' meters'
445	IF(lwp) WRITE(numout,*) ' zdta_1 = ', zdta(1,1) , ' meters'
446	IF(lwp) WRITE(numout,*) ' zdta_201 = ', zdta(1,201) , ' meters'
447	IF(lwp) WRITE(numout,*) ' zdta_202 = ', zdta(1,202) , ' meters'
448	! ! idta :
449	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
450	idta(:,:) = jpkm1
451	ELSE ! z-coordinate (zco or zps): step-like topography
452	idta(:,:) = jpkm1
453	DO jk = 1, jpkm1
454	WHERE( gdept_1d(jk) < zdta(:,:) .AND. zdta(:,:) <= gdept_1d(jk+1) ) idta(:,:) = jk
455	END DO
456	ENDIF
457	!
458	ENDIF
459	! ! set GLOBAL boundary conditions
460	! ! Caution : idta on the global domain: use of jperio, not nperio
461	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
462	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1._wp
463	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp
464	ELSEIF( jperio == 2 ) THEN
465	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
466	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0._wp
467	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0._wp
468	idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0._wp
469	!SF add periodicity North-South , nperio = 7
470	!SF ELSEIF( jperio == 7 ) THEN
471	!SF IF(lwp) WRITE(numout,*) ' NORTH SOUTH periodicity '
472	!SF idta( 1 , : ) = 0
473	!SF idta( jpidta, : ) = 0
474	ELSE
475	ih = 0 ; zh = 0._wp
476	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
477	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
478	idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh
479	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
480	idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh
481	ENDIF
482
483	! ! local domain level and meter bathymetries (mbathy,bathy)
484	mbathy(:,:) = 0 ! set to zero extra halo points
485	bathy (:,:) = 0._wp ! (require for mpp case)
486	DO jj = 1, nlcj ! interior values
487	DO ji = 1, nlci
488	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
489	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
490	END DO
491	END DO
492	!
493	! ! ================ !
494	ELSEIF( ntopo == 1 ) THEN ! read in file ! (over the local domain)
495	! ! ================ !
496	!
497	IF( ln_zco ) THEN ! zco : read level bathymetry
498	CALL iom_open ( 'bathy_level.nc', inum )
499	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
500	CALL iom_close( inum )
501	mbathy(:,:) = INT( bathy(:,:) )
502	! ! =====================
503	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
504	! ! =====================
505	IF( nn_cla == 0 ) THEN
506	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
507	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
508	DO ji = mi0(ii0), mi1(ii1)
509	DO jj = mj0(ij0), mj1(ij1)
510	mbathy(ji,jj) = 15
511	END DO
512	END DO
513	IF(lwp) WRITE(numout,*)
514	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
515	!
516	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
517	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
518	DO ji = mi0(ii0), mi1(ii1)
519	DO jj = mj0(ij0), mj1(ij1)
520	mbathy(ji,jj) = 12
521	END DO
522	END DO
523	IF(lwp) WRITE(numout,*)
524	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
525	ENDIF
526	!
527	ENDIF
528	!
529	ENDIF
530	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
531	CALL iom_open ( 'bathy_meter.nc', inum )
532	CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
533	CALL iom_close( inum )
534	!
535	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
536	!
537	IF( nn_cla == 0 ) THEN
538	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
539	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
540	DO ji = mi0(ii0), mi1(ii1)
541	DO jj = mj0(ij0), mj1(ij1)
542	bathy(ji,jj) = 284._wp
543	END DO
544	END DO
545	IF(lwp) WRITE(numout,*)
546	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
547	!
548	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
549	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
550	DO ji = mi0(ii0), mi1(ii1)
551	DO jj = mj0(ij0), mj1(ij1)
552	bathy(ji,jj) = 137._wp
553	END DO
554	END DO
555	IF(lwp) WRITE(numout,*)
556	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
557	ENDIF
558	!
559	ENDIF
560	!
561	ENDIF
562	! ! =============== !
563	ELSE ! error !
564	! ! =============== !
565	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
566	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
567	ENDIF
568	!
569	IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==!
570	!
571	IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
572	IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
573	ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth
574	ENDIF
575	zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels
576	WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
577	ELSE WHERE ; bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
578	END WHERE
579	IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
580	ENDIF
581	!
582	CALL wrk_dealloc( jpidta, jpjdta, idta )
583	CALL wrk_dealloc( jpidta, jpjdta, zdta )
584	!
585	IF( nn_timing == 1 ) CALL timing_stop('zgr_bat')
586	!
587	END SUBROUTINE zgr_bat
588
589
590	SUBROUTINE zgr_bat_zoom
591	!!----------------------------------------------------------------------
592	!! * ROUTINE zgr_bat_zoom *
593	!!
594	!! ** Purpose : - Close zoom domain boundary if necessary
595	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
596	!!
597	!! ** Method :
598	!!
599	!! ** Action : - update mbathy: level bathymetry (in level index)
600	!!----------------------------------------------------------------------
601	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
602	!!----------------------------------------------------------------------
603	!
604	IF(lwp) WRITE(numout,*)
605	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
606	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
607	!
608	! Zoom domain
609	! ===========
610	!
611	! Forced closed boundary if required
612	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
613	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
614	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
615	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
616	!
617	! Configuration specific domain modifications
618	! (here, ORCA arctic configuration: suppress Med Sea)
619	IF( cp_cfg == "orca" .AND. cp_cfz == "arctic" ) THEN
620	SELECT CASE ( jp_cfg )
621	! ! =======================
622	CASE ( 2 ) ! ORCA_R2 configuration
623	! ! =======================
624	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
625	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
626	ij0 = 98 ; ij1 = 110
627	! ! =======================
628	CASE ( 05 ) ! ORCA_R05 configuration
629	! ! =======================
630	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
631	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
632	ij0 = 314 ; ij1 = 370
633	END SELECT
634	!
635	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
636	!
637	ENDIF
638	!
639	END SUBROUTINE zgr_bat_zoom
640
641
642	SUBROUTINE zgr_bat_ctl
643	!!----------------------------------------------------------------------
644	!! * ROUTINE zgr_bat_ctl *
645	!!
646	!! ** Purpose : check the bathymetry in levels
647	!!
648	!! ** Method : The array mbathy is checked to verified its consistency
649	!! with the model options. in particular:
650	!! mbathy must have at least 1 land grid-points (mbathy<=0)
651	!! along closed boundary.
652	!! mbathy must be cyclic IF jperio=1.
653	!! mbathy must be lower or equal to jpk-1.
654	!! isolated ocean grid points are suppressed from mbathy
655	!! since they are only connected to remaining
656	!! ocean through vertical diffusion.
657	!! C A U T I O N : mbathy will be modified during the initializa-
658	!! tion phase to become the number of non-zero w-levels of a water
659	!! column, with a minimum value of 1.
660	!!
661	!! ** Action : - update mbathy: level bathymetry (in level index)
662	!! - update bathy : meter bathymetry (in meters)
663	!!----------------------------------------------------------------------
664	!!
665	INTEGER :: ji, jj, jl ! dummy loop indices
666	INTEGER :: icompt, ibtest, ikmax ! temporary integers
667	REAL(wp), POINTER, DIMENSION(:,:) :: zbathy
668
669	!!----------------------------------------------------------------------
670	!
671	IF( nn_timing == 1 ) CALL timing_start('zgr_bat_ctl')
672	!
673	CALL wrk_alloc( jpi, jpj, zbathy )
674	!
675	IF(lwp) WRITE(numout,*)
676	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
677	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
678	! ! Suppress isolated ocean grid points
679	IF(lwp) WRITE(numout,*)
680	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
681	IF(lwp) WRITE(numout,*)' -----------------------------------'
682	icompt = 0
683	DO jl = 1, 2
684	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
685	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
686	mbathy(jpi,:) = mbathy( 2 ,:)
687	ENDIF
688	!SF old case mono periodicity NORTH/SOUTH
689	!SF ELSEIF( nperio == 7 ) THEN
690	!SF mbathy(:, 1 ) = mbathy(:,jpjm1) ! local domain is cyclic north-south
691	!SF mbathy(:,jpj) = mbathy(:, 2 )
692	!SF ENDIF
693	DO jj = 2, jpjm1
694	DO ji = 2, jpim1
695	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
696	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
697	IF( ibtest < mbathy(ji,jj) ) THEN
698	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
699	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
700	mbathy(ji,jj) = ibtest
701	icompt = icompt + 1
702	ENDIF
703	END DO
704	END DO
705	END DO
706	IF( lk_mpp ) CALL mpp_sum( icompt )
707	IF( icompt == 0 ) THEN
708	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
709	ELSE
710	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
711	ENDIF
712	IF( lk_mpp ) THEN
713	zbathy(:,:) = FLOAT( mbathy(:,:) )
714	CALL lbc_lnk( zbathy, 'T', 1._wp )
715	mbathy(:,:) = INT( zbathy(:,:) )
716	!SF
717	IF(lwp) WRITE(numout,*) ' caso mpp dopo lk_mpp lbc_lnk valore di mbathy(1,:) = ', mbathy(1,:)
718	IF(lwp) WRITE(numout,*) ' caso mpp dopo lk_mpp lbc_lnk valore di mbathy(jpim1,:) = ', mbathy(jpim1,:)
719	IF(lwp) WRITE(numout,*) ' ======================================== '
720	IF(lwp) WRITE(numout,*) ' caso mpp dopo lk_mpp lbc_lnk valore di mbathy(jpi,:) = ', mbathy(jpi,:)
721	IF(lwp) WRITE(numout,*) ' caso mpp dopo lk_mpp lbc_lnk valore di mbathy(2,:) = ', mbathy(2,:)
722	!SF
723	ENDIF
724	! ! East-west cyclic boundary conditions
725	!SF
726	IF(lwp) WRITE(numout,*) ' valore di nperio = ', nperio
727	!SF
728	IF( nperio == 0 ) THEN
729	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: nperio = ', nperio
730	IF( lk_mpp ) THEN
731	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
732	IF( jperio /= 1 ) mbathy(1,:) = 0
733	ENDIF
734	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
735	IF( jperio /= 1 ) mbathy(nlci,:) = 0
736	ENDIF
737	ELSE
738	IF( ln_zco .OR. ln_zps ) THEN
739	mbathy( 1 ,:) = 0
740	mbathy(jpi,:) = 0
741	ELSE
742	mbathy( 1 ,:) = jpkm1
743	mbathy(jpi,:) = jpkm1
744	ENDIF
745	ENDIF
746	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
747	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: nperio = ', nperio
748	mbathy( 1 ,:) = mbathy(jpim1,:)
749	mbathy(jpi,:) = mbathy( 2 ,:)
750	ELSEIF( nperio == 2 ) THEN
751	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: nperio = ', nperio
752	!SF case periodicity NORTH/SOUTH
753	!SF ELSEIF( nperio == 7 ) THEN
754	!SF IF(lwp) WRITE(numout,*) ' new case : nperio = ', nperio
755	ELSE
756	IF(lwp) WRITE(numout,*) ' e r r o r'
757	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
758	! STOP 'dom_mba'
759	ENDIF
760	! Boundary condition on mbathy
761	IF( .NOT.lk_mpp ) THEN
762	!!gm !!bug ??? think about it !
763	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
764	zbathy(:,:) = FLOAT( mbathy(:,:) )
765	CALL lbc_lnk( zbathy, 'T', 1._wp )
766	mbathy(:,:) = INT( zbathy(:,:) )
767	ENDIF
768	! Number of ocean level inferior or equal to jpkm1
769	ikmax = 0
770	DO jj = 1, jpj
771	DO ji = 1, jpi
772	ikmax = MAX( ikmax, mbathy(ji,jj) )
773	END DO
774	END DO
775	!!gm !!! test to do: ikmax = MAX( mbathy(:,:) ) ???
776	IF( ikmax > jpkm1 ) THEN
777	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
778	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
779	ELSE IF( ikmax < jpkm1 ) THEN
780	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
781	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
782	ENDIF
783
784	IF( lwp .AND. nprint == 1 ) THEN ! control print
785	WRITE(numout,*)
786	WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '
787	WRITE(numout,*) ' ------------------'
788	CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )
789	WRITE(numout,*)
790	ENDIF
791	!
792	CALL wrk_dealloc( jpi, jpj, zbathy )
793	!
794	IF( nn_timing == 1 ) CALL timing_stop('zgr_bat_ctl')
795	!
796	END SUBROUTINE zgr_bat_ctl
797
798
799	SUBROUTINE zgr_bot_level
800	!!----------------------------------------------------------------------
801	!! * ROUTINE zgr_bot_level *
802	!!
803	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
804	!!
805	!! ** Method : computes from mbathy with a minimum value of 1 over land
806	!!
807	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
808	!! ocean level at t-, u- & v-points
809	!! (min value = 1 over land)
810	!!----------------------------------------------------------------------
811	!!
812	INTEGER :: ji, jj ! dummy loop indices
813	REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
814	!!----------------------------------------------------------------------
815	!
816	IF( nn_timing == 1 ) CALL timing_start('zgr_bot_level')
817	!
818	CALL wrk_alloc( jpi, jpj, zmbk )
819	!
820	IF(lwp) WRITE(numout,*)
821	IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
822	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
823	!
824	mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
825
826	! ! bottom k-index of W-level = mbkt+1
827	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
828	DO ji = 1, jpim1
829	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
830	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
831	END DO
832	END DO
833	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
834	zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk(zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
835	zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk(zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
836	!
837	CALL wrk_dealloc( jpi, jpj, zmbk )
838	!
839	IF( nn_timing == 1 ) CALL timing_stop('zgr_bot_level')
840	!
841	END SUBROUTINE zgr_bot_level
842
843
844	SUBROUTINE zgr_zco
845	!!----------------------------------------------------------------------
846	!! * ROUTINE zgr_zco *
847	!!
848	!! ** Purpose : define the z-coordinate system
849	!!
850	!! ** Method : set 3D coord. arrays to reference 1D array
851	!!----------------------------------------------------------------------
852	INTEGER :: jk
853	!!----------------------------------------------------------------------
854	!
855	IF( nn_timing == 1 ) CALL timing_start('zgr_zco')
856	!
857	DO jk = 1, jpk
858	gdept_0 (:,:,jk) = gdept_1d(jk)
859	gdepw_0 (:,:,jk) = gdepw_1d(jk)
860	gdep3w_0(:,:,jk) = gdepw_1d(jk)
861	e3t_0 (:,:,jk) = e3t_1d (jk)
862	e3u_0 (:,:,jk) = e3t_1d (jk)
863	e3v_0 (:,:,jk) = e3t_1d (jk)
864	e3f_0 (:,:,jk) = e3t_1d (jk)
865	e3w_0 (:,:,jk) = e3w_1d (jk)
866	e3uw_0 (:,:,jk) = e3w_1d (jk)
867	e3vw_0 (:,:,jk) = e3w_1d (jk)
868	END DO
869	!
870	IF( nn_timing == 1 ) CALL timing_stop('zgr_zco')
871	!
872	END SUBROUTINE zgr_zco
873
874
875	SUBROUTINE zgr_zps
876	!!----------------------------------------------------------------------
877	!! * ROUTINE zgr_zps *
878	!!
879	!! ** Purpose : the depth and vertical scale factor in partial step
880	!! z-coordinate case
881	!!
882	!! ** Method : Partial steps : computes the 3D vertical scale factors
883	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
884	!! a partial step representation of bottom topography.
885	!!
886	!! The reference depth of model levels is defined from an analytical
887	!! function the derivative of which gives the reference vertical
888	!! scale factors.
889	!! From depth and scale factors reference, we compute there new value
890	!! with partial steps on 3d arrays ( i, j, k ).
891	!!
892	!! w-level: gdepw_0(i,j,k) = gdep(k)
893	!! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k)
894	!! t-level: gdept_0(i,j,k) = gdep(k+0.5)
895	!! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
896	!!
897	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
898	!! we find the mbathy index of the depth at each grid point.
899	!! This leads us to three cases:
900	!!
901	!! - bathy = 0 => mbathy = 0
902	!! - 1 < mbathy < jpkm1
903	!! - bathy > gdepw_0(jpk) => mbathy = jpkm1
904	!!
905	!! Then, for each case, we find the new depth at t- and w- levels
906	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
907	!! and f-points.
908	!!
909	!! This routine is given as an example, it must be modified
910	!! following the user s desiderata. nevertheless, the output as
911	!! well as the way to compute the model levels and scale factors
912	!! must be respected in order to insure second order accuracy
913	!! schemes.
914	!!
915	!! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
916	!! - - - - - - - gdept_0, gdepw_0 and e3. are positives
917	!!
918	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
919	!!----------------------------------------------------------------------
920	!!
921	INTEGER :: ji, jj, jk ! dummy loop indices
922	INTEGER :: ik, it ! temporary integers
923	LOGICAL :: ll_print ! Allow control print for debugging
924	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
925	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
926	REAL(wp) :: zmax ! Maximum depth
927	REAL(wp) :: zdiff ! temporary scalar
928	REAL(wp) :: zrefdep ! temporary scalar
929	REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt
930	!!---------------------------------------------------------------------
931	!
932	IF( nn_timing == 1 ) CALL timing_start('zgr_zps')
933	!
934	CALL wrk_alloc( jpi, jpj, jpk, zprt )
935	!
936	IF(lwp) WRITE(numout,*)
937	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
938	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
939	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
940
941	ll_print = .FALSE. ! Local variable for debugging
942
943	IF(lwp .AND. ll_print) THEN ! control print of the ocean depth
944	WRITE(numout,*)
945	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
946	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
947	ENDIF
948
949
950	! bathymetry in level (from bathy_meter)
951	! ===================
952	zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
953	bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
954	WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
955	ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
956	END WHERE
957
958	! Compute mbathy for ocean points (i.e. the number of ocean levels)
959	! find the number of ocean levels such that the last level thickness
960	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
961	! e3t_1d is the reference level thickness
962	DO jk = jpkm1, 1, -1
963	zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
964	WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
965	END DO
966
967	! Scale factors and depth at T- and W-points
968	DO jk = 1, jpk ! intitialization to the reference z-coordinate
969	gdept_0(:,:,jk) = gdept_1d(jk)
970	gdepw_0(:,:,jk) = gdepw_1d(jk)
971	e3t_0 (:,:,jk) = e3t_1d (jk)
972	e3w_0 (:,:,jk) = e3w_1d (jk)
973	END DO
974	!
975	DO jj = 1, jpj
976	DO ji = 1, jpi
977	ik = mbathy(ji,jj)
978	IF( ik > 0 ) THEN ! ocean point only
979	! max ocean level case
980	IF( ik == jpkm1 ) THEN
981	zdepwp = bathy(ji,jj)
982	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
983	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
984	e3t_0(ji,jj,ik ) = ze3tp
985	e3t_0(ji,jj,ik+1) = ze3tp
986	e3w_0(ji,jj,ik ) = ze3wp
987	e3w_0(ji,jj,ik+1) = ze3tp
988	gdepw_0(ji,jj,ik+1) = zdepwp
989	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
990	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
991	!
992	ELSE ! standard case
993	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
994	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
995	ENDIF
996	!gm Bug? check the gdepw_1d
997	! ... on ik
998	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
999	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
1000	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
1001	e3t_0(ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
1002	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
1003	e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
1004	& * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1005	! ... on ik+1
1006	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1007	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1008	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1009	ENDIF
1010	ENDIF
1011	END DO
1012	END DO
1013	!
1014	it = 0
1015	DO jj = 1, jpj
1016	DO ji = 1, jpi
1017	ik = mbathy(ji,jj)
1018	IF( ik > 0 ) THEN ! ocean point only
1019	e3tp (ji,jj) = e3t_0(ji,jj,ik)
1020	e3wp (ji,jj) = e3w_0(ji,jj,ik)
1021	! test
1022	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
1023	IF( zdiff <= 0._wp .AND. lwp ) THEN
1024	it = it + 1
1025	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1026	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1027	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1028	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
1029	ENDIF
1030	ENDIF
1031	END DO
1032	END DO
1033
1034	! Scale factors and depth at U-, V-, UW and VW-points
1035	DO jk = 1, jpk ! initialisation to z-scale factors
1036	e3u_0 (:,:,jk) = e3t_1d(jk)
1037	e3v_0 (:,:,jk) = e3t_1d(jk)
1038	e3uw_0(:,:,jk) = e3w_1d(jk)
1039	e3vw_0(:,:,jk) = e3w_1d(jk)
1040	END DO
1041	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
1042	DO jj = 1, jpjm1
1043	DO ji = 1, fs_jpim1 ! vector opt.
1044	e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1045	e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1046	e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1047	e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1048	END DO
1049	END DO
1050	END DO
1051	CALL lbc_lnk( e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk( e3uw_0, 'U', 1._wp ) ! lateral boundary conditions
1052	CALL lbc_lnk( e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk( e3vw_0, 'V', 1._wp )
1053	!
1054	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1055	WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)
1056	WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)
1057	WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)
1058	WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)
1059	END DO
1060
1061	! Scale factor at F-point
1062	DO jk = 1, jpk ! initialisation to z-scale factors
1063	e3f_0(:,:,jk) = e3t_1d(jk)
1064	END DO
1065	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
1066	DO jj = 1, jpjm1
1067	DO ji = 1, fs_jpim1 ! vector opt.
1068	e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1069	END DO
1070	END DO
1071	END DO
1072	CALL lbc_lnk( e3f_0, 'F', 1._wp ) ! Lateral boundary conditions
1073	!
1074	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1075	WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)
1076	END DO
1077	!!gm bug ? : must be a do loop with mj0,mj1
1078	!
1079	e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
1080	e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)
1081	e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)
1082	e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)
1083	e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)
1084
1085	! Control of the sign
1086	IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )
1087	IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )
1088	IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )
1089	IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )
1090
1091	! Compute gdep3w_0 (vertical sum of e3w)
1092	gdep3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
1093	DO jk = 2, jpk
1094	gdep3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)
1095	END DO
1096
1097	! ! ================= !
1098	IF(lwp .AND. ll_print) THEN ! Control print !
1099	! ! ================= !
1100	DO jj = 1,jpj
1101	DO ji = 1, jpi
1102	ik = MAX( mbathy(ji,jj), 1 )
1103	zprt(ji,jj,1) = e3t_0 (ji,jj,ik)
1104	zprt(ji,jj,2) = e3w_0 (ji,jj,ik)
1105	zprt(ji,jj,3) = e3u_0 (ji,jj,ik)
1106	zprt(ji,jj,4) = e3v_0 (ji,jj,ik)
1107	zprt(ji,jj,5) = e3f_0 (ji,jj,ik)
1108	zprt(ji,jj,6) = gdep3w_0(ji,jj,ik)
1109	END DO
1110	END DO
1111	WRITE(numout,*)
1112	WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1113	WRITE(numout,*)
1114	WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1115	WRITE(numout,*)
1116	WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1117	WRITE(numout,*)
1118	WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1119	WRITE(numout,*)
1120	WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1121	WRITE(numout,*)
1122	WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1123	ENDIF
1124	!
1125	CALL wrk_dealloc( jpi, jpj, jpk, zprt )
1126	!
1127	IF( nn_timing == 1 ) CALL timing_stop('zgr_zps')
1128	!
1129	END SUBROUTINE zgr_zps
1130
1131	SUBROUTINE zgr_sco
1132	!!----------------------------------------------------------------------
1133	!! * ROUTINE zgr_sco *
1134	!!
1135	!! ** Purpose : define the s-coordinate system
1136	!!
1137	!! ** Method : s-coordinate
1138	!! The depth of model levels is defined as the product of an
1139	!! analytical function by the local bathymetry, while the vertical
1140	!! scale factors are defined as the product of the first derivative
1141	!! of the analytical function by the bathymetry.
1142	!! (this solution save memory as depth and scale factors are not
1143	!! 3d fields)
1144	!! - Read bathymetry (in meters) at t-point and compute the
1145	!! bathymetry at u-, v-, and f-points.
1146	!! hbatu = mi( hbatt )
1147	!! hbatv = mj( hbatt )
1148	!! hbatf = mi( mj( hbatt ) )
1149	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1150	!! function and its derivative given as function.
1151	!! z_gsigt(k) = fssig (k )
1152	!! z_gsigw(k) = fssig (k-0.5)
1153	!! z_esigt(k) = fsdsig(k )
1154	!! z_esigw(k) = fsdsig(k-0.5)
1155	!! Three options for stretching are give, and they can be modified
1156	!! following the users requirements. Nevertheless, the output as
1157	!! well as the way to compute the model levels and scale factors
1158	!! must be respected in order to insure second order accuracy
1159	!! schemes.
1160	!!
1161	!! The three methods for stretching available are:
1162	!!
1163	!! s_sh94 (Song and Haidvogel 1994)
1164	!! a sinh/tanh function that allows sigma and stretched sigma
1165	!!
1166	!! s_sf12 (Siddorn and Furner 2012?)
1167	!! allows the maintenance of fixed surface and or
1168	!! bottom cell resolutions (cf. geopotential coordinates)
1169	!! within an analytically derived stretched S-coordinate framework.
1170	!!
1171	!! s_tanh (Madec et al 1996)
1172	!! a cosh/tanh function that gives stretched coordinates
1173	!!
1174	!!----------------------------------------------------------------------
1175	!
1176	INTEGER :: ji, jj, jk, jl ! dummy loop argument
1177	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
1178	INTEGER :: ios ! Local integer output status for namelist read
1179	REAL(wp) :: zrmax, ztaper ! temporary scalars
1180	REAL(wp) :: zrfact
1181	REAL(wp) :: h_oce ! temporary scalars
1182	!
1183	REAL(wp), SAVE :: r1_jpkm1
1184	REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1185	REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1186
1187	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1188	rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1189	!!----------------------------------------------------------------------
1190	!
1191	IF( nn_timing == 1 ) CALL timing_start('zgr_sco')
1192	!
1193	CALL wrk_alloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1194	!
1195	REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1196	READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1197	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1198
1199	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1200	READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1201	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1202	WRITE ( numond, namzgr_sco )
1203
1204	IF(lwp) THEN ! control print
1205	WRITE(numout,*)
1206	WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1207	WRITE(numout,*) '~~~~~~~~~~~'
1208	WRITE(numout,*) ' Namelist namzgr_sco'
1209	WRITE(numout,*) ' stretching coeffs '
1210	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
1211	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
1212	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
1213	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
1214	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
1215	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
1216	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
1217	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
1218	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
1219	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
1220	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
1221	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
1222	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
1223	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
1224	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
1225	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
1226	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
1227	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
1228	ENDIF
1229
1230	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1231	hifu(:,:) = rn_sbot_min
1232	hifv(:,:) = rn_sbot_min
1233	hiff(:,:) = rn_sbot_min
1234
1235	r1_jpkm1 = 1._wp / REAL( jpkm1, wp )
1236
1237	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) ) ! set the maximum ocean depth
1238
1239	DO jj = 1, jpj
1240	DO ji = 1, jpi
1241	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1242	END DO
1243	END DO
1244	! ! ==============================
1245	! ! hbatt field
1246	! ! ==============================
1247	!SF old periodicity north-south
1248	!SF hbatt( : ,:) = bathy( : ,:) ! Special case "Upwelling 2-dimensions"
1249	!SF hbatt( 1 ,:) = hbatt( 2 ,:) ! NO zero hbatt in the domain
1250	!SF hbatt(jpi,:) = hbatt(jpim1,:)
1251	hbatt( :,: ) = bathy( :,: ) ! Special case "Upwelling 2-dimensions"
1252	hbatt( :,1 ) = hbatt( :,2 ) ! NO zero hbatt in the domain
1253	hbatt( :,jpj) = hbatt( :,jpjm1)
1254
1255	IF(lwp) THEN ! Control print
1256	WRITE(numout,*)
1257	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1258	WRITE(numout,*)
1259	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )
1260	IF( nprint == 1 ) THEN
1261	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1262	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1263	ENDIF
1264	ENDIF
1265
1266	! ! ==============================
1267	! ! hbatu, hbatv, hbatf fields
1268	! ! ==============================
1269	IF(lwp) THEN
1270	WRITE(numout,*)
1271	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1272	ENDIF
1273
1274	!SF special case "Upwelling 2-dimensions"
1275	h_oce = gdepw_1d(jpk)
1276	WRITE(numout,*) 'h_oce value : ', h_oce
1277	!SF old periodicity north-south
1278	!SF DO ji = 1, jpi
1279	!SF hbatt(ji,:) = 1.e-3 * ( h_oce - REAL(ji,wp) + 1.5_wp ) * ppe1_m
1280	!SF hbatu(ji,:) = 1.e-3 * ( h_oce - REAL(ji,wp) + 1._wp ) * ppe1_m
1281	!SF hbatv(ji,:) = 1.e-3 * ( h_oce - REAL(ji,wp) + 1.5_wp ) * ppe1_m
1282	!SF hbatf(ji,:) = 1.e-3 * ( h_oce - REAL(ji,wp) + 1._wp ) * ppe1_m
1283	!SF END DO
1284	DO jj = 1, jpj
1285	hbatt(:,jj) = 1.e-3 * ( h_oce - REAL(jj,wp) + 1.5_wp ) * ppe2_m
1286	hbatu(:,jj) = 1.e-3 * ( h_oce - REAL(jj,wp) + 1._wp ) * ppe2_m
1287	hbatv(:,jj) = 1.e-3 * ( h_oce - REAL(jj,wp) + 1.5_wp ) * ppe2_m
1288	hbatf(:,jj) = 1.e-3 * ( h_oce - REAL(jj,wp) + 1._wp ) * ppe2_m
1289	END DO
1290	!
1291
1292	DO jk = 1, jpk ! uniform vertical grid spacing: e3 independent of jk
1293	e3t_0(:,:,jk) = hbatt(:,:) * r1_jpkm1
1294	e3u_0(:,:,jk) = hbatu(:,:) * r1_jpkm1
1295	e3v_0(:,:,jk) = hbatv(:,:) * r1_jpkm1
1296	e3f_0(:,:,jk) = hbatf(:,:) * r1_jpkm1
1297	!
1298	e3w_0(:,:,jk) = hbatt(:,:) * r1_jpkm1
1299	e3uw_0(:,:,jk) = hbatu(:,:) * r1_jpkm1
1300	e3vw_0(:,:,jk) = hbatv(:,:) * r1_jpkm1
1301	END DO
1302	! ! uniform vertical grid spacing: gdep increases linearly with jk
1303	DO jk = 1, jpk ! use jk-1 ==>> gdepw(jk=1) = 0. and gdepw(k=jpkm1) = hbatt
1304	gdepw_0(:,:,jk) = hbatt(:,:) * REAL( jk-1,wp ) * r1_jpkm1
1305	gdept_0(:,:,jk) = hbatt(:,:) * ( (REAL( jk-1,wp ) )+ 0.5 ) * r1_jpkm1
1306	gdep3w_0(:,:,jk) = gdepw_0(:,:,jk) ! here special case: linear variation gdepw=gdep3w
1307	END DO
1308	!
1309	fsdept(:,:,:) = gdept_0 (:,:,:)
1310	fsdepw(:,:,:) = gdepw_0 (:,:,:)
1311	fsde3w(:,:,:) = gdep3w_0(:,:,:)
1312	fse3t (:,:,:) = e3t_0 (:,:,:)
1313	fse3u (:,:,:) = e3u_0 (:,:,:)
1314	fse3v (:,:,:) = e3v_0 (:,:,:)
1315	fse3f (:,:,:) = e3f_0 (:,:,:)
1316	fse3w (:,:,:) = e3w_0 (:,:,:)
1317	fse3uw(:,:,:) = e3uw_0 (:,:,:)
1318	fse3vw(:,:,:) = e3vw_0 (:,:,:)
1319	!!gm
1320	nprint = 1 ! force the control print below
1321	!!gm
1322
1323	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
1324	& ' MAX ', MAXVAL( mbathy(:,:) )
1325
1326	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
1327	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1328	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
1329	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gdep3w_0(:,:,:) )
1330	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
1331	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
1332	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
1333	& ' w ', MINVAL( e3w_0 (:,:,:) )
1334
1335	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
1336	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gdep3w_0(:,:,:) )
1337	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
1338	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
1339	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
1340	& ' w ', MAXVAL( e3w_0 (:,:,:) )
1341	ENDIF
1342	! END DO
1343	IF(lwp) THEN ! selected vertical profiles
1344	WRITE(numout,*)
1345	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1346	WRITE(numout,*) ' ~~~~~~ --------------------'
1347	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
1348	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
1349	& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
1350	DO jj = mj0(20), mj1(20)
1351	DO ji = mi0(20), mi1(20)
1352	WRITE(numout,*)
1353	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1354	WRITE(numout,*) ' ~~~~~~ --------------------'
1355	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
1356	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
1357	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1358	END DO
1359	END DO
1360	DO jj = mj0(74), mj1(74)
1361	DO ji = mi0(100), mi1(100)
1362	WRITE(numout,*)
1363	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1364	WRITE(numout,*) ' ~~~~~~ --------------------'
1365	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
1366	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
1367	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1368	END DO
1369	END DO
1370	ENDIF
1371
1372	!================================================================================
1373	! check the coordinate makes sense
1374	!================================================================================
1375	DO ji = 1, jpi
1376	DO jj = 1, jpj
1377
1378	IF( hbatt(ji,jj) > 0._wp) THEN
1379	DO jk = 1, mbathy(ji,jj)
1380	! check coordinate is monotonically increasing
1381	IF (fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN
1382	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1383	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1384	WRITE(numout,*) 'e3w',fse3w(ji,jj,:)
1385	WRITE(numout,*) 'e3t',fse3t(ji,jj,:)
1386	CALL ctl_stop( ctmp1 )
1387	ENDIF
1388	! and check it has never gone negative
1389	IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN
1390	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1391	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1392	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
1393	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
1394	CALL ctl_stop( ctmp1 )
1395	ENDIF
1396	! and check it never exceeds the total depth
1397	IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN
1398	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
1399	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
1400	WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)
1401	CALL ctl_stop( ctmp1 )
1402	ENDIF
1403	END DO
1404
1405	DO jk = 1, mbathy(ji,jj)-1
1406	! and check it never exceeds the total depth
1407	IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN
1408	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
1409	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
1410	WRITE(numout,*) 'gdept',fsdept(ji,jj,:)
1411	CALL ctl_stop( ctmp1 )
1412	ENDIF
1413	END DO
1414
1415	ENDIF
1416
1417	END DO
1418	END DO
1419	!
1420	CALL wrk_dealloc( jpi, jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1421	!
1422	IF( nn_timing == 1 ) CALL timing_stop('zgr_sco')
1423	!
1424	END SUBROUTINE zgr_sco
1425
1426	!!======================================================================
1427	SUBROUTINE s_sh94()
1428
1429	!!----------------------------------------------------------------------
1430	!! * ROUTINE s_sh94 *
1431	!!
1432	!! ** Purpose : stretch the s-coordinate system
1433	!!
1434	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
1435	!! mixed S/sigma coordinate
1436	!!
1437	!! Reference : Song and Haidvogel 1994.
1438	!!----------------------------------------------------------------------
1439	!
1440	INTEGER :: ji, jj, jk ! dummy loop argument
1441	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
1442	!
1443	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
1444	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
1445
1446	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1447	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1448
1449	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
1450	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
1451	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
1452	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
1453
1454	DO ji = 1, jpi
1455	DO jj = 1, jpj
1456
1457	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
1458	DO jk = 1, jpk
1459	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1460	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
1461	END DO
1462	ELSE ! shallow water, uniform sigma
1463	DO jk = 1, jpk
1464	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
1465	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1466	END DO
1467	ENDIF
1468	!
1469	DO jk = 1, jpkm1
1470	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1471	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1472	END DO
1473	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
1474	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
1475	!
1476	! Coefficients for vertical depth as the sum of e3w scale factors
1477	z_gsi3w3(ji,jj,1) = 0.5_wp * z_esigw3(ji,jj,1)
1478	DO jk = 2, jpk
1479	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
1480	END DO
1481	!
1482	DO jk = 1, jpk
1483	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1484	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1485	gdept_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
1486	gdepw_0 (ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
1487	gdep3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsi3w3(ji,jj,jk)+rn_hczcoeft )
1488	END DO
1489	!
1490	END DO ! for all jj's
1491	END DO ! for all ji's
1492
1493	DO ji = 1, jpim1
1494	DO jj = 1, jpjm1
1495	DO jk = 1, jpk
1496	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1497	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1498	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1499	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1500	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) &
1501	& + hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) &
1502	& / ( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1503	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1504	& / ( hbatt(ji,jj)+hbatt(ji+1,jj) )
1505	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1506	& / ( hbatt(ji,jj)+hbatt(ji,jj+1) )
1507	!
1508	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1509	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1510	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1511	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1512	!
1513	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1514	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1515	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1516	END DO
1517	END DO
1518	END DO
1519
1520	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1521	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1522
1523	END SUBROUTINE s_sh94
1524
1525	SUBROUTINE s_sf12
1526
1527	!!----------------------------------------------------------------------
1528	!! * ROUTINE s_sf12 *
1529	!!
1530	!! ** Purpose : stretch the s-coordinate system
1531	!!
1532	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
1533	!! mixed S/sigma/Z coordinate
1534	!!
1535	!! This method allows the maintenance of fixed surface and or
1536	!! bottom cell resolutions (cf. geopotential coordinates)
1537	!! within an analytically derived stretched S-coordinate framework.
1538	!!
1539	!!
1540	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1541	!!----------------------------------------------------------------------
1542	!
1543	INTEGER :: ji, jj, jk ! dummy loop argument
1544	REAL(wp) :: zsmth ! smoothing around critical depth
1545	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
1546	!
1547	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3
1548	REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
1549
1550	!
1551	CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1552	CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1553
1554	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp
1555	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
1556	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
1557	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
1558
1559	DO ji = 1, jpi
1560	DO jj = 1, jpj
1561
1562	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1563
1564	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
1565	! could be changed by users but care must be taken to do so carefully
1566	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
1567
1568	zzs = rn_zs / hbatt(ji,jj)
1569
1570	IF (rn_efold /= 0.0_wp) THEN
1571	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1572	ELSE
1573	zsmth = 1.0_wp
1574	ENDIF
1575
1576	DO jk = 1, jpk
1577	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1578	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1579	ENDDO
1580	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
1581	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
1582
1583	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1584
1585	DO jk = 1, jpk
1586	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1587	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1588	END DO
1589
1590	ELSE ! shallow water, z coordinates
1591
1592	DO jk = 1, jpk
1593	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1594	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1595	END DO
1596
1597	ENDIF
1598
1599	DO jk = 1, jpkm1
1600	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1601	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1602	END DO
1603	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
1604	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
1605
1606	! Coefficients for vertical depth as the sum of e3w scale factors
1607	z_gsi3w3(ji,jj,1) = 0.5 * z_esigw3(ji,jj,1)
1608	DO jk = 2, jpk
1609	z_gsi3w3(ji,jj,jk) = z_gsi3w3(ji,jj,jk-1) + z_esigw3(ji,jj,jk)
1610	END DO
1611
1612	DO jk = 1, jpk
1613	gdept_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
1614	gdepw_0 (ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
1615	gdep3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)
1616	END DO
1617
1618	ENDDO ! for all jj's
1619	ENDDO ! for all ji's
1620
1621	DO ji=1,jpi-1
1622	DO jj=1,jpj-1
1623
1624	DO jk = 1, jpk
1625	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) / &
1626	( hbatt(ji,jj)+hbatt(ji+1,jj) )
1627	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) / &
1628	( hbatt(ji,jj)+hbatt(ji,jj+1) )
1629	z_esigtf3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) + &
1630	hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk)+hbatt(ji+1,jj+1)z_esigt3(ji+1,jj+1,jk) ) / &
1631	( hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1632	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) / &
1633	( hbatt(ji,jj)+hbatt(ji+1,jj) )
1634	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) / &
1635	( hbatt(ji,jj)+hbatt(ji,jj+1) )
1636
1637	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
1638	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
1639	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
1640	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
1641	!
1642	e3w_0(ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
1643	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
1644	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
1645	END DO
1646
1647	ENDDO
1648	ENDDO
1649	!
1650	CALL lbc_lnk(e3t_0 ,'T',1.) ; CALL lbc_lnk(e3u_0 ,'T',1.)
1651	CALL lbc_lnk(e3v_0 ,'T',1.) ; CALL lbc_lnk(e3f_0 ,'T',1.)
1652	CALL lbc_lnk(e3w_0 ,'T',1.)
1653	CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.)
1654	!
1655	! ! =============
1656
1657	CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 )
1658	CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1659
1660	END SUBROUTINE s_sf12
1661
1662	SUBROUTINE s_tanh()
1663
1664	!!----------------------------------------------------------------------
1665	!! * ROUTINE s_tanh*
1666	!!
1667	!! ** Purpose : stretch the s-coordinate system
1668	!!
1669	!! ** Method : s-coordinate stretch
1670	!!
1671	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1672	!!----------------------------------------------------------------------
1673
1674	INTEGER :: ji, jj, jk ! dummy loop argument
1675	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
1676
1677	REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w
1678	REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw
1679
1680	CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
1681	CALL wrk_alloc( jpk, z_esigt, z_esigw )
1682
1683	z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp
1684	z_esigt = 0._wp ; z_esigw = 0._wp
1685
1686	DO jk = 1, jpk
1687	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1688	z_gsigt(jk) = -fssig( REAL(jk,wp) )
1689	END DO
1690	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
1691	!
1692	! Coefficients for vertical scale factors at w-, t- levels
1693	!!gm bug : define it from analytical function, not like juste bellow....
1694	!!gm or betteroffer the 2 possibilities....
1695	DO jk = 1, jpkm1
1696	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
1697	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
1698	END DO
1699	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
1700	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
1701	!
1702	! Coefficients for vertical depth as the sum of e3w scale factors
1703	z_gsi3w(1) = 0.5_wp * z_esigw(1)
1704	DO jk = 2, jpk
1705	z_gsi3w(jk) = z_gsi3w(jk-1) + z_esigw(jk)
1706	END DO
1707	!!gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1708	DO jk = 1, jpk
1709	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1710	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1711	gdept_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
1712	gdepw_0 (:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
1713	gdep3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsi3w(jk) + hift(:,:)zcoeft )
1714	END DO
1715	!!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1716	DO jj = 1, jpj
1717	DO ji = 1, jpi
1718	DO jk = 1, jpk
1719	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1720	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1721	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1722	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1723	!
1724	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1725	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1726	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1727	END DO
1728	END DO
1729	END DO
1730
1731	CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w )
1732	CALL wrk_dealloc( jpk, z_esigt, z_esigw )
1733
1734	END SUBROUTINE s_tanh
1735
1736	FUNCTION fssig( pk ) RESULT( pf )
1737	!!----------------------------------------------------------------------
1738	!! * ROUTINE fssig *
1739	!!
1740	!! ** Purpose : provide the analytical function in s-coordinate
1741	!!
1742	!! ** Method : the function provide the non-dimensional position of
1743	!! T and W (i.e. between 0 and 1)
1744	!! T-points at integer values (between 1 and jpk)
1745	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1746	!!----------------------------------------------------------------------
1747	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
1748	REAL(wp) :: pf ! sigma value
1749	!!----------------------------------------------------------------------
1750	!
1751	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
1752	& - TANH( rn_thetb * rn_theta ) ) &
1753	& * ( COSH( rn_theta ) &
1754	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
1755	& / ( 2._wp * SINH( rn_theta ) )
1756	!
1757	END FUNCTION fssig
1758
1759
1760	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
1761	!!----------------------------------------------------------------------
1762	!! * ROUTINE fssig1 *
1763	!!
1764	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
1765	!!
1766	!! ** Method : the function provides the non-dimensional position of
1767	!! T and W (i.e. between 0 and 1)
1768	!! T-points at integer values (between 1 and jpk)
1769	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1770	!!----------------------------------------------------------------------
1771	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
1772	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
1773	REAL(wp) :: pf1 ! sigma value
1774	!!----------------------------------------------------------------------
1775	!
1776	IF ( rn_theta == 0 ) then ! uniform sigma
1777	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
1778	ELSE ! stretched sigma
1779	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
1780	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
1781	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
1782	ENDIF
1783	!
1784	END FUNCTION fssig1
1785
1786
1787	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
1788	!!----------------------------------------------------------------------
1789	!! * ROUTINE fgamma *
1790	!!
1791	!! ** Purpose : provide analytical function for the s-coordinate
1792	!!
1793	!! ** Method : the function provides the non-dimensional position of
1794	!! T and W (i.e. between 0 and 1)
1795	!! T-points at integer values (between 1 and jpk)
1796	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1797	!!
1798	!! This method allows the maintenance of fixed surface and or
1799	!! bottom cell resolutions (cf. geopotential coordinates)
1800	!! within an analytically derived stretched S-coordinate framework.
1801	!!
1802	!! Reference : Siddorn and Furner, in prep
1803	!!----------------------------------------------------------------------
1804	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
1805	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
1806	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
1807	REAL(wp), INTENT(in ) :: pzs ! surface box depth
1808	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
1809	REAL(wp) :: za1,za2,za3 ! local variables
1810	REAL(wp) :: zn1,zn2 ! local variables
1811	REAL(wp) :: za,zb,zx ! local variables
1812	integer :: jk
1813	!!----------------------------------------------------------------------
1814	!
1815
1816	zn1 = 1./(jpk-1.)
1817	zn2 = 1. - zn1
1818
1819	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
1820	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
1821	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
1822
1823	za = pzb - za3*(pzs-za1)-za2
1824	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
1825	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
1826	zx = 1.0_wp-za/2.0_wp-zb
1827
1828	DO jk = 1, jpk
1829	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
1830	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
1831	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
1832	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
1833	ENDDO
1834
1835	!
1836	END FUNCTION fgamma
1837
1838	!!======================================================================
1839	END MODULE domzgr

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

Download in other formats: