source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domvvl.F90 @ 10978

Last change on this file since 10978 was 10978, checked in by davestorkey, 18 months ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : DOM and sshwzv.F90. (sshwzv.F90 will be rewritten somewhat when we rewrite the time-level swapping but I've done it anyway). Passes all non-AGRIF SETTE tests.

  • 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_swp   : 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_swp     ! 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_n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
184      hu_b(:,:) = e3u(:,:,1,Kbb) * umask(:,:,1)
185      hu_n(:,:) = e3u(:,:,1,Kmm) * umask(:,:,1)
186      hv_b(:,:) = e3v(:,:,1,Kbb) * vmask(:,:,1)
187      hv_n(:,:) = e3v(:,:,1,Kmm) * vmask(:,:,1)
188      DO jk = 2, jpkm1
189         ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
190         hu_b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
191         hu_n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
192         hv_b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
193         hv_n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
194      END DO
195      !
196      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
197      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
198      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )
199      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
200      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 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_a(:,:) = e3u(:,:,1,Kaa) * umask(:,:,1)
553      hv_a(:,:) = e3v(:,:,1,Kaa) * vmask(:,:,1)
554      DO jk = 2, jpkm1
555         hu_a(:,:) = hu_a(:,:) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
556         hv_a(:,:) = hv_a(:,:) + 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_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
561      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 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_swp( kt, Kbb, Kmm, Kaa )
569      !!----------------------------------------------------------------------
570      !!                ***  ROUTINE dom_vvl_sf_swp  ***
571      !!                   
572      !! ** Purpose :  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
577      !!               - reconstruct scale factor at other grid points (interpolate)
578      !!               - recompute depths and water height fields
579      !!
580      !! ** Action  :  - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_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_swp')
602      !
603      IF( kt == nit000 )   THEN
604         IF(lwp) WRITE(numout,*)
605         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
606         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
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      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
622      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
623
624      e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)
625      e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa)
626      e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa)
627
628      ! Compute all missing vertical scale factor and depths
629      ! ====================================================
630      ! Horizontal scale factor interpolations
631      ! --------------------------------------
632      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are allready computed in dynnxt
633      ! - JC - hu_b, hv_b, hur_b, hvr_b also
634     
635      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
636     
637      ! Vertical scale factor interpolations
638      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
639      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
640      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
641      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
642      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
643      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
644
645      ! t- and w- points depth (set the isf depth as it is in the initial step)
646      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
647      gdepw(:,:,1,Kmm) = 0.0_wp
648      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
649      DO jk = 2, jpk
650         DO jj = 1,jpj
651            DO ji = 1,jpi
652              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
653                                                                 ! 1 for jk = mikt
654               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
655               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
656               gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
657                   &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
658               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
659            END DO
660         END DO
661      END DO
662
663      ! Local depth and Inverse of the local depth of the water
664      ! -------------------------------------------------------
665      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
666      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
667      !
668      ht_n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
669      DO jk = 2, jpkm1
670         ht_n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
671      END DO
672
673      ! write restart file
674      ! ==================
675      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
676      !
677      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
678      !
679   END SUBROUTINE dom_vvl_sf_swp
680
681
682   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
683      !!---------------------------------------------------------------------
684      !!                  ***  ROUTINE dom_vvl__interpol  ***
685      !!
686      !! ** Purpose :   interpolate scale factors from one grid point to another
687      !!
688      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
689      !!                - horizontal interpolation: grid cell surface averaging
690      !!                - vertical interpolation: simple averaging
691      !!----------------------------------------------------------------------
692      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
693      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
694      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
695      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
696      !
697      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
698      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
699      !!----------------------------------------------------------------------
700      !
701      IF(ln_wd_il) THEN
702        zlnwd = 1.0_wp
703      ELSE
704        zlnwd = 0.0_wp
705      END IF
706      !
707      SELECT CASE ( pout )    !==  type of interpolation  ==!
708         !
709      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
710         DO jk = 1, jpk
711            DO jj = 1, jpjm1
712               DO ji = 1, fs_jpim1   ! vector opt.
713                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
714                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
715                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
716               END DO
717            END DO
718         END DO
719         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
720         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
721         !
722      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
723         DO jk = 1, jpk
724            DO jj = 1, jpjm1
725               DO ji = 1, fs_jpim1   ! vector opt.
726                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
727                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
728                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
729               END DO
730            END DO
731         END DO
732         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
733         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
734         !
735      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
736         DO jk = 1, jpk
737            DO jj = 1, jpjm1
738               DO ji = 1, fs_jpim1   ! vector opt.
739                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
740                     &                       *    r1_e1e2f(ji,jj)                                                  &
741                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
742                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
743               END DO
744            END DO
745         END DO
746         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
747         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
748         !
749      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
750         !
751         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
752         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
753!!gm BUG? use here wmask in case of ISF ?  to be checked
754         DO jk = 2, jpk
755            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
756               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
757               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
758               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
759         END DO
760         !
761      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
762         !
763         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
764         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
765!!gm BUG? use here wumask in case of ISF ?  to be checked
766         DO jk = 2, jpk
767            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
768               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
769               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
770               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
771         END DO
772         !
773      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
774         !
775         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
776         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
777!!gm BUG? use here wvmask in case of ISF ?  to be checked
778         DO jk = 2, jpk
779            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
780               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
781               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
782               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
783         END DO
784      END SELECT
785      !
786   END SUBROUTINE dom_vvl_interpol
787
788
789   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
790      !!---------------------------------------------------------------------
791      !!                   ***  ROUTINE dom_vvl_rst  ***
792      !!                     
793      !! ** Purpose :   Read or write VVL file in restart file
794      !!
795      !! ** Method  :   use of IOM library
796      !!                if the restart does not contain vertical scale factors,
797      !!                they are set to the _0 values
798      !!                if the restart does not contain vertical scale factors increments (z_tilde),
799      !!                they are set to 0.
800      !!----------------------------------------------------------------------
801      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
802      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
803      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
804      !
805      INTEGER ::   ji, jj, jk
806      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
807      !!----------------------------------------------------------------------
808      !
809      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
810         !                                   ! ===============
811         IF( ln_rstart ) THEN                   !* Read the restart file
812            CALL rst_read_open                  !  open the restart file if necessary
813            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
814            !
815            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
816            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
817            id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. )
818            id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. )
819            id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. )
820            !                             ! --------- !
821            !                             ! all cases !
822            !                             ! --------- !
823            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
824               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
825               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
826               ! needed to restart if land processor not computed
827               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
828               WHERE ( tmask(:,:,:) == 0.0_wp ) 
829                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
830                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
831               END WHERE
832               IF( neuler == 0 ) THEN
833                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
834               ENDIF
835            ELSE IF( id1 > 0 ) THEN
836               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
837               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
838               IF(lwp) write(numout,*) 'neuler is forced to 0'
839               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
840               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
841               neuler = 0
842            ELSE IF( id2 > 0 ) THEN
843               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
844               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
845               IF(lwp) write(numout,*) 'neuler is forced to 0'
846               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
847               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
848               neuler = 0
849            ELSE
850               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
851               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
852               IF(lwp) write(numout,*) 'neuler is forced to 0'
853               DO jk = 1, jpk
854                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
855                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
856                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
857               END DO
858               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
859               neuler = 0
860            ENDIF
861            !                             ! ----------- !
862            IF( ln_vvl_zstar ) THEN       ! z_star case !
863               !                          ! ----------- !
864               IF( MIN( id3, id4 ) > 0 ) THEN
865                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
866               ENDIF
867               !                          ! ----------------------- !
868            ELSE                          ! z_tilde and layer cases !
869               !                          ! ----------------------- !
870               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
871                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lrxios )
872                  CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lrxios )
873               ELSE                            ! one at least array is missing
874                  tilde_e3t_b(:,:,:) = 0.0_wp
875                  tilde_e3t_n(:,:,:) = 0.0_wp
876               ENDIF
877               !                          ! ------------ !
878               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
879                  !                       ! ------------ !
880                  IF( id5 > 0 ) THEN  ! required array exists
881                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lrxios )
882                  ELSE                ! array is missing
883                     hdiv_lf(:,:,:) = 0.0_wp
884                  ENDIF
885               ENDIF
886            ENDIF
887            !
888         ELSE                                   !* Initialize at "rest"
889            !
890
891            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
892               !
893               IF( cn_cfg == 'wad' ) THEN
894                  ! Wetting and drying test case
895                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
896                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
897                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
898                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
899                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
900               ELSE
901                  ! if not test case
902                  ssh(:,:,Kmm) = -ssh_ref
903                  ssh(:,:,Kbb) = -ssh_ref
904
905                  DO jj = 1, jpj
906                     DO ji = 1, jpi
907                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
908                           ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
909                           ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
910                        ENDIF
911                     ENDDO
912                  ENDDO
913               ENDIF !If test case else
914
915               ! Adjust vertical metrics for all wad
916               DO jk = 1, jpk
917                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
918                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
919                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
920               END DO
921               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
922
923               DO ji = 1, jpi
924                  DO jj = 1, jpj
925                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
926                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
927                     ENDIF
928                  END DO
929               END DO 
930               !
931            ELSE
932               !
933               ! Just to read set ssh in fact, called latter once vertical grid
934               ! is set up:
935!               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
936!               !
937!               DO jk=1,jpk
938!                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
939!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
940!               END DO
941!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
942                ssh(:,:,Kmm)=0._wp
943                e3t(:,:,:,Kmm)=e3t_0(:,:,:)
944                e3t(:,:,:,Kbb)=e3t_0(:,:,:)
945               !
946            END IF           ! end of ll_wd edits
947
948            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
949               tilde_e3t_b(:,:,:) = 0._wp
950               tilde_e3t_n(:,:,:) = 0._wp
951               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
952            END IF
953         ENDIF
954         !
955      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
956         !                                   ! ===================
957         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
958         IF( lwxios ) CALL iom_swap(      cwxios_context          )
959         !                                           ! --------- !
960         !                                           ! all cases !
961         !                                           ! --------- !
962         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
963         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
964         !                                           ! ----------------------- !
965         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
966            !                                        ! ----------------------- !
967            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
968            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
969         END IF
970         !                                           ! -------------!   
971         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
972            !                                        ! ------------ !
973            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
974         ENDIF
975         !
976         IF( lwxios ) CALL iom_swap(      cxios_context          )
977      ENDIF
978      !
979   END SUBROUTINE dom_vvl_rst
980
981
982   SUBROUTINE dom_vvl_ctl
983      !!---------------------------------------------------------------------
984      !!                  ***  ROUTINE dom_vvl_ctl  ***
985      !!               
986      !! ** Purpose :   Control the consistency between namelist options
987      !!                for vertical coordinate
988      !!----------------------------------------------------------------------
989      INTEGER ::   ioptio, ios
990      !!
991      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
992         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
993         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
994      !!----------------------------------------------------------------------
995      !
996      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
997      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
998901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
999      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
1000      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
1001902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
1002      IF(lwm) WRITE ( numond, nam_vvl )
1003      !
1004      IF(lwp) THEN                    ! Namelist print
1005         WRITE(numout,*)
1006         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1007         WRITE(numout,*) '~~~~~~~~~~~'
1008         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
1009         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1010         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1011         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
1012         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1013         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1014         WRITE(numout,*) '      !'
1015         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
1016         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
1017         IF( ln_vvl_ztilde_as_zstar ) THEN
1018            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1019            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
1020            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
1021            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
1022            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1023            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
1024         ELSE
1025            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
1026            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
1027         ENDIF
1028         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
1029      ENDIF
1030      !
1031      ioptio = 0                      ! Parameter control
1032      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1033      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1034      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1035      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1036      !
1037      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1038      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1039      !
1040      IF(lwp) THEN                   ! Print the choice
1041         WRITE(numout,*)
1042         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1043         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1044         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1045         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1046      ENDIF
1047      !
1048#if defined key_agrif
1049      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1050#endif
1051      !
1052   END SUBROUTINE dom_vvl_ctl
1053
1054   !!======================================================================
1055END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.