source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 7698

Last change on this file since 7698 was 7698, checked in by mocavero, 4 years ago

update trunk with OpenMP parallelization

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