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

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

Last change on this file since 10759 was 10759, checked in by andmirek, 5 years ago

GMED 450 write output.namelist.dyn only for nprint > 2

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