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.
mes.F90 in NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src – NEMO

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

Last change on this file since 15271 was 15271, checked in by dbruciaferri, 17 months ago

rewriting saw tooth dignostic

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