Context Navigation

mes.F90 @ 15733

Last change on this file since 15733 was 15278, checked in by dbruciaferri, 3 years ago
reorganising code in more general structure and adding partial-steps for MEs
File size: 64.3 KB

Line
1	MODULE mes
2	!!==============================================================================
3	!! * MODULE mes *
4	!! Ocean initialization : Multiple Enveloped s coordinate (MES)
5	!!==============================================================================
6	!! NEMO 3.6 ! 2017-05 D. Bruciaferri
7	!! NEMO 4.0.4 ! 2021-07 D. Bruciaferri
8	!!----------------------------------------------------------------------
9	!
10	USE oce ! ocean variables
11	USE dom_oce ! ocean domain
12	USE closea ! closed seas
13	!
14	USE in_out_manager ! I/O manager
15	USE iom ! I/O library
16	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
17	USE lib_mpp ! distributed memory computing library
18	USE wrk_nemo ! Memory allocation
19	USE timing ! Timing
20
21	IMPLICIT NONE
22	PRIVATE
23
24	PUBLIC mes_build ! called by zgrmes.F90
25	!
26	! Envelopes and stretching parameters
27	CHARACTER(lc) :: ctlmes ! control message error (lc=256)
28	INTEGER, PARAMETER :: max_nn_env = 5 ! Maximum allowed number of envelopes.
29	INTEGER :: tot_env ! Tot number of requested envelopes
30	LOGICAL :: ln_envl(max_nn_env) ! Array with flags (T/F) specifing which envelope is used
31	INTEGER :: nn_strt(max_nn_env) ! Array specifing the stretching function for each envelope:
32	! Madec 1996 (0), Song and Haidvogel 1994 (1)
33	REAL(wp) :: max_dep_env(max_nn_env) ! global maximum depth of envelopes
34	REAL(wp) :: min_dep_env(max_nn_env) ! global minimum depth of envelopes
35	INTEGER :: nn_slev(max_nn_env) ! Array specifing number of levels of each enveloped vertical zone
36	REAL(wp) :: rn_e_hc(max_nn_env) ! Array specifing critical depth for transition to stretched
37	! coordinates of each envelope
38	REAL(wp) :: rn_e_th(max_nn_env) ! Array specifing surface control parameter (0<=th<=20) of SH94
39	! or rn_zs of SF12 for each vertical sub-zone
40	REAL(wp) :: rn_e_bb(max_nn_env) ! Array specifing bottom control parameter (0<=bb<=1) of SH94
41	! or rn_zb_b of SF12 for each vertical sub-zone
42	REAL(wp) :: rn_e_ba(max_nn_env) ! Array specifing bottom control parameter rn_zb_a of SF12 for
43	! each vertical sub-zone
44	REAL(wp) :: rn_e_al(max_nn_env) ! Array specifing stretching parameter rn_alpha of SF12 for
45	! each vertical sub-zone
46	LOGICAL, PUBLIC :: ln_pst_mes ! To apply partial steps to MEs (.TRUE.) or not (.FALSE.).
47	LOGICAL, PUBLIC :: ln_pst_l2g ! To apply partial steps to the transition zone (.TRUE.)
48	! or not (.FALSE.) when ln_loc_zgr = .TRUE.
49	REAL, PUBLIC :: rn_e3pst_min ! partial steps parameters (require ln_pst_mes = .TRUE.)
50	REAL, PUBLIC :: rn_e3pst_rat ! partial steps parameters (require ln_pst_mes = .TRUE.)
51	!
52	REAL(wp), POINTER, DIMENSION(:,:,:) :: envlt ! array for the envelopes
53
54	!! * Substitutions
55	# include "vectopt_loop_substitute.h90"
56
57	CONTAINS
58
59	! =====================================================================================================
60
61	SUBROUTINE mes_build
62	!!-----------------------------------------------------------------------------
63	!! * ROUTINE mes_build *
64	!!
65	!! ** Purpose : define the Multi Enveloped S-coordinate (MES) system
66	!!
67	!! ** Method : s-coordinate defined respect to multiple
68	!! externally defined envelope as detailed in
69	!!
70	!! Bruciaferri, Shapiro, Wobus, 2018. Oce. Dyn.
71	!! https://doi.org/10.1007/s10236-018-1189-x
72	!!
73	!! Three options for stretching are given:
74	!!
75	!! *) nn_strt = 0 : Madec et al 1996 cosh/tanh function
76	!! *) nn_strt = 1 : Song and Haidvogel 1994 sinh/tanh function
77	!! *) nn_strt = 2 : Siddorn and Furner gamma function
78	!!
79	!!----------------------------------------------------------------------
80	!! SKETCH of the GEOMETRY OF A MEs-COORDINATE SYSTEM
81	!!
82	!! Lines represent W-levels:
83	!!
84	!! 0: First s-level part. The total number
85	!! of levels in this part is controlled
86	!! by the nn_slev(1) namelist parameter.
87	!!
88	!! Depth first W-lev: 0 m (surfcace)
89	!! Depth last W-lev: depth of 1st envelope
90	!!
91	!! o: Second s-level part. The total number
92	!! of levels in this part is controlled
93	!! by the nn_slev(2) namelist parameter.
94	!!
95	!! Depth last W-lev: depth of 2nd envelope
96	!!
97	!! @: Third s-level part. The total number
98	!! of levels in this part is controlled
99	!! by the nn_slev(3) namelist parameter.
100	!!
101	!! Depth last W-lev: depth of 3rd envelope
102	!!
103	!! z \|~~~~~~~~~~~~~~~~~~~~~~~0~~~~~~~~~~~~~~~~~~~~~~~ SURFACE
104	!! \|
105	!! \|-----------------------0----------------------- ln_envl(1) = true, nn_slev(1) = 3
106	!! \|
107	!! \|=======================0======================= ENVELOPE 1
108	!! \|
109	!! \|-----------------------o-----------------------
110	!! \| ln_envl(2) = true, nn_slev(2) = 3
111	!! \|-----------------------o-----------------------
112	!! \|
113	!! \|¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬o¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ ENVELOPE 2
114	!! \|
115	!! \|-----------------------@----------------------- ln_envl(3) = true, nn_slev(3) = 2
116	!! \|
117	!! \|_______________________@_______________________ ENVELOPE 3
118	!! \\|/
119	!!
120	!!
121	!!-----------------------------------------------------------------------------
122	!! Author: Diego Bruciaferri (University of Plymouth), September 2017
123	!!-----------------------------------------------------------------------------
124
125	!
126	INTEGER :: ji, jj, jk, jl, je ! dummy loop argument
127	INTEGER :: eii, eij, irnk ! dummy loop argument
128	INTEGER :: num_s, s_1st, ind ! for loops over envelopes
129	INTEGER :: num_s_up, num_s_dw ! for loops over envelopes
130	REAL(wp) :: D1s_up, D1s_dw ! for loops over envelopes
131	INTEGER :: ibcbeg, ibcend ! boundary condition for cubic splines
132	INTEGER, PARAMETER :: n = 2 ! number of considered points for each
133	! ! interval
134	REAL(wp) :: c(4,n) ! matrix for cubic splines coefficents
135	REAL(wp) :: d(2,2) ! matrix for depth values and depth
136	! ! derivative at intervals' boundaries
137	REAL(wp) :: tau(n) ! abscissas of intervals' boundaries
138	REAL(wp) :: coefa, coefb ! for loops over envelopes
139	REAL(wp) :: alpha, env_hc ! for loops over envelopes
140	REAL(wp) :: zcoeft, zcoefw ! temporary scalars
141	REAL(wp) :: max_env_up, min_env_up ! to identify geopotential levels
142	REAL(wp) :: max_env_dw, min_env_dw ! to identify geopotential levels
143	REAL(wp) :: mindep
144	!
145	REAL(wp), DIMENSION(max_nn_env) :: rn_ebot_max
146	REAL(wp), DIMENSION(jpk) :: z_gsigw1, z_gsigt1 ! 1D arrays for debugging
147	REAL(wp), DIMENSION(jpk) :: gdepw1, gdept1 ! 1D arrays for debugging
148	REAL(wp), DIMENSION(jpk) :: e3t1, e3w1 ! 1D arrays for debugging
149	!
150	REAL(wp), DIMENSION(jpi,jpj) :: env_up, env_dw ! for loops over envelopes
151	REAL(wp), DIMENSION(jpi,jpj) :: env0, env1, env2, env3 ! for loops over envelopes
152	! ! for cubic splines
153	REAL(wp), DIMENSION(jpi,jpj) :: pmsk ! working array
154	!
155	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_gsigw3, z_gsigt3 ! 3D arrays for debugging
156	!!----------------------------------------------------------------------
157	!
158	IF( nn_timing == 1 ) CALL timing_start('mes_build')
159	!
160	CALL mes_init
161	!
162	! Initialize mbathy to the maximum ocean level available
163	mbathy(:,:) = jpkm1
164
165	! Defining scosrf and scobot
166	scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea)
167	scobot(:,:) = bathy(:,:) ! ocean bottom depth
168
169	!========================
170	! Compute 3D MEs mesh
171	!========================
172
173	DO je = 1, tot_env ! LOOP over all the requested envelopes
174
175	IF (je == 1) THEN
176	s_1st = 1
177	num_s = nn_slev(je)
178	env_up = scosrf ! surface
179	ELSE
180	s_1st = s_1st + num_s - 1
181	num_s = nn_slev(je) + 1
182	env_up = envlt(:,:,je-1)
183	ENDIF
184	env_dw = envlt(:,:,je)
185
186	IF ( isodd(je) ) THEN
187	!
188	env_hc = rn_e_hc(je)
189	coefa = rn_e_th(je)
190	coefb = rn_e_bb(je)
191	alpha = rn_e_al(je)
192	!
193	DO jk = 1, num_s
194	ind = s_1st+jk-1
195	!
196	DO jj = 1,jpj
197	DO ji = 1,jpi
198	!
199	! -------------------------------------------
200	! Computing MEs-coordinates (-1 <= MEs <= 0)
201	!
202	! shallow water, uniform sigma
203	IF (env_dw(ji,jj) < env_hc) THEN
204	! W-GRID
205	zcoefw = -sigma(1, jk, num_s)
206	! T-GRID
207	zcoeft = -sigma(0, jk, num_s)
208	!IF (nn_strt(je) == 2) THEN
209	! ! shallow water, z coordinates
210	! ! SF12 in AMM* configurations uses
211	! ! this (ln_sigcrit=.false.)
212	! zcoefw = zcoefw*(env_hc/(env_dw(ji,jj)-env_up(ji,jj)))
213	! zcoeft = zcoeft*(env_hc/(env_dw(ji,jj)-env_up(ji,jj)))
214	!ENDIF
215	! deep water, stretched s
216	ELSE
217	IF (nn_strt(je) == 2) THEN
218	coefa = rn_e_th(je) / (env_dw(ji,jj)-env_up(ji,jj))
219	coefb = rn_e_bb(je) + ((env_dw(ji,jj)-env_up(ji,jj))*rn_e_ba(je))
220	coefb = 1.0_wp-(coefb/(env_dw(ji,jj)-env_up(ji,jj)))
221	END IF
222	! W-GRID
223	zcoefw = -s_coord( sigma(1,jk,num_s), coefa, &
224	coefb, alpha, num_s, nn_strt(je) )
225	! T-GRID
226	zcoeft = -s_coord( sigma(0,jk,num_s), coefa, &
227	coefb, alpha, num_s, nn_strt(je) )
228	ENDIF
229	!
230	z_gsigw3(ji,jj,ind) = zcoefw
231	z_gsigt3(ji,jj,ind) = zcoeft
232	!
233	! --------------------------------------------------
234	! Computing model levels depths
235	IF (nn_strt(je) /= 2) THEN
236	gdept_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoeft, &
237	-sigma(0, jk, num_s), env_hc)
238
239	gdepw_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoefw, &
240	-sigma(1, jk, num_s), env_hc)
241	ELSE
242	gdept_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoeft, &
243	-sigma(0, jk, num_s), 0._wp)
244
245	gdepw_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoefw, &
246	-sigma(1, jk, num_s), 0._wp)
247	END IF
248	END DO
249	END DO
250	END DO
251	!
252	ELSE
253	!
254	! Interval's endpoints TAU are 0 and 1.
255	tau(1) = 0._wp
256	tau(n) = 1._wp
257	!
258	IF ( je == 2 ) THEN
259	env0 = scosrf
260	ELSE
261	env0 = envlt(:,:,je-2)
262	END IF
263	env1 = env_up
264	env2 = env_dw
265	IF ( je < tot_env ) env3 = envlt(:,:,je+1)
266	!
267	DO jj = 1,jpj
268	DO ji = 1,jpi
269	d(1,1) = env_up(ji,jj)
270	d(1,2) = env_dw(ji,jj)
271	!
272	! Computing derivative analytically
273	!
274	! 1. Derivative for upper envelope
275	ibcbeg = 1
276	IF (nn_strt(je-1) == 2) THEN
277	alpha = rn_e_al(je-1)
278	num_s_up = nn_slev(je-1)
279	coefa = rn_e_th(je-1) / &
280	(env1(ji,jj)-env0(ji,jj))
281	coefb = rn_e_bb(je-1) + &
282	((env1(ji,jj)-env0(ji,jj))*rn_e_ba(je-1))
283	coefb = 1.0_wp-(coefb/(env1(ji,jj)-env0(ji,jj)))
284	ELSE
285	alpha = 0._wp
286	num_s_up = 1
287	coefa = rn_e_th(je-1)
288	coefb = rn_e_bb(je-1)
289	END IF
290	D1s_up = D1s_coord(-1._wp, coefa, coefb, &
291	alpha, num_s_up, nn_strt(je-1))
292	IF (nn_strt(je-1) == 2) THEN
293	d(2,1) = D1z_mes(env0(ji,jj), env1(ji,jj), D1s_up, &
294	0._wp, 1)
295	ELSE
296	d(2,1) = D1z_mes(env0(ji,jj), env1(ji,jj), D1s_up, &
297	rn_e_hc(je-1), 1)
298	END IF
299	!
300	! 2. Derivative for lower envelope
301	IF ( je < tot_env ) THEN
302	ibcend = 1 ! Bound. cond for the 1st derivative
303	IF (nn_strt(je+1) == 2) THEN
304	alpha = rn_e_al(je+1)
305	num_s_dw = nn_slev(je+1)
306	coefa = rn_e_th(je+1) / &
307	(env3(ji,jj)-env2(ji,jj))
308	coefb = rn_e_bb(je+1) + &
309	((env3(ji,jj)-env2(ji,jj))*rn_e_ba(je+1))
310	coefb = 1.0_wp-(coefb/(env3(ji,jj)-env2(ji,jj)))
311	ELSE
312	alpha = 0._wp
313	num_s_dw = 1
314	coefa = rn_e_th(je+1)
315	coefb = rn_e_bb(je+1)
316	END IF
317	D1s_dw = D1s_coord(0._wp, coefa, coefb, &
318	alpha, num_s_dw, nn_strt(je+1))
319
320	IF (nn_strt(je+1) == 2) THEN
321	d(2,2) = D1z_mes(env2(ji,jj), env3(ji,jj), D1s_dw, &
322	0._wp, 1)
323	ELSE
324	d(2,2) = D1z_mes(env2(ji,jj), env3(ji,jj), D1s_dw, &
325	rn_e_hc(je+1), 1)
326	END IF
327	ELSE
328	ibcend = 2 ! Bound. cond for the 2nd derivative
329	d(2,2) = 0._wp ! 2nd der = 0
330	ENDIF
331
332	! Computing spline's coefficients
333	!IF ( lwp ) THEN
334	! WRITE(numout,*) "BOUNDARY COND. to CONSTRAIN SPLINES:"
335	! WRITE(numout,*) "ibcbeg: ", ibcbeg
336	! WRITE(numout,*) "ibcend: ", ibcend
337	!ENDIF
338
339	c = cub_spl(tau, d, n, ibcbeg, ibcend )
340
341	env_hc = rn_e_hc(je)
342	coefa = rn_e_th(je)
343	coefb = rn_e_bb(je)
344	alpha = rn_e_al(je)
345
346	DO jk = 1, num_s
347	ind = s_1st+jk-1
348	!
349	! Computing stretched distribution of interpolation points
350	! W-GRID
351	zcoefw = -s_coord( sigma(1,jk,num_s), coefa, &
352	coefb, alpha, num_s, nn_strt(je) )
353	! T-GRID
354	zcoeft = -s_coord( sigma(0,jk,num_s), coefa, &
355	coefb, alpha, num_s, nn_strt(je) )
356	!
357	z_gsigw3(ji,jj,ind) = zcoefw
358	z_gsigt3(ji,jj,ind) = zcoeft
359	!
360	gdept_0(ji,jj,ind) = z_cub_spl(zcoeft, tau, c)
361	gdepw_0(ji,jj,ind) = z_cub_spl(zcoefw, tau, c)
362	END DO
363
364	END DO
365	END DO
366
367	END IF
368
369	END DO
370
371	! =============================================
372	! 4. COMPUTE e3t, e3w for ALL general levels
373	! as finite differences
374	! =============================================
375	!
376	DO jk=1,jpkm1
377	e3t_0(:,:,jk) = gdepw_0(:,:,jk+1) - gdepw_0(:,:,jk)
378	e3w_0(:,:,jk+1) = gdept_0(:,:,jk+1) - gdept_0(:,:,jk)
379	ENDDO
380	! Surface
381	jk = 1
382	e3w_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,1) - gdepw_0(:,:,1))
383	!
384	! Bottom
385	jk = jpk
386	e3t_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,jk) - gdepw_0(:,:,jk))
387
388	!==============================
389	! Computing mbathy
390	!==============================
391
392	IF( lwp ) WRITE(numout,*) ''
393	IF( lwp ) WRITE(numout,*) 'MES mbathy:'
394	IF( lwp ) WRITE(numout,*) ''
395
396	DO jj = 1, jpj
397	DO ji = 1, jpi
398	DO jk = 1, jpkm1
399	!IF( lwp ) WRITE(numout,*) 'scobot: ', scobot(ji,jj), 'gdept_0: ', gdept_0(ji,jj,jk)
400	IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
401	IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0
402	IF( scobot(ji,jj) < 0.e0 ) mbathy(ji,jj) = int( scobot(ji,jj)) !???? do we need it?
403	END DO
404	!IF( lwp ) WRITE(numout,*) 'mbathy(',ji,',',jj,') = ', mbathy(ji,jj)
405	END DO
406	END DO
407
408	!==============================================
409	! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0
410	!==============================================
411
412	DO jj = 1, jpjm1
413	DO ji = 1, jpim1
414	DO jk = 1, jpk
415
416	e3u_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + &
417	MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) ) / &
418	MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)) )
419
420	e3v_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + &
421	MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) / &
422	MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) )
423
424	e3uw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk) + &
425	MIN(1,mbathy(ji+1,jj))*e3w_0(ji+1,jj,jk) ) / &
426	MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)) )
427
428	e3vw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk) + &
429	MIN(1,mbathy(ji,jj+1))*e3w_0(ji,jj+1,jk) ) / &
430	MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) )
431
432	e3f_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + &
433	MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) + &
434	MIN(1,mbathy(ji+1,jj+1))* e3t_0(ji+1,jj+1,jk) + &
435	MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) / &
436	MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) &
437	+ MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1)) )
438	ENDDO
439	ENDDO
440	ENDDO
441
442	! Adjusting for geopotential levels, if any
443
444	DO je = 1, tot_env ! LOOP over all the requested envelopes
445
446	IF (je == 1) THEN
447	s_1st = 1
448	num_s = nn_slev(je)
449	max_env_up = 0._wp ! surface
450	min_env_up = 0._wp ! surface
451	ELSE
452	s_1st = s_1st + num_s - 1
453	num_s = nn_slev(je) + 1
454	max_env_up = max_dep_env(je-1)
455	min_env_up = min_dep_env(je-1)
456	ENDIF
457
458	max_env_dw = max_dep_env(je)
459	min_env_dw = min_dep_env(je)
460
461	IF ( max_env_up == min_env_up .AND. max_env_dw == min_env_dw ) THEN
462
463	DO jj = 1,jpjm1
464	DO ji = 1,jpim1
465	DO jk = 1, num_s
466
467	ind = s_1st+jk-1
468	e3u_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji+1,jj,ind))
469	e3uw_0(ji,jj,ind) = MIN( e3w_0(ji,jj,ind), e3w_0(ji+1,jj,ind))
470	e3v_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji,jj+1,ind))
471	e3vw_0(ji,jj,ind) = MIN( e3w_0(ji,jj,ind), e3w_0(ji,jj+1,ind))
472	e3f_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji+1,jj,ind))
473
474	ENDDO
475	ENDDO
476	ENDDO
477
478	ENDIF
479
480	ENDDO
481
482	!==============================
483	! Computing gdept3w
484	!==============================
485
486	gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1)
487	DO jk = 2, jpk
488	gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
489	END DO
490
491	!=================================================================================
492	! From here equal to sco code - domzgr.F90 line 2079
493	! Lateral B.C.
494
495	CALL lbc_lnk( e3t_0 , 'T', 1._wp )
496	CALL lbc_lnk( e3u_0 , 'U', 1._wp )
497	CALL lbc_lnk( e3v_0 , 'V', 1._wp )
498	CALL lbc_lnk( e3f_0 , 'F', 1._wp )
499	CALL lbc_lnk( e3w_0 , 'W', 1._wp )
500	CALL lbc_lnk( e3uw_0, 'U', 1._wp )
501	CALL lbc_lnk( e3vw_0, 'V', 1._wp )
502
503	where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0
504	where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0
505	where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0
506	where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0
507	where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0
508	where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0
509	where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0
510
511	IF( lwp ) WRITE(numout,*) 'Refine mbathy'
512	DO jj = 1, jpj
513	DO ji = 1, jpi
514	DO jk = 1, jpkm1
515	IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )
516	END DO
517	IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0
518	END DO
519	END DO
520
521	IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), &
522	& ' MAX ', MAXVAL( mbathy(:,:) )
523
524	IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain
525	WRITE(numout,*) 'zgr_sco : mbathy'
526	WRITE(numout,*) '~~~~~~~'
527	WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )
528	WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
529	& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) )
530	WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), &
531	& ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), &
532	& ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), &
533	& ' w ', MINVAL( e3w_0 (:,:,:) )
534
535	WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
536	& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) )
537	WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), &
538	& ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), &
539	& ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), &
540	& ' w ', MAXVAL( e3w_0 (:,:,:) )
541	ENDIF
542
543	!================================================================================
544	! Check coordinates makes sense
545	!================================================================================
546
547	! Extracting MEs depth profile in the shallowest point of the deepest
548	! envelope for a first check of monotonicity of transformation
549	! (also useful for debugging)
550	pmsk(:,:) = 0.0
551	WHERE ( bathy > 0.0 ) pmsk = 1.0
552	CALL mpp_minloc( envlt(:,:,tot_env), pmsk(:,:), mindep, eii, eij )
553	IF ((mi0(eii)>=1 .AND. mi0(eii)<=jpi) .AND. (mj0(eij)>=1 .AND. mj0(eij)<=jpj)) THEN
554	irnk = mpprank
555	DO je = 1, tot_env
556	rn_ebot_max(je) = envlt(mi0(eii),mj0(eij),je)
557	END DO
558	z_gsigt1(:) = z_gsigt3(mi0(eii),mj0(eij),:)
559	z_gsigw1(:) = z_gsigw3(mi0(eii),mj0(eij),:)
560	gdept1(:) = gdept_0(mi0(eii),mj0(eij),:)
561	gdepw1(:) = gdepw_0(mi0(eii),mj0(eij),:)
562	e3t1(:) = e3t_0(mi0(eii),mj0(eij),:)
563	e3w1(:) = e3w_0(mi0(eii),mj0(eij),:)
564	ELSE
565	irnk = -1
566	END IF
567	IF( lk_mpp ) CALL mppsync
568	IF( lk_mpp ) CALL mpp_max(irnk)
569	IF( lk_mpp ) CALL mppbcast_a_real(rn_ebot_max, max_nn_env, irnk )
570	IF( lk_mpp ) CALL mppbcast_a_real(z_gsigt1 , jpk , irnk )
571	IF( lk_mpp ) CALL mppbcast_a_real(z_gsigw1 , jpk , irnk )
572	IF( lk_mpp ) CALL mppbcast_a_real(gdept1 , jpk , irnk )
573	IF( lk_mpp ) CALL mppbcast_a_real(gdepw1 , jpk , irnk )
574	IF( lk_mpp ) CALL mppbcast_a_real(e3t1 , jpk , irnk )
575	IF( lk_mpp ) CALL mppbcast_a_real(e3w1 , jpk , irnk )
576
577	IF( lwp ) THEN
578	WRITE(numout,*) ""
579	WRITE(numout,*) "mes_build:"
580	WRITE(numout,*) "~~~~~~~~~"
581	WRITE(numout,*) ""
582	WRITE(numout,*) " FIRST CHECK: Checking MEs-levels profile in the shallowest point of the last envelope:"
583	WRITE(numout,*) " it is the most likely point where monotonicty of splines may be violeted. "
584	WRITE(numout,*) ""
585	DO je = 1, tot_env
586	WRITE(numout,) ' depth of envelope ', je, ' at point (',eii,',',eij,') is ', rn_ebot_max(je)
587	END DO
588	WRITE(numout,*) ""
589	WRITE(numout,*) " MEs-coordinates"
590	WRITE(numout,*) ""
591	WRITE(numout,*) " k z_gsigw1 z_gsigt1"
592	WRITE(numout,*) ""
593	DO jk = 1, jpk
594	WRITE(numout,*) ' ', jk, z_gsigw1(jk), z_gsigt1(jk)
595	ENDDO
596	WRITE(numout,*) ""
597	WRITE(numout,*) "-----------------------------------------------------------------"
598	WRITE(numout,*) ""
599	WRITE(numout,*) " MEs-levels depths and scale factors"
600	WRITE(numout,*) ""
601	WRITE(numout,*) " k gdepw1 e3w1 gdept1 e3t1"
602	WRITE(numout,*) ""
603	DO jk = 1, jpk
604	WRITE(numout,*) ' ', jk, gdepw1(jk), e3w1(jk), gdept1(jk), e3t1(jk)
605	ENDDO
606	WRITE(numout,*) "-----------------------------------------------------------------"
607	ENDIF
608
609	! Checking monotonicity for cubic splines
610	DO jk = 1, jpk-1
611	IF ( gdept1(jk+1) < gdept1(jk) ) THEN
612	WRITE(ctlmes,*) 'NOT MONOTONIC gdept_0: change envelopes'
613	CALL ctl_stop( ctlmes )
614	END IF
615	IF ( gdepw1(jk+1) < gdepw1(jk) ) THEN
616	WRITE(ctlmes,*) 'NOT MONOTONIC gdepw_0: change envelopes'
617	CALL ctl_stop( ctlmes )
618	END IF
619	IF ( e3t1(jk) < 0.0 ) THEN
620	WRITE(ctlmes,*) 'NEGATIVE e3t_0: change envelopes'
621	CALL ctl_stop( ctlmes )
622	END IF
623	IF ( e3w1(jk) < 0.0 ) THEN
624	WRITE(ctlmes,*) 'NEGATIVE e3w_0: change envelopes'
625	CALL ctl_stop( ctlmes )
626	END IF
627	END DO
628
629	! CHECKING THE WHOLE DOMAIN
630	DO ji = 1, jpi
631	DO jj = 1, jpj
632	DO jk = 1, jpk !mbathy(ji,jj)
633	! check coordinate is monotonically increasing
634	IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN
635	WRITE(ctmp1,*) 'ERROR mes_build: e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
636	WRITE(numout,*) 'ERROR mes_build: e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk
637	WRITE(numout,*) 'e3w',e3w_0(ji,jj,:)
638	WRITE(numout,*) 'e3t',e3t_0(ji,jj,:)
639	CALL ctl_stop( ctmp1 )
640	ENDIF
641	! and check it has never gone negative
642	IF ( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN
643	WRITE(ctmp1,*) 'ERROR mes_build: gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
644	WRITE(numout,*) 'ERROR mes_build: gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk
645	WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:)
646	WRITE(numout,*) 'gdept',gdept_0(ji,jj,:)
647	CALL ctl_stop( ctmp1 )
648	ENDIF
649	END DO
650	END DO
651	END DO
652	!
653	CALL wrk_dealloc( jpi, jpj, max_nn_env , envlt )
654	!
655	IF( nn_timing == 1 ) CALL timing_stop('mes_build')
656
657	END SUBROUTINE mes_build
658
659	! =====================================================================================================
660
661	SUBROUTINE mes_init
662	!!-----------------------------------------------------------------------------
663	!! * ROUTINE mes_init *
664	!!
665	!! ** Purpose : Initilaise a Multi Enveloped s-coordinate (MEs) system
666	!!
667	!!-----------------------------------------------------------------------------
668
669	CHARACTER(lc) :: env_name ! name of the externally defined envelope
670	INTEGER :: ji, jj, je, inum
671	INTEGER :: cor_env, ios
672	REAL(wp) :: rn_bot_min ! min depth of ocean bottom (>0) (m)
673	REAL(wp) :: rn_bot_max ! max depth of ocean bottom (= ocean depth) (>0) (m)
674
675	NAMELIST/namzgr_mes/rn_bot_min , rn_bot_max , ln_envl , &
676	& nn_strt , nn_slev , rn_e_hc , &
677	& rn_e_th , rn_e_bb , rn_e_ba , &
678	& rn_e_al , ln_pst_mes , ln_pst_l2g, &
679	& rn_e3pst_min, rn_e3pst_rat
680	!!----------------------------------------------------------------------
681	!
682	IF( nn_timing == 1 ) CALL timing_start('mes_init')
683	!
684	CALL wrk_alloc( jpi, jpj, max_nn_env , envlt )
685	!
686	! Initialising some arrays
687	max_dep_env(:) = 0.0
688	min_dep_env(:) = 0.0
689	!
690	IF(lwp) THEN ! control print
691	WRITE(numout,*) ''
692	WRITE(numout,*) 'mes_init: Setting a Multi-Envelope s-ccordinate system (Bruciaferri et al. 2018)'
693	WRITE(numout,*) '~~~~~~~~'
694	END IF
695	!
696	! Namelist namzgr_mes in reference namelist : envelopes and
697	! sigma-stretching parameters
698	REWIND( numnam_ref )
699	READ ( numnam_ref, namzgr_mes, IOSTAT = ios, ERR = 901)
700	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_mes in reference namelist', lwp )
701	! Namelist namzgr_mes in configuration namelist : envelopes and
702	! sigma-stretching parameters
703	REWIND( numnam_cfg )
704	READ ( numnam_cfg, namzgr_mes, IOSTAT = ios, ERR = 902 )
705	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_mes in configuration namelist', lwp )
706	IF(lwm) WRITE ( numond, namzgr_mes )
707	!
708	! 1) Check if namelist is defined correctly.
709	! Not strictly needed but we force the user
710	! to define the namelist correctly.
711	tot_env = 0 ! total number of requested envelopes
712	cor_env = 0 ! total number of correctly defined envelopes
713	DO je = 1, max_nn_env
714	IF ( ln_envl(je) ) tot_env = tot_env + 1
715	IF ( ln_envl(je) .AND. cor_env == (je-1)) cor_env = cor_env + 1
716	ENDDO
717	WRITE(ctlmes,*) 'num. of REQUESTED env. and num. of CORRECTLY defined env. are DIFFERENT'
718	IF ( tot_env /= cor_env ) CALL ctl_stop( ctlmes )
719	IF(lwp) THEN
720	WRITE(numout,*) ' Num. or requested envelopes: ', tot_env
721	WRITE(numout,*) ''
722	END IF
723	!
724	! 2) Checking consistency of user defined parameters
725	WRITE(ctlmes,*) 'TOT number of levels (jpk) IS DIFFERENT from sum over nn_slev(:)'
726	IF ( SUM(nn_slev(:)) /= jpk ) CALL ctl_stop( ctlmes )
727	!
728	IF( .NOT.(ln_loc_zgr) .AND. ln_pst_l2g) ln_pst_l2g = .FALSE.
729	!
730	! 3) Reading Bathymetry and envelopes
731	IF( ntopo == 1 ) THEN
732	IF(lwp) THEN
733	WRITE(numout,*) ' Reading bathymetry and envelopes ...'
734	WRITE(numout,*) ''
735	END IF
736
737	CALL iom_open ( 'bathy_meter.nc', inum )
738	CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr )
739	DO je = 1, tot_env
740	WRITE (env_name, '(A6,I0)') 'hbatt_', je
741	CALL iom_get ( inum, jpdom_data, TRIM(env_name) , envlt(:,:,je), lrowattr=ln_use_jattr )
742	ENDDO
743	CALL iom_close( inum )
744	ELSE
745	WRITE(ctlmes,*) 'parameter , ntopo = ', ntopo
746	CALL ctl_stop( ctlmes )
747	ENDIF
748	IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==!
749	!
750	! 4) Checking consistency of envelopes
751	DO je = 1, tot_env-1
752	WRITE(ctlmes,*) 'Envelope ', je+1, ' is shallower that Envelope ', je
753	IF (MAXVAL(envlt(:,:,je+1)) < MAXVAL(envlt(:,:,je))) CALL ctl_stop( ctlmes )
754	ENDDO
755	! 5) Computing max and min depths of envelopes
756	DO je = 1, tot_env
757	max_dep_env(je) = MAXVAL(envlt(:,:,je))
758	min_dep_env(je) = MINVAL(envlt(:,:,je))
759	IF( lk_mpp ) CALL mpp_max( max_dep_env(je) )
760	IF( lk_mpp ) CALL mpp_min( min_dep_env(je) )
761	END DO
762	IF( lk_mpp ) CALL mppsync
763	!
764	! 6) Set maximum and minimum ocean depth
765	bathy(:,:) = MIN( rn_bot_max, bathy(:,:) )
766	DO jj = 1, jpj
767	DO ji = 1, jpi
768	IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_bot_min, bathy(ji,jj) )
769	END DO
770	END DO
771	!
772	IF(lwp) THEN ! control print
773	WRITE(numout,*)
774	WRITE(numout,*) ' GEOMETRICAL FEATURES OF THE REQUESTED MEs GRID:'
775	WRITE(numout,*) ' -----------------------------------------------'
776	WRITE(numout,*) ' Minimum depth of the ocean rn_bot_min = ', rn_bot_min
777	WRITE(numout,*) ' Maximum depth of the ocean rn_bot_max = ', rn_bot_max
778	WRITE(numout,*) ''
779	WRITE(numout,*) ' ENVELOPES:'
780	WRITE(numout,*) ' ========= '
781	WRITE(numout,*) ''
782	DO je = 1, tot_env
783	WRITE(numout,) ' envelope ', je, ': min depth = ', min_dep_env(je)
784	WRITE(numout,*) ' max depth = ', max_dep_env(je)
785	WRITE(numout,*) ''
786	END DO
787	WRITE(numout,*) ''
788	WRITE(numout,*) ' SUBDOMAINS:'
789	WRITE(numout,*) ' ========== '
790	WRITE(numout,*) ''
791	DO je = 1, tot_env
792	WRITE(numout,) ' subdomain ', je,':'
793	IF ( je == 1) THEN
794	WRITE(numout,*) ' Shallower envelope: free surface'
795	ELSE
796	WRITE(numout,*) ' Shallower envelope: envelope ',je-1
797	END IF
798	WRITE(numout,*) ' Deeper envelope: envelope ',je
799	WRITE(numout,*) ' Number of MEs-levs: nn_slev(',je,') = ',nn_slev(je)
800	IF ( isodd(je) ) THEN
801	WRITE(numout,*) ' Stretched s-coordinates: '
802	ELSE
803	WRITE(numout,*) ' Stretched CUBIC SPLINES: '
804	END IF
805	IF (nn_strt(je) == 0) WRITE(numout,*) ' M96 stretching function'
806	IF (nn_strt(je) == 1) WRITE(numout,*) ' SH94 stretching function'
807	IF (nn_strt(je) == 2) WRITE(numout,*) ' SF12 stretching function'
808	WRITE(numout,*) ' critical depth rn_e_hc(',je,') = ',rn_e_hc(je)
809	WRITE(numout,*) ' surface stretc. coef. rn_e_th(',je,') = ',rn_e_th(je)
810	IF (nn_strt(je) == 2) THEN
811	WRITE(numout,*) ' bottom stretc. coef. rn_e_ba(',je,') = ',rn_e_ba(je)
812	END IF
813	WRITE(numout,*) ' bottom stretc. coef. rn_e_bb(',je,') = ',rn_e_bb(je)
814	IF (nn_strt(je) == 2) THEN
815	WRITE(numout,*) ' bottom stretc. coef. rn_e_al(',je,') = ',rn_e_al(je)
816	END IF
817	WRITE(numout,*) ' ------------------------------------------------------------------'
818	ENDDO
819	IF( ln_loc_zgr) THEN
820	WRITE(numout,*) ''
821	WRITE(numout,*) ' LOCALISED MEs '
822	WRITE(numout,*) ' ============= '
823	ENDIF
824	IF( ln_pst_mes .OR. ln_pst_l2g) THEN
825	WRITE(numout,*) ''
826	WRITE(numout,*) ' APPLYING PARTIAL-STEPS to '
827	WRITE(numout,*) ' ========================= '
828	WRITE(numout,*) ''
829	WRITE(numout,*) ' a) MEs area : ', ln_pst_mes
830	IF( ln_loc_zgr) WRITE(numout,*) ' b) Transition area: ', ln_pst_l2g
831	WRITE(numout,*) ''
832	WRITE(numout,*) ' Partial step thickness is set larger than the minimum of'
833	WRITE(numout,*) ''
834	WRITE(numout,*) ' rn_e3pst_min = ', rn_e3pst_min
835	WRITE(numout,*) ' and'
836	WRITE(numout,*) ' rn_e3pst_rat = ', rn_e3pst_rat
837	WRITE(numout,*) ''
838	WRITE(numout,*) ' with 0 < rn_e3pst_rat < 1'
839	ENDIF
840	ENDIF
841
842	END SUBROUTINE mes_init
843
844	! =====================================================================================================
845
846	FUNCTION isodd( a ) RESULT(odd)
847	! Notice that this method uses a btest intrinsic function. If you look at any number in binary
848	! notation, you'll note that the Least Significant Bit (LSB) will always be set for odd numbers
849	! and cleared for even numbers. Hence, all that is necessary to determine if a number is even or
850	! odd is to check the LSB. Since bitwise operations like AND, OR, XOR etc. are usually much faster
851	! than the arithmetic operations, this method will run a lot faster than the one with modulo
852	! operation.
853	! http://forums.devshed.com/software-design-43/quick-algorithm-determine-odd-29843.html
854
855	INTEGER, INTENT (in) :: a
856	LOGICAL :: odd
857	odd = btest(a, 0)
858	END FUNCTION isodd
859
860	! =====================================================================================================
861
862	FUNCTION iseven( a ) RESULT(even)
863	INTEGER, INTENT (in) :: a
864	LOGICAL :: even
865	even = .NOT. isodd(a)
866	END FUNCTION iseven
867
868	! =====================================================================================================
869
870	FUNCTION sech( x ) RESULT( seh )
871
872	REAL(wp), INTENT(in ) :: x
873	REAL(wp) :: seh
874
875	seh = 1._wp / COSH(x)
876
877	END FUNCTION sech
878
879	! =====================================================================================================
880
881	FUNCTION sigma( vgrid, pk, kmax ) RESULT( ps1 )
882	!!----------------------------------------------------------------------
883	!! * ROUTINE sigma *
884	!!
885	!! ** Purpose : provide the analytical function for sigma-coordinate
886	!! (not stretched s-coordinate).
887	!!
888	!! ** Method : the function provide the non-dimensional position of
889	!! T and W points (i.e. between 0 and 1).
890	!! T-points at integer values (between 1 and jpk)
891	!! W-points at integer values - 1/2 (between 0.5 and
892	!! jpk-0.5)
893	!!
894	!!----------------------------------------------------------------------
895	INTEGER, INTENT (in) :: pk ! continuous k coordinate
896	INTEGER, INTENT (in) :: kmax ! number of levels
897	INTEGER, INTENT (in) :: vgrid ! type of vertical grid: 0 = T, 1 = W
898	REAL(wp) :: kindx ! index of T or W points
899	REAL(wp) :: ps1 ! value of sigma coordinate (-1 <= ps1 <=0)
900	!
901	SELECT CASE (vgrid)
902
903	CASE (0) ! T-points at integer values (between 1 and jpk)
904	kindx = REAL(pk,wp)
905	CASE (1) ! W-points at integer values - 1/2 (between 0.5 and jpk-0.5)
906	kindx = REAL(pk,wp) - 0.5_wp
907	CASE DEFAULT
908	WRITE(ctlmes,*) 'No valid vertical grid option selected'
909
910	END SELECT
911
912
913	ps1 = -(kindx - 0.5) / REAL(kmax-1) ! The surface is at the first W level
914	! while the bottom coincides with the
915	! last W level
916	END FUNCTION sigma
917
918	! =====================================================================================================
919
920	FUNCTION s_coord( s, ca, cb, alpha, kmax, ftype ) RESULT ( pf1 )
921	!!----------------------------------------------------------------------
922	!! * ROUTINE s_coord *
923	!!
924	!! ** Purpose : provide the analytical stretching function
925	!! for s-coordinate.
926	!!
927	!! ** Method : if ftype = 2: Siddorn and Furner 2012
928	!! analytical function is used
929	!! if ftype = 1: Song and Haidvogel 1994
930	!! analytical function is used
931	!! if ftype = 0: Madec et al. 1996
932	!! analytical function is used
933	!! (pag 65 of NEMO Manual)
934	!! Reference : Madec, Lott, Delecluse and
935	!! Crepon, 1996. JPO, 26, 1393-1408
936	!!
937	!! s MUST be NEGATIVE: -1 <= s <= 0
938	!!----------------------------------------------------------------------
939	REAL(wp), INTENT(in) :: s ! continuous not stretched
940	! "s" coordinate
941	REAL(wp), INTENT(in) :: ca ! surface stretch. coeff
942	REAL(wp), INTENT(in) :: cb ! bottom stretch. coeff
943	REAL(wp), INTENT(in) :: alpha ! alpha stretch. coeff for SF12
944	INTEGER, INTENT (in) :: kmax ! number of levels
945	INTEGER, INTENT(in) :: ftype ! type of stretching function
946	REAL(wp) :: pf1 ! "s" stretched value
947	! SF12 stretch. funct. parameters
948	REAL(wp) :: psmth ! smoothing parameter for transition
949	! between shallow and deep areas
950	REAL(wp) :: za1,za2,za3 ! local variables
951	REAL(wp) :: zn1,zn2,ps ! local variables
952	REAL(wp) :: za,zb,zx ! local variables
953	!!----------------------------------------------------------------------
954	SELECT CASE (ftype)
955
956	CASE (0) ! M96 stretching function
957	pf1 = ( TANH( ca * ( s + cb ) ) - TANH( cb * ca ) ) &
958	& * ( COSH( ca ) &
959	& + COSH( ca * ( 2.e0 * cb - 1.e0 ) ) ) &
960	& / ( 2._wp * SINH( ca ) )
961	CASE (1) ! SH94 stretching function
962	IF ( ca == 0 ) THEN ! uniform sigma
963	pf1 = s
964	ELSE ! stretched sigma
965	pf1 = (1._wp - cb) * SINH(ca*s) / SINH(ca) + &
966	& cb * ( ( TANH(ca(s + 0.5_wp)) - TANH(0.5_wpca) ) / &
967	& (2._wpTANH(0.5_wpca)) )
968	END IF
969	CASE (2) ! SF12 stretching function
970	psmth = 1.0_wp ! We consider only the case for efold = 0
971	ps = -s
972	zn1 = 1. / REAL(kmax-1)
973	zn2 = 1. - zn1
974
975	za1 = (alpha+2.0_wp)zn1(alpha+1.0_wp)-(alpha+1.0_wp)zn1**(alpha+2.0_wp)
976	za2 = (alpha+2.0_wp)zn2(alpha+1.0_wp)-(alpha+1.0_wp)zn2**(alpha+2.0_wp)
977	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
978
979	za = cb - za3*(ca-za1)-za2
980	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp)) )
981	zb = (ca - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
982	zx = 1.0_wp-za/2.0_wp-zb
983
984	pf1 = za(ps(1.0_wp-ps/2.0_wp))+zbps*3.0_wp + &
985	& zx( (alpha+2.0_wp)ps**(alpha+1.0_wp) - &
986	& (alpha+1.0_wp)ps*(alpha+2.0_wp) )
987	pf1 = -pf1psmth+ps(1.0_wp-psmth)
988	CASE (3) ! stretching function for deeper part of the domain
989	IF ( ca == 0 ) THEN ! uniform sigma
990	pf1 = s
991	ELSE ! stretched sigma
992	pf1 = (1._wp - cb) * SINH(ca*s) / SINH(ca)
993	END IF
994	END SELECT
995
996	END FUNCTION s_coord
997
998	! =====================================================================================================
999
1000	FUNCTION D1s_coord( s, ca, cb, alpha, kmax, ftype) RESULT( pf1 )
1001	!!----------------------------------------------------------------------
1002	!! * ROUTINE D1s_coord *
1003	!!
1004	!! ** Purpose : provide the 1st derivative of the analytical
1005	!! stretching function for s-coordinate.
1006	!!
1007	!! ** Method : if ftype = 2: Siddorn and Furner 2012
1008	!! analytical function is used
1009	!! if ftype = 1: Song and Haidvogel 1994
1010	!! analytical function is used
1011	!! if ftype = 0: Madec et al. 1996
1012	!! analytical function is used
1013	!! (pag 65 of NEMO Manual)
1014	!! Reference : Madec, Lott, Delecluse and
1015	!! Crepon, 1996. JPO, 26, 1393-1408
1016	!!
1017	!! s MUST be NEGATIVE: -1 <= s <= 0
1018	!!----------------------------------------------------------------------
1019	REAL(wp), INTENT(in ) :: s ! continuous not stretched
1020	! "s" coordinate
1021	REAL(wp), INTENT(in) :: ca ! surface stretch. coeff
1022	REAL(wp), INTENT(in) :: cb ! bottom stretch. coeff
1023	REAL(wp), INTENT(in) :: alpha ! alpha stretch. coeff for SF12
1024	INTEGER, INTENT (in) :: kmax ! number of levels
1025	INTEGER, INTENT(in ) :: ftype ! type of stretching function
1026	REAL(wp) :: pf1 ! "s" stretched value
1027	! SF12 stretch. funct. parameters
1028	REAL(wp) :: psmth ! smoothing parameter for transition
1029	! between shallow and deep areas
1030	REAL(wp) :: za1,za2,za3 ! local variables
1031	REAL(wp) :: zn1,zn2,ps ! local variables
1032	REAL(wp) :: za,zb,zx ! local variables
1033	REAL(wp) :: zt1,zt2,zt3 ! local variables
1034	!!----------------------------------------------------------------------
1035	SELECT CASE (ftype)
1036
1037	CASE (0) ! M96 stretching function
1038	pf1 = ( ca * ( COSH(ca(2._wp cb - 1._wp)) + COSH(ca)) * &
1039	sech(ca * (s+cb))*2 ) / (2._wp SINH(ca))
1040	CASE (1) ! SH94 stretching function
1041	IF ( ca == 0 ) then ! uniform sigma
1042	pf1 = 1._wp / REAL(kmax,wp)
1043	ELSE ! stretched sigma
1044	pf1 = (1._wp - cb) * ca * COSH(cas) / (SINH(ca) REAL(kmax,wp)) + &
1045	cb * ca * &
1046	(sech((s+0.5_wp)ca)2) / (2._wp TANH(0.5_wpca) REAL(kmax,wp))
1047	END IF
1048	CASE (2) ! SF12 stretching function
1049	ps = -s
1050	zn1 = 1. / REAL(kmax-1)
1051	zn2 = 1. - zn1
1052
1053	za1 = (alpha+2.0_wp)zn1(alpha+1.0_wp)-(alpha+1.0_wp)zn1**(alpha+2.0_wp)
1054	za2 = (alpha+2.0_wp)zn2(alpha+1.0_wp)-(alpha+1.0_wp)zn2**(alpha+2.0_wp)
1055	za3 = (zn23.0_wp - za2)/( zn13.0_wp - za1)
1056
1057	za = cb - za3*(ca-za1)-za2
1058	za = za/( zn2-0.5_wp(za2+zn22.0_wp) - za3(zn1-0.5_wp(za1+zn1*2.0_wp)) )
1059	zb = (ca - za1 - za( zn1-0.5_wp(za1+zn12.0_wp ) ) ) / (zn13.0_wp - za1)
1060	zx = (alpha+2.0_wp)*(alpha+1.0_wp)
1061
1062	zt1 = 0.5_wpza(ps-1._wp)((ps2.0_wp + 3._wpps + 2._wp)ps*alpha - 2._wp)
1063	zt2 = zb(zxps*(alpha+1.0_wp) - zxps*(alpha) + 3._wpps**2._wp)
1064	zt3 = zx(ps-1._wp)ps**(alpha)
1065
1066	pf1 = zt1 + zt2 - zt3
1067
1068	END SELECT
1069
1070	END FUNCTION D1s_coord
1071
1072	! =====================================================================================================
1073
1074	FUNCTION z_mes(dep_up, dep_dw, Cs, s, hc) RESULT( z )
1075	!!----------------------------------------------------------------------
1076	!! * ROUTINE z_mes *
1077	!!
1078	!! ** Purpose : provide the analytical trasformation from the
1079	!! computational space (mes-coordinate) to the
1080	!! physical space (depth z)
1081	!!
1082	!! N.B.: z is downward positive defined as well as envelope surfaces.
1083	!! Therefore, Cs MUST be positive, meaning 0. <= Cs <= 1.
1084	!!----------------------------------------------------------------------
1085	REAL(wp), INTENT(in ) :: dep_up ! shallower envelope
1086	REAL(wp), INTENT(in ) :: dep_dw ! deeper envelope
1087	REAL(wp), INTENT(in ) :: Cs ! stretched s coordinate
1088	REAL(wp), INTENT(in ) :: s ! not stretched s coordinate
1089	REAL(wp), INTENT(in ) :: hc ! critical depth
1090	REAL :: z ! downward positive depth
1091	!!----------------------------------------------------------------------
1092	!
1093
1094	z = dep_up + hc * s + Cs * (dep_dw - hc - dep_up)
1095
1096	END FUNCTION z_mes
1097
1098	! =====================================================================================================
1099
1100	FUNCTION D1z_mes(dep_up, dep_dw, d1Cs, hc, kmax) RESULT( d1z )
1101	!!----------------------------------------------------------------------
1102	!! * ROUTINE D1z_mes *
1103	!!
1104	!! ** Purpose : provide the 1st derivative of the analytical
1105	!! trasformation from the computational space
1106	!! (mes-coordinate) to the physical space (depth z)
1107	!!
1108	!! N.B.: z is downward positive defined as well as envelope surfaces.
1109	!! Therefore, Cs MUST be positive, meaning 0. <= Cs <= 1.
1110	!!----------------------------------------------------------------------
1111	REAL(wp), INTENT(in ) :: dep_up ! shallower envelope
1112	REAL(wp), INTENT(in ) :: dep_dw ! deeper envelope
1113	REAL(wp), INTENT(in ) :: d1Cs ! 1st derivative of the
1114	! stretched s coordinate
1115	REAL(wp), INTENT(in ) :: hc ! critical depth
1116	INTEGER, INTENT(in ) :: kmax ! number of levels
1117	REAL(wp) :: d1z ! 1st derivative of downward
1118	! positive depth
1119	!!----------------------------------------------------------------------
1120	!
1121
1122	d1z = hc / REAL(kmax,wp) + d1Cs * (dep_dw - hc - dep_up)
1123	!d1z = hc + d1Cs * (dep_dw - hc - dep_up)
1124
1125	END FUNCTION D1z_mes
1126
1127	! =====================================================================================================
1128
1129	FUNCTION cub_spl(tau, d, n, ibcbeg, ibcend) RESULT ( c )
1130	!!----------------------------------------------------------------------
1131	!!
1132	!! CUBSPL defines an interpolatory cubic spline.
1133	!!
1134	!! Discussion:
1135	!!
1136	!! A tridiagonal linear system for the unknown slopes S(I) of
1137	!! F at TAU(I), I=1,..., N, is generated and then solved by Gauss
1138	!! elimination, with S(I) ending up in C(2,I), for all I.
1139	!!
1140	!! Author: Carl de Boor
1141	!!
1142	!! Reference: Carl de Boor, Practical Guide to Splines,
1143	!! Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27.
1144	!!
1145	!! Parameters:
1146	!!
1147	!! Input:
1148	!! TAU(N) : the abscissas or X values of the data points.
1149	!! The entries of TAU are assumed to be strictly
1150	!! increasing.
1151	!!
1152	!! N : the number of data points. N is assumed to be
1153	!! at least 2.
1154	!!
1155	!! IBCBEG : boundary condition indicator at TAU(1).
1156	!! = 0 no boundary condition at TAU(1) is given.
1157	!! In this case, the "not-a-knot condition"
1158	!! is used.
1159	!! = 1 the 1st derivative at TAU(1) is equal to
1160	!! the input value D(2,1).
1161	!! = 2 the 2nd derivative at TAU(1) is equal to
1162	!! the input value D(2,1).
1163	!!
1164	!! IBCEND : boundary condition indicator at TAU(N).
1165	!! = 0 no boundary condition at TAU(N) is given.
1166	!! In this case, the "not-a-knot condition"
1167	!! is used.
1168	!! = 1 the 1st derivative at TAU(N) is equal to
1169	!! the input value D(2,2).
1170	!! = 2 the 2nd derivative at TAU(N) is equal to
1171	!! the input value D(2,2).
1172	!!
1173	!! D(1,1): value of the function at TAU(1)
1174	!! D(1,1): value of the function at TAU(N)
1175	!! D(2,1): if IBCBEG is 1 (2) it is the value of the
1176	!! 1st (2nd) derivative at TAU(1)
1177	!! D(2,2): if IBCBEG is 1 (2) it is the value of the
1178	!! 1st (2nd) derivative at TAU(N)
1179	!!
1180	!! Output:
1181	!! C(4,N) : contains the polynomial coefficients of the
1182	!! cubic interpolating spline.
1183	!! In the interval interval (TAU(I), TAU(I+1)),
1184	!! the spline F is given by
1185	!!
1186	!! F(X) =
1187	!! C(1,I) +
1188	!! C(2,I) * H +
1189	!! C(3,I) * H^2 / 2 +
1190	!! C(4,I) * H^3 / 6.
1191	!!
1192	!! where H = X - TAU(I).
1193	!!
1194	!!----------------------------------------------------------------------
1195	INTEGER, INTENT(in ) :: n
1196	REAL(wp), INTENT(in ) :: tau(n)
1197	REAL(wp), INTENT(in ) :: d(2,2)
1198	INTEGER, INTENT(in ) :: ibcbeg, ibcend
1199	REAL(wp) :: c(4,n)
1200	REAL(wp) :: divdf1
1201	REAL(wp) :: divdf3
1202	REAL(wp) :: dtau
1203	REAL(wp) :: g
1204	INTEGER(wp) :: i
1205	!!----------------------------------------------------------------------
1206	! Initialise c and copy d values to c
1207	c(:,:) = 0.0D+00
1208	c(1,1) = d(1,1)
1209	c(1,n) = d(1,2)
1210	c(2,1) = d(2,1)
1211	c(2,n) = d(2,2)
1212	!
1213	! C(3,) and C(4,) are used initially for temporary storage.
1214	! Store first differences of the TAU sequence in C(3,*).
1215	! Store first divided difference of data in C(4,*).
1216	DO i = 2, n
1217	c(3,i) = tau(i) - tau(i-1)
1218	c(4,i) = ( c(1,i) - c(1,i-1) ) / c(3,i)
1219	END DO
1220
1221	! Construct the first equation from the boundary condition
1222	! at the left endpoint, of the form:
1223	!
1224	! C(4,1) * S(1) + C(3,1) * S(2) = C(2,1)
1225	!
1226	! IBCBEG = 0: Not-a-knot
1227	IF ( ibcbeg == 0 ) THEN
1228
1229	IF ( n <= 2 ) THEN
1230	c(4,1) = 1._wp
1231	c(3,1) = 1._wp
1232	c(2,1) = 2._wp * c(4,2)
1233	ELSE
1234	c(4,1) = c(3,3)
1235	c(3,1) = c(3,2) + c(3,3)
1236	c(2,1) = ( ( c(3,2) + 2._wp * c(3,1) ) * c(4,2) &
1237	* c(3,3) + c(3,2)*2 c(4,3) ) / c(3,1)
1238	END IF
1239
1240	! IBCBEG = 1: derivative specified.
1241	ELSE IF ( ibcbeg == 1 ) then
1242
1243	c(4,1) = 1._wp
1244	c(3,1) = 0._wp
1245
1246	! IBCBEG = 2: Second derivative prescribed at left end.
1247	ELSE IF ( ibcbeg == 2 ) then
1248
1249	c(4,1) = 2._wp
1250	c(3,1) = 1._wp
1251	c(2,1) = 3._wp * c(4,2) - c(3,2) / 2._wp * c(2,1)
1252
1253	ELSE
1254	WRITE(ctlmes,*) 'CUBSPL - Error, invalid IBCBEG input option!'
1255	CALL ctl_stop( ctlmes )
1256	END IF
1257
1258	! If there are interior knots, generate the corresponding
1259	! equations and carry out the forward pass of Gauss,
1260	! elimination after which the I-th equation reads:
1261	!
1262	! C(4,I) * S(I) + C(3,I) * S(I+1) = C(2,I).
1263
1264	IF ( n > 2 ) THEN
1265
1266	DO i = 2, n-1
1267	g = -c(3,i+1) / c(4,i-1)
1268	c(2,i) = g * c(2,i-1) + 3._wp * &
1269	( c(3,i) * c(4,i+1) + c(3,i+1) * c(4,i) )
1270	c(4,i) = g * c(3,i-1) + 2._wp * ( c(3,i) + c(3,i+1))
1271	END DO
1272
1273	! Construct the last equation from the second boundary,
1274	! condition of the form
1275	!
1276	! -G * C(4,N-1) * S(N-1) + C(4,N) * S(N) = C(2,N)
1277	!
1278	! If 1st der. is prescribed at right end ( ibcend == 1 ),
1279	! one can go directly to back-substitution, since the C
1280	! array happens to be set up just right for it at this
1281	! point.
1282
1283	IF ( ibcend < 1 ) THEN
1284	! Not-a-knot and 3 <= N, and either 3 < N or also
1285	! not-a-knot at left end point.
1286	IF ( n /= 3 .OR. ibcbeg /= 0 ) THEN
1287	g = c(3,n-1) + c(3,n)
1288	c(2,n) = ( ( c(3,n) + 2._wp * g ) * c(4,n) * c(3,n-1) &
1289	+ c(3,n)*2 ( c(1,n-1) - c(1,n-2) ) / &
1290	c(3,n-1) ) / g
1291	g = - g / c(4,n-1)
1292	c(4,n) = c(3,n-1)
1293	c(4,n) = c(4,n) + g * c(3,n-1)
1294	c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
1295	ELSE
1296	! N = 3 and not-a-knot also at left.
1297	c(2,n) = 2._wp * c(4,n)
1298	c(4,n) = 1._wp
1299	g = -1._wp / c(4,n-1)
1300	c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
1301	c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
1302	END IF
1303
1304	ELSE IF ( ibcend == 2 ) THEN
1305	! IBCEND = 2: Second derivative prescribed at right endpoint.
1306	c(2,n) = 3._wp * c(4,n) + c(3,n) / 2._wp * c(2,n)
1307	c(4,n) = 2._wp
1308	g = -1._wp / c(4,n-1)
1309	c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
1310	c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
1311	END IF
1312
1313	ELSE
1314	! N = 2 (assumed to be at least equal to 2!).
1315
1316	IF ( ibcend == 2 ) THEN
1317
1318	c(2,n) = 3._wp * c(4,n) + c(3,n) / 2._wp * c(2,n)
1319	c(4,n) = 2._wp
1320	g = -1._wp / c(4,n-1)
1321	c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
1322	c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
1323
1324	ELSE IF ( ibcend == 0 .AND. ibcbeg /= 0 ) THEN
1325
1326	c(2,n) = 2._wp * c(4,n)
1327	c(4,n) = 1._wp
1328	g = -1._wp / c(4,n-1)
1329	c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1)
1330	c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n)
1331
1332	ELSE IF ( ibcend == 0 .AND. ibcbeg == 0 ) THEN
1333
1334	c(2,n) = c(4,n)
1335
1336	END IF
1337
1338	END IF
1339
1340	! Back solve the upper triangular system
1341	!
1342	! C(4,I) * S(I) + C(3,I) * S(I+1) = B(I)
1343	!
1344	! for the slopes C(2,I), given that S(N) is already known.
1345	!
1346	DO i = n-1, 1, -1
1347	c(2,i) = ( c(2,i) - c(3,i) * c(2,i+1) ) / c(4,i)
1348	END DO
1349
1350	! Generate cubic coefficients in each interval, that is, the
1351	! derivatives at its left endpoint, from value and slope at its
1352	! endpoints.
1353	DO i = 2, n
1354	dtau = c(3,i)
1355	divdf1 = ( c(1,i) - c(1,i-1) ) / dtau
1356	divdf3 = c(2,i-1) + c(2,i) - 2._wp * divdf1
1357	c(3,i-1) = 2._wp * ( divdf1 - c(2,i-1) - divdf3 ) / dtau
1358	c(4,i-1) = 6._wp * divdf3 / dtau**2
1359	END DO
1360
1361	END FUNCTION cub_spl
1362
1363	! =====================================================================================================
1364
1365	FUNCTION z_cub_spl(s, tau, c) RESULT ( z_cs )
1366	!----------------------------------------------------------------------
1367	! * function z_cub_spl *
1368	!
1369	! ** Purpose : evaluate the cubic spline F in the interval
1370	! (TAU(I), TAU(I+1)). F is given by
1371	!
1372	!
1373	! F(X) = C1 + C2H + (C3H^2)/2 + (C4*H^3)/6
1374	!
1375	! where H = S-TAU(I).
1376	!
1377	! s MUST be positive: 0 <= s <= 1
1378	!---------------------------------------------------------------------
1379	INTEGER, PARAMETER :: n = 2
1380	REAL(wp), INTENT(in ) :: s
1381	REAL(wp), INTENT(in ) :: tau(n)
1382	REAL(wp), INTENT(in ) :: c(4,n)
1383	REAL(wp) :: z_cs
1384	REAL(wp) :: c1, c2, c3, c4
1385	REAL(wp) :: H
1386	!---------------------------------------------------------------------
1387	c1 = c(1,1)
1388	c2 = c(2,1)
1389	c3 = c(3,1)
1390	c4 = c(4,1)
1391
1392	H = s - tau(1)
1393
1394	z_cs = c1 + c2 * H + (c3 * H*2._wp)/2._wp + (c4 H**3._wp)/6._wp
1395
1396	END FUNCTION z_cub_spl
1397
1398	! =====================================================================================================
1399
1400	FUNCTION D1z_cub_spl(s, tau, c, kmax) RESULT ( d1z_cs )
1401	!----------------------------------------------------------------------
1402	! * function D1z_cub_spl *
1403	!
1404	! ** Purpose : evaluate the 1st derivative of cubic spline F
1405	! in the interval (TAU(I), TAU(I+1)).
1406	! The 1st derivative D1F is given by
1407	!
1408	!
1409	! D1F(S) = C1 + C2H + (C3H^2)/2 + (C4*H^3)/6
1410	!
1411	! where H = S-TAU(I).
1412	!
1413	! s MUST be positive: 0 <= s <= 1
1414	!---------------------------------------------------------------------
1415	INTEGER, PARAMETER :: n = 2
1416	REAL(wp), INTENT(in ) :: s
1417	REAL(wp), INTENT(in ) :: tau(n)
1418	REAL(wp), INTENT(in ) :: c(4,n)
1419	INTEGER, INTENT(in ) :: kmax
1420	REAL(wp) :: d1z_cs
1421	REAL(wp) :: c2, c3, c4
1422	REAL(wp) :: H
1423	!---------------------------------------------------------------------
1424	c2 = c(2,1) / REAL(kmax,wp)
1425	c3 = c(3,1) / REAL(kmax,wp)
1426	c4 = c(4,1) / REAL(kmax,wp)
1427
1428	H = s - tau(1)
1429
1430	d1z_cs = c2 + c3 * H + (c4 * H**2._wp)/2._wp
1431
1432	END FUNCTION D1z_cub_spl
1433
1434	! =====================================================================================================
1435
1436
1437	END MODULE mes

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/mes.F90 @ 15733

Download in other formats: