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 utils/tools/DOMAINcfg/src – NEMO

Context Navigation

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

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

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

Download in other formats: