source: trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 @ 7753

Last change on this file since 7753 was 7753, checked in by mocavero, 4 years ago

Reverting trunk to remove OpenMP

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