source: NEMO/trunk/src/OCE/DOM/domvvl.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 14 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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