New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domvvl.F90 in branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago

#2050 first version. Compiled OK in moci test suite

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