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/2019/UKMO_MERGE_2019/src/OCE/DOM – NEMO

source: NEMO/branches/2019/UKMO_MERGE_2019/src/OCE/DOM/domvvl.F90 @ 12079

Last change on this file since 12079 was 12068, checked in by davestorkey, 4 years ago

2019/UKMO_MERGE_2019 : Merging in changes from ENHANCE-02_ISF_nemo.

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