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

source: branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 6069

Last change on this file since 6069 was 6069, checked in by timgraham, 8 years ago

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