New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
domvvl.F90 in NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM – NEMO

source: NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/domvvl.F90 @ 9976

Last change on this file since 9976 was 9976, checked in by acc, 6 years ago

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