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/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/tests/CANAL/MY_SRC/domvvl.F90 @ 12424

Last change on this file since 12424 was 12424, checked in by davestorkey, 4 years ago
  1. Rename r2dt -> rDt
  2. Rename r1_2dt -> r1_Dt
  3. Reorganise management of initial Euler timestep for leapfrogging.

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