New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domvvl.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 9367

Last change on this file since 9367 was 9367, checked in by mathiot, 6 years ago

Add restart read/write via XIOS capability (#1953 and #1962 and twiki: 2017WP/Met_Office-1_Mirek_XIOSread). WARNING: need to upgrade XIOS to r1296 to compile
  • Property svn:keywords set to Id
File size: 54.9 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/OPA 3.7 , NEMO-Consortium (2015)
67   !! $Id$
68   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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      IF(lwxios) THEN
241! define variables in restart file when writing with XIOS
242         CALL iom_set_rstw_var_active('e3t_b')
243         CALL iom_set_rstw_var_active('e3t_n')
244         !                                           ! ----------------------- !
245         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
246            !                                        ! ----------------------- !
247            CALL iom_set_rstw_var_active('tilde_e3t_b')
248            CALL iom_set_rstw_var_active('tilde_e3t_n')
249         END IF
250         !                                           ! -------------!   
251         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
252            !                                        ! ------------ !
253            CALL iom_set_rstw_var_active('hdiv_lf')
254         ENDIF
255         !
256      ENDIF
257
258      IF( ln_timing )  CALL timing_stop('dom_vvl_init')
259      !
260   END SUBROUTINE dom_vvl_init
261
262
263   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
264      !!----------------------------------------------------------------------
265      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
266      !!                   
267      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
268      !!                 tranxt and dynspg routines
269      !!
270      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
271      !!               - z_tilde_case: after scale factor increment =
272      !!                                    high frequency part of horizontal divergence
273      !!                                  + retsoring towards the background grid
274      !!                                  + thickness difusion
275      !!                               Then repartition of ssh INCREMENT proportionnaly
276      !!                               to the "baroclinic" level thickness.
277      !!
278      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
279      !!               - tilde_e3t_a: after increment of vertical scale factor
280      !!                              in z_tilde case
281      !!               - e3(t/u/v)_a
282      !!
283      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
284      !!----------------------------------------------------------------------
285      INTEGER, INTENT( in )           ::   kt      ! time step
286      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
287      !
288      INTEGER                ::   ji, jj, jk            ! dummy loop indices
289      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
290      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
291      LOGICAL                ::   ll_do_bclinic         ! local logical
292      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
293      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
294      !!----------------------------------------------------------------------
295      !
296      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
297      !
298      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
299      !
300      IF( kt == nit000 ) THEN
301         IF(lwp) WRITE(numout,*)
302         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
303         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
304      ENDIF
305
306      ll_do_bclinic = .TRUE.
307      IF( PRESENT(kcall) ) THEN
308         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
309      ENDIF
310
311      ! ******************************* !
312      ! After acale factors at t-points !
313      ! ******************************* !
314      !                                                ! --------------------------------------------- !
315      !                                                ! z_star coordinate and barotropic z-tilde part !
316      !                                                ! --------------------------------------------- !
317      !
318      z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
319      DO jk = 1, jpkm1
320         ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0)
321         e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
322      END DO
323      !
324      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
325         !                                                            ! ------baroclinic part------ !
326         ! I - initialization
327         ! ==================
328
329         ! 1 - barotropic divergence
330         ! -------------------------
331         zhdiv(:,:) = 0._wp
332         zht(:,:)   = 0._wp
333         DO jk = 1, jpkm1
334            zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
335            zht  (:,:) = zht  (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
336         END DO
337         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
338
339         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
340         ! --------------------------------------------------
341         IF( ln_vvl_ztilde ) THEN
342            IF( kt > nit000 ) THEN
343               DO jk = 1, jpkm1
344                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
345                     &          * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
346               END DO
347            ENDIF
348         ENDIF
349
350         ! II - after z_tilde increments of vertical scale factors
351         ! =======================================================
352         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
353
354         ! 1 - High frequency divergence term
355         ! ----------------------------------
356         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
357            DO jk = 1, jpkm1
358               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
359            END DO
360         ELSE                         ! layer case
361            DO jk = 1, jpkm1
362               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
363            END DO
364         ENDIF
365
366         ! 2 - Restoring term (z-tilde case only)
367         ! ------------------
368         IF( ln_vvl_ztilde ) THEN
369            DO jk = 1, jpk
370               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
371            END DO
372         ENDIF
373
374         ! 3 - Thickness diffusion term
375         ! ----------------------------
376         zwu(:,:) = 0._wp
377         zwv(:,:) = 0._wp
378         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
379            DO jj = 1, jpjm1
380               DO ji = 1, fs_jpim1   ! vector opt.
381                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
382                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
383                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
384                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
385                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
386                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
387               END DO
388            END DO
389         END DO
390         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
391            DO ji = 1, jpi
392               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
393               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
394            END DO
395         END DO
396         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
397            DO jj = 2, jpjm1
398               DO ji = fs_2, fs_jpim1   ! vector opt.
399                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
400                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
401                     &                                            ) * r1_e1e2t(ji,jj)
402               END DO
403            END DO
404         END DO
405         !                       ! d - thickness diffusion transport: boundary conditions
406         !                             (stored for tracer advction and continuity equation)
407         CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
408
409         ! 4 - Time stepping of baroclinic scale factors
410         ! ---------------------------------------------
411         ! Leapfrog time stepping
412         ! ~~~~~~~~~~~~~~~~~~~~~~
413         IF( neuler == 0 .AND. kt == nit000 ) THEN
414            z2dt =  rdt
415         ELSE
416            z2dt = 2.0_wp * rdt
417         ENDIF
418         CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp )
419         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
420
421         ! Maximum deformation control
422         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
423         ze3t(:,:,jpk) = 0._wp
424         DO jk = 1, jpkm1
425            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
426         END DO
427         z_tmax = MAXVAL( ze3t(:,:,:) )
428         IF( lk_mpp )   CALL mpp_max( z_tmax )                 ! max over the global domain
429         z_tmin = MINVAL( ze3t(:,:,:) )
430         IF( lk_mpp )   CALL mpp_min( z_tmin )                 ! min over the global domain
431         ! - ML - test: for the moment, stop simulation for too large e3_t variations
432         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
433            IF( lk_mpp ) THEN
434               CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) )
435               CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) )
436            ELSE
437               ijk_max = MAXLOC( ze3t(:,:,:) )
438               ijk_max(1) = ijk_max(1) + nimpp - 1
439               ijk_max(2) = ijk_max(2) + njmpp - 1
440               ijk_min = MINLOC( ze3t(:,:,:) )
441               ijk_min(1) = ijk_min(1) + nimpp - 1
442               ijk_min(2) = ijk_min(2) + njmpp - 1
443            ENDIF
444            IF (lwp) THEN
445               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
446               WRITE(numout, *) 'at i, j, k=', ijk_max
447               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
448               WRITE(numout, *) 'at i, j, k=', ijk_min           
449               CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
450            ENDIF
451         ENDIF
452         ! - ML - end test
453         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
454         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
455         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
456
457         !
458         ! "tilda" change in the after scale factor
459         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
460         DO jk = 1, jpkm1
461            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
462         END DO
463         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
464         ! ==================================================================================
465         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
466         ! - ML - baroclinicity error should be better treated in the future
467         !        i.e. locally and not spread over the water column.
468         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
469         zht(:,:) = 0.
470         DO jk = 1, jpkm1
471            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
472         END DO
473         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )
474         DO jk = 1, jpkm1
475            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
476         END DO
477
478      ENDIF
479
480      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
481      !                                           ! ---baroclinic part--------- !
482         DO jk = 1, jpkm1
483            e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
484         END DO
485      ENDIF
486
487      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
488         !
489         IF( lwp ) WRITE(numout, *) 'kt =', kt
490         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
491            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
492            IF( lk_mpp ) CALL mpp_max( z_tmax )                             ! max over the global domain
493            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
494         END IF
495         !
496         zht(:,:) = 0.0_wp
497         DO jk = 1, jpkm1
498            zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
499         END DO
500         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
501         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
502         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax
503         !
504         zht(:,:) = 0.0_wp
505         DO jk = 1, jpkm1
506            zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk)
507         END DO
508         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
509         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
510         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax
511         !
512         zht(:,:) = 0.0_wp
513         DO jk = 1, jpkm1
514            zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk)
515         END DO
516         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) )
517         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
518         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax
519         !
520         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshb(:,:) ) )
521         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
522         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax
523         !
524         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( sshn(:,:) ) )
525         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
526         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax
527         !
528         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssha(:,:) ) )
529         IF( lk_mpp ) CALL mpp_max( z_tmax )                                ! max over the global domain
530         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax
531      END IF
532
533      ! *********************************** !
534      ! After scale factors at u- v- points !
535      ! *********************************** !
536
537      CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' )
538      CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' )
539
540      ! *********************************** !
541      ! After depths at u- v points         !
542      ! *********************************** !
543
544      hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1)
545      hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1)
546      DO jk = 2, jpkm1
547         hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk)
548         hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk)
549      END DO
550      !                                        ! Inverse of the local depth
551!!gm BUG ?  don't understand the use of umask_i here .....
552      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
553      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
554      !
555      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
556      !
557   END SUBROUTINE dom_vvl_sf_nxt
558
559
560   SUBROUTINE dom_vvl_sf_swp( kt )
561      !!----------------------------------------------------------------------
562      !!                ***  ROUTINE dom_vvl_sf_swp  ***
563      !!                   
564      !! ** Purpose :  compute time filter and swap of scale factors
565      !!               compute all depths and related variables for next time step
566      !!               write outputs and restart file
567      !!
568      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
569      !!               - reconstruct scale factor at other grid points (interpolate)
570      !!               - recompute depths and water height fields
571      !!
572      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_n ready for next time step
573      !!               - Recompute:
574      !!                    e3(u/v)_b       
575      !!                    e3w_n           
576      !!                    e3(u/v)w_b     
577      !!                    e3(u/v)w_n     
578      !!                    gdept_n, gdepw_n  and gde3w_n
579      !!                    h(u/v) and h(u/v)r
580      !!
581      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
582      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
583      !!----------------------------------------------------------------------
584      INTEGER, INTENT( in ) ::   kt   ! time step
585      !
586      INTEGER  ::   ji, jj, jk   ! dummy loop indices
587      REAL(wp) ::   zcoef        ! local scalar
588      !!----------------------------------------------------------------------
589      !
590      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
591      !
592      IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp')
593      !
594      IF( kt == nit000 )   THEN
595         IF(lwp) WRITE(numout,*)
596         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
597         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
598      ENDIF
599      !
600      ! Time filter and swap of scale factors
601      ! =====================================
602      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
603      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
604         IF( neuler == 0 .AND. kt == nit000 ) THEN
605            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
606         ELSE
607            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
608            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
609         ENDIF
610         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
611      ENDIF
612      gdept_b(:,:,:) = gdept_n(:,:,:)
613      gdepw_b(:,:,:) = gdepw_n(:,:,:)
614
615      e3t_n(:,:,:) = e3t_a(:,:,:)
616      e3u_n(:,:,:) = e3u_a(:,:,:)
617      e3v_n(:,:,:) = e3v_a(:,:,:)
618
619      ! Compute all missing vertical scale factor and depths
620      ! ====================================================
621      ! Horizontal scale factor interpolations
622      ! --------------------------------------
623      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
624      ! - JC - hu_b, hv_b, hur_b, hvr_b also
625     
626      CALL dom_vvl_interpol( e3u_n(:,:,:), e3f_n(:,:,:), 'F'  )
627     
628      ! Vertical scale factor interpolations
629      CALL dom_vvl_interpol( e3t_n(:,:,:),  e3w_n(:,:,:), 'W'  )
630      CALL dom_vvl_interpol( e3u_n(:,:,:), e3uw_n(:,:,:), 'UW' )
631      CALL dom_vvl_interpol( e3v_n(:,:,:), e3vw_n(:,:,:), 'VW' )
632      CALL dom_vvl_interpol( e3t_b(:,:,:),  e3w_b(:,:,:), 'W'  )
633      CALL dom_vvl_interpol( e3u_b(:,:,:), e3uw_b(:,:,:), 'UW' )
634      CALL dom_vvl_interpol( e3v_b(:,:,:), e3vw_b(:,:,:), 'VW' )
635
636      ! t- and w- points depth (set the isf depth as it is in the initial step)
637      gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1)
638      gdepw_n(:,:,1) = 0.0_wp
639      gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:)
640      DO jk = 2, jpk
641         DO jj = 1,jpj
642            DO ji = 1,jpi
643              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
644                                                                 ! 1 for jk = mikt
645               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
646               gdepw_n(ji,jj,jk) = gdepw_n(ji,jj,jk-1) + e3t_n(ji,jj,jk-1)
647               gdept_n(ji,jj,jk) =    zcoef  * ( gdepw_n(ji,jj,jk  ) + 0.5 * e3w_n(ji,jj,jk) )  &
648                   &             + (1-zcoef) * ( gdept_n(ji,jj,jk-1) +       e3w_n(ji,jj,jk) ) 
649               gde3w_n(ji,jj,jk) = gdept_n(ji,jj,jk) - sshn(ji,jj)
650            END DO
651         END DO
652      END DO
653
654      ! Local depth and Inverse of the local depth of the water
655      ! -------------------------------------------------------
656      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
657      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
658      !
659      ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1)
660      DO jk = 2, jpkm1
661         ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk)
662      END DO
663
664      ! write restart file
665      ! ==================
666      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
667      !
668      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
669      !
670   END SUBROUTINE dom_vvl_sf_swp
671
672
673   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
674      !!---------------------------------------------------------------------
675      !!                  ***  ROUTINE dom_vvl__interpol  ***
676      !!
677      !! ** Purpose :   interpolate scale factors from one grid point to another
678      !!
679      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
680      !!                - horizontal interpolation: grid cell surface averaging
681      !!                - vertical interpolation: simple averaging
682      !!----------------------------------------------------------------------
683      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
684      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
685      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
686      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
687      !
688      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
689      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
690      !!----------------------------------------------------------------------
691      !
692      IF(ln_wd_il) THEN
693        zlnwd = 1.0_wp
694      ELSE
695        zlnwd = 0.0_wp
696      END IF
697      !
698      SELECT CASE ( pout )    !==  type of interpolation  ==!
699         !
700      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
701         DO jk = 1, jpk
702            DO jj = 1, jpjm1
703               DO ji = 1, fs_jpim1   ! vector opt.
704                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
705                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
706                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
707               END DO
708            END DO
709         END DO
710         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
711         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
712         !
713      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
714         DO jk = 1, jpk
715            DO jj = 1, jpjm1
716               DO ji = 1, fs_jpim1   ! vector opt.
717                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
718                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
719                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
720               END DO
721            END DO
722         END DO
723         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
724         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
725         !
726      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
727         DO jk = 1, jpk
728            DO jj = 1, jpjm1
729               DO ji = 1, fs_jpim1   ! vector opt.
730                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
731                     &                       *    r1_e1e2f(ji,jj)                                                  &
732                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
733                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
734               END DO
735            END DO
736         END DO
737         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp )
738         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
739         !
740      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
741         !
742         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
743         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
744!!gm BUG? use here wmask in case of ISF ?  to be checked
745         DO jk = 2, jpk
746            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
747               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
748               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
749               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
750         END DO
751         !
752      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
753         !
754         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
755         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
756!!gm BUG? use here wumask in case of ISF ?  to be checked
757         DO jk = 2, jpk
758            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
759               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
760               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
761               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
762         END DO
763         !
764      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
765         !
766         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
767         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
768!!gm BUG? use here wvmask in case of ISF ?  to be checked
769         DO jk = 2, jpk
770            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
771               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
772               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
773               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
774         END DO
775      END SELECT
776      !
777   END SUBROUTINE dom_vvl_interpol
778
779
780   SUBROUTINE dom_vvl_rst( kt, cdrw )
781      !!---------------------------------------------------------------------
782      !!                   ***  ROUTINE dom_vvl_rst  ***
783      !!                     
784      !! ** Purpose :   Read or write VVL file in restart file
785      !!
786      !! ** Method  :   use of IOM library
787      !!                if the restart does not contain vertical scale factors,
788      !!                they are set to the _0 values
789      !!                if the restart does not contain vertical scale factors increments (z_tilde),
790      !!                they are set to 0.
791      !!----------------------------------------------------------------------
792      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
793      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
794      !
795      INTEGER ::   ji, jj, jk
796      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
797      !!----------------------------------------------------------------------
798      !
799      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
800         !                                   ! ===============
801         IF( ln_rstart ) THEN                   !* Read the restart file
802            CALL rst_read_open                  !  open the restart file if necessary
803            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , sshn, ldxios = lrxios    )
804            !
805            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
806            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
807            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
808            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
809            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
810            !                             ! --------- !
811            !                             ! all cases !
812            !                             ! --------- !
813            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
814               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
815               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
816               ! needed to restart if land processor not computed
817               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t_b and e3t_n found in restart files'
818               WHERE ( tmask(:,:,:) == 0.0_wp ) 
819                  e3t_n(:,:,:) = e3t_0(:,:,:)
820                  e3t_b(:,:,:) = e3t_0(:,:,:)
821               END WHERE
822               IF( neuler == 0 ) THEN
823                  e3t_b(:,:,:) = e3t_n(:,:,:)
824               ENDIF
825            ELSE IF( id1 > 0 ) THEN
826               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart files'
827               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
828               IF(lwp) write(numout,*) 'neuler is forced to 0'
829               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
830               e3t_n(:,:,:) = e3t_b(:,:,:)
831               neuler = 0
832            ELSE IF( id2 > 0 ) THEN
833               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_b not found in restart files'
834               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
835               IF(lwp) write(numout,*) 'neuler is forced to 0'
836               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t_n(:,:,:), ldxios = lrxios )
837               e3t_b(:,:,:) = e3t_n(:,:,:)
838               neuler = 0
839            ELSE
840               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t_n not found in restart file'
841               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
842               IF(lwp) write(numout,*) 'neuler is forced to 0'
843               DO jk = 1, jpk
844                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:) ) &
845                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
846                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
847               END DO
848               e3t_b(:,:,:) = e3t_n(:,:,:)
849               neuler = 0
850            ENDIF
851            !                             ! ----------- !
852            IF( ln_vvl_zstar ) THEN       ! z_star case !
853               !                          ! ----------- !
854               IF( MIN( id3, id4 ) > 0 ) THEN
855                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
856               ENDIF
857               !                          ! ----------------------- !
858            ELSE                          ! z_tilde and layer cases !
859               !                          ! ----------------------- !
860               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
861                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
862                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
863               ELSE                            ! one at least array is missing
864                  tilde_e3t_b(:,:,:) = 0.0_wp
865                  tilde_e3t_n(:,:,:) = 0.0_wp
866               ENDIF
867               !                          ! ------------ !
868               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
869                  !                       ! ------------ !
870                  IF( id5 > 0 ) THEN  ! required array exists
871                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
872                  ELSE                ! array is missing
873                     hdiv_lf(:,:,:) = 0.0_wp
874                  ENDIF
875               ENDIF
876            ENDIF
877            !
878         ELSE                                   !* Initialize at "rest"
879            !
880
881            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
882               !
883               IF( cn_cfg == 'wad' ) THEN
884                  ! Wetting and drying test case
885                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, sshb  )
886                  tsn  (:,:,:,:) = tsb (:,:,:,:)       ! set now values from to before ones
887                  sshn (:,:)     = sshb(:,:)
888                  un   (:,:,:)   = ub  (:,:,:)
889                  vn   (:,:,:)   = vb  (:,:,:)
890               ELSE
891                  ! if not test case
892                  sshn(:,:) = -ssh_ref
893                  sshb(:,:) = -ssh_ref
894
895                  DO jj = 1, jpj
896                     DO ji = 1, jpi
897                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
898
899                           sshb(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
900                           sshn(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
901                           ssha(ji,jj) = rn_wdmin1 - (ht_0(ji,jj) )
902                        ENDIF
903                     ENDDO
904                  ENDDO
905               ENDIF !If test case else
906
907               ! Adjust vertical metrics for all wad
908               DO jk = 1, jpk
909                  e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshn(:,:)  ) &
910                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
911                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
912               END DO
913               e3t_b(:,:,:) = e3t_n(:,:,:)
914
915               DO ji = 1, jpi
916                  DO jj = 1, jpj
917                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
918                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
919                     ENDIF
920                  END DO
921               END DO 
922               !
923            ELSE
924               !
925               ! Just to read set ssh in fact, called latter once vertical grid
926               ! is set up:
927!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, sshb  )
928!               !
929!               DO jk=1,jpk
930!                  e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( ht_0(:,:) + sshb(:,:) ) &
931!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
932!               END DO
933!               e3t_n(:,:,:) = e3t_b(:,:,:)
934                sshn(:,:)=0._wp
935                e3t_n(:,:,:)=e3t_0(:,:,:)
936                e3t_b(:,:,:)=e3t_0(:,:,:)
937               !
938            END IF           ! end of ll_wd edits
939
940            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
941               tilde_e3t_b(:,:,:) = 0._wp
942               tilde_e3t_n(:,:,:) = 0._wp
943               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
944            END IF
945         ENDIF
946         !
947      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
948         !                                   ! ===================
949         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
950         IF( lwxios ) CALL iom_swap(      cwxios_context          )
951         !                                           ! --------- !
952         !                                           ! all cases !
953         !                                           ! --------- !
954         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t_b(:,:,:), ldxios = lwxios )
955         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t_n(:,:,:), ldxios = lwxios )
956         !                                           ! ----------------------- !
957         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
958            !                                        ! ----------------------- !
959            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
960            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
961         END IF
962         !                                           ! -------------!   
963         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
964            !                                        ! ------------ !
965            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
966         ENDIF
967         !
968         IF( lwxios ) CALL iom_swap(      cxios_context          )
969      ENDIF
970      !
971   END SUBROUTINE dom_vvl_rst
972
973
974   SUBROUTINE dom_vvl_ctl
975      !!---------------------------------------------------------------------
976      !!                  ***  ROUTINE dom_vvl_ctl  ***
977      !!               
978      !! ** Purpose :   Control the consistency between namelist options
979      !!                for vertical coordinate
980      !!----------------------------------------------------------------------
981      INTEGER ::   ioptio, ios
982      !!
983      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
984         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
985         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
986      !!----------------------------------------------------------------------
987      !
988      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
989      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
990901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
991      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
992      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
993902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
994      IF(lwm) WRITE ( numond, nam_vvl )
995      !
996      IF(lwp) THEN                    ! Namelist print
997         WRITE(numout,*)
998         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
999         WRITE(numout,*) '~~~~~~~~~~~'
1000         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
1001         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1002         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1003         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
1004         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1005         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1006         WRITE(numout,*) '      !'
1007         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
1008         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
1009         IF( ln_vvl_ztilde_as_zstar ) THEN
1010            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1011            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
1012            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
1013            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
1014            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1015            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
1016         ELSE
1017            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
1018            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
1019         ENDIF
1020         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
1021      ENDIF
1022      !
1023      ioptio = 0                      ! Parameter control
1024      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1025      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1026      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1027      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1028      !
1029      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1030      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1031      !
1032      IF(lwp) THEN                   ! Print the choice
1033         WRITE(numout,*)
1034         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1035         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1036         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1037         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1038      ENDIF
1039      !
1040#if defined key_agrif
1041      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1042#endif
1043      !
1044   END SUBROUTINE dom_vvl_ctl
1045
1046   !!======================================================================
1047END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.