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

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

Branch 2013/dev_r3858_NOC_ZTC, #863. Changes to make new vvl code restartable. Vertical scale factors are read in dom_vvl_init which is called from domain before istate. Introduced rst_read_open in restart.F90 so that the restart file can be optionally opened in advance of the call to rst_read
  • Property svn:keywords set to Id
File size: 65.3 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            CALL rst_read_open                  !  open the restart file if necessary
705            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
706            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
707            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
708            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
709            id5 = iom_varid( numror, 'hdif_lf', ldstop = .FALSE. )
710            !                             ! --------- !
711            !                             ! all cases !
712            !                             ! --------- !
713            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
714               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
715               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
716               IF( neuler == 0 ) THEN
717                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
718               ENDIF
719            ELSE                                 ! one at least array is missing
720               CALL ctl_stop( 'dom_vvl_rst: vvl cannot restart from a non vvl run' )
721            ENDIF
722            !                             ! ----------- !
723            IF( ln_vvl_zstar ) THEN       ! z_star case !
724               !                          ! ----------- !
725               IF( MIN( id3, id4 ) > 0 ) THEN
726                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
727               ENDIF
728               !                          ! ----------------------- !
729            ELSE                          ! z_tilde and layer cases !
730               !                          ! ----------------------- !
731               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
732                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
733                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
734               ELSE                            ! one at least array is missing
735                  tilde_e3t_b(:,:,:) = 0.e0
736                  tilde_e3t_n(:,:,:) = 0.e0
737               ENDIF
738               !                          ! ------------ !
739               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
740                  !                       ! ------------ !
741                  IF( id5 > 0 ) THEN  ! required array exists
742                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
743                  ELSE                ! array is missing
744                     hdiv_lf(:,:,:) = 0.e0
745                  ENDIF
746               ENDIF
747            ENDIF
748            !
749         ELSE                                   !* Initialize at "rest"
750            fse3t_b(:,:,:) = e3t_0(:,:,:)
751            fse3t_n(:,:,:) = e3t_0(:,:,:)
752            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
753               tilde_e3t_b(:,:,:) = 0.e0
754               tilde_e3t_n(:,:,:) = 0.e0
755               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.e0
756            END IF
757         ENDIF
758
759      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
760         !                                   ! ===================
761         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
762         !                                           ! --------- !
763         !                                           ! all cases !
764         !                                           ! --------- !
765         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
766         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
767         !                                           ! ----------------------- !
768         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
769            !                                        ! ----------------------- !
770            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
771            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
772         END IF
773         !                                           ! -------------!   
774         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
775            !                                        ! ------------ !
776            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
777         ENDIF
778
779      ENDIF
780      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
781
782   END SUBROUTINE dom_vvl_rst
783
784
785   SUBROUTINE dom_vvl_ctl
786      !!---------------------------------------------------------------------
787      !!                  ***  ROUTINE dom_vvl_ctl  ***
788      !!               
789      !! ** Purpose :   Control the consistency between namelist options
790      !!                for vertical coordinate
791      !!----------------------------------------------------------------------
792      INTEGER ::   ioptio
793
794      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
795                      &      rn_ahe3,    rn_rst_e3t, rn_lf_cutoff, rn_zdef_max            , &
796                      &   ln_vvl_dbg               ! not yet implemented: ln_vvl_kepe
797      !!----------------------------------------------------------------------
798
799      REWIND ( numnam )               ! Read Namelist nam_vvl : vertical coordinate
800      READ   ( numnam, nam_vvl )
801
802      IF(lwp) THEN                    ! Namelist print
803         WRITE(numout,*)
804         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
805         WRITE(numout,*) '~~~~~~~~~~~'
806         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
807         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
808         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
809         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
810         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
811         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
812         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
813         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
814         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
815         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
816         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
817         IF( ln_vvl_ztilde_as_zstar ) THEN
818            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
819            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
820            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
821            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
822            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
823            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
824         ELSE
825            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
826            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
827            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
828            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
829         ENDIF
830         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
831         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
832      ENDIF
833
834      ioptio = 0                      ! Parameter control
835      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
836      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
837      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
838      IF( ln_vvl_layer           )        ioptio = ioptio + 1
839
840      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
841
842      IF(lwp) THEN                   ! Print the choice
843         WRITE(numout,*)
844         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
845         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
846         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
847         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
848         ! - ML - Option not developed yet
849         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
850         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
851      ENDIF
852
853   END SUBROUTINE dom_vvl_ctl
854
855   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
856      !!---------------------------------------------------------------------
857      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
858      !!                     
859      !! ** Purpose :   Correct surface weighted, horizontally interpolated,
860      !!                scale factors at locations that have been individually
861      !!                modified in domhgr. Such modifications break the
862      !!                relationship between e12t and e1u*e2u etc.
863      !!                Recompute some scale factors ignoring the modified metric.
864      !!----------------------------------------------------------------------
865      !! * Arguments
866      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
867      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
868      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
869      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
870      !! * Local declarations
871      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
872      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
873      !                                                ! =====================
874      IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN    ! ORCA R2 configuration
875         !                                             ! =====================
876         IF( nn_cla == 0 ) THEN
877            !
878            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
879            ij0 = 102   ;   ij1 = 102
880            DO jk = 1, jpkm1
881               DO jj = mj0(ij0), mj1(ij1)
882                  DO ji = mi0(ii0), mi1(ii1)
883                     SELECT CASE ( pout )
884                     CASE( 'U' )
885                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
886                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
887                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
888                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
889                     CASE( 'F' )
890                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
891                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
892                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
893                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
894                     END SELECT
895                  END DO
896               END DO
897            END DO
898            !
899            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
900            ij0 =  88   ;   ij1 =  88
901            DO jk = 1, jpkm1
902               DO jj = mj0(ij0), mj1(ij1)
903                  DO ji = mi0(ii0), mi1(ii1)
904                     SELECT CASE ( pout )
905                     CASE( 'U' )
906                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
907                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
908                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
909                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
910                     CASE( 'V' )
911                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
912                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
913                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
914                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
915                     CASE( 'F' )
916                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
917                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
918                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
919                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
920                     END SELECT
921                  END DO
922               END DO
923            END DO
924         ENDIF
925
926         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
927         ij0 = 116   ;   ij1 = 116
928         DO jk = 1, jpkm1
929            DO jj = mj0(ij0), mj1(ij1)
930               DO ji = mi0(ii0), mi1(ii1)
931                  SELECT CASE ( pout )
932                  CASE( 'U' )
933                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
934                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
935                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
936                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
937                  CASE( 'F' )
938                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
939                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
940                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
941                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
942                  END SELECT
943               END DO
944            END DO
945         END DO
946      ENDIF
947      !
948         !                                             ! =====================
949      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
950         !                                             ! =====================
951         !
952         ii0 = 281   ;   ii1 = 282        ! Gibraltar Strait (e2u was modified)
953         ij0 = 200   ;   ij1 = 200
954         DO jk = 1, jpkm1
955            DO jj = mj0(ij0), mj1(ij1)
956               DO ji = mi0(ii0), mi1(ii1)
957                  SELECT CASE ( pout )
958                  CASE( 'U' )
959                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
960                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
961                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
962                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
963                  CASE( 'F' )
964                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
965                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
966                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
967                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
968                  END SELECT
969               END DO
970            END DO
971         END DO
972         !
973         ii0 = 314   ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
974         ij0 = 208   ;   ij1 = 208
975         DO jk = 1, jpkm1
976            DO jj = mj0(ij0), mj1(ij1)
977               DO ji = mi0(ii0), mi1(ii1)
978                  SELECT CASE ( pout )
979                  CASE( 'U' )
980                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
981                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
982                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
983                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
984                  CASE( 'F' )
985                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
986                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
987                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
988                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
989                  END SELECT
990               END DO
991            END DO
992         END DO
993         !
994         ii0 =  44   ;   ii1 =  44        ! Lombok Strait (e1v was modified)
995         ij0 = 124   ;   ij1 = 125
996         DO jk = 1, jpkm1
997            DO jj = mj0(ij0), mj1(ij1)
998               DO ji = mi0(ii0), mi1(ii1)
999                  SELECT CASE ( pout )
1000                  CASE( 'V' )
1001                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1002                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1003                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1004                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1005                  END SELECT
1006               END DO
1007            END DO
1008         END DO
1009         !
1010         ii0 =  48   ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
1011         ij0 = 124   ;   ij1 = 125
1012         DO jk = 1, jpkm1
1013            DO jj = mj0(ij0), mj1(ij1)
1014               DO ji = mi0(ii0), mi1(ii1)
1015                  SELECT CASE ( pout )
1016                  CASE( 'V' )
1017                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1018                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1019                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1020                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1021                  END SELECT
1022               END DO
1023            END DO
1024         END DO
1025         !
1026         ii0 =  53   ;   ii1 =  53        ! Ombai Strait (e1v was modified)
1027         ij0 = 124   ;   ij1 = 125
1028         DO jk = 1, jpkm1
1029            DO jj = mj0(ij0), mj1(ij1)
1030               DO ji = mi0(ii0), mi1(ii1)
1031                  SELECT CASE ( pout )
1032                  CASE( 'V' )
1033                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1034                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1035                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1036                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1037                  END SELECT
1038               END DO
1039            END DO
1040         END DO
1041         !
1042         ii0 =  56   ;   ii1 =  56        ! Timor Passage (e1v was modified)
1043         ij0 = 124   ;   ij1 = 125
1044         DO jk = 1, jpkm1
1045            DO jj = mj0(ij0), mj1(ij1)
1046               DO ji = mi0(ii0), mi1(ii1)
1047                  SELECT CASE ( pout )
1048                  CASE( 'V' )
1049                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1050                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1051                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1052                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1053                  END SELECT
1054               END DO
1055            END DO
1056         END DO
1057         !
1058         ii0 =  55   ;   ii1 =  55        ! West Halmahera Strait (e1v was modified)
1059         ij0 = 141   ;   ij1 = 142
1060         DO jk = 1, jpkm1
1061            DO jj = mj0(ij0), mj1(ij1)
1062               DO ji = mi0(ii0), mi1(ii1)
1063                  SELECT CASE ( pout )
1064                  CASE( 'V' )
1065                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1066                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1067                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1068                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1069                  END SELECT
1070               END DO
1071            END DO
1072         END DO
1073         !
1074         ii0 =  58   ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
1075         ij0 = 141   ;   ij1 = 142
1076         DO jk = 1, jpkm1
1077            DO jj = mj0(ij0), mj1(ij1)
1078               DO ji = mi0(ii0), mi1(ii1)
1079                  SELECT CASE ( pout )
1080                  CASE( 'V' )
1081                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1082                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1083                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1084                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1085                  END SELECT
1086               END DO
1087            END DO
1088         END DO
1089      ENDIF
1090         !                                             ! =====================
1091      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
1092         !                                             ! =====================
1093         !
1094         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
1095         ij0 = 327   ;   ij1 = 327
1096         DO jk = 1, jpkm1
1097            DO jj = mj0(ij0), mj1(ij1)
1098               DO ji = mi0(ii0), mi1(ii1)
1099                  SELECT CASE ( pout )
1100                  CASE( 'U' )
1101                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1102                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1103                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1104                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1105                  CASE( 'F' )
1106                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1107                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1108                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1109                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1110                  END SELECT
1111               END DO
1112            END DO
1113         END DO
1114         !
1115         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
1116         ij0 = 343   ;   ij1 = 343
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         !
1136         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
1137         ij0 = 232   ;   ij1 = 232
1138         DO jk = 1, jpkm1
1139            DO jj = mj0(ij0), mj1(ij1)
1140               DO ji = mi0(ii0), mi1(ii1)
1141                  SELECT CASE ( pout )
1142                  CASE( 'U' )
1143                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1144                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1145                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1146                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1147                  CASE( 'F' )
1148                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1149                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1150                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1151                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1152                  END SELECT
1153               END DO
1154            END DO
1155         END DO
1156         !
1157         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
1158         ij0 = 232   ;   ij1 = 232
1159         DO jk = 1, jpkm1
1160            DO jj = mj0(ij0), mj1(ij1)
1161               DO ji = mi0(ii0), mi1(ii1)
1162                  SELECT CASE ( pout )
1163                  CASE( 'U' )
1164                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1165                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1166                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1167                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1168                  CASE( 'F' )
1169                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1170                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1171                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1172                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1173                  END SELECT
1174               END DO
1175            END DO
1176         END DO
1177         !
1178         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
1179         ij0 = 270   ;   ij1 = 270
1180         DO jk = 1, jpkm1
1181            DO jj = mj0(ij0), mj1(ij1)
1182               DO ji = mi0(ii0), mi1(ii1)
1183                  SELECT CASE ( pout )
1184                  CASE( 'U' )
1185                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1186                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1187                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1188                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1189                  CASE( 'F' )
1190                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1191                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1192                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1193                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1194                  END SELECT
1195               END DO
1196            END DO
1197         END DO
1198         !
1199         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
1200         ij0 = 232   ;   ij1 = 233
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 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
1216         ij0 = 276   ;   ij1 = 276
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      ENDIF
1231   END SUBROUTINE dom_vvl_orca_fix
1232
1233   !!======================================================================
1234END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.