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/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 3887

Last change on this file since 3887 was 3887, checked in by acc, 11 years ago

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