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

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

#2050 fixes and changes

File size: 76.6 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
7   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
8   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl:
9   !!                                          vvl option includes z_star and z_tilde coordinates
10   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
11   !!----------------------------------------------------------------------
12   !!   'key_vvl'                              variable volume
13   !!----------------------------------------------------------------------
14   !!----------------------------------------------------------------------
15   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
16   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
17   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
18   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
19   !!   dom_vvl_rst      : read/write restart file
20   !!   dom_vvl_ctl      : Check the vvl options
21   !!   dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors
22   !!                    : to account for manual changes to e[1,2][u,v] in some Straits
23   !!----------------------------------------------------------------------
24   !! * Modules used
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE sbc_oce         ! ocean surface boundary condition
28   USE in_out_manager  ! I/O manager
29   USE iom             ! I/O manager library
30   USE restart         ! ocean restart
31   USE lib_mpp         ! distributed memory computing library
32   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
33   USE wrk_nemo        ! Memory allocation
34   USE timing          ! Timing
35   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      CONTINUE
955      ENDIF
956      call mpp_bcast(ios)
957      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
958      IF(lwm) THEN
959         REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
960         READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
961902      CONTINUE
962      ENDIF
963      call mpp_bcast(ios)
964      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
965
966      IF(lwm) WRITE ( numond, nam_vvl )
967
968      CALL dom_namelist()
969
970      IF(lwp) THEN                    ! Namelist print
971         WRITE(numout,*)
972         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
973         WRITE(numout,*) '~~~~~~~~~~~'
974         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
975         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
976         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
977         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
978         WRITE(numout,*) '              ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
979         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
980         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
981         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
982         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
983         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
984         WRITE(numout,*) '           Namelist nam_vvl : maximum e3t deformation fractional change'
985         WRITE(numout,*) '                                         rn_zdef_max    = ', rn_zdef_max
986         IF( ln_vvl_ztilde_as_zstar ) THEN
987            WRITE(numout,*) '           ztilde running in zstar emulation mode; '
988            WRITE(numout,*) '           ignoring namelist timescale parameters and using:'
989            WRITE(numout,*) '                 hard-wired : z-tilde to zstar restoration timescale (days)'
990            WRITE(numout,*) '                                         rn_rst_e3t     =    0.0'
991            WRITE(numout,*) '                 hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
992            WRITE(numout,*) '                                         rn_lf_cutoff   =    1.0/rdt'
993         ELSE
994            WRITE(numout,*) '           Namelist nam_vvl : z-tilde to zstar restoration timescale (days)'
995            WRITE(numout,*) '                                         rn_rst_e3t     = ', rn_rst_e3t
996            WRITE(numout,*) '           Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)'
997            WRITE(numout,*) '                                         rn_lf_cutoff   = ', rn_lf_cutoff
998         ENDIF
999         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
1000         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
1001      ENDIF
1002
1003      ioptio = 0                      ! Parameter control
1004      IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true.
1005      IF( ln_vvl_zstar           )        ioptio = ioptio + 1
1006      IF( ln_vvl_ztilde          )        ioptio = ioptio + 1
1007      IF( ln_vvl_layer           )        ioptio = ioptio + 1
1008
1009      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1010      IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1011
1012      IF(lwp) THEN                   ! Print the choice
1013         WRITE(numout,*)
1014         IF( ln_vvl_zstar           ) WRITE(numout,*) '              zstar vertical coordinate is used'
1015         IF( ln_vvl_ztilde          ) WRITE(numout,*) '              ztilde vertical coordinate is used'
1016         IF( ln_vvl_layer           ) WRITE(numout,*) '              layer vertical coordinate is used'
1017         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '              to emulate a zstar coordinate'
1018         ! - ML - Option not developed yet
1019         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
1020         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
1021      ENDIF
1022
1023#if defined key_agrif
1024      IF (.NOT.Agrif_Root()) CALL ctl_stop( 'AGRIF not implemented with non-linear free surface (key_vvl)' )
1025#endif
1026
1027   END SUBROUTINE dom_vvl_ctl
1028
1029   SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout )
1030      !!---------------------------------------------------------------------
1031      !!                   ***  ROUTINE dom_vvl_orca_fix  ***
1032      !!                     
1033      !! ** Purpose :   Correct surface weighted, horizontally interpolated,
1034      !!                scale factors at locations that have been individually
1035      !!                modified in domhgr. Such modifications break the
1036      !!                relationship between e12t and e1u*e2u etc.
1037      !!                Recompute some scale factors ignoring the modified metric.
1038      !!----------------------------------------------------------------------
1039      !! * Arguments
1040      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
1041      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
1042      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
1043      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
1044      !! * Local declarations
1045      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
1046      INTEGER ::   ij0, ij1, ii0, ii1                                  ! dummy loop indices
1047      INTEGER ::   isrow                                               ! index for ORCA1 starting row
1048      !! acc
1049      !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for
1050      !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations
1051      !!
1052      !                                                ! =====================
1053      IF( cp_cfg == "orca" .AND. jp_cfg == 3 ) THEN    ! ORCA R2 configuration
1054         !                                             ! =====================
1055      !! acc
1056         IF( nn_cla == 0 ) THEN
1057            !
1058            ii0 = 139   ;   ii1 = 140        ! Gibraltar Strait (e2u was modified)
1059            ij0 = 102   ;   ij1 = 102
1060            DO jk = 1, jpkm1
1061               DO jj = mj0(ij0), mj1(ij1)
1062                  DO ji = mi0(ii0), mi1(ii1)
1063                     SELECT CASE ( pout )
1064                     CASE( 'U' )
1065                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1066                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1067                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1068                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1069                     CASE( 'F' )
1070                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1071                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1072                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1073                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1074                     END SELECT
1075                  END DO
1076               END DO
1077            END DO
1078            !
1079            ii0 = 160   ;   ii1 = 160        ! Bab el Mandeb (e2u and e1v were modified)
1080            ij0 =  88   ;   ij1 =  88
1081            DO jk = 1, jpkm1
1082               DO jj = mj0(ij0), mj1(ij1)
1083                  DO ji = mi0(ii0), mi1(ii1)
1084                     SELECT CASE ( pout )
1085                     CASE( 'U' )
1086                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1087                       &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1088                       &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1089                       &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1090                     CASE( 'V' )
1091                        pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1092                       &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1093                       &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1094                       &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1095                     CASE( 'F' )
1096                        pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1097                       &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1098                       &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1099                       &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1100                     END SELECT
1101                  END DO
1102               END DO
1103            END DO
1104         ENDIF
1105
1106         ii0 = 145   ;   ii1 = 146        ! Danish Straits (e2u was modified)
1107         ij0 = 116   ;   ij1 = 116
1108         DO jk = 1, jpkm1
1109            DO jj = mj0(ij0), mj1(ij1)
1110               DO ji = mi0(ii0), mi1(ii1)
1111                  SELECT CASE ( pout )
1112                  CASE( 'U' )
1113                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1114                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1115                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1116                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1117                  CASE( 'F' )
1118                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1119                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1120                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1121                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1122                  END SELECT
1123               END DO
1124            END DO
1125         END DO
1126      ENDIF
1127      !
1128         !                                             ! =====================
1129      IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN    ! ORCA R1 configuration
1130         !                                             ! =====================
1131         ! This dirty section will be suppressed by simplification process:
1132         ! all this will come back in input files
1133         ! Currently these hard-wired indices relate to configuration with
1134         ! extend grid (jpjglo=332)
1135         ! which had a grid-size of 362x292.
1136         isrow = 332 - jpjglo
1137         !
1138         ii0 = 282           ;   ii1 = 283        ! Gibraltar Strait (e2u was modified)
1139         ij0 = 241 - isrow   ;   ij1 = 241 - isrow
1140         DO jk = 1, jpkm1
1141            DO jj = mj0(ij0), mj1(ij1)
1142               DO ji = mi0(ii0), mi1(ii1)
1143                  SELECT CASE ( pout )
1144                  CASE( 'U' )
1145                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1146                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1147                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1148                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1149                  CASE( 'F' )
1150                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1151                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1152                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1153                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1154                  END SELECT
1155               END DO
1156            END DO
1157         END DO
1158         !
1159         ii0 = 314           ;   ii1 = 315        ! Bhosporus Strait (e2u was modified)
1160         ij0 = 248 - isrow   ;   ij1 = 248 - isrow
1161         DO jk = 1, jpkm1
1162            DO jj = mj0(ij0), mj1(ij1)
1163               DO ji = mi0(ii0), mi1(ii1)
1164                  SELECT CASE ( pout )
1165                  CASE( 'U' )
1166                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1167                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1168                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1169                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1170                  CASE( 'F' )
1171                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1172                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1173                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1174                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1175                  END SELECT
1176               END DO
1177            END DO
1178         END DO
1179         !
1180         ii0 =  44           ;   ii1 =  44        ! Lombok Strait (e1v was modified)
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 =  48           ;   ii1 =  48        ! Sumba Strait (e1v was modified) [closed from bathy_11 on]
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 =  53          ;   ii1 =  53        ! Ombai Strait (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 =  56            ;   ii1 =  56        ! Timor Passage (e1v was modified)
1229         ij0 = 164 - isrow    ;   ij1 = 165  - 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 =  55            ;   ii1 =  55        ! West 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         !
1260         ii0 =  58            ;   ii1 =  58        ! East Halmahera Strait (e1v was modified)
1261         ij0 = 181 - isrow    ;   ij1 = 182 - isrow 
1262         DO jk = 1, jpkm1
1263            DO jj = mj0(ij0), mj1(ij1)
1264               DO ji = mi0(ii0), mi1(ii1)
1265                  SELECT CASE ( pout )
1266                  CASE( 'V' )
1267                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1268                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1269                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1270                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1271                  END SELECT
1272               END DO
1273            END DO
1274         END DO
1275      ENDIF
1276         !                                             ! =====================
1277      IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN   ! ORCA R05 configuration
1278         !                                             ! =====================
1279         !
1280         ii0 = 563   ;   ii1 = 564        ! Gibraltar Strait (e2u was modified)
1281         ij0 = 327   ;   ij1 = 327
1282         DO jk = 1, jpkm1
1283            DO jj = mj0(ij0), mj1(ij1)
1284               DO ji = mi0(ii0), mi1(ii1)
1285                  SELECT CASE ( pout )
1286                  CASE( 'U' )
1287                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1288                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1289                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1290                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1291                  CASE( 'F' )
1292                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1293                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1294                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1295                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1296                  END SELECT
1297               END DO
1298            END DO
1299         END DO
1300         !
1301         ii0 = 627   ;   ii1 = 628        ! Bosphorus Strait (e2u was modified)
1302         ij0 = 343   ;   ij1 = 343
1303         DO jk = 1, jpkm1
1304            DO jj = mj0(ij0), mj1(ij1)
1305               DO ji = mi0(ii0), mi1(ii1)
1306                  SELECT CASE ( pout )
1307                  CASE( 'U' )
1308                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        & 
1309                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1310                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1311                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1312                  CASE( 'F' )
1313                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    & 
1314                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1315                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1316                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1317                  END SELECT
1318               END DO
1319            END DO
1320         END DO
1321         !
1322         ii0 =  93   ;   ii1 =  94        ! Sumba Strait (e2u was modified)
1323         ij0 = 232   ;   ij1 = 232
1324         DO jk = 1, jpkm1
1325            DO jj = mj0(ij0), mj1(ij1)
1326               DO ji = mi0(ii0), mi1(ii1)
1327                  SELECT CASE ( pout )
1328                  CASE( 'U' )
1329                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1330                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1331                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1332                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1333                  CASE( 'F' )
1334                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1335                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1336                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1337                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1338                  END SELECT
1339               END DO
1340            END DO
1341         END DO
1342         !
1343         ii0 = 103   ;   ii1 = 103        ! Ombai Strait (e2u was modified)
1344         ij0 = 232   ;   ij1 = 232
1345         DO jk = 1, jpkm1
1346            DO jj = mj0(ij0), mj1(ij1)
1347               DO ji = mi0(ii0), mi1(ii1)
1348                  SELECT CASE ( pout )
1349                  CASE( 'U' )
1350                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1351                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1352                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1353                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1354                  CASE( 'F' )
1355                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1356                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1357                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1358                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1359                  END SELECT
1360               END DO
1361            END DO
1362         END DO
1363         !
1364         ii0 =  15   ;   ii1 =  15        ! Palk Strait (e2u was modified)
1365         ij0 = 270   ;   ij1 = 270
1366         DO jk = 1, jpkm1
1367            DO jj = mj0(ij0), mj1(ij1)
1368               DO ji = mi0(ii0), mi1(ii1)
1369                  SELECT CASE ( pout )
1370                  CASE( 'U' )
1371                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk)                                        &
1372                    &                    * (   e1t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) ) &
1373                    &                    +     e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) &
1374                    &                      ) / e1u(ji,jj)   +   e3u_0(ji,jj,jk)
1375                  CASE( 'F' )
1376                     pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk)                    &
1377                    &                    * (   e1u(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3u_0(ji  ,jj,jk) ) &
1378                    &                    +     e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) &
1379                    &                      ) / e1f(ji,jj)   +   e3f_0(ji,jj,jk)
1380                  END SELECT
1381               END DO
1382            END DO
1383         END DO
1384         !
1385         ii0 =  87   ;   ii1 =  87        ! Lombok Strait (e1v was modified)
1386         ij0 = 232   ;   ij1 = 233
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         !
1401         ii0 = 662   ;   ii1 = 662        ! Bab el Mandeb (e1v was modified)
1402         ij0 = 276   ;   ij1 = 276
1403         DO jk = 1, jpkm1
1404            DO jj = mj0(ij0), mj1(ij1)
1405               DO ji = mi0(ii0), mi1(ii1)
1406                  SELECT CASE ( pout )
1407                  CASE( 'V' )
1408                     pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk)                                        &
1409                    &                    * (   e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) ) &
1410                    &                    +     e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) &
1411                    &                      ) / e2v(ji,jj)   +   e3v_0(ji,jj,jk)
1412                  END SELECT
1413               END DO
1414            END DO
1415         END DO
1416      ENDIF
1417   END SUBROUTINE dom_vvl_orca_fix
1418
1419   SUBROUTINE dom_namelist()
1420     !!---------------------------------------------------------------------
1421     !!                   ***  ROUTINE dom_namelist  ***
1422     !!                     
1423     !! ** Purpose :   Broadcast namelist variables read by procesor lwm
1424     !!
1425     !! ** Method  :   use lib_mpp
1426     !!----------------------------------------------------------------------
1427#if defined key_mpp_mpi
1428      CALL mpp_bcast(ln_vvl_zstar)
1429      CALL mpp_bcast(ln_vvl_ztilde)
1430      CALL mpp_bcast(ln_vvl_layer)
1431      CALL mpp_bcast(ln_vvl_ztilde_as_zstar)
1432      CALL mpp_bcast(ln_vvl_zstar_at_eqtor)
1433      CALL mpp_bcast(rn_ahe3)
1434      CALL mpp_bcast(rn_rst_e3t)
1435      CALL mpp_bcast(rn_lf_cutoff)
1436      CALL mpp_bcast(rn_zdef_max)
1437      CALL mpp_bcast(ln_vvl_dbg)
1438#endif
1439   END SUBROUTINE dom_namelist
1440   !!======================================================================
1441END MODULE domvvl
1442
1443
1444
Note: See TracBrowser for help on using the repository browser.