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

source: NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DOM/domvvl.F90 @ 9923

Last change on this file since 9923 was 9923, checked in by gm, 6 years ago

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