source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 6348

Last change on this file since 6348 was 6348, checked in by cetlod, 4 years ago

bugfix: move the output of scale factor before time swapping and output some variables needs for offline, see ticket #1682

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