source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM/domvvl.F90 @ 11823

Last change on this file since 11823 was 11823, checked in by mathiot, 12 months ago

rm useless USE statement, option compatibility test + minor changes

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