source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domvvl.F90 @ 12680

Last change on this file since 12680 was 12680, checked in by techene, 8 months ago

dynatfQCO.F90, stepLF.F90 : fixed (remove pe3. from dyn_atf_qco input arguments), all : remove e3. tables and include gurvan's feedbacks

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