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

Context Navigation

source: utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/domzgr.F90 @ 13024

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

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

Download in other formats: