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

Last change on this file since 584 was 528, checked in by opalod, 18 years ago
RB:nemo_v1_bugfix_063: change ln_zco to lk_zco
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 51.0 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	!!---------------------------------------------------------------------
17	!! * Modules used
18	USE oce ! ocean dynamics and tracers
19	USE dom_oce ! ocean space and time domain
20	USE in_out_manager ! I/O manager
21	USE lib_mpp ! distributed memory computing library
22	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
23	USE closea
24	USE solisl
25	USE ini1d ! initialization of the 1D configuration
26
27	IMPLICIT NONE
28	PRIVATE
29
30	!! * Routine accessibility
31	PUBLIC dom_zgr ! called by dom_init.F90
32
33	!! * Module variables
34	REAL(wp) :: & !!: Namelist nam_zgr_sco
35	sbot_min = 300. , & !: minimum depth of s-bottom surface (>0) (m)
36	sbot_max = 5250. , & !: maximum depth of s-bottom surface (= ocean depth) (>0) (m)
37	theta = 6.0 , & !: surface control parameter (0<=theta<=20)
38	thetb = 0.75, & !: bottom control parameter (0<=thetb<= 1)
39	r_max = 0.15 !: maximum cut-off r-value allowed (0<r_max<1)
40
41
42
43	!! * Substitutions
44	# include "domzgr_substitute.h90"
45	# include "vectopt_loop_substitute.h90"
46	!!----------------------------------------------------------------------
47	!! OPA 9.0 , LOCEAN-IPSL (2005)
48	!!----------------------------------------------------------------------
49
50	CONTAINS
51
52	SUBROUTINE dom_zgr
53	!!----------------------------------------------------------------------
54	!! * ROUTINE dom_zgr *
55	!!
56	!! ** Purpose : set the depth of model levels and the resulting
57	!! vertical scale factors.
58	!!
59	!! ** Method reference vertical coordinate
60	!!
61	!! ** Action :
62	!!
63	!! History :
64	!! 9.0 ! 03-08 (G. Madec) original code
65	!! 9.0 ! 05-10 (A. Beckmann) modifications for hybrid s-ccordinates
66	!!----------------------------------------------------------------------
67	INTEGER :: ioptio = 0 ! temporary integer
68
69	NAMELIST/nam_zgr/ ln_zco, ln_zps, ln_sco
70	!!----------------------------------------------------------------------
71
72	! Read Namelist nam_zgr : vertical coordinate'
73	! ---------------------
74	REWIND ( numnam )
75	READ ( numnam, nam_zgr )
76
77	! Parameter control and print
78	! ---------------------------
79	! Control print
80	IF(lwp) THEN
81	WRITE(numout,*)
82	WRITE(numout,*) 'dom_zgr : vertical coordinate'
83	WRITE(numout,*) '~~~~~~~'
84	WRITE(numout,*) ' Namelist nam_zgr : set vertical coordinate'
85	WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco
86	WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps
87	WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco
88	ENDIF
89
90	! Check Vertical coordinate options
91	ioptio = 0
92	IF( ln_zco ) ioptio = ioptio + 1
93	IF( ln_zps ) ioptio = ioptio + 1
94	IF( ln_sco ) ioptio = ioptio + 1
95	IF ( ioptio /= 1 ) CALL ctl_stop( ' none or several vertical coordinate options used' )
96
97	IF( lk_zco ) THEN
98	IF(lwp) WRITE(numout,*) ' z-coordinate with reduced incore memory requirement'
99	IF( ln_zps .OR. ln_sco ) CALL ctl_stop( ' reduced memory with zps or sco option is impossible' )
100	ENDIF
101
102	! Build the vertical coordinate system
103	! ------------------------------------
104
105	CALL zgr_z ! Reference z-coordinate system (always called)
106
107	CALL zgr_bat ! Bathymetry fields (levels and meters)
108
109	IF( ln_zco ) CALL zgr_zco ! z-coordinate
110
111	IF( ln_zps ) CALL zgr_zps ! Partial step z-coordinate
112
113	IF( ln_sco ) CALL zgr_sco ! s-coordinate or hybrid z-s coordinate
114
115	!!bug gm control print:
116	IF( nprint == 1 .AND. lwp ) THEN
117	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
118	WRITE(numout,*) ' MIN val depth t ', MINVAL( fsdept(:,:,:) ), &
119	& ' w ', MINVAL( fsdepw(:,:,:) ), '3w ', MINVAL( fsde3w(:,:,:) )
120	WRITE(numout,*) ' MIN val e3 t ', MINVAL( fse3t(:,:,:) ), ' f ', MINVAL( fse3f(:,:,:) ), &
121	& ' u ', MINVAL( fse3u(:,:,:) ), ' u ', MINVAL( fse3v(:,:,:) ), &
122	& ' uw', MINVAL( fse3uw(:,:,:)), ' vw', MINVAL( fse3vw(:,:,:)), &
123	& ' w ', MINVAL( fse3w(:,:,:) )
124
125	WRITE(numout,*) ' MAX val depth t ', MAXVAL( fsdept(:,:,:) ), &
126	& ' w ', MAXVAL( fsdepw(:,:,:) ), '3w ', MAXVAL( fsde3w(:,:,:) )
127	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( fse3t(:,:,:) ), ' f ', MAXVAL( fse3f(:,:,:) ), &
128	& ' u ', MAXVAL( fse3u(:,:,:) ), ' u ', MAXVAL( fse3v(:,:,:) ), &
129	& ' uw', MAXVAL( fse3uw(:,:,:)), ' vw', MAXVAL( fse3vw(:,:,:)), &
130	& ' w ', MAXVAL( fse3w(:,:,:) )
131	ENDIF
132	!!!bug gm
133
134	END SUBROUTINE dom_zgr
135
136
137	SUBROUTINE zgr_z
138	!!----------------------------------------------------------------------
139	!! * ROUTINE zgr_z *
140	!!
141	!! ** Purpose : set the depth of model levels and the resulting
142	!! vertical scale factors.
143	!!
144	!! ** Method : z-coordinate system (use in all type of coordinate)
145	!! The depth of model levels is defined from an analytical
146	!! function the derivative of which gives the scale factors.
147	!! both depth and scale factors only depend on k (1d arrays).
148	!! w-level: gdepw_0 = fsdep(k)
149	!! e3w_0(k) = dk(fsdep)(k) = fse3(k)
150	!! t-level: gdept_0 = fsdep(k+0.5)
151	!! e3t_0(k) = dk(fsdep)(k+0.5) = fse3(k+0.5)
152	!!
153	!! ** Action : - gdept_0, gdepw_0 : depth of T- and W-point (m)
154	!! - e3t_0, e3w_0 : scale factors at T- and W-levels (m)
155	!!
156	!! Reference :
157	!! Marti, Madec & Delecluse, 1992, JGR, 97, No8, 12,763-12,766.
158	!!
159	!! History :
160	!! 9.0 ! 03-08 (G. Madec) F90: Free form and module
161	!!----------------------------------------------------------------------
162	!! * Local declarations
163	INTEGER :: jk ! dummy loop indices
164	REAL(wp) :: zt, zw ! temporary scalars
165	REAL(wp) :: &
166	zsur , za0, za1, zkth, zacr, & ! Values set from parameters in
167	zdzmin, zhmax ! par_CONFIG_Rxx.h90
168	!!----------------------------------------------------------------------
169
170	! Set variables from parameters
171	! ------------------------------
172	zkth = ppkth ; zacr = ppacr
173	zdzmin = ppdzmin ; zhmax = pphmax
174
175	! If ppa1 and ppa0 and ppsur are et to pp_to_be_computed
176	! za0, za1, zsur are computed from ppdzmin , pphmax, ppkth, ppacr
177	!
178	IF( ppa1 == pp_to_be_computed .AND. &
179	& ppa0 == pp_to_be_computed .AND. &
180	& ppsur == pp_to_be_computed ) THEN
181	za1 = ( ppdzmin - pphmax / FLOAT(jpk-1) ) &
182	/ ( TANH((1-ppkth)/ppacr) - ppacr/FLOAT(jpk-1) &
183	& * ( LOG( COSH( (jpk - ppkth) / ppacr) ) &
184	& - LOG( COSH( ( 1 - ppkth) / ppacr) ) ) )
185
186	za0 = ppdzmin - za1 * TANH( (1-ppkth) / ppacr )
187	zsur = - za0 - za1 * ppacr * LOG( COSH( (1-ppkth) / ppacr ) )
188
189	ELSE
190	za1 = ppa1 ; za0 = ppa0 ; zsur = ppsur
191	ENDIF
192
193
194	! Parameter print
195	! ---------------
196	IF(lwp) THEN
197	WRITE(numout,*)
198	WRITE(numout,*) ' zgr_z : Reference vertical z-coordinates'
199	WRITE(numout,*) ' ~~~~~~~'
200	IF( ppkth == 0. ) THEN
201	WRITE(numout,*) ' Uniform grid with ',jpk-1,' layers'
202	WRITE(numout,*) ' Total depth :', zhmax
203	WRITE(numout,*) ' Layer thickness:', zhmax/(jpk-1)
204	ELSE
205	IF( ppa1 == 0. .AND. ppa0 == 0. .AND. ppsur == 0. ) THEN
206	WRITE(numout,*) ' zsur, za0, za1 computed from '
207	WRITE(numout,*) ' zdzmin = ', zdzmin
208	WRITE(numout,*) ' zhmax = ', zhmax
209	ENDIF
210	WRITE(numout,*) ' Value of coefficients for vertical mesh:'
211	WRITE(numout,*) ' zsur = ', zsur
212	WRITE(numout,*) ' za0 = ', za0
213	WRITE(numout,*) ' za1 = ', za1
214	WRITE(numout,*) ' zkth = ', zkth
215	WRITE(numout,*) ' zacr = ', zacr
216	ENDIF
217	ENDIF
218
219
220	! Reference z-coordinate (depth - scale factor at T- and W-points)
221	! ======================
222	IF( ppkth == 0.e0 ) THEN ! uniform vertical grid
223
224	za1 = zhmax / FLOAT(jpk-1)
225	DO jk = 1, jpk
226	zw = FLOAT( jk )
227	zt = FLOAT( jk ) + 0.5
228	gdepw_0(jk) = ( zw - 1 ) * za1
229	gdept_0(jk) = ( zt - 1 ) * za1
230	e3w_0 (jk) = za1
231	e3t_0 (jk) = za1
232	END DO
233
234	ELSE
235
236	DO jk = 1, jpk
237	zw = FLOAT( jk )
238	zt = FLOAT( jk ) + 0.5
239	gdepw_0(jk) = ( zsur + za0 * zw + za1 * zacr * LOG( COSH( (zw-zkth)/zacr ) ) )
240	gdept_0(jk) = ( zsur + za0 * zt + za1 * zacr * LOG( COSH( (zt-zkth)/zacr ) ) )
241	e3w_0 (jk) = za0 + za1 * TANH( (zw-zkth)/zacr )
242	e3t_0 (jk) = za0 + za1 * TANH( (zt-zkth)/zacr )
243	END DO
244	gdepw_0(1) = 0.e0 ! force first w-level to be exactly at zero
245
246	ENDIF
247
248	! Control and print
249	! ==================
250
251	IF(lwp) THEN
252	WRITE(numout,*)
253	WRITE(numout,*) ' Reference z-coordinate depth and scale factors:'
254	WRITE(numout, "(9x,' level gdept gdepw e3t e3w ')" )
255	WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_0(jk), gdepw_0(jk), e3t_0(jk), e3w_0(jk), jk = 1, jpk )
256	ENDIF
257
258	DO jk = 1, jpk
259	IF( e3w_0(jk) <= 0. .OR. e3t_0(jk) <= 0. ) CALL ctl_stop( ' e3w or e3t =< 0 ' )
260	IF( gdepw_0(jk) < 0. .OR. gdept_0(jk) < 0. ) CALL ctl_stop( ' gdepw or gdept < 0 ' )
261	END DO
262
263	END SUBROUTINE zgr_z
264
265
266	SUBROUTINE zgr_bat
267	!!----------------------------------------------------------------------
268	!! * ROUTINE zgr_bat *
269	!!
270	!! ** Purpose : set bathymetry both in levels and meters
271	!!
272	!! ** Method : read or define mbathy and bathy arrays
273	!! * level bathymetry:
274	!! The ocean basin geometry is given by a two-dimensional array,
275	!! mbathy, which is defined as follow :
276	!! mbathy(ji,jj) = 1, ..., jpk-1, the number of ocean level
277	!! at t-point (ji,jj).
278	!! = 0 over the continental t-point.
279	!! = -n over the nth island t-point.
280	!! The array mbathy is checked to verified its consistency with
281	!! model option. in particular:
282	!! mbathy must have at least 1 land grid-points (mbathy<=0)
283	!! along closed boundary.
284	!! mbathy must be cyclic IF jperio=1.
285	!! mbathy must be lower or equal to jpk-1.
286	!! isolated ocean grid points are suppressed from mbathy
287	!! since they are only connected to remaining
288	!! ocean through vertical diffusion.
289	!! ntopo=-1 : rectangular channel or bassin with a bump
290	!! ntopo= 0 : flat rectangular channel or basin
291	!! ntopo= 1 : mbathy is read in 'bathy_level.nc' NetCDF file
292	!! bathy is read in 'bathy_meter.nc' NetCDF file
293	!! C A U T I O N : mbathy will be modified during the initializa-
294	!! tion phase to become the number of non-zero w-levels of a water
295	!! column, with a minimum value of 1.
296	!!
297	!! ** Action : - mbathy: level bathymetry (in level index)
298	!! - bathy : meter bathymetry (in meters)
299	!!
300	!! History :
301	!! 9.0 ! 03-08 (G. Madec) Original code
302	!! 9.0 ! 05-10 (A. Beckmann) modifications for s-ccordinates
303	!!----------------------------------------------------------------------
304	!! * Modules used
305	USE iom
306
307	!! * Local declarations
308	INTEGER :: ji, jj, jl, jk ! dummy loop indices
309	INTEGER :: inum ! temporary logical unit
310	INTEGER :: &
311	ii_bump, ij_bump, ih ! bump center position
312	INTEGER , DIMENSION(jpidta,jpjdta) :: &
313	idta ! global domain integer data
314	REAL(wp) :: &
315	r_bump, h_bump, h_oce, & ! bump characteristics
316	zi, zj, zh ! temporary scalars
317	REAL(wp), DIMENSION(jpidta,jpjdta) :: &
318	zdta ! global domain scalar data
319	!!----------------------------------------------------------------------
320
321	IF(lwp) WRITE(numout,*)
322	IF(lwp) WRITE(numout,*) ' zgr_bat : defines level and meter bathymetry'
323	IF(lwp) WRITE(numout,*) ' ~~~~~~~'
324
325	! ========================================
326	! global domain level and meter bathymetry (idta,zdta)
327	! ========================================
328	! ! =============== !
329	IF( ntopo == 0 .OR. ntopo == -1 ) THEN ! defined by hand !
330	! ! =============== !
331
332	IF( ntopo == 0 ) THEN ! flat basin
333
334	IF(lwp) WRITE(numout,*)
335	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin'
336
337	idta(:,:) = jpkm1 ! flat basin
338	zdta(:,:) = gdepw_0(jpk)
339
340	ELSE ! bump
341	IF(lwp) WRITE(numout,*)
342	IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump'
343
344	ii_bump = jpidta / 3 + 3 ! i-index of the bump center
345	ij_bump = jpjdta / 2 ! j-index of the bump center
346	r_bump = 6 ! bump radius (index)
347	h_bump = 240.e0 ! bump height (meters)
348	h_oce = gdepw_0(jpk) ! background ocean depth (meters)
349	IF(lwp) WRITE(numout,*) ' bump characteristics: '
350	IF(lwp) WRITE(numout,*) ' bump center (i,j) = ', ii_bump, ii_bump
351	IF(lwp) WRITE(numout,*) ' bump height = ', h_bump , ' meters'
352	IF(lwp) WRITE(numout,*) ' bump radius = ', r_bump , ' index'
353	IF(lwp) WRITE(numout,*) ' background ocean depth = ', h_oce , ' meters'
354	! zdta :
355	DO jj = 1, jpjdta
356	DO ji = 1, jpidta
357	zi = FLOAT( ji - ii_bump ) / r_bump
358	zj = FLOAT( jj - ij_bump ) / r_bump
359	zdta(ji,jj) = h_oce - h_bump * EXP( -( zizi + zjzj ) )
360	END DO
361	END DO
362
363	! idta :
364	IF( ln_sco ) THEN ! s-coordinate (zsc ): idta()=jpk
365	idta(:,:) = jpkm1
366	ELSE ! z-coordinate (zco or zps): step-like topography
367	idta(:,:) = jpkm1
368	DO jk = 1, jpkm1
369	DO jj = 1, jpjdta
370	DO ji = 1, jpidta
371	IF( gdept_0(jk) < zdta(ji,jj) .AND. zdta(ji,jj) <= gdept_0(jk+1) ) idta(ji,jj) = jk
372	END DO
373	END DO
374	END DO
375	ENDIF
376	ENDIF
377
378	! set boundary conditions (caution, idta on the global domain: use of jperio, not nperio)
379	IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 ) THEN
380	idta( : , 1 ) = -1 ; zdta( : , 1 ) = -1.e0
381	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0
382	ELSEIF( jperio == 2 ) THEN
383	idta( : , 1 ) = idta( : , 3 ) ; zdta( : , 1 ) = zdta( : , 3 )
384	idta( : ,jpjdta) = 0 ; zdta( : ,jpjdta) = 0.e0
385	idta( 1 , : ) = 0 ; zdta( 1 , : ) = 0.e0
386	idta(jpidta, : ) = 0 ; zdta(jpidta, : ) = 0.e0
387	ELSE
388	ih = 0 ; zh = 0.e0
389	IF( ln_sco ) ih = jpkm1 ; IF( ln_sco ) zh = h_oce
390	idta( : , 1 ) = ih ; zdta( : , 1 ) = zh
391	idta( : ,jpjdta) = ih ; zdta( : ,jpjdta) = zh
392	idta( 1 , : ) = ih ; zdta( 1 , : ) = zh
393	idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh
394	ENDIF
395
396	! EEL R5 configuration with east and west open boundaries.
397	! Two rows of zeroes are needed at the south and north for OBCs
398
399	IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN
400	ih = 0 ; zh = 0.e0
401	IF( ln_sco ) zh = jpkm1 ; IF( ln_sco ) zh = h_oce
402	idta( : , 2 ) = jpkm1 ; zdta( : , 2 ) = h_oce
403	idta( : ,jpjdta-1) = jpkm1 ; zdta( : ,jpjdta-1) = h_oce
404	ENDIF
405
406	! =======================================
407	! local domain level and meter bathymetry (mbathy,bathy)
408	! =======================================
409
410	mbathy(:,:) = 0 ! set to zero extra halo points
411	bathy (:,:) = 0.e0 ! (require for mpp case)
412
413	DO jj = 1, nlcj ! interior values
414	DO ji = 1, nlci
415	mbathy(ji,jj) = idta( mig(ji), mjg(jj) )
416	bathy (ji,jj) = zdta( mig(ji), mjg(jj) )
417	END DO
418	END DO
419
420	! ! =============== !
421	ELSEIF( ntopo == 1 ) THEN ! read in file !
422	! ! =============== !
423
424	IF( ln_zco ) THEN
425	CALL iom_open ( 'bathy_level.nc', inum ) ! Level bathymetry
426	CALL iom_get ( inum, jpdom_data, 'Bathy_level', bathy )
427	CALL iom_close (inum)
428	mbathy(:,:) = INT( bathy(:,:) )
429	ENDIF
430
431	IF( ln_zps .OR. ln_sco ) THEN
432	CALL iom_open ( 'bathy_meter.nc', inum ) ! meter bathymetry
433	CALL iom_get ( inum, jpdom_data, 'Bathymetry', bathy )
434	CALL iom_close (inum)
435	ENDIF
436	! ! =============== !
437	ELSE ! error !
438	! ! =============== !
439	WRITE(ctmp1,*) ' parameter , ntopo = ', ntopo
440	CALL ctl_stop( ' zgr_bat : '//trim(ctmp1) )
441	ENDIF
442
443	! =======================
444	! NO closed seas or lakes
445	! =======================
446
447	IF( nclosea == 0 ) THEN
448	DO jl = 1, jpncs
449	DO jj = ncsj1(jl), ncsj2(jl)
450	DO ji = ncsi1(jl), ncsi2(jl)
451	mbathy(ji,jj) = 0 ! suppress closed seas
452	bathy (ji,jj) = 0.e0 ! and lakes
453	!!bug gm bathy (ji,jj) = 300.e0 ! and lakes
454	END DO
455	END DO
456	END DO
457	ENDIF
458
459	#if defined key_orca_lev10
460	! 10 time the vertical resolution
461	mbathy(:,:) = 10 * mbathy(:,:)
462	IF(lwp) WRITE(numout,*) ' ATTENTION: 300 niveaux avec bathy levels "vraie?"'
463	#endif
464	! ===========
465	! Zoom domain
466	! ===========
467
468	IF( lzoom ) CALL zgr_bat_zoom
469
470	! ================
471	! Bathymetry check
472	! ================
473
474	CALL zgr_bat_ctl
475
476	END SUBROUTINE zgr_bat
477
478
479	SUBROUTINE zgr_bat_zoom
480	!!----------------------------------------------------------------------
481	!! * ROUTINE zgr_bat_zoom *
482	!!
483	!! ** Purpose : - Close zoom domain boundary if necessary
484	!! - Suppress Med Sea from ORCA R2 and R05 arctic zoom
485	!!
486	!! ** Method :
487	!!
488	!! ** Action : - update mbathy: level bathymetry (in level index)
489	!!
490	!! History :
491	!! 9.0 ! 03-08 (G. Madec) Original code
492	!!----------------------------------------------------------------------
493	!! * Local variables
494	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
495	!!----------------------------------------------------------------------
496
497	IF(lwp) WRITE(numout,*)
498	IF(lwp) WRITE(numout,*) ' zgr_bat_zoom : modify the level bathymetry for zoom domain'
499	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~'
500
501	! Zoom domain
502	! ===========
503
504	! Forced closed boundary if required
505	IF( lzoom_w ) mbathy( mi0(jpizoom):mi1(jpizoom) , : ) = 0
506	IF( lzoom_s ) mbathy( : , mj0(jpjzoom):mj1(jpjzoom) ) = 0
507	IF( lzoom_e ) mbathy( mi0(jpiglo+jpizoom-1):mi1(jpiglo+jpizoom-1) , : ) = 0
508	IF( lzoom_n ) mbathy( : , mj0(jpjglo+jpjzoom-1):mj1(jpjglo+jpjzoom-1) ) = 0
509
510	! Configuration specific domain modifications
511	! (here, ORCA arctic configuration: suppress Med Sea)
512	IF( cp_cfg == "orca" .AND. lzoom_arct ) THEN
513	SELECT CASE ( jp_cfg )
514	! ! =======================
515	CASE ( 2 ) ! ORCA_R2 configuration
516	! ! =======================
517	IF(lwp) WRITE(numout,*) ' ORCA R2 arctic zoom: suppress the Med Sea'
518
519	ii0 = 141 ; ii1 = 162 ! Sea box i,j indices
520	ij0 = 98 ; ij1 = 110
521	! ! =======================
522	CASE ( 05 ) ! ORCA_R05 configuration
523	! ! =======================
524	IF(lwp) WRITE(numout,*) ' ORCA R05 arctic zoom: suppress the Med Sea'
525
526	ii0 = 563 ; ii1 = 642 ! zero over the Med Sea boxe
527	ij0 = 314 ; ij1 = 370
528
529	END SELECT
530	!
531	mbathy( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0 ! zero over the Med Sea boxe
532	!
533	ENDIF
534
535	END SUBROUTINE zgr_bat_zoom
536
537
538	SUBROUTINE zgr_bat_ctl
539	!!----------------------------------------------------------------------
540	!! * ROUTINE zgr_bat_ctl *
541	!!
542	!! ** Purpose : check the bathymetry in levels
543	!!
544	!! ** Method : The array mbathy is checked to verified its consistency
545	!! with the model options. in particular:
546	!! mbathy must have at least 1 land grid-points (mbathy<=0)
547	!! along closed boundary.
548	!! mbathy must be cyclic IF jperio=1.
549	!! mbathy must be lower or equal to jpk-1.
550	!! isolated ocean grid points are suppressed from mbathy
551	!! since they are only connected to remaining
552	!! ocean through vertical diffusion.
553	!! C A U T I O N : mbathy will be modified during the initializa-
554	!! tion phase to become the number of non-zero w-levels of a water
555	!! column, with a minimum value of 1.
556	!!
557	!! ** Action : - update mbathy: level bathymetry (in level index)
558	!! - update bathy : meter bathymetry (in meters)
559	!!
560	!! History :
561	!! 9.0 ! 03-08 (G. Madec) Original code
562	!! 9.0 ! 05-10 (A. Beckmann) modifications for s-ccordinates
563	!!----------------------------------------------------------------------
564	!! * Local declarations
565	INTEGER :: ji, jj, jl ! dummy loop indices
566	INTEGER :: &
567	icompt, ibtest, ikmax ! temporary integers
568	REAL(wp), DIMENSION(jpi,jpj) :: &
569	zbathy ! temporary workspace
570	!!----------------------------------------------------------------------
571
572	IF(lwp) WRITE(numout,*)
573	IF(lwp) WRITE(numout,*) ' zgr_bat_ctl : check the bathymetry'
574	IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~'
575
576	! ================
577	! Bathymetry check
578	! ================
579
580	IF( .NOT. lk_cfg_1d ) THEN
581
582	! Suppress isolated ocean grid points
583
584	IF(lwp) WRITE(numout,*)
585	IF(lwp) WRITE(numout,*)' suppress isolated ocean grid points'
586	IF(lwp) WRITE(numout,*)' -----------------------------------'
587
588	icompt = 0
589
590	DO jl = 1, 2
591
592	IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
593	mbathy( 1 ,:) = mbathy(jpim1,:)
594	mbathy(jpi,:) = mbathy( 2 ,:)
595	ENDIF
596	DO jj = 2, jpjm1
597	DO ji = 2, jpim1
598	ibtest = MAX( mbathy(ji-1,jj), mbathy(ji+1,jj), &
599	& mbathy(ji,jj-1),mbathy(ji,jj+1) )
600	IF( ibtest < mbathy(ji,jj) ) THEN
601	IF(lwp) WRITE(numout,*) ' the number of ocean level at ', &
602	& 'grid-point (i,j) = ',ji,jj,' is changed from ', mbathy(ji,jj),' to ', ibtest
603	mbathy(ji,jj) = ibtest
604	icompt = icompt + 1
605	ENDIF
606	END DO
607	END DO
608
609	END DO
610	IF( icompt == 0 ) THEN
611	IF(lwp) WRITE(numout,*)' no isolated ocean grid points'
612	ELSE
613	IF(lwp) WRITE(numout,*)' ',icompt,' ocean grid points suppressed'
614	ENDIF
615	IF( lk_mpp ) THEN
616	zbathy(:,:) = FLOAT( mbathy(:,:) )
617	CALL lbc_lnk( zbathy, 'T', 1. )
618	mbathy(:,:) = INT( zbathy(:,:) )
619	ENDIF
620
621	! 3.2 East-west cyclic boundary conditions
622
623	IF( nperio == 0 ) THEN
624	IF(lwp) WRITE(numout,*) ' mbathy set to 0 along east and west', &
625	' boundary: nperio = ', nperio
626	IF( lk_mpp ) THEN
627	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
628	IF( jperio /= 1 ) mbathy(1,:) = 0
629	ENDIF
630	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
631	IF( jperio /= 1 ) mbathy(nlci,:) = 0
632	ENDIF
633	ELSE
634	IF( ln_zco .OR. ln_zps ) THEN
635	mbathy( 1 ,:) = 0
636	mbathy(jpi,:) = 0
637	ELSE
638	mbathy( 1 ,:) = jpkm1
639	mbathy(jpi,:) = jpkm1
640	ENDIF
641	ENDIF
642	ELSEIF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
643	IF(lwp) WRITE(numout,*)' east-west cyclic boundary conditions', &
644	' on mbathy: nperio = ', nperio
645	mbathy( 1 ,:) = mbathy(jpim1,:)
646	mbathy(jpi,:) = mbathy( 2 ,:)
647	ELSEIF( nperio == 2 ) THEN
648	IF(lwp) WRITE(numout,*) ' equatorial boundary conditions', &
649	' on mbathy: nperio = ', nperio
650	ELSE
651	IF(lwp) WRITE(numout,*) ' e r r o r'
652	IF(lwp) WRITE(numout,*) ' parameter , nperio = ', nperio
653	! STOP 'dom_mba'
654	ENDIF
655
656	! Set to zero mbathy over islands if necessary (lk_isl=F)
657	IF( .NOT. lk_isl ) THEN ! No island
658	IF(lwp) WRITE(numout,*)
659	IF(lwp) WRITE(numout,*) ' mbathy set to 0 over islands'
660	IF(lwp) WRITE(numout,*) ' ----------------------------'
661
662	mbathy(:,:) = MAX( 0, mbathy(:,:) )
663
664	! Boundary condition on mbathy
665	IF( .NOT.lk_mpp ) THEN
666	!!bug ??? y reflechir!
667	! ... mono- or macro-tasking: T-point, >0, 2D array, no slab
668	zbathy(:,:) = FLOAT( mbathy(:,:) )
669	CALL lbc_lnk( zbathy, 'T', 1. )
670	mbathy(:,:) = INT( zbathy(:,:) )
671	ENDIF
672
673	ENDIF
674
675	ENDIF
676
677	! Number of ocean level inferior or equal to jpkm1
678
679	ikmax = 0
680	DO jj = 1, jpj
681	DO ji = 1, jpi
682	ikmax = MAX( ikmax, mbathy(ji,jj) )
683	END DO
684	END DO
685	!!! test a faire: ikmax = MAX( mbathy(:,:) ) ???
686
687	IF( ikmax > jpkm1 ) THEN
688	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' > jpk-1'
689	IF(lwp) WRITE(numout,*) ' change jpk to ',ikmax+1,' to use the exact ead bathymetry'
690	ELSE IF( ikmax < jpkm1 ) THEN
691	IF(lwp) WRITE(numout,*) ' maximum number of ocean level = ', ikmax,' < jpk-1'
692	IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1
693	ENDIF
694
695	IF( lwp .AND. nprint == 1 ) THEN
696	WRITE(numout,*)
697	WRITE(numout,*) ' bathymetric field '
698	WRITE(numout,*) ' ------------------'
699	WRITE(numout,*) ' number of non-zero T-levels '
700	CALL prihin( mbathy, jpi, jpj, 1, jpi, &
701	1 , 1 , jpj, 1, 3 , &
702	numout )
703	WRITE(numout,*)
704	ENDIF
705
706	END SUBROUTINE zgr_bat_ctl
707
708
709	SUBROUTINE zgr_zco
710	!!----------------------------------------------------------------------
711	!! * ROUTINE zgr_zco *
712	!!
713	!! ** Purpose : define the z-coordinate system
714	!!
715	!! ** Method : set 3D coord. arrays to reference 1D array if lk_zco=T
716	!!
717	!! History :
718	!! ! 06-04 (R. Benshila, G. Madec) Original code
719	!!----------------------------------------------------------------------
720	INTEGER :: jk
721	!!----------------------------------------------------------------------
722
723	IF( .NOT.lk_zco ) THEN
724	DO jk = 1, jpk
725	fsdept(:,:,jk) = gdept_0(jk)
726	fsdepw(:,:,jk) = gdepw_0(jk)
727	fsde3w(:,:,jk) = gdepw_0(jk)
728	fse3t (:,:,jk) = e3t_0(jk)
729	fse3u (:,:,jk) = e3t_0(jk)
730	fse3v (:,:,jk) = e3t_0(jk)
731	fse3f (:,:,jk) = e3t_0(jk)
732	fse3w (:,:,jk) = e3w_0(jk)
733	fse3uw(:,:,jk) = e3w_0(jk)
734	fse3vw(:,:,jk) = e3w_0(jk)
735	END DO
736	ENDIF
737
738	END SUBROUTINE zgr_zco
739
740
741
742	# include "domzgr_zps.h90"
743
744
745	#if ! defined key_zco
746	!!----------------------------------------------------------------------
747	!! NOT 'key_zco' : 3D gdep. and e3.
748	!!----------------------------------------------------------------------
749
750	SUBROUTINE zgr_sco
751	!!----------------------------------------------------------------------
752	!! * ROUTINE zgr_sco *
753	!!
754	!! ** Purpose : define the s-coordinate system
755	!!
756	!! ** Method : s-coordinate
757	!! The depth of model levels is defined as the product of an
758	!! analytical function by the local bathymetry, while the vertical
759	!! scale factors are defined as the product of the first derivative
760	!! of the analytical function by the bathymetry.
761	!! (this solution save memory as depth and scale factors are not
762	!! 3d fields)
763	!! - Read bathymetry (in meters) at t-point and compute the
764	!! bathymetry at u-, v-, and f-points.
765	!! hbatu = mi( hbatt )
766	!! hbatv = mj( hbatt )
767	!! hbatf = mi( mj( hbatt ) )
768	!! - Compute gsigt, gsigw, esigt, esigw from an analytical
769	!! function and its derivative given as statement function.
770	!! gsigt(k) = fssig (k+0.5)
771	!! gsigw(k) = fssig (k )
772	!! esigt(k) = fsdsig(k+0.5)
773	!! esigw(k) = fsdsig(k )
774	!! This routine is given as an example, it must be modified
775	!! following the user s desiderata. nevertheless, the output as
776	!! well as the way to compute the model levels and scale factors
777	!! must be respected in order to insure second order a!!uracy
778	!! schemes.
779	!!
780	!! Reference :
781	!! Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408.
782	!!
783	!! History :
784	!! ! 95-12 (G. Madec) Original code : s vertical coordinate
785	!! ! 97-07 (G. Madec) lbc_lnk call
786	!! ! 97-04 (J.-O. Beismann)
787	!! 8.5 ! 02-06 (G. Madec) F90: Free form and module
788	!! 9.0 ! 05-10 (A. Beckmann) new stretching function
789	!!----------------------------------------------------------------------
790	!! * Local declarations
791	INTEGER :: ji, jj, jk, jl
792	REAL(wp) :: fssig, fsdsig, pfkk
793
794	INTEGER :: iip1, ijp1, iim1, ijm1
795	REAL(wp) :: &
796	fssigt, fssigw, zcoeft, zcoefw, &
797	zrmax, ztaper
798
799	REAL(wp), DIMENSION(jpi,jpj) :: &
800	zenv, ztmp, zmsk, zri, zrj, zhbat
801
802	NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max
803	!!----------------------------------------------------------------------
804
805	! a. Sigma stretching coefficient
806	fssigt(pfkk) = (tanh(theta(-(pfkk-0.5)/FLOAT(jpkm1)+thetb))-tanh(thetbtheta)) &
807	* (cosh(theta)+cosh(theta(2.thetb-1.)))/(2.*sinh(theta))
808	fssigw(pfkk) = (tanh(theta(-(pfkk-1.0)/FLOAT(jpkm1)+thetb))-tanh(thetbtheta)) &
809	* (cosh(theta)+cosh(theta(2.thetb-1.)))/(2.*sinh(theta))
810
811	! b. Vertical derivative of sigma stretching coefficient
812
813	!!bug fsdsig(pfkk)= put here the analytical derivative of fssigt and w ...
814	!!----------------------------------------------------------------------
815
816	! Read Namelist nam_zgr_sco : sigma-stretching parameters
817	! -------------------------
818	REWIND( numnam )
819	READ ( numnam, nam_zgr_sco )
820
821	IF(lwp) THEN
822	WRITE(numout,*)
823	WRITE(numout,*) 'dom:zgr_sco : s-coordinate or hybrid z-s-coordinate'
824	WRITE(numout,*) '~~~~~~~~~~~'
825	WRITE(numout,*) ' Namelist nam_zgr_sco'
826	WRITE(numout,*) ' sigma-stretching coeffs '
827	WRITE(numout,*) ' maximum depth of s-bottom surface (>0) sbot_max = ' ,sbot_max
828	WRITE(numout,*) ' minimum depth of s-bottom surface (>0) sbot_min = ' ,sbot_min
829	WRITE(numout,*) ' surface control parameter (0<=theta<=20) theta = ', theta
830	WRITE(numout,*) ' bottom control parameter (0<=thetb<= 1) thetb = ', thetb
831	WRITE(numout,*) ' maximum cut-off r-value allowed r_max = ' , r_max
832	ENDIF
833
834	! ???
835	! ------------------
836	hift(:,:) = sbot_min
837	hifu(:,:) = sbot_min
838	hifv(:,:) = sbot_min
839	hiff(:,:) = sbot_min
840
841
842	! set maximum ocean depth
843	! -----------------------
844	DO jj = 1, jpj
845	DO ji = 1, jpi
846	bathy(ji,jj) = MIN( sbot_max, bathy(ji,jj) )
847	END DO
848	END DO
849
850	! Define the envelop bathymetry
851	! =============================
852	! Smooth the bathymetry (if required)
853
854	scosrf(:,:) = 0.e0 ! ocean surface depth (here zero: no under ice-shelf sea)
855	scobot(:,:) = bathy(:,:) ! ocean bottom depth
856
857
858	! use r-value to create hybrid coordinates
859
860	DO jj = 1, jpj
861	DO ji = 1, jpi
862	zenv(ji,jj) = MAX( bathy(ji,jj), sbot_min )
863	END DO
864	END DO
865
866	jl = 0
867	zrmax = 1.e0
868	! ! ================ !
869	DO WHILE ( jl <= 10000 .AND. zrmax > r_max ) ! Iterative loop !
870	! ! ================ !
871	jl = jl + 1
872
873	zrmax = 0.e0
874	zmsk(:,:) = 0.e0
875
876	DO jj = 1, nlcj
877	DO ji = 1, nlci
878	iip1 = MIN( ji+1, nlci ) ! force zri = 0 on last line (ji=ncli+1 to jpi)
879	ijp1 = MIN( jj+1, nlcj ) ! force zrj = 0 on last raw (jj=nclj+1 to jpj)
880	zri(ji,jj) = ABS( zenv(iip1,jj ) - zenv(ji,jj) ) / ( zenv(iip1,jj ) + zenv(ji,jj) )
881	zrj(ji,jj) = ABS( zenv(ji ,ijp1) - zenv(ji,jj) ) / ( zenv(ji ,ijp1) + zenv(ji,jj) )
882	zrmax = MAX( zrmax, zri(ji,jj), zrj(ji,jj) )
883	IF( zri(ji,jj) > r_max ) zmsk(ji ,jj ) = 1.0
884	IF( zri(ji,jj) > r_max ) zmsk(iip1,jj ) = 1.0
885	IF( zrj(ji,jj) > r_max ) zmsk(ji ,jj ) = 1.0
886	IF( zrj(ji,jj) > r_max ) zmsk(ji ,ijp1) = 1.0
887	END DO
888	END DO
889	IF( lk_mpp ) CALL mpp_max( zrmax ) ! max over the global domain
890	! lateral boundary condition on zmsk: keep 1 along closed boundary (use of MAX)
891	ztmp(:,:) = zmsk(:,:)
892	CALL lbc_lnk( zmsk, 'T', 1. )
893	DO jj = 1, nlcj
894	DO ji = 1, nlci
895	zmsk(ji,jj) = MAX( zmsk(ji,jj), ztmp(ji,jj) )
896	END DO
897	END DO
898
899	WRITE(numout,*) 'zgr_sco : iter= ',jl, ' rmax= ', zrmax, ' nb of pt= ', INT( SUM(zmsk(:,:) ) )
900
901	DO jj = 1, nlcj
902	DO ji = 1, nlci
903	iip1 = MIN( ji+1, nlci ) ! last line (ji=nlci)
904	ijp1 = MIN( jj+1, nlcj ) ! last raw (jj=nlcj)
905	iim1 = MAX( ji-1, 1 ) ! first line (ji=nlci)
906	ijm1 = MAX( jj-1, 1 ) ! first raw (jj=nlcj)
907	IF( zmsk(ji,jj) == 1.0 ) THEN
908	ztmp(ji,jj) = ( &
909	& zenv(iim1,ijp1)zmsk(iim1,ijp1) + zenv(ji,ijp1)zmsk(ji,ijp1) + zenv(iip1,ijp1)*zmsk(iip1,ijp1) &
910	& + zenv(iim1,jj )zmsk(iim1,jj ) + zenv(ji,jj ) 2.e0 + zenv(iip1,jj )*zmsk(iip1,jj ) &
911	& + zenv(iim1,ijm1)zmsk(iim1,ijm1) + zenv(ji,ijm1)zmsk(ji,ijm1) + zenv(iip1,ijm1)*zmsk(iip1,ijm1) &
912	& ) / ( &
913	& zmsk(iim1,ijp1) + zmsk(ji,ijp1) + zmsk(iip1,ijp1) &
914	& + zmsk(iim1,jj ) + 2.e0 + zmsk(iip1,jj ) &
915	& + zmsk(iim1,ijm1) + zmsk(ji,ijm1) + zmsk(iip1,ijm1) &
916	& )
917	ENDIF
918	END DO
919	END DO
920
921	DO jj = 1, nlcj
922	DO ji = 1, nlci
923	IF( zmsk(ji,jj) == 1.0 ) zenv(ji,jj) = MAX( ztmp(ji,jj), bathy(ji,jj) )
924	END DO
925	END DO
926
927	! Apply lateral boundary condition CAUTION: kept the value when the lbc field is zero
928	ztmp(:,:) = zenv(:,:)
929	CALL lbc_lnk( zenv, 'T', 1. )
930	DO jj = 1, nlcj
931	DO ji = 1, nlci
932	IF( zenv(ji,jj) == 0.e0 ) zenv(ji,jj) = ztmp(ji,jj)
933	END DO
934	END DO
935	! ! ================ !
936	END DO ! End loop !
937	! ! ================ !
938
939	! save envelop bathymetry in hbatt
940	hbatt(:,:) = zenv(:,:)
941
942	!!bug gm : CAUTION: tapering hard coded !!!! orca2 only
943	!!bug gm NOT valid in mpp ===> taper have been changed
944
945	DO jj = 1, jpj
946	DO ji = 1, jpi
947	!!bug mpp ztaper = EXP( -(FLOAT(jj-74)/10.)**2 )
948	ztaper = EXP( -(gphit(ji,jj)/8)**2 )
949	hbatt(ji,jj) = sbot_max * ztaper + hbatt(ji,jj) * (1.-ztaper)
950	END DO
951	END DO
952
953	DO jj = 1, jpj
954	zrmax = EXP( -(gphit(10,jj)/8)**2 )
955	ztaper = EXP( -(FLOAT(jj-74)/10.)**2 )
956	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) jj, (FLOAT(jj-74)/10.), ztaper,(gphit(10,jj)/8), zrmax
957	END DO
958
959	IF( nprint == 1 .AND. lwp ) THEN
960	WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )
961	WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )
962	ENDIF
963
964	! Control print
965	IF(lwp) THEN
966	WRITE(numout,*)
967	WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'
968	WRITE(numout,*)
969	CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0., numout )
970	ENDIF
971
972	! 1. computation of hbatu, hbatv, hbatf fields
973	! --------------------------------------------
974
975	hbatu(:,:) = sbot_min
976	hbatv(:,:) = sbot_min
977	hbatf(:,:) = sbot_min
978
979	IF(lwp) THEN
980	WRITE(numout,*)
981	WRITE(numout,*) ' zgr_sco: minimum depth of the envelop topography set to : ', sbot_min
982	WRITE(numout,*)
983	ENDIF
984
985	DO jj = 1, jpjm1
986	DO ji = 1, jpim1
987	hbatu(ji,jj)= 0.5 *( hbatt(ji ,jj)+hbatt(ji+1,jj ) )
988	hbatv(ji,jj)= 0.5 *( hbatt(ji ,jj)+hbatt(ji ,jj+1) )
989	hbatf(ji,jj)= 0.25*( hbatt(ji ,jj)+hbatt(ji ,jj+1) &
990	+hbatt(ji+1,jj)+hbatt(ji+1,jj+1) )
991	END DO
992	END DO
993
994	! Apply lateral boundary condition
995	! CAUTION: retain non zero value in the initial file
996	! this should be OK for orca cfg, not for EEL
997	zhbat(:,:) = hbatu(:,:)
998	CALL lbc_lnk( hbatu, 'U', 1. )
999	DO jj = 1, jpj
1000	DO ji = 1, jpi
1001	IF( hbatu(ji,jj) == 0.e0 ) THEN
1002	IF( zhbat(ji,jj) == 0.e0 ) hbatu(ji,jj) = sbot_min
1003	IF( zhbat(ji,jj) /= 0.e0 ) hbatu(ji,jj) = zhbat(ji,jj)
1004	ENDIF
1005	END DO
1006	END DO
1007
1008	zhbat(:,:) = hbatv(:,:)
1009	CALL lbc_lnk( hbatv, 'V', 1. )
1010	DO jj = 1, jpj
1011	DO ji = 1, jpi
1012	IF( hbatv(ji,jj) == 0.e0 ) THEN
1013	IF( zhbat(ji,jj) == 0.e0 ) hbatv(ji,jj) = sbot_min
1014	IF( zhbat(ji,jj) /= 0.e0 ) hbatv(ji,jj) = zhbat(ji,jj)
1015	ENDIF
1016	END DO
1017	END DO
1018
1019	zhbat(:,:) = hbatf(:,:)
1020	CALL lbc_lnk( hbatf, 'F', 1. )
1021	DO jj = 1, jpj
1022	DO ji = 1, jpi
1023	IF( hbatf(ji,jj) == 0.e0 ) THEN
1024	IF( zhbat(ji,jj) == 0.e0 ) hbatf(ji,jj) = sbot_min
1025	IF( zhbat(ji,jj) /= 0.e0 ) hbatf(ji,jj) = zhbat(ji,jj)
1026	ENDIF
1027	END DO
1028	END DO
1029
1030
1031	!!bug: key_helsinki a verifer
1032	hift(:,:) = MIN( hift(:,:), hbatt(:,:) )
1033	hifu(:,:) = MIN( hifu(:,:), hbatu(:,:) )
1034	hifv(:,:) = MIN( hifv(:,:), hbatv(:,:) )
1035	hiff(:,:) = MIN( hiff(:,:), hbatf(:,:) )
1036
1037	IF( nprint == 1 .AND. lwp ) THEN
1038	WRITE(numout,*) ' MAX val hif t ', MAXVAL( hift(:,:) ), ' f ', MAXVAL( hiff(:,:) ), &
1039	& ' u ', MAXVAL( hifu(:,:) ), ' v ', MAXVAL( hifv(:,:) )
1040	WRITE(numout,*) ' MIN val hif t ', MINVAL( hift(:,:) ), ' f ', MINVAL( hiff(:,:) ), &
1041	& ' u ', MINVAL( hifu(:,:) ), ' v ', MINVAL( hifv(:,:) )
1042	WRITE(numout,*) ' MAX val hbat t ', MAXVAL( hbatt(:,:) ), ' f ', MAXVAL( hbatf(:,:) ), &
1043	& ' u ', MAXVAL( hbatu(:,:) ), ' v ', MAXVAL( hbatv(:,:) )
1044	WRITE(numout,*) ' MIN val hbat t ', MINVAL( hbatt(:,:) ), ' f ', MINVAL( hbatf(:,:) ), &
1045	& ' u ', MINVAL( hbatu(:,:) ), ' v ', MINVAL( hbatv(:,:) )
1046	ENDIF
1047	!! helsinki
1048
1049	! 2. Computation of gsig and esig fields
1050	! --------------------------------------
1051
1052	! 2.1 Coefficients for model level depth at w- and t-levels
1053
1054	!!bug gm : change the def of statment function....
1055	DO jk = 1, jpk
1056	gsigw(jk) = -fssigw(FLOAT(jk))
1057	gsigt(jk) = -fssigt(FLOAT(jk))
1058	END DO
1059
1060	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk)
1061
1062	!!org DO jk = 1, jpk
1063	!!org gsigw(jk) = -fssig ( FLOAT(jk) )
1064	!!org gsigt(jk) = -fssig ( FLOAT(jk)+0.5)
1065	!!org END DO
1066
1067	! 2.2 Coefficients for vertical scale factors at w-, t- levels
1068
1069	!!bug gm : define it from analytical function, not like juste bellow....
1070	!! or betteroffer the 2 possibilities....
1071	DO jk=2,jpk
1072	esigw(jk)=gsigt(jk)-gsigt(jk-1)
1073	END DO
1074	DO jk=1,jpkm1
1075	esigt(jk)=gsigw(jk+1)-gsigw(jk)
1076	END DO
1077	esigw(1)=esigw(2)
1078	esigt(jpk)=esigt(jpkm1)
1079
1080	!!org DO jk = 1, jpk
1081	!!org esigw(jk)=fsdsig( FLOAT(jk))
1082	!!org esigt(jk)=fsdsig( FLOAT(jk)+0.5)
1083	!!org END DO
1084
1085	! 2.3 Coefficients for vertical depth as the sum of e3w scale factors
1086
1087	gsi3w(1) = 0.5 * esigw(1)
1088	DO jk = 2, jpk
1089	gsi3w(jk) = gsi3w(jk-1)+ esigw(jk)
1090	END DO
1091
1092	!!bug gm: depuw, depvw can be suppressed (modif in ldfslp) and depw=dep3w can be set (save 3 3D arrays)
1093	DO jk = 1, jpk
1094	!! change the solver stat....
1095	!! zcoeft = ( FLOAT(jk) - 0.5 ) / FLOAT(jpkm1)
1096	!! gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigt(jk)+hift(:,:)zcoeft)
1097	gdept (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigt(jk)+hift(:,:)(FLOAT(jk)-0.5)/FLOAT(jpkm1))
1098	!! zcoefw = ( FLOAT(jk) - 1.0 ) / FLOAT(jpkm1)
1099	!! gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigw(jk)+hift(:,:)zcoefw)
1100	!! gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsi3w(jk)+hift(:,:)zcoefw)
1101	gdepw (:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsigw(jk)+hift(:,:)FLOAT(jk-1)/FLOAT(jpkm1))
1102	gdep3w(:,:,jk) = (scosrf(:,:)+(hbatt(:,:)-hift(:,:))gsi3w(jk)+hift(:,:)FLOAT(jk-1)/FLOAT(jpkm1))
1103	END DO
1104
1105	!!bug gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays)
1106	DO jj = 1, jpj
1107	DO ji = 1, jpi
1108	DO jk = 1, jpk
1109	e3t(ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigt(jk) + hift(ji,jj)/FLOAT(jpkm1))
1110	e3u(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigt(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1111	e3v(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigt(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1112	e3f(ji,jj,jk)=((hbatf(ji,jj)-hiff(ji,jj))*esigt(jk) + hiff(ji,jj)/FLOAT(jpkm1))
1113
1114	e3w (ji,jj,jk)=((hbatt(ji,jj)-hift(ji,jj))*esigw(jk) + hift(ji,jj)/FLOAT(jpkm1))
1115	e3uw(ji,jj,jk)=((hbatu(ji,jj)-hifu(ji,jj))*esigw(jk) + hifu(ji,jj)/FLOAT(jpkm1))
1116	e3vw(ji,jj,jk)=((hbatv(ji,jj)-hifv(ji,jj))*esigw(jk) + hifv(ji,jj)/FLOAT(jpkm1))
1117	END DO
1118	END DO
1119	END DO
1120
1121	! HYBRID
1122	DO jj = 1, jpj
1123	DO ji = 1, jpi
1124	DO jk = 1, jpkm1
1125	IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
1126	IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0
1127	END DO
1128	END DO
1129	END DO
1130
1131	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
1132
1133
1134	! ===========
1135	! Zoom domain
1136	! ===========
1137
1138	IF( lzoom ) CALL zgr_bat_zoom
1139
1140	! 2.4 Control print
1141
1142	IF(lwp) THEN
1143	WRITE(numout,*)
1144	WRITE(numout,*) ' domzgr: vertical coefficients for model level'
1145	WRITE(numout,9400)
1146	WRITE(numout,9410) (jk,gsigt(jk),gsigw(jk),esigt(jk),esigw(jk),gsi3w(jk),jk=1,jpk)
1147	ENDIF
1148	9400 FORMAT(9x,' level gsigt gsigw esigt esigw gsi3w')
1149	9410 FORMAT(10x,i4,5f11.4)
1150
1151	IF(lwp) THEN
1152	WRITE(numout,*)
1153	WRITE(numout,*) ' domzgr: vertical coordinates : point (1,1,k) bathy = ', bathy(1,1), hbatt(1,1)
1154	WRITE(numout,*) ' ~~~~~~ --------------------'
1155	WRITE(numout,9420)
1156	WRITE(numout,9430) (jk,fsdept(1,1,jk),fsdepw(1,1,jk), &
1157	fse3t (1,1,jk),fse3w (1,1,jk),jk=1,jpk)
1158	DO jj = mj0(20), mj1(20)
1159	DO ji = mi0(20), mi1(20)
1160	WRITE(numout,*)
1161	WRITE(numout,*) ' domzgr: vertical coordinates : point (20,20,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1162	WRITE(numout,*) ' ~~~~~~ --------------------'
1163	WRITE(numout,9420)
1164	WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk), &
1165	& fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk)
1166	END DO
1167	END DO
1168	DO jj = mj0(74), mj1(74)
1169	DO ji = mi0(100), mi1(100)
1170	WRITE(numout,*)
1171	WRITE(numout,*) ' domzgr: vertical coordinates : point (100,74,k) bathy = ', bathy(ji,jj), hbatt(ji,jj)
1172	WRITE(numout,*) ' ~~~~~~ --------------------'
1173	WRITE(numout,9420)
1174	WRITE(numout,9430) (jk,fsdept(ji,jj,jk),fsdepw(ji,jj,jk), &
1175	& fse3t (ji,jj,jk),fse3w (ji,jj,jk),jk=1,jpk)
1176	END DO
1177	END DO
1178	ENDIF
1179
1180	9420 FORMAT(9x,' level gdept gdepw gde3w e3t e3w ')
1181	9430 FORMAT(10x,i4,4f9.2)
1182
1183	!!bug gm no more necessary? if ! defined key_helsinki
1184	DO jk = 1, jpk
1185	DO jj = 1, jpj
1186	DO ji = 1, jpi
1187	IF( fse3w(ji,jj,jk) <= 0. .OR. fse3t(ji,jj,jk) <= 0. ) THEN
1188	IF(lwp) THEN
1189	WRITE(numout,*)
1190	WRITE(numout,*) ' e r r o r : e3w or e3t =< 0 '
1191	WRITE(numout,*) ' ========= --------------- '
1192	WRITE(numout,*)
1193	WRITE(numout,*) ' point (i,j,k)= ',ji,jj,jk
1194	WRITE(numout,*)
1195	ENDIF
1196	STOP 'domzgr'
1197	ENDIF
1198	IF( fsdepw(ji,jj,jk) < 0. .OR. fsdept(ji,jj,jk) < 0. ) THEN
1199	IF(lwp) THEN
1200	WRITE(numout,*)
1201	WRITE(numout,*) ' e r r o r : gdepw or gdept < 0 '
1202	WRITE(numout,*) ' ========= ------------------ '
1203	WRITE(numout,*)
1204	WRITE(numout,*) ' point (i,j,k)= ', ji, jj, jk
1205	WRITE(numout,*)
1206	ENDIF
1207	STOP 'domzgr'
1208	ENDIF
1209	END DO
1210	END DO
1211	END DO
1212	!!bug gm #endif
1213
1214	END SUBROUTINE zgr_sco
1215
1216	#else
1217	!!----------------------------------------------------------------------
1218	!! Default option : Empty routine
1219	!!----------------------------------------------------------------------
1220	SUBROUTINE zgr_sco
1221	END SUBROUTINE zgr_sco
1222	#endif
1223
1224
1225	!!======================================================================
1226	END MODULE domzgr

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

Download in other formats: