Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

domzgr.F90

Last change on this file was 15348, checked in by jchanut, 3 years ago
#2638, set bathymetry minimum on child grids prior boundary update from parent.
File size: 106.1 KB

Line
1	MODULE domzgr
2	!!==============================================================================
3	!! * MODULE domzgr *
4	!! Ocean domain : definition of the vertical coordinate system
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	!! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye
20	!! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying
21	!!----------------------------------------------------------------------
22
23	!!----------------------------------------------------------------------
24	!! dom_zgr : defined the ocean vertical coordinate system
25	!! zgr_bat : bathymetry fields (levels and meters)
26	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
27	!! zgr_bat_ctl : check the bathymetry files
28	!! zgr_bot_level: deepest ocean level for t-, u, and v-points
29	!! zgr_z : reference z-coordinate
30	!! zgr_zco : z-coordinate
31	!! zgr_zps : z-coordinate with partial steps
32	!! zgr_sco : s-coordinate
33	!! fssig : tanh stretch function
34	!! fssig1 : Song and Haidvogel 1994 stretch function
35	!! fgamma : Siddorn and Furner 2012 stretching function
36	!!---------------------------------------------------------------------
37	USE dom_oce ! ocean domain
38	USE depth_e3 ! depth <=> e3
39	!
40	USE in_out_manager ! I/O manager
41	USE iom ! I/O library
42	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
43	USE lib_mpp ! distributed memory computing library
44	USE lib_fortran
45	USE dombat
46	USE domisf
47	USE agrif_connect
48	USE agrif_domzgr
49
50	IMPLICIT NONE
51	PRIVATE
52
53	PUBLIC dom_zgr ! called by dom_init.F90
54	!
55	! !!* Namelist namzgr_sco *
56	LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
57	LOGICAL :: ln_s_sf12 ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
58	!
59	REAL(wp) :: rn_sbot_min ! minimum depth of s-bottom surface (>0) (m)
60	REAL(wp) :: rn_sbot_max ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
61	REAL(wp) :: rn_rmax ! maximum cut-off r-value allowed (0<rn_rmax<1)
62	REAL(wp) :: rn_hc ! Critical depth for transition from sigma to stretched coordinates
63
64	! Song and Haidvogel 1994 stretching parameters
65	REAL(wp) :: rn_theta ! surface control parameter (0<=rn_theta<=20)
66	REAL(wp) :: rn_thetb ! bottom control parameter (0<=rn_thetb<= 1)
67	REAL(wp) :: rn_bb ! stretching parameter
68	! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
69	! Siddorn and Furner stretching parameters
70	LOGICAL :: ln_sigcrit ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
71	REAL(wp) :: rn_alpha ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
72	REAL(wp) :: rn_efold ! efold length scale for transition to stretched coord
73	REAL(wp) :: rn_zs ! depth of surface grid box
74	! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
75	REAL(wp) :: rn_zb_a ! bathymetry scaling factor for calculating Zb
76	REAL(wp) :: rn_zb_b ! offset for calculating Zb
77
78	!! * Substitutions
79	# include "do_loop_substitute.h90"
80	!!----------------------------------------------------------------------
81	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
82	!! $Id: dommsk.F90 13305 2020-07-14 17:12:25Z acc $
83	!! Software governed by the CeCILL license (see ./LICENSE)
84	!!----------------------------------------------------------------------
85	CONTAINS
86
87	SUBROUTINE dom_zgr( k_top, k_bot )
88	!!----------------------------------------------------------------------
89	!! * ROUTINE dom_zgr *
90	!!
91	!! ** Purpose : set the depth of model levels and the resulting
92	!! vertical scale factors.
93	!!
94	!! ** Method : - reference 1D vertical coordinate (gdep._1d, e3._1d)
95	!! - read/set ocean depth and ocean levels (bathy, mbathy)
96	!! - vertical coordinate (gdep., e3.) depending on the
97	!! coordinate chosen :
98	!! ln_zco=T z-coordinate
99	!! ln_zps=T z-coordinate with partial steps
100	!! ln_zco=T s-coordinate
101	!!
102	!! ** Action : define gdep., e3., mbathy and bathy
103	!!----------------------------------------------------------------------
104	INTEGER, DIMENSION(:,:), INTENT(out) :: k_top, k_bot ! ocean first and last level indices
105	!
106	INTEGER :: ioptio, ibat ! local integer
107	INTEGER :: ios
108	!
109	INTEGER :: jk
110	REAL(wp) :: zrefdep ! depth of the reference level (~10m)
111
112
113	NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav
114	!!----------------------------------------------------------------------
115	!
116	!
117	! REWIND( numnam_ref ) ! Namelist namzgr in reference namelist : Vertical coordinate
118	READ ( numnam_ref, namzgr, IOSTAT = ios, ERR = 901 )
119	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in reference namelist')
120
121	! REWIND( numnam_cfg ) ! Namelist namzgr in configuration namelist : Vertical coordinate
122	READ ( numnam_cfg, namzgr, IOSTAT = ios, ERR = 902 )
123	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr in configuration namelist')
124	IF(lwm) WRITE ( numond, namzgr )
125
126	IF(ln_read_cfg) THEN
127	IF(lwp) WRITE(numout,*)
128	IF(lwp) WRITE(numout,*) ' ==>>> Read vertical mesh in ', TRIM( cn_domcfg ), ' file'
129	!
130	CALL zgr_read ( ln_zco , ln_zps , ln_sco, ln_isfcav, &
131	& gdept_1d, gdepw_1d, e3t_1d, e3w_1d , & ! 1D gridpoints depth
132	& gdept_0 , gdepw_0 , & ! gridpoints depth
133	& e3t_0 , e3u_0 , e3v_0 , e3f_0 , & ! vertical scale factors
134	& e3w_0 , e3uw_0 , e3vw_0 , & ! vertical scale factors
135	& k_top , k_bot ) ! 1st & last ocean level
136	!
137	!!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears
138	! ! Compute gde3w_0 (vertical sum of e3w)
139	! gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
140	! DO jk = 2, jpk
141	! gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
142	! END DO
143	!
144
145	! ! top/bottom ocean level indices for t-, u- and v-points (f-point also for top)
146	CALL zgr_top_bot( k_top, k_bot ) ! with a minimum value set to 1
147
148	! ! deepest/shallowest W level Above/Below ~10m
149	!!gm BUG in s-coordinate this does not work!
150	zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
151	nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
152	nla10 = nlb10 - 1 ! deepest W level Above ~10m
153	!!gm end bug
154
155	ENDIF
156
157	IF(lwp) THEN ! Control print
158	WRITE(numout,*)
159	WRITE(numout,*) 'dom_zgr : vertical coordinate'
160	WRITE(numout,*) '~~~~~~~'
161	WRITE(numout,*) ' Namelist namzgr : set vertical coordinate'
162	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
163	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
164	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
165	WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav
166	ENDIF
167
168	ioptio = 0 ! Check Vertical coordinate options
169	IF( ln_zco ) ioptio = ioptio + 1
170	IF( ln_zps ) ioptio = ioptio + 1
171	IF( ln_sco ) ioptio = ioptio + 1
172	IF( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
173	!
174	IF ( ln_isfcav ) CALL zgr_isf_nam
175	ioptio = 0
176	IF ( ln_zco .AND. ln_isfcav ) ioptio = ioptio + 1
177	IF ( ln_sco .AND. ln_isfcav ) ioptio = ioptio + 1
178	IF( ioptio > 0 ) CALL ctl_stop( ' Cavity not tested/compatible with full step (zco) and sigma (ln_sco) ' )
179
180	#if defined key_agrif
181	! IF ( (.NOT.Agrif_Root()).AND.((.NOT.ln_zps).OR.(.NOT.Agrif_parent(ln_zps)))) THEN
182	! CALL ctl_stop( 'STOP', 'AGRIF zooms require ln_zps=T for both Child and Parent')
183	! ENDIF
184	#endif
185
186	IF(.NOT.ln_read_cfg) THEN
187	!
188	! Build the vertical coordinate system
189	! ------------------------------------
190	CALL zgr_z ! Reference z-coordinate system (always called)
191	CALL zgr_bat ! Bathymetry fields (levels and meters)
192	IF( ln_zco ) CALL zgr_zco ! z-coordinate
193	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
194	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
195	CALL zgr_bat_ctl ! check bathymetry (mbathy) and suppress isolated ocean points
196	!
197	! final adjustment of mbathy & check
198	! -----------------------------------
199	CALL zgr_bot_level ! deepest ocean level for t-, u- and v-points
200	k_bot = mbkt
201	CALL zgr_top_level ! shallowest ocean level for T-, U-, V- points
202	k_top = mikt
203	WHERE( bathy(:,:) <= 0._wp ) k_top(:,:) = 0 ! set k_top to zero over land
204	ENDIF
205	!
206	IF( lwp ) THEN
207	WRITE(numout,*) ' MIN val k_top ', MINVAL( k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) )
208	WRITE(numout,*) ' MIN val k_bot ', MINVAL( k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) )
209	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
210	& ' w ', MINVAL( gdepw_0(:,:,:) )
211	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), &
212	& ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), &
213	& ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), &
214	& ' w ', MINVAL( e3w_0(:,:,:) )
215
216	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
217	& ' w ', MAXVAL( gdepw_0(:,:,:) )
218	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), &
219	& ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), &
220	& ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), &
221	& ' w ', MAXVAL( e3w_0(:,:,:) )
222	ENDIF
223
224	!
225	END SUBROUTINE dom_zgr
226
227	SUBROUTINE zgr_read( ld_zco , ld_zps , ld_sco , ld_isfcav, & ! type of vertical coordinate
228	& pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d , & ! 1D reference vertical coordinate
229	& pdept , pdepw , & ! 3D t & w-points depth
230	& pe3t , pe3u , pe3v , pe3f , & ! vertical scale factors
231	& pe3w , pe3uw , pe3vw , & ! - - -
232	& k_top , k_bot ) ! top & bottom ocean level
233	!!---------------------------------------------------------------------
234	!! * ROUTINE zgr_read *
235	!!
236	!! ** Purpose : Read the vertical information in the domain configuration file
237	!!
238	!!----------------------------------------------------------------------
239	LOGICAL , INTENT(out) :: ld_zco, ld_zps, ld_sco ! vertical coordinate flags
240	LOGICAL , INTENT(out) :: ld_isfcav ! under iceshelf cavity flag
241	REAL(wp), DIMENSION(:) , INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m]
242	REAL(wp), DIMENSION(:) , INTENT(out) :: pe3t_1d , pe3w_1d ! 1D vertical scale factors [m]
243	REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pdept, pdepw ! grid-point depth [m]
244	REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3t , pe3u , pe3v , pe3f ! vertical scale factors [m]
245	REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3w , pe3uw, pe3vw ! - - -
246	INTEGER , DIMENSION(:,:) , INTENT(out) :: k_top , k_bot ! first & last ocean level
247	!
248	INTEGER :: jk ! dummy loop index
249	INTEGER :: inum ! local logical unit
250	REAL(WP) :: z_zco, z_zps, z_sco, z_cav
251	REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
252	!!----------------------------------------------------------------------
253	!
254	IF(lwp) THEN
255	WRITE(numout,*)
256	WRITE(numout,*) ' zgr_read : read the vertical coordinates in ', TRIM( cn_domcfg ), ' file'
257	WRITE(numout,*) ' ~~~~~~~~'
258	ENDIF
259	!
260	CALL iom_open( cn_domcfg, inum )
261	!
262	! !* type of vertical coordinate
263	CALL iom_get( inum, 'ln_zco' , z_zco )
264	CALL iom_get( inum, 'ln_zps' , z_zps )
265	CALL iom_get( inum, 'ln_sco' , z_sco )
266	IF( z_zco == 0._wp ) THEN ; ld_zco = .false. ; ELSE ; ld_zco = .true. ; ENDIF
267	IF( z_zps == 0._wp ) THEN ; ld_zps = .false. ; ELSE ; ld_zps = .true. ; ENDIF
268	IF( z_sco == 0._wp ) THEN ; ld_sco = .false. ; ELSE ; ld_sco = .true. ; ENDIF
269	!
270	! !* ocean cavities under iceshelves
271	CALL iom_get( inum, 'ln_isfcav', z_cav )
272	IF( z_cav == 0._wp ) THEN ; ld_isfcav = .false. ; ELSE ; ld_isfcav = .true. ; ENDIF
273	!
274	! !* vertical scale factors
275	CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , pe3t_1d ) ! 1D reference coordinate
276	CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , pe3w_1d )
277	!
278	CALL iom_get( inum, jpdom_global, 'e3t_0' , pe3t , cd_type = 'T', psgn = 1._wp, kfill = jpfillcopy ) ! 3D coordinate
279	CALL iom_get( inum, jpdom_global, 'e3u_0' , pe3u , cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
280	CALL iom_get( inum, jpdom_global, 'e3v_0' , pe3v , cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
281	CALL iom_get( inum, jpdom_global, 'e3f_0' , pe3f , cd_type = 'F', psgn = 1._wp, kfill = jpfillcopy )
282	CALL iom_get( inum, jpdom_global, 'e3w_0' , pe3w , cd_type = 'W', psgn = 1._wp, kfill = jpfillcopy )
283	CALL iom_get( inum, jpdom_global, 'e3uw_0' , pe3uw, cd_type = 'U', psgn = 1._wp, kfill = jpfillcopy )
284	CALL iom_get( inum, jpdom_global, 'e3vw_0' , pe3vw, cd_type = 'V', psgn = 1._wp, kfill = jpfillcopy )
285	!
286	! !* depths
287	! !- old depth definition (obsolescent feature)
288	IF( iom_varid( inum, 'gdept_1d', ldstop = .FALSE. ) > 0 .AND. &
289	& iom_varid( inum, 'gdepw_1d', ldstop = .FALSE. ) > 0 .AND. &
290	& iom_varid( inum, 'gdept_0' , ldstop = .FALSE. ) > 0 .AND. &
291	& iom_varid( inum, 'gdepw_0' , ldstop = .FALSE. ) > 0 ) THEN
292	CALL ctl_warn( 'zgr_read : old definition of depths and scale factors used ', &
293	& ' depths at t- and w-points read in the domain configuration file')
294	CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d )
295	CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d )
296	CALL iom_get( inum, jpdom_global , 'gdept_0' , pdept, kfill = jpfillcopy )
297	CALL iom_get( inum, jpdom_global , 'gdepw_0' , pdepw, kfill = jpfillcopy )
298	!
299	ELSE !- depths computed from e3. scale factors
300	CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d ) ! 1D reference depth
301	CALL e3_to_depth( pe3t , pe3w , pdept , pdepw ) ! 3D depths
302	IF(lwp) THEN
303	WRITE(numout,*)
304	WRITE(numout,*) ' Reference 1D z-coordinate depth and scale factors:'
305	WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
306	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )
307	ENDIF
308	ENDIF
309	!
310	! !* ocean top and bottom level
311	CALL iom_get( inum, jpdom_global, 'top_level' , z2d ) ! 1st wet T-points (ISF)
312	k_top(:,:) = NINT( z2d(:,:) )
313	CALL iom_get( inum, jpdom_global, 'bottom_level' , z2d ) ! last wet T-points
314	k_bot(:,:) = NINT( z2d(:,:) )
315	!
316	! reference depth for negative bathy (wetting and drying only)
317	! IF( ll_wd ) CALL iom_get( inum, 'rn_wd_ref_depth' , ssh_ref )
318	!
319	CALL iom_close( inum )
320	!
321	END SUBROUTINE zgr_read
322
323	SUBROUTINE zgr_top_bot( k_top, k_bot )
324	!!----------------------------------------------------------------------
325	!! * ROUTINE zgr_top_bot *
326	!!
327	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
328	!!
329	!! ** Method : computes from k_top and k_bot with a minimum value of 1 over land
330	!!
331	!! ** Action : mikt, miku, mikv : vertical indices of the shallowest
332	!! ocean level at t-, u- & v-points
333	!! (min value = 1)
334	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
335	!! ocean level at t-, u- & v-points
336	!! (min value = 1 over land)
337	!!----------------------------------------------------------------------
338	INTEGER , DIMENSION(:,:), INTENT(in) :: k_top, k_bot ! top & bottom ocean level indices
339	!
340	INTEGER :: ji, jj ! dummy loop indices
341	REAL(wp), DIMENSION(jpi,jpj) :: zk ! workspace
342	!!----------------------------------------------------------------------
343	!
344	IF(lwp) WRITE(numout,*)
345	IF(lwp) WRITE(numout,*) ' zgr_top_bot : ocean top and bottom k-index of T-, U-, V- and W-levels '
346	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
347	!
348	mikt(:,:) = MAX( k_top(:,:) , 1 ) ! top ocean k-index of T-level (=1 over land)
349	!
350	mbkt(:,:) = MAX( k_bot(:,:) , 1 ) ! bottom ocean k-index of T-level (=1 over land)
351
352	! ! N.B. top k-index of W-level = mikt
353	! ! bottom k-index of W-level = mbkt+1
354	DO jj = 1, jpjm1
355	DO ji = 1, jpim1
356	miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) )
357	mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) )
358	mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) )
359	!
360	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
361	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
362	END DO
363	END DO
364	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
365	zk(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'U', 1. ) ; miku(:,:) = MAX( NINT( zk(:,:) ), 1 )
366	zk(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'V', 1. ) ; mikv(:,:) = MAX( NINT( zk(:,:) ), 1 )
367	zk(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'F', 1. ) ; mikf(:,:) = MAX( NINT( zk(:,:) ), 1 )
368	!
369	zk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'U', 1. ) ; mbku(:,:) = MAX( NINT( zk(:,:) ), 1 )
370	zk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk( 'domzgr', zk, 'V', 1. ) ; mbkv(:,:) = MAX( NINT( zk(:,:) ), 1 )
371	!
372	END SUBROUTINE zgr_top_bot
373
374	SUBROUTINE zgr_z
375	!!----------------------------------------------------------------------
376	!! * ROUTINE zgr_z *
377	!!
378	!! ** Purpose : set the depth of model levels and the resulting
379	!! vertical scale factors.
380	!!
381	!! ** Method : z-coordinate system (use in all type of coordinate)
382	!! The depth of model levels is defined from an analytical
383	!! function the derivative of which gives the scale factors.
384	!! both depth and scale factors only depend on k (1d arrays).
385	!! w-level: gdepw_1d = gdep(k)
386	!! e3w_1d(k) = dk(gdep)(k) = e3(k)
387	!! t-level: gdept_1d = gdep(k+0.5)
388	!! e3t_1d(k) = dk(gdep)(k+0.5) = e3(k+0.5)
389	!!
390	!! ** Action : - gdept_1d, gdepw_1d : depth of T- and W-point (m)
391	!! - e3t_1d , e3w_1d : scale factors at T- and W-levels (m)
392	!!
393	!! Reference : Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
394	!!----------------------------------------------------------------------
395	INTEGER :: jk ! dummy loop indices
396	REAL(wp) :: zt, zw ! temporary scalars
397	REAL(wp) :: zsur, za0, za1, zkth ! Values set from parameters in
398	REAL(wp) :: zacr, zdzmin, zhmax ! par_CONFIG_Rxx.h90
399	REAL(wp) :: zrefdep ! depth of the reference level (~10m)
400	REAL(wp) :: za2, zkth2, zacr2 ! Values for optional double tanh function set from parameters
401	!!----------------------------------------------------------------------
402	!
403	! Set variables from parameters
404	! ------------------------------
405	zkth = ppkth ; zacr = ppacr
406	zdzmin = ppdzmin ; zhmax = pphmax
407	zkth2 = ppkth2 ; zacr2 = ppacr2 ! optional (ldbletanh=T) double tanh parameters
408
409	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
410	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
411	IF( ppa1 == pp_to_be_computed .AND. &
412	& ppa0 == pp_to_be_computed .AND. &
413	& ppsur == pp_to_be_computed ) THEN
414	!
415	za1 = ( ppdzmin - pphmax / FLOAT(jpkm1) ) &
416	& / ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
417	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
418	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
419	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
420	ELSE
421	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
422	za2 = ppa2 ! optional (ldbletanh=T) double tanh parameter
423	ENDIF
424
425	IF(lwp) THEN ! Parameter print
426	WRITE(numout,*)
427	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
428	WRITE(numout,*) ' ~~~~~~~'
429	IF( ppkth == 0._wp ) THEN
430	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
431	WRITE(numout,*) ' Total depth :', zhmax
432	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
433	ELSE
434	IF( ppa1 == 0._wp .AND. ppa0 == 0._wp .AND. ppsur == 0._wp ) THEN
435	WRITE(numout,*) ' zsur, za0, za1 computed from '
436	WRITE(numout,*) ' zdzmin = ', zdzmin
437	WRITE(numout,*) ' zhmax = ', zhmax
438	ENDIF
439	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
440	WRITE(numout,*) ' zsur = ', zsur
441	WRITE(numout,*) ' za0 = ', za0
442	WRITE(numout,*) ' za1 = ', za1
443	WRITE(numout,*) ' zkth = ', zkth
444	WRITE(numout,*) ' zacr = ', zacr
445	IF( ldbletanh ) THEN
446	WRITE(numout,*) ' (Double tanh za2 = ', za2
447	WRITE(numout,*) ' parameters) zkth2= ', zkth2
448	WRITE(numout,*) ' zacr2= ', zacr2
449	ENDIF
450	ENDIF
451	ENDIF
452
453
454	! Reference z-coordinate (depth - scale factor at T- and W-points)
455	! ======================
456	IF( ppkth == 0._wp ) THEN ! uniform vertical grid
457
458
459
460	za1 = zhmax / FLOAT(jpk-1)
461
462	DO jk = 1, jpk
463	zw = FLOAT( jk )
464	zt = FLOAT( jk ) + 0.5_wp
465	gdepw_1d(jk) = ( zw - 1 ) * za1
466	gdept_1d(jk) = ( zt - 1 ) * za1
467	e3w_1d (jk) = za1
468	e3t_1d (jk) = za1
469	END DO
470	ELSE ! Madec & Imbard 1996 function
471	IF( .NOT. ldbletanh ) THEN
472	DO jk = 1, jpk
473	zw = REAL( jk , wp )
474	zt = REAL( jk , wp ) + 0.5_wp
475	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth) / zacr ) ) )
476	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth) / zacr ) ) )
477	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth) / zacr )
478	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth) / zacr )
479	END DO
480	ELSE
481	DO jk = 1, jpk
482	zw = FLOAT( jk )
483	zt = FLOAT( jk ) + 0.5_wp
484	! Double tanh function
485	gdepw_1d(jk) = ( zsur + za0 * zw + za1 * zacr * LOG ( COSH( (zw-zkth ) / zacr ) ) &
486	& + za2 * zacr2* LOG ( COSH( (zw-zkth2) / zacr2 ) ) )
487	gdept_1d(jk) = ( zsur + za0 * zt + za1 * zacr * LOG ( COSH( (zt-zkth ) / zacr ) ) &
488	& + za2 * zacr2* LOG ( COSH( (zt-zkth2) / zacr2 ) ) )
489	e3w_1d (jk) = za0 + za1 * TANH( (zw-zkth ) / zacr ) &
490	& + za2 * TANH( (zw-zkth2) / zacr2 )
491	e3t_1d (jk) = za0 + za1 * TANH( (zt-zkth ) / zacr ) &
492	& + za2 * TANH( (zt-zkth2) / zacr2 )
493	END DO
494	ENDIF
495	gdepw_1d(1) = 0._wp ! force first w-level to be exactly at zero
496	ENDIF
497
498	IF ( ln_isfcav .OR. ln_e3_dep ) THEN ! e3. = dk[gdep]
499	!
500	DO jk = 1, jpkm1
501	e3t_1d(jk) = gdepw_1d(jk+1)-gdepw_1d(jk)
502	END DO
503	e3t_1d(jpk) = e3t_1d(jpk-1) ! we don't care because this level is masked in NEMO
504
505	DO jk = 2, jpk
506	e3w_1d(jk) = gdept_1d(jk) - gdept_1d(jk-1)
507	END DO
508	e3w_1d(1 ) = 2._wp * (gdept_1d(1) - gdepw_1d(1))
509	END IF
510
511	!!gm BUG in s-coordinate this does not work!
512	! deepest/shallowest W level Above/Below ~10m
513	zrefdep = 10._wp - 0.1_wp * MINVAL( e3w_1d ) ! ref. depth with tolerance (10% of minimum layer thickness)
514	nlb10 = MINLOC( gdepw_1d, mask = gdepw_1d > zrefdep, dim = 1 ) ! shallowest W level Below ~10m
515	nla10 = nlb10 - 1 ! deepest W level Above ~10m
516	!!gm end bug
517
518	IF(lwp) THEN ! control print
519	WRITE(numout,*)
520	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
521	WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" )
522	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d(jk), gdepw_1d(jk), e3t_1d(jk), e3w_1d(jk), jk = 1, jpk )
523	ENDIF
524	DO jk = 1, jpk ! control positivity
525	IF( e3w_1d (jk) <= 0._wp .OR. e3t_1d (jk) <= 0._wp ) CALL ctl_stop( 'dom:zgr_z: e3w_1d or e3t_1d =< 0 ' )
526	IF( gdepw_1d(jk) < 0._wp .OR. gdept_1d(jk) < 0._wp ) CALL ctl_stop( 'dom:zgr_z: gdepw_1d or gdept_1d < 0 ' )
527	END DO
528	!
529	END SUBROUTINE zgr_z
530
531
532	SUBROUTINE zgr_bat
533	!!----------------------------------------------------------------------
534	!! * ROUTINE zgr_bat *
535	!!
536	!! ** Purpose : set bathymetry both in levels and meters
537	!!
538	!! ** Method : read or define mbathy and bathy arrays
539	!! * level bathymetry:
540	!! The ocean basin geometry is given by a two-dimensional array,
541	!! mbathy, which is defined as follow :
542	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
543	!! at t-point (ji,jj).
544	!! = 0 over the continental t-point.
545	!! The array mbathy is checked to verified its consistency with
546	!! model option. in particular:
547	!! mbathy must have at least 1 land grid-points (mbathy<=0)
548	!! along closed boundary.
549	!! mbathy must be cyclic IF jperio=1.
550	!! mbathy must be lower or equal to jpk-1.
551	!! isolated ocean grid points are suppressed from mbathy
552	!! since they are only connected to remaining
553	!! ocean through vertical diffusion.
554	!! ntopo=-1 : rectangular channel or bassin with a bump
555	!! ntopo= 0 : flat rectangular channel or basin
556	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
557	!! bathy is read in 'bathy_meter.nc' NetCDF file
558	!!
559	!! ** Action : - mbathy: level bathymetry (in level index)
560	!! - bathy : meter bathymetry (in meters)
561	!!----------------------------------------------------------------------
562	INTEGER :: ji, jj, jk ! dummy loop indices
563	INTEGER :: inum ! temporary logical unit
564	INTEGER :: ierror ! error flag
565	INTEGER :: ii_bump, ij_bump, ih ! bump center position
566	INTEGER :: ii0, ii1, ij0, ij1, ik ! local indices
567	REAL(wp) :: r_bump , h_bump , h_oce ! bump characteristics
568	REAL(wp) :: zi, zj, zh, zhmin ! local scalars
569	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrand, zbatv
570	!!----------------------------------------------------------------------
571	!
572	IF(lwp) WRITE(numout,*)
573	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
574	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
575	! ! ================== !
576	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
577	! ! ================== !
578	!
579	!
580	IF( ntopo == 0 ) THEN ! flat basin
581	IF(lwp) WRITE(numout,*)
582	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
583	IF( rn_bathy > 0.01 ) THEN
584	IF(lwp) WRITE(numout,*) ' Depth = rn_bathy read in namelist'
585	bathy(:,:) = rn_bathy
586	IF( ln_sco ) THEN ! s-coordinate (zsc ):
587	mbathy(:,:) = jpkm1
588	ELSE ! z-coordinate (zco or zps): step-like topography
589	mbathy(:,:) = jpkm1
590	DO jk = 1, jpkm1
591	WHERE( gdept_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= gdept_1d(jk+1) ) mbathy(:,:) = jk
592	END DO
593	ENDIF
594	ELSE
595	IF(lwp) WRITE(numout,*) ' Depth = depthw(jpkm1)'
596	mbathy(:,:) = jpkm1 ! before last level
597	bathy(:,:) = gdepw_1d(jpk) ! last w-point depth
598	h_oce = gdepw_1d(jpk)
599	ENDIF
600	ELSE ! flat basin with random noise
601	IF(lwp) WRITE(numout,*)
602	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with random noise'
603	ALLOCATE( zrand(jpiglo,jpjglo), STAT=ierror )
604	IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'zgr_bat: unable to allocate zrand array' )
605	h_oce = gdepw_1d(jpk) ! background ocean depth (meters)
606	!
607	CALL RANDOM_SEED()
608	CALL RANDOM_NUMBER(zrand)
609	DO_2D( 0, 0, 0, 0 )
610	bathy(ji,jj) = h_oce + 0.1_wp h_oce (zrand(mig(ji),mjg(jj))-1._wp)
611	END_2D
612	IF ( cp_cfg=='OVERFLOW' ) THEN
613	DO jj=1,jpj
614	bathy(:,jj) = +( 500. + 0.5 * 1500. * ( 1.0 + tanh( (glamt(:,3) - 40.) / 7. ) ))
615	END DO
616	ENDIF
617	IF ( cp_cfg=='DOME' ) THEN
618	ALLOCATE(zbatv(jpi,jpj))
619	zbatv(:,:) = MIN(3600._wp, MAX( 600._wp, 600._wp -gphiv(:,:)1.e30.01 ))
620	bathy(:,1) = 0._wp
621	DO jj =2,jpj
622	bathy(:,jj) = 0.5_wp*(zbatv(:,jj) + zbatv(:,jj-1))
623	END DO
624	CALL lbc_lnk( 'zgr_bat', bathy, 'T', 1._wp )
625	WHERE (gphit(:,:) >0._wp) bathy(:,:) = 0._wp
626	! Dig inlet:
627	WHERE ((gphit(:,:)>0._wp).AND.(glamt(:,:)>-50._wp).AND.(glamt(:,:)<50._wp)) bathy(:,:) = 600._wp
628	DEALLOCATE(zbatv)
629	ENDIF
630	! CALL lbc_lnk( 'zgr_bat', bathy, 'T', 1._wp )
631	!
632	DEALLOCATE(zrand)
633	!
634	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
635	mbathy(:,:) = jpkm1
636	ELSE ! z-coordinate (zco or zps): step-like topography
637	mbathy(:,:) = jpkm1
638	DO jk = 1, jpkm1
639	WHERE( gdept_1d(jk) < bathy(:,:) .AND. bathy(:,:) <= gdept_1d(jk+1) ) mbathy(:,:) = jk
640	END DO
641	ENDIF
642	ENDIF
643	!
644	risfdep(:,:)=0.e0
645	misfdep(:,:)=1
646	!
647	! ! ================ !
648	ELSEIF( ntopo == 1 .OR. ntopo ==2 ) THEN ! read in file ! (over the local domain)
649	! ! ================ !
650	!
651	IF( ln_zco ) THEN ! zco : read level bathymetry
652	CALL iom_open ( cn_topolvl, inum )
653	CALL iom_get ( inum, jpdom_auto, cn_bathlvl, bathy )
654	CALL iom_close( inum )
655	mbathy(:,:) = INT( bathy(:,:) )
656	! initialisation isf variables
657	risfdep(:,:)=0._wp ; misfdep(:,:)=1
658	! ! =====================
659	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
660	! ! =====================
661	!
662	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
663	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
664	DO ji = mi0(ii0), mi1(ii1)
665	DO jj = mj0(ij0), mj1(ij1)
666	mbathy(ji,jj) = 15
667	END DO
668	END DO
669	IF(lwp) WRITE(numout,*)
670	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
671	!
672	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
673	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
674	DO ji = mi0(ii0), mi1(ii1)
675	DO jj = mj0(ij0), mj1(ij1)
676	mbathy(ji,jj) = 12
677	END DO
678	END DO
679	IF(lwp) WRITE(numout,*)
680	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
681	!
682	ENDIF
683	!
684	ENDIF
685	IF( ln_zps .OR. ln_sco ) THEN ! zps or sco : read meter bathymetry
686	#if defined key_agrif
687	IF (agrif_root()) THEN
688	#endif
689	IF( ntopo == 1) THEN
690	CALL iom_open ( cn_topo, inum )
691	CALL iom_get ( inum, jpdom_auto, cn_bath, bathy )
692	CALL iom_close( inum )
693	ELSE
694	CALL dom_bat
695	ENDIF
696	#if defined key_agrif
697	ELSE
698	IF( ntopo == 1) THEN
699	CALL agrif_create_bathy_meter()
700	ELSE
701	CALL dom_bat
702	ENDIF
703	ENDIF
704	#endif
705	!
706	! initialisation isf variables
707	risfdep(:,:)=0._wp ; misfdep(:,:)=1
708	!
709	IF ( ln_isfcav ) THEN
710	CALL iom_open ( cn_fisfd, inum )
711	CALL iom_get ( inum, jpdom_auto, cn_visfd, risfdep )
712	CALL iom_close( inum )
713	END IF
714	!
715	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
716	!
717	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
718	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
719	DO ji = mi0(ii0), mi1(ii1)
720	DO jj = mj0(ij0), mj1(ij1)
721	bathy(ji,jj) = 284._wp
722	END DO
723	END DO
724	IF(lwp) WRITE(numout,*)
725	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
726	!
727	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
728	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
729	DO ji = mi0(ii0), mi1(ii1)
730	DO jj = mj0(ij0), mj1(ij1)
731	bathy(ji,jj) = 137._wp
732	END DO
733	END DO
734	IF(lwp) WRITE(numout,*)
735	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
736	!
737	ENDIF
738	!
739	ENDIF
740	! ! =============== !
741	ELSE ! error !
742	! ! =============== !
743	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
744	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
745	ENDIF
746	!
747	IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
748	IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
749	ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth
750	ENDIF
751	zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels
752	WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
753	ELSE WHERE ( risfdep == 0._wp ); bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
754	END WHERE
755	IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
756	ENDIF
757	#if defined key_agrif
758	IF ( .NOT.Agrif_Root() ) CALL agrif_bathymetry_connect
759	#endif
760	!
761	END SUBROUTINE zgr_bat
762
763	SUBROUTINE zgr_bat_ctl
764	!!----------------------------------------------------------------------
765	!! * ROUTINE zgr_bat_ctl *
766	!!
767	!! ** Purpose : check the bathymetry in levels
768	!!
769	!! ** Method : The array mbathy is checked to verified its consistency
770	!! with the model options. in particular:
771	!! mbathy must have at least 1 land grid-points (mbathy<=0)
772	!! along closed boundary.
773	!! mbathy must be cyclic IF jperio=1.
774	!! mbathy must be lower or equal to jpk-1.
775	!! isolated ocean grid points are suppressed from mbathy
776	!! since they are only connected to remaining
777	!! ocean through vertical diffusion.
778	!! C A U T I O N : mbathy will be modified during the initializa-
779	!! tion phase to become the number of non-zero w-levels of a water
780	!! column, with a minimum value of 1.
781	!!
782	!! ** Action : - update mbathy: level bathymetry (in level index)
783	!! - update bathy : meter bathymetry (in meters)
784	!!----------------------------------------------------------------------
785	INTEGER :: ji, jj, jl ! dummy loop indices
786	INTEGER :: icompt, ibtest, ikmax ! temporary integers
787	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zbathy
788	!!----------------------------------------------------------------------
789	!
790	ALLOCATE(zbathy(jpi,jpj))
791	!
792	IF(lwp) WRITE(numout,*)
793	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
794	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
795	! ! Suppress isolated ocean grid points
796	IF(lwp) WRITE(numout,*)
797	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
798	IF(lwp) WRITE(numout,*)' -----------------------------------'
799	icompt = 0
800	DO jl = 1, 2
801	IF( l_Iperio ) THEN
802	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
803	mbathy(jpi,:) = mbathy( 2 ,:)
804	ENDIF
805	zbathy(:,:) = FLOAT( mbathy(:,:) )
806	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
807	mbathy(:,:) = INT( zbathy(:,:) )
808
809	DO jj = 2, jpjm1
810	DO ji = 2, jpim1
811	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
812	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
813	IF( ibtest < mbathy(ji,jj) ) THEN
814	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
815	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
816	mbathy(ji,jj) = ibtest
817	icompt = icompt + 1
818	ENDIF
819	END DO
820	END DO
821	END DO
822
823	IF( lk_mpp ) CALL mpp_sum( 'domzgr', icompt )
824	IF( icompt == 0 ) THEN
825	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
826	ELSE
827	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
828	ENDIF
829
830	IF( lk_mpp ) THEN
831	zbathy(:,:) = FLOAT( mbathy(:,:) )
832	CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp )
833	mbathy(:,:) = INT( zbathy(:,:) )
834	ENDIF
835	! ! East-west cyclic boundary conditions
836
837	IF( jperio == 0 ) THEN
838	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: jperio = ', jperio
839	IF( ln_zco .OR. ln_zps ) THEN
840	mbathy( mi0( 1+nn_hls):mi1( 1+nn_hls),:) = 0
841	mbathy( mi0(jpiglo-nn_hls):mi1(jpiglo-nn_hls),:) = 0
842	ELSE
843	mbathy( mi0( 1+nn_hls):mi1( 1+nn_hls),:) = jpkm1
844	mbathy( mi0(jpiglo-nn_hls):mi1(jpiglo-nn_hls),:) = jpkm1
845	ENDIF
846	ELSEIF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
847	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio
848	! mbathy( 1 ,:) = mbathy(jpim1,:)
849	! mbathy(jpi,:) = mbathy( 2 ,:)
850	ELSEIF( jperio == 2 ) THEN
851	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: jperio = ', jperio
852	ELSE
853	IF(lwp) WRITE(numout,*) ' e r r o r'
854	IF(lwp) WRITE(numout,*) ' parameter , jperio = ', jperio
855	! STOP 'dom_mba'
856	ENDIF
857
858	! Boundary condition on mbathy
859	IF( .NOT.lk_mpp ) THEN
860	!!gm !!bug ??? think about it !
861	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
862	zbathy(:,:) = FLOAT( mbathy(:,:) )
863	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
864	mbathy(:,:) = INT( zbathy(:,:) )
865	ENDIF
866
867	! Number of ocean level inferior or equal to jpkm1
868	zbathy(:,:) = FLOAT( mbathy(:,:) )
869	ikmax = MAXVAL(zbathy(:,:))
870	CALL mpp_max( 'domzgr',ikmax)
871
872	IF( ikmax > jpkm1 ) THEN
873	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
874	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
875	ELSE IF( ikmax < jpkm1 ) THEN
876	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
877	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
878	ENDIF
879	!
880	DEALLOCATE( zbathy )
881	!
882	END SUBROUTINE zgr_bat_ctl
883
884
885	SUBROUTINE zgr_bot_level
886	!!----------------------------------------------------------------------
887	!! * ROUTINE zgr_bot_level *
888	!!
889	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
890	!!
891	!! ** Method : computes from mbathy with a minimum value of 1 over land
892	!!
893	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
894	!! ocean level at t-, u- & v-points
895	!! (min value = 1 over land)
896	!!----------------------------------------------------------------------
897	INTEGER :: ji, jj ! dummy loop indices
898	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmbk
899	!!----------------------------------------------------------------------
900	!
901	ALLOCATE( zmbk(jpi,jpj) )
902	!
903	IF(lwp) WRITE(numout,*)
904	IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
905	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
906	!
907	mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
908
909	! ! bottom k-index of W-level = mbkt+1
910	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
911	DO ji = 1, jpim1
912	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
913	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
914	END DO
915	END DO
916	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
917	zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
918	zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
919	!
920	DEALLOCATE( zmbk )
921	!
922	END SUBROUTINE zgr_bot_level
923
924
925	SUBROUTINE zgr_top_level
926	!!----------------------------------------------------------------------
927	!! * ROUTINE zgr_top_level *
928	!!
929	!! ** Purpose : defines the vertical index of ocean top (mik. arrays)
930	!!
931	!! ** Method : computes from misfdep with a minimum value of 1
932	!!
933	!! ** Action : mikt, miku, mikv : vertical indices of the shallowest
934	!! ocean level at t-, u- & v-points
935	!! (min value = 1)
936	!!----------------------------------------------------------------------
937	INTEGER :: ji, jj ! dummy loop indices
938	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmik
939	!!----------------------------------------------------------------------
940	!
941	ALLOCATE( zmik(jpi,jpj) )
942	!
943	IF(lwp) WRITE(numout,*)
944	IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
945	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
946	!
947	mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1)
948	! ! top k-index of W-level (=mikt)
949	DO jj = 1, jpjm1 ! top k-index of U- (U-) level
950	DO ji = 1, jpim1
951	miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) )
952	mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) )
953	mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) )
954	END DO
955	END DO
956
957	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
958	zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 )
959	zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 )
960	zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 )
961	!
962	DEALLOCATE( zmik )
963	!
964	END SUBROUTINE zgr_top_level
965
966
967	SUBROUTINE zgr_zco
968	!!----------------------------------------------------------------------
969	!! * ROUTINE zgr_zco *
970	!!
971	!! ** Purpose : define the reference z-coordinate system
972	!!
973	!! ** Method : set 3D coord. arrays to reference 1D array
974	!!----------------------------------------------------------------------
975	INTEGER :: jk
976	!!----------------------------------------------------------------------
977	!
978	DO jk = 1, jpk
979	gdept_0(:,:,jk) = gdept_1d(jk)
980	gdepw_0(:,:,jk) = gdepw_1d(jk)
981	e3t_0 (:,:,jk) = e3t_1d (jk)
982	e3u_0 (:,:,jk) = e3t_1d (jk)
983	e3v_0 (:,:,jk) = e3t_1d (jk)
984	e3f_0 (:,:,jk) = e3t_1d (jk)
985	e3w_0 (:,:,jk) = e3w_1d (jk)
986	e3uw_0 (:,:,jk) = e3w_1d (jk)
987	e3vw_0 (:,:,jk) = e3w_1d (jk)
988	END DO
989	!
990	END SUBROUTINE zgr_zco
991
992
993	SUBROUTINE zgr_zps
994	!!----------------------------------------------------------------------
995	!! * ROUTINE zgr_zps *
996	!!
997	!! ** Purpose : the depth and vertical scale factor in partial step
998	!! reference z-coordinate case
999	!!
1000	!! ** Method : Partial steps : computes the 3D vertical scale factors
1001	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
1002	!! a partial step representation of bottom topography.
1003	!!
1004	!! The reference depth of model levels is defined from an analytical
1005	!! function the derivative of which gives the reference vertical
1006	!! scale factors.
1007	!! From depth and scale factors reference, we compute there new value
1008	!! with partial steps on 3d arrays ( i, j, k ).
1009	!!
1010	!! w-level: gdepw_0(i,j,k) = gdep(k)
1011	!! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k)
1012	!! t-level: gdept_0(i,j,k) = gdep(k+0.5)
1013	!! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
1014	!!
1015	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
1016	!! we find the mbathy index of the depth at each grid point.
1017	!! This leads us to three cases:
1018	!!
1019	!! - bathy = 0 => mbathy = 0
1020	!! - 1 < mbathy < jpkm1
1021	!! - bathy > gdepw_0(jpk) => mbathy = jpkm1
1022	!!
1023	!! Then, for each case, we find the new depth at t- and w- levels
1024	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
1025	!! and f-points.
1026	!!
1027	!! This routine is given as an example, it must be modified
1028	!! following the user s desiderata. nevertheless, the output as
1029	!! well as the way to compute the model levels and scale factors
1030	!! must be respected in order to insure second order accuracy
1031	!! schemes.
1032	!!
1033	!! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
1034	!! - - - - - - - gdept_0, gdepw_0 and e3. are positives
1035	!!
1036	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
1037	!!----------------------------------------------------------------------
1038	INTEGER :: ji, jj, jk ! dummy loop indices
1039	INTEGER :: ik, it, ikb, ikt ! temporary integers
1040	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
1041	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
1042	REAL(wp) :: zdiff ! temporary scalar
1043	REAL(wp) :: zmax ! temporary scalar
1044	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zprt
1045	!!---------------------------------------------------------------------
1046	!
1047	ALLOCATE( zprt(jpi,jpj,jpk) )
1048	!
1049	IF(lwp) WRITE(numout,*)
1050	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
1051	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
1052	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
1053
1054	! compute position of the ice shelf grounding line
1055	! set bathy and isfdraft to 0 where grounded
1056	IF ( ln_isfcav ) CALL zgr_isf_zspace
1057
1058	! bathymetry in level (from bathy_meter)
1059	! ===================
1060	zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
1061	bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
1062	WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
1063	ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
1064	END WHERE
1065
1066	! Compute mbathy for ocean points (i.e. the number of ocean levels)
1067	! find the number of ocean levels such that the last level thickness
1068	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
1069	! e3t_1d is the reference level thickness
1070	DO jk = jpkm1, 1, -1
1071	zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1072	WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
1073	END DO
1074
1075	! Check compatibility between bathy and iceshelf draft
1076	! insure at least 2 wet level on the vertical under an ice shelf
1077	! compute misfdep and adjust isf draft if needed
1078	IF ( ln_isfcav ) CALL zgr_isf_kspace
1079
1080	! Scale factors and depth at T- and W-points
1081	DO jk = 1, jpk ! intitialization to the reference z-coordinate
1082	gdept_0(:,:,jk) = gdept_1d(jk)
1083	gdepw_0(:,:,jk) = gdepw_1d(jk)
1084	e3t_0 (:,:,jk) = e3t_1d (jk)
1085	e3w_0 (:,:,jk) = e3w_1d (jk)
1086	END DO
1087
1088	! Scale factors and depth at T- and W-points
1089	DO jj = 1, jpj
1090	DO ji = 1, jpi
1091	ik = mbathy(ji,jj)
1092	IF( ik > 0 ) THEN ! ocean point only
1093	! max ocean level case
1094	IF( ik == jpkm1 ) THEN
1095	zdepwp = bathy(ji,jj)
1096	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
1097	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1098	e3t_0(ji,jj,ik ) = ze3tp
1099	e3t_0(ji,jj,ik+1) = ze3tp
1100	e3w_0(ji,jj,ik ) = ze3wp
1101	e3w_0(ji,jj,ik+1) = ze3tp
1102	gdepw_0(ji,jj,ik+1) = zdepwp
1103	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
1104	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1105	!
1106	ELSE ! standard case
1107	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1108	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1109	ENDIF
1110	!gm Bug? check the gdepw_1d
1111	! ... on ik
1112	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
1113	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
1114	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
1115	e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
1116	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
1117	e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
1118	& * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1119	! ... on ik+1
1120	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1121	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1122	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1123	ENDIF
1124	ENDIF
1125	END DO
1126	END DO
1127	!
1128	it = 0
1129	DO jj = 1, jpj
1130	DO ji = 1, jpi
1131	ik = mbathy(ji,jj)
1132	IF( ik > 0 ) THEN ! ocean point only
1133	e3tp (ji,jj) = e3t_0(ji,jj,ik)
1134	e3wp (ji,jj) = e3w_0(ji,jj,ik)
1135	! test
1136	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
1137	IF( zdiff <= 0._wp .AND. lwp ) THEN
1138	it = it + 1
1139	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1140	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1141	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1142	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
1143	ENDIF
1144	ENDIF
1145	END DO
1146	END DO
1147	!
1148	! compute top scale factor if ice shelf
1149	IF (ln_isfcav) CALL zps_isf
1150	!
1151	! Scale factors and depth at U-, V-, UW and VW-points
1152	DO jk = 1, jpk ! initialisation to z-scale factors
1153	e3u_0 (:,:,jk) = e3t_1d(jk)
1154	e3v_0 (:,:,jk) = e3t_1d(jk)
1155	e3uw_0(:,:,jk) = e3w_1d(jk)
1156	e3vw_0(:,:,jk) = e3w_1d(jk)
1157	END DO
1158
1159	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
1160	DO jj = 1, jpjm1
1161	DO ji = 1, jpim1 ! vector opt.
1162	e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1163	e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1164	e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1165	e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1166	END DO
1167	END DO
1168	END DO
1169
1170	! update e3uw in case only 2 cells in the water column
1171	IF ( ln_isfcav ) CALL zps_isf_e3uv_w
1172	!
1173	CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp ) ! lateral boundary conditions
1174	CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1175	!
1176	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1177	WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)
1178	WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)
1179	WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)
1180	WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)
1181	END DO
1182
1183	! Scale factor at F-point
1184	DO jk = 1, jpk ! initialisation to z-scale factors
1185	e3f_0(:,:,jk) = e3t_1d(jk)
1186	END DO
1187	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
1188	DO jj = 1, jpjm1
1189	DO ji = 1, jpim1 ! vector opt.
1190	e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1191	END DO
1192	END DO
1193	END DO
1194	CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp ) ! Lateral boundary conditions
1195	!
1196	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1197	WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)
1198	END DO
1199	!!gm bug ? : must be a do loop with mj0,mj1
1200	!
1201	e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
1202	e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)
1203	e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)
1204	e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)
1205	e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)
1206
1207	! Control of the sign
1208	IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )
1209	IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )
1210	IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )
1211	IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )
1212	!
1213	! if in the future gde3w_0 need to be compute, use the function defined in NEMO
1214	! for now gde3w_0 computation is removed as not an output of domcfg
1215
1216	DEALLOCATE( zprt )
1217	!
1218	END SUBROUTINE zgr_zps
1219
1220
1221	SUBROUTINE zgr_sco
1222	!!----------------------------------------------------------------------
1223	!! * ROUTINE zgr_sco *
1224	!!
1225	!! ** Purpose : define the s-coordinate system
1226	!!
1227	!! ** Method : s-coordinate
1228	!! The depth of model levels is defined as the product of an
1229	!! analytical function by the local bathymetry, while the vertical
1230	!! scale factors are defined as the product of the first derivative
1231	!! of the analytical function by the bathymetry.
1232	!! (this solution save memory as depth and scale factors are not
1233	!! 3d fields)
1234	!! - Read bathymetry (in meters) at t-point and compute the
1235	!! bathymetry at u-, v-, and f-points.
1236	!! hbatu = mi( hbatt )
1237	!! hbatv = mj( hbatt )
1238	!! hbatf = mi( mj( hbatt ) )
1239	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1240	!! function and its derivative given as function.
1241	!! z_gsigt(k) = fssig (k )
1242	!! z_gsigw(k) = fssig (k-0.5)
1243	!! z_esigt(k) = fsdsig(k )
1244	!! z_esigw(k) = fsdsig(k-0.5)
1245	!! Three options for stretching are give, and they can be modified
1246	!! following the users requirements. Nevertheless, the output as
1247	!! well as the way to compute the model levels and scale factors
1248	!! must be respected in order to insure second order accuracy
1249	!! schemes.
1250	!!
1251	!! The three methods for stretching available are:
1252	!!
1253	!! s_sh94 (Song and Haidvogel 1994)
1254	!! a sinh/tanh function that allows sigma and stretched sigma
1255	!!
1256	!! s_sf12 (Siddorn and Furner 2012?)
1257	!! allows the maintenance of fixed surface and or
1258	!! bottom cell resolutions (cf. geopotential coordinates)
1259	!! within an analytically derived stretched S-coordinate framework.
1260	!!
1261	!! s_tanh (Madec et al 1996)
1262	!! a cosh/tanh function that gives stretched coordinates
1263	!!
1264	!!----------------------------------------------------------------------
1265	INTEGER :: ji, jj, jk, jl ! dummy loop argument
1266	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
1267	INTEGER :: ios ! Local integer output status for namelist read
1268	REAL(wp) :: zrmax, ztaper ! temporary scalars
1269	REAL(wp) :: zrfact
1270	!
1271	REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1272	REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1273	!!
1274	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1275	& rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1276	!!----------------------------------------------------------------------
1277	!
1278	ALLOCATE( zenv(jpi,jpj), ztmp(jpi,jpj), zmsk(jpi,jpj), zri(jpi,jpj), zrj(jpi,jpj), zhbat(jpi,jpj) , ztmpi1(jpi,jpj), ztmpi2(jpi,jpj), ztmpj1(jpi,jpj), ztmpj2(jpi,jpj) )
1279	!
1280	!REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1281	READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1282	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist')
1283
1284	!REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1285	READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1286	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist')
1287	IF(lwm) WRITE ( numond, namzgr_sco )
1288
1289	IF(lwp) THEN ! control print
1290	WRITE(numout,*)
1291	WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1292	WRITE(numout,*) '~~~~~~~~~~~'
1293	WRITE(numout,*) ' Namelist namzgr_sco'
1294	WRITE(numout,*) ' stretching coeffs '
1295	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
1296	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
1297	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
1298	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
1299	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
1300	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
1301	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
1302	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
1303	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
1304	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
1305	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
1306	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
1307	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
1308	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
1309	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
1310	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
1311	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
1312	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
1313	ENDIF
1314
1315	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1316	hifu(:,:) = rn_sbot_min
1317	hifv(:,:) = rn_sbot_min
1318	hiff(:,:) = rn_sbot_min
1319
1320	! ! set maximum ocean depth
1321	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1322
1323	DO jj = 1, jpj
1324	DO ji = 1, jpi
1325	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1326	END DO
1327	END DO
1328	! ! =============================
1329	! ! Define the envelop bathymetry (hbatt)
1330	! ! =============================
1331	! use r-value to create hybrid coordinates
1332	zenv(:,:) = bathy(:,:)
1333	!
1334	! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1335	DO jj = 1, jpj
1336	DO ji = 1, jpi
1337	IF( bathy(ji,jj) == 0._wp ) THEN
1338	iip1 = MIN( ji+1, jpi )
1339	ijp1 = MIN( jj+1, jpj )
1340	iim1 = MAX( ji-1, 1 )
1341	ijm1 = MAX( jj-1, 1 )
1342	!!gm BUG fix see ticket #1617
1343	IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) &
1344	& + bathy(iim1,jj ) + bathy(iip1,jj ) &
1345	& + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) &
1346	& zenv(ji,jj) = rn_sbot_min
1347	!!gm
1348	!!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + &
1349	!!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1350	!!gm zenv(ji,jj) = rn_sbot_min
1351	!!gm ENDIF
1352	!!gm end
1353	ENDIF
1354	END DO
1355	END DO
1356
1357	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1358	CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, kfillmode=jpfillnothing )
1359	!
1360	! smooth the bathymetry (if required)
1361	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
1362	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1363	!
1364	jl = 0
1365	zrmax = 1._wp
1366	!
1367	!
1368	! set scaling factor used in reducing vertical gradients
1369	zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1370	!
1371	! initialise temporary evelope depth arrays
1372	ztmpi1(:,:) = zenv(:,:)
1373	ztmpi2(:,:) = zenv(:,:)
1374	ztmpj1(:,:) = zenv(:,:)
1375	ztmpj2(:,:) = zenv(:,:)
1376	!
1377	! initialise temporary r-value arrays
1378	zri(:,:) = 1._wp
1379	zrj(:,:) = 1._wp
1380	! ! ================ !
1381	DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop !
1382	! ! ================ !
1383	jl = jl + 1
1384	zrmax = 0._wp
1385	! we set zrmax from previous r-values (zri and zrj) first
1386	! if set after current r-value calculation (as previously)
1387	! we could exit DO WHILE prematurely before checking r-value
1388	! of current zenv
1389	DO_2D( 0, 0, 0, 0 )
1390	zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1391	END_2D
1392	zri(:,:) = 0._wp
1393	zrj(:,:) = 0._wp
1394	DO_2D( 0, 0, 0, 0 )
1395	iip1 = MIN( ji+1, jpi ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1396	ijp1 = MIN( jj+1, jpj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1397	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1398	zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1399	END IF
1400	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1401	zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1402	END IF
1403	IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact
1404	IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact
1405	IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact
1406	IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact
1407	END_2D
1408	! IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1409	!
1410	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax
1411	!
1412	DO_2D( 0, 0, 0, 0 )
1413	zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1414	END_2D
1415	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1416	CALL lbc_lnk( 'toto',zenv, 'T', 1._wp, kfillmode=jpfillnothing)
1417	! ! ================ !
1418	END DO ! End loop !
1419	! ! ================ !
1420	DO jj = 1, jpj
1421	DO ji = 1, jpi
1422	zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1423	END DO
1424	END DO
1425	!
1426	! Envelope bathymetry saved in hbatt
1427	hbatt(:,:) = zenv(:,:)
1428	IF ((ntopo>0).AND.MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1429	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1430	DO jj = 1, jpj
1431	DO ji = 1, jpi
1432	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1433	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1434	END DO
1435	END DO
1436	ENDIF
1437	!
1438	! ! ==============================
1439	! ! hbatu, hbatv, hbatf fields
1440	! ! ==============================
1441	IF(lwp) THEN
1442	WRITE(numout,*)
1443	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1444	ENDIF
1445	hbatu(:,:) = rn_sbot_min
1446	hbatv(:,:) = rn_sbot_min
1447	hbatf(:,:) = rn_sbot_min
1448	DO jj = 1, jpjm1
1449	DO ji = 1, jpim1 ! NO vector opt.
1450	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
1451	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
1452	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
1453	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1454	END DO
1455	END DO
1456
1457	!
1458	! Apply lateral boundary condition
1459	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1460	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp )
1461	DO jj = 1, jpj
1462	DO ji = 1, jpi
1463	IF( hbatu(ji,jj) == 0._wp ) THEN
1464	!No worries about the following line when ln_wd == .true.
1465	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
1466	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
1467	ENDIF
1468	END DO
1469	END DO
1470	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp )
1471	DO jj = 1, jpj
1472	DO ji = 1, jpi
1473	IF( hbatv(ji,jj) == 0._wp ) THEN
1474	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
1475	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
1476	ENDIF
1477	END DO
1478	END DO
1479	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp )
1480	DO jj = 1, jpj
1481	DO ji = 1, jpi
1482	IF( hbatf(ji,jj) == 0._wp ) THEN
1483	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
1484	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
1485	ENDIF
1486	END DO
1487	END DO
1488
1489	!!bug: key_helsinki a verifer
1490	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1491	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1492	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1493	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1494
1495	IF( lwp ) THEN
1496	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
1497	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1498	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
1499	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1500	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1501	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1502	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1503	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1504	ENDIF
1505	!! helsinki
1506
1507	! ! =======================
1508	! ! s-ccordinate fields (gdep., e3.)
1509	! ! =======================
1510	!
1511	! non-dimensional "sigma" for model level depth at w- and t-levels
1512
1513
1514	!========================================================================
1515	! Song and Haidvogel 1994 (ln_s_sh94=T)
1516	! Siddorn and Furner 2012 (ln_sf12=T)
1517	! or tanh function (both false)
1518	!========================================================================
1519	IF ( ln_s_sh94 ) THEN
1520	CALL s_sh94()
1521	ELSE IF ( ln_s_sf12 ) THEN
1522	CALL s_sf12()
1523	ELSE
1524	CALL s_tanh()
1525	ENDIF
1526
1527	CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp )
1528	CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp )
1529	CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp )
1530	CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp )
1531	CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp )
1532	CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp )
1533	CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1534	!
1535	WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp
1536	WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp
1537	WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp
1538	WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp
1539	WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp
1540	WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp
1541	WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp
1542	!!
1543	! HYBRID :
1544	DO jj = 1, jpj
1545	DO ji = 1, jpi
1546	DO jk = 1, jpkm1
1547	IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1548	END DO
1549	END DO
1550	END DO
1551
1552	WHERE (bathy(:,:)<=0) mbathy(:,:) = 0
1553
1554	IF(lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
1555	& ' MAX ', MAXVAL( mbathy(:,:) )
1556
1557	IF( lwp ) THEN ! min max values over the local domain
1558	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1559	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
1560	& ' w ', MINVAL( gdepw_0(:,:,:) )
1561	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
1562	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
1563	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
1564	& ' w ', MINVAL( e3w_0 (:,:,:) )
1565
1566	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
1567	& ' w ', MAXVAL( gdepw_0(:,:,:) )
1568	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
1569	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
1570	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
1571	& ' w ', MAXVAL( e3w_0 (:,:,:) )
1572	ENDIF
1573	! END DO
1574	IF(lwp) THEN ! selected vertical profiles
1575	WRITE(numout,*)
1576	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1577	WRITE(numout,*) ' ~~~~~~ --------------------'
1578	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
1579	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
1580	& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
1581	DO jj = mj0(20), mj1(20)
1582	DO ji = mi0(20), mi1(20)
1583	WRITE(numout,*)
1584	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1585	WRITE(numout,*) ' ~~~~~~ --------------------'
1586	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
1587	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
1588	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1589	END DO
1590	END DO
1591	DO jj = mj0(74), mj1(74)
1592	DO ji = mi0(10), mi1(10)
1593	WRITE(numout,*)
1594	WRITE(numout,*) ' domzgr: vertical coordinates : point (10,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1595	WRITE(numout,*) ' ~~~~~~ --------------------'
1596	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
1597	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
1598	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1599	END DO
1600	END DO
1601	ENDIF
1602	!
1603	!================================================================================
1604	! check the coordinate makes sense
1605	!================================================================================
1606	DO ji = 1, jpi
1607	DO jj = 1, jpj
1608	!
1609	IF( hbatt(ji,jj) > 0._wp) THEN
1610	DO jk = 1, mbathy(ji,jj)
1611	! check coordinate is monotonically increasing
1612	IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
1613	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1614	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1615	WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
1616	WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
1617	CALL ctl_stop( ctmp1 )
1618	ENDIF
1619	! and check it has never gone negative
1620	IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
1621	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1622	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1623	WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1624	WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1625	CALL ctl_stop( ctmp1 )
1626	ENDIF
1627	! and check it never exceeds the total depth
1628	IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1629	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
1630	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
1631	WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1632	CALL ctl_stop( ctmp1 )
1633	ENDIF
1634	END DO
1635	!
1636	DO jk = 1, mbathy(ji,jj)-1
1637	! and check it never exceeds the total depth
1638	IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1639	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
1640	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
1641	WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1642	CALL ctl_stop( ctmp1 )
1643	ENDIF
1644	END DO
1645	ENDIF
1646	END DO
1647	END DO
1648	!
1649	DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1650	!
1651	END SUBROUTINE zgr_sco
1652
1653
1654	SUBROUTINE s_sh94()
1655	!!----------------------------------------------------------------------
1656	!! * ROUTINE s_sh94 *
1657	!!
1658	!! ** Purpose : stretch the s-coordinate system
1659	!!
1660	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
1661	!! mixed S/sigma coordinate
1662	!!
1663	!! Reference : Song and Haidvogel 1994.
1664	!!----------------------------------------------------------------------
1665	INTEGER :: ji, jj, jk ! dummy loop argument
1666	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
1667	REAL(wp) :: ztmpu, ztmpv, ztmpf
1668	REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
1669	!
1670	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1671	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
1672	!!----------------------------------------------------------------------
1673
1674	ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1675	ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) )
1676	ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1677
1678	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp
1679	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
1680	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
1681	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
1682	!
1683	DO ji = 1, jpi
1684	DO jj = 1, jpj
1685	!
1686	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
1687	DO jk = 1, jpk
1688	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1689	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
1690	END DO
1691	ELSE ! shallow water, uniform sigma
1692	DO jk = 1, jpk
1693	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
1694	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1695	END DO
1696	ENDIF
1697	!
1698	DO jk = 1, jpkm1
1699	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1700	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1701	END DO
1702	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
1703	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
1704	!
1705	DO jk = 1, jpk
1706	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1707	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1708	gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
1709	gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
1710	END DO
1711	!
1712	END DO ! for all jj's
1713	END DO ! for all ji's
1714
1715	DO ji = 1, jpim1
1716	DO jj = 1, jpjm1
1717	! extended for Wetting/Drying case
1718	ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
1719	ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
1720	ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1721	ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1722	ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1723	ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1724	& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1725	DO jk = 1, jpk
1726	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1727	& / ztmpu
1728	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1729	& / ztmpu
1730
1731	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1732	& / ztmpv
1733	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1734	& / ztmpv
1735
1736	z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
1737	& + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
1738	& + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
1739	& + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1740
1741	!
1742	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1743	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1744	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1745	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1746	!
1747	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1748	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1749	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1750	END DO
1751	END DO
1752	END DO
1753	!
1754	DEALLOCATE( z_gsigw3, z_gsigt3 )
1755	DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1756	!
1757	END SUBROUTINE s_sh94
1758
1759
1760	SUBROUTINE s_sf12
1761	!!----------------------------------------------------------------------
1762	!! * ROUTINE s_sf12 *
1763	!!
1764	!! ** Purpose : stretch the s-coordinate system
1765	!!
1766	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
1767	!! mixed S/sigma/Z coordinate
1768	!!
1769	!! This method allows the maintenance of fixed surface and or
1770	!! bottom cell resolutions (cf. geopotential coordinates)
1771	!! within an analytically derived stretched S-coordinate framework.
1772	!!
1773	!!
1774	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1775	!!----------------------------------------------------------------------
1776	INTEGER :: ji, jj, jk ! dummy loop argument
1777	REAL(wp) :: zsmth ! smoothing around critical depth
1778	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
1779	REAL(wp) :: ztmpu, ztmpv, ztmpf
1780	REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
1781	!
1782	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1783	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
1784	!!----------------------------------------------------------------------
1785	!
1786	ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1787	ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk))
1788	ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1789
1790	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp
1791	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
1792	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
1793	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
1794
1795	DO ji = 1, jpi
1796	DO jj = 1, jpj
1797
1798	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1799
1800	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
1801	! could be changed by users but care must be taken to do so carefully
1802	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
1803
1804	zzs = rn_zs / hbatt(ji,jj)
1805
1806	IF (rn_efold /= 0.0_wp) THEN
1807	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1808	ELSE
1809	zsmth = 1.0_wp
1810	ENDIF
1811
1812	DO jk = 1, jpk
1813	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1814	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1815	ENDDO
1816	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
1817	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
1818
1819	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1820
1821	DO jk = 1, jpk
1822	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1823	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1824	END DO
1825
1826	ELSE ! shallow water, z coordinates
1827
1828	DO jk = 1, jpk
1829	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1830	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1831	END DO
1832
1833	ENDIF
1834
1835	DO jk = 1, jpkm1
1836	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1837	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1838	END DO
1839	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
1840	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
1841
1842	DO jk = 1, jpk
1843	gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
1844	gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
1845	END DO
1846
1847	ENDDO ! for all jj's
1848	ENDDO ! for all ji's
1849
1850	DO ji=1,jpi-1
1851	DO jj=1,jpj-1
1852
1853	! extend to suit for Wetting/Drying case
1854	ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
1855	ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
1856	ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1857	ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1858	ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1859	ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1860	& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1861	DO jk = 1, jpk
1862	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1863	& / ztmpu
1864	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1865	& / ztmpu
1866	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1867	& / ztmpv
1868	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1869	& / ztmpv
1870
1871	z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
1872	& + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
1873	& + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
1874	& + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1875
1876	! Code prior to wetting and drying option (for reference)
1877	!z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1878	! /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1879	!
1880	!z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1881	! /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1882	!
1883	!z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1884	! /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1885	!
1886	!z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1887	! /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1888	!
1889	!z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
1890	! & +hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
1891	! +hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
1892	! & +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) &
1893	! /( hbatt(ji ,jj )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1894
1895	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
1896	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
1897	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
1898	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
1899	!
1900	e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
1901	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
1902	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
1903	END DO
1904
1905	ENDDO
1906	ENDDO
1907	!
1908	CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.)
1909	CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.)
1910	CALL lbc_lnk('domzgr',e3w_0 ,'T',1.)
1911	CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.)
1912	!
1913	DEALLOCATE( z_gsigw3, z_gsigt3 )
1914	DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1915	!
1916	END SUBROUTINE s_sf12
1917
1918
1919	SUBROUTINE s_tanh()
1920	!!----------------------------------------------------------------------
1921	!! * ROUTINE s_tanh*
1922	!!
1923	!! ** Purpose : stretch the s-coordinate system
1924	!!
1925	!! ** Method : s-coordinate stretch
1926	!!
1927	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1928	!!----------------------------------------------------------------------
1929	INTEGER :: ji, jj, jk ! dummy loop argument
1930	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
1931	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt
1932	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
1933	!!----------------------------------------------------------------------
1934
1935	ALLOCATE( z_gsigw(jpk), z_gsigt(jpk) )
1936	ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
1937
1938	z_gsigw = 0._wp ; z_gsigt = 0._wp
1939	z_esigt = 0._wp ; z_esigw = 0._wp
1940
1941	DO jk = 1, jpk
1942	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1943	z_gsigt(jk) = -fssig( REAL(jk,wp) )
1944	END DO
1945	IF( lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
1946	!
1947	! Coefficients for vertical scale factors at w-, t- levels
1948	!!gm bug : define it from analytical function, not like juste bellow....
1949	!!gm or betteroffer the 2 possibilities....
1950	DO jk = 1, jpkm1
1951	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
1952	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
1953	END DO
1954	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
1955	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
1956	!
1957	DO jk = 1, jpk
1958	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1959	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1960	gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
1961	gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
1962	END DO
1963
1964	DO jj = 1, jpj
1965	DO ji = 1, jpi
1966	DO jk = 1, jpk
1967	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1968	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1969	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1970	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
1971	!
1972	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
1973	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
1974	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
1975	END DO
1976	END DO
1977	END DO
1978	!
1979	DEALLOCATE( z_gsigw, z_gsigt )
1980	DEALLOCATE( z_esigt, z_esigw )
1981	!
1982	END SUBROUTINE s_tanh
1983
1984
1985	FUNCTION fssig( pk ) RESULT( pf )
1986	!!----------------------------------------------------------------------
1987	!! * ROUTINE fssig *
1988	!!
1989	!! ** Purpose : provide the analytical function in s-coordinate
1990	!!
1991	!! ** Method : the function provide the non-dimensional position of
1992	!! T and W (i.e. between 0 and 1)
1993	!! T-points at integer values (between 1 and jpk)
1994	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1995	!!----------------------------------------------------------------------
1996	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
1997	REAL(wp) :: pf ! sigma value
1998	!!----------------------------------------------------------------------
1999	!
2000	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
2001	& - TANH( rn_thetb * rn_theta ) ) &
2002	& * ( COSH( rn_theta ) &
2003	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
2004	& / ( 2._wp * SINH( rn_theta ) )
2005	!
2006	END FUNCTION fssig
2007
2008
2009	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2010	!!----------------------------------------------------------------------
2011	!! * ROUTINE fssig1 *
2012	!!
2013	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
2014	!!
2015	!! ** Method : the function provides the non-dimensional position of
2016	!! T and W (i.e. between 0 and 1)
2017	!! T-points at integer values (between 1 and jpk)
2018	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2019	!!----------------------------------------------------------------------
2020	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
2021	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
2022	REAL(wp) :: pf1 ! sigma value
2023	!!----------------------------------------------------------------------
2024	!
2025	IF ( rn_theta == 0 ) then ! uniform sigma
2026	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2027	ELSE ! stretched sigma
2028	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
2029	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
2030	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
2031	ENDIF
2032	!
2033	END FUNCTION fssig1
2034
2035
2036	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2037	!!----------------------------------------------------------------------
2038	!! * ROUTINE fgamma *
2039	!!
2040	!! ** Purpose : provide analytical function for the s-coordinate
2041	!!
2042	!! ** Method : the function provides the non-dimensional position of
2043	!! T and W (i.e. between 0 and 1)
2044	!! T-points at integer values (between 1 and jpk)
2045	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2046	!!
2047	!! This method allows the maintenance of fixed surface and or
2048	!! bottom cell resolutions (cf. geopotential coordinates)
2049	!! within an analytically derived stretched S-coordinate framework.
2050	!!
2051	!! Reference : Siddorn and Furner, in prep
2052	!!----------------------------------------------------------------------
2053	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
2054	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
2055	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
2056	REAL(wp), INTENT(in ) :: pzs ! surface box depth
2057	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
2058	!
2059	INTEGER :: jk ! dummy loop index
2060	REAL(wp) :: za1,za2,za3 ! local scalar
2061	REAL(wp) :: zn1,zn2 ! - -
2062	REAL(wp) :: za,zb,zx ! - -
2063	!!----------------------------------------------------------------------
2064	!
2065	zn1 = 1._wp / REAL( jpkm1, wp )
2066	zn2 = 1._wp - zn1
2067	!
2068	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
2069	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
2070	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
2071	!
2072	za = pzb - za3*(pzs-za1)-za2
2073	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
2074	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
2075	zx = 1.0_wp-za/2.0_wp-zb
2076	!
2077	DO jk = 1, jpk
2078	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
2079	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
2080	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
2081	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
2082	END DO
2083	!
2084	END FUNCTION fgamma
2085
2086	!!======================================================================
2087	END MODULE domzgr

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: utils/tools/DOMAINcfg/src/domzgr.F90

Download in other formats: