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

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

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

Download in other formats: