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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

domzgr.F90 in utils/tools/DOMAINcfg/src – NEMO

Context Navigation

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

Last change on this file since 14449 was 13390, checked in by mathiot, 4 years ago
ticket #2502: merge ticket branch into trunk. DOMAIN_cfg namelist contains now fields to specify input files names (bathy meter and level files, coord file, isf draft meter and level files), save it into the netcdf (dom_doc.exe) and re-generate the namelist if needed (xtrac_namelist.bash). The usage is documented in the DOMAIN_cfg README.rst.
File size: 108.7 KB

Line
1	MODULE domzgr
2	!!==============================================================================
3	!! * MODULE domzgr *
4	!! Ocean domain : definition of the vertical coordinate system
5	!!==============================================================================
6	!! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate
7	!! ! 1997-07 (G. Madec) lbc_lnk call
8	!! ! 1997-04 (J.-O. Beismann)
9	!! 8.5 ! 2002-09 (A. Bozec, G. Madec) F90: Free form and module
10	!! - ! 2002-09 (A. de Miranda) rigid-lid + islands
11	!! NEMO 1.0 ! 2003-08 (G. Madec) F90: Free form and module
12	!! - ! 2005-10 (A. Beckmann) modifications for hybrid s-ccordinates & new stretching function
13	!! 2.0 ! 2006-04 (R. Benshila, G. Madec) add zgr_zco
14	!! 3.0 ! 2008-06 (G. Madec) insertion of domzgr_zps.h90 & conding style
15	!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
16	!! 3.3 ! 2010-11 (G. Madec) add mbk. arrays associated to the deepest ocean level
17	!! 3.4 ! 2012-08 (J. Siddorn) added Siddorn and Furner stretching function
18	!! 3.4 ! 2012-12 (R. Bourdalle-Badie and G. Reffray) modify C1D case
19	!! 3.6 ! 2014-11 (P. Mathiot and C. Harris) add ice shelf capabilitye
20	!! 3.? ! 2015-11 (H. Liu) Modifications for Wetting/Drying
21	!!----------------------------------------------------------------------
22
23	!!----------------------------------------------------------------------
24	!! dom_zgr : defined the ocean vertical coordinate system
25	!! zgr_bat : bathymetry fields (levels and meters)
26	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
27	!! zgr_bat_ctl : check the bathymetry files
28	!! zgr_bot_level: deepest ocean level for t-, u, and v-points
29	!! zgr_z : reference z-coordinate
30	!! zgr_zco : z-coordinate
31	!! zgr_zps : z-coordinate with partial steps
32	!! zgr_sco : s-coordinate
33	!! fssig : tanh stretch function
34	!! fssig1 : Song and Haidvogel 1994 stretch function
35	!! fgamma : Siddorn and Furner 2012 stretching function
36	!!---------------------------------------------------------------------
37	USE dom_oce ! ocean domain
38	USE depth_e3 ! depth <=> e3
39	!
40	USE in_out_manager ! I/O manager
41	USE iom ! I/O library
42	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
43	USE lib_mpp ! distributed memory computing library
44	USE lib_fortran
45	USE dombat
46	USE domisf
47	USE agrif_domzgr
48
49	IMPLICIT NONE
50	PRIVATE
51
52	PUBLIC dom_zgr ! called by dom_init.F90
53	!
54	! !!* Namelist namzgr_sco *
55	LOGICAL :: ln_s_sh94 ! use hybrid s-sig Song and Haidvogel 1994 stretching function fssig1 (ln_sco=T)
56	LOGICAL :: ln_s_sf12 ! use hybrid s-z-sig Siddorn and Furner 2012 stretching function fgamma (ln_sco=T)
57	!
58	REAL(wp) :: rn_sbot_min ! minimum depth of s-bottom surface (>0) (m)
59	REAL(wp) :: rn_sbot_max ! maximum depth of s-bottom surface (= ocean depth) (>0) (m)
60	REAL(wp) :: rn_rmax ! maximum cut-off r-value allowed (0<rn_rmax<1)
61	REAL(wp) :: rn_hc ! Critical depth for transition from sigma to stretched coordinates
62
63	! Song and Haidvogel 1994 stretching parameters
64	REAL(wp) :: rn_theta ! surface control parameter (0<=rn_theta<=20)
65	REAL(wp) :: rn_thetb ! bottom control parameter (0<=rn_thetb<= 1)
66	REAL(wp) :: rn_bb ! stretching parameter
67	! ! ( rn_bb=0; top only, rn_bb =1; top and bottom)
68	! Siddorn and Furner stretching parameters
69	LOGICAL :: ln_sigcrit ! use sigma coordinates below critical depth (T) or Z coordinates (F) for Siddorn & Furner stretch
70	REAL(wp) :: rn_alpha ! control parameter ( > 1 stretch towards surface, < 1 towards seabed)
71	REAL(wp) :: rn_efold ! efold length scale for transition to stretched coord
72	REAL(wp) :: rn_zs ! depth of surface grid box
73	! bottom cell depth (Zb) is a linear function of water depth Zb = H*a + b
74	REAL(wp) :: rn_zb_a ! bathymetry scaling factor for calculating Zb
75	REAL(wp) :: rn_zb_b ! offset for calculating Zb
76
77	!! * Substitutions
78	!!----------------------------------------------------------------------
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	WHERE( bathy(:,:) <= 0._wp ) k_top(:,:) = 0 ! set k_top to zero over land
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 ( cn_topolvl, inum )
679	CALL iom_get ( inum, jpdom_data, cn_bathlvl, 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 ( cn_topo, inum )
717	CALL iom_get ( inum, jpdom_data, cn_bath, bathy, lrowattr=ln_use_jattr )
718	CALL iom_close( inum )
719	ELSE
720	CALL dom_bat
721	ENDIF
722	#if defined key_agrif
723	ELSE
724	IF( ntopo == 1) THEN
725	CALL agrif_create_bathy_meter()
726	ELSE
727	CALL dom_bat
728	ENDIF
729	ENDIF
730	#endif
731	!
732	! initialisation isf variables
733	risfdep(:,:)=0._wp ; misfdep(:,:)=1
734	!
735	IF ( ln_isfcav ) THEN
736	CALL iom_open ( cn_fisfd, inum )
737	CALL iom_get ( inum, jpdom_data, cn_visfd, risfdep )
738	CALL iom_close( inum )
739	END IF
740	!
741	IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration
742	!
743	ii0 = 140 ; ii1 = 140 ! Gibraltar Strait open
744	ij0 = 102 ; ij1 = 102 ! (Thomson, Ocean Modelling, 1995)
745	DO ji = mi0(ii0), mi1(ii1)
746	DO jj = mj0(ij0), mj1(ij1)
747	bathy(ji,jj) = 284._wp
748	END DO
749	END DO
750	IF(lwp) WRITE(numout,*)
751	IF(lwp) WRITE(numout,*) ' orca_r2: Gibraltar strait open at i=',ii0,' j=',ij0
752	!
753	ii0 = 160 ; ii1 = 160 ! Bab el mandeb Strait open
754	ij0 = 88 ; ij1 = 88 ! (Thomson, Ocean Modelling, 1995)
755	DO ji = mi0(ii0), mi1(ii1)
756	DO jj = mj0(ij0), mj1(ij1)
757	bathy(ji,jj) = 137._wp
758	END DO
759	END DO
760	IF(lwp) WRITE(numout,*)
761	IF(lwp) WRITE(numout,*) ' orca_r2: Bab el Mandeb strait open at i=',ii0,' j=',ij0
762	!
763	ENDIF
764	!
765	ENDIF
766	! ! =============== !
767	ELSE ! error !
768	! ! =============== !
769	WRITE(ctmp1,*) 'parameter , ntopo = ', ntopo
770	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
771	ENDIF
772	!
773	IF ( .not. ln_sco ) THEN !== set a minimum depth ==!
774	IF( rn_hmin < 0._wp ) THEN ; ik = - INT( rn_hmin ) ! from a nb of level
775	ELSE ; ik = MINLOC( gdepw_1d, mask = gdepw_1d > rn_hmin, dim = 1 ) ! from a depth
776	ENDIF
777	zhmin = gdepw_1d(ik+1) ! minimum depth = ik+1 w-levels
778	WHERE( bathy(:,:) <= 0._wp ) ; bathy(:,:) = 0._wp ! min=0 over the lands
779	ELSE WHERE ( risfdep == 0._wp ); bathy(:,:) = MAX( zhmin , bathy(:,:) ) ! min=zhmin over the oceans
780	END WHERE
781	IF(lwp) write(numout,*) 'Minimum ocean depth: ', zhmin, ' minimum number of ocean levels : ', ik
782	ENDIF
783	!
784	END SUBROUTINE zgr_bat
785
786	SUBROUTINE zgr_bat_ctl
787	!!----------------------------------------------------------------------
788	!! * ROUTINE zgr_bat_ctl *
789	!!
790	!! ** Purpose : check the bathymetry in levels
791	!!
792	!! ** Method : The array mbathy is checked to verified its consistency
793	!! with the model options. in particular:
794	!! mbathy must have at least 1 land grid-points (mbathy<=0)
795	!! along closed boundary.
796	!! mbathy must be cyclic IF jperio=1.
797	!! mbathy must be lower or equal to jpk-1.
798	!! isolated ocean grid points are suppressed from mbathy
799	!! since they are only connected to remaining
800	!! ocean through vertical diffusion.
801	!! C A U T I O N : mbathy will be modified during the initializa-
802	!! tion phase to become the number of non-zero w-levels of a water
803	!! column, with a minimum value of 1.
804	!!
805	!! ** Action : - update mbathy: level bathymetry (in level index)
806	!! - update bathy : meter bathymetry (in meters)
807	!!----------------------------------------------------------------------
808	INTEGER :: ji, jj, jl ! dummy loop indices
809	INTEGER :: icompt, ibtest, ikmax ! temporary integers
810	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zbathy
811	!!----------------------------------------------------------------------
812	!
813	ALLOCATE(zbathy(jpi,jpj))
814	!
815	IF(lwp) WRITE(numout,*)
816	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
817	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
818	! ! Suppress isolated ocean grid points
819	IF(lwp) WRITE(numout,*)
820	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
821	IF(lwp) WRITE(numout,*)' -----------------------------------'
822	icompt = 0
823	DO jl = 1, 2
824	IF( l_Iperio ) THEN
825	mbathy( 1 ,:) = mbathy(jpim1,:) ! local domain is cyclic east-west
826	mbathy(jpi,:) = mbathy( 2 ,:)
827	ENDIF
828	zbathy(:,:) = FLOAT( mbathy(:,:) )
829	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
830	mbathy(:,:) = INT( zbathy(:,:) )
831
832	DO jj = 2, jpjm1
833	DO ji = 2, jpim1
834	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
835	& mbathy(ji,jj-1), mbathy(ji,jj+1) )
836	IF( ibtest < mbathy(ji,jj) ) THEN
837	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
838	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
839	mbathy(ji,jj) = ibtest
840	icompt = icompt + 1
841	ENDIF
842	END DO
843	END DO
844	END DO
845
846	IF( lk_mpp ) CALL mpp_sum( 'domzgr', icompt )
847	IF( icompt == 0 ) THEN
848	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
849	ELSE
850	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
851	ENDIF
852
853	IF( lk_mpp ) THEN
854	zbathy(:,:) = FLOAT( mbathy(:,:) )
855	CALL lbc_lnk( 'toto',zbathy, 'T', 1._wp )
856	mbathy(:,:) = INT( zbathy(:,:) )
857	ENDIF
858	! ! East-west cyclic boundary conditions
859	IF( jperio == 0 ) THEN
860	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west boundary: jperio = ', jperio
861	IF( lk_mpp ) THEN
862	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
863	IF( jperio /= 1 ) mbathy(1,:) = 0
864	ENDIF
865	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
866	IF( jperio /= 1 ) mbathy(nlci,:) = 0
867	ENDIF
868	ELSE
869	IF( ln_zco .OR. ln_zps ) THEN
870	mbathy( 1 ,:) = 0
871	mbathy(jpi,:) = 0
872	ELSE
873	mbathy( 1 ,:) = jpkm1
874	mbathy(jpi,:) = jpkm1
875	ENDIF
876	ENDIF
877	ELSEIF( l_Iperio ) THEN
878	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions on mbathy: jperio = ', jperio
879	mbathy( 1 ,:) = mbathy(jpim1,:)
880	mbathy(jpi,:) = mbathy( 2 ,:)
881	ELSEIF( jperio == 2 ) THEN
882	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions on mbathy: jperio = ', jperio
883	ELSE
884	IF(lwp) WRITE(numout,*) ' e r r o r'
885	IF(lwp) WRITE(numout,*) ' parameter , jperio = ', jperio
886	! STOP 'dom_mba'
887	ENDIF
888
889	! Boundary condition on mbathy
890	IF( .NOT.lk_mpp ) THEN
891	!!gm !!bug ??? think about it !
892	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
893	zbathy(:,:) = FLOAT( mbathy(:,:) )
894	CALL lbc_lnk( 'domzgr',zbathy, 'T', 1._wp )
895	mbathy(:,:) = INT( zbathy(:,:) )
896	ENDIF
897
898	! Number of ocean level inferior or equal to jpkm1
899	zbathy(:,:) = FLOAT( mbathy(:,:) )
900	ikmax = glob_max( 'domzgr', zbathy(:,:) )
901
902	IF( ikmax > jpkm1 ) THEN
903	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
904	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
905	ELSE IF( ikmax < jpkm1 ) THEN
906	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
907	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
908	ENDIF
909	!
910	DEALLOCATE( zbathy )
911	!
912	END SUBROUTINE zgr_bat_ctl
913
914
915	SUBROUTINE zgr_bot_level
916	!!----------------------------------------------------------------------
917	!! * ROUTINE zgr_bot_level *
918	!!
919	!! ** Purpose : defines the vertical index of ocean bottom (mbk. arrays)
920	!!
921	!! ** Method : computes from mbathy with a minimum value of 1 over land
922	!!
923	!! ** Action : mbkt, mbku, mbkv : vertical indices of the deeptest
924	!! ocean level at t-, u- & v-points
925	!! (min value = 1 over land)
926	!!----------------------------------------------------------------------
927	INTEGER :: ji, jj ! dummy loop indices
928	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmbk
929	!!----------------------------------------------------------------------
930	!
931	ALLOCATE( zmbk(jpi,jpj) )
932	!
933	IF(lwp) WRITE(numout,*)
934	IF(lwp) WRITE(numout,*) ' zgr_bot_level : ocean bottom k-index of T-, U-, V- and W-levels '
935	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
936	!
937	mbkt(:,:) = MAX( mbathy(:,:) , 1 ) ! bottom k-index of T-level (=1 over land)
938
939	! ! bottom k-index of W-level = mbkt+1
940	DO jj = 1, jpjm1 ! bottom k-index of u- (v-) level
941	DO ji = 1, jpim1
942	mbku(ji,jj) = MIN( mbkt(ji+1,jj ) , mbkt(ji,jj) )
943	mbkv(ji,jj) = MIN( mbkt(ji ,jj+1) , mbkt(ji,jj) )
944	END DO
945	END DO
946	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
947	zmbk(:,:) = REAL( mbku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'U',1.) ; mbku (:,:) = MAX( INT( zmbk(:,:) ), 1 )
948	zmbk(:,:) = REAL( mbkv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmbk,'V',1.) ; mbkv (:,:) = MAX( INT( zmbk(:,:) ), 1 )
949	!
950	DEALLOCATE( zmbk )
951	!
952	END SUBROUTINE zgr_bot_level
953
954
955	SUBROUTINE zgr_top_level
956	!!----------------------------------------------------------------------
957	!! * ROUTINE zgr_top_level *
958	!!
959	!! ** Purpose : defines the vertical index of ocean top (mik. arrays)
960	!!
961	!! ** Method : computes from misfdep with a minimum value of 1
962	!!
963	!! ** Action : mikt, miku, mikv : vertical indices of the shallowest
964	!! ocean level at t-, u- & v-points
965	!! (min value = 1)
966	!!----------------------------------------------------------------------
967	INTEGER :: ji, jj ! dummy loop indices
968	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zmik
969	!!----------------------------------------------------------------------
970	!
971	ALLOCATE( zmik(jpi,jpj) )
972	!
973	IF(lwp) WRITE(numout,*)
974	IF(lwp) WRITE(numout,*) ' zgr_top_level : ocean top k-index of T-, U-, V- and W-levels '
975	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~'
976	!
977	mikt(:,:) = MAX( misfdep(:,:) , 1 ) ! top k-index of T-level (=1)
978	! ! top k-index of W-level (=mikt)
979	DO jj = 1, jpjm1 ! top k-index of U- (U-) level
980	DO ji = 1, jpim1
981	miku(ji,jj) = MAX( mikt(ji+1,jj ) , mikt(ji,jj) )
982	mikv(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj) )
983	mikf(ji,jj) = MAX( mikt(ji ,jj+1) , mikt(ji,jj), mikt(ji+1,jj ), mikt(ji+1,jj+1) )
984	END DO
985	END DO
986
987	! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
988	zmik(:,:) = REAL( miku(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'U',1.) ; miku (:,:) = MAX( INT( zmik(:,:) ), 1 )
989	zmik(:,:) = REAL( mikv(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'V',1.) ; mikv (:,:) = MAX( INT( zmik(:,:) ), 1 )
990	zmik(:,:) = REAL( mikf(:,:), wp ) ; CALL lbc_lnk('domzgr',zmik,'F',1.) ; mikf (:,:) = MAX( INT( zmik(:,:) ), 1 )
991	!
992	DEALLOCATE( zmik )
993	!
994	END SUBROUTINE zgr_top_level
995
996
997	SUBROUTINE zgr_zco
998	!!----------------------------------------------------------------------
999	!! * ROUTINE zgr_zco *
1000	!!
1001	!! ** Purpose : define the reference z-coordinate system
1002	!!
1003	!! ** Method : set 3D coord. arrays to reference 1D array
1004	!!----------------------------------------------------------------------
1005	INTEGER :: jk
1006	!!----------------------------------------------------------------------
1007	!
1008	DO jk = 1, jpk
1009	gdept_0(:,:,jk) = gdept_1d(jk)
1010	gdepw_0(:,:,jk) = gdepw_1d(jk)
1011	e3t_0 (:,:,jk) = e3t_1d (jk)
1012	e3u_0 (:,:,jk) = e3t_1d (jk)
1013	e3v_0 (:,:,jk) = e3t_1d (jk)
1014	e3f_0 (:,:,jk) = e3t_1d (jk)
1015	e3w_0 (:,:,jk) = e3w_1d (jk)
1016	e3uw_0 (:,:,jk) = e3w_1d (jk)
1017	e3vw_0 (:,:,jk) = e3w_1d (jk)
1018	END DO
1019	!
1020	END SUBROUTINE zgr_zco
1021
1022
1023	SUBROUTINE zgr_zps
1024	!!----------------------------------------------------------------------
1025	!! * ROUTINE zgr_zps *
1026	!!
1027	!! ** Purpose : the depth and vertical scale factor in partial step
1028	!! reference z-coordinate case
1029	!!
1030	!! ** Method : Partial steps : computes the 3D vertical scale factors
1031	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
1032	!! a partial step representation of bottom topography.
1033	!!
1034	!! The reference depth of model levels is defined from an analytical
1035	!! function the derivative of which gives the reference vertical
1036	!! scale factors.
1037	!! From depth and scale factors reference, we compute there new value
1038	!! with partial steps on 3d arrays ( i, j, k ).
1039	!!
1040	!! w-level: gdepw_0(i,j,k) = gdep(k)
1041	!! e3w_0(i,j,k) = dk(gdep)(k) = e3(i,j,k)
1042	!! t-level: gdept_0(i,j,k) = gdep(k+0.5)
1043	!! e3t_0(i,j,k) = dk(gdep)(k+0.5) = e3(i,j,k+0.5)
1044	!!
1045	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
1046	!! we find the mbathy index of the depth at each grid point.
1047	!! This leads us to three cases:
1048	!!
1049	!! - bathy = 0 => mbathy = 0
1050	!! - 1 < mbathy < jpkm1
1051	!! - bathy > gdepw_0(jpk) => mbathy = jpkm1
1052	!!
1053	!! Then, for each case, we find the new depth at t- and w- levels
1054	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
1055	!! and f-points.
1056	!!
1057	!! This routine is given as an example, it must be modified
1058	!! following the user s desiderata. nevertheless, the output as
1059	!! well as the way to compute the model levels and scale factors
1060	!! must be respected in order to insure second order accuracy
1061	!! schemes.
1062	!!
1063	!! c a u t i o n : gdept_1d, gdepw_1d and e3._1d are positives
1064	!! - - - - - - - gdept_0, gdepw_0 and e3. are positives
1065	!!
1066	!! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
1067	!!----------------------------------------------------------------------
1068	INTEGER :: ji, jj, jk ! dummy loop indices
1069	INTEGER :: ik, it, ikb, ikt ! temporary integers
1070	REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points
1071	REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t
1072	REAL(wp) :: zdiff ! temporary scalar
1073	REAL(wp) :: zmax ! temporary scalar
1074	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zprt
1075	!!---------------------------------------------------------------------
1076	!
1077	ALLOCATE( zprt(jpi,jpj,jpk) )
1078	!
1079	IF(lwp) WRITE(numout,*)
1080	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
1081	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
1082	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
1083
1084	! compute position of the ice shelf grounding line
1085	! set bathy and isfdraft to 0 where grounded
1086	IF ( ln_isfcav ) CALL zgr_isf_zspace
1087
1088	! bathymetry in level (from bathy_meter)
1089	! ===================
1090	zmax = gdepw_1d(jpk) + e3t_1d(jpk) ! maximum depth (i.e. the last ocean level thickness <= 2*e3t_1d(jpkm1) )
1091	bathy(:,:) = MIN( zmax , bathy(:,:) ) ! bounded value of bathy (min already set at the end of zgr_bat)
1092	WHERE( bathy(:,:) == 0._wp ) ; mbathy(:,:) = 0 ! land : set mbathy to 0
1093	ELSE WHERE ; mbathy(:,:) = jpkm1 ! ocean : initialize mbathy to the max ocean level
1094	END WHERE
1095
1096	! Compute mbathy for ocean points (i.e. the number of ocean levels)
1097	! find the number of ocean levels such that the last level thickness
1098	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_1d (where
1099	! e3t_1d is the reference level thickness
1100	DO jk = jpkm1, 1, -1
1101	zdepth = gdepw_1d(jk) + MIN( e3zps_min, e3t_1d(jk)*e3zps_rat )
1102	WHERE( 0._wp < bathy(:,:) .AND. bathy(:,:) <= zdepth ) mbathy(:,:) = jk-1
1103	END DO
1104
1105	! Check compatibility between bathy and iceshelf draft
1106	! insure at least 2 wet level on the vertical under an ice shelf
1107	! compute misfdep and adjust isf draft if needed
1108	IF ( ln_isfcav ) CALL zgr_isf_kspace
1109
1110	! Scale factors and depth at T- and W-points
1111	DO jk = 1, jpk ! intitialization to the reference z-coordinate
1112	gdept_0(:,:,jk) = gdept_1d(jk)
1113	gdepw_0(:,:,jk) = gdepw_1d(jk)
1114	e3t_0 (:,:,jk) = e3t_1d (jk)
1115	e3w_0 (:,:,jk) = e3w_1d (jk)
1116	END DO
1117
1118	! Scale factors and depth at T- and W-points
1119	DO jj = 1, jpj
1120	DO ji = 1, jpi
1121	ik = mbathy(ji,jj)
1122	IF( ik > 0 ) THEN ! ocean point only
1123	! max ocean level case
1124	IF( ik == jpkm1 ) THEN
1125	zdepwp = bathy(ji,jj)
1126	ze3tp = bathy(ji,jj) - gdepw_1d(ik)
1127	ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) )
1128	e3t_0(ji,jj,ik ) = ze3tp
1129	e3t_0(ji,jj,ik+1) = ze3tp
1130	e3w_0(ji,jj,ik ) = ze3wp
1131	e3w_0(ji,jj,ik+1) = ze3tp
1132	gdepw_0(ji,jj,ik+1) = zdepwp
1133	gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp
1134	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp
1135	!
1136	ELSE ! standard case
1137	IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj)
1138	ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1)
1139	ENDIF
1140	!gm Bug? check the gdepw_1d
1141	! ... on ik
1142	gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &
1143	& * ((gdept_1d( ik ) - gdepw_1d(ik) ) &
1144	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))
1145	e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &
1146	& / ( gdepw_1d( ik+1) - gdepw_1d(ik) )
1147	e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &
1148	& * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )
1149	! ... on ik+1
1150	e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1151	e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)
1152	gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)
1153	ENDIF
1154	ENDIF
1155	END DO
1156	END DO
1157	!
1158	it = 0
1159	DO jj = 1, jpj
1160	DO ji = 1, jpi
1161	ik = mbathy(ji,jj)
1162	IF( ik > 0 ) THEN ! ocean point only
1163	e3tp (ji,jj) = e3t_0(ji,jj,ik)
1164	e3wp (ji,jj) = e3w_0(ji,jj,ik)
1165	! test
1166	zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )
1167	IF( zdiff <= 0._wp .AND. lwp ) THEN
1168	it = it + 1
1169	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
1170	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
1171	WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff
1172	WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )
1173	ENDIF
1174	ENDIF
1175	END DO
1176	END DO
1177	!
1178	! compute top scale factor if ice shelf
1179	IF (ln_isfcav) CALL zps_isf
1180	!
1181	! Scale factors and depth at U-, V-, UW and VW-points
1182	DO jk = 1, jpk ! initialisation to z-scale factors
1183	e3u_0 (:,:,jk) = e3t_1d(jk)
1184	e3v_0 (:,:,jk) = e3t_1d(jk)
1185	e3uw_0(:,:,jk) = e3w_1d(jk)
1186	e3vw_0(:,:,jk) = e3w_1d(jk)
1187	END DO
1188
1189	DO jk = 1,jpk ! Computed as the minimum of neighbooring scale factors
1190	DO jj = 1, jpjm1
1191	DO ji = 1, jpim1 ! vector opt.
1192	e3u_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji+1,jj,jk) )
1193	e3v_0 (ji,jj,jk) = MIN( e3t_0(ji,jj,jk), e3t_0(ji,jj+1,jk) )
1194	e3uw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji+1,jj,jk) )
1195	e3vw_0(ji,jj,jk) = MIN( e3w_0(ji,jj,jk), e3w_0(ji,jj+1,jk) )
1196	END DO
1197	END DO
1198	END DO
1199
1200	! update e3uw in case only 2 cells in the water column
1201	IF ( ln_isfcav ) CALL zps_isf_e3uv_w
1202	!
1203	CALL lbc_lnk('domzgr', e3u_0 , 'U', 1._wp ) ; CALL lbc_lnk('domzgr', e3uw_0, 'U', 1._wp ) ! lateral boundary conditions
1204	CALL lbc_lnk('domzgr', e3v_0 , 'V', 1._wp ) ; CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1205	!
1206	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1207	WHERE( e3u_0 (:,:,jk) == 0._wp ) e3u_0 (:,:,jk) = e3t_1d(jk)
1208	WHERE( e3v_0 (:,:,jk) == 0._wp ) e3v_0 (:,:,jk) = e3t_1d(jk)
1209	WHERE( e3uw_0(:,:,jk) == 0._wp ) e3uw_0(:,:,jk) = e3w_1d(jk)
1210	WHERE( e3vw_0(:,:,jk) == 0._wp ) e3vw_0(:,:,jk) = e3w_1d(jk)
1211	END DO
1212
1213	! Scale factor at F-point
1214	DO jk = 1, jpk ! initialisation to z-scale factors
1215	e3f_0(:,:,jk) = e3t_1d(jk)
1216	END DO
1217	DO jk = 1, jpk ! Computed as the minimum of neighbooring V-scale factors
1218	DO jj = 1, jpjm1
1219	DO ji = 1, jpim1 ! vector opt.
1220	e3f_0(ji,jj,jk) = MIN( e3v_0(ji,jj,jk), e3v_0(ji+1,jj,jk) )
1221	END DO
1222	END DO
1223	END DO
1224	CALL lbc_lnk('domzgr', e3f_0, 'F', 1._wp ) ! Lateral boundary conditions
1225	!
1226	DO jk = 1, jpk ! set to z-scale factor if zero (i.e. along closed boundaries)
1227	WHERE( e3f_0(:,:,jk) == 0._wp ) e3f_0(:,:,jk) = e3t_1d(jk)
1228	END DO
1229	!!gm bug ? : must be a do loop with mj0,mj1
1230	!
1231	e3t_0(:,mj0(1),:) = e3t_0(:,mj0(2),:) ! we duplicate factor scales for jj = 1 and jj = 2
1232	e3w_0(:,mj0(1),:) = e3w_0(:,mj0(2),:)
1233	e3u_0(:,mj0(1),:) = e3u_0(:,mj0(2),:)
1234	e3v_0(:,mj0(1),:) = e3v_0(:,mj0(2),:)
1235	e3f_0(:,mj0(1),:) = e3f_0(:,mj0(2),:)
1236
1237	! Control of the sign
1238	IF( MINVAL( e3t_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3t_0 <= 0' )
1239	IF( MINVAL( e3w_0 (:,:,:) ) <= 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r e3w_0 <= 0' )
1240	IF( MINVAL( gdept_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdept_0 < 0' )
1241	IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' )
1242	!
1243	! if in the future gde3w_0 need to be compute, use the function defined in NEMO
1244	! for now gde3w_0 computation is removed as not an output of domcfg
1245
1246	DEALLOCATE( zprt )
1247	!
1248	END SUBROUTINE zgr_zps
1249
1250
1251	SUBROUTINE zgr_sco
1252	!!----------------------------------------------------------------------
1253	!! * ROUTINE zgr_sco *
1254	!!
1255	!! ** Purpose : define the s-coordinate system
1256	!!
1257	!! ** Method : s-coordinate
1258	!! The depth of model levels is defined as the product of an
1259	!! analytical function by the local bathymetry, while the vertical
1260	!! scale factors are defined as the product of the first derivative
1261	!! of the analytical function by the bathymetry.
1262	!! (this solution save memory as depth and scale factors are not
1263	!! 3d fields)
1264	!! - Read bathymetry (in meters) at t-point and compute the
1265	!! bathymetry at u-, v-, and f-points.
1266	!! hbatu = mi( hbatt )
1267	!! hbatv = mj( hbatt )
1268	!! hbatf = mi( mj( hbatt ) )
1269	!! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical
1270	!! function and its derivative given as function.
1271	!! z_gsigt(k) = fssig (k )
1272	!! z_gsigw(k) = fssig (k-0.5)
1273	!! z_esigt(k) = fsdsig(k )
1274	!! z_esigw(k) = fsdsig(k-0.5)
1275	!! Three options for stretching are give, and they can be modified
1276	!! following the users requirements. Nevertheless, the output as
1277	!! well as the way to compute the model levels and scale factors
1278	!! must be respected in order to insure second order accuracy
1279	!! schemes.
1280	!!
1281	!! The three methods for stretching available are:
1282	!!
1283	!! s_sh94 (Song and Haidvogel 1994)
1284	!! a sinh/tanh function that allows sigma and stretched sigma
1285	!!
1286	!! s_sf12 (Siddorn and Furner 2012?)
1287	!! allows the maintenance of fixed surface and or
1288	!! bottom cell resolutions (cf. geopotential coordinates)
1289	!! within an analytically derived stretched S-coordinate framework.
1290	!!
1291	!! s_tanh (Madec et al 1996)
1292	!! a cosh/tanh function that gives stretched coordinates
1293	!!
1294	!!----------------------------------------------------------------------
1295	INTEGER :: ji, jj, jk, jl ! dummy loop argument
1296	INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers
1297	INTEGER :: ios ! Local integer output status for namelist read
1298	REAL(wp) :: zrmax, ztaper ! temporary scalars
1299	REAL(wp) :: zrfact
1300	!
1301	REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2
1302	REAL(wp), ALLOCATABLE, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat
1303	!!
1304	NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, &
1305	& rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b
1306	!!----------------------------------------------------------------------
1307	!
1308	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) )
1309	!
1310	REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters
1311	READ ( numnam_ref, namzgr_sco, IOSTAT = ios, ERR = 901)
1312	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in reference namelist', lwp )
1313
1314	REWIND( numnam_cfg ) ! Namelist namzgr_sco in configuration namelist : Sigma-stretching parameters
1315	READ ( numnam_cfg, namzgr_sco, IOSTAT = ios, ERR = 902 )
1316	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_sco in configuration namelist', lwp )
1317	IF(lwm) WRITE ( numond, namzgr_sco )
1318
1319	IF(lwp) THEN ! control print
1320	WRITE(numout,*)
1321	WRITE(numout,*) 'domzgr_sco : s-coordinate or hybrid z-s-coordinate'
1322	WRITE(numout,*) '~~~~~~~~~~~'
1323	WRITE(numout,*) ' Namelist namzgr_sco'
1324	WRITE(numout,*) ' stretching coeffs '
1325	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) rn_sbot_max = ',rn_sbot_max
1326	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) rn_sbot_min = ',rn_sbot_min
1327	WRITE(numout,*) ' Critical depth rn_hc = ',rn_hc
1328	WRITE(numout,*) ' maximum cut-off r-value allowed rn_rmax = ',rn_rmax
1329	WRITE(numout,*) ' Song and Haidvogel 1994 stretching ln_s_sh94 = ',ln_s_sh94
1330	WRITE(numout,*) ' Song and Haidvogel 1994 stretching coefficients'
1331	WRITE(numout,*) ' surface control parameter (0<=rn_theta<=20) rn_theta = ',rn_theta
1332	WRITE(numout,*) ' bottom control parameter (0<=rn_thetb<= 1) rn_thetb = ',rn_thetb
1333	WRITE(numout,*) ' stretching parameter (song and haidvogel) rn_bb = ',rn_bb
1334	WRITE(numout,*) ' Siddorn and Furner 2012 stretching ln_s_sf12 = ',ln_s_sf12
1335	WRITE(numout,*) ' switching to sigma (T) or Z (F) at H<Hc ln_sigcrit = ',ln_sigcrit
1336	WRITE(numout,*) ' Siddorn and Furner 2012 stretching coefficients'
1337	WRITE(numout,*) ' stretchin parameter ( >1 surface; <1 bottom) rn_alpha = ',rn_alpha
1338	WRITE(numout,*) ' e-fold length scale for transition region rn_efold = ',rn_efold
1339	WRITE(numout,*) ' Surface cell depth (Zs) (m) rn_zs = ',rn_zs
1340	WRITE(numout,*) ' Bathymetry multiplier for Zb rn_zb_a = ',rn_zb_a
1341	WRITE(numout,*) ' Offset for Zb rn_zb_b = ',rn_zb_b
1342	WRITE(numout,) ' Bottom cell (Zb) (m) = Hrn_zb_a + rn_zb_b'
1343	ENDIF
1344
1345	hift(:,:) = rn_sbot_min ! set the minimum depth for the s-coordinate
1346	hifu(:,:) = rn_sbot_min
1347	hifv(:,:) = rn_sbot_min
1348	hiff(:,:) = rn_sbot_min
1349
1350	! ! set maximum ocean depth
1351	bathy(:,:) = MIN( rn_sbot_max, bathy(:,:) )
1352
1353	DO jj = 1, jpj
1354	DO ji = 1, jpi
1355	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_sbot_min, bathy(ji,jj) )
1356	END DO
1357	END DO
1358	! ! =============================
1359	! ! Define the envelop bathymetry (hbatt)
1360	! ! =============================
1361	! use r-value to create hybrid coordinates
1362	zenv(:,:) = bathy(:,:)
1363	!
1364	! set first land point adjacent to a wet cell to sbot_min as this needs to be included in smoothing
1365	DO jj = 1, jpj
1366	DO ji = 1, jpi
1367	IF( bathy(ji,jj) == 0._wp ) THEN
1368	iip1 = MIN( ji+1, jpi )
1369	ijp1 = MIN( jj+1, jpj )
1370	iim1 = MAX( ji-1, 1 )
1371	ijm1 = MAX( jj-1, 1 )
1372	!!gm BUG fix see ticket #1617
1373	IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) &
1374	& + bathy(iim1,jj ) + bathy(iip1,jj ) &
1375	& + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) &
1376	& zenv(ji,jj) = rn_sbot_min
1377	!!gm
1378	!!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + &
1379	!!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN
1380	!!gm zenv(ji,jj) = rn_sbot_min
1381	!!gm ENDIF
1382	!!gm end
1383	ENDIF
1384	END DO
1385	END DO
1386
1387	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1388	CALL lbc_lnk( 'domzgr',zenv, 'T', 1._wp, 'no0' )
1389	!
1390	! smooth the bathymetry (if required)
1391	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
1392	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1393	!
1394	jl = 0
1395	zrmax = 1._wp
1396	!
1397	!
1398	! set scaling factor used in reducing vertical gradients
1399	zrfact = ( 1._wp - rn_rmax ) / ( 1._wp + rn_rmax )
1400	!
1401	! initialise temporary evelope depth arrays
1402	ztmpi1(:,:) = zenv(:,:)
1403	ztmpi2(:,:) = zenv(:,:)
1404	ztmpj1(:,:) = zenv(:,:)
1405	ztmpj2(:,:) = zenv(:,:)
1406	!
1407	! initialise temporary r-value arrays
1408	zri(:,:) = 1._wp
1409	zrj(:,:) = 1._wp
1410	! ! ================ !
1411	DO WHILE( jl <= 10000 .AND. ( zrmax - rn_rmax ) > 1.e-8_wp ) ! Iterative loop !
1412	! ! ================ !
1413	jl = jl + 1
1414	zrmax = 0._wp
1415	! we set zrmax from previous r-values (zri and zrj) first
1416	! if set after current r-value calculation (as previously)
1417	! we could exit DO WHILE prematurely before checking r-value
1418	! of current zenv
1419	DO jj = 1, nlcj
1420	DO ji = 1, nlci
1421	zrmax = MAX( zrmax, ABS(zri(ji,jj)), ABS(zrj(ji,jj)) )
1422	END DO
1423	END DO
1424	zri(:,:) = 0._wp
1425	zrj(:,:) = 0._wp
1426	DO jj = 1, nlcj
1427	DO ji = 1, nlci
1428	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1429	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1430	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(iip1,jj) > 0._wp)) THEN
1431	zri(ji,jj) = ( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1432	END IF
1433	IF( (zenv(ji,jj) > 0._wp) .AND. (zenv(ji,ijp1) > 0._wp)) THEN
1434	zrj(ji,jj) = ( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1435	END IF
1436	IF( zri(ji,jj) > rn_rmax ) ztmpi1(ji ,jj ) = zenv(iip1,jj ) * zrfact
1437	IF( zri(ji,jj) < -rn_rmax ) ztmpi2(iip1,jj ) = zenv(ji ,jj ) * zrfact
1438	IF( zrj(ji,jj) > rn_rmax ) ztmpj1(ji ,jj ) = zenv(ji ,ijp1) * zrfact
1439	IF( zrj(ji,jj) < -rn_rmax ) ztmpj2(ji ,ijp1) = zenv(ji ,jj ) * zrfact
1440	END DO
1441	END DO
1442	! IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1443	!
1444	IF(lwp)WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax
1445	!
1446	DO jj = 1, nlcj
1447	DO ji = 1, nlci
1448	zenv(ji,jj) = MAX(zenv(ji,jj), ztmpi1(ji,jj), ztmpi2(ji,jj), ztmpj1(ji,jj), ztmpj2(ji,jj) )
1449	END DO
1450	END DO
1451	! apply lateral boundary condition CAUTION: keep the value when the lbc field is zero
1452	CALL lbc_lnk( 'toto',zenv, 'T', 1._wp, 'no0' )
1453	! ! ================ !
1454	END DO ! End loop !
1455	! ! ================ !
1456	DO jj = 1, jpj
1457	DO ji = 1, jpi
1458	zenv(ji,jj) = MAX( zenv(ji,jj), rn_sbot_min ) ! set all points to avoid undefined scale value warnings
1459	END DO
1460	END DO
1461	!
1462	! Envelope bathymetry saved in hbatt
1463	hbatt(:,:) = zenv(:,:)
1464	IF( MINVAL( gphit(:,:) ) * MAXVAL( gphit(:,:) ) <= 0._wp ) THEN
1465	CALL ctl_warn( ' s-coordinates are tapered in vicinity of the Equator' )
1466	DO jj = 1, jpj
1467	DO ji = 1, jpi
1468	ztaper = EXP( -(gphit(ji,jj)/8._wp)**2._wp )
1469	hbatt(ji,jj) = rn_sbot_max * ztaper + hbatt(ji,jj) * ( 1._wp - ztaper )
1470	END DO
1471	END DO
1472	ENDIF
1473	!
1474	! ! ==============================
1475	! ! hbatu, hbatv, hbatf fields
1476	! ! ==============================
1477	IF(lwp) THEN
1478	WRITE(numout,*)
1479	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', rn_sbot_min
1480	ENDIF
1481	hbatu(:,:) = rn_sbot_min
1482	hbatv(:,:) = rn_sbot_min
1483	hbatf(:,:) = rn_sbot_min
1484	DO jj = 1, jpjm1
1485	DO ji = 1, jpim1 ! NO vector opt.
1486	hbatu(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji+1,jj ) )
1487	hbatv(ji,jj) = 0.50_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) )
1488	hbatf(ji,jj) = 0.25_wp * ( hbatt(ji ,jj) + hbatt(ji ,jj+1) &
1489	& + hbatt(ji+1,jj) + hbatt(ji+1,jj+1) )
1490	END DO
1491	END DO
1492
1493	!
1494	! Apply lateral boundary condition
1495	!!gm ! CAUTION: retain non zero value in the initial file this should be OK for orca cfg, not for EEL
1496	zhbat(:,:) = hbatu(:,:) ; CALL lbc_lnk('domzgr', hbatu, 'U', 1._wp )
1497	DO jj = 1, jpj
1498	DO ji = 1, jpi
1499	IF( hbatu(ji,jj) == 0._wp ) THEN
1500	!No worries about the following line when ln_wd == .true.
1501	IF( zhbat(ji,jj) == 0._wp ) hbatu(ji,jj) = rn_sbot_min
1502	IF( zhbat(ji,jj) /= 0._wp ) hbatu(ji,jj) = zhbat(ji,jj)
1503	ENDIF
1504	END DO
1505	END DO
1506	zhbat(:,:) = hbatv(:,:) ; CALL lbc_lnk('domzgr', hbatv, 'V', 1._wp )
1507	DO jj = 1, jpj
1508	DO ji = 1, jpi
1509	IF( hbatv(ji,jj) == 0._wp ) THEN
1510	IF( zhbat(ji,jj) == 0._wp ) hbatv(ji,jj) = rn_sbot_min
1511	IF( zhbat(ji,jj) /= 0._wp ) hbatv(ji,jj) = zhbat(ji,jj)
1512	ENDIF
1513	END DO
1514	END DO
1515	zhbat(:,:) = hbatf(:,:) ; CALL lbc_lnk('domzgr', hbatf, 'F', 1._wp )
1516	DO jj = 1, jpj
1517	DO ji = 1, jpi
1518	IF( hbatf(ji,jj) == 0._wp ) THEN
1519	IF( zhbat(ji,jj) == 0._wp ) hbatf(ji,jj) = rn_sbot_min
1520	IF( zhbat(ji,jj) /= 0._wp ) hbatf(ji,jj) = zhbat(ji,jj)
1521	ENDIF
1522	END DO
1523	END DO
1524
1525	!!bug: key_helsinki a verifer
1526	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1527	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1528	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1529	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1530
1531	IF( lwp ) THEN
1532	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift (:,:) ), ' f ', MAXVAL( hiff (:,:) ), &
1533	& ' u ', MAXVAL( hifu (:,:) ), ' v ', MAXVAL( hifv (:,:) )
1534	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift (:,:) ), ' f ', MINVAL( hiff (:,:) ), &
1535	& ' u ', MINVAL( hifu (:,:) ), ' v ', MINVAL( hifv (:,:) )
1536	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1537	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1538	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1539	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1540	ENDIF
1541	!! helsinki
1542
1543	! ! =======================
1544	! ! s-ccordinate fields (gdep., e3.)
1545	! ! =======================
1546	!
1547	! non-dimensional "sigma" for model level depth at w- and t-levels
1548
1549
1550	!========================================================================
1551	! Song and Haidvogel 1994 (ln_s_sh94=T)
1552	! Siddorn and Furner 2012 (ln_sf12=T)
1553	! or tanh function (both false)
1554	!========================================================================
1555	IF ( ln_s_sh94 ) THEN
1556	CALL s_sh94()
1557	ELSE IF ( ln_s_sf12 ) THEN
1558	CALL s_sf12()
1559	ELSE
1560	CALL s_tanh()
1561	ENDIF
1562
1563	CALL lbc_lnk( 'domzgr',e3t_0 , 'T', 1._wp )
1564	CALL lbc_lnk( 'domzgr',e3u_0 , 'U', 1._wp )
1565	CALL lbc_lnk( 'domzgr',e3v_0 , 'V', 1._wp )
1566	CALL lbc_lnk( 'domzgr',e3f_0 , 'F', 1._wp )
1567	CALL lbc_lnk( 'domzgr',e3w_0 , 'W', 1._wp )
1568	CALL lbc_lnk( 'domzgr',e3uw_0, 'U', 1._wp )
1569	CALL lbc_lnk('domzgr', e3vw_0, 'V', 1._wp )
1570	!
1571	WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp
1572	WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp
1573	WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp
1574	WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp
1575	WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp
1576	WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp
1577	WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp
1578	!!
1579	! HYBRID :
1580	DO jj = 1, jpj
1581	DO ji = 1, jpi
1582	DO jk = 1, jpkm1
1583	IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1584	END DO
1585	END DO
1586	END DO
1587	IF(lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
1588	& ' MAX ', MAXVAL( mbathy(:,:) )
1589
1590	IF( lwp ) THEN ! min max values over the local domain
1591	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1592	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
1593	& ' w ', MINVAL( gdepw_0(:,:,:) )
1594	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
1595	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
1596	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
1597	& ' w ', MINVAL( e3w_0 (:,:,:) )
1598
1599	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
1600	& ' w ', MAXVAL( gdepw_0(:,:,:) )
1601	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
1602	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
1603	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
1604	& ' w ', MAXVAL( e3w_0 (:,:,:) )
1605	ENDIF
1606	! END DO
1607	IF(lwp) THEN ! selected vertical profiles
1608	WRITE(numout,*)
1609	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1610	WRITE(numout,*) ' ~~~~~~ --------------------'
1611	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
1612	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(1,1,jk), gdepw_0(1,1,jk), &
1613	& e3t_0 (1,1,jk) , e3w_0 (1,1,jk) , jk=1,jpk )
1614	DO jj = mj0(20), mj1(20)
1615	DO ji = mi0(20), mi1(20)
1616	WRITE(numout,*)
1617	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1618	WRITE(numout,*) ' ~~~~~~ --------------------'
1619	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
1620	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
1621	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1622	END DO
1623	END DO
1624	DO jj = mj0(74), mj1(74)
1625	DO ji = mi0(100), mi1(100)
1626	WRITE(numout,*)
1627	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1628	WRITE(numout,*) ' ~~~~~~ --------------------'
1629	WRITE(numout,"(9x,' level gdept_0 gdepw_0 e3t_0 e3w_0')")
1630	WRITE(numout,"(10x,i4,4f9.2)") ( jk, gdept_0(ji,jj,jk), gdepw_0(ji,jj,jk), &
1631	& e3t_0 (ji,jj,jk) , e3w_0 (ji,jj,jk) , jk=1,jpk )
1632	END DO
1633	END DO
1634	ENDIF
1635	!
1636	!================================================================================
1637	! check the coordinate makes sense
1638	!================================================================================
1639	DO ji = 1, jpi
1640	DO jj = 1, jpj
1641	!
1642	IF( hbatt(ji,jj) > 0._wp) THEN
1643	DO jk = 1, mbathy(ji,jj)
1644	! check coordinate is monotonically increasing
1645	IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
1646	WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1647	WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
1648	WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
1649	WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
1650	CALL ctl_stop( ctmp1 )
1651	ENDIF
1652	! and check it has never gone negative
1653	IF( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
1654	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1655	WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
1656	WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1657	WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1658	CALL ctl_stop( ctmp1 )
1659	ENDIF
1660	! and check it never exceeds the total depth
1661	IF( gdepw_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1662	WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
1663	WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk
1664	WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
1665	CALL ctl_stop( ctmp1 )
1666	ENDIF
1667	END DO
1668	!
1669	DO jk = 1, mbathy(ji,jj)-1
1670	! and check it never exceeds the total depth
1671	IF( gdept_0(ji,jj,jk) > hbatt(ji,jj) ) THEN
1672	WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
1673	WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk
1674	WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
1675	CALL ctl_stop( ctmp1 )
1676	ENDIF
1677	END DO
1678	ENDIF
1679	END DO
1680	END DO
1681	!
1682	DEALLOCATE( zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )
1683	!
1684	END SUBROUTINE zgr_sco
1685
1686
1687	SUBROUTINE s_sh94()
1688	!!----------------------------------------------------------------------
1689	!! * ROUTINE s_sh94 *
1690	!!
1691	!! ** Purpose : stretch the s-coordinate system
1692	!!
1693	!! ** Method : s-coordinate stretch using the Song and Haidvogel 1994
1694	!! mixed S/sigma coordinate
1695	!!
1696	!! Reference : Song and Haidvogel 1994.
1697	!!----------------------------------------------------------------------
1698	INTEGER :: ji, jj, jk ! dummy loop argument
1699	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
1700	REAL(wp) :: ztmpu, ztmpv, ztmpf
1701	REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
1702	!
1703	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1704	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
1705	!!----------------------------------------------------------------------
1706
1707	ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1708	ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk) )
1709	ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1710
1711	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp
1712	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
1713	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
1714	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
1715	!
1716	DO ji = 1, jpi
1717	DO jj = 1, jpj
1718	!
1719	IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma
1720	DO jk = 1, jpk
1721	z_gsigw3(ji,jj,jk) = -fssig1( REAL(jk,wp)-0.5_wp, rn_bb )
1722	z_gsigt3(ji,jj,jk) = -fssig1( REAL(jk,wp) , rn_bb )
1723	END DO
1724	ELSE ! shallow water, uniform sigma
1725	DO jk = 1, jpk
1726	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) / REAL(jpk-1,wp)
1727	z_gsigt3(ji,jj,jk) = ( REAL(jk-1,wp) + 0.5_wp ) / REAL(jpk-1,wp)
1728	END DO
1729	ENDIF
1730	!
1731	DO jk = 1, jpkm1
1732	z_esigt3(ji,jj,jk ) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1733	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1734	END DO
1735	z_esigw3(ji,jj,1 ) = 2._wp * ( z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ) )
1736	z_esigt3(ji,jj,jpk) = 2._wp * ( z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk) )
1737	!
1738	DO jk = 1, jpk
1739	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1740	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1741	gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigt3(ji,jj,jk)+rn_hczcoeft )
1742	gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)z_gsigw3(ji,jj,jk)+rn_hczcoefw )
1743	END DO
1744	!
1745	END DO ! for all jj's
1746	END DO ! for all ji's
1747
1748	DO ji = 1, jpim1
1749	DO jj = 1, jpjm1
1750	! extended for Wetting/Drying case
1751	ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
1752	ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
1753	ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1754	ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1755	ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1756	ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1757	& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1758	DO jk = 1, jpk
1759	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1760	& / ztmpu
1761	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1762	& / ztmpu
1763
1764	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1765	& / ztmpv
1766	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1767	& / ztmpv
1768
1769	z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
1770	& + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
1771	& + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
1772	& + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1773
1774	!
1775	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigt3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1776	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigtu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1777	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigtv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1778	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-rn_hc)*z_esigtf3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1779	!
1780	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-rn_hc)*z_esigw3 (ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1781	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-rn_hc)*z_esigwu3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1782	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-rn_hc)*z_esigwv3(ji,jj,jk) + rn_hc/REAL(jpkm1,wp) )
1783	END DO
1784	END DO
1785	END DO
1786	!
1787	DEALLOCATE( z_gsigw3, z_gsigt3 )
1788	DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1789	!
1790	END SUBROUTINE s_sh94
1791
1792
1793	SUBROUTINE s_sf12
1794	!!----------------------------------------------------------------------
1795	!! * ROUTINE s_sf12 *
1796	!!
1797	!! ** Purpose : stretch the s-coordinate system
1798	!!
1799	!! ** Method : s-coordinate stretch using the Siddorn and Furner 2012?
1800	!! mixed S/sigma/Z coordinate
1801	!!
1802	!! This method allows the maintenance of fixed surface and or
1803	!! bottom cell resolutions (cf. geopotential coordinates)
1804	!! within an analytically derived stretched S-coordinate framework.
1805	!!
1806	!!
1807	!! Reference : Siddorn and Furner 2012 (submitted Ocean modelling).
1808	!!----------------------------------------------------------------------
1809	INTEGER :: ji, jj, jk ! dummy loop argument
1810	REAL(wp) :: zsmth ! smoothing around critical depth
1811	REAL(wp) :: zzs, zzb ! Surface and bottom cell thickness in sigma space
1812	REAL(wp) :: ztmpu, ztmpv, ztmpf
1813	REAL(wp) :: ztmpu1, ztmpv1, ztmpf1
1814	!
1815	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3
1816	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3
1817	!!----------------------------------------------------------------------
1818	!
1819	ALLOCATE( z_gsigw3 (jpi,jpj,jpk), z_gsigt3 (jpi,jpj,jpk) )
1820	ALLOCATE( z_esigt3 (jpi,jpj,jpk), z_esigw3 (jpi,jpj,jpk), z_esigtu3(jpi,jpj,jpk), z_esigtv3(jpi,jpj,jpk))
1821	ALLOCATE( z_esigtf3(jpi,jpj,jpk), z_esigwu3(jpi,jpj,jpk), z_esigwv3(jpi,jpj,jpk) )
1822
1823	z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp
1824	z_esigt3 = 0._wp ; z_esigw3 = 0._wp
1825	z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp
1826	z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp
1827
1828	DO ji = 1, jpi
1829	DO jj = 1, jpj
1830
1831	IF (hbatt(ji,jj)>rn_hc) THEN !deep water, stretched sigma
1832
1833	zzb = hbatt(ji,jj)*rn_zb_a + rn_zb_b ! this forces a linear bottom cell depth relationship with H,.
1834	! could be changed by users but care must be taken to do so carefully
1835	zzb = 1.0_wp-(zzb/hbatt(ji,jj))
1836
1837	zzs = rn_zs / hbatt(ji,jj)
1838
1839	IF (rn_efold /= 0.0_wp) THEN
1840	zsmth = tanh( (hbatt(ji,jj)- rn_hc ) / rn_efold )
1841	ELSE
1842	zsmth = 1.0_wp
1843	ENDIF
1844
1845	DO jk = 1, jpk
1846	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1847	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)
1848	ENDDO
1849	z_gsigw3(ji,jj,:) = fgamma( z_gsigw3(ji,jj,:), zzb, zzs, zsmth )
1850	z_gsigt3(ji,jj,:) = fgamma( z_gsigt3(ji,jj,:), zzb, zzs, zsmth )
1851
1852	ELSE IF (ln_sigcrit) THEN ! shallow water, uniform sigma
1853
1854	DO jk = 1, jpk
1855	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)
1856	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5)/REAL(jpk-1,wp)
1857	END DO
1858
1859	ELSE ! shallow water, z coordinates
1860
1861	DO jk = 1, jpk
1862	z_gsigw3(ji,jj,jk) = REAL(jk-1,wp) /REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1863	z_gsigt3(ji,jj,jk) = (REAL(jk-1,wp)+0.5_wp)/REAL(jpk-1,wp)*(rn_hc/hbatt(ji,jj))
1864	END DO
1865
1866	ENDIF
1867
1868	DO jk = 1, jpkm1
1869	z_esigt3(ji,jj,jk) = z_gsigw3(ji,jj,jk+1) - z_gsigw3(ji,jj,jk)
1870	z_esigw3(ji,jj,jk+1) = z_gsigt3(ji,jj,jk+1) - z_gsigt3(ji,jj,jk)
1871	END DO
1872	z_esigw3(ji,jj,1 ) = 2.0_wp * (z_gsigt3(ji,jj,1 ) - z_gsigw3(ji,jj,1 ))
1873	z_esigt3(ji,jj,jpk) = 2.0_wp * (z_gsigt3(ji,jj,jpk) - z_gsigw3(ji,jj,jpk))
1874
1875	DO jk = 1, jpk
1876	gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk)
1877	gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk)
1878	END DO
1879
1880	ENDDO ! for all jj's
1881	ENDDO ! for all ji's
1882
1883	DO ji=1,jpi-1
1884	DO jj=1,jpj-1
1885
1886	! extend to suit for Wetting/Drying case
1887	ztmpu = hbatt(ji,jj)+hbatt(ji+1,jj)
1888	ztmpv = hbatt(ji,jj)+hbatt(ji,jj+1)
1889	ztmpf = hbatt(ji,jj)+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1)
1890	ztmpu1 = hbatt(ji,jj)*hbatt(ji+1,jj)
1891	ztmpv1 = hbatt(ji,jj)*hbatt(ji,jj+1)
1892	ztmpf1 = MIN(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1)) * &
1893	& MAX(hbatt(ji,jj), hbatt(ji+1,jj), hbatt(ji,jj+1), hbatt(ji+1,jj+1))
1894	DO jk = 1, jpk
1895	z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1896	& / ztmpu
1897	z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1898	& / ztmpu
1899	z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1900	& / ztmpv
1901	z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1902	& / ztmpv
1903
1904	z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
1905	& + hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
1906	& + hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
1907	& + hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) / ztmpf
1908
1909	! Code prior to wetting and drying option (for reference)
1910	!z_esigtu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji+1,jj)z_esigt3(ji+1,jj,jk) ) &
1911	! /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1912	!
1913	!z_esigwu3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji+1,jj)z_esigw3(ji+1,jj,jk) ) &
1914	! /( hbatt(ji,jj)+hbatt(ji+1,jj) )
1915	!
1916	!z_esigtv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigt3(ji,jj,jk)+hbatt(ji,jj+1)z_esigt3(ji,jj+1,jk) ) &
1917	! /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1918	!
1919	!z_esigwv3(ji,jj,jk) = ( hbatt(ji,jj)z_esigw3(ji,jj,jk)+hbatt(ji,jj+1)z_esigw3(ji,jj+1,jk) ) &
1920	! /( hbatt(ji,jj)+hbatt(ji,jj+1) )
1921	!
1922	!z_esigtf3(ji,jj,jk) = ( hbatt(ji ,jj )*z_esigt3(ji ,jj ,jk) &
1923	! & +hbatt(ji+1,jj )*z_esigt3(ji+1,jj ,jk) &
1924	! +hbatt(ji ,jj+1)*z_esigt3(ji ,jj+1,jk) &
1925	! & +hbatt(ji+1,jj+1)*z_esigt3(ji+1,jj+1,jk) ) &
1926	! /( hbatt(ji ,jj )+hbatt(ji+1,jj)+hbatt(ji,jj+1)+hbatt(ji+1,jj+1) )
1927
1928	e3t_0(ji,jj,jk)=(scosrf(ji,jj)+hbatt(ji,jj))*z_esigt3(ji,jj,jk)
1929	e3u_0(ji,jj,jk)=(scosrf(ji,jj)+hbatu(ji,jj))*z_esigtu3(ji,jj,jk)
1930	e3v_0(ji,jj,jk)=(scosrf(ji,jj)+hbatv(ji,jj))*z_esigtv3(ji,jj,jk)
1931	e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk)
1932	!
1933	e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)
1934	e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk)
1935	e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk)
1936	END DO
1937
1938	ENDDO
1939	ENDDO
1940	!
1941	CALL lbc_lnk('domzgr',e3t_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3u_0 ,'T',1.)
1942	CALL lbc_lnk('domzgr',e3v_0 ,'T',1.) ; CALL lbc_lnk('domzgr',e3f_0 ,'T',1.)
1943	CALL lbc_lnk('domzgr',e3w_0 ,'T',1.)
1944	CALL lbc_lnk('domzgr',e3uw_0,'T',1.) ; CALL lbc_lnk('domzgr',e3vw_0,'T',1.)
1945	!
1946	DEALLOCATE( z_gsigw3, z_gsigt3 )
1947	DEALLOCATE( z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )
1948	!
1949	END SUBROUTINE s_sf12
1950
1951
1952	SUBROUTINE s_tanh()
1953	!!----------------------------------------------------------------------
1954	!! * ROUTINE s_tanh*
1955	!!
1956	!! ** Purpose : stretch the s-coordinate system
1957	!!
1958	!! ** Method : s-coordinate stretch
1959	!!
1960	!! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1961	!!----------------------------------------------------------------------
1962	INTEGER :: ji, jj, jk ! dummy loop argument
1963	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
1964	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_gsigw, z_gsigt
1965	REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_esigt, z_esigw
1966	!!----------------------------------------------------------------------
1967
1968	ALLOCATE( z_gsigw(jpk), z_gsigt(jpk) )
1969	ALLOCATE( z_esigt(jpk), z_esigw(jpk) )
1970
1971	z_gsigw = 0._wp ; z_gsigt = 0._wp
1972	z_esigt = 0._wp ; z_esigw = 0._wp
1973
1974	DO jk = 1, jpk
1975	z_gsigw(jk) = -fssig( REAL(jk,wp)-0.5_wp )
1976	z_gsigt(jk) = -fssig( REAL(jk,wp) )
1977	END DO
1978	IF( lwp ) WRITE(numout,*) 'z_gsigw 1 jpk ', z_gsigw(1), z_gsigw(jpk)
1979	!
1980	! Coefficients for vertical scale factors at w-, t- levels
1981	!!gm bug : define it from analytical function, not like juste bellow....
1982	!!gm or betteroffer the 2 possibilities....
1983	DO jk = 1, jpkm1
1984	z_esigt(jk ) = z_gsigw(jk+1) - z_gsigw(jk)
1985	z_esigw(jk+1) = z_gsigt(jk+1) - z_gsigt(jk)
1986	END DO
1987	z_esigw( 1 ) = 2._wp * ( z_gsigt(1 ) - z_gsigw(1 ) )
1988	z_esigt(jpk) = 2._wp * ( z_gsigt(jpk) - z_gsigw(jpk) )
1989	!
1990	DO jk = 1, jpk
1991	zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp)
1992	zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp)
1993	gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigt(jk) + hift(:,:)zcoeft )
1994	gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))z_gsigw(jk) + hift(:,:)zcoefw )
1995	END DO
1996
1997	DO jj = 1, jpj
1998	DO ji = 1, jpi
1999	DO jk = 1, jpk
2000	e3t_0(ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigt(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2001	e3u_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigt(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2002	e3v_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigt(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2003	e3f_0(ji,jj,jk) = ( (hbatf(ji,jj)-hiff(ji,jj))*z_esigt(jk) + hiff(ji,jj)/REAL(jpkm1,wp) )
2004	!
2005	e3w_0 (ji,jj,jk) = ( (hbatt(ji,jj)-hift(ji,jj))*z_esigw(jk) + hift(ji,jj)/REAL(jpkm1,wp) )
2006	e3uw_0(ji,jj,jk) = ( (hbatu(ji,jj)-hifu(ji,jj))*z_esigw(jk) + hifu(ji,jj)/REAL(jpkm1,wp) )
2007	e3vw_0(ji,jj,jk) = ( (hbatv(ji,jj)-hifv(ji,jj))*z_esigw(jk) + hifv(ji,jj)/REAL(jpkm1,wp) )
2008	END DO
2009	END DO
2010	END DO
2011	!
2012	DEALLOCATE( z_gsigw, z_gsigt )
2013	DEALLOCATE( z_esigt, z_esigw )
2014	!
2015	END SUBROUTINE s_tanh
2016
2017
2018	FUNCTION fssig( pk ) RESULT( pf )
2019	!!----------------------------------------------------------------------
2020	!! * ROUTINE fssig *
2021	!!
2022	!! ** Purpose : provide the analytical function in s-coordinate
2023	!!
2024	!! ** Method : the function provide the non-dimensional position of
2025	!! T and W (i.e. between 0 and 1)
2026	!! T-points at integer values (between 1 and jpk)
2027	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2028	!!----------------------------------------------------------------------
2029	REAL(wp), INTENT(in) :: pk ! continuous "k" coordinate
2030	REAL(wp) :: pf ! sigma value
2031	!!----------------------------------------------------------------------
2032	!
2033	pf = ( TANH( rn_theta * ( -(pk-0.5_wp) / REAL(jpkm1) + rn_thetb ) ) &
2034	& - TANH( rn_thetb * rn_theta ) ) &
2035	& * ( COSH( rn_theta ) &
2036	& + COSH( rn_theta * ( 2._wp * rn_thetb - 1._wp ) ) ) &
2037	& / ( 2._wp * SINH( rn_theta ) )
2038	!
2039	END FUNCTION fssig
2040
2041
2042	FUNCTION fssig1( pk1, pbb ) RESULT( pf1 )
2043	!!----------------------------------------------------------------------
2044	!! * ROUTINE fssig1 *
2045	!!
2046	!! ** Purpose : provide the Song and Haidvogel version of the analytical function in s-coordinate
2047	!!
2048	!! ** Method : the function provides the non-dimensional position of
2049	!! T and W (i.e. between 0 and 1)
2050	!! T-points at integer values (between 1 and jpk)
2051	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2052	!!----------------------------------------------------------------------
2053	REAL(wp), INTENT(in) :: pk1 ! continuous "k" coordinate
2054	REAL(wp), INTENT(in) :: pbb ! Stretching coefficient
2055	REAL(wp) :: pf1 ! sigma value
2056	!!----------------------------------------------------------------------
2057	!
2058	IF ( rn_theta == 0 ) then ! uniform sigma
2059	pf1 = - ( pk1 - 0.5_wp ) / REAL( jpkm1 )
2060	ELSE ! stretched sigma
2061	pf1 = ( 1._wp - pbb ) * ( SINH( rn_theta*(-(pk1-0.5_wp)/REAL(jpkm1)) ) ) / SINH( rn_theta ) &
2062	& + pbb * ( (TANH( rn_theta( (-(pk1-0.5_wp)/REAL(jpkm1)) + 0.5_wp) ) - TANH( 0.5_wp rn_theta ) ) &
2063	& / ( 2._wp * TANH( 0.5_wp * rn_theta ) ) )
2064	ENDIF
2065	!
2066	END FUNCTION fssig1
2067
2068
2069	FUNCTION fgamma( pk1, pzb, pzs, psmth) RESULT( p_gamma )
2070	!!----------------------------------------------------------------------
2071	!! * ROUTINE fgamma *
2072	!!
2073	!! ** Purpose : provide analytical function for the s-coordinate
2074	!!
2075	!! ** Method : the function provides the non-dimensional position of
2076	!! T and W (i.e. between 0 and 1)
2077	!! T-points at integer values (between 1 and jpk)
2078	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
2079	!!
2080	!! This method allows the maintenance of fixed surface and or
2081	!! bottom cell resolutions (cf. geopotential coordinates)
2082	!! within an analytically derived stretched S-coordinate framework.
2083	!!
2084	!! Reference : Siddorn and Furner, in prep
2085	!!----------------------------------------------------------------------
2086	REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate
2087	REAL(wp) :: p_gamma(jpk) ! stretched coordinate
2088	REAL(wp), INTENT(in ) :: pzb ! Bottom box depth
2089	REAL(wp), INTENT(in ) :: pzs ! surface box depth
2090	REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter
2091	!
2092	INTEGER :: jk ! dummy loop index
2093	REAL(wp) :: za1,za2,za3 ! local scalar
2094	REAL(wp) :: zn1,zn2 ! - -
2095	REAL(wp) :: za,zb,zx ! - -
2096	!!----------------------------------------------------------------------
2097	!
2098	zn1 = 1._wp / REAL( jpkm1, wp )
2099	zn2 = 1._wp - zn1
2100	!
2101	za1 = (rn_alpha+2.0_wp)zn1(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn1**(rn_alpha+2.0_wp)
2102	za2 = (rn_alpha+2.0_wp)zn2(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)zn2**(rn_alpha+2.0_wp)
2103	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
2104	!
2105	za = pzb - za3*(pzs-za1)-za2
2106	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp) ) )
2107	zb = (pzs - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
2108	zx = 1.0_wp-za/2.0_wp-zb
2109	!
2110	DO jk = 1, jpk
2111	p_gamma(jk) = za(pk1(jk)(1.0_wp-pk1(jk)/2.0_wp))+zbpk1(jk)*3.0_wp + &
2112	& zx( (rn_alpha+2.0_wp)pk1(jk)**(rn_alpha+1.0_wp)- &
2113	& (rn_alpha+1.0_wp)pk1(jk)*(rn_alpha+2.0_wp) )
2114	p_gamma(jk) = p_gamma(jk)psmth+pk1(jk)(1.0_wp-psmth)
2115	END DO
2116	!
2117	END FUNCTION fgamma
2118
2119	!!======================================================================
2120	END MODULE domzgr

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

Download in other formats: