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/VORTEX/MY_SRC – NEMO

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