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_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/OCE/DOM – NEMO

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/OCE/DOM/domvvl.F90 @ 13334

Last change on this file since 13334 was 13334, checked in by jchanut, 4 years ago

finish bypassing ocean/ice initialization with AGRIF, #2222, #2129

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