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/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM – NEMO

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 3865

Last change on this file since 3865 was 3865, checked in by acc, 11 years ago

Branch 2013/dev_r3858_NOC_ZTC, #863. Nearly complete port of 2011/dev_r2739_LOCEAN8_ZTC development branch into v3.5aplha base. Compiles and runs but currently unstable after 8 timesteps with ORCA2_LIM reference configuration.
  • Property svn:keywords set to Id
File size: 39.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   !!----------------------------------------------------------------------
11   !!   'key_vvl'                              variable volume
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   !! * Modules used
22   USE oce             ! ocean dynamics and tracers
23   USE dom_oce         ! ocean space and time domain
24   USE sbc_oce         ! ocean surface boundary condition
25   USE in_out_manager  ! I/O manager
26   USE iom             ! I/O manager library
27   USE restart         ! ocean restart
28   USE lib_mpp         ! distributed memory computing library
29   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
30   USE wrk_nemo        ! Memory allocation
31   USE timing          ! Timing
32
33   IMPLICIT NONE
34   PRIVATE
35
36   !! * Routine accessibility
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_kepe   = .FALSE.   ! kinetic/potential energy transfer
47   !                                                                                  ! conservation: not used yet
48   REAL(wp)                                              :: rn_ahe3       =  0.e0     ! thickness diffusion coefficient
49   LOGICAL , PUBLIC                                      :: ln_vvl_dbg    = .FALSE.   ! debug control prints
50
51   !! * Module variables
52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td              ! thickness diffusion transport
53   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                   ! low frequency part of hz divergence
54   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3t_t_b, e3t_t_n, e3t_t_a ! baroclinic scale factors
55   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t               ! retoring period for scale factors
56   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv               ! retoring period for low freq. divergence
57
58   REAL(wp),         ALLOCATABLE, SAVE, DIMENSION(:)     ::   r2dt   ! vertical profile time-step, = 2 rdttra
59      !                                                              ! except at nit000 (=rdttra) if neuler=0
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( e3t_t_b(jpi,jpj,jpk) , e3t_t_n(jpi,jpj,jpk) , e3t_t_a(jpi,jpj,jpk) ,   &
79            &      un_td  (jpi,jpj,jpk) , vn_td  (jpi,jpj,jpk) , STAT = dom_vvl_alloc        )
80         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
81         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
82      ENDIF
83      IF( ln_vvl_ztilde ) THEN
84         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
85         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
86         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
87      ENDIF
88
89   END FUNCTION dom_vvl_alloc
90
91
92   SUBROUTINE dom_vvl_init
93      !!----------------------------------------------------------------------
94      !!                ***  ROUTINE dom_vvl_init  ***
95      !!                   
96      !! ** Purpose :  Initialization of all scale factors, depths
97      !!               and water column heights
98      !!
99      !! ** Method  :  - use restart file and/or initialize
100      !!               - interpolate scale factors
101      !!
102      !! ** Action  : - fse3t_(n/b) and e3t_t_(n/b)
103      !!              - Regrid: fse3(u/v)_n
104      !!                        fse3(u/v)_b       
105      !!                        fse3w_n           
106      !!                        fse3(u/v)w_b     
107      !!                        fse3(u/v)w_n     
108      !!                        fsdept_n, fsdepw_n and fsde3w_n
109      !!              - h(t/u/v)_0
110      !!              - frq_rst_e3t and frq_rst_hdv
111      !!
112      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
113      !!----------------------------------------------------------------------
114      USE phycst,  ONLY : rpi
115      !! * Local declarations
116      INTEGER ::   jk
117      !!----------------------------------------------------------------------
118      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_init')
119
120      IF(lwp) WRITE(numout,*)
121      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
122      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
123
124      ! choose vertical coordinate (z_star, z_tilde or layer)
125      ! ==========================
126      CALL dom_vvl_ctl
127
128      ! Allocate module arrays
129      ! ======================
130      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
131
132      ! Read or initialize fse3t_(b/n), e3t_t_(b/n) and hdiv_lf (and e3t_a(jpk))
133      ! ========================================================================
134      CALL dom_vvl_rst( nit000, 'READ' )
135      fse3t_a(:,:,jpk) = e3t_0(:,:,jpk)
136
137      ! Reconstruction of all vertical scale factors at now and before time steps
138      ! =========================================================================
139      ! Horizontal scale factor interpolations
140      ! --------------------------------------
141      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' )
142      CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' )
143      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' )
144      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' )
145      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' )
146      ! Vertical scale factor interpolations
147      ! ------------------------------------
148      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
149      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
150      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
151      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
152      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
153      ! t- and w- points depth
154      ! ----------------------
155      fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1)
156      fsdepw_n(:,:,1) = 0.e0
157      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
158      DO jk = 2, jpk
159         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
160         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
161         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
162      END DO
163      ! Reference water column height at t-, u- and v- point
164      ! ----------------------------------------------------
165      ht_0(:,:) = 0.e0
166      hu_0(:,:) = 0.e0
167      hv_0(:,:) = 0.e0
168      DO jk = 1, jpk
169         ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk)
170         hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk)
171         hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk)
172      END DO
173
174      ! Restoring frequencies for z_tilde coordinate
175      ! ============================================
176      IF( ln_vvl_ztilde ) THEN
177         ! - ML - In the future, this should be tunable by the user (namelist)
178         frq_rst_e3t(:,:) = 2.e0_wp * rpi / ( 30.e0_wp * 86400.e0_wp )
179         frq_rst_hdv(:,:) = 2.e0_wp * rpi / (  5.e0_wp * 86400.e0_wp )
180      ENDIF
181
182      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_init')
183
184   END SUBROUTINE dom_vvl_init
185
186
187   SUBROUTINE dom_vvl_sf_nxt( kt ) 
188      !!----------------------------------------------------------------------
189      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
190      !!                   
191      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
192      !!                 tranxt and dynspg routines
193      !!
194      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
195      !!               - z_tilde_case: after scale factor increment =
196      !!                                    high frequency part of horizontal divergence
197      !!                                  + retsoring towards the background grid
198      !!                                  + thickness difusion
199      !!                               Then repartition of ssh INCREMENT proportionnaly
200      !!                               to the "baroclinic" level thickness.
201      !!
202      !! ** Action  :  - hdiv_lf: restoring towards full baroclinic divergence in z_tilde case
203      !!               - e3t_t_a: after increment of vertical scale factor
204      !!                          in z_tilde case
205      !!               - fse3(t/u/v)_a
206      !!
207      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
208      !!----------------------------------------------------------------------
209      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t
210      REAL(wp), POINTER, DIMENSION(:,:  ) :: zht, z_scale, zwu, zwv, zhdiv
211      !! * Arguments
212      INTEGER, INTENT( in )                  :: kt                    ! time step
213      !! * Local declarations
214      INTEGER                                :: ji, jj, jk            ! dummy loop indices
215      INTEGER , DIMENSION(3)                 :: ijk_max, ijk_min      ! temporary integers
216      REAL(wp)                               :: z2dt                  ! temporary scalars
217      REAL(wp)                               :: z_def_max             ! temporary scalar
218      REAL(wp)                               :: z_tmin, z_tmax        ! temporary scalars
219      !!----------------------------------------------------------------------
220      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
221      CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
222      CALL wrk_alloc( jpi, jpj, jpk, ze3t                     )
223
224      IF(kt == nit000)   THEN
225         IF(lwp) WRITE(numout,*)
226         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
227         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
228      ENDIF
229
230      ! ******************************* !
231      ! After acale factors at t-points !
232      ! ******************************* !
233
234      !                                                ! ----------------- !
235      IF( ln_vvl_zstar ) THEN                          ! z_star coordinate !
236         !                                             ! ----------------- !
237
238         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
239         DO jk = 1, jpkm1
240            fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
241         END DO
242
243      !                                                ! --------------------------- !
244      ELSEIF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN   ! z_tilde or layer coordinate !
245         !                                             ! --------------------------- !
246
247         ! I - initialization
248         ! ==================
249
250         ! 1 - barotropic divergence
251         ! -------------------------
252         zhdiv(:,:) = 0.
253         zht(:,:)   = 0.
254         DO jk = 1, jpkm1
255            zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)
256            zht  (:,:) = zht  (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
257         END DO
258         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) )
259
260         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
261         ! --------------------------------------------------
262         IF( ln_vvl_ztilde ) THEN
263            IF( kt .GT. nit000 ) THEN
264               DO jk = 1, jpkm1
265                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
266                     &          * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) )
267               END DO
268            ENDIF
269         END IF
270
271         ! II - after z_tilde increments of vertical scale factors
272         ! =======================================================
273         e3t_t_a(:,:,:) = 0.e0   ! e3t_t_a used to store tendency terms
274
275         ! 1 - High frequency divergence term
276         ! ----------------------------------
277         IF( ln_vvl_ztilde ) THEN   ! z_tilde case
278            DO jk = 1, jpkm1
279               e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
280            END DO
281         ELSE                       ! layer case
282            DO jk = 1, jpkm1
283               e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) )
284            END DO
285         END IF
286
287         ! 2 - Restoring term (z-tilde case only)
288         ! ------------------
289         IF( ln_vvl_ztilde ) THEN
290            DO jk = 1, jpk
291               e3t_t_a(:,:,jk) = e3t_t_a(:,:,jk) - frq_rst_e3t(:,:) * e3t_t_b(:,:,jk)
292            END DO
293         END IF
294
295         ! 3 - Thickness diffusion term
296         ! ----------------------------
297         zwu(:,:) = 0.e0
298         zwv(:,:) = 0.e0
299         ! a - first derivative: diffusive fluxes
300         DO jk = 1, jpkm1
301            DO jj = 1, jpjm1
302               DO ji = 1, fs_jpim1   ! vector opt.
303                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji+1,jj  ,jk) )
304                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) * ( e3t_t_b(ji,jj,jk) - e3t_t_b(ji  ,jj+1,jk) )
305                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
306                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
307               END DO
308            END DO
309         END DO
310         ! b - correction for last oceanic u-v points
311         DO jj = 1, jpj
312            DO ji = 1, jpi
313               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
314               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
315            END DO
316         END DO
317         ! c - second derivative: divergence of diffusive fluxes
318         DO jk = 1, jpkm1
319            DO jj = 2, jpjm1
320               DO ji = fs_2, fs_jpim1   ! vector opt.
321                  e3t_t_a(ji,jj,jk) = e3t_t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
322                     &                                      + vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk) )  &
323                     &                                    * r1_e12t(ji,jj)
324               END DO
325            END DO
326         END DO
327         ! d - thickness diffusion transport: boundary conditions
328         !     (stored for tracer advction and continuity equation)
329         CALL lbc_lnk( un_td , 'U' , -1.)
330         CALL lbc_lnk( vn_td , 'V' , -1.)
331
332         ! 4 - Time stepping of baroclinic scale factors
333         ! ---------------------------------------------
334         ! Leapfrog time stepping
335         ! ~~~~~~~~~~~~~~~~~~~~~~
336         IF( neuler == 0 .AND. kt == nit000 ) THEN
337            z2dt =  rdt
338         ELSE
339            z2dt = 2.e0 * rdt
340         ENDIF
341         CALL lbc_lnk( e3t_t_a(:,:,:), 'T', 1. )
342         e3t_t_a(:,:,:) = e3t_t_b(:,:,:) + z2dt * tmask(:,:,:) * e3t_t_a(:,:,:)
343
344         ! Maximum deformation control
345         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
346         ! - ML - Should perhaps be put in the namelist
347         z_def_max = 0.9e0
348         ze3t(:,:,jpk) = 0.e0
349         DO jk = 1, jpkm1
350            ze3t(:,:,jk) = e3t_t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
351         END DO
352         z_tmax = MAXVAL( ze3t(:,:,:) )
353         z_tmin = MINVAL( ze3t(:,:,:) )
354         ! - ML - test: for the moment, stop simulation for too large e3_t variations
355         IF( ( z_tmax .GT. z_def_max ) .OR. ( z_tmin .LT. - z_def_max ) ) THEN
356            ijk_max = MAXLOC( ze3t(:,:,:) )
357            ijk_min = MINLOC( ze3t(:,:,:) )
358            WRITE(numout, *) 'MAX( e3t_t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
359            WRITE(numout, *) 'at i, j, k=', ijk_max
360            WRITE(numout, *) 'MIN( e3t_t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
361            WRITE(numout, *) 'at i, j, k=', ijk_min           
362            CALL ctl_stop('MAX( ABS( e3t_t_a(:,:,:) ) / e3t_0(:,:,:) ) too high')
363         ENDIF
364         ! - ML - end test
365         ! - ML - This will cause a baroclinicity error if the ctl_stop above is not used
366         e3t_t_a(:,:,:) = MIN( e3t_t_a(:,:,:),   z_def_max * e3t_0(:,:,:) )
367         e3t_t_a(:,:,:) = MAX( e3t_t_a(:,:,:), - z_def_max * e3t_0(:,:,:) )
368
369         ! Add "tilda" part to the after scale factor
370         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
371         fse3t_a(:,:,:) = e3t_0(:,:,:) + e3t_t_a(:,:,:)
372
373         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
374         ! ==================================================================================
375         ! add e3t(n-1) "star" Asselin-filtered
376         DO jk = 1, jpkm1
377            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - e3t_0(:,:,jk) - e3t_t_b(:,:,jk)
378         END DO
379         ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n)
380         ! - ML - baroclinicity error should be better treated in the future
381         !        i.e. locally and not spread over the water column.
382         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
383         zht(:,:) = 0.
384         DO jk = 1, jpkm1
385            zht(:,:) = zht(:,:) + e3t_t_a(:,:,jk) * tmask(:,:,jk)
386         END DO
387         z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) )
388         DO jk = 1, jpkm1
389            fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)
390         END DO
391
392      ENDIF
393
394      IF( ln_vvl_dbg ) THEN   ! - ML - test: control prints for debuging
395         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
396            WRITE(numout, *) 'kt =', kt
397            WRITE(numout, *) 'MAXVAL(abs(SUM(e3t_t_a))) =',   &
398               &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
399         END IF
400         zht(:,:) = 0.e0
401         DO jk = 1, jpkm1
402            zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk)
403         END DO
404         WRITE(numout, *) 'MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =',   &
405            &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) )
406         zht(:,:) = 0.e0
407         DO jk = 1, jpkm1
408            zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk)
409         END DO
410         WRITE(numout, *) 'MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =',   &
411            &              MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) )
412      END IF
413
414      ! *********************************** !
415      ! After scale factors at u- v- points !
416      ! *********************************** !
417
418      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' )
419      CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' )
420
421      CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv )
422      CALL wrk_dealloc( jpi, jpj, jpk, ze3t                     )
423
424      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_nxt')
425
426   END SUBROUTINE dom_vvl_sf_nxt
427
428
429   SUBROUTINE dom_vvl_sf_swp( kt )
430      !!----------------------------------------------------------------------
431      !!                ***  ROUTINE dom_vvl_sf_swp  ***
432      !!                   
433      !! ** Purpose :  compute time filter and swap of scale factors
434      !!               compute all depths and related variables for next time step
435      !!               write outputs and restart file
436      !!
437      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
438      !!               - reconstruct scale factor at other grid points (interpolate)
439      !!               - recompute depths and water height fields
440      !!
441      !! ** Action  :  - fse3t_(b/n), e3t_t_(b/n) and fse3(u/v)_n ready for next time step
442      !!               - Recompute:
443      !!                    fse3(u/v)_b       
444      !!                    fse3w_n           
445      !!                    fse3(u/v)w_b     
446      !!                    fse3(u/v)w_n     
447      !!                    fsdept_n, fsdepw_n  and fsde3w_n
448      !!                    h(u/v) and h(u/v)r
449      !!
450      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
451      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
452      !!----------------------------------------------------------------------
453      !! * Arguments
454      INTEGER, INTENT( in )               :: kt       ! time step
455      !! * Local declarations
456      REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def
457      INTEGER                             :: jk       ! dummy loop indices
458      !!----------------------------------------------------------------------
459
460      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_sf_swp')
461      !
462      CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def                )
463      !
464      IF( kt == nit000 )   THEN
465         IF(lwp) WRITE(numout,*)
466         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
467         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
468      ENDIF
469      !
470      ! Time filter and swap of scale factors
471      ! =====================================
472      ! - ML - fse3(t/u/v)_b are allready computed in dynnxt.
473      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
474         IF( neuler == 0 .AND. kt == nit000 ) THEN
475            e3t_t_b(:,:,:) = e3t_t_n(:,:,:)
476         ELSE
477            e3t_t_b(:,:,:) = e3t_t_n(:,:,:) + atfp * ( e3t_t_b(:,:,:) - 2.e0 * e3t_t_n(:,:,:) + e3t_t_a(:,:,:) )
478         ENDIF
479         e3t_t_n(:,:,:) = e3t_t_a(:,:,:)
480      ENDIF
481      fse3t_n(:,:,:) = fse3t_a(:,:,:)
482      fse3u_n(:,:,:) = fse3u_a(:,:,:)
483      fse3v_n(:,:,:) = fse3v_a(:,:,:)
484
485      ! Compute all missing vertical scale factor and depths
486      ! ====================================================
487      ! Horizontal scale factor interpolations
488      ! --------------------------------------
489      ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt
490      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F'  )
491      ! Vertical scale factor interpolations
492      ! ------------------------------------
493      CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W'  )
494      CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' )
495      CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' )
496      CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' )
497      CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' )
498      ! t- and w- points depth
499      ! ----------------------
500      fsdept_n(:,:,1) = 0.5 * fse3w_n(:,:,1)
501      fsdepw_n(:,:,1) = 0.e0
502      fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:)
503      DO jk = 2, jpk
504         fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk)
505         fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1)
506         fsde3w_n(:,:,jk) = fsdept_n(:,:,jk  ) - sshn   (:,:)
507      END DO
508      ! Local depth and Inverse of the local depth of the water column at u- and v- points
509      ! ----------------------------------------------------------------------------------
510      hu(:,:) = 0.
511      hv(:,:) = 0.
512      DO jk = 1, jpk
513         hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk)
514         hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk)
515      END DO
516      ! Inverse of the local depth
517      hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1. - umask(:,:,1) )
518      hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1. - vmask(:,:,1) )
519
520      ! Write outputs
521      ! =============
522      z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
523      CALL iom_put( "e3t_n"  , fse3t_n  (:,:,:) )
524      CALL iom_put( "dept_n" , fsde3w_n (:,:,:) )
525      CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) )
526
527      ! write restart file
528      ! ==================
529      IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' )
530      !
531      CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def )
532      !
533      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_sf_swp')
534
535   END SUBROUTINE dom_vvl_sf_swp
536
537
538   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
539      !!---------------------------------------------------------------------
540      !!                  ***  ROUTINE dom_vvl__interpol  ***
541      !!
542      !! ** Purpose :   interpolate scale factors from one grid point to another
543      !!
544      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
545      !!                - horizontal interpolation: grid cell surface averaging
546      !!                - vertical interpolation: simple averaging
547      !!----------------------------------------------------------------------
548      !! * Arguments
549      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
550      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
551      CHARACTER(LEN=1), INTENT( in )                    ::  pout       ! grid point of out scale factors
552      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
553      !! * Local declarations
554      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
555      !!----------------------------------------------------------------------
556      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
557      SELECT CASE ( pout )
558         !               ! ------------------------------------- !
559      CASE( 'U' )        ! interpolation from T-point to U-point !
560         !               ! ------------------------------------- !
561         ! horizontal surface weighted interpolation
562         DO jk = 1, jpk
563            DO jj = 1, jpjm1
564               DO ji = 1, fs_jpim1   ! vector opt.
565                  pe3_out(ji,jj,jk) = 0.5 * umask(ji,jj,jk) * r1_e12u(ji,jj)                                      &
566                     &                    * (   e12t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
567                     &                        + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
568               END DO
569            END DO
570         END DO
571         ! boundary conditions
572         CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. )
573         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
574         !               ! ------------------------------------- !
575      CASE( 'V' )        ! interpolation from T-point to V-point !
576         !               ! ------------------------------------- !
577         ! horizontal surface weighted interpolation
578         DO jk = 1, jpk
579            DO jj = 1, jpjm1
580               DO ji = 1, fs_jpim1   ! vector opt.
581                  pe3_out(ji,jj,jk) = 0.5 * vmask(ji,jj,jk) * r1_e12v(ji,jj)                                      &
582                     &                    * (   e12t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
583                     &                        + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
584               END DO
585            END DO
586         END DO
587         ! boundary conditions
588         CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. )
589         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
590         !               ! ------------------------------------- !
591      CASE( 'F' )        ! interpolation from U-point to F-point !
592         !               ! ------------------------------------- !
593         ! horizontal surface weighted interpolation
594         DO jk = 1, jpk
595            DO jj = 1, jpjm1
596               DO ji = 1, fs_jpim1   ! vector opt.
597                  pe3_out(ji,jj,jk) = 0.5 * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj)                  &
598                     &                    * (   e12u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
599                     &                        + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
600               END DO
601            END DO
602         END DO
603         ! boundary conditions
604         CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. )
605         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
606         !               ! ------------------------------------- !
607      CASE( 'W' )        ! interpolation from T-point to W-point !
608         !               ! ------------------------------------- !
609         ! vertical simple interpolation
610         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
611         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
612         DO jk = 2, jpk
613            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1. - 0.5 * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )   &
614               &                              +        0.5 * tmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
615         END DO
616         !               ! -------------------------------------- !
617      CASE( 'UW' )       ! interpolation from U-point to UW-point !
618         !               ! -------------------------------------- !
619         ! vertical simple interpolation
620         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
621         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
622         DO jk = 2, jpk
623            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1. - 0.5 * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )   &
624               &                               +        0.5 * umask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
625         END DO
626         !               ! -------------------------------------- !
627      CASE( 'VW' )       ! interpolation from V-point to VW-point !
628         !               ! -------------------------------------- !
629         ! vertical simple interpolation
630         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
631         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
632         DO jk = 2, jpk
633            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1. - 0.5 * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )   &
634               &                               +        0.5 * vmask(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
635         END DO
636      END SELECT
637
638      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
639
640   END SUBROUTINE dom_vvl_interpol
641
642
643   SUBROUTINE dom_vvl_rst( kt, cdrw )
644      !!---------------------------------------------------------------------
645      !!                   ***  ROUTINE dom_vvl_rst  ***
646      !!                     
647      !! ** Purpose :   Read or write VVL file in restart file
648      !!
649      !! ** Method  :   use of IOM library
650      !!                if the restart does not contain vertical scale factors,
651      !!                they are set to the _0 values
652      !!                if the restart does not contain vertical scale factors increments (z_tilde),
653      !!                they are set to 0.
654      !!----------------------------------------------------------------------
655      !! * Arguments
656      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
657      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
658      !! * Local declarations
659      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
660      !!----------------------------------------------------------------------
661      !
662      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_rst')
663      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
664         !                                   ! ===============
665         IF( ln_rstart ) THEN                   !* Read the restart file
666            id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. )
667            id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. )
668            id3 = iom_varid( numror, 'e3t_t_b', ldstop = .FALSE. )
669            id4 = iom_varid( numror, 'e3t_t_n', ldstop = .FALSE. )
670            id5 = iom_varid( numror, 'hdif_lf', ldstop = .FALSE. )
671            !                             ! --------- !
672            !                             ! all cases !
673            !                             ! --------- !
674            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
675               CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) )
676               CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) )
677               IF( neuler == 0 ) THEN
678                  fse3t_b(:,:,:) = fse3t_n(:,:,:)
679               ENDIF
680            ELSE                                 ! one at least array is missing
681               CALL ctl_stop( 'dom_vvl_rst: vvl cannot restart from a non vvl run' )
682            ENDIF
683            !                             ! ----------- !
684            IF( ln_vvl_zstar ) THEN       ! z_star case !
685               !                          ! ----------- !
686               IF( MIN( id3, id4 ) > 0 ) THEN
687                  CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' )
688               ENDIF
689               !                          ! ----------------------- !
690            ELSE                          ! z_tilde and layer cases !
691               !                          ! ----------------------- !
692               IF( MIN( id3, id4 ) > 0 ) THEN  ! all required arrays exist
693                  CALL iom_get( numror, jpdom_autoglo, 'e3t_t_b', e3t_t_b(:,:,:) )
694                  CALL iom_get( numror, jpdom_autoglo, 'e3t_t_n', e3t_t_n(:,:,:) )
695               ELSE                            ! one at least array is missing
696                  e3t_t_b(:,:,:) = 0.e0
697                  e3t_t_n(:,:,:) = 0.e0
698               ENDIF
699               !                          ! ------------ !
700               IF( ln_vvl_ztilde ) THEN   ! z_tilde case !
701                  !                       ! ------------ !
702                  IF( id5 > 0 ) THEN  ! required array exists
703                     CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) )
704                  ELSE                ! array is missing
705                     hdiv_lf(:,:,:) = 0.e0
706                  ENDIF
707               ENDIF
708            ENDIF
709            !
710         ELSE                                   !* Initialize at "rest"
711            fse3t_b(:,:,:) = e3t_0(:,:,:)
712            fse3t_n(:,:,:) = e3t_0(:,:,:)
713            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
714               e3t_t_b(:,:,:) = 0.e0
715               e3t_t_n(:,:,:) = 0.e0
716               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.e0
717            END IF
718         ENDIF
719
720      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
721         !                                   ! ===================
722         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
723         !                                           ! --------- !
724         !                                           ! all cases !
725         !                                           ! --------- !
726         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) )
727         CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) )
728         !                                           ! ----------------------- !
729         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
730            !                                        ! ----------------------- !
731            CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_b', e3t_t_b(:,:,:) )
732            CALL iom_rstput( kt, nitrst, numrow, 'e3t_t_n', e3t_t_n(:,:,:) )
733         END IF
734         !                                           ! -------------!   
735         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
736            !                                        ! ------------ !
737            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) )
738         ENDIF
739
740      ENDIF
741      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_rst')
742
743   END SUBROUTINE dom_vvl_rst
744
745
746   SUBROUTINE dom_vvl_ctl
747      !!---------------------------------------------------------------------
748      !!                  ***  ROUTINE dom_vvl_ctl  ***
749      !!               
750      !! ** Purpose :   Control the consistency between namelist options
751      !!                for vertical coordinate
752      !!----------------------------------------------------------------------
753      INTEGER ::   ioptio
754
755      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, rn_ahe3, ln_vvl_dbg! , ln_vvl_kepe
756      !!----------------------------------------------------------------------
757
758      REWIND ( numnam )               ! Read Namelist nam_vvl : vertical coordinate
759      READ   ( numnam, nam_vvl )
760
761      IF(lwp) THEN                    ! Namelist print
762         WRITE(numout,*)
763         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
764         WRITE(numout,*) '~~~~~~~~~~~'
765         WRITE(numout,*) '           Namelist nam_vvl : chose a vertical coordinate'
766         WRITE(numout,*) '              zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
767         WRITE(numout,*) '              ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
768         WRITE(numout,*) '              layer                      ln_vvl_layer   = ', ln_vvl_layer
769         ! WRITE(numout,*) '           Namelist nam_vvl : chose kinetic-to-potential energy conservation'
770         ! WRITE(numout,*) '                                         ln_vvl_kepe    = ', ln_vvl_kepe
771         WRITE(numout,*) '           Namelist nam_vvl : thickness diffusion coefficient'
772         WRITE(numout,*) '                                         rn_ahe3        = ', rn_ahe3
773         WRITE(numout,*) '           Namelist nam_vvl : debug prints'
774         WRITE(numout,*) '                                         ln_vvl_dbg     = ', ln_vvl_dbg
775      ENDIF
776
777      ioptio = 0                      ! Parameter control
778      IF( ln_vvl_zstar     )   ioptio = ioptio + 1
779      IF( ln_vvl_ztilde    )   ioptio = ioptio + 1
780      IF( ln_vvl_layer     )   ioptio = ioptio + 1
781
782      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
783
784      IF(lwp) THEN                   ! Print the choice
785         WRITE(numout,*)
786         IF( ln_vvl_zstar      ) WRITE(numout,*) '              zstar vertical coordinate is used'
787         IF( ln_vvl_ztilde     ) WRITE(numout,*) '              ztilde vertical coordinate is used'
788         IF( ln_vvl_layer      ) WRITE(numout,*) '              layer vertical coordinate is used'
789         ! - ML - Option not developed yet
790         ! IF(       ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option used'
791         ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) '              kinetic to potential energy transfer : option not used'
792      ENDIF
793
794   END SUBROUTINE dom_vvl_ctl
795
796   !!======================================================================
797END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.