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

source: NEMO/branches/2020/dev_r13787_doc_latex_recovery/src/OCE/DOM/domvvl.F90 @ 14098

Last change on this file since 14098 was 14098, checked in by nicolasmartin, 3 years ago

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