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.
domvvl.F90 in branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5441

Last change on this file since 5441 was 5441, checked in by hliu, 9 years ago

new vertical length scale updating/interpolation for wetting and drying
  • Property svn:keywords set to Id
File size: 75.6 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
8   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl:
9   !!                                          vvl option includes z_star and z_tilde coordinates
10   !!----------------------------------------------------------------------
11   !!   'key_vvl'                              variable volume
12   !!----------------------------------------------------------------------
13   !!----------------------------------------------------------------------
14   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
15   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
16   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
17   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
18   !!   dom_vvl_rst      : read/write restart file
19   !!   dom_vvl_ctl      : Check the vvl options
20   !!   dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors
21   !!                    : to account for manual changes to e[1,2][u,v] in some Straits
22   !!----------------------------------------------------------------------
23   !! * Modules used
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
26   USE sbc_oce         ! ocean surface boundary condition
27   USE in_out_manager  ! I/O manager
28   USE iom             ! I/O manager library
29   USE restart         ! ocean restart
30   USE lib_mpp         ! distributed memory computing library
31   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
32   USE wrk_nemo        ! Memory allocation
33   USE timing          ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   !! * Routine accessibility
39   PUBLIC  dom_vvl_init       ! called by domain.F90
40   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
41   PUBLIC  dom_vvl_sf_swp     ! called by step.F90
42   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
43   PRIVATE dom_vvl_orca_fix   ! called by dom_vvl_interpol
44
45   !!* Namelist nam_vvl
46   LOGICAL , PUBLIC                                      :: ln_vvl_zstar              ! zstar  vertical coordinate
47   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde             ! ztilde vertical coordinate
48   LOGICAL , PUBLIC                                      :: ln_vvl_layer              ! level  vertical coordinate
49   LOGICAL , PUBLIC                                      :: ln_vvl_ztilde_as_zstar    ! ztilde vertical coordinate
50   LOGICAL , PUBLIC                                      :: ln_vvl_zstar_at_eqtor     ! ztilde vertical coordinate
51   LOGICAL , PUBLIC                                      :: ln_vvl_kepe               ! kinetic/potential energy transfer
52   !                                                                                           ! conservation: not used yet
53   REAL(wp)                                              :: rn_ahe3                   ! thickness diffusion coefficient
54   REAL(wp)                                              :: rn_rst_e3t                ! ztilde to zstar restoration timescale [days]
55   REAL(wp)                                              :: rn_lf_cutoff              ! cutoff frequency for low-pass filter  [days]
56   REAL(wp)                                              :: rn_zdef_max               ! maximum fractional e3t deformation
57   LOGICAL , PUBLIC                                      :: ln_vvl_dbg                ! debug control prints
58
59   !! * Module variables
60   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                       ! thickness diffusion transport
61   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                            ! low frequency part of hz divergence
62   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n           ! baroclinic scale factors
63   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a          ! baroclinic scale factors
64   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                        ! retoring period for scale factors
65   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                        ! retoring period for low freq. divergence
66
67   !! * Substitutions
68#  include "domzgr_substitute.h90"
69#  include "vectopt_loop_substitute.h90"
70   !!----------------------------------------------------------------------
71   !! NEMO/OPA 3.3 , NEMO-Consortium (2010)
72   !! $Id$
73   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
74   !!----------------------------------------------------------------------
75
76CONTAINS
77
78   INTEGER FUNCTION dom_vvl_alloc()
79      !!----------------------------------------------------------------------
80      !!                ***  FUNCTION dom_vvl_alloc  ***
81      !!----------------------------------------------------------------------
82      IF( ln_vvl_zstar ) dom_vvl_alloc = 0
83      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
84         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
85            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
86            &      STAT = dom_vvl_alloc        )
87         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
88         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
89         un_td = 0.0_wp
90         vn_td = 0.0_wp
91      ENDIF
92      IF( ln_vvl_ztilde ) THEN
93         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
94         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
95         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
96      ENDIF
97
98   END FUNCTION dom_vvl_alloc
99
100
101   SUBROUTINE dom_vvl_init
102      !!----------------------------------------------------------------------
103      !!                ***  ROUTINE dom_vvl_init  ***
104      !!                   
105      !! ** Purpose :  Initialization of all scale factors, depths
106      !!               and water column heights
107      !!
108      !! ** Method  :  - use restart file and/or initialize
109      !!               - interpolate scale factors
110      !!
111      !! ** Action  : - fse3t_(n/b) and tilde_e3t_(n/b)
112      !!              - Regrid: fse3(u/v)_n
113      !!                        fse3(u/v)_b       
114      !!                        fse3w_n           
115      !!                        fse3(u/v)w_b     
116      !!                        fse3(u/v)w_n     
117      !!                        fsdept_n, fsdepw_n and fsde3w_n
118      !!              - h(t/u/v)_0
119      !!              - frq_rst_e3t and frq_rst_hdv
120      !!
121      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
122      !!----------------------------------------------------------------------
123      USE phycst,  ONLY : rpi, rsmall, rad
124      !! * Local declarations
125      INTEGER ::   ji,jj,jk
126      INTEGER ::   ii0, ii1, ij0, ij1
127      !!----------------------------------------------------------------------
128      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init')
129
130      IF(lwp) WRITE(numout,*)
131      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
132      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
133
134      ! choose vertical coordinate (z_star, z_tilde or layer)
135      ! ==========================
136      CALL dom_vvl_ctl
137
138      ! Allocate module arrays
139      ! ======================
140      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
141
142      ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk))
143      ! ============================================================================
144      CALL dom_vvl_rst( nit000, 'READ' )
145      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
146
147      !IF(ln_wd) THEN
148      !  DO jj = 1, jpj
149      !    DO ji = 1, jpi
150      !      IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN
151      !       fse3t_a(ji,jj,1:2) = 0.5_wp * rn_wdmin1
152      !       fse3t_n(ji,jj,1:2) = 0.5_wp * rn_wdmin1
153      !       fse3t_b(ji,jj,1:2) = 0.5_wp * rn_wdmin1
154      !      END IF
155      !    ENDDO
156      !  ENDDO
157      !END IF
158
159      ! Reconstruction of all vertical scale factors at now and before time steps
160      ! =============================================================================
161      ! Horizontal scale factor interpolations
162      ! --------------------------------------
163      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
164      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
165      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
166      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
167      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
168      ! Vertical scale factor interpolations
169      ! ------------------------------------
170      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
171      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
172      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
173      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
174      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
175      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
176      ! t- and w- points depth
177      ! ----------------------
178      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
179      fsdepw_n(:,:,1) = 0.0_wp
180      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
181      fsdept_b(:,:,1) = 0.5_wp * fse3w_b(:,:,1)
182      fsdepw_b(:,:,1) = 0.0_wp
183      DO jk = 2, jpk
184         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
185         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
186         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
187         fsdept_b(:,:,jk) = fsdept_b(:,:,jk-1) + fse3w_b(:,:,jk)
188         fsdepw_b(:,:,jk) = fsdepw_b(:,:,jk-1) + fse3t_b(:,:,jk-1)
189      END DO
190
191      ! Before depth and Inverse of the local depth of the water column at u- and v- points
192      ! -----------------------------------------------------------------------------------
193      hu_b(:,:) = 0.
194      hv_b(:,:) = 0.
195      DO jk = 1, jpkm1
196         hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk)
197         hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk)
198      END DO
199      hur_b(:,:) = umask(:,:,1) / ( hu_b(:,:) + 1. - umask(:,:,1) )
200      hvr_b(:,:) = vmask(:,:,1) / ( hv_b(:,:) + 1. - vmask(:,:,1) )
201
202      ! Restoring frequencies for z_tilde coordinate
203      ! ============================================
204      IF( ln_vvl_ztilde ) THEN
205         ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings
206         frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
207         frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
208         IF( ln_vvl_ztilde_as_zstar ) THEN
209            ! Ignore namelist settings and use these next two to emulate z-star using z-tilde
210            frq_rst_e3t(:,:) = 0.0_wp 
211            frq_rst_hdv(:,:) = 1.0_wp / rdt
212         ENDIF
213         IF ( ln_vvl_zstar_at_eqtor ) THEN
214            DO jj = 1, jpj
215               DO ji = 1, jpi
216                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
217                     ! values outside the equatorial band and transition zone (ztilde)
218                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
219                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
220                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN
221                     ! values inside the equatorial band (ztilde as zstar)
222                     frq_rst_e3t(ji,jj) =  0.0_wp
223                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
224                  ELSE
225                     ! values in the transition band (linearly vary from ztilde to ztilde as zstar values)
226                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
227                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
228                        &                                          * 180._wp / 3.5_wp ) )
229                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
230                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
231                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
232                        &                                          * 180._wp / 3.5_wp ) )
233                  ENDIF
234               END DO
235            END DO
236            IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN
237               ii0 = 103   ;   ii1 = 111        ! Suppress ztilde in the Foxe Basin for ORCA2
238               ij0 = 128   ;   ij1 = 135   ;   
239               frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
240               frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rdt
241            ENDIF
242         ENDIF
243      ENDIF
244
245      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_init')
246
247   END SUBROUTINE dom_vvl_init
248
249
250   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
251      !!----------------------------------------------------------------------
252      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
253      !!                   
254      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
255      !!                 tranxt and dynspg routines
256      !!
257      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
258      !!               - z_tilde_case: after scale factor increment =
259      !!                                    high frequency part of horizontal divergence
260      !!                                  + retsoring towards the background grid
261      !!                                  + thickness difusion
262      !!                               Then repartition of ssh INCREMENT proportionnaly
263      !!                               to the "baroclinic" level thickness.
264      !!
265      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
266      !!               - tilde_e3t_a: after increment of vertical scale factor
267      !!                              in z_tilde case
268      !!               - fse3(t/u/v)_a
269      !!
270      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
271      !!----------------------------------------------------------------------
272      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t
273      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv
274      !! * Arguments
275      INTEGER, INTENT( in )                  :: kt                    ! time step
276      INTEGER, INTENT( in ), OPTIONAL        :: kcall                 ! optional argument indicating call sequence
277      !! * Local declarations
278      INTEGER                                :: ji, jj, jk            ! dummy loop indices
279      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
280      REAL(wp)                               :: z2dt                  ! temporary scalars
281      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
282      LOGICAL                                :: ll_do_bclinic         ! temporary logical
283      !!----------------------------------------------------------------------
284      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
285      CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
286      CALL wrk_alloc( jpi, jpj, jpk, ze3t                     )
287
288      IF(kt == nit000)   THEN
289         IF(lwp) WRITE(numout,*)
290         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
291         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
292      ENDIF
293
294      ll_do_bclinic = .TRUE.
295      IF( PRESENT(kcall) ) THEN
296         IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.
297      ENDIF
298
299      ! ******************************* !
300      ! After acale factors at t-points !
301      ! ******************************* !
302
303      !                                                ! --------------------------------------------- !
304                                                       ! z_star coordinate and barotropic z-tilde part !
305      !                                                ! --------------------------------------------- !
306
307      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
308      DO jk = 1, jpkm1
309         ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0)
310         fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
311      END DO
312
313      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
314         !                                                            ! ------baroclinic part------ !
315
316         ! I - initialization
317         ! ==================
318
319         ! 1 - barotropic divergence
320         ! -------------------------
321         zhdiv(:,:) = 0.
322         zht(:,:)   = 0.
323         DO jk = 1, jpkm1
324            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
325            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
326         END DO
327         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) )
328
329         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
330         ! --------------------------------------------------
331         IF( ln_vvl_ztilde ) THEN
332            IF( kt .GT. nit000 ) THEN
333               DO jk = 1, jpkm1
334                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
335                     &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
336               END DO
337            ENDIF
338         END IF
339
340         ! II - after z_tilde increments of vertical scale factors
341         ! =======================================================
342         tilde_e3t_a(:,:,:) = 0.0_wp  ! tilde_e3t_a used to store tendency terms
343
344         ! 1 - High frequency divergence term
345         ! ----------------------------------
346         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
347            DO jk = 1, jpkm1
348               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
349            END DO
350         ELSE                         ! layer case
351            DO jk = 1, jpkm1
352               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) )
353            END DO
354         END IF
355
356         ! 2 - Restoring term (z-tilde case only)
357         ! ------------------
358         IF( ln_vvl_ztilde ) THEN
359            DO jk = 1, jpk
360               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
361            END DO
362         END IF
363
364         ! 3 - Thickness diffusion term
365         ! ----------------------------
366         zwu(:,:) = 0.0_wp
367         zwv(:,:) = 0.0_wp
368         ! a - first derivative: diffusive fluxes
369         DO jk = 1, jpkm1
370            DO jj = 1, jpjm1
371               DO ji = 1, fs_jpim1   ! vector opt.
372                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) &
373                                  & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
374                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 
375                                  & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
376                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
377                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
378               END DO
379            END DO
380         END DO
381         ! b - correction for last oceanic u-v points
382         DO jj = 1, jpj
383            DO ji = 1, jpi
384               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
385               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
386            END DO
387         END DO
388         ! c - second derivative: divergence of diffusive fluxes
389         DO jk = 1, jpkm1
390            DO jj = 2, jpjm1
391               DO ji = fs_2, fs_jpim1   ! vector opt.
392                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
393                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
394                     &                                            ) * r1_e12t(ji,jj)
395               END DO
396            END DO
397         END DO
398         ! d - thickness diffusion transport: boundary conditions
399         !     (stored for tracer advction and continuity equation)
400         CALL lbc_lnk( un_td , 'U' , -1.)
401         CALL lbc_lnk( vn_td , 'V' , -1.)
402
403         ! 4 - Time stepping of baroclinic scale factors
404         ! ---------------------------------------------
405         ! Leapfrog time stepping
406         ! ~~~~~~~~~~~~~~~~~~~~~~
407         IF( neuler == 0 .AND. kt == nit000 ) THEN
408            z2dt =  rdt
409         ELSE
410            z2dt = 2.0_wp * rdt
411         ENDIF
412         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. )
413         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
414
415         ! Maximum deformation control
416         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
417         ze3t(:,:,jpk) = 0.0_wp
418         DO jk = 1, jpkm1
419            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
420         END DO
421         z_tmax = MAXVAL( ze3t(:,:,:) )
422         IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
423         z_tmin = MINVAL( ze3t(:,:,:) )
424         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
425         ! - ML - test: for the moment, stop simulation for too large e3_t variations
426         IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN
427            IF( lk_mpp ) THEN
428               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
429               CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
430            ELSE
431               ijk_max = MAXLOC( ze3t(:,:,:) )
432               ijk_max(1) = ijk_max(1) + nimpp - 1
433               ijk_max(2) = ijk_max(2) + njmpp - 1
434               ijk_min = MINLOC( ze3t(:,:,:) )
435               ijk_min(1) = ijk_min(1) + nimpp - 1
436               ijk_min(2) = ijk_min(2) + njmpp - 1
437            ENDIF
438            IF (lwp) THEN
439               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
440               WRITE(numout, *) 'at i, j, k=', ijk_max
441               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
442               WRITE(numout, *) 'at i, j, k=', ijk_min           
443               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
444            ENDIF
445         ENDIF
446         ! - ML - end test
447         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
448         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
449         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
450
451         !
452         ! "tilda" change in the after scale factor
453         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
454         DO jk = 1, jpkm1
455            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
456         END DO
457         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
458         ! ==================================================================================
459         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
460         ! - ML - baroclinicity error should be better treated in the future
461         !        i.e. locally and not spread over the water column.
462         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
463         zht(:,:) = 0.
464         DO jk = 1, jpkm1
465            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
466         END DO
467         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
468         DO jk = 1, jpkm1
469            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
470         END DO
471
472      ENDIF
473
474      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
475      !                                           ! ---baroclinic part--------- !
476         DO jk = 1, jpkm1
477            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk)
478         END DO
479      ENDIF
480
481      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
482         !
483         IF( lwp ) WRITE(numout, *) 'kt =', kt
484         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
485            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
486            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
487            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
488         END IF
489         !
490         zht(:,:) = 0.0_wp
491         DO jk = 1, jpkm1
492            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
493         END DO
494         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
495         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
496         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax
497         !
498         zht(:,:) = 0.0_wp
499         DO jk = 1, jpkm1
500            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
501         END DO
502         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
503         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
504         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax
505         !
506         zht(:,:) = 0.0_wp
507         DO jk = 1, jpkm1
508            zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk)
509         END DO
510         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
511         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
512         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax
513         !
514         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
515         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
516         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
517         !
518         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
519         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
520         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
521         !
522         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
523         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
524         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
525      END IF
526
527      ! *********************************** !
528      ! After scale factors at u- v- points !
529      ! *********************************** !
530
531      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
532      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
533
534      ! *********************************** !
535      ! After depths at u- v points         !
536      ! *********************************** !
537
538      hu_a(:,:) = 0._wp                        ! Ocean depth at U-points
539      hv_a(:,:) = 0._wp                        ! Ocean depth at V-points
540      DO jk = 1, jpkm1
541         hu_a(:,:) = hu_a(:,:) + fse3u_a(:,:,jk) * umask(:,:,jk)
542         hv_a(:,:) = hv_a(:,:) + fse3v_a(:,:,jk) * vmask(:,:,jk)
543      END DO
544      !                                        ! Inverse of the local depth
545      hur_a(:,:) = 1._wp / ( hu_a(:,:) + 1._wp - umask(:,:,1) ) * umask(:,:,1)
546      hvr_a(:,:) = 1._wp / ( hv_a(:,:) + 1._wp - vmask(:,:,1) ) * vmask(:,:,1)
547
548      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
549      CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
550
551      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_nxt')
552
553   END SUBROUTINE dom_vvl_sf_nxt
554
555
556   SUBROUTINE dom_vvl_sf_swp( kt )
557      !!----------------------------------------------------------------------
558      !!                ***  ROUTINE dom_vvl_sf_swp  ***
559      !!                   
560      !! ** Purpose :  compute time filter and swap of scale factors
561      !!               compute all depths and related variables for next time step
562      !!               write outputs and restart file
563      !!
564      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
565      !!               - reconstruct scale factor at other grid points (interpolate)
566      !!               - recompute depths and water height fields
567      !!
568      !! ** Action  :  - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step
569      !!               - Recompute:
570      !!                    fse3(u/v)_b       
571      !!                    fse3w_n           
572      !!                    fse3(u/v)w_b     
573      !!                    fse3(u/v)w_n     
574      !!                    fsdept_n, fsdepw_n  and fsde3w_n
575      !!                    h(u/v) and h(u/v)r
576      !!
577      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
578      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
579      !!----------------------------------------------------------------------
580      !! * Arguments
581      INTEGER, INTENT( in )               :: kt       ! time step
582      !! * Local declarations
583      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def
584      INTEGER                             :: jk       ! dummy loop indices
585      !!----------------------------------------------------------------------
586
587      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
588      !
589      CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def                )
590      !
591      IF( kt == nit000 )   THEN
592         IF(lwp) WRITE(numout,*)
593         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
594         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
595      ENDIF
596      !
597      ! Time filter and swap of scale factors
598      ! =====================================
599      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
600      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
601         IF( neuler == 0 .AND. kt == nit000 ) THEN
602            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
603         ELSE
604            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
605            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
606         ENDIF
607         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
608      ENDIF
609      fsdept_b(:,:,:) = fsdept_n(:,:,:)
610      fsdepw_b(:,:,:) = fsdepw_n(:,:,:)
611
612      fse3t_n(:,:,:) = fse3t_a(:,:,:)
613      fse3u_n(:,:,:) = fse3u_a(:,:,:)
614      fse3v_n(:,:,:) = fse3v_a(:,:,:)
615
616      ! Compute all missing vertical scale factor and depths
617      ! ====================================================
618      ! Horizontal scale factor interpolations
619      ! --------------------------------------
620      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
621      ! - JC - hu_b, hv_b, hur_b, hvr_b also
622      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
623      ! Vertical scale factor interpolations
624      ! ------------------------------------
625      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
626      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
627      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
628      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3w_b (:,:,:), 'W'  )
629      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
630      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
631      ! t- and w- points depth
632      ! ----------------------
633      fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1)
634      fsdepw_n(:,:,1) = 0.0_wp
635      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
636      DO jk = 2, jpk
637         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
638         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
639         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
640      END DO
641      ! Local depth and Inverse of the local depth of the water column at u- and v- points
642      ! ----------------------------------------------------------------------------------
643      hu (:,:) = hu_a (:,:)
644      hv (:,:) = hv_a (:,:)
645
646      ! Inverse of the local depth
647      hur(:,:) = hur_a(:,:)
648      hvr(:,:) = hvr_a(:,:)
649
650      ! Local depth of the water column at t- points
651      ! --------------------------------------------
652      ht(:,:) = 0.
653      DO jk = 1, jpkm1
654         ht(:,:) = ht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
655      END DO
656
657      ! Write outputs
658      ! =============
659      z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
660      CALL iom_put( "cellthc" , fse3t_n  (:,:,:) )
661      CALL iom_put( "tpt_dep" , fsde3w_n (:,:,:) )
662      CALL iom_put( "e3tdef"  , z_e3t_def(:,:,:) )
663
664      ! write restart file
665      ! ==================
666      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
667      !
668      CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def )
669      !
670      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
671
672   END SUBROUTINE dom_vvl_sf_swp
673
674
675   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
676      !!---------------------------------------------------------------------
677      !!                  ***  ROUTINE dom_vvl__interpol  ***
678      !!
679      !! ** Purpose :   interpolate scale factors from one grid point to another
680      !!
681      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
682      !!                - horizontal interpolation: grid cell surface averaging
683      !!                - vertical interpolation: simple averaging
684      !!----------------------------------------------------------------------
685      !! * Arguments
686      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
687      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
688      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
689      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
690      !! * Local declarations
691      REAL(wp) ::  zwad                                                ! = 1.0 when ln_wd = .true.
692                                                                       ! = 0.0 when ln_wd = .false.
693                                                                       !
694      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
695      LOGICAL ::   l_is_orca                                           ! local logical
696      !!----------------------------------------------------------------------
697      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
698         !
699      IF(ln_wd) THEN
700        zwad = 1.0_wp
701      ELSE
702        zwad = 0.0_wp
703      END IF
704
705      l_is_orca = .FALSE.
706      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE.      ! ORCA R2 configuration - will need to correct some locations
707
708      SELECT CASE ( pout )
709         !               ! ------------------------------------- !
710      CASE( 'U' )        ! interpolation from T-point to U-point !
711         !               ! ------------------------------------- !
712         ! horizontal surface weighted interpolation
713         DO jk = 1, jpk
714            DO jj = 1, jpjm1
715               DO ji = 1, fs_jpim1   ! vector opt.
716                  !pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj)                                   &
717                  pe3_out(ji,jj,jk) = 0.5_wp * (umask(ji,jj,jk) * (1.0_wp - zwad) + zwad) * r1_e12u(ji,jj)        &
718                     &                       * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
719                     &                           + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
720               END DO
721            END DO
722         END DO
723         !
724         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
725         ! boundary conditions
726         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )
727         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
728         !               ! ------------------------------------- !
729      CASE( 'V' )        ! interpolation from T-point to V-point !
730         !               ! ------------------------------------- !
731         ! horizontal surface weighted interpolation
732         DO jk = 1, jpk
733            DO jj = 1, jpjm1
734               DO ji = 1, fs_jpim1   ! vector opt.
735                  !pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                   &
736                  pe3_out(ji,jj,jk) = 0.5_wp * (vmask(ji,jj,jk) * (1.0_wp - zwad) + zwad) * r1_e12v(ji,jj)        &
737                     &                       * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
738                     &                           + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
739               END DO
740            END DO
741         END DO
742         !
743         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
744         ! boundary conditions
745         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )
746         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
747         !               ! ------------------------------------- !
748      CASE( 'F' )        ! interpolation from U-point to F-point !
749         !               ! ------------------------------------- !
750         ! horizontal surface weighted interpolation
751         DO jk = 1, jpk
752            DO jj = 1, jpjm1
753               DO ji = 1, fs_jpim1   ! vector opt.
754                  !pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)               &
755                  pe3_out(ji,jj,jk) = 0.5_wp * (umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zwad) + zwad)     &
756                     &                       * r1_e12f(ji,jj)                                                     &
757                     &                       * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
758                     &                           + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
759               END DO
760            END DO
761         END DO
762         !
763         IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout )
764         ! boundary conditions
765         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )
766         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
767         !               ! ------------------------------------- !
768      CASE( 'W' )        ! interpolation from T-point to W-point !
769         !               ! ------------------------------------- !
770         ! vertical simple interpolation
771         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
772         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
773         DO jk = 2, jpk
774            !pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
775            !   &                            +            0.5_wp * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
776            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (tmask(:,:,jk) * (1.0_wp - zwad) + zwad) ) &
777               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                         &
778               &                            +            0.5_wp * (tmask(:,:,jk) * (1.0_wp - zwad) + zwad)   &
779               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
780         END DO
781         !               ! -------------------------------------- !
782      CASE( 'UW' )       ! interpolation from U-point to UW-point !
783         !               ! -------------------------------------- !
784         ! vertical simple interpolation
785         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
786         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
787         DO jk = 2, jpk
788            !pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
789            !   &                             +            0.5_wp * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
790            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (umask(:,:,jk) * (1.0_wp - zwad) + zwad) ) &
791               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                         &
792               &                             +            0.5_wp * (umask(:,:,jk) * (1.0_wp - zwad) + zwad)   &
793               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
794         END DO
795         !               ! -------------------------------------- !
796      CASE( 'VW' )       ! interpolation from V-point to VW-point !
797         !               ! -------------------------------------- !
798         ! vertical simple interpolation
799         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
800         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
801         DO jk = 2, jpk
802            !pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
803            !   &                             +            0.5_wp * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
804            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * (vmask(:,:,jk) * (1.0_wp - zwad) + zwad) ) &
805               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                         &
806               &                             +            0.5_wp * (vmask(:,:,jk) * (1.0_wp - zwad) + zwad)   &
807               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
808         END DO
809      END SELECT
810      !
811
812      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
813
814   END SUBROUTINE dom_vvl_interpol
815
816   SUBROUTINE dom_vvl_rst( kt, cdrw )
817      !!---------------------------------------------------------------------
818      !!                   ***  ROUTINE dom_vvl_rst  ***
819      !!                     
820      !! ** Purpose :   Read or write VVL file in restart file
821      !!
822      !! ** Method  :   use of IOM library
823      !!                if the restart does not contain vertical scale factors,
824      !!                they are set to the _0 values
825      !!                if the restart does not contain vertical scale factors increments (z_tilde),
826      !!                they are set to 0.
827      !!----------------------------------------------------------------------
828      !! * Arguments
829      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
830      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
831      !! * Local declarations
832      INTEGER ::   ji, jj, jk
833      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
834      !!----------------------------------------------------------------------
835      !
836      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
837      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
838         !                                   ! ===============
839         IF( ln_rstart ) THEN                   !* Read the restart file
840            CALL rst_read_open                  !  open the restart file if necessary
841            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn    )
842            !
843            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
844            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
845            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
846            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
847            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
848            !                             ! --------- !
849            !                             ! all cases !
850            !                             ! --------- !
851            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
852               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
853               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
854               IF( neuler == 0 ) THEN
855                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
856               ENDIF
857            ELSE IF( id1 > 0 ) THEN
858               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files'
859               IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.'
860               IF(lwp) write(numout,*) 'neuler is forced to 0'
861               fse3t_b(:,:,:) = fse3t_n(:,:,:)
862               neuler = 0
863            ELSE
864               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart file'
865               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
866               IF(lwp) write(numout,*) 'neuler is forced to 0'
867               DO jk=1,jpk
868                  fse3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
869                      &                            / ( ht_0(:,:) + 1._wp - tmask(:,:,1) ) * tmask(:,:,jk) &
870                      &            + e3t_0(:,:,jk) * (1._wp -tmask(:,:,jk))
871               END DO
872               fse3t_b(:,:,:) = fse3t_n(:,:,:)
873               neuler = 0
874            ENDIF
875            !                             ! ----------- !
876            IF( ln_vvl_zstar ) THEN       ! z_star case !
877               !                          ! ----------- !
878               IF( MIN( id3, id4 ) > 0 ) THEN
879                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
880               ENDIF
881               !                          ! ----------------------- !
882            ELSE                          ! z_tilde and layer cases !
883               !                          ! ----------------------- !
884               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
885                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
886                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
887               ELSE                            ! one at least array is missing
888                  tilde_e3t_b(:,:,:) = 0.0_wp
889                  tilde_e3t_n(:,:,:) = 0.0_wp
890               ENDIF
891               !                          ! ------------ !
892               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
893                  !                       ! ------------ !
894                  IF( id5 > 0 ) THEN  ! required array exists
895                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
896                  ELSE                ! array is missing
897                     hdiv_lf(:,:,:) = 0.0_wp
898                  ENDIF
899               ENDIF
900            ENDIF
901            !
902         ELSE                                   !* Initialize at "rest"
903            fse3t_b(:,:,:) = e3t_0(:,:,:)
904            fse3t_n(:,:,:) = e3t_0(:,:,:)
905            sshn(:,:) = 0.0_wp
906
907            IF(ln_wd) THEN
908              DO jj = 1, jpj
909                DO ji = 1, jpi
910                  !IF(e3t_0(ji,jj,1) < 0._wp) THEN
911                  !IF(mbathy(ji,jj) == 2 .AND. e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN
912                  IF( e3t_0(ji,jj,1) <= 0.5_wp * rn_wdmin1) THEN
913                    fse3t_b(ji,jj,:) = 0.5_wp * rn_wdmin1 
914                    fse3t_n(ji,jj,:) = 0.5_wp * rn_wdmin1 
915                    fse3t_a(ji,jj,:) = 0.5_wp * rn_wdmin1 
916                    sshb(ji,jj) = rn_wdmin1 - bathy(ji,jj)
917                    sshn(ji,jj) = rn_wdmin1 - bathy(ji,jj)
918                    ssha(ji,jj) = rn_wdmin1 - bathy(ji,jj)
919                  ENDIF
920                ENDDO
921              ENDDO
922            END IF
923
924            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
925               tilde_e3t_b(:,:,:) = 0.0_wp
926               tilde_e3t_n(:,:,:) = 0.0_wp
927               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp
928            END IF
929         ENDIF
930
931      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
932         !                                   ! ===================
933         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
934         !                                           ! --------- !
935         !                                           ! all cases !
936         !                                           ! --------- !
937         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
938         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
939         !                                           ! ----------------------- !
940         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
941            !                                        ! ----------------------- !
942            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
943            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
944         END IF
945         !                                           ! -------------!   
946         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
947            !                                        ! ------------ !
948            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
949         ENDIF
950
951      ENDIF
952      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
953
954   END SUBROUTINE dom_vvl_rst
955
956
957   SUBROUTINE dom_vvl_ctl
958      !!---------------------------------------------------------------------
959      !!                  ***  ROUTINE dom_vvl_ctl  ***
960      !!               
961      !! ** Purpose :   Control the consistency between namelist options
962      !!                for vertical coordinate
963      !!----------------------------------------------------------------------
964      INTEGER ::   ioptio
965      INTEGER ::   ios
966
967      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
968                      & ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
969                      & rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
970      !!----------------------------------------------------------------------
971
972      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
973      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
974901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
975
976      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
977      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
978902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
979      IF(lwm) WRITE ( numond, nam_vvl )
980
981      IF(lwp) THEN                    ! Namelist print
982         WRITE(numout,*)
983         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
984         WRITE(numout,*) '~~~~~~~~~~~'
985         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
986         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
987         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
988         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
989         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
990         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
991         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
992         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
993         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
994         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
995         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
996         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
997         IF( ln_vvl_ztilde_as_zstar ) THEN
998            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
999            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
1000            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
1001            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
1002            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1003            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
1004         ELSE
1005            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
1006            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
1007            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
1008            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
1009         ENDIF
1010         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
1011         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
1012      ENDIF
1013
1014      ioptio = 0                      ! Parameter control
1015      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
1016      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
1017      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
1018      IF( ln_vvl_layer           )        ioptio = ioptio + 1
1019
1020      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1021
1022      IF(lwp) THEN                   ! Print the choice
1023         WRITE(numout,*)
1024         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
1025         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
1026         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
1027         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
1028         ! - ML - Option not developed yet
1029         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
1030         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
1031      ENDIF
1032
1033#if defined key_agrif
1034      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
1035#endif
1036
1037   END SUBROUTINE dom_vvl_ctl
1038
1039   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1040      !!---------------------------------------------------------------------
1041      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
1042      !!                     
1043      !! ** Purpose :   Correct surface weighted, horizontally interpolated,
1044      !!                scale factors at locations that have been individually
1045      !!                modified in domhgr. Such modifications break the
1046      !!                relationship between e12t and e1u*e2u etc.
1047      !!                Recompute some scale factors ignoring the modified metric.
1048      !!----------------------------------------------------------------------
1049      !! * Arguments
1050      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
1051      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
1052      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
1053      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
1054      !! * Local declarations
1055      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
1056      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
1057      !! acc
1058      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
1059      !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
1060      !!
1061      !                                                ! =====================
1062      IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration
1063         !                                             ! =====================
1064      !! acc
1065         IF( nn_cla == 0 ) THEN
1066            !
1067            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
1068            ij0 = 102   ;   ij1 = 102
1069            DO jk = 1, jpkm1
1070               DO jj = mj0(ij0), mj1(ij1)
1071                  DO ji = mi0(ii0), mi1(ii1)
1072                     SELECT CASE ( pout )
1073                     CASE( 'U' )
1074                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1075                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1076                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1077                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1078                     CASE( 'F' )
1079                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1080                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1081                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1082                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1083                     END SELECT
1084                  END DO
1085               END DO
1086            END DO
1087            !
1088            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
1089            ij0 =  88   ;   ij1 =  88
1090            DO jk = 1, jpkm1
1091               DO jj = mj0(ij0), mj1(ij1)
1092                  DO ji = mi0(ii0), mi1(ii1)
1093                     SELECT CASE ( pout )
1094                     CASE( 'U' )
1095                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1096                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1097                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1098                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1099                     CASE( 'V' )
1100                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1101                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1102                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1103                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1104                     CASE( 'F' )
1105                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1106                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1107                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1108                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1109                     END SELECT
1110                  END DO
1111               END DO
1112            END DO
1113         ENDIF
1114
1115         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
1116         ij0 = 116   ;   ij1 = 116
1117         DO jk = 1, jpkm1
1118            DO jj = mj0(ij0), mj1(ij1)
1119               DO ji = mi0(ii0), mi1(ii1)
1120                  SELECT CASE ( pout )
1121                  CASE( 'U' )
1122                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1123                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1124                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1125                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1126                  CASE( 'F' )
1127                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1128                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1129                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1130                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1131                  END SELECT
1132               END DO
1133            END DO
1134         END DO
1135      ENDIF
1136      !
1137         !                                             ! =====================
1138      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
1139         !                                             ! =====================
1140         !
1141         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
1142         ij0 = 200   ;   ij1 = 200
1143         DO jk = 1, jpkm1
1144            DO jj = mj0(ij0), mj1(ij1)
1145               DO ji = mi0(ii0), mi1(ii1)
1146                  SELECT CASE ( pout )
1147                  CASE( 'U' )
1148                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1149                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1150                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1151                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1152                  CASE( 'F' )
1153                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1154                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1155                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1156                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1157                  END SELECT
1158               END DO
1159            END DO
1160         END DO
1161         !
1162         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
1163         ij0 = 208   ;   ij1 = 208
1164         DO jk = 1, jpkm1
1165            DO jj = mj0(ij0), mj1(ij1)
1166               DO ji = mi0(ii0), mi1(ii1)
1167                  SELECT CASE ( pout )
1168                  CASE( 'U' )
1169                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1170                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1171                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1172                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1173                  CASE( 'F' )
1174                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1175                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1176                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1177                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1178                  END SELECT
1179               END DO
1180            END DO
1181         END DO
1182         !
1183         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
1184         ij0 = 124   ;   ij1 = 125
1185         DO jk = 1, jpkm1
1186            DO jj = mj0(ij0), mj1(ij1)
1187               DO ji = mi0(ii0), mi1(ii1)
1188                  SELECT CASE ( pout )
1189                  CASE( 'V' )
1190                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1191                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1192                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1193                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1194                  END SELECT
1195               END DO
1196            END DO
1197         END DO
1198         !
1199         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1200         ij0 = 124   ;   ij1 = 125
1201         DO jk = 1, jpkm1
1202            DO jj = mj0(ij0), mj1(ij1)
1203               DO ji = mi0(ii0), mi1(ii1)
1204                  SELECT CASE ( pout )
1205                  CASE( 'V' )
1206                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1207                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1208                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1209                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1210                  END SELECT
1211               END DO
1212            END DO
1213         END DO
1214         !
1215         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
1216         ij0 = 124   ;   ij1 = 125
1217         DO jk = 1, jpkm1
1218            DO jj = mj0(ij0), mj1(ij1)
1219               DO ji = mi0(ii0), mi1(ii1)
1220                  SELECT CASE ( pout )
1221                  CASE( 'V' )
1222                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1223                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1224                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1225                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1226                  END SELECT
1227               END DO
1228            END DO
1229         END DO
1230         !
1231         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
1232         ij0 = 124   ;   ij1 = 125
1233         DO jk = 1, jpkm1
1234            DO jj = mj0(ij0), mj1(ij1)
1235               DO ji = mi0(ii0), mi1(ii1)
1236                  SELECT CASE ( pout )
1237                  CASE( 'V' )
1238                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1239                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1240                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1241                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1242                  END SELECT
1243               END DO
1244            END DO
1245         END DO
1246         !
1247         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
1248         ij0 = 141   ;   ij1 = 142
1249         DO jk = 1, jpkm1
1250            DO jj = mj0(ij0), mj1(ij1)
1251               DO ji = mi0(ii0), mi1(ii1)
1252                  SELECT CASE ( pout )
1253                  CASE( 'V' )
1254                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1255                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1256                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1257                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1258                  END SELECT
1259               END DO
1260            END DO
1261         END DO
1262         !
1263         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
1264         ij0 = 141   ;   ij1 = 142
1265         DO jk = 1, jpkm1
1266            DO jj = mj0(ij0), mj1(ij1)
1267               DO ji = mi0(ii0), mi1(ii1)
1268                  SELECT CASE ( pout )
1269                  CASE( 'V' )
1270                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1271                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1272                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1273                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1274                  END SELECT
1275               END DO
1276            END DO
1277         END DO
1278      ENDIF
1279         !                                             ! =====================
1280      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
1281         !                                             ! =====================
1282         !
1283         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
1284         ij0 = 327   ;   ij1 = 327
1285         DO jk = 1, jpkm1
1286            DO jj = mj0(ij0), mj1(ij1)
1287               DO ji = mi0(ii0), mi1(ii1)
1288                  SELECT CASE ( pout )
1289                  CASE( 'U' )
1290                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1291                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1292                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1293                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1294                  CASE( 'F' )
1295                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1296                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1297                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1298                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1299                  END SELECT
1300               END DO
1301            END DO
1302         END DO
1303         !
1304         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
1305         ij0 = 343   ;   ij1 = 343
1306         DO jk = 1, jpkm1
1307            DO jj = mj0(ij0), mj1(ij1)
1308               DO ji = mi0(ii0), mi1(ii1)
1309                  SELECT CASE ( pout )
1310                  CASE( 'U' )
1311                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1312                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1313                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1314                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1315                  CASE( 'F' )
1316                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1317                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1318                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1319                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1320                  END SELECT
1321               END DO
1322            END DO
1323         END DO
1324         !
1325         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
1326         ij0 = 232   ;   ij1 = 232
1327         DO jk = 1, jpkm1
1328            DO jj = mj0(ij0), mj1(ij1)
1329               DO ji = mi0(ii0), mi1(ii1)
1330                  SELECT CASE ( pout )
1331                  CASE( 'U' )
1332                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1333                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1334                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1335                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1336                  CASE( 'F' )
1337                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1338                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1339                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1340                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1341                  END SELECT
1342               END DO
1343            END DO
1344         END DO
1345         !
1346         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
1347         ij0 = 232   ;   ij1 = 232
1348         DO jk = 1, jpkm1
1349            DO jj = mj0(ij0), mj1(ij1)
1350               DO ji = mi0(ii0), mi1(ii1)
1351                  SELECT CASE ( pout )
1352                  CASE( 'U' )
1353                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1354                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1355                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1356                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1357                  CASE( 'F' )
1358                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1359                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1360                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1361                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1362                  END SELECT
1363               END DO
1364            END DO
1365         END DO
1366         !
1367         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
1368         ij0 = 270   ;   ij1 = 270
1369         DO jk = 1, jpkm1
1370            DO jj = mj0(ij0), mj1(ij1)
1371               DO ji = mi0(ii0), mi1(ii1)
1372                  SELECT CASE ( pout )
1373                  CASE( 'U' )
1374                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1375                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1376                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1377                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1378                  CASE( 'F' )
1379                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1380                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1381                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1382                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1383                  END SELECT
1384               END DO
1385            END DO
1386         END DO
1387         !
1388         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
1389         ij0 = 232   ;   ij1 = 233
1390         DO jk = 1, jpkm1
1391            DO jj = mj0(ij0), mj1(ij1)
1392               DO ji = mi0(ii0), mi1(ii1)
1393                  SELECT CASE ( pout )
1394                  CASE( 'V' )
1395                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1396                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1397                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1398                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1399                  END SELECT
1400               END DO
1401            END DO
1402         END DO
1403         !
1404         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
1405         ij0 = 276   ;   ij1 = 276
1406         DO jk = 1, jpkm1
1407            DO jj = mj0(ij0), mj1(ij1)
1408               DO ji = mi0(ii0), mi1(ii1)
1409                  SELECT CASE ( pout )
1410                  CASE( 'V' )
1411                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1412                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1413                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1414                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1415                  END SELECT
1416               END DO
1417            END DO
1418         END DO
1419      ENDIF
1420   END SUBROUTINE dom_vvl_orca_fix
1421
1422   !!======================================================================
1423END MODULE domvvl
1424
1425
1426
Note: See TracBrowser for help on using the repository browser.