source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/tests/CANAL/MY_SRC/domvvl.F90 @ 10170

Last change on this file since 10170 was 10170, checked in by smasson, 2 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 2a: add report calls of lbc_lnk, see #2133

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