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

Last change on this file since 5014 was 5014, checked in by hliu, 6 years ago

upload the modifications for W/D based on r:4826

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