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 branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2015/dev_r5721_CNRS9_NOC3_LDF/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 5737

Last change on this file since 5737 was 5737, checked in by gm, 9 years ago

#1593: LDF-ADV, step I: Phasing of horizontal scale factors correct 2

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