New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domvvl.F90 in NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/tests/VORTEX/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/tests/VORTEX/MY_SRC/domvvl.F90 @ 13689

Last change on this file since 13689 was 13689, checked in by jchanut, 3 years ago

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