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/VORTEX/MY_SRC – NEMO

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