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

Last change on this file since 14698 was 14630, checked in by jchanut, 3 years ago
AGFdomcfg: 1) restore use of analytical grids (i.e. case where jphgr_msh>0). That may be useful to set up test cases with AGRIF. 2) Add random topography over a flat bottom if nn_bathy = -1 (in place of a Gaussian bump). This illustrates well where the interface matching and update are done. #2638
File size: 108.0 KB

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

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

Download in other formats: