source: NEMO/trunk/tests/VORTEX/MY_SRC/domvvl.F90 @ 12740

Last change on this file since 12740 was 12740, checked in by smasson, 7 months ago

trunk: update/debug of tests and C1D, see #2442

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