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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/tests/CANAL/MY_SRC/domvvl.F90 @ 11960

Last change on this file since 11960 was 11960, checked in by acc, 4 years ago

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