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_r13898_Tiling_Cleanup_MPI3/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2020/dev_r13898_Tiling_Cleanup_MPI3/tests/CANAL/MY_SRC/domvvl.F90 @ 13906

Last change on this file since 13906 was 13906, checked in by mocavero, 3 years ago

Merge with dev_r13296_HPC-07_mocavero_mpi3
  • Property svn:keywords set to Id
File size: 56.9 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(:,:,:), ALLOCATABLE ::   ze3t
337      LOGICAL , DIMENSION(:,:,:), ALLOCATABLE ::   llmsk
338      !!----------------------------------------------------------------------
339      !
340      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
341      !
342      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
343      !
344      IF( kt == nit000 ) THEN
345         IF(lwp) WRITE(numout,*)
346         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
347         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
348      ENDIF
349
350      ll_do_bclinic = .TRUE.
351      IF( PRESENT(kcall) ) THEN
352         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
353      ENDIF
354
355      ! ******************************* !
356      ! After acale factors at t-points !
357      ! ******************************* !
358      !                                                ! --------------------------------------------- !
359      !                                                ! z_star coordinate and barotropic z-tilde part !
360      !                                                ! --------------------------------------------- !
361      !
362      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
363      DO jk = 1, jpkm1
364         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
365         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
366      END DO
367      !
368      IF( (ln_vvl_ztilde .OR. ln_vvl_layer) .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
369         !                                                               ! ------baroclinic part------ !
370         ! I - initialization
371         ! ==================
372
373         ! 1 - barotropic divergence
374         ! -------------------------
375         zhdiv(:,:) = 0._wp
376         zht(:,:)   = 0._wp
377         DO jk = 1, jpkm1
378            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
379            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
380         END DO
381         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
382
383         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
384         ! --------------------------------------------------
385         IF( ln_vvl_ztilde ) THEN
386            IF( kt > nit000 ) THEN
387               DO jk = 1, jpkm1
388                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rn_Dt * frq_rst_hdv(:,:)   &
389                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
390               END DO
391            ENDIF
392         ENDIF
393
394         ! II - after z_tilde increments of vertical scale factors
395         ! =======================================================
396         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
397
398         ! 1 - High frequency divergence term
399         ! ----------------------------------
400         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
401            DO jk = 1, jpkm1
402               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
403            END DO
404         ELSE                         ! layer case
405            DO jk = 1, jpkm1
406               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
407            END DO
408         ENDIF
409
410         ! 2 - Restoring term (z-tilde case only)
411         ! ------------------
412         IF( ln_vvl_ztilde ) THEN
413            DO jk = 1, jpk
414               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
415            END DO
416         ENDIF
417
418         ! 3 - Thickness diffusion term
419         ! ----------------------------
420         zwu(:,:) = 0._wp
421         zwv(:,:) = 0._wp
422         DO_3D( 1, 0, 1, 0, 1, jpkm1 )
423            un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
424               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
425            vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
426               &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
427            zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
428            zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
429         END_3D
430         DO_2D( 1, 1, 1, 1 )
431            un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
432            vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
433         END_2D
434         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
435            tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
436               &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
437               &                                            ) * r1_e1e2t(ji,jj)
438         END_3D
439         !                       ! d - thickness diffusion transport: boundary conditions
440         !                             (stored for tracer advction and continuity equation)
441#if defined key_mpi3
442         CALL lbc_lnk_nc_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
443#else
444         CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
445#endif
446
447         ! 4 - Time stepping of baroclinic scale factors
448         ! ---------------------------------------------
449#if defined key_mpi3
450         CALL lbc_lnk_nc_multi( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
451#else
452         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
453#endif
454         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + rDt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
455
456         ! Maximum deformation control
457         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
458         ALLOCATE( ze3t(jpi,jpj,jpk), llmsk(jpi,jpj,jpk) )
459         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
460            ze3t(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) / e3t_0(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
461         END_3D
462         !
463         llmsk(   1:Nis1,:,:) = .FALSE.   ! exclude halos from the checked region
464         llmsk(Nie1: jpi,:,:) = .FALSE.
465         llmsk(:,   1:Njs1,:) = .FALSE.
466         llmsk(:,Nje1: jpj,:) = .FALSE.
467         !
468         llmsk(Nis0:Nie0,Njs0:Nje0,:) = tmask(Nis0:Nie0,Njs0:Nje0,:) == 1._wp                  ! define only the inner domain
469         z_tmax = MAXVAL( ze3t(:,:,:), mask = llmsk )   ;   CALL mpp_max( 'domvvl', z_tmax )   ! max over the global domain
470         z_tmin = MINVAL( ze3t(:,:,:), mask = llmsk )   ;   CALL mpp_min( 'domvvl', z_tmin )   ! min over the global domain
471         ! - ML - test: for the moment, stop simulation for too large e3_t variations
472         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
473            CALL mpp_maxloc( 'domvvl', ze3t, llmsk, z_tmax, ijk_max )
474            CALL mpp_minloc( 'domvvl', ze3t, llmsk, z_tmin, ijk_min )
475            IF (lwp) THEN
476               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
477               WRITE(numout, *) 'at i, j, k=', ijk_max
478               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
479               WRITE(numout, *) 'at i, j, k=', ijk_min           
480               CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high')
481            ENDIF
482         ENDIF
483         DEALLOCATE( ze3t, llmsk )
484         ! - ML - end test
485         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
486         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
487         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
488
489         !
490         ! "tilda" change in the after scale factor
491         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
492         DO jk = 1, jpkm1
493            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
494         END DO
495         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
496         ! ==================================================================================
497         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
498         ! - ML - baroclinicity error should be better treated in the future
499         !        i.e. locally and not spread over the water column.
500         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
501         zht(:,:) = 0.
502         DO jk = 1, jpkm1
503            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
504         END DO
505         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
506         DO jk = 1, jpkm1
507            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
508         END DO
509
510      ENDIF
511
512      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
513      !                                           ! ---baroclinic part--------- !
514         DO jk = 1, jpkm1
515            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
516         END DO
517      ENDIF
518
519      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
520         !
521         IF( lwp ) WRITE(numout, *) 'kt =', kt
522         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
523            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
524            CALL mpp_max( 'domvvl', z_tmax )                             ! max over the global domain
525            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
526         END IF
527         !
528         zht(:,:) = 0.0_wp
529         DO jk = 1, jpkm1
530            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
531         END DO
532         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
533         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
534         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
535         !
536         zht(:,:) = 0.0_wp
537         DO jk = 1, jpkm1
538            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
539         END DO
540         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
541         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
542         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
543         !
544         zht(:,:) = 0.0_wp
545         DO jk = 1, jpkm1
546            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
547         END DO
548         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
549         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
550         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
551         !
552         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) )
553         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
554         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
555         !
556         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) )
557         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
558         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
559         !
560         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) )
561         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
562         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
563      END IF
564
565      ! *********************************** !
566      ! After scale factors at u- v- points !
567      ! *********************************** !
568
569      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
570      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
571
572      ! *********************************** !
573      ! After depths at u- v points         !
574      ! *********************************** !
575
576      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
577      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
578      DO jk = 2, jpkm1
579         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
580         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
581      END DO
582      !                                        ! Inverse of the local depth
583!!gm BUG ?  don't understand the use of umask_i here .....
584      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
585      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
586      !
587      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
588      !
589   END SUBROUTINE dom_vvl_sf_nxt
590
591
592   SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa )
593      !!----------------------------------------------------------------------
594      !!                ***  ROUTINE dom_vvl_sf_update  ***
595      !!                   
596      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
597      !!               compute all depths and related variables for next time step
598      !!               write outputs and restart file
599      !!
600      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE)
601      !!               - reconstruct scale factor at other grid points (interpolate)
602      !!               - recompute depths and water height fields
603      !!
604      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step
605      !!               - Recompute:
606      !!                    e3(u/v)_b       
607      !!                    e3w(:,:,:,Kmm)           
608      !!                    e3(u/v)w_b     
609      !!                    e3(u/v)w_n     
610      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
611      !!                    h(u/v) and h(u/v)r
612      !!
613      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
614      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
615      !!----------------------------------------------------------------------
616      INTEGER, INTENT( in ) ::   kt              ! time step
617      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
618      !
619      INTEGER  ::   ji, jj, jk   ! dummy loop indices
620      REAL(wp) ::   zcoef        ! local scalar
621      !!----------------------------------------------------------------------
622      !
623      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
624      !
625      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update')
626      !
627      IF( kt == nit000 )   THEN
628         IF(lwp) WRITE(numout,*)
629         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
630         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
631      ENDIF
632      !
633      ! Time filter and swap of scale factors
634      ! =====================================
635      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
636      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
637         IF( l_1st_euler ) THEN
638            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
639         ELSE
640            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
641            &         + rn_atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
642         ENDIF
643         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
644      ENDIF
645
646      ! Compute all missing vertical scale factor and depths
647      ! ====================================================
648      ! Horizontal scale factor interpolations
649      ! --------------------------------------
650      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
651      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
652     
653      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
654     
655      ! Vertical scale factor interpolations
656      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
657      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
658      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
659      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
660      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
661      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
662
663      ! t- and w- points depth (set the isf depth as it is in the initial step)
664      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
665      gdepw(:,:,1,Kmm) = 0.0_wp
666      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
667      DO_3D( 1, 1, 1, 1, 2, jpk )
668        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
669                                                           ! 1 for jk = mikt
670         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
671         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
672         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
673             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
674         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
675      END_3D
676
677      ! Local depth and Inverse of the local depth of the water
678      ! -------------------------------------------------------
679      !
680      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
681      DO jk = 2, jpkm1
682         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
683      END DO
684
685      ! write restart file
686      ! ==================
687      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
688      !
689      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update')
690      !
691   END SUBROUTINE dom_vvl_sf_update
692
693
694   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
695      !!---------------------------------------------------------------------
696      !!                  ***  ROUTINE dom_vvl__interpol  ***
697      !!
698      !! ** Purpose :   interpolate scale factors from one grid point to another
699      !!
700      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
701      !!                - horizontal interpolation: grid cell surface averaging
702      !!                - vertical interpolation: simple averaging
703      !!----------------------------------------------------------------------
704      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
705      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
706      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
707      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
708      !
709      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
710      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
711      !!----------------------------------------------------------------------
712      !
713      IF(ln_wd_il) THEN
714        zlnwd = 1.0_wp
715      ELSE
716        zlnwd = 0.0_wp
717      END IF
718      !
719      SELECT CASE ( pout )    !==  type of interpolation  ==!
720         !
721      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
722         DO_3D( 1, 0, 1, 0, 1, jpk )
723            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
724               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
725               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
726         END_3D
727#if defined key_mpi3
728         CALL lbc_lnk_nc_multi( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
729#else
730         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
731#endif
732         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
733         !
734      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
735         DO_3D( 1, 0, 1, 0, 1, jpk )
736            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
737               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
738               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
739         END_3D
740#if defined key_mpi3
741         CALL lbc_lnk_nc_multi( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
742#else
743         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
744#endif
745         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
746         !
747      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
748         DO_3D( 1, 0, 1, 0, 1, jpk )
749            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
750               &                       *    r1_e1e2f(ji,jj)                                                  &
751               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
752               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
753         END_3D
754#if defined key_mpi3
755         CALL lbc_lnk_nc_multi( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
756#else
757         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
758#endif
759         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
760         !
761      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
762         !
763         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
764         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
765!!gm BUG? use here wmask in case of ISF ?  to be checked
766         DO jk = 2, jpk
767            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
768               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
769               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
770               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
771         END DO
772         !
773      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
774         !
775         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
776         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
777!!gm BUG? use here wumask in case of ISF ?  to be checked
778         DO jk = 2, jpk
779            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
780               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
781               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
782               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
783         END DO
784         !
785      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
786         !
787         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
788         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
789!!gm BUG? use here wvmask in case of ISF ?  to be checked
790         DO jk = 2, jpk
791            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
792               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
793               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
794               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
795         END DO
796      END SELECT
797      !
798   END SUBROUTINE dom_vvl_interpol
799
800
801   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
802      !!---------------------------------------------------------------------
803      !!                   ***  ROUTINE dom_vvl_rst  ***
804      !!                     
805      !! ** Purpose :   Read or write VVL file in restart file
806      !!
807      !! ** Method  :   use of IOM library
808      !!                if the restart does not contain vertical scale factors,
809      !!                they are set to the _0 values
810      !!                if the restart does not contain vertical scale factors increments (z_tilde),
811      !!                they are set to 0.
812      !!----------------------------------------------------------------------
813      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
814      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
815      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
816      !
817      INTEGER ::   ji, jj, jk
818      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
819      !!----------------------------------------------------------------------
820      !
821      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
822         !                                   ! ===============
823         IF( ln_rstart ) THEN                   !* Read the restart file
824            CALL rst_read_open                  !  open the restart file if necessary
825            CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
826            !
827            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
828            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
829            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
830            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
831            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
832            !
833            !                             ! --------- !
834            !                             ! all cases !
835            !                             ! --------- !
836            !
837            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
838               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
839               CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
840               ! needed to restart if land processor not computed
841               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
842               WHERE ( tmask(:,:,:) == 0.0_wp ) 
843                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
844                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
845               END WHERE
846               IF( l_1st_euler ) THEN
847                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
848               ENDIF
849            ELSE IF( id1 > 0 ) THEN
850               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
851               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
852               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
853               CALL iom_get( numror, jpdom_auto, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
854               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
855               l_1st_euler = .true.
856            ELSE IF( id2 > 0 ) THEN
857               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
858               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
859               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
860               CALL iom_get( numror, jpdom_auto, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
861               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
862               l_1st_euler = .true.
863            ELSE
864               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
865               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
866               IF(lwp) write(numout,*) 'l_1st_euler is forced to true'
867               DO jk = 1, jpk
868                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
869                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
870                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
871               END DO
872               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
873               l_1st_euler = .true.
874            ENDIF
875            !                             ! ----------- !
876            IF( ln_vvl_zstar ) THEN       ! z_star case !
877               !                          ! ----------- !
878               IF( MIN( id3, id4 ) > 0 ) THEN
879                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
880               ENDIF
881               !                          ! ----------------------- !
882            ELSE                          ! z_tilde and layer cases !
883               !                          ! ----------------------- !
884               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
885                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
886                  CALL iom_get( numror, jpdom_auto, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
887               ELSE                            ! one at least array is missing
888                  tilde_e3t_b(:,:,:) = 0.0_wp
889                  tilde_e3t_n(:,:,:) = 0.0_wp
890               ENDIF
891               !                          ! ------------ !
892               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
893                  !                       ! ------------ !
894                  IF( id5 > 0 ) THEN  ! required array exists
895                     CALL iom_get( numror, jpdom_auto, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
896                  ELSE                ! array is missing
897                     hdiv_lf(:,:,:) = 0.0_wp
898                  ENDIF
899               ENDIF
900            ENDIF
901            !
902         ELSE                                   !* Initialize at "rest"
903            !
904
905            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
906               !
907               IF( cn_cfg == 'wad' ) THEN
908                  ! Wetting and drying test case
909                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
910                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
911                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
912                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
913                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
914               ELSE
915                  ! if not test case
916                  ssh(:,:,Kmm) = -ssh_ref
917                  ssh(:,:,Kbb) = -ssh_ref
918
919                  DO_2D( 1, 1, 1, 1 )
920                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
921                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
922                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
923                     ENDIF
924                  END_2D
925               ENDIF !If test case else
926
927               ! Adjust vertical metrics for all wad
928               DO jk = 1, jpk
929                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
930                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
931                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
932               END DO
933               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
934
935               DO_2D( 1, 1, 1, 1 )
936                  IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
937                     CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
938                  ENDIF
939               END_2D
940               !
941            ELSE
942               !
943               ! usr_def_istate called here only to get ssh(Kbb) needed to initialize e3t(Kbb) and e3t(Kmm)
944               !
945               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  ) 
946               !
947               ! usr_def_istate will be called again in istate_init to initialize ts, ssh, u and v
948               !
949               DO jk=1,jpk
950                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
951                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
952                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t(:,:,:,Kbb) != 0 on land points
953               END DO
954               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
955               ssh(:,:,Kmm) = ssh(:,:,Kbb)                                     ! needed later for gde3w
956               !
957            END IF           ! end of ll_wd edits
958
959            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
960               tilde_e3t_b(:,:,:) = 0._wp
961               tilde_e3t_n(:,:,:) = 0._wp
962               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
963            END IF
964         ENDIF
965         !
966      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
967         !                                   ! ===================
968         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
969         IF( lwxios ) CALL iom_swap(      cwxios_context          )
970         !                                           ! --------- !
971         !                                           ! all cases !
972         !                                           ! --------- !
973         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
974         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
975         !                                           ! ----------------------- !
976         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
977            !                                        ! ----------------------- !
978            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
979            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
980         END IF
981         !                                           ! -------------!   
982         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
983            !                                        ! ------------ !
984            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
985         ENDIF
986         !
987         IF( lwxios ) CALL iom_swap(      cxios_context          )
988      ENDIF
989      !
990   END SUBROUTINE dom_vvl_rst
991
992
993   SUBROUTINE dom_vvl_ctl
994      !!---------------------------------------------------------------------
995      !!                  ***  ROUTINE dom_vvl_ctl  ***
996      !!               
997      !! ** Purpose :   Control the consistency between namelist options
998      !!                for vertical coordinate
999      !!----------------------------------------------------------------------
1000      INTEGER ::   ioptio, ios
1001      !!
1002      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
1003         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
1004         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
1005      !!----------------------------------------------------------------------
1006      !
1007      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
1008901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
1009      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1010902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
1011      IF(lwm) WRITE ( numond, nam_vvl )
1012      !
1013      IF(lwp) THEN                    ! Namelist print
1014         WRITE(numout,*)
1015         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1016         WRITE(numout,*) '~~~~~~~~~~~'
1017         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
1018         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1019         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1020         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
1021         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1022         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1023         WRITE(numout,*) '      !'
1024         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
1025         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
1026         IF( ln_vvl_ztilde_as_zstar ) THEN
1027            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1028            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
1029            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
1030            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
1031            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1032            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rn_Dt'
1033         ELSE
1034            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
1035            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
1036         ENDIF
1037         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
1038      ENDIF
1039      !
1040      ioptio = 0                      ! Parameter control
1041      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1042      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1043      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1044      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1045      !
1046      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1047      !
1048      IF(lwp) THEN                   ! Print the choice
1049         WRITE(numout,*)
1050         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1051         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1052         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1053         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1054      ENDIF
1055      !
1056#if defined key_agrif
1057      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1058#endif
1059      !
1060   END SUBROUTINE dom_vvl_ctl
1061
1062#endif
1063
1064   !!======================================================================
1065END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.