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 @ 13056

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

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

Download in other formats: