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_r13296_HPC-07_mocavero_mpi3/src/SWE – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/SWE/domvvl.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

Add neighborhood collectives calls in the NEMO src - ticket #2496
File size: 69.6 KB
Line 
1
2MODULE domvvl
3   !!======================================================================
4   !!                       ***  MODULE domvvl   ***
5   !! Ocean :
6   !!======================================================================
7   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
8   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
9   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates
10   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
11   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping
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!!stoops
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!!st
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   PUBLIC  dom_vvl_interpol_st! called by dynnxt.F90
76   PUBLIC  dom_vvl_sf_nxt_st  ! called by step.F90
77   PUBLIC  dom_vvl_sf_update_st
78!!st
79   
80   !! * Substitutions
81#  include "do_loop_substitute.h90"
82   !!----------------------------------------------------------------------
83   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
84   !! $Id: domvvl.F90 12614 2020-03-26 14:59:52Z gm $
85   !! Software governed by the CeCILL license (see ./LICENSE)
86   !!----------------------------------------------------------------------
87CONTAINS
88
89   INTEGER FUNCTION dom_vvl_alloc()
90      !!----------------------------------------------------------------------
91      !!                ***  FUNCTION dom_vvl_alloc  ***
92      !!----------------------------------------------------------------------
93      IF( ln_vvl_zstar )   dom_vvl_alloc = 0
94      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
95         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
96            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
97            &      STAT = dom_vvl_alloc        )
98         CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
99         IF( dom_vvl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
100         un_td = 0._wp
101         vn_td = 0._wp
102      ENDIF
103      IF( ln_vvl_ztilde ) THEN
104         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
105         CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
106         IF( dom_vvl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
107      ENDIF
108      !
109   END FUNCTION dom_vvl_alloc
110
111
112   SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa )
113      !!----------------------------------------------------------------------
114      !!                ***  ROUTINE dom_vvl_init  ***
115      !!                   
116      !! ** Purpose :  Initialization of all scale factors, depths
117      !!               and water column heights
118      !!
119      !! ** Method  :  - use restart file and/or initialize
120      !!               - interpolate scale factors
121      !!
122      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
123      !!              - Regrid: e3[u/v](:,:,:,Kmm)
124      !!                        e3[u/v](:,:,:,Kmm)       
125      !!                        e3w(:,:,:,Kmm)           
126      !!                        e3[u/v]w_b
127      !!                        e3[u/v]w_n     
128      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
129      !!              - h(t/u/v)_0
130      !!              - frq_rst_e3t and frq_rst_hdv
131      !!
132      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
133      !!----------------------------------------------------------------------
134      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
135      !
136      IF(lwp) WRITE(numout,*)
137      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
138      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
139      !
140      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
141      !
142      !                    ! Allocate module arrays
143      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
144      !
145      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
146      CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' )
147      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
148      !
149      CALL dom_vvl_zgr_st(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column
150      !
151   END SUBROUTINE dom_vvl_init
152   !
153   SUBROUTINE dom_vvl_zgr(Kbb, Kmm, Kaa)
154      !!----------------------------------------------------------------------
155      !!                ***  ROUTINE dom_vvl_init  ***
156      !!                   
157      !! ** Purpose :  Interpolation of all scale factors,
158      !!               depths and water column heights
159      !!
160      !! ** Method  :  - interpolate scale factors
161      !!
162      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
163      !!              - Regrid: e3(u/v)_n
164      !!                        e3(u/v)_b       
165      !!                        e3w_n           
166      !!                        e3(u/v)w_b     
167      !!                        e3(u/v)w_n     
168      !!                        gdept_n, gdepw_n and gde3w_n
169      !!              - h(t/u/v)_0
170      !!              - frq_rst_e3t and frq_rst_hdv
171      !!
172      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
173      !!----------------------------------------------------------------------
174      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
175      !!----------------------------------------------------------------------
176      INTEGER ::   ji, jj, jk
177      INTEGER ::   ii0, ii1, ij0, ij1
178      REAL(wp)::   zcoef
179      !!----------------------------------------------------------------------
180      !
181      !                    !== Set of all other vertical scale factors  ==!  (now and before)
182      !                                ! Horizontal interpolation of e3t
183      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U
184      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
185      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V
186      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
187      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F
188      !                                ! Vertical interpolation of e3t,u,v
189      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W
190      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3w (:,:,:,Kbb), 'W'  )
191      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW
192      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
193      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW
194      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
195
196      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
197      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
198      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
199      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
200      !
201      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
202      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)       ! reference to the ocean surface (used for MLD and light penetration)
203      gdepw(:,:,1,Kmm) = 0.0_wp
204      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
205      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
206      gdepw(:,:,1,Kbb) = 0.0_wp
207      DO_3D( 1, 1, 1, 1, 2, jpk )
208         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
209         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
210         !                             ! 0.5 where jk = mikt     
211!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ??
212         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
213         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
214         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
215            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm)) 
216         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
217         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
218         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  &
219            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb)) 
220      END_3D
221      !
222      !                    !==  thickness of the water column  !!   (ocean portion only)
223      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
224      hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1)
225      hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1)
226      hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1)
227      hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1)
228      DO jk = 2, jpkm1
229         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
230         hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
231         hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
232         hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
233         hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
234      END DO
235      !
236      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
237      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
238      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
239      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
240      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
241
242      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
243      IF( ln_vvl_ztilde ) THEN
244!!gm : idea: add here a READ in a file of custumized restoring frequency
245         !                                   ! Values in days provided via the namelist
246         !                                   ! use rsmall to avoid possible division by zero errors with faulty settings
247         frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
248         frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
249         !
250         IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile
251            frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings
252            frq_rst_hdv(:,:) = 1._wp / rn_Dt
253         ENDIF
254         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
255            DO_2D( 1, 1, 1, 1 )
256!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default
257               IF( ABS(gphit(ji,jj)) >= 6.) THEN
258                  ! values outside the equatorial band and transition zone (ztilde)
259                  frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
260                  frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
261               ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star
262                  ! values inside the equatorial band (ztilde as zstar)
263                  frq_rst_e3t(ji,jj) =  0.0_wp
264                  frq_rst_hdv(ji,jj) =  1.0_wp / rn_Dt
265               ELSE                                      ! transition band (2.5 to 6 degrees N/S)
266                  !                                      ! (linearly transition from z-tilde to z-star)
267                  frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
268                     &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
269                     &                                          * 180._wp / 3.5_wp ) )
270                  frq_rst_hdv(ji,jj) = (1.0_wp / rn_Dt)                                &
271                     &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rn_Dt) )*0.5_wp   &
272                     &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
273                     &                                          * 180._wp / 3.5_wp ) )
274               ENDIF
275            END_2D
276            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN
277               IF( nn_cfg == 3 ) THEN   ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2
278                  ii0 = 103 + nn_hls - 1   ;   ii1 = 111 + nn_hls - 1     
279                  ij0 = 128 + nn_hls       ;   ij1 = 135 + nn_hls
280                  frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  0.0_wp
281                  frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) =  1.e0_wp / rn_Dt
282               ENDIF
283            ENDIF
284         ENDIF
285      ENDIF
286      !
287      IF(lwxios) THEN
288! define variables in restart file when writing with XIOS
289         CALL iom_set_rstw_var_active('e3t_b')
290         CALL iom_set_rstw_var_active('e3t_n')
291         !                                           ! ----------------------- !
292         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
293            !                                        ! ----------------------- !
294            CALL iom_set_rstw_var_active('tilde_e3t_b')
295            CALL iom_set_rstw_var_active('tilde_e3t_n')
296         END IF
297         !                                           ! -------------!   
298         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
299            !                                        ! ------------ !
300            CALL iom_set_rstw_var_active('hdiv_lf')
301         ENDIF
302         !
303      ENDIF
304      !
305   END SUBROUTINE dom_vvl_zgr
306
307   
308   SUBROUTINE dom_vvl_zgr_st(Kbb, Kmm, Kaa)
309      !!----------------------------------------------------------------------
310      !!                ***  ROUTINE dom_vvl_init  ***
311      !!                   
312      !! ** Purpose :  Interpolation of all scale factors,
313      !!               depths and water column heights
314      !!
315      !! ** Method  :  - interpolate scale factors
316      !!
317      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
318      !!              - Regrid: e3(u/v)_n
319      !!                        e3(u/v)_b       
320      !!                        e3w_n           
321      !!                        e3(u/v)w_b     
322      !!                        e3(u/v)w_n     
323      !!                        gdept_n, gdepw_n and gde3w_n
324      !!              - h(t/u/v)_0
325      !!              - frq_rst_e3t and frq_rst_hdv
326      !!
327      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
328      !!----------------------------------------------------------------------
329      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
330      !!----------------------------------------------------------------------
331      INTEGER ::   ji, jj, jk
332      INTEGER ::   ii0, ii1, ij0, ij1
333      REAL(wp)::   zcoef
334      !!----------------------------------------------------------------------
335      !
336      !                    !== Set of all other vertical scale factors  ==!  (now and before)
337      !                                ! Horizontal interpolation of e3t
338      CALL dom_vvl_interpol_st( r3u(:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U
339      CALL dom_vvl_interpol_st( r3u(:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
340      CALL dom_vvl_interpol_st( r3v(:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V
341      CALL dom_vvl_interpol_st( r3v(:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
342      CALL dom_vvl_interpol_st( r3f(:,:), e3f(:,:,:), 'F' )    ! from U to F
343      !                                ! Vertical interpolation of e3t,u,v
344      CALL dom_vvl_interpol_st( r3t(:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W
345      CALL dom_vvl_interpol_st( r3t(:,:,Kbb), e3w (:,:,:,Kbb), 'W'  )
346      CALL dom_vvl_interpol_st( r3u(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW
347      CALL dom_vvl_interpol_st( r3u(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
348      CALL dom_vvl_interpol_st( r3v(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW
349      CALL dom_vvl_interpol_st( r3v(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
350
351      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
352      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
353      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
354      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
355      !
356      DO_3D( 1, 1, 1, 1, 1, jpk )
357         gdepw(ji,jj,jk,Kmm) =  gdepw_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm))
358         gdept(ji,jj,jk,Kmm) =  gdept_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm)) 
359         gde3w(ji,jj,jk    ) =  gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
360         gdepw(ji,jj,jk,Kbb) =  gdepw_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kbb))
361         gdept(ji,jj,jk,Kbb) =  gdept_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kbb)) 
362      END_3D     
363      !
364      !                    !==  thickness of the water column  !!   (ocean portion only)
365      ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
366      hu(:,:,Kbb) = (hu_0(:,:)*(1._wp+r3u(:,:,Kbb)))
367      hv(:,:,Kbb) = (hv_0(:,:)*(1._wp+r3v(:,:,Kbb)))
368      hu(:,:,Kbb) = (hu_0(:,:)*(1._wp+r3u(:,:,Kmm)))
369      hv(:,:,Kbb) = (hv_0(:,:)*(1._wp+r3v(:,:,Kmm)))
370      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
371      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
372      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
373      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
374      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 
375      !
376      IF(lwxios) THEN
377! define variables in restart file when writing with XIOS
378         CALL iom_set_rstw_var_active('e3t_b')
379         CALL iom_set_rstw_var_active('e3t_n')
380         !
381      ENDIF
382      !
383   END SUBROUTINE dom_vvl_zgr_st
384   
385
386   SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 
387      !!----------------------------------------------------------------------
388      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
389      !!                   
390      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
391      !!                 tranxt and dynspg routines
392      !!
393      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
394      !!               - z_tilde_case: after scale factor increment =
395      !!                                    high frequency part of horizontal divergence
396      !!                                  + retsoring towards the background grid
397      !!                                  + thickness difusion
398      !!                               Then repartition of ssh INCREMENT proportionnaly
399      !!                               to the "baroclinic" level thickness.
400      !!
401      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
402      !!               - tilde_e3t_a: after increment of vertical scale factor
403      !!                              in z_tilde case
404      !!               - e3(t/u/v)_a
405      !!
406      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
407      !!----------------------------------------------------------------------
408      INTEGER, INTENT( in )           ::   kt             ! time step
409      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
410      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
411      !
412      INTEGER                ::   ji, jj, jk            ! dummy loop indices
413      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
414      REAL(wp)               ::   z_tmin, z_tmax        ! local scalars
415      LOGICAL                ::   ll_do_bclinic         ! local logical
416      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
417      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ze3t
418      LOGICAL , DIMENSION(:,:,:), ALLOCATABLE ::   llmsk
419      !!----------------------------------------------------------------------
420      !
421      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
422      !
423      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
424      !
425      IF( kt == nit000 ) THEN
426         IF(lwp) WRITE(numout,*)
427         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
428         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
429      ENDIF
430
431      ll_do_bclinic = .TRUE.
432      IF( PRESENT(kcall) ) THEN
433         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
434      ENDIF
435
436      ! ******************************* !
437      ! After acale factors at t-points !
438      ! ******************************* !
439      !                                                ! --------------------------------------------- !
440      !                                                ! z_star coordinate and barotropic z-tilde part !
441      !                                                ! --------------------------------------------- !
442      !
443      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
444      DO jk = 1, jpkm1
445         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
446         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
447      END DO
448      !
449      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
450         !                                                               ! ------baroclinic part------ !
451         ! I - initialization
452         ! ==================
453
454         ! 1 - barotropic divergence
455         ! -------------------------
456         zhdiv(:,:) = 0._wp
457         zht(:,:)   = 0._wp
458         DO jk = 1, jpkm1
459            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
460            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
461         END DO
462         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
463
464         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
465         ! --------------------------------------------------
466         IF( ln_vvl_ztilde ) THEN
467            IF( kt > nit000 ) THEN
468               DO jk = 1, jpkm1
469                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:)   &
470                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
471               END DO
472            ENDIF
473         ENDIF
474
475         ! II - after z_tilde increments of vertical scale factors
476         ! =======================================================
477         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
478
479         ! 1 - High frequency divergence term
480         ! ----------------------------------
481         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
482            DO jk = 1, jpkm1
483               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
484            END DO
485         ELSE                         ! layer case
486            DO jk = 1, jpkm1
487               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
488            END DO
489         ENDIF
490
491         ! 2 - Restoring term (z-tilde case only)
492         ! ------------------
493         IF( ln_vvl_ztilde ) THEN
494            DO jk = 1, jpk
495               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
496            END DO
497         ENDIF
498
499         ! 3 - Thickness diffusion term
500         ! ----------------------------
501         zwu(:,:) = 0._wp
502         zwv(:,:) = 0._wp
503         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
504            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
505               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
506            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
507               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
508            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
509            zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
510         END_3D
511         DO_2D( 1, 1, 1, 1 )
512            un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
513            vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
514         END_2D
515         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
516            tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
517               &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
518               &                                            ) * r1_e1e2t(ji,jj)
519         END_3D
520         !                       ! d - thickness diffusion transport: boundary conditions
521         !                             (stored for tracer advction and continuity equation)
522#if defined key_mpi3
523         CALL lbc_lnk_nc_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
524#else
525         CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
526#endif
527
528         ! 4 - Time stepping of baroclinic scale factors
529         ! ---------------------------------------------
530#if defined key_mpi3
531         CALL lbc_lnk_nc_multi( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
532#else
533         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
534#endif
535         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
536
537         ! Maximum deformation control
538         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
539         ALLOCATE( ze3t(jpi,jpj,jpk), llmsk(jpi,jpj,jpk) )
540         ze3t(:,:,jpk) = 0._wp
541         DO jk = 1, jpkm1
542            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
543         END DO
544         !
545         llmsk(   1:Nis1,:,:) = .FALSE.   ! exclude halos from the checked region
546         llmsk(Nie1: jpi,:,:) = .FALSE.
547         llmsk(:,   1:Njs1,:) = .FALSE.
548         llmsk(:,Nje1: jpj,:) = .FALSE.
549         !
550         llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp                  ! define only the inner domain
551         z_tmax = MAXVAL( ze3t(:,:,:), mask = llmsk )   ;   CALL mpp_max( 'domvvl', z_tmax )   ! max over the global domain
552         z_tmin = MINVAL( ze3t(:,:,:), mask = llmsk )   ;   CALL mpp_min( 'domvvl', z_tmin )   ! min over the global domain
553         ! - ML - test: for the moment, stop simulation for too large e3_t variations
554         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
555            CALL mpp_maxloc( 'domvvl', ze3t, llmsk, z_tmax, ijk_max )
556            CALL mpp_minloc( 'domvvl', ze3t, llmsk, z_tmin, ijk_min )
557            IF (lwp) THEN
558               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
559               WRITE(numout, *) 'at i, j, k=', ijk_max
560               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
561               WRITE(numout, *) 'at i, j, k=', ijk_min           
562               CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high')
563            ENDIF
564         ENDIF
565         DEALLOCATE( ze3t, llmsk )
566         ! - ML - end test
567         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
568         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
569         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
570
571         !
572         ! "tilda" change in the after scale factor
573         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
574         DO jk = 1, jpkm1
575            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
576         END DO
577         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
578         ! ==================================================================================
579         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
580         ! - ML - baroclinicity error should be better treated in the future
581         !        i.e. locally and not spread over the water column.
582         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
583         zht(:,:) = 0.
584         DO jk = 1, jpkm1
585            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
586         END DO
587         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
588         DO jk = 1, jpkm1
589            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
590         END DO
591
592      ENDIF
593
594      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
595      !                                           ! ---baroclinic part--------- !
596         DO jk = 1, jpkm1
597            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
598         END DO
599      ENDIF
600
601      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
602         !
603         IF( lwp ) WRITE(numout, *) 'kt =', kt
604         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
605            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
606            CALL mpp_max( 'domvvl', z_tmax )                             ! max over the global domain
607            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
608         END IF
609         !
610         zht(:,:) = 0.0_wp
611         DO jk = 1, jpkm1
612            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
613         END DO
614         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
615         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
616         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
617         !
618         zht(:,:) = 0.0_wp
619         DO jk = 1, jpkm1
620            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
621         END DO
622         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
623         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
624         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
625         !
626         zht(:,:) = 0.0_wp
627         DO jk = 1, jpkm1
628            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
629         END DO
630         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
631         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
632         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
633         !
634         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) )
635         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
636         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
637         !
638         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) )
639         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
640         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
641         !
642         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) )
643         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
644         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
645      END IF
646
647      ! *********************************** !
648      ! After scale factors at u- v- points !
649      ! *********************************** !
650
651      CALL dom_vvl_interpol( ssh(:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
652      CALL dom_vvl_interpol( ssh(:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
653
654      ! *********************************** !
655      ! After depths at u- v points         !
656      ! *********************************** !
657
658      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
659      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
660      DO jk = 2, jpkm1
661         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
662         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
663      END DO
664      !                                        ! Inverse of the local depth
665!!gm BUG ?  don't understand the use of umask_i here .....
666      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
667      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
668      !
669      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
670      !
671   END SUBROUTINE dom_vvl_sf_nxt
672
673
674
675   SUBROUTINE dom_vvl_sf_nxt_st( kt, Kbb, Kmm, Kaa, kcall ) 
676      !!----------------------------------------------------------------------
677      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
678      !!                   
679      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
680      !!                 tranxt and dynspg routines
681      !!
682      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
683      !!               - z_tilde_case: after scale factor increment =
684      !!                                    high frequency part of horizontal divergence
685      !!                                  + retsoring towards the background grid
686      !!                                  + thickness difusion
687      !!                               Then repartition of ssh INCREMENT proportionnaly
688      !!                               to the "baroclinic" level thickness.
689      !!
690      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
691      !!               - tilde_e3t_a: after increment of vertical scale factor
692      !!                              in z_tilde case
693      !!               - e3(t/u/v)_a
694      !!
695      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
696      !!----------------------------------------------------------------------
697      INTEGER, INTENT( in )           ::   kt             ! time step
698      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
699      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
700      !
701      INTEGER                ::   ji, jj, jk            ! dummy loop indices
702      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
703      REAL(wp)               ::   z_tmin, z_tmax        ! local scalars
704      LOGICAL                ::   ll_do_bclinic         ! local logical
705      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
706      !!----------------------------------------------------------------------
707      !
708      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
709      !
710      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
711      !
712      IF( kt == nit000 ) THEN
713         IF(lwp) WRITE(numout,*)
714         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
715         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
716      ENDIF
717
718      ll_do_bclinic = .TRUE.
719      IF( PRESENT(kcall) ) THEN
720         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
721      ENDIF
722
723      ! ******************************* !
724      ! After acale factors at t-points !
725      ! ******************************* !
726      !
727      DO jk = 1, jpkm1
728         e3t(:,:,jk,Kaa) = e3t_0(:,:,jk) * ( 1._wp + r3t(:,:,Kaa) )
729         e3u(:,:,jk,Kaa) = e3u_0(:,:,jk) * ( 1._wp + r3u(:,:,Kaa) )
730         e3v(:,:,jk,Kaa) = e3v_0(:,:,jk) * ( 1._wp + r3v(:,:,Kaa) )
731      END DO
732      !
733      ! *********************************** !
734      ! After scale factors at u- v- points !
735      ! *********************************** !
736     
737      !!st CALL dom_vvl_interpol_st( r3u(:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
738      !!st CALL dom_vvl_interpol_st( r3v(:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
739
740      ! *********************************** !
741      ! After depths at u- v points         !
742      ! *********************************** !
743
744      !!st hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
745      !!st hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
746      !!st DO jk = 2, jpkm1
747      !!st    hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
748      !!st    hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
749      !!st   
750      !!st END DO
751      hu(:,:,Kaa) = (hu_0(:,:)*(1._wp+r3u(:,:,Kaa)))
752      hv(:,:,Kaa) = (hv_0(:,:)*(1._wp+r3v(:,:,Kaa)))
753      !                                        ! Inverse of the local depth
754!!gm BUG ?  don't understand the use of umask_i here .....
755      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
756      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
757      !
758      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
759      !
760   END SUBROUTINE dom_vvl_sf_nxt_st
761   
762
763
764   SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa )
765      !!----------------------------------------------------------------------
766      !!                ***  ROUTINE dom_vvl_sf_update  ***
767      !!                   
768      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
769      !!               compute all depths and related variables for next time step
770      !!               write outputs and restart file
771      !!
772      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
773      !!               - reconstruct scale factor at other grid points (interpolate)
774      !!               - recompute depths and water height fields
775      !!
776      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step
777      !!               - Recompute:
778      !!                    e3(u/v)_b       
779      !!                    e3w(:,:,:,Kmm)           
780      !!                    e3(u/v)w_b     
781      !!                    e3(u/v)w_n     
782      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
783      !!                    h(u/v) and h(u/v)r
784      !!
785      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
786      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
787      !!----------------------------------------------------------------------
788      INTEGER, INTENT( in ) ::   kt              ! time step
789      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
790      !
791      INTEGER  ::   ji, jj, jk   ! dummy loop indices
792      REAL(wp) ::   zcoef        ! local scalar
793      !!----------------------------------------------------------------------
794      !
795      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
796      !
797      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update')
798      !
799      IF( kt == nit000 )   THEN
800         IF(lwp) WRITE(numout,*)
801         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
802         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
803      ENDIF
804      !
805      ! Time filter and swap of scale factors
806      ! =====================================
807      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
808      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
809         IF( l_1st_euler ) THEN
810            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
811         ELSE
812            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
813            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
814         ENDIF
815         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
816      ENDIF
817
818      ! Compute all missing vertical scale factor and depths
819      ! ====================================================
820      ! Horizontal scale factor interpolations
821      ! --------------------------------------
822      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
823      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
824     
825      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3f(:,:,:), 'F'  )
826     
827      ! Vertical scale factor interpolations
828      CALL dom_vvl_interpol( ssh(:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
829      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
830      CALL dom_vvl_interpol( ssh(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
831      CALL dom_vvl_interpol( ssh(:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
832      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
833      CALL dom_vvl_interpol( ssh(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
834
835      ! t- and w- points depth (set the isf depth as it is in the initial step)
836      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
837      gdepw(:,:,1,Kmm) = 0.0_wp
838      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
839      DO_3D( 1, 1, 1, 1, 2, jpk )
840        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
841                                                           ! 1 for jk = mikt
842         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
843         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
844         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
845             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
846         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
847      END_3D
848
849      ! Local depth and Inverse of the local depth of the water
850      ! -------------------------------------------------------
851      !
852      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
853      DO jk = 2, jpkm1
854         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
855      END DO
856
857      ! write restart file
858      ! ==================
859      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
860      !
861      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update')
862      !
863   END SUBROUTINE dom_vvl_sf_update
864   
865
866   SUBROUTINE dom_vvl_sf_update_st( kt, Kbb, Kmm, Kaa )
867      !!----------------------------------------------------------------------
868      !!                ***  ROUTINE dom_vvl_sf_update  ***
869      !!                   
870      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
871      !!               compute all depths and related variables for next time step
872      !!               write outputs and restart file
873      !!
874      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
875      !!               - reconstruct scale factor at other grid points (interpolate)
876      !!               - recompute depths and water height fields
877      !!
878      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step
879      !!               - Recompute:
880      !!                    e3(u/v)_b       
881      !!                    e3w(:,:,:,Kmm)           
882      !!                    e3(u/v)w_b     
883      !!                    e3(u/v)w_n     
884      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
885      !!                    h(u/v) and h(u/v)r
886      !!
887      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
888      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
889      !!----------------------------------------------------------------------
890      INTEGER, INTENT( in ) ::   kt              ! time step
891      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
892      !
893      INTEGER  ::   ji, jj, jk   ! dummy loop indices
894      REAL(wp) ::   zcoef        ! local scalar
895      !!----------------------------------------------------------------------
896      !
897      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
898      !
899      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update')
900      !
901      IF( kt == nit000 )   THEN
902         IF(lwp) WRITE(numout,*)
903         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
904         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
905      ENDIF
906      !
907
908      ! Compute all missing vertical scale factor and depths
909      ! ====================================================
910      ! Horizontal scale factor interpolations
911      ! --------------------------------------
912      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
913      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
914     
915      CALL dom_vvl_interpol_st( r3f(:,:), e3f(:,:,:), 'F'  )
916
917      ! Vertical scale factor interpolations
918      CALL dom_vvl_interpol_st( r3t(:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
919      CALL dom_vvl_interpol_st( r3u(:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
920      CALL dom_vvl_interpol_st( r3v(:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
921      CALL dom_vvl_interpol_st( r3t(:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
922      CALL dom_vvl_interpol_st( r3u(:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
923      CALL dom_vvl_interpol_st( r3v(:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
924
925      ! t- and w- points depth (set the isf depth as it is in the initial step)
926      DO_3D( 1, 1, 1, 1, 1, jpk )
927         gdepw(ji,jj,jk,Kmm) =  gdepw_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm))
928         gdept(ji,jj,jk,Kmm) =  gdept_0(ji,jj,jk) * (1._wp + r3t(ji,jj,Kmm)) 
929         gde3w(ji,jj,jk    ) =  gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
930      END_3D
931
932      ! Local depth and Inverse of the local depth of the water
933      ! -------------------------------------------------------
934      !
935      ht(:,:) = ht_0(:,:) + ssh(:,:,Kmm)
936
937      ! write restart file
938      ! ==================
939      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
940      !
941      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update')
942      !
943   END SUBROUTINE dom_vvl_sf_update_st
944   
945
946   
947   SUBROUTINE dom_vvl_interpol_st( rc3, pe3, cdp )
948      !!---------------------------------------------------------------------
949      !!                  ***  ROUTINE dom_vvl__interpol  ***
950      !!
951      !! ** Purpose :   interpolate scale factors from one grid point to another
952      !!
953      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
954      !!                - horizontal interpolation: grid cell surface averaging
955      !!                - vertical interpolation: simple averaging
956      !!----------------------------------------------------------------------
957      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   rc3     ! input e3   NOT used here (ssh is used instead)
958      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3     ! scale factor e3 to be updated   [m]
959      CHARACTER(LEN=*)          , INTENT(in   ) ::   cdp     ! grid point of the scale factor ( 'U', 'V', 'W, 'F', 'UW' or 'VW' )
960      !
961      INTEGER ::   ji, jj, jk                 ! dummy loop indices
962      REAL(wp), DIMENSION(jpi,jpj) ::   zc3   ! 2D workspace
963      !!----------------------------------------------------------------------
964      !
965      SELECT CASE ( cdp )     !==  type of interpolation  ==!
966         !
967      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
968         DO jk = 1, jpkm1
969            pe3(:,:,jk) = e3u_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
970         END DO
971         !
972      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
973         DO jk = 1, jpkm1
974            pe3(:,:,jk) = e3v_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
975         END DO
976         !
977      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
978         DO jk = 1, jpkm1                    ! Horizontal interpolation of e3f from ssh
979            e3f(:,:,jk) = e3f_0(:,:,jk) * ( 1._wp + rc3(:,:) )
980         END DO
981         !
982      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
983         DO jk = 1, jpk
984            pe3(:,:,jk) = e3w_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
985         END DO
986         !
987      CASE( 'UW' )                  !* from U- to UW-point
988         DO jk = 1, jpk
989            pe3(:,:,jk) = e3uw_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
990         END DO
991      CASE( 'VW' )                  !* from U- to UW-point : vertical simple mean
992         DO jk = 1, jpk
993            pe3(:,:,jk) = e3vw_0(:,:,jk) * ( 1.0_wp + rc3(:,:) )
994         END DO
995         !
996      END SELECT
997      !
998   END SUBROUTINE dom_vvl_interpol_st
999   
1000
1001   SUBROUTINE dom_vvl_interpol( pssh, pe3, cdp )
1002      !!---------------------------------------------------------------------
1003      !!                  ***  ROUTINE dom_vvl__interpol  ***
1004      !!
1005      !! ** Purpose :   interpolate scale factors from one grid point to another
1006      !!
1007      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
1008      !!                - horizontal interpolation: grid cell surface averaging
1009      !!                - vertical interpolation: simple averaging
1010      !!----------------------------------------------------------------------
1011      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pssh    ! input e3   NOT used here (ssh is used instead)
1012      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe3     ! scale factor e3 to be updated   [m]
1013      CHARACTER(LEN=*)          , INTENT(in   ) ::   cdp     ! grid point of the scale factor ( 'U', 'V', 'W, 'F', 'UW' or 'VW' )
1014      !
1015      INTEGER ::   ji, jj, jk                 ! dummy loop indices
1016      REAL(wp), DIMENSION(jpi,jpj) ::   zc3   ! 2D workspace
1017      !!----------------------------------------------------------------------
1018      !
1019      SELECT CASE ( cdp )     !==  type of interpolation  ==!
1020         !
1021      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
1022         DO_2D( 0, 0, 0, 0 )
1023            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
1024               &                   + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
1025         END_2D
1026#if defined key_mpi3
1027         CALL lbc_lnk_nc_multi( 'domvvl', zc3(:,:), 'U', 1._wp )
1028#else
1029         CALL lbc_lnk( 'domvvl', zc3(:,:), 'U', 1._wp )
1030#endif
1031         !
1032         DO jk = 1, jpkm1
1033            pe3(:,:,jk) = e3u_0(:,:,jk) * ( 1.0_wp + zc3(:,:) )
1034         END DO
1035         !
1036      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
1037         DO_2D( 0, 0, 0, 0 )
1038            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
1039               &                   + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
1040         END_2D
1041#if defined key_mpi3
1042         CALL lbc_lnk_nc_multi( 'domvvl', zc3(:,:), 'V', 1._wp )
1043#else
1044         CALL lbc_lnk( 'domvvl', zc3(:,:), 'V', 1._wp )
1045#endif
1046         !
1047         DO jk = 1, jpkm1
1048            pe3(:,:,jk) = e3v_0(:,:,jk) * ( 1.0_wp + zc3(:,:) )
1049         END DO
1050         !
1051      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
1052         DO_2D( 1, 0, 1, 0 )
1053            zc3(ji,jj) = 0.25_wp * (  e1e2t(ji  ,jj  ) * pssh(ji  ,jj  )  &
1054               &                    + e1e2t(ji+1,jj  ) * pssh(ji+1,jj  )  &
1055               &                    + e1e2t(ji  ,jj+1) * pssh(ji  ,jj+1)  &
1056               &                    + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1)  ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
1057         END_2D
1058#if defined key_mpi3
1059         CALL lbc_lnk_nc_multi( 'domvvl', zc3(:,:), 'F', 1._wp )
1060#else
1061         CALL lbc_lnk( 'domvvl', zc3(:,:), 'F', 1._wp )
1062#endif
1063         !
1064         DO jk = 1, jpkm1                    ! Horizontal interpolation of e3f from ssh
1065            e3f(:,:,jk)     = e3f_0(:,:,jk) * ( 1._wp + zc3(:,:) )
1066         END DO
1067         !
1068      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
1069         zc3(:,:) = pssh(:,:) * r1_ht_0(:,:)
1070         !
1071         DO jk = 1, jpk
1072            pe3(:,:,jk) = e3w_0(:,:,jk) * ( 1.0_wp + zc3(:,:) )
1073         END DO
1074         !
1075      CASE( 'UW' )                  !* from U- to UW-point
1076         !
1077         DO_2D( 0, 0, 0, 0 )
1078            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji  ,jj) * pssh(ji  ,jj)  &
1079               &                   + e1e2t(ji+1,jj) * pssh(ji+1,jj)  ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
1080         END_2D
1081#if defined key_mpi3
1082         CALL lbc_lnk_nc_multi( 'domvvl', zc3(:,:), 'U', 1._wp )
1083#else
1084         CALL lbc_lnk( 'domvvl', zc3(:,:), 'U', 1._wp )
1085#endif
1086         !
1087         DO jk = 1, jpk
1088            pe3(:,:,jk) = e3uw_0(:,:,jk) * ( 1.0_wp + zc3(:,:) )
1089         END DO
1090      CASE( 'VW' )                  !* from U- to UW-point : vertical simple mean
1091         !
1092         DO_2D( 0, 0, 0, 0 )
1093            zc3(ji,jj) = 0.5_wp * (  e1e2t(ji,jj  ) * pssh(ji,jj  )  &
1094               &                   + e1e2t(ji,jj+1) * pssh(ji,jj+1)  ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
1095         END_2D
1096#if defined key_mpi3
1097         CALL lbc_lnk_nc_multi( 'domvvl', zc3(:,:), 'V', 1._wp )
1098#else
1099         CALL lbc_lnk( 'domvvl', zc3(:,:), 'V', 1._wp )
1100#endif
1101          !
1102         DO jk = 1, jpk
1103            pe3(:,:,jk) = e3vw_0(:,:,jk) * ( 1.0_wp + zc3(:,:) )
1104         END DO
1105         !
1106      END SELECT
1107      !
1108   END SUBROUTINE dom_vvl_interpol
1109   
1110
1111   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
1112      !!---------------------------------------------------------------------
1113      !!                   ***  ROUTINE dom_vvl_rst  ***
1114      !!                     
1115      !! ** Purpose :   Read or write VVL file in restart file
1116      !!
1117      !! ** Method  :   use of IOM library
1118      !!                if the restart does not contain vertical scale factors,
1119      !!                they are set to the _0 values
1120      !!                if the restart does not contain vertical scale factors increments (z_tilde),
1121      !!                they are set to 0.
1122      !!----------------------------------------------------------------------
1123      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
1124      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
1125      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
1126      !
1127      INTEGER ::   ji, jj, jk
1128      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
1129      !!----------------------------------------------------------------------
1130      !
1131      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
1132         !                                   ! ===============
1133         IF( ln_rstart ) THEN                   !* Read the restart file
1134            CALL rst_read_open                  !  open the restart file if necessary
1135            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
1136            !
1137            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
1138            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
1139            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
1140            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
1141            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
1142            !
1143            !                             ! --------- !
1144            !                             ! all cases !
1145            !                             ! --------- !
1146            !
1147            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
1148               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
1149               CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
1150               ! needed to restart if land processor not computed
1151               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
1152               WHERE ( tmask(:,:,:) == 0.0_wp ) 
1153                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
1154                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
1155               END WHERE
1156               IF( l_1st_euler ) THEN
1157                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
1158               ENDIF
1159            ELSE IF( id1 > 0 ) THEN
1160               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
1161               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
1162               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
1163               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
1164               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
1165               l_1st_euler = .true.
1166            ELSE IF( id2 > 0 ) THEN
1167               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
1168               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
1169               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
1170               CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
1171               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
1172               l_1st_euler = .true.
1173            ELSE
1174               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
1175               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
1176               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
1177               DO jk = 1, jpk
1178                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
1179                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
1180                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
1181               END DO
1182               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
1183               l_1st_euler = .true.
1184            ENDIF
1185            !                             ! ----------- !
1186            IF( ln_vvl_zstar ) THEN       ! z_star case !
1187               !                          ! ----------- !
1188               IF( MIN( id3, id4 ) > 0 ) THEN
1189                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
1190               ENDIF
1191               !                          ! ----------------------- !
1192            ELSE                          ! z_tilde and layer cases !
1193               !                          ! ----------------------- !
1194               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
1195                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
1196                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
1197               ELSE                            ! one at least array is missing
1198                  tilde_e3t_b(:,:,:) = 0.0_wp
1199                  tilde_e3t_n(:,:,:) = 0.0_wp
1200               ENDIF
1201               !                          ! ------------ !
1202               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
1203                  !                       ! ------------ !
1204                  IF( id5 > 0 ) THEN  ! required array exists
1205                     CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
1206                  ELSE                ! array is missing
1207                     hdiv_lf(:,:,:) = 0.0_wp
1208                  ENDIF
1209               ENDIF
1210            ENDIF
1211            !
1212         ELSE                                   !* Initialize at "rest"
1213            !
1214
1215            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
1216               !
1217               IF( cn_cfg == 'wad' ) THEN
1218                  ! Wetting and drying test case
1219                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
1220!!an                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
1221                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
1222                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
1223                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
1224               ELSE
1225                  ! if not test case
1226                  ssh(:,:,Kmm) = -ssh_ref
1227                  ssh(:,:,Kbb) = -ssh_ref
1228
1229                  DO_2D( 1, 1, 1, 1 )
1230                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
1231                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
1232                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
1233                     ENDIF
1234                  END_2D
1235               ENDIF !If test case else
1236
1237               ! Adjust vertical metrics for all wad
1238               DO jk = 1, jpk
1239                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
1240                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
1241                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
1242               END DO
1243               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
1244
1245               DO ji = 1, jpi
1246                  DO jj = 1, jpj
1247                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
1248                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
1249                     ENDIF
1250                  END DO
1251               END DO 
1252               !
1253            ELSE
1254               !
1255               ! Just to read set ssh in fact, called latter once vertical grid
1256               ! is set up:
1257!               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
1258!               !
1259!               DO jk=1,jpk
1260!                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
1261!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
1262!               END DO
1263!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
1264                ssh(:,:,Kmm)=0._wp
1265                ssh(:,:,Kbb)=0._wp
1266                e3t(:,:,:,Kmm)=e3t_0(:,:,:)
1267                e3t(:,:,:,Kbb)=e3t_0(:,:,:)
1268               !
1269            END IF           ! end of ll_wd edits
1270
1271            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
1272               tilde_e3t_b(:,:,:) = 0._wp
1273               tilde_e3t_n(:,:,:) = 0._wp
1274               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
1275            END IF
1276         ENDIF
1277         !
1278      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1279         !                                   ! ===================
1280         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
1281         IF( lwxios ) CALL iom_swap(      cwxios_context          )
1282         !                                           ! --------- !
1283         !                                           ! all cases !
1284         !                                           ! --------- !
1285         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
1286         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
1287         !                                           ! ----------------------- !
1288         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
1289            !                                        ! ----------------------- !
1290            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
1291            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
1292         END IF
1293         !                                           ! -------------!   
1294         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
1295            !                                        ! ------------ !
1296            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
1297         ENDIF
1298         !
1299         IF( lwxios ) CALL iom_swap(      cxios_context          )
1300      ENDIF
1301      !
1302   END SUBROUTINE dom_vvl_rst
1303
1304
1305   SUBROUTINE dom_vvl_ctl
1306      !!---------------------------------------------------------------------
1307      !!                  ***  ROUTINE dom_vvl_ctl  ***
1308      !!               
1309      !! ** Purpose :   Control the consistency between namelist options
1310      !!                for vertical coordinate
1311      !!----------------------------------------------------------------------
1312      INTEGER ::   ioptio, ios
1313      !!
1314      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
1315         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
1316         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
1317      !!----------------------------------------------------------------------
1318      !
1319      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
1320901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
1321      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1322902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
1323      IF(lwm) WRITE ( numond, nam_vvl )
1324      !
1325      IF(lwp) THEN                    ! Namelist print
1326         WRITE(numout,*)
1327         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1328         WRITE(numout,*) '~~~~~~~~~~~'
1329         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
1330         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1331         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1332         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
1333         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1334         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1335         WRITE(numout,*) '      !'
1336         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
1337         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
1338         IF( ln_vvl_ztilde_as_zstar ) THEN
1339            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1340            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
1341            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
1342            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
1343            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1344            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rn_Dt'
1345         ELSE
1346            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
1347            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
1348         ENDIF
1349         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
1350      ENDIF
1351      !
1352      ioptio = 0                      ! Parameter control
1353      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1354      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1355      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1356      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1357      !
1358      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1359      !
1360      IF(lwp) THEN                   ! Print the choice
1361         WRITE(numout,*)
1362         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1363         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1364         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1365         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1366      ENDIF
1367      !
1368#if defined key_agrif
1369      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1370#endif
1371      !
1372   END SUBROUTINE dom_vvl_ctl
1373
1374#endif
1375!!stoops
1376
1377   !!======================================================================
1378END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.