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

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

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

Download in other formats: