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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/VORTEX/MY_SRC – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/VORTEX/MY_SRC/domvvl.F90 @ 11480

Last change on this file since 11480 was 11480, checked in by davestorkey, 5 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Merge in changes from branch of branch.
Main changes:

  1. "nxt" modules renamed as "atf" and now just do Asselin time filtering. The time level swapping is achieved by swapping indices.
  2. Some additional prognostic grid variables changed to use a time dimension.

Notes:

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