source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 20 months ago

GMED 450 add flush after prints

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