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_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM – NEMO

source: NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/DOM/domvvl.F90 @ 15475

Last change on this file since 15475 was 14058, checked in by ayoung, 4 years ago

Updating to trunk revision 14057. No conflicts since last sette test. Ticket #2480.
  • Property svn:keywords set to Id
File size: 52.8 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   !!             -   !  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 )                     ! vertical sum
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 )   ! a - first derivative: diffusive fluxes
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 )             ! b - correction for last oceanic u-v points
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 )   ! c - second derivative: divergence of diffusive fluxes
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         ! 4 - Time stepping of baroclinic scale factors
425         ! ---------------------------------------------
426         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
427         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
428
429         ! Maximum deformation control
430         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
431         ALLOCATE( ze3t(jpi,jpj,jpk), llmsk(jpi,jpj,jpk) )
432         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
433            ze3t(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) / e3t_0(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
434         END_3D
435         !
436         llmsk(   1:Nis1,:,:) = .FALSE.   ! exclude halos from the checked region
437         llmsk(Nie1: jpi,:,:) = .FALSE.
438         llmsk(:,   1:Njs1,:) = .FALSE.
439         llmsk(:,Nje1: jpj,:) = .FALSE.
440         !
441         llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp                  ! define only the inner domain
442         z_tmax = MAXVAL( ze3t(:,:,:), mask = llmsk )   ;   CALL mpp_max( 'domvvl', z_tmax )   ! max over the global domain
443         z_tmin = MINVAL( ze3t(:,:,:), mask = llmsk )   ;   CALL mpp_min( 'domvvl', z_tmin )   ! min over the global domain
444         ! - ML - test: for the moment, stop simulation for too large e3_t variations
445         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
446            CALL mpp_maxloc( 'domvvl', ze3t, llmsk, z_tmax, ijk_max )
447            CALL mpp_minloc( 'domvvl', ze3t, llmsk, z_tmin, ijk_min )
448            IF (lwp) THEN
449               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
450               WRITE(numout, *) 'at i, j, k=', ijk_max
451               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
452               WRITE(numout, *) 'at i, j, k=', ijk_min           
453               CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high')
454            ENDIF
455         ENDIF
456         DEALLOCATE( ze3t, llmsk )
457         ! - ML - end test
458         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
459         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
460         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
461
462         !
463         ! "tilda" change in the after scale factor
464         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
465         DO jk = 1, jpkm1
466            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
467         END DO
468         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
469         ! ==================================================================================
470         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
471         ! - ML - baroclinicity error should be better treated in the future
472         !        i.e. locally and not spread over the water column.
473         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
474         zht(:,:) = 0.
475         DO jk = 1, jpkm1
476            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
477         END DO
478         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
479         DO jk = 1, jpkm1
480            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
481         END DO
482
483      ENDIF
484
485      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
486      !                                           ! ---baroclinic part--------- !
487         DO jk = 1, jpkm1
488            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
489         END DO
490      ENDIF
491
492      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
493         !
494         IF( lwp ) WRITE(numout, *) 'kt =', kt
495         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
496            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
497            CALL mpp_max( 'domvvl', z_tmax )                             ! max over the global domain
498            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
499         END IF
500         !
501         zht(:,:) = 0.0_wp
502         DO jk = 1, jpkm1
503            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
504         END DO
505         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
506         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
507         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
508         !
509         zht(:,:) = 0.0_wp
510         DO jk = 1, jpkm1
511            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
512         END DO
513         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
514         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
515         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
516         !
517         zht(:,:) = 0.0_wp
518         DO jk = 1, jpkm1
519            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
520         END DO
521         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
522         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
523         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
524         !
525         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) )
526         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
527         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
528         !
529         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) )
530         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
531         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
532         !
533         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) )
534         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
535         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
536      END IF
537
538      ! *********************************** !
539      ! After scale factors at u- v- points !
540      ! *********************************** !
541
542      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
543      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
544
545      ! *********************************** !
546      ! After depths at u- v points         !
547      ! *********************************** !
548
549      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
550      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
551      DO jk = 2, jpkm1
552         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
553         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
554      END DO
555      !                                        ! Inverse of the local depth
556!!gm BUG ?  don't understand the use of umask_i here .....
557      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
558      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
559      !
560      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
561      !
562   END SUBROUTINE dom_vvl_sf_nxt
563
564
565   SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa )
566      !!----------------------------------------------------------------------
567      !!                ***  ROUTINE dom_vvl_sf_update  ***
568      !!                   
569      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
570      !!               compute all depths and related variables for next time step
571      !!               write outputs and restart file
572      !!
573      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
574      !!               - reconstruct scale factor at other grid points (interpolate)
575      !!               - recompute depths and water height fields
576      !!
577      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step
578      !!               - Recompute:
579      !!                    e3(u/v)_b       
580      !!                    e3w(:,:,:,Kmm)           
581      !!                    e3(u/v)w_b     
582      !!                    e3(u/v)w_n     
583      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
584      !!                    h(u/v) and h(u/v)r
585      !!
586      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
587      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
588      !!----------------------------------------------------------------------
589      INTEGER, INTENT( in ) ::   kt              ! time step
590      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
591      !
592      INTEGER  ::   ji, jj, jk   ! dummy loop indices
593      REAL(wp) ::   zcoef        ! local scalar
594      !!----------------------------------------------------------------------
595      !
596      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
597      !
598      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update')
599      !
600      IF( kt == nit000 )   THEN
601         IF(lwp) WRITE(numout,*)
602         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
603         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
604      ENDIF
605      !
606      ! Time filter and swap of scale factors
607      ! =====================================
608      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
609      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
610         IF( l_1st_euler ) THEN
611            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
612         ELSE
613            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
614            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
615         ENDIF
616         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
617      ENDIF
618
619      ! Compute all missing vertical scale factor and depths
620      ! ====================================================
621      ! Horizontal scale factor interpolations
622      ! --------------------------------------
623      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
624      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
625     
626      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
627     
628      ! Vertical scale factor interpolations
629      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
630      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
631      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
632      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
633      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
634      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
635
636      ! t- and w- points depth (set the isf depth as it is in the initial step)
637      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
638      gdepw(:,:,1,Kmm) = 0.0_wp
639      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
640      DO_3D( 1, 1, 1, 1, 2, jpk )
641        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
642                                                           ! 1 for jk = mikt
643         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
644         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
645         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
646             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
647         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
648      END_3D
649
650      ! Local depth and Inverse of the local depth of the water
651      ! -------------------------------------------------------
652      !
653      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
654      DO jk = 2, jpkm1
655         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
656      END DO
657
658      ! write restart file
659      ! ==================
660      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
661      !
662      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update')
663      !
664   END SUBROUTINE dom_vvl_sf_update
665
666
667   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
668      !!---------------------------------------------------------------------
669      !!                  ***  ROUTINE dom_vvl__interpol  ***
670      !!
671      !! ** Purpose :   interpolate scale factors from one grid point to another
672      !!
673      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
674      !!                - horizontal interpolation: grid cell surface averaging
675      !!                - vertical interpolation: simple averaging
676      !!----------------------------------------------------------------------
677      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
678      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
679      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
680      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
681      !
682      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
683      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
684      !!----------------------------------------------------------------------
685      !
686      IF(ln_wd_il) THEN
687        zlnwd = 1.0_wp
688      ELSE
689        zlnwd = 0.0_wp
690      END IF
691      !
692      SELECT CASE ( pout )    !==  type of interpolation  ==!
693         !
694      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
695         DO_3D( 1, 0, 1, 0, 1, jpk )
696            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
697               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
698               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
699         END_3D
700         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
701         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
702         !
703      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
704         DO_3D( 1, 0, 1, 0, 1, jpk )
705            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
706               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
707               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
708         END_3D
709         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
710         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
711         !
712      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
713         DO_3D( 1, 0, 1, 0, 1, jpk )
714            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
715               &                       *    r1_e1e2f(ji,jj)                                                  &
716               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
717               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
718         END_3D
719         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
720         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
721         !
722      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
723         !
724         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
725         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
726!!gm BUG? use here wmask in case of ISF ?  to be checked
727         DO jk = 2, jpk
728            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
729               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
730               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
731               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
732         END DO
733         !
734      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
735         !
736         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
737         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
738!!gm BUG? use here wumask in case of ISF ?  to be checked
739         DO jk = 2, jpk
740            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
741               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
742               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
743               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
744         END DO
745         !
746      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
747         !
748         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
749         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
750!!gm BUG? use here wvmask in case of ISF ?  to be checked
751         DO jk = 2, jpk
752            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
753               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
754               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
755               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
756         END DO
757      END SELECT
758      !
759   END SUBROUTINE dom_vvl_interpol
760
761
762   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
763      !!---------------------------------------------------------------------
764      !!                   ***  ROUTINE dom_vvl_rst  ***
765      !!                     
766      !! ** Purpose :   Read or write VVL file in restart file
767      !!
768      !! ** Method  : * restart comes from a linear ssh simulation :
769      !!                   an attempt to read e3t_n stops simulation
770      !!              * restart comes from a z-star, z-tilde, or layer :
771      !!                   read e3t_n and e3t_b
772      !!              * restart comes from a z-star :
773      !!                   set tilde_e3t_n, tilde_e3t_n, and hdiv_lf to 0
774      !!              * restart comes from layer :
775      !!                   read tilde_e3t_n and tilde_e3t_b
776      !!                   set hdiv_lf to 0
777      !!              * restart comes from a z-tilde:
778      !!                   read tilde_e3t_n, tilde_e3t_b, and hdiv_lf
779      !!
780      !!              NB: if l_1st_euler = T (ln_1st_euler or ssh_b not found)
781      !!                   Kbb fields set to Kmm ones
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      ! dummy loop indices
788      INTEGER ::   id3, id4, id5   ! local integers
789      !!----------------------------------------------------------------------
790      !
791      !                                      !=====================!
792      IF( TRIM(cdrw) == 'READ' ) THEN        !  Read / initialise  !
793         !                                   !=====================!
794         !
795         IF( ln_rstart ) THEN                   !==  Read the restart file  ==!
796            !
797            CALL rst_read_open                                          !*  open the restart file if necessary
798            !                                         ! --------- !
799            !                                         ! all cases !
800            !                                         ! --------- !
801            !
802            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )  !*  check presence
803            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
804            id5 = iom_varid( numror, 'hdiv_lf'    , ldstop = .FALSE. )
805            !
806            !                                                           !*  scale factors
807            IF(lwp) WRITE(numout,*)    '          Kmm scale factor read in the restart file'
808            CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm) )
809            WHERE ( tmask(:,:,:) == 0.0_wp ) 
810               e3t(:,:,:,Kmm) = e3t_0(:,:,:)
811            END WHERE
812            IF( l_1st_euler ) THEN                       ! euler
813               IF(lwp) WRITE(numout,*) '          Euler first time step : e3t(Kbb) = e3t(Kmm)'
814               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
815            ELSE                                         ! leap frog
816               IF(lwp) WRITE(numout,*) '          Kbb scale factor read in the restart file'
817               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb) )
818               WHERE ( tmask(:,:,:) == 0.0_wp ) 
819                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
820               END WHERE
821            ENDIF
822            !                                         ! ------------ !
823            IF( ln_vvl_zstar ) THEN                   !  z_star case !
824               !                                      ! ------------ !
825               IF( MIN( id3, id4 ) > 0 ) THEN
826                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
827               ENDIF
828               !                                      ! ------------------------ !
829            ELSE                                      !  z_tilde and layer cases !
830               !                                      ! ------------------------ !
831               !
832               IF( id4 > 0 ) THEN                                       !*  scale factor increments
833                  IF(lwp) WRITE(numout,*)    '          Kmm scale factor increments read in the restart file'
834                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_n', tilde_e3t_n(:,:,:) )
835                  IF( l_1st_euler ) THEN                 ! euler
836                     IF(lwp) WRITE(numout,*) '          Euler first time step : tilde_e3t(Kbb) = tilde_e3t(Kmm)'
837                     tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
838                  ELSE                                   ! leap frog
839                     IF(lwp) WRITE(numout,*) '          Kbb scale factor increments read in the restart file'
840                     CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:) )
841                  ENDIF
842               ELSE
843                  tilde_e3t_b(:,:,:) = 0.0_wp
844                  tilde_e3t_n(:,:,:) = 0.0_wp
845               ENDIF
846               !                                      ! ------------ !
847               IF( ln_vvl_ztilde ) THEN               ! z_tilde case !
848                  !                                   ! ------------ !
849                  IF( id5 > 0 ) THEN  ! required array exists
850                     CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:) )
851                  ELSE                ! array is missing
852                     hdiv_lf(:,:,:) = 0.0_wp 
853                  ENDIF
854               ENDIF
855            ENDIF
856            !
857         ELSE                                   !==  Initialize at "rest" with ssh  ==!
858            !
859            DO jk = 1, jpk
860               e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( 1._wp + ssh(:,:,Kmm) * r1_ht_0(:,:) * tmask(:,:,jk) )
861            END DO
862            e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
863            !
864            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
865               tilde_e3t_b(:,:,:) = 0._wp
866               tilde_e3t_n(:,:,:) = 0._wp
867               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
868            ENDIF
869         ENDIF
870         !                                       !=======================!
871      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN       !  Create restart file  !
872         !                                       !=======================!
873         !
874         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
875         !                                           ! --------- !
876         !                                           ! all cases !
877         !                                           ! --------- !
878         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb) )
879         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm) )
880         !                                           ! ----------------------- !
881         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
882            !                                        ! ----------------------- !
883            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:))
884            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:))
885         END IF
886         !                                           ! -------------!   
887         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
888            !                                        ! ------------ !
889            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:))
890         ENDIF
891         !
892      ENDIF
893      !
894   END SUBROUTINE dom_vvl_rst
895
896
897   SUBROUTINE dom_vvl_ctl
898      !!---------------------------------------------------------------------
899      !!                  ***  ROUTINE dom_vvl_ctl  ***
900      !!               
901      !! ** Purpose :   Control the consistency between namelist options
902      !!                for vertical coordinate
903      !!----------------------------------------------------------------------
904      INTEGER ::   ioptio, ios
905      !!
906      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
907         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
908         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
909      !!----------------------------------------------------------------------
910      !
911      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
912901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
913      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
914902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
915      IF(lwm) WRITE ( numond, nam_vvl )
916      !
917      IF(lwp) THEN                    ! Namelist print
918         WRITE(numout,*)
919         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
920         WRITE(numout,*) '~~~~~~~~~~~'
921         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
922         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
923         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
924         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
925         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
926         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
927         WRITE(numout,*) '      !'
928         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
929         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
930         IF( ln_vvl_ztilde_as_zstar ) THEN
931            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
932            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
933            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
934            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
935            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
936            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rn_Dt'
937         ELSE
938            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
939            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
940         ENDIF
941         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
942      ENDIF
943      !
944      ioptio = 0                      ! Parameter control
945      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
946      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
947      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
948      IF( ln_vvl_layer           )   ioptio = ioptio + 1
949      !
950      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
951      !
952      IF(lwp) THEN                   ! Print the choice
953         WRITE(numout,*)
954         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
955         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
956         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
957         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
958      ENDIF
959      !
960#if defined key_agrif
961      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
962#endif
963      !
964   END SUBROUTINE dom_vvl_ctl
965
966#endif
967
968   !!======================================================================
969END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.