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 trunk/NEMO/OPA_SRC/DOM – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/DOM/domzgr.F90 @ 1083

Last change on this file since 1083 was 1083, checked in by ctlod, 16 years ago
trunk: small reorgansition of domzgr.F90 module, see ticket: #191
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 66.3 KB

Line
1	MODULE domzgr
2	!!==============================================================================
3	!! * MODULE domzgr *
4	!! Ocean initialization : domain initialization
5	!!==============================================================================
6
7	!!----------------------------------------------------------------------
8	!! dom_zgr : defined the ocean vertical coordinate system
9	!! zgr_bat : bathymetry fields (levels and meters)
10	!! zgr_bat_zoom : modify the bathymetry field if zoom domain
11	!! zgr_bat_ctl : check the bathymetry files
12	!! zgr_z : reference z-coordinate
13	!! zgr_zco : z-coordinate
14	!! zgr_zps : z-coordinate with partial steps
15	!! zgr_sco : s-coordinate
16	!! fssig : sigma coordinate non-dimensional function
17	!!---------------------------------------------------------------------
18	!! * Modules used
19	USE oce ! ocean dynamics and tracers
20	USE dom_oce ! ocean space and time domain
21	USE in_out_manager ! I/O manager
22	USE lib_mpp ! distributed memory computing library
23	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
24	USE closea
25	USE solisl
26	USE c1d ! 1D configuration
27
28	IMPLICIT NONE
29	PRIVATE
30
31	!! * Routine accessibility
32	PUBLIC dom_zgr ! called by dom_init.F90
33
34	!! * Module variables
35	REAL(wp) :: & !!: Namelist nam_zgr_sco
36	sbot_min = 300. , & !: minimum depth of s-bottom surface (>0) (m)
37	sbot_max = 5250. , & !: maximum depth of s-bottom surface (= ocean depth) (>0) (m)
38	theta = 6.0 , & !: surface control parameter (0<=theta<=20)
39	thetb = 0.75, & !: bottom control parameter (0<=thetb<= 1)
40	r_max = 0.15 !: maximum cut-off r-value allowed (0<r_max<1)
41
42
43
44	!! * Substitutions
45	# include "domzgr_substitute.h90"
46	# include "vectopt_loop_substitute.h90"
47	!!----------------------------------------------------------------------
48	!! OPA 9.0 , LOCEAN-IPSL (2005)
49	!!----------------------------------------------------------------------
50
51	CONTAINS
52
53	SUBROUTINE dom_zgr
54	!!----------------------------------------------------------------------
55	!! * ROUTINE dom_zgr *
56	!!
57	!! ** Purpose : set the depth of model levels and the resulting
58	!! vertical scale factors.
59	!!
60	!! ** Method reference vertical coordinate
61	!!
62	!! ** Action :
63	!!
64	!! History :
65	!! 9.0 ! 03-08 (G. Madec) original code
66	!! 9.0 ! 05-10 (A. Beckmann) modifications for hybrid s-ccordinates
67	!!----------------------------------------------------------------------
68	INTEGER :: ioptio = 0 ! temporary integer
69
70	NAMELIST/nam_zgr/ ln_zco, ln_zps, ln_sco
71	!!----------------------------------------------------------------------
72
73	! Read Namelist nam_zgr : vertical coordinate'
74	! ---------------------
75	REWIND ( numnam )
76	READ ( numnam, nam_zgr )
77
78	! Parameter control and print
79	! ---------------------------
80	! Control print
81	IF(lwp) THEN
82	WRITE(numout,*)
83	WRITE(numout,*) 'dom_zgr : vertical coordinate'
84	WRITE(numout,*) '~~~~~~~'
85	WRITE(numout,*) ' Namelist nam_zgr : set vertical coordinate'
86	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
87	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
88	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
89	ENDIF
90
91	! Check Vertical coordinate options
92	ioptio = 0
93	IF( ln_zco ) ioptio = ioptio + 1
94	IF( ln_zps ) ioptio = ioptio + 1
95	IF( ln_sco ) ioptio = ioptio + 1
96	IF ( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
97
98	IF( lk_zco ) THEN
99	IF(lwp) WRITE(numout,*) ' z-coordinate with reduced incore memory requirement'
100	IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' reduced memory with zps or sco option is impossible' )
101	ENDIF
102
103	! Build the vertical coordinate system
104	! ------------------------------------
105
106	CALL zgr_z ! Reference z-coordinate system (always called)
107
108	CALL zgr_bat ! Bathymetry fields (levels and meters)
109
110	IF( ln_zco ) CALL zgr_zco ! z-coordinate
111
112	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
113
114	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
115
116	!!bug gm control print:
117	IF( nprint == 1 .AND. lwp ) THEN
118	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
119	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
120	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
121	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ), &
122	& ' u ', MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ), &
123	& ' uw', MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)), &
124	& ' w ', MINVAL( fse3w(:,:,:) )
125
126	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
127	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
128	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ), &
129	& ' u ', MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ), &
130	& ' uw', MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)), &
131	& ' w ', MAXVAL( fse3w(:,:,:) )
132	ENDIF
133	!!!bug gm
134
135	END SUBROUTINE dom_zgr
136
137
138	SUBROUTINE zgr_z
139	!!----------------------------------------------------------------------
140	!! * ROUTINE zgr_z *
141	!!
142	!! ** Purpose : set the depth of model levels and the resulting
143	!! vertical scale factors.
144	!!
145	!! ** Method : z-coordinate system (use in all type of coordinate)
146	!! The depth of model levels is defined from an analytical
147	!! function the derivative of which gives the scale factors.
148	!! both depth and scale factors only depend on k (1d arrays).
149	!! w-level: gdepw_0 = fsdep(k)
150	!! e3w_0(k) = dk(fsdep)(k) = fse3(k)
151	!! t-level: gdept_0 = fsdep(k+0.5)
152	!! e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
153	!!
154	!! ** Action : - gdept_0, gdepw_0 : depth of T- and W-point (m)
155	!! - e3t_0, e3w_0 : scale factors at T- and W-levels (m)
156	!!
157	!! Reference :
158	!! Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
159	!!
160	!! History :
161	!! 9.0 ! 03-08 (G. Madec) F90: Free form and module
162	!!----------------------------------------------------------------------
163	!! * Local declarations
164	INTEGER :: jk ! dummy loop indices
165	REAL(wp) :: zt, zw ! temporary scalars
166	REAL(wp) :: &
167	zsur , za0, za1, zkth, zacr, & ! Values set from parameters in
168	zdzmin, zhmax ! par_CONFIG_Rxx.h90
169	!!----------------------------------------------------------------------
170
171	! Set variables from parameters
172	! ------------------------------
173	zkth = ppkth ; zacr = ppacr
174	zdzmin = ppdzmin ; zhmax = pphmax
175
176	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
177	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
178	!
179	IF( ppa1 == pp_to_be_computed .AND. &
180	& ppa0 == pp_to_be_computed .AND. &
181	& ppsur == pp_to_be_computed ) THEN
182	za1 = ( ppdzmin - pphmax / FLOAT(jpk-1) ) &
183	/ ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) &
184	& * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
185	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
186
187	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
188	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
189
190	ELSE
191	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
192	ENDIF
193
194
195	! Parameter print
196	! ---------------
197	IF(lwp) THEN
198	WRITE(numout,*)
199	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
200	WRITE(numout,*) ' ~~~~~~~'
201	IF( ppkth == 0. ) THEN
202	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
203	WRITE(numout,*) ' Total depth :', zhmax
204	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
205	ELSE
206	IF( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0. ) THEN
207	WRITE(numout,*) ' zsur, za0, za1 computed from '
208	WRITE(numout,*) ' zdzmin = ', zdzmin
209	WRITE(numout,*) ' zhmax = ', zhmax
210	ENDIF
211	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
212	WRITE(numout,*) ' zsur = ', zsur
213	WRITE(numout,*) ' za0 = ', za0
214	WRITE(numout,*) ' za1 = ', za1
215	WRITE(numout,*) ' zkth = ', zkth
216	WRITE(numout,*) ' zacr = ', zacr
217	ENDIF
218	ENDIF
219
220
221	! Reference z-coordinate (depth - scale factor at T- and W-points)
222	! ======================
223	IF( ppkth == 0.e0 ) THEN ! uniform vertical grid
224
225	za1 = zhmax / FLOAT(jpk-1)
226	DO jk = 1, jpk
227	zw = FLOAT( jk )
228	zt = FLOAT( jk ) + 0.5
229	gdepw_0(jk) = ( zw - 1 ) * za1
230	gdept_0(jk) = ( zt - 1 ) * za1
231	e3w_0 (jk) = za1
232	e3t_0 (jk) = za1
233	END DO
234
235	ELSE
236
237	DO jk = 1, jpk
238	zw = FLOAT( jk )
239	zt = FLOAT( jk ) + 0.5
240	gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG( COSH( (zw-zkth)/zacr ) ) )
241	gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG( COSH( (zt-zkth)/zacr ) ) )
242	e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth)/zacr )
243	e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth)/zacr )
244	END DO
245	gdepw_0(1) = 0.e0 ! force first w-level to be exactly at zero
246
247	ENDIF
248
249	! Control and print
250	! ==================
251
252	IF(lwp) THEN
253	WRITE(numout,*)
254	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
255	WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" )
256	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
257	ENDIF
258
259	DO jk = 1, jpk
260	IF( e3w_0(jk) <= 0. .OR. e3t_0(jk) <= 0. ) CALL ctl_stop( ' e3w or e3t =< 0 ' )
261	IF( gdepw_0(jk) < 0. .OR. gdept_0(jk) < 0. ) CALL ctl_stop( ' gdepw or gdept < 0 ' )
262	END DO
263
264	END SUBROUTINE zgr_z
265
266
267	SUBROUTINE zgr_bat
268	!!----------------------------------------------------------------------
269	!! * ROUTINE zgr_bat *
270	!!
271	!! ** Purpose : set bathymetry both in levels and meters
272	!!
273	!! ** Method : read or define mbathy and bathy arrays
274	!! * level bathymetry:
275	!! The ocean basin geometry is given by a two-dimensional array,
276	!! mbathy, which is defined as follow :
277	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
278	!! at t-point (ji,jj).
279	!! = 0 over the continental t-point.
280	!! = -n over the nth island t-point.
281	!! The array mbathy is checked to verified its consistency with
282	!! model option. in particular:
283	!! mbathy must have at least 1 land grid-points (mbathy<=0)
284	!! along closed boundary.
285	!! mbathy must be cyclic IF jperio=1.
286	!! mbathy must be lower or equal to jpk-1.
287	!! isolated ocean grid points are suppressed from mbathy
288	!! since they are only connected to remaining
289	!! ocean through vertical diffusion.
290	!! ntopo=-1 : rectangular channel or bassin with a bump
291	!! ntopo= 0 : flat rectangular channel or basin
292	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
293	!! bathy is read in 'bathy_meter.nc' NetCDF file
294	!! C A U T I O N : mbathy will be modified during the initializa-
295	!! tion phase to become the number of non-zero w-levels of a water
296	!! column, with a minimum value of 1.
297	!!
298	!! ** Action : - mbathy: level bathymetry (in level index)
299	!! - bathy : meter bathymetry (in meters)
300	!!
301	!! History :
302	!! 9.0 ! 03-08 (G. Madec) Original code
303	!! 9.0 ! 05-10 (A. Beckmann) modifications for s-ccordinates
304	!!----------------------------------------------------------------------
305	!! * Modules used
306	USE iom
307
308	!! * Local declarations
309	INTEGER :: ji, jj, jl, jk ! dummy loop indices
310	INTEGER :: inum ! temporary logical unit
311	INTEGER :: &
312	ii_bump, ij_bump, ih ! bump center position
313	INTEGER , DIMENSION(jpidta,jpjdta) :: &
314	idta ! global domain integer data
315	REAL(wp) :: &
316	r_bump, h_bump, h_oce, & ! bump characteristics
317	zi, zj, zh ! temporary scalars
318	REAL(wp), DIMENSION(jpidta,jpjdta) :: &
319	zdta ! global domain scalar data
320	!!----------------------------------------------------------------------
321
322	IF(lwp) WRITE(numout,*)
323	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
324	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
325
326	! ========================================
327	! global domain level and meter bathymetry (idta,zdta)
328	! ========================================
329	! ! =============== !
330	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
331	! ! =============== !
332
333	IF( ntopo == 0 ) THEN ! flat basin
334
335	IF(lwp) WRITE(numout,*)
336	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
337
338	idta(:,:) = jpkm1 ! flat basin
339	zdta(:,:) = gdepw_0(jpk)
340	h_oce = gdepw_0(jpk)
341
342	ELSE ! bump
343	IF(lwp) WRITE(numout,*)
344	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
345
346	ii_bump = jpidta / 2 ! i-index of the bump center
347	ij_bump = jpjdta / 2 ! j-index of the bump center
348	r_bump = 50000.e0 ! bump radius (meters)
349	h_bump = 2700.e0 ! bump height (meters)
350	h_oce = gdepw_0(jpk) ! background ocean depth (meters)
351	IF(lwp) WRITE(numout,*) ' bump characteristics: '
352	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
353	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
354	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
355	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
356	! zdta :
357	DO jj = 1, jpjdta
358	DO ji = 1, jpidta
359	zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump
360	zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump
361	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
362	END DO
363	END DO
364
365	! idta :
366	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
367	idta(:,:) = jpkm1
368	ELSE ! z-coordinate (zco or zps): step-like topography
369	idta(:,:) = jpkm1
370	DO jk = 1, jpkm1
371	DO jj = 1, jpjdta
372	DO ji = 1, jpidta
373	IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) ) idta(ji,jj) = jk
374	END DO
375	END DO
376	END DO
377	ENDIF
378	ENDIF
379
380	! set boundary conditions (caution, idta on the global domain: use of jperio, not nperio)
381	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
382	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1.e0
383	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0
384	ELSEIF( jperio == 2 ) THEN
385	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
386	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0
387	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0.e0
388	idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0.e0
389	ELSE
390	ih = 0 ; zh = 0.e0
391	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
392	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
393	idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh
394	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
395	idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh
396	ENDIF
397
398	! =======================================
399	! local domain level and meter bathymetry (mbathy,bathy)
400	! =======================================
401
402	mbathy(:,:) = 0 ! set to zero extra halo points
403	bathy (:,:) = 0.e0 ! (require for mpp case)
404
405	DO jj = 1, nlcj ! interior values
406	DO ji = 1, nlci
407	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
408	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
409	END DO
410	END DO
411
412	! ! =============== !
413	ELSEIF( ntopo == 1 ) THEN ! read in file !
414	! ! =============== !
415
416	IF( ln_zco ) THEN
417	CALL iom_open ( 'bathy_level.nc', inum ) ! Level bathymetry
418	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
419	CALL iom_close (inum)
420	mbathy(:,:) = INT( bathy(:,:) )
421	ENDIF
422
423	IF( ln_zps .OR. ln_sco ) THEN
424	CALL iom_open ( 'bathy_meter.nc', inum ) ! meter bathymetry
425	CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
426	CALL iom_close (inum)
427	ENDIF
428	! ! =============== !
429	ELSE ! error !
430	! ! =============== !
431	WRITE(ctmp1,*) ' parameter , ntopo = ', ntopo
432	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
433	ENDIF
434
435	! =======================
436	! NO closed seas or lakes
437	! =======================
438
439	IF( nclosea == 0 ) THEN
440	DO jl = 1, jpncs
441	DO jj = ncsj1(jl), ncsj2(jl)
442	DO ji = ncsi1(jl), ncsi2(jl)
443	mbathy(ji,jj) = 0 ! suppress closed seas
444	bathy (ji,jj) = 0.e0 ! and lakes
445	!!bug gm bathy (ji,jj) = 300.e0 ! and lakes
446	END DO
447	END DO
448	END DO
449	ENDIF
450
451	#if defined key_orca_lev10
452	! 10 time the vertical resolution
453	mbathy(:,:) = 10 * mbathy(:,:)
454	IF(lwp) WRITE(numout,*) ' ATTENTION: 300 niveaux avec bathy levels "vraie?"'
455	#endif
456	! ===========
457	! Zoom domain
458	! ===========
459
460	IF( lzoom ) CALL zgr_bat_zoom
461
462	! ================
463	! Bathymetry check
464	! ================
465
466	CALL zgr_bat_ctl
467
468	END SUBROUTINE zgr_bat
469
470
471	SUBROUTINE zgr_bat_zoom
472	!!----------------------------------------------------------------------
473	!! * ROUTINE zgr_bat_zoom *
474	!!
475	!! ** Purpose : - Close zoom domain boundary if necessary
476	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
477	!!
478	!! ** Method :
479	!!
480	!! ** Action : - update mbathy: level bathymetry (in level index)
481	!!
482	!! History :
483	!! 9.0 ! 03-08 (G. Madec) Original code
484	!!----------------------------------------------------------------------
485	!! * Local variables
486	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
487	!!----------------------------------------------------------------------
488
489	IF(lwp) WRITE(numout,*)
490	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
491	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
492
493	! Zoom domain
494	! ===========
495
496	! Forced closed boundary if required
497	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
498	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
499	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
500	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
501
502	! Configuration specific domain modifications
503	! (here, ORCA arctic configuration: suppress Med Sea)
504	IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
505	SELECT CASE ( jp_cfg )
506	! ! =======================
507	CASE ( 2 ) ! ORCA_R2 configuration
508	! ! =======================
509	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
510
511	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
512	ij0 = 98 ; ij1 = 110
513	! ! =======================
514	CASE ( 05 ) ! ORCA_R05 configuration
515	! ! =======================
516	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
517
518	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
519	ij0 = 314 ; ij1 = 370
520
521	END SELECT
522	!
523	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
524	!
525	ENDIF
526
527	END SUBROUTINE zgr_bat_zoom
528
529
530	SUBROUTINE zgr_bat_ctl
531	!!----------------------------------------------------------------------
532	!! * ROUTINE zgr_bat_ctl *
533	!!
534	!! ** Purpose : check the bathymetry in levels
535	!!
536	!! ** Method : The array mbathy is checked to verified its consistency
537	!! with the model options. in particular:
538	!! mbathy must have at least 1 land grid-points (mbathy<=0)
539	!! along closed boundary.
540	!! mbathy must be cyclic IF jperio=1.
541	!! mbathy must be lower or equal to jpk-1.
542	!! isolated ocean grid points are suppressed from mbathy
543	!! since they are only connected to remaining
544	!! ocean through vertical diffusion.
545	!! C A U T I O N : mbathy will be modified during the initializa-
546	!! tion phase to become the number of non-zero w-levels of a water
547	!! column, with a minimum value of 1.
548	!!
549	!! ** Action : - update mbathy: level bathymetry (in level index)
550	!! - update bathy : meter bathymetry (in meters)
551	!!
552	!! History :
553	!! 9.0 ! 03-08 (G. Madec) Original code
554	!! 9.0 ! 05-10 (A. Beckmann) modifications for s-ccordinates
555	!!----------------------------------------------------------------------
556	!! * Local declarations
557	INTEGER :: ji, jj, jl ! dummy loop indices
558	INTEGER :: &
559	icompt, ibtest, ikmax ! temporary integers
560	REAL(wp), DIMENSION(jpi,jpj) :: &
561	zbathy ! temporary workspace
562	!!----------------------------------------------------------------------
563
564	IF(lwp) WRITE(numout,*)
565	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
566	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
567
568	! ================
569	! Bathymetry check
570	! ================
571
572	IF( .NOT. lk_c1d ) THEN
573
574	! Suppress isolated ocean grid points
575
576	IF(lwp) WRITE(numout,*)
577	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
578	IF(lwp) WRITE(numout,*)' -----------------------------------'
579
580	icompt = 0
581
582	DO jl = 1, 2
583
584	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
585	mbathy( 1 ,:) = mbathy(jpim1,:)
586	mbathy(jpi,:) = mbathy( 2 ,:)
587	ENDIF
588	DO jj = 2, jpjm1
589	DO ji = 2, jpim1
590	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
591	& mbathy(ji,jj-1),mbathy(ji,jj+1) )
592	IF( ibtest < mbathy(ji,jj) ) THEN
593	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
594	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
595	mbathy(ji,jj) = ibtest
596	icompt = icompt + 1
597	ENDIF
598	END DO
599	END DO
600
601	END DO
602	IF( icompt == 0 ) THEN
603	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
604	ELSE
605	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
606	ENDIF
607	IF( lk_mpp ) THEN
608	zbathy(:,:) = FLOAT( mbathy(:,:) )
609	CALL lbc_lnk( zbathy, 'T', 1. )
610	mbathy(:,:) = INT( zbathy(:,:) )
611	ENDIF
612
613	! 3.2 East-west cyclic boundary conditions
614
615	IF( nperio == 0 ) THEN
616	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west', &
617	' boundary: nperio = ', nperio
618	IF( lk_mpp ) THEN
619	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
620	IF( jperio /= 1 ) mbathy(1,:) = 0
621	ENDIF
622	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
623	IF( jperio /= 1 ) mbathy(nlci,:) = 0
624	ENDIF
625	ELSE
626	IF( ln_zco .OR. ln_zps ) THEN
627	mbathy( 1 ,:) = 0
628	mbathy(jpi,:) = 0
629	ELSE
630	mbathy( 1 ,:) = jpkm1
631	mbathy(jpi,:) = jpkm1
632	ENDIF
633	ENDIF
634	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
635	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions', &
636	' on mbathy: nperio = ', nperio
637	mbathy( 1 ,:) = mbathy(jpim1,:)
638	mbathy(jpi,:) = mbathy( 2 ,:)
639	ELSEIF( nperio == 2 ) THEN
640	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions', &
641	' on mbathy: nperio = ', nperio
642	ELSE
643	IF(lwp) WRITE(numout,*) ' e r r o r'
644	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
645	! STOP 'dom_mba'
646	ENDIF
647
648	! Set to zero mbathy over islands if necessary (lk_isl=F)
649	IF( .NOT. lk_isl ) THEN ! No island
650	IF(lwp) WRITE(numout,*)
651	IF(lwp) WRITE(numout,*) ' mbathy set to 0 over islands'
652	IF(lwp) WRITE(numout,*) ' ----------------------------'
653
654	mbathy(:,:) = MAX( 0, mbathy(:,:) )
655
656	! Boundary condition on mbathy
657	IF( .NOT.lk_mpp ) THEN
658	!!bug ??? y reflechir!
659	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
660	zbathy(:,:) = FLOAT( mbathy(:,:) )
661	CALL lbc_lnk( zbathy, 'T', 1. )
662	mbathy(:,:) = INT( zbathy(:,:) )
663	ENDIF
664
665	ENDIF
666
667	ENDIF
668
669	! Number of ocean level inferior or equal to jpkm1
670
671	ikmax = 0
672	DO jj = 1, jpj
673	DO ji = 1, jpi
674	ikmax = MAX( ikmax, mbathy(ji,jj) )
675	END DO
676	END DO
677	!!! test a faire: ikmax = MAX( mbathy(:,:) ) ???
678
679	IF( ikmax > jpkm1 ) THEN
680	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
681	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
682	ELSE IF( ikmax < jpkm1 ) THEN
683	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
684	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
685	ENDIF
686
687	IF( lwp .AND. nprint == 1 ) THEN
688	WRITE(numout,*)
689	WRITE(numout,*) ' bathymetric field '
690	WRITE(numout,*) ' ------------------'
691	WRITE(numout,*) ' number of non-zero T-levels '
692	CALL prihin( mbathy, jpi, jpj, 1, jpi, &
693	1 , 1 , jpj, 1, 3 , &
694	numout )
695	WRITE(numout,*)
696	ENDIF
697
698	END SUBROUTINE zgr_bat_ctl
699
700
701	SUBROUTINE zgr_zco
702	!!----------------------------------------------------------------------
703	!! * ROUTINE zgr_zco *
704	!!
705	!! ** Purpose : define the z-coordinate system
706	!!
707	!! ** Method : set 3D coord. arrays to reference 1D array if lk_zco=T
708	!!
709	!! History :
710	!! ! 06-04 (R. Benshila, G. Madec) Original code
711	!!----------------------------------------------------------------------
712	INTEGER :: jk
713	!!----------------------------------------------------------------------
714
715	IF( .NOT.lk_zco ) THEN
716	DO jk = 1, jpk
717	fsdept(:,:,jk) = gdept_0(jk)
718	fsdepw(:,:,jk) = gdepw_0(jk)
719	fsde3w(:,:,jk) = gdepw_0(jk)
720	fse3t (:,:,jk) = e3t_0(jk)
721	fse3u (:,:,jk) = e3t_0(jk)
722	fse3v (:,:,jk) = e3t_0(jk)
723	fse3f (:,:,jk) = e3t_0(jk)
724	fse3w (:,:,jk) = e3w_0(jk)
725	fse3uw(:,:,jk) = e3w_0(jk)
726	fse3vw(:,:,jk) = e3w_0(jk)
727	END DO
728	ENDIF
729
730	END SUBROUTINE zgr_zco
731
732	#if defined key_zco
733	!!----------------------------------------------------------------------
734	!! 'key_zco' : "pure" zco (gdep & e3 are 1D arrays)
735	!!----------------------------------------------------------------------
736	SUBROUTINE zgr_zps ! Empty routine
737	END SUBROUTINE zgr_zps
738	SUBROUTINE zgr_sco ! Empty routine
739	END SUBROUTINE zgr_sco
740
741	#else
742	!!----------------------------------------------------------------------
743	!! Default option : zco, zps and/or sco available (gedp & e3 are 3D arrays)
744	!!----------------------------------------------------------------------
745
746	SUBROUTINE zgr_zps
747	!!----------------------------------------------------------------------
748	!! * ROUTINE zgr_zps *
749	!!
750	!! ** Purpose : the depth and vertical scale factor in partial step
751	!! z-coordinate case
752	!!
753	!! ** Method : Partial steps : computes the 3D vertical scale factors
754	!! of T-, U-, V-, W-, UW-, VW and F-points that are associated with
755	!! a partial step representation of bottom topography.
756	!!
757	!! The reference depth of model levels is defined from an analytical
758	!! function the derivative of which gives the reference vertical
759	!! scale factors.
760	!! From depth and scale factors reference, we compute there new value
761	!! with partial steps on 3d arrays ( i, j, k ).
762	!!
763	!! w-level: gdepw(i,j,k) = fsdep(k)
764	!! e3w(i,j,k) = dk(fsdep)(k) = fse3(i,j,k)
765	!! t-level: gdept(i,j,k) = fsdep(k+0.5)
766	!! e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5)
767	!!
768	!! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc),
769	!! we find the mbathy index of the depth at each grid point.
770	!! This leads us to three cases:
771	!!
772	!! - bathy = 0 => mbathy = 0
773	!! - 1 < mbathy < jpkm1
774	!! - bathy > gdepw(jpk) => mbathy = jpkm1
775	!!
776	!! Then, for each case, we find the new depth at t- and w- levels
777	!! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw-
778	!! and f-points.
779	!!
780	!! This routine is given as an example, it must be modified
781	!! following the user s desiderata. nevertheless, the output as
782	!! well as the way to compute the model levels and scale factors
783	!! must be respected in order to insure second order accuracy
784	!! schemes.
785	!!
786	!! c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives
787	!! - - - - - - - gdept, gdepw and e3. are positives
788	!!
789	!! Reference :
790	!! Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270.
791	!!
792	!! History :
793	!! 8.5 ! 02-09 (A. Bozec, G. Madec) F90: Free form and module
794	!! 9.0 ! 02-09 (A. de Miranda) rigid-lid + islands
795	!!----------------------------------------------------------------------
796	!! * Local declarations
797	INTEGER :: ji, jj, jk ! dummy loop indices
798	INTEGER :: &
799	ik, it ! temporary integers
800
801	REAL(wp) :: &
802	ze3tp, ze3wp, & ! Last ocean level thickness at T- and W-points
803	zdepwp, & ! Ajusted ocean depth to avoid too small e3t
804	zdepth, & ! " "
805	zmax, zmin, & ! Maximum and minimum depth
806	zdiff ! temporary scalar
807
808	REAL(wp), DIMENSION(jpi,jpj) :: &
809	zprt ! " "
810
811	LOGICAL :: ll_print ! Allow control print for debugging
812
813	!!---------------------------------------------------------------------
814	!! OPA8.5, LODYC-IPSL (2002)
815	!!---------------------------------------------------------------------
816
817	! Local variable for debugging
818	ll_print=.FALSE.
819	!!! ll_print=.TRUE.
820
821	! Initialization of constant
822	zmax = gdepw_0(jpk) + e3t_0(jpk)
823	zmin = gdepw_0(4)
824
825	! Ocean depth
826	IF(lwp .AND. ll_print) THEN
827	WRITE(numout,*)
828	WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'
829	CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )
830	ENDIF
831
832	IF(lwp) WRITE(numout,*)
833	IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps'
834	IF(lwp) WRITE(numout,*) ' ~~~~~~~ '
835	IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used'
836
837
838	! bathymetry in level (from bathy_meter)
839	! ===================
840
841	! initialize mbathy to the maximum ocean level available
842	mbathy(:,:) = jpkm1
843
844	! storage of land and island's number (zera and negative values) in mbathy
845	DO jj = 1, jpj
846	DO ji= 1, jpi
847	IF( bathy(ji,jj) <= 0. ) mbathy(ji,jj) = INT( bathy(ji,jj) )
848	END DO
849	END DO
850
851	! bounded value of bathy
852	! minimum depth == 3 levels
853	! maximum depth == gdepw_0(jpk)+e3t_0(jpk)
854	! i.e. the last ocean level thickness cannot exceed e3t_0(jpkm1)+e3t_0(jpk)
855	DO jj = 1, jpj
856	DO ji= 1, jpi
857	IF( bathy(ji,jj) <= 0. ) THEN
858	bathy(ji,jj) = 0.e0
859	ELSE
860	bathy(ji,jj) = MAX( bathy(ji,jj), zmin )
861	bathy(ji,jj) = MIN( bathy(ji,jj), zmax )
862	ENDIF
863	END DO
864	END DO
865
866	! Compute mbathy for ocean points (i.e. the number of ocean levels)
867	! find the number of ocean levels such that the last level thickness
868	! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where
869	! e3t_0 is the reference level thickness
870	DO jk = jpkm1, 1, -1
871	zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat )
872	DO jj = 1, jpj
873	DO ji = 1, jpi
874	IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth ) mbathy(ji,jj) = jk-1
875	END DO
876	END DO
877	END DO
878
879	! Scale factors and depth at T- and W-points
880
881	! intitialization to the reference z-coordinate
882	DO jk = 1, jpk
883	gdept(:,:,jk) = gdept_0(jk)
884	gdepw(:,:,jk) = gdepw_0(jk)
885	e3t (:,:,jk) = e3t_0 (jk)
886	e3w (:,:,jk) = e3w_0 (jk)
887	END DO
888	hdept(:,:) = gdept(:,:,2 )
889	hdepw(:,:) = gdepw(:,:,3 )
890
891	!
892	DO jj = 1, jpj
893	DO ji = 1, jpi
894	ik = mbathy(ji,jj)
895	! ocean point only
896	IF( ik > 0 ) THEN
897	! max ocean level case
898	IF( ik == jpkm1 ) THEN
899	zdepwp = bathy(ji,jj)
900	ze3tp = bathy(ji,jj) - gdepw_0(ik)
901	ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) )
902	e3t(ji,jj,ik ) = ze3tp
903	e3t(ji,jj,ik+1) = ze3tp
904	e3w(ji,jj,ik ) = ze3wp
905	e3w(ji,jj,ik+1) = ze3tp
906	gdepw(ji,jj,ik+1) = zdepwp
907	gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp
908	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp
909	! standard case
910	ELSE
911	!!alex
912	IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN
913	gdepw(ji,jj,ik+1) = bathy(ji,jj)
914	ELSE
915	!!alex ctl write(,) 'zps',ji,jj,'bathy', bathy(ji,jj), 'depw ',gdepw_0(ik+1)
916	gdepw(ji,jj,ik+1) = gdepw_0(ik+1)
917	ENDIF
918	!!Alex
919	!!Alex gdepw(ji,jj,ik+1) = bathy(ji,jj)
920	!!bug gm verifier les gdepw_0
921	gdept(ji,jj,ik ) = gdepw_0(ik) + ( gdepw(ji,jj,ik+1) - gdepw_0(ik) ) &
922	& * ((gdept_0( ik ) - gdepw_0(ik) ) &
923	& / ( gdepw_0( ik+1) - gdepw_0(ik) ))
924	e3t(ji,jj,ik) = e3t_0(ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) &
925	& / ( gdepw_0( ik+1) - gdepw_0(ik) )
926	e3w(ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) &
927	& * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) )
928	! ... on ik+1
929	e3w(ji,jj,ik+1) = e3t(ji,jj,ik)
930	e3t(ji,jj,ik+1) = e3t(ji,jj,ik)
931	gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik)
932	ENDIF
933	ENDIF
934	END DO
935	END DO
936
937	it = 0
938	DO jj = 1, jpj
939	DO ji = 1, jpi
940	ik = mbathy(ji,jj)
941	! ocean point only
942	IF( ik > 0 ) THEN
943	! bathymetry output
944	hdept(ji,jj) = gdept(ji,jj,ik )
945	hdepw(ji,jj) = gdepw(ji,jj,ik+1)
946	e3tp (ji,jj) = e3t(ji,jj,ik )
947	e3wp (ji,jj) = e3w(ji,jj,ik )
948	! test
949	zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik )
950	IF( zdiff <= 0. .AND. lwp ) THEN
951	it=it+1
952	WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj
953	WRITE(numout,*) ' bathy = ', bathy(ji,jj)
954	WRITE(numout,*) ' gdept= ', gdept(ji,jj,ik), ' gdepw= ', gdepw(ji,jj,ik+1), &
955	' zdiff = ', zdiff
956	WRITE(numout,*) ' e3tp = ', e3t(ji,jj,ik ), ' e3wp = ', e3w(ji,jj,ik )
957	ENDIF
958	ENDIF
959	END DO
960	END DO
961
962	! Scale factors and depth at U-, V-, UW and VW-points
963
964	! initialisation to z-scale factors
965	DO jk = 1, jpk
966	e3u (:,:,jk) = e3t_0(jk)
967	e3v (:,:,jk) = e3t_0(jk)
968	e3uw(:,:,jk) = e3w_0(jk)
969	e3vw(:,:,jk) = e3w_0(jk)
970	END DO
971
972	! Computed as the minimum of neighbooring scale factors
973	DO jk = 1,jpk
974	DO jj = 1, jpjm1
975	DO ji = 1, fs_jpim1 ! vector opt.
976	e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk))
977	e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk))
978	e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) )
979	e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) )
980	END DO
981	END DO
982	END DO
983
984	! Boundary conditions
985	CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. )
986	CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. )
987
988	! set to z-scale factor if zero (i.e. along closed boundaries)
989	DO jk = 1, jpk
990	DO jj = 1, jpj
991	DO ji = 1, jpi
992	IF( e3u (ji,jj,jk) == 0.e0 ) e3u (ji,jj,jk) = e3t_0(jk)
993	IF( e3v (ji,jj,jk) == 0.e0 ) e3v (ji,jj,jk) = e3t_0(jk)
994	IF( e3uw(ji,jj,jk) == 0.e0 ) e3uw(ji,jj,jk) = e3w_0(jk)
995	IF( e3vw(ji,jj,jk) == 0.e0 ) e3vw(ji,jj,jk) = e3w_0(jk)
996	END DO
997	END DO
998	END DO
999
1000	! Scale factor at F-point
1001
1002	! initialisation to z-scale factors
1003	DO jk = 1, jpk
1004	e3f (:,:,jk) = e3t_0(jk)
1005	END DO
1006
1007	! Computed as the minimum of neighbooring V-scale factors
1008	DO jk = 1, jpk
1009	DO jj = 1, jpjm1
1010	DO ji = 1, fs_jpim1 ! vector opt.
1011	e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) )
1012	END DO
1013	END DO
1014	END DO
1015	! Boundary conditions
1016	CALL lbc_lnk( e3f, 'F', 1. )
1017
1018	! set to z-scale factor if zero (i.e. along closed boundaries)
1019	DO jk = 1, jpk
1020	DO jj = 1, jpj
1021	DO ji = 1, jpi
1022	IF( e3f(ji,jj,jk) == 0.e0 ) e3f(ji,jj,jk) = e3t_0(jk)
1023	END DO
1024	END DO
1025	END DO
1026	!!bug ? gm: must be a do loop with mj0,mj1
1027
1028	! we duplicate factor scales for jj = 1 and jj = 2
1029	e3t(:,mj0(1),:) = e3t(:,mj0(2),:)
1030	e3w(:,mj0(1),:) = e3w(:,mj0(2),:)
1031	e3u(:,mj0(1),:) = e3u(:,mj0(2),:)
1032	e3v(:,mj0(1),:) = e3v(:,mj0(2),:)
1033	e3f(:,mj0(1),:) = e3f(:,mj0(2),:)
1034
1035
1036
1037	! Compute gdep3w (vertical sum of e3w)
1038
1039	gdep3w (:,:,1) = 0.5 * e3w (:,:,1)
1040	DO jk = 2, jpk
1041	gdep3w (:,:,jk) = gdep3w (:,:,jk-1) + e3w (:,:,jk)
1042	END DO
1043
1044	! Control print
1045	9600 FORMAT(9x,' level gdept gdepw e3t e3w ')
1046	9610 FORMAT(10x,i4,4f9.2)
1047	IF(lwp .AND. ll_print) THEN
1048	DO jj = 1,jpj
1049	DO ji = 1, jpi
1050	ik = MAX(mbathy(ji,jj),1)
1051	zprt(ji,jj) = e3t(ji,jj,ik)
1052	END DO
1053	END DO
1054	WRITE(numout,*)
1055	WRITE(numout,*) 'domzgr e3t(mbathy)'
1056	CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1057	DO jj = 1,jpj
1058	DO ji = 1, jpi
1059	ik = MAX(mbathy(ji,jj),1)
1060	zprt(ji,jj) = e3w(ji,jj,ik)
1061	END DO
1062	END DO
1063	WRITE(numout,*)
1064	WRITE(numout,*) 'domzgr e3w(mbathy)'
1065	CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1066	DO jj = 1,jpj
1067	DO ji = 1, jpi
1068	ik = MAX(mbathy(ji,jj),1)
1069	zprt(ji,jj) = e3u(ji,jj,ik)
1070	END DO
1071	END DO
1072
1073	WRITE(numout,*)
1074	WRITE(numout,*) 'domzgr e3u(mbathy)'
1075	CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1076	DO jj = 1,jpj
1077	DO ji = 1, jpi
1078	ik = MAX(mbathy(ji,jj),1)
1079	zprt(ji,jj) = e3v(ji,jj,ik)
1080	END DO
1081	END DO
1082	WRITE(numout,*)
1083	WRITE(numout,*) 'domzgr e3v(mbathy)'
1084	CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1085	DO jj = 1,jpj
1086	DO ji = 1, jpi
1087	ik = MAX(mbathy(ji,jj),1)
1088	zprt(ji,jj) = e3f(ji,jj,ik)
1089	END DO
1090	END DO
1091
1092	WRITE(numout,*)
1093	WRITE(numout,*) 'domzgr e3f(mbathy)'
1094	CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1095	DO jj = 1,jpj
1096	DO ji = 1, jpi
1097	ik = MAX(mbathy(ji,jj),1)
1098	zprt(ji,jj) = gdep3w(ji,jj,ik)
1099	END DO
1100	END DO
1101	WRITE(numout,*)
1102	WRITE(numout,*) 'domzgr gdep3w(mbathy)'
1103	CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
1104
1105	ENDIF
1106
1107
1108	DO jk = 1,jpk
1109	DO jj = 1, jpj
1110	DO ji = 1, jpi
1111	IF( e3w(ji,jj,jk) <= 0. .or. e3t(ji,jj,jk) <= 0. ) THEN
1112	IF(lwp) THEN
1113	WRITE(numout,*) ' e r r o r : e3w or e3t =< 0 '
1114	WRITE(numout,*) ' ========= --------------- '
1115	WRITE(numout,*)
1116	ENDIF
1117	STOP 'domzgr.psteps'
1118	ENDIF
1119	IF( gdepw(ji,jj,jk) < 0. .or. gdept(ji,jj,jk) < 0. ) THEN
1120	IF(lwp) THEN
1121	WRITE(numout,*) ' e r r o r : gdepw or gdept < 0 '
1122	WRITE(numout,*) ' ========= ------------------ '
1123	WRITE(numout,*)
1124	ENDIF
1125	STOP 'domzgr.psteps'
1126	ENDIF
1127	END DO
1128	END DO
1129	END DO
1130
1131	IF(lwp) THEN
1132	WRITE(numout,*) ' e3t lev 21 '
1133	CALL prihre(e3t(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
1134	WRITE(numout,*) ' e3w lev 21 '
1135	CALL prihre(e3w(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
1136	WRITE(numout,*) ' e3u lev 21 '
1137	CALL prihre(e3u(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
1138	WRITE(numout,*) ' e3v lev 21 '
1139	CALL prihre(e3v(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
1140	WRITE(numout,*) ' e3f lev 21 '
1141	CALL prihre(e3f(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout)
1142	WRITE(numout,*) ' e3t lev 22 '
1143	CALL prihre(e3t(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
1144	WRITE(numout,*) ' e3w lev 22 '
1145	CALL prihre(e3w(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
1146	WRITE(numout,*) ' e3u lev 22 '
1147	CALL prihre(e3u(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
1148	WRITE(numout,*) ' e3v lev 22 '
1149	CALL prihre(e3v(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
1150	WRITE(numout,*) ' e3f lev 22 '
1151	CALL prihre(e3f(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout)
1152	ENDIF
1153
1154	! ===========
1155	! Zoom domain
1156	! ===========
1157
1158	IF( lzoom ) CALL zgr_bat_zoom
1159
1160	! ================
1161	! Bathymetry check
1162	! ================
1163
1164	CALL zgr_bat_ctl
1165
1166	END SUBROUTINE zgr_zps
1167
1168
1169	FUNCTION fssig( pk ) RESULT( pf )
1170	!!----------------------------------------------------------------------
1171	!! * ROUTINE eos_init *
1172	!!
1173	!! ** Purpose : provide the analytical function in s-coordinate
1174	!!
1175	!! ** Method : the function provide the non-dimensional position of
1176	!! T and W (i.e. between 0 and 1)
1177	!! T-points at integer values (between 1 and jpk)
1178	!! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
1179	!!
1180	!! Reference : ???
1181	!!----------------------------------------------------------------------
1182	REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate
1183	REAL(wp) :: pf ! sigma value
1184	!!----------------------------------------------------------------------
1185	!
1186	pf = ( TANH( theta * ( -(pk-0.5) / REAL(jpkm1) + thetb ) ) &
1187	& - TANH( thetb * theta ) ) &
1188	& * ( COSH( theta ) &
1189	& + COSH( theta * ( 2.e0 * thetb - 1.e0 ) ) ) &
1190	& / ( 2.e0 * SINH( theta ) )
1191	!
1192	END FUNCTION fssig
1193
1194
1195	SUBROUTINE zgr_sco
1196	!!----------------------------------------------------------------------
1197	!! * ROUTINE zgr_sco *
1198	!!
1199	!! ** Purpose : define the s-coordinate system
1200	!!
1201	!! ** Method : s-coordinate
1202	!! The depth of model levels is defined as the product of an
1203	!! analytical function by the local bathymetry, while the vertical
1204	!! scale factors are defined as the product of the first derivative
1205	!! of the analytical function by the bathymetry.
1206	!! (this solution save memory as depth and scale factors are not
1207	!! 3d fields)
1208	!! - Read bathymetry (in meters) at t-point and compute the
1209	!! bathymetry at u-, v-, and f-points.
1210	!! hbatu = mi( hbatt )
1211	!! hbatv = mj( hbatt )
1212	!! hbatf = mi( mj( hbatt ) )
1213	!! - Compute gsigt, gsigw, esigt, esigw from an analytical
1214	!! function and its derivative given as function.
1215	!! gsigt(k) = fssig (k )
1216	!! gsigw(k) = fssig (k-0.5)
1217	!! esigt(k) = fsdsig(k )
1218	!! esigw(k) = fsdsig(k-0.5)
1219	!! This routine is given as an example, it must be modified
1220	!! following the user s desiderata. nevertheless, the output as
1221	!! well as the way to compute the model levels and scale factors
1222	!! must be respected in order to insure second order a!!uracy
1223	!! schemes.
1224	!!
1225	!! Reference :
1226	!! Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
1227	!!
1228	!! History :
1229	!! ! 95-12 (G. Madec) Original code : s vertical coordinate
1230	!! ! 97-07 (G. Madec) lbc_lnk call
1231	!! ! 97-04 (J.-O. Beismann)
1232	!! 8.5 ! 02-06 (G. Madec) F90: Free form and module
1233	!! 9.0 ! 05-10 (A. Beckmann) new stretching function
1234	!!----------------------------------------------------------------------
1235	!! * Local declarations
1236	INTEGER :: ji, jj, jk, jl
1237	INTEGER :: iip1, ijp1, iim1, ijm1
1238	REAL(wp) :: zrmax, ztaper
1239
1240	REAL(wp), DIMENSION(jpi,jpj) :: &
1241	zenv, ztmp, zmsk, zri, zrj, zhbat
1242
1243	NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max
1244	!!----------------------------------------------------------------------
1245
1246	! Read Namelist nam_zgr_sco : sigma-stretching parameters
1247	! -------------------------
1248	REWIND( numnam )
1249	READ ( numnam, nam_zgr_sco )
1250
1251	IF(lwp) THEN
1252	WRITE(numout,*)
1253	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
1254	WRITE(numout,*) '~~~~~~~~~~~'
1255	WRITE(numout,*) ' Namelist nam_zgr_sco'
1256	WRITE(numout,*) ' sigma-stretching coeffs '
1257	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) sbot_max = ' ,sbot_max
1258	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) sbot_min = ' ,sbot_min
1259	WRITE(numout,*) ' surface control parameter (0<=theta<=20) theta = ', theta
1260	WRITE(numout,*) ' bottom control parameter (0<=thetb<= 1) thetb = ', thetb
1261	WRITE(numout,*) ' maximum cut-off r-value allowed r_max = ' , r_max
1262	ENDIF
1263
1264	! ???
1265	! ------------------
1266	hift(:,:) = sbot_min
1267	hifu(:,:) = sbot_min
1268	hifv(:,:) = sbot_min
1269	hiff(:,:) = sbot_min
1270
1271
1272	! set maximum ocean depth
1273	! -----------------------
1274	DO jj = 1, jpj
1275	DO ji = 1, jpi
1276	bathy(ji,jj) = MIN( sbot_max, bathy(ji,jj) )
1277	END DO
1278	END DO
1279
1280	! Define the envelop bathymetry
1281	! =============================
1282	! Smooth the bathymetry (if required)
1283
1284	scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea)
1285	scobot(:,:) = bathy(:,:) ! ocean bottom depth
1286
1287
1288	! use r-value to create hybrid coordinates
1289
1290	DO jj = 1, jpj
1291	DO ji = 1, jpi
1292	zenv(ji,jj) = MAX( bathy(ji,jj), sbot_min )
1293	END DO
1294	END DO
1295
1296	jl = 0
1297	zrmax = 1.e0
1298	! ! ================ !
1299	DO WHILE ( jl <= 10000 .AND. zrmax > r_max ) ! Iterative loop !
1300	! ! ================ !
1301	jl = jl + 1
1302
1303	zrmax = 0.e0
1304	zmsk(:,:) = 0.e0
1305
1306	DO jj = 1, nlcj
1307	DO ji = 1, nlci
1308	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
1309	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
1310	zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
1311	zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
1312	zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
1313	IF( zri(ji,jj) > r_max ) zmsk(ji ,jj ) = 1.0
1314	IF( zri(ji,jj) > r_max ) zmsk(iip1,jj ) = 1.0
1315	IF( zrj(ji,jj) > r_max ) zmsk(ji ,jj ) = 1.0
1316	IF( zrj(ji,jj) > r_max ) zmsk(ji ,ijp1) = 1.0
1317	END DO
1318	END DO
1319	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
1320	! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
1321	ztmp(:,:) = zmsk(:,:)
1322	CALL lbc_lnk( zmsk, 'T', 1. )
1323	DO jj = 1, nlcj
1324	DO ji = 1, nlci
1325	zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
1326	END DO
1327	END DO
1328
1329	WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
1330
1331	DO jj = 1, nlcj
1332	DO ji = 1, nlci
1333	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
1334	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
1335	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
1336	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
1337	IF( zmsk(ji,jj) == 1.0 ) THEN
1338	ztmp(ji,jj) = ( &
1339	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
1340	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2.e0 + zenv(iip1,jj )*zmsk(iip1,jj ) &
1341	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
1342	& ) / ( &
1343	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
1344	& + zmsk(iim1,jj ) + 2.e0 + zmsk(iip1,jj ) &
1345	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
1346	& )
1347	ENDIF
1348	END DO
1349	END DO
1350
1351	DO jj = 1, nlcj
1352	DO ji = 1, nlci
1353	IF( zmsk(ji,jj) == 1.0 ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
1354	END DO
1355	END DO
1356
1357	! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero
1358	ztmp(:,:) = zenv(:,:)
1359	CALL lbc_lnk( zenv, 'T', 1. )
1360	DO jj = 1, nlcj
1361	DO ji = 1, nlci
1362	IF( zenv(ji,jj) == 0.e0 ) zenv(ji,jj) = ztmp(ji,jj)
1363	END DO
1364	END DO
1365	! ! ================ !
1366	END DO ! End loop !
1367	! ! ================ !
1368
1369	! save envelop bathymetry in hbatt
1370	hbatt(:,:) = zenv(:,:)
1371
1372	!!bug gm : CAUTION: tapering hard coded !!!! orca2 only
1373	!!bug gm NOT valid in mpp ===> taper have been changed
1374
1375	DO jj = 1, jpj
1376	DO ji = 1, jpi
1377	!!bug mpp ztaper = EXP( -(FLOAT(jj-74)/10.)**2 )
1378	ztaper = EXP( -(gphit(ji,jj)/8)**2 )
1379	hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
1380	END DO
1381	END DO
1382
1383	DO jj = 1, jpj
1384	zrmax = EXP( -(gphit(10,jj)/8)**2 )
1385	ztaper = EXP( -(FLOAT(jj-74)/10.)**2 )
1386	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) jj, (FLOAT(jj-74)/10.), ztaper,(gphit(10,jj)/8), zrmax
1387	END DO
1388
1389	IF( nprint == 1 .AND. lwp ) THEN
1390	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
1391	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
1392	ENDIF
1393
1394	! Control print
1395	IF(lwp) THEN
1396	WRITE(numout,*)
1397	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
1398	WRITE(numout,*)
1399	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
1400	ENDIF
1401
1402	! 1. computation of hbatu, hbatv, hbatf fields
1403	! --------------------------------------------
1404
1405	hbatu(:,:) = sbot_min
1406	hbatv(:,:) = sbot_min
1407	hbatf(:,:) = sbot_min
1408
1409	IF(lwp) THEN
1410	WRITE(numout,*)
1411	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min
1412	WRITE(numout,*)
1413	ENDIF
1414
1415	DO jj = 1, jpjm1
1416	DO ji = 1, jpim1
1417	hbatu(ji,jj)= 0.5 *( hbatt(ji ,jj)+hbatt(ji+1,jj ) )
1418	hbatv(ji,jj)= 0.5 *( hbatt(ji ,jj)+hbatt(ji ,jj+1) )
1419	hbatf(ji,jj)= 0.25*( hbatt(ji ,jj)+hbatt(ji ,jj+1) &
1420	+hbatt(ji+1,jj)+hbatt(ji+1,jj+1) )
1421	END DO
1422	END DO
1423
1424	! Apply lateral boundary condition
1425	! CAUTION: retain non zero value in the initial file
1426	! this should be OK for orca cfg, not for EEL
1427	zhbat(:,:) = hbatu(:,:)
1428	CALL lbc_lnk( hbatu, 'U', 1. )
1429	DO jj = 1, jpj
1430	DO ji = 1, jpi
1431	IF( hbatu(ji,jj) == 0.e0 ) THEN
1432	IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = sbot_min
1433	IF( zhbat(ji,jj) /= 0.e0 ) hbatu(ji,jj) = zhbat(ji,jj)
1434	ENDIF
1435	END DO
1436	END DO
1437
1438	zhbat(:,:) = hbatv(:,:)
1439	CALL lbc_lnk( hbatv, 'V', 1. )
1440	DO jj = 1, jpj
1441	DO ji = 1, jpi
1442	IF( hbatv(ji,jj) == 0.e0 ) THEN
1443	IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = sbot_min
1444	IF( zhbat(ji,jj) /= 0.e0 ) hbatv(ji,jj) = zhbat(ji,jj)
1445	ENDIF
1446	END DO
1447	END DO
1448
1449	zhbat(:,:) = hbatf(:,:)
1450	CALL lbc_lnk( hbatf, 'F', 1. )
1451	DO jj = 1, jpj
1452	DO ji = 1, jpi
1453	IF( hbatf(ji,jj) == 0.e0 ) THEN
1454	IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = sbot_min
1455	IF( zhbat(ji,jj) /= 0.e0 ) hbatf(ji,jj) = zhbat(ji,jj)
1456	ENDIF
1457	END DO
1458	END DO
1459
1460
1461	!!bug: key_helsinki a verifer
1462	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1463	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1464	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1465	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1466
1467	IF( nprint == 1 .AND. lwp ) THEN
1468	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift(:,:) ), ' f ', MAXVAL( hiff(:,:) ), &
1469	& ' u ', MAXVAL( hifu(:,:) ), ' v ', MAXVAL( hifv(:,:) )
1470	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift(:,:) ), ' f ', MINVAL( hiff(:,:) ), &
1471	& ' u ', MINVAL( hifu(:,:) ), ' v ', MINVAL( hifv(:,:) )
1472	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1473	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1474	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1475	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1476	ENDIF
1477	!! helsinki
1478
1479	! 2. Computation of gsig and esig fields
1480	! --------------------------------------
1481
1482	! 2.1 Coefficients for model level depth at w- and t-levels
1483
1484	DO jk = 1, jpk
1485	gsigw(jk) = -fssig( FLOAT(jk)-0.5 )
1486	gsigt(jk) = -fssig( FLOAT(jk) )
1487	END DO
1488
1489	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)
1490
1491	! 2.2 Coefficients for vertical scale factors at w-, t- levels
1492
1493	!!bug gm : define it from analytical function, not like juste bellow....
1494	!! or betteroffer the 2 possibilities....
1495	DO jk=2,jpk
1496	esigw(jk)=gsigt(jk)-gsigt(jk-1)
1497	END DO
1498	DO jk=1,jpkm1
1499	esigt(jk)=gsigw(jk+1)-gsigw(jk)
1500	END DO
1501	esigw(1)=esigw(2)
1502	esigt(jpk)=esigt(jpkm1)
1503
1504	! 2.3 Coefficients for vertical depth as the sum of e3w scale factors
1505
1506	gsi3w(1) = 0.5 * esigw(1)
1507	DO jk = 2, jpk
1508	gsi3w(jk) = gsi3w(jk-1)+ esigw(jk)
1509	END DO
1510
1511	!!bug gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1512	DO jk = 1, jpk
1513	!! change the solver stat....
1514	!! zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
1515	!! gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigt(jk)+hift(:,:)zcoeft)
1516	gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigt(jk)+hift(:,:)(FLOAT(jk)-0.5)/FLOAT(jpkm1))
1517	!! zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
1518	!! gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigw(jk)+hift(:,:)zcoefw)
1519	!! gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsi3w(jk)+hift(:,:)zcoefw)
1520	gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigw(jk)+hift(:,:)FLOAT(jk-1)/FLOAT(jpkm1))
1521	gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsi3w(jk)+hift(:,:)FLOAT(jk-1)/FLOAT(jpkm1))
1522	END DO
1523
1524	!!bug gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1525	DO jj = 1, jpj
1526	DO ji = 1, jpi
1527	DO jk = 1, jpk
1528	e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
1529	e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1530	e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1531	e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
1532
1533	e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
1534	e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1535	e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1536	END DO
1537	END DO
1538	END DO
1539
1540	! HYBRID
1541	DO jj = 1, jpj
1542	DO ji = 1, jpi
1543	DO jk = 1, jpkm1
1544	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1545	IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0
1546	END DO
1547	END DO
1548	END DO
1549
1550	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1551
1552
1553	! ===========
1554	! Zoom domain
1555	! ===========
1556
1557	IF( lzoom ) CALL zgr_bat_zoom
1558
1559	! 2.4 Control print
1560
1561	IF(lwp) THEN
1562	WRITE(numout,*)
1563	WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1564	WRITE(numout,9400)
1565	WRITE(numout,9410) (jk,gsigt(jk),gsigw(jk),esigt(jk),esigw(jk),gsi3w(jk),jk=1,jpk)
1566	ENDIF
1567	9400 FORMAT(9x,' level gsigt gsigw esigt esigw gsi3w')
1568	9410 FORMAT(10x,i4,5f11.4)
1569
1570	IF(lwp) THEN
1571	WRITE(numout,*)
1572	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1573	WRITE(numout,*) ' ~~~~~~ --------------------'
1574	WRITE(numout,9420)
1575	WRITE(numout,9430) (jk,fsdept(1,1,jk),fsdepw(1,1,jk), &
1576	fse3t (1,1,jk),fse3w (1,1,jk),jk=1,jpk)
1577	DO jj = mj0(20), mj1(20)
1578	DO ji = mi0(20), mi1(20)
1579	WRITE(numout,*)
1580	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1581	WRITE(numout,*) ' ~~~~~~ --------------------'
1582	WRITE(numout,9420)
1583	WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk), &
1584	& fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk)
1585	END DO
1586	END DO
1587	DO jj = mj0(74), mj1(74)
1588	DO ji = mi0(100), mi1(100)
1589	WRITE(numout,*)
1590	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1591	WRITE(numout,*) ' ~~~~~~ --------------------'
1592	WRITE(numout,9420)
1593	WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk), &
1594	& fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk)
1595	END DO
1596	END DO
1597	ENDIF
1598
1599	9420 FORMAT(9x,' level gdept gdepw gde3w e3t e3w ')
1600	9430 FORMAT(10x,i4,4f9.2)
1601
1602	!!bug gm no more necessary? if ! defined key_helsinki
1603	DO jk = 1, jpk
1604	DO jj = 1, jpj
1605	DO ji = 1, jpi
1606	IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
1607	IF(lwp) THEN
1608	WRITE(numout,*)
1609	WRITE(numout,*) ' e r r o r : e3w or e3t =< 0 '
1610	WRITE(numout,*) ' ========= --------------- '
1611	WRITE(numout,*)
1612	WRITE(numout,*) ' point (i,j,k)= ',ji,jj,jk
1613	WRITE(numout,*)
1614	ENDIF
1615	STOP 'domzgr'
1616	ENDIF
1617	IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
1618	IF(lwp) THEN
1619	WRITE(numout,*)
1620	WRITE(numout,*) ' e r r o r : gdepw or gdept < 0 '
1621	WRITE(numout,*) ' ========= ------------------ '
1622	WRITE(numout,*)
1623	WRITE(numout,*) ' point (i,j,k)= ', ji, jj, jk
1624	WRITE(numout,*)
1625	ENDIF
1626	STOP 'domzgr'
1627	ENDIF
1628	END DO
1629	END DO
1630	END DO
1631	!!bug gm #endif
1632
1633	END SUBROUTINE zgr_sco
1634
1635	#endif
1636
1637
1638	!!======================================================================
1639	END MODULE domzgr

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

Download in other formats: