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/2020/dev_12905_xios_restart/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2020/dev_12905_xios_restart/tests/CANAL/MY_SRC/domvvl.F90 @ 13727

Last change on this file since 13727 was 13727, checked in by andmirek, 4 years ago

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