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/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 4386

Last change on this file since 4386 was 4386, checked in by jchanut, 10 years ago

Typo error for timing

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