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_r13312_AGRIF-03-04_jchanut_vinterp_tstep/tests/VORTEX/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/tests/VORTEX/MY_SRC/domvvl.F90 @ 13338

Last change on this file since 13338 was 13338, checked in by jchanut, 4 years ago

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