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 NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/DOM – NEMO

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

Last change on this file since 11553 was 11553, checked in by mathiot, 5 years ago

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