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/2019/dev_r11943_MERGE_2019/tests/CANAL/MY_SRC – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/tests/CANAL/MY_SRC/domvvl.F90 @ 12353

Last change on this file since 12353 was 12353, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. Additions to the do loop macro implementation: converted a few loops previously missed because they used jpi-1 instead of jpim1 etc.; changed internal macro names in do_loop_substitute.h90 to strings that are much more unlikely to appear in any future code elsewhere and removed the key_vectopt_loop option (and all related code) since the do loop macros have suppressed this option. These changes have been fully SETTE-tested and this branch should now be ready to go back to the trunk.
  • Property svn:keywords set to Id
File size: 56.1 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   !!            4.1  !  2019-08  (A. Coward, D. Storkey) rename dom_vvl_sf_swp -> dom_vvl_sf_update for new timestepping
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_update   : 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 timing          ! Timing
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC  dom_vvl_init       ! called by domain.F90
39   PUBLIC  dom_vvl_sf_nxt     ! called by step.F90
40   PUBLIC  dom_vvl_sf_update  ! called by step.F90
41   PUBLIC  dom_vvl_interpol   ! called by dynnxt.F90
42
43   !                                                      !!* Namelist nam_vvl
44   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
45   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
46   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
47   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
48   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
49   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
50   !                                                       ! conservation: not used yet
51   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
52   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days]
53   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days]
54   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation
55   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
56
57   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
58   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf                     ! low frequency part of hz divergence
59   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n    ! baroclinic scale factors
60   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a   ! baroclinic scale factors
61   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_e3t                 ! retoring period for scale factors
62   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   :: frq_rst_hdv                 ! retoring period for low freq. divergence
63
64   !!----------------------------------------------------------------------
65   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
66   !! $Id$
67   !! Software governed by the CeCILL license (see ./LICENSE)
68   !!----------------------------------------------------------------------
69CONTAINS
70
71   INTEGER FUNCTION dom_vvl_alloc()
72      !!----------------------------------------------------------------------
73      !!                ***  FUNCTION dom_vvl_alloc  ***
74      !!----------------------------------------------------------------------
75      IF( ln_vvl_zstar )   dom_vvl_alloc = 0
76      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
77         ALLOCATE( tilde_e3t_b(jpi,jpj,jpk)  , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) ,   &
78            &      dtilde_e3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
79            &      STAT = dom_vvl_alloc        )
80         CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
81         IF( dom_vvl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
82         un_td = 0._wp
83         vn_td = 0._wp
84      ENDIF
85      IF( ln_vvl_ztilde ) THEN
86         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
87         CALL mpp_sum ( 'domvvl', dom_vvl_alloc )
88         IF( dom_vvl_alloc /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_alloc: failed to allocate arrays' )
89      ENDIF
90      !
91   END FUNCTION dom_vvl_alloc
92
93
94   SUBROUTINE dom_vvl_init( Kbb, Kmm, Kaa )
95      !!----------------------------------------------------------------------
96      !!                ***  ROUTINE dom_vvl_init  ***
97      !!                   
98      !! ** Purpose :  Initialization of all scale factors, depths
99      !!               and water column heights
100      !!
101      !! ** Method  :  - use restart file and/or initialize
102      !!               - interpolate scale factors
103      !!
104      !! ** Action  : - e3t_(n/b) and tilde_e3t_(n/b)
105      !!              - Regrid: e3[u/v](:,:,:,Kmm)
106      !!                        e3[u/v](:,:,:,Kmm)       
107      !!                        e3w(:,:,:,Kmm)           
108      !!                        e3[u/v]w_b
109      !!                        e3[u/v]w_n     
110      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
111      !!              - h(t/u/v)_0
112      !!              - frq_rst_e3t and frq_rst_hdv
113      !!
114      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
115      !!----------------------------------------------------------------------
116      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
117      !
118      INTEGER ::   ji, jj, jk
119      INTEGER ::   ii0, ii1, ij0, ij1
120      REAL(wp)::   zcoef
121      !!----------------------------------------------------------------------
122      !
123      IF(lwp) WRITE(numout,*)
124      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
125      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
126      !
127      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
128      !
129      !                    ! Allocate module arrays
130      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
131      !
132      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
133      CALL dom_vvl_rst( nit000, Kbb, Kmm, 'READ' )
134      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
135      !
136      !                    !== Set of all other vertical scale factors  ==!  (now and before)
137      !                                ! Horizontal interpolation of e3t
138      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U
139      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
140      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V
141      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
142      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F
143      !                                ! Vertical interpolation of e3t,u,v
144      CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W
145      CALL dom_vvl_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  )
146      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW
147      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
148      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW
149      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
150
151      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
152      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
153      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
154      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
155      !
156      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
157      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)       ! reference to the ocean surface (used for MLD and light penetration)
158      gdepw(:,:,1,Kmm) = 0.0_wp
159      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
160      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
161      gdepw(:,:,1,Kbb) = 0.0_wp
162      DO jk = 2, jpk                               ! vertical sum
163         DO jj = 1,jpj
164            DO ji = 1,jpi
165               !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
166               !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
167               !                             ! 0.5 where jk = mikt     
168!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ??
169               zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
170               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
171               gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
172                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm)) 
173               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
174               gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
175               gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  &
176                  &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb)) 
177            END DO
178         END DO
179      END DO
180      !
181      !                    !==  thickness of the water column  !!   (ocean portion only)
182      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
183      hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1)
184      hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1)
185      hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1)
186      hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1)
187      DO jk = 2, jpkm1
188         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
189         hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
190         hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
191         hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
192         hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
193      END DO
194      !
195      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
196      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
197      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
198      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
199      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
200
201      !                    !==   z_tilde coordinate case  ==!   (Restoring frequencies)
202      IF( ln_vvl_ztilde ) THEN
203!!gm : idea: add here a READ in a file of custumized restoring frequency
204         !                                   ! Values in days provided via the namelist
205         !                                   ! use rsmall to avoid possible division by zero errors with faulty settings
206         frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.0_wp )
207         frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp )
208         !
209         IF( ln_vvl_ztilde_as_zstar ) THEN   ! z-star emulation using z-tile
210            frq_rst_e3t(:,:) = 0._wp               !Ignore namelist settings
211            frq_rst_hdv(:,:) = 1._wp / rdt
212         ENDIF
213         IF ( ln_vvl_zstar_at_eqtor ) THEN   ! use z-star in vicinity of the Equator
214            DO jj = 1, jpj
215               DO ji = 1, jpi
216!!gm  case |gphi| >= 6 degrees is useless   initialized just above by default
217                  IF( ABS(gphit(ji,jj)) >= 6.) THEN
218                     ! values outside the equatorial band and transition zone (ztilde)
219                     frq_rst_e3t(ji,jj) =  2.0_wp * rpi / ( MAX( rn_rst_e3t  , rsmall ) * 86400.e0_wp )
220                     frq_rst_hdv(ji,jj) =  2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp )
221                  ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN    ! Equator strip ==> z-star
222                     ! values inside the equatorial band (ztilde as zstar)
223                     frq_rst_e3t(ji,jj) =  0.0_wp
224                     frq_rst_hdv(ji,jj) =  1.0_wp / rdt
225                  ELSE                                      ! transition band (2.5 to 6 degrees N/S)
226                     !                                      ! (linearly transition from z-tilde to z-star)
227                     frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp   &
228                        &            * (  1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
229                        &                                          * 180._wp / 3.5_wp ) )
230                     frq_rst_hdv(ji,jj) = (1.0_wp / rdt)                                &
231                        &            + (  frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp   &
232                        &            * (  1._wp  - COS( rad*(ABS(gphit(ji,jj))-2.5_wp)  &
233                        &                                          * 180._wp / 3.5_wp ) )
234                  ENDIF
235               END DO
236            END DO
237            IF( cn_cfg == "orca" .OR. cn_cfg == "ORCA" ) THEN
238               IF( 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      ENDIF
247      !
248      IF(lwxios) THEN
249! define variables in restart file when writing with XIOS
250         CALL iom_set_rstw_var_active('e3t_b')
251         CALL iom_set_rstw_var_active('e3t_n')
252         !                                           ! ----------------------- !
253         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
254            !                                        ! ----------------------- !
255            CALL iom_set_rstw_var_active('tilde_e3t_b')
256            CALL iom_set_rstw_var_active('tilde_e3t_n')
257         END IF
258         !                                           ! -------------!   
259         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
260            !                                        ! ------------ !
261            CALL iom_set_rstw_var_active('hdiv_lf')
262         ENDIF
263         !
264      ENDIF
265      !
266   END SUBROUTINE dom_vvl_init
267
268
269   SUBROUTINE dom_vvl_sf_nxt( kt, Kbb, Kmm, Kaa, kcall ) 
270      !!----------------------------------------------------------------------
271      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
272      !!                   
273      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
274      !!                 tranxt and dynspg routines
275      !!
276      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
277      !!               - z_tilde_case: after scale factor increment =
278      !!                                    high frequency part of horizontal divergence
279      !!                                  + retsoring towards the background grid
280      !!                                  + thickness difusion
281      !!                               Then repartition of ssh INCREMENT proportionnaly
282      !!                               to the "baroclinic" level thickness.
283      !!
284      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
285      !!               - tilde_e3t_a: after increment of vertical scale factor
286      !!                              in z_tilde case
287      !!               - e3(t/u/v)_a
288      !!
289      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
290      !!----------------------------------------------------------------------
291      INTEGER, INTENT( in )           ::   kt             ! time step
292      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
293      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
294      !
295      INTEGER                ::   ji, jj, jk            ! dummy loop indices
296      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
297      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
298      LOGICAL                ::   ll_do_bclinic         ! local logical
299      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
300      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
301      !!----------------------------------------------------------------------
302      !
303      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
304      !
305      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
306      !
307      IF( kt == nit000 ) THEN
308         IF(lwp) WRITE(numout,*)
309         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
310         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
311      ENDIF
312
313      ll_do_bclinic = .TRUE.
314      IF( PRESENT(kcall) ) THEN
315         IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
316      ENDIF
317
318      ! ******************************* !
319      ! After acale factors at t-points !
320      ! ******************************* !
321      !                                                ! --------------------------------------------- !
322      !                                                ! z_star coordinate and barotropic z-tilde part !
323      !                                                ! --------------------------------------------- !
324      !
325      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
326      DO jk = 1, jpkm1
327         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
328         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
329      END DO
330      !
331      IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN   ! z_tilde or layer coordinate !
332         !                                                            ! ------baroclinic part------ !
333         ! I - initialization
334         ! ==================
335
336         ! 1 - barotropic divergence
337         ! -------------------------
338         zhdiv(:,:) = 0._wp
339         zht(:,:)   = 0._wp
340         DO jk = 1, jpkm1
341            zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
342            zht  (:,:) = zht  (:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
343         END DO
344         zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) )
345
346         ! 2 - Low frequency baroclinic horizontal divergence  (z-tilde case only)
347         ! --------------------------------------------------
348         IF( ln_vvl_ztilde ) THEN
349            IF( kt > nit000 ) THEN
350               DO jk = 1, jpkm1
351                  hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:)   &
352                     &          * ( hdiv_lf(:,:,jk) - e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) )
353               END DO
354            ENDIF
355         ENDIF
356
357         ! II - after z_tilde increments of vertical scale factors
358         ! =======================================================
359         tilde_e3t_a(:,:,:) = 0._wp  ! tilde_e3t_a used to store tendency terms
360
361         ! 1 - High frequency divergence term
362         ! ----------------------------------
363         IF( ln_vvl_ztilde ) THEN     ! z_tilde case
364            DO jk = 1, jpkm1
365               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) )
366            END DO
367         ELSE                         ! layer case
368            DO jk = 1, jpkm1
369               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) -   e3t(:,:,jk,Kmm) * ( hdiv(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk)
370            END DO
371         ENDIF
372
373         ! 2 - Restoring term (z-tilde case only)
374         ! ------------------
375         IF( ln_vvl_ztilde ) THEN
376            DO jk = 1, jpk
377               tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk)
378            END DO
379         ENDIF
380
381         ! 3 - Thickness diffusion term
382         ! ----------------------------
383         zwu(:,:) = 0._wp
384         zwv(:,:) = 0._wp
385         DO jk = 1, jpkm1        ! a - first derivative: diffusive fluxes
386            DO jj = 1, jpjm1
387               DO ji = 1, fs_jpim1   ! vector opt.
388                  un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj)           &
389                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj  ,jk) )
390                  vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj)           & 
391                     &            * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji  ,jj+1,jk) )
392                  zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk)
393                  zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk)
394               END DO
395            END DO
396         END DO
397         DO jj = 1, jpj          ! b - correction for last oceanic u-v points
398            DO ji = 1, jpi
399               un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj)
400               vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj)
401            END DO
402         END DO
403         DO jk = 1, jpkm1        ! c - second derivative: divergence of diffusive fluxes
404            DO jj = 2, jpjm1
405               DO ji = fs_2, fs_jpim1   ! vector opt.
406                  tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + (   un_td(ji-1,jj  ,jk) - un_td(ji,jj,jk)    &
407                     &                                          +     vn_td(ji  ,jj-1,jk) - vn_td(ji,jj,jk)    &
408                     &                                            ) * r1_e1e2t(ji,jj)
409               END DO
410            END DO
411         END DO
412         !                       ! d - thickness diffusion transport: boundary conditions
413         !                             (stored for tracer advction and continuity equation)
414         CALL lbc_lnk_multi( 'domvvl', un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp)
415
416         ! 4 - Time stepping of baroclinic scale factors
417         ! ---------------------------------------------
418         ! Leapfrog time stepping
419         ! ~~~~~~~~~~~~~~~~~~~~~~
420         IF( neuler == 0 .AND. kt == nit000 ) THEN
421            z2dt =  rdt
422         ELSE
423            z2dt = 2.0_wp * rdt
424         ENDIF
425         CALL lbc_lnk( 'domvvl', tilde_e3t_a(:,:,:), 'T', 1._wp )
426         tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:)
427
428         ! Maximum deformation control
429         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~
430         ze3t(:,:,jpk) = 0._wp
431         DO jk = 1, jpkm1
432            ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:)
433         END DO
434         z_tmax = MAXVAL( ze3t(:,:,:) )
435         CALL mpp_max( 'domvvl', z_tmax )                 ! max over the global domain
436         z_tmin = MINVAL( ze3t(:,:,:) )
437         CALL mpp_min( 'domvvl', z_tmin )                 ! min over the global domain
438         ! - ML - test: for the moment, stop simulation for too large e3_t variations
439         IF( ( z_tmax >  rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max ) ) THEN
440            IF( lk_mpp ) THEN
441               CALL mpp_maxloc( 'domvvl', ze3t, tmask, z_tmax, ijk_max )
442               CALL mpp_minloc( 'domvvl', ze3t, tmask, z_tmin, ijk_min )
443            ELSE
444               ijk_max = MAXLOC( ze3t(:,:,:) )
445               ijk_max(1) = ijk_max(1) + nimpp - 1
446               ijk_max(2) = ijk_max(2) + njmpp - 1
447               ijk_min = MINLOC( ze3t(:,:,:) )
448               ijk_min(1) = ijk_min(1) + nimpp - 1
449               ijk_min(2) = ijk_min(2) + njmpp - 1
450            ENDIF
451            IF (lwp) THEN
452               WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax
453               WRITE(numout, *) 'at i, j, k=', ijk_max
454               WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin
455               WRITE(numout, *) 'at i, j, k=', ijk_min           
456               CALL ctl_stop( 'STOP', 'MAX( ABS( tilde_e3t_a(:,:,: ) ) / e3t_0(:,:,:) ) too high')
457            ENDIF
458         ENDIF
459         ! - ML - end test
460         ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below
461         tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:),   rn_zdef_max * e3t_0(:,:,:) )
462         tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )
463
464         !
465         ! "tilda" change in the after scale factor
466         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
467         DO jk = 1, jpkm1
468            dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)
469         END DO
470         ! III - Barotropic repartition of the sea surface height over the baroclinic profile
471         ! ==================================================================================
472         ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)
473         ! - ML - baroclinicity error should be better treated in the future
474         !        i.e. locally and not spread over the water column.
475         !        (keep in mind that the idea is to reduce Eulerian velocity as much as possible)
476         zht(:,:) = 0.
477         DO jk = 1, jpkm1
478            zht(:,:)  = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)
479         END DO
480         z_scale(:,:) =  - zht(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
481         DO jk = 1, jpkm1
482            dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
483         END DO
484
485      ENDIF
486
487      IF( ln_vvl_ztilde .OR. ln_vvl_layer )  THEN   ! z_tilde or layer coordinate !
488      !                                           ! ---baroclinic part--------- !
489         DO jk = 1, jpkm1
490            e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kaa) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)
491         END DO
492      ENDIF
493
494      IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN   ! - ML - test: control prints for debuging
495         !
496         IF( lwp ) WRITE(numout, *) 'kt =', kt
497         IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
498            z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )
499            CALL mpp_max( 'domvvl', z_tmax )                             ! max over the global domain
500            IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax
501         END IF
502         !
503         zht(:,:) = 0.0_wp
504         DO jk = 1, jpkm1
505            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
506         END DO
507         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kmm) - zht(:,:) ) )
508         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
509         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t(:,:,:,Kmm)))) =', z_tmax
510         !
511         zht(:,:) = 0.0_wp
512         DO jk = 1, jpkm1
513            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kaa) * tmask(:,:,jk)
514         END DO
515         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kaa) - zht(:,:) ) )
516         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
517         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t(:,:,:,Kaa)))) =', z_tmax
518         !
519         zht(:,:) = 0.0_wp
520         DO jk = 1, jpkm1
521            zht(:,:) = zht(:,:) + e3t(:,:,jk,Kbb) * tmask(:,:,jk)
522         END DO
523         z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssh(:,:,Kbb) - zht(:,:) ) )
524         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
525         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t(:,:,:,Kbb)))) =', z_tmax
526         !
527         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kbb) ) )
528         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
529         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kbb)))) =', z_tmax
530         !
531         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kmm) ) )
532         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
533         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kmm)))) =', z_tmax
534         !
535         z_tmax = MAXVAL( tmask(:,:,1) *  ABS( ssh(:,:,Kaa) ) )
536         CALL mpp_max( 'domvvl', z_tmax )                                ! max over the global domain
537         IF( lwp    ) WRITE(numout, *) kt,' MAXVAL(abs(ssh(:,:,Kaa)))) =', z_tmax
538      END IF
539
540      ! *********************************** !
541      ! After scale factors at u- v- points !
542      ! *********************************** !
543
544      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
545      CALL dom_vvl_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
546
547      ! *********************************** !
548      ! After depths at u- v points         !
549      ! *********************************** !
550
551      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
552      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
553      DO jk = 2, jpkm1
554         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
555         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * 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(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
560      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 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_update( kt, Kbb, Kmm, Kaa )
568      !!----------------------------------------------------------------------
569      !!                ***  ROUTINE dom_vvl_sf_update  ***
570      !!                   
571      !! ** Purpose :  for z tilde case: 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 (ONLY FOR Z TILDE CASE)
576      !!               - reconstruct scale factor at other grid points (interpolate)
577      !!               - recompute depths and water height fields
578      !!
579      !! ** Action  :  - tilde_e3t_(b/n) ready for next time step
580      !!               - Recompute:
581      !!                    e3(u/v)_b       
582      !!                    e3w(:,:,:,Kmm)           
583      !!                    e3(u/v)w_b     
584      !!                    e3(u/v)w_n     
585      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
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      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
593      !
594      INTEGER  ::   ji, jj, jk   ! dummy loop indices
595      REAL(wp) ::   zcoef        ! local scalar
596      !!----------------------------------------------------------------------
597      !
598      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
599      !
600      IF( ln_timing )   CALL timing_start('dom_vvl_sf_update')
601      !
602      IF( kt == nit000 )   THEN
603         IF(lwp) WRITE(numout,*)
604         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step'
605         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
606      ENDIF
607      !
608      ! Time filter and swap of scale factors
609      ! =====================================
610      ! - ML - e3(t/u/v)_b are allready computed in dynnxt.
611      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
612         IF( neuler == 0 .AND. kt == nit000 ) THEN
613            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:)
614         ELSE
615            tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 
616            &         + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) )
617         ENDIF
618         tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:)
619      ENDIF
620      gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)
621      gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)
622
623      e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)
624      e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa)
625      e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa)
626
627      ! Compute all missing vertical scale factor and depths
628      ! ====================================================
629      ! Horizontal scale factor interpolations
630      ! --------------------------------------
631      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
632      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
633     
634      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
635     
636      ! Vertical scale factor interpolations
637      CALL dom_vvl_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
638      CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
639      CALL dom_vvl_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
640      CALL dom_vvl_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
641      CALL dom_vvl_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
642      CALL dom_vvl_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
643
644      ! t- and w- points depth (set the isf depth as it is in the initial step)
645      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
646      gdepw(:,:,1,Kmm) = 0.0_wp
647      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
648      DO jk = 2, jpk
649         DO jj = 1,jpj
650            DO ji = 1,jpi
651              !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
652                                                                 ! 1 for jk = mikt
653               zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
654               gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
655               gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
656                   &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) ) 
657               gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
658            END DO
659         END DO
660      END DO
661
662      ! Local depth and Inverse of the local depth of the water
663      ! -------------------------------------------------------
664      !
665      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
666      DO jk = 2, jpkm1
667         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
668      END DO
669
670      ! write restart file
671      ! ==================
672      IF( lrst_oce  )   CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' )
673      !
674      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_update')
675      !
676   END SUBROUTINE dom_vvl_sf_update
677
678
679   SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout )
680      !!---------------------------------------------------------------------
681      !!                  ***  ROUTINE dom_vvl__interpol  ***
682      !!
683      !! ** Purpose :   interpolate scale factors from one grid point to another
684      !!
685      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
686      !!                - horizontal interpolation: grid cell surface averaging
687      !!                - vertical interpolation: simple averaging
688      !!----------------------------------------------------------------------
689      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
690      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
691      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
692      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
693      !
694      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
695      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
696      !!----------------------------------------------------------------------
697      !
698      IF(ln_wd_il) THEN
699        zlnwd = 1.0_wp
700      ELSE
701        zlnwd = 0.0_wp
702      END IF
703      !
704      SELECT CASE ( pout )    !==  type of interpolation  ==!
705         !
706      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
707         DO jk = 1, jpk
708            DO jj = 1, jpjm1
709               DO ji = 1, fs_jpim1   ! vector opt.
710                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
711                     &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
712                     &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
713               END DO
714            END DO
715         END DO
716         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'U', 1._wp )
717         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
718         !
719      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
720         DO jk = 1, jpk
721            DO jj = 1, jpjm1
722               DO ji = 1, fs_jpim1   ! vector opt.
723                  pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
724                     &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
725                     &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
726               END DO
727            END DO
728         END DO
729         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'V', 1._wp )
730         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
731         !
732      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
733         DO jk = 1, jpk
734            DO jj = 1, jpjm1
735               DO ji = 1, fs_jpim1   ! vector opt.
736                  pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
737                     &                       *    r1_e1e2f(ji,jj)                                                  &
738                     &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
739                     &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
740               END DO
741            END DO
742         END DO
743         CALL lbc_lnk( 'domvvl', pe3_out(:,:,:), 'F', 1._wp )
744         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
745         !
746      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
747         !
748         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
749         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
750!!gm BUG? use here wmask in case of ISF ?  to be checked
751         DO jk = 2, jpk
752            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
753               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
754               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
755               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
756         END DO
757         !
758      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
759         !
760         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
761         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
762!!gm BUG? use here wumask in case of ISF ?  to be checked
763         DO jk = 2, jpk
764            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
765               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
766               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
767               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
768         END DO
769         !
770      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
771         !
772         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
773         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
774!!gm BUG? use here wvmask in case of ISF ?  to be checked
775         DO jk = 2, jpk
776            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
777               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
778               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
779               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
780         END DO
781      END SELECT
782      !
783   END SUBROUTINE dom_vvl_interpol
784
785
786   SUBROUTINE dom_vvl_rst( kt, Kbb, Kmm, cdrw )
787      !!---------------------------------------------------------------------
788      !!                   ***  ROUTINE dom_vvl_rst  ***
789      !!                     
790      !! ** Purpose :   Read or write VVL file in restart file
791      !!
792      !! ** Method  :   use of IOM library
793      !!                if the restart does not contain vertical scale factors,
794      !!                they are set to the _0 values
795      !!                if the restart does not contain vertical scale factors increments (z_tilde),
796      !!                they are set to 0.
797      !!----------------------------------------------------------------------
798      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
799      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
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'   , ssh(:,:,Kmm), 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(:,:,:,Kbb), ldxios = lrxios )
822               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
823               ! needed to restart if land processor not computed
824               IF(lwp) write(numout,*) 'dom_vvl_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
825               WHERE ( tmask(:,:,:) == 0.0_wp ) 
826                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
827                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
828               END WHERE
829               IF( neuler == 0 ) THEN
830                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
831               ENDIF
832            ELSE IF( id1 > 0 ) THEN
833               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) 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(:,:,:,Kbb), ldxios = lrxios )
837               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
838               neuler = 0
839            ELSE IF( id2 > 0 ) THEN
840               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kbb) 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(:,:,:,Kmm), ldxios = lrxios )
844               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
845               neuler = 0
846            ELSE
847               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : e3t(:,:,:,Kmm) 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(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
852                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
853                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
854               END DO
855               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
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(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
893                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
894                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
895                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
896                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
897               ELSE
898                  ! if not test case
899                  ssh(:,:,Kmm) = -ssh_ref
900                  ssh(:,:,Kbb) = -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                           ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
906                           ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
907                        ENDIF
908                     ENDDO
909                  ENDDO
910               ENDIF !If test case else
911
912               ! Adjust vertical metrics for all wad
913               DO jk = 1, jpk
914                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
915                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
916                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
917               END DO
918               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
919
920               DO ji = 1, jpi
921                  DO jj = 1, jpj
922                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
923                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
924                     ENDIF
925                  END DO
926               END DO 
927               !
928            ELSE
929               !
930               ! usr_def_istate called here only to get sshb, that is needed to initialize e3t(Kbb) and e3t(Kmm)
931               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
932               ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn)
933               !
934               DO jk=1,jpk
935                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
936                    &                              / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
937                    &              + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )   ! make sure e3t_b != 0 on land points
938               END DO
939               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
940               ssh(:,:  ,Kmm) = ssh(:,:  ,Kbb)   ! needed later for gde3w
941               !
942            END IF           ! end of ll_wd edits
943
944            IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN
945               tilde_e3t_b(:,:,:) = 0._wp
946               tilde_e3t_n(:,:,:) = 0._wp
947               IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp
948            END IF
949         ENDIF
950         !
951      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
952         !                                   ! ===================
953         IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----'
954         IF( lwxios ) CALL iom_swap(      cwxios_context          )
955         !                                           ! --------- !
956         !                                           ! all cases !
957         !                                           ! --------- !
958         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
959         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
960         !                                           ! ----------------------- !
961         IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN  ! z_tilde and layer cases !
962            !                                        ! ----------------------- !
963            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:), ldxios = lwxios)
964            CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:), ldxios = lwxios)
965         END IF
966         !                                           ! -------------!   
967         IF( ln_vvl_ztilde ) THEN                    ! z_tilde case !
968            !                                        ! ------------ !
969            CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:), ldxios = lwxios)
970         ENDIF
971         !
972         IF( lwxios ) CALL iom_swap(      cxios_context          )
973      ENDIF
974      !
975   END SUBROUTINE dom_vvl_rst
976
977
978   SUBROUTINE dom_vvl_ctl
979      !!---------------------------------------------------------------------
980      !!                  ***  ROUTINE dom_vvl_ctl  ***
981      !!               
982      !! ** Purpose :   Control the consistency between namelist options
983      !!                for vertical coordinate
984      !!----------------------------------------------------------------------
985      INTEGER ::   ioptio, ios
986      !!
987      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
988         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
989         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
990      !!----------------------------------------------------------------------
991      !
992      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
993901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
994      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
995902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
996      IF(lwm) WRITE ( numond, nam_vvl )
997      !
998      IF(lwp) THEN                    ! Namelist print
999         WRITE(numout,*)
1000         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
1001         WRITE(numout,*) '~~~~~~~~~~~'
1002         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
1003         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
1004         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
1005         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
1006         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
1007         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
1008         WRITE(numout,*) '      !'
1009         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
1010         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
1011         IF( ln_vvl_ztilde_as_zstar ) THEN
1012            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
1013            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
1014            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
1015            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
1016            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
1017            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
1018         ELSE
1019            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
1020            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
1021         ENDIF
1022         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
1023      ENDIF
1024      !
1025      ioptio = 0                      ! Parameter control
1026      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
1027      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
1028      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
1029      IF( ln_vvl_layer           )   ioptio = ioptio + 1
1030      !
1031      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
1032      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
1033      !
1034      IF(lwp) THEN                   ! Print the choice
1035         WRITE(numout,*)
1036         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
1037         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
1038         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
1039         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
1040      ENDIF
1041      !
1042#if defined key_agrif
1043      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
1044#endif
1045      !
1046   END SUBROUTINE dom_vvl_ctl
1047
1048   !!======================================================================
1049END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.