source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DOM/domvvl.F90 @ 10009

Last change on this file since 10009 was 10009, checked in by gm, 2 years ago

#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh

  • Property svn:keywords set to Id
File size: 36.0 KB
Line 
1MODULE domvvl
2   !!======================================================================
3   !!                       ***  MODULE domvvl   ***
4   !! Ocean :
5   !!======================================================================
6   !! History :  5.0  !  2018-07  (G. Madec)  Flux Form with Kinetic Energy conservation
7   !!                             ==>>>  here z* and s* only (no z-tilde)
8   
9   ! 1- remove   z-tilde          ==>>>  pure z-star (or s-star)
10   ! 2- remove   dom_vvl_interpol 
11   
12   !!----------------------------------------------------------------------
13
14   !!----------------------------------------------------------------------
15   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
16   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
17   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
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_swp     ! called by step.F90
41   PUBLIC  ssh2e3_before      ! ...
42   PUBLIC  ssh2e3_now         ! ...
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)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_b, te3t_n    ! baroclinic scale factors
61   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: te3t_a, dte3t_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!!gm add
67!!gm
68
69
70
71   !! * Substitutions
72#  include "vectopt_loop_substitute.h90"
73   !!----------------------------------------------------------------------
74   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
75   !! $Id$
76   !! Software governed by the CeCILL licence     (./LICENSE)
77   !!----------------------------------------------------------------------
78CONTAINS
79
80   INTEGER FUNCTION dom_vvl_alloc()
81      !!----------------------------------------------------------------------
82      !!                ***  FUNCTION dom_vvl_alloc  ***
83      !!----------------------------------------------------------------------
84      !
85      IF( ln_vvl_zstar )   dom_vvl_alloc = 0
86      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
87         ALLOCATE( te3t_b(jpi,jpj,jpk)  , te3t_n(jpi,jpj,jpk) , te3t_a(jpi,jpj,jpk) ,   &
88            &      dte3t_a(jpi,jpj,jpk) , un_td  (jpi,jpj,jpk)     , vn_td  (jpi,jpj,jpk)     ,   &
89            &      STAT = dom_vvl_alloc        )
90         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
91         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
92         un_td = 0._wp
93         vn_td = 0._wp
94      ENDIF
95      IF( ln_vvl_ztilde ) THEN
96         ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc )
97         IF( lk_mpp             )   CALL mpp_sum ( dom_vvl_alloc )
98         IF( dom_vvl_alloc /= 0 )   CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays')
99      ENDIF
100      !
101   END FUNCTION dom_vvl_alloc
102
103
104   SUBROUTINE dom_vvl_init
105      !!----------------------------------------------------------------------
106      !!                ***  ROUTINE dom_vvl_init  ***
107      !!                   
108      !! ** Purpose :  Initialization of all scale factors, depths
109      !!               and water column heights
110      !!
111      !! ** Method  :  - use restart file and/or initialize
112      !!               - interpolate scale factors
113      !!
114      !! ** Action  : - e3t_(n/b) and te3t_(n/b)
115      !!              - Regrid: e3(u/v)_n
116      !!                        e3(u/v)_b       
117      !!                        e3w_n           
118      !!                        e3(u/v)w_b     
119      !!                        e3(u/v)w_n     
120      !!                        gdept_n, gdepw_n and gde3w_n
121      !!              - h(t/u/v)_0
122      !!              - frq_rst_e3t and frq_rst_hdv
123      !!
124      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
125      !!----------------------------------------------------------------------
126      !!----------------------------------------------------------------------
127      !
128      IF(lwp) WRITE(numout,*)
129      IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated'
130      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
131      !
132      CALL dom_vvl_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
133      !
134      !                    ! Allocate module arrays
135      IF( dom_vvl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' )
136      !
137      !                    ! Read or initialize ssh(Nbb) & ssh(Nnn)
138      CALL dom_vvl_rst( nit000, 'READ' )
139      !
140      !                    !== Set of all other vertical mesh fields  ==!  (now and before)     
141      !
142      !                          !* BEFORE fields :
143      CALL ssh2e3_before               ! set:      hu , hv , r1_hu, r1_hv
144      !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw        (from 1 to jpkm1)
145      !
146      !                                ! set jpk level one to the e3._0 values
147      e3t_b(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_b(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_b(:,:,jpk) =  e3v_0(:,:,jpk)
148      e3w_b(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_b(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_b(:,:,jpk) = e3vw_0(:,:,jpk)
149      !
150      !                          !* NOW fields :
151      CALL ssh2e3_now                  ! set: ht , hu , hv , r1_hu, r1_hv
152      !                                !      e3t, e3w, e3u, e3uw, e3v, e3vw, e3f   (from 1 to jpkm1)
153      !                                !      gdept_n, gdepw_n, gde3w_n
154!!gm issue?   gdept_n, gdepw_n, gde3w_n never defined at jpk
155      !
156      !                                ! set one for all last level to the e3._0 value
157      e3t_n(:,:,jpk) = e3t_0(:,:,jpk)  ;   e3u_n(:,:,jpk) =  e3w_0(:,:,jpk)  ;   e3v_n(:,:,jpk) =  e3v_0(:,:,jpk)
158      e3w_n(:,:,jpk) = e3w_0(:,:,jpk)  ;  e3uw_n(:,:,jpk) = e3uw_0(:,:,jpk)  ;  e3vw_n(:,:,jpk) = e3vw_0(:,:,jpk)
159      e3f_n(:,:,jpk) = e3f_0(:,:,jpk)
160      !
161      !                          !* AFTER fields : (last level for OPA, 3D required for AGRIF initialisation)
162      e3t_a(:,:,:) = e3t_n(:,:,:)   ;   e3u_a(:,:,:) = e3u_n(:,:,:)   ;   e3v_a(:,:,:) = e3v_n(:,:,:)
163     
164!!gm           
165!!gm        ===>>>>   below:  some issues to think about !!!
166!!gm 
167!!gm   fmask definition checked (0 or 1 nothing else)
168!!        in z-tilde or ALE    e3._0  should be the time varying fields step forward with an euler scheme
169!!gm   e3w_b  & gdept_b are not used   except its update in agrif 
170!!     and used to compute before slope of surface in dynldf_iso ==>>>  remove it !!!!
171!!         NB: in triads on TRA, gdept_n  is used !!!!   BUG?
172!!gm   e3f_n  almost not used  ===>>>>  verify whether it can be removed or not...
173!!       verify the use of wumask & wvmask   mau be replaced by a multiplication by umask(k)*umask(k+1)  ???
174!!
175!!gm  ISF case to be checked by Pierre Mathiot
176!!
177!!gm  setting of e3._a  for agrif....  TO BE CHECKED
178
179      !
180      IF(lwxios) THEN      ! define variables in restart file when writing with XIOS
181         CALL iom_set_rstw_var_active('e3t_b')
182         CALL iom_set_rstw_var_active('e3t_n')
183      ENDIF
184      !
185   END SUBROUTINE dom_vvl_init
186
187
188   SUBROUTINE dom_vvl_sf_nxt( kt, kcall ) 
189      !!----------------------------------------------------------------------
190      !!                ***  ROUTINE dom_vvl_sf_nxt  ***
191      !!                   
192      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
193      !!                 tranxt and dynspg routines
194      !!
195      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
196      !!               - z_tilde_case: after scale factor increment =
197      !!                                    high frequency part of horizontal divergence
198      !!                                  + retsoring towards the background grid
199      !!                                  + thickness difusion
200      !!                               Then repartition of ssh INCREMENT proportionnaly
201      !!                               to the "baroclinic" level thickness.
202      !!
203      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
204      !!               - te3t_a: after increment of vertical scale factor
205      !!                              in z_tilde case
206      !!               - e3(t/u/v)_a
207      !!
208      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
209      !!----------------------------------------------------------------------
210      INTEGER, INTENT( in )           ::   kt      ! time step
211      INTEGER, INTENT( in ), OPTIONAL ::   kcall   ! optional argument indicating call sequence
212      !
213      INTEGER ::   ji, jj, jk   ! dummy loop indices
214      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
215      !!----------------------------------------------------------------------
216      !
217      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
218      !
219      IF( ln_timing )   CALL timing_start('dom_vvl_sf_nxt')
220      !
221      IF( kt == nit000 ) THEN
222         IF(lwp) WRITE(numout,*)
223         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors'
224         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
225      ENDIF
226
227      !                             ! ------------------!
228      !                             ! z_star coordinate !     (and barotropic z-tilde part)
229      !                             ! ------------------!
230      !
231      !                                   !==  after ssh  ==!  (u- and v-points)
232      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
233         zsshu_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji+1,jj,Naa) ) * ssumask(ji,jj)
234         zsshv_h(ji,jj) = 0.5_wp * ( ssh(ji,jj,Naa) + ssh(ji,jj+1,Naa) ) * ssvmask(ji,jj)
235      END DO             ;   END DO     
236      CALL lbc_lnk_multi( zsshu_h(:,:), 'U', 1._wp , zsshv_h(:,:), 'V', 1._wp )
237      !
238      !                                   !==  after depths and its inverse  ==!
239         hu_a(:,:) = hu_0(:,:) + zsshu_h(:,:)
240         hv_a(:,:) = hv_0(:,:) + zsshv_h(:,:)
241      r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )
242      r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )
243      !
244      !                                   !==  after scale factors  ==!  (e3t , e3u , e3v)
245      zssht_h(:,:) = ssh    (:,:,Naa) * r1_ht_0(:,:)     ! t-point
246      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! u-point
247      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! v-point
248      DO jk = 1, jpkm1
249         e3t_a(:,:,jk) = e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * tmask(:,:,jk) )
250         e3u_a(:,:,jk) = e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
251         e3v_a(:,:,jk) = e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
252      END DO
253      !
254      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_nxt')
255      !
256   END SUBROUTINE dom_vvl_sf_nxt
257
258
259   SUBROUTINE dom_vvl_sf_swp( kt )
260      !!----------------------------------------------------------------------
261      !!                ***  ROUTINE dom_vvl_sf_swp  ***
262      !!                   
263      !! ** Purpose :  compute time filter and swap of scale factors
264      !!               compute all depths and related variables for next time step
265      !!               write outputs and restart file
266      !!
267      !! ** Method  :  - swap of e3t with trick for volume/tracer conservation
268      !!               - reconstruct scale factor at other grid points (interpolate)
269      !!               - recompute depths and water height fields
270      !!
271      !! ** Action  :  - e3t_(b/n), te3t_(b/n) and e3(u/v)_n ready for next time step
272      !!               - Recompute:
273      !!                    e3(u/v)_b       
274      !!                    e3w_n           
275      !!                    e3(u/v)w_b     
276      !!                    e3(u/v)w_n     
277      !!                    gdept_n, gdepw_n  and gde3w_n
278      !!                    h(u/v) and h(u/v)r
279      !!
280      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
281      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
282      !!----------------------------------------------------------------------
283      INTEGER, INTENT( in ) ::   kt   ! time step
284      !
285      INTEGER  ::   ji, jj, jk    ! dummy loop indices
286      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h
287      !!----------------------------------------------------------------------
288      !
289      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
290      !
291      IF( ln_timing )   CALL timing_start('dom_vvl_sf_swp')
292      !
293      IF( kt == nit000 )   THEN
294         IF(lwp) WRITE(numout,*)
295         IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors'
296         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~   - interpolate scale factors and compute depths for next time step'
297      ENDIF
298      !
299      ! Swap and Compute all missing vertical scale factor and depths
300      ! =============================================================
301      ! - ML - e3u_b and e3v_b are allready computed in dynnxt
302      ! - JC - hu_b, hv_b, hur_b, hvr_b also
303      !
304      ! - GM - to be updated :   e3f_n,  e3w_n , e3uw_n , e3vw_n
305      !                                  e3w_b , e3uw_b , e3vw_b
306      !                          gdept_n , gdepw_n , gde3w_n
307      !                          ht_n
308      !        to be swap    :   hu_a , hv_a  , r1_hu_a , r1_hv_a
309      !
310      ! Local depth and Inverse of the local depth of the water
311      ! -------------------------------------------------------
312      !                       ! swap of depth and scale factors
313      !                       ! ===============================
314      DO jk = 1, jpkm1                          ! swap n--> a 
315         gdept_b(:,:,jk) = gdept_n(:,:,jk)         ! depth at t and w
316         gdepw_b(:,:,jk) = gdepw_n(:,:,jk)
317         e3t_n  (:,:,jk) = e3t_a  (:,:,jk)         ! e3t, e3u, e3v
318         e3u_n  (:,:,jk) = e3u_a  (:,:,jk)
319         e3v_n  (:,:,jk) = e3v_a  (:,:,jk)
320      END DO
321      ht_n(:,:) = ht_0(:,:) + ssh(:,:,Nnn)            ! ocean thickness
322      !
323      hu_n(:,:) = hu_a(:,:)   ;   r1_hu_n(:,:) = r1_hu_a(:,:)
324      hv_n(:,:) = hv_a(:,:)   ;   r1_hv_n(:,:) = r1_hv_a(:,:)
325      !
326      !                    !==  before :
327      !                                            !* ssh at u- and v-points)
328      DO jj = 2, jpjm1   ;   DO ji = 2, jpim1
329         zsshu_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
330         zsshv_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
331      END DO             ;   END DO     
332      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp )
333      !
334      !                                            !*  e3w_b , e3uw_b , e3vw_b
335      zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)     ! w-point
336      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! uw-point
337      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! vw-point
338      DO jk = 1, jpkm1
339          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
340         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * umask(:,:,jk) )
341         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * vmask(:,:,jk) )
342      END DO
343      !     
344      !                    !==  now    :
345      !                                            !* ssh at u- and v-points)
346      DO jj = 1, jpjm1   ;   DO ji = 1, jpim1            ! start from 1 for f-point
347         zsshu_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
348         zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
349         zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
350            &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
351      END DO             ;   END DO     
352      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp )     
353      !
354      !                                            !* e3w_n , e3uw_n , e3vw_n, e3f_n
355      zssht_h(:,:) = ssh    (:,:,Nnn) * r1_ht_0(:,:)     ! t- & w-point
356      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)     ! uw-point
357      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)     ! vw-point
358      zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)     ! f-point
359      DO jk = 1, jpkm1
360          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
361         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) )
362         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) )
363          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) *  fmask(:,:,jk) )
364      END DO     
365      !
366      zssht_h(:,:) = 1._wp + zssht_h(:,:)               ! t-point
367      !
368      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
369         DO jk = 1, jpkm1
370            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
371            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
372            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh    (:,:,Nnn)
373         END DO
374      ELSE                    ! no ISF cavities
375         DO jk = 1, jpkm1
376            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
377            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
378            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh    (:,:,Nnn)
379         END DO
380      ENDIF
381      !
382      ! write restart file
383      ! ==================
384      IF( lrst_oce  )   CALL dom_vvl_rst( kt, 'WRITE' )
385      !
386      IF( ln_timing )   CALL timing_stop('dom_vvl_sf_swp')
387      !
388   END SUBROUTINE dom_vvl_sf_swp
389
390
391   SUBROUTINE dom_vvl_rst( kt, cdrw )
392      !!---------------------------------------------------------------------
393      !!                   ***  ROUTINE dom_vvl_rst  ***
394      !!                     
395      !! ** Purpose :   Read or write VVL file in restart file
396      !!
397      !! ** Method  :   use of IOM library
398      !!                if the restart does not contain vertical scale factors,
399      !!                they are set to the _0 values
400      !!                if the restart does not contain vertical scale factors increments (z_tilde),
401      !!                they are set to 0.
402      !!----------------------------------------------------------------------
403      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
404      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
405      !
406      INTEGER ::   ji, jj, jk
407      INTEGER ::   id1, id2, id3, id4, id5     ! local integers
408      !!----------------------------------------------------------------------
409      !
410      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
411         !                                   ! ===============
412         IF( ln_rstart ) THEN                   !* Read the restart file
413            CALL rst_read_open                  !  open the restart file if necessary
414            !
415            id1 = iom_varid( numror, 'sshb'      , ldstop = .FALSE. )
416            id2 = iom_varid( numror, 'sshn'      , ldstop = .FALSE. )
417            !
418            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
419               IF(lwp) write(numout,*) 'dom_vvl_rst : both sshb and sshn found in restart files'
420               !
421!!gm  Question: use jpdom_data above to read data over jpi x jpj    (like is dom_hgr_read and dom_zgr_read)
422!!              so that it will work with land processor suppression
423!               CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Nnn), ldxios = lrxios    )
424!               CALL iom_get( numror, jpdom_autoglo, 'sshb'   , ssh(:,:,Nbb), ldxios = lrxios    )
425!!gm
426               CALL iom_get( numror, jpdom_data, 'sshn'   , ssh(:,:,Nnn), ldxios = lrxios    )
427               CALL iom_get( numror, jpdom_data, 'sshb'   , ssh(:,:,Nbb), ldxios = lrxios    )
428!!gm end
429               IF( l_1st_euler )   ssh(:,:,Nbb) = ssh(:,:,Nnn)
430            ELSE IF( id1 > 0 ) THEN
431               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshn not found in restart files'
432               IF(lwp) write(numout,*) '   set ssh(Nnn) = ssh(Nbb)  and force l_1st_euler = true'
433!!gm               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t_b(:,:,:), ldxios = lrxios )
434               CALL iom_get( numror, jpdom_data, 'sshb', ssh(:,:,Nbb), ldxios = lrxios )
435               ssh(:,:,Nnn) = ssh(:,:,Nbb)
436               l_1st_euler = .TRUE.
437            ELSE IF( id2 > 0 ) THEN
438               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb not found in restart files'
439               IF(lwp) write(numout,*) 'set ssh(Nbb) = ssh(Nnn)  and force l_1st_euler = true'
440               CALL iom_get( numror, jpdom_data, 'sshn', ssh(:,:,Nnn), ldxios = lrxios )
441               ssh(:,:,Nbb) = ssh(:,:,Nnn)
442               l_1st_euler = .TRUE.
443            ELSE
444               IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : sshb and sshn not found in restart file'
445               IF(lwp) write(numout,*) 'set ssh(Nbb) = ssh(Nnn) = 0  and force l_1st_euler = true'
446               ssh(:,:,Nbb) = 0._wp
447               ssh(:,:,Nnn) = 0._wp
448               l_1st_euler = .TRUE.
449            ENDIF
450         ELSE                                   !* Initialize at "rest"
451            !
452            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
453               !
454               IF( cn_cfg == 'wad' ) THEN             ! Wetting and drying test case
455                  CALL usr_def_istate( gdept_b, tmask, tsb, ub, vb, ssh(:,:,Nbb)  )
456                  tsn(:,:,:,:) = tsb(:,:,:,:)            ! set now values from to before ones
457                  ssh(:,:,Nnn) = ssh(:,:,Nbb)
458                  un (:,:,:)   = ub (:,:,:)
459                  vn (:,:,:)   = vb (:,:,:)
460               ELSE                                   ! Not the test case
461                  ssh(:,:,Nnn) = -ssh_ref
462                  ssh(:,:,Nbb) = -ssh_ref
463                  !
464                  DO jj = 1, jpj
465                     DO ji = 1, jpi
466                        IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN  ! if total depth is less than min depth
467                           ssh(ji,jj,Nbb) = rn_wdmin1 - (ht_0(ji,jj) )
468                           ssh(ji,jj,Nnn) = rn_wdmin1 - (ht_0(ji,jj) )
469                           ssh(ji,jj,Naa) = rn_wdmin1 - (ht_0(ji,jj) )
470                        ENDIF
471                     END DO
472                  END DO
473               ENDIF     ! If test case else
474               !
475               DO ji = 1, jpi
476                  DO jj = 1, jpj
477                     IF ( ht_0(ji,jj) /= 0._wp .AND. NINT( ssmask(ji,jj) ) == 1 ) THEN
478                       CALL ctl_stop( 'dom_vvl_rst: ht_0 must be positive at potentially wet points' )
479                     ENDIF
480                  END DO
481               END DO 
482               !
483            ELSE
484               !
485               ! Just to read set ssh in fact, called latter once vertical grid is set up:
486!               CALL usr_def_istate( gdept_0, tmask, tsb, ub, vb, ssh(:,:,Nbb)  )
487               ssh(:,:,Nnn) = 0._wp
488               ssh(:,:,Nbb) = 0._wp
489               !
490            END IF
491            !
492         ENDIF
493         !
494      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
495         !                                   ! ===================
496
497!!gm      DO NOTHING,   ssh(Nbb) and ssh(Nnn)  are written in restart.F90
498
499      ENDIF
500      !
501   END SUBROUTINE dom_vvl_rst
502
503
504   SUBROUTINE dom_vvl_ctl
505      !!---------------------------------------------------------------------
506      !!                  ***  ROUTINE dom_vvl_ctl  ***
507      !!               
508      !! ** Purpose :   Control the consistency between namelist options
509      !!                for vertical coordinate
510      !!----------------------------------------------------------------------
511      INTEGER ::   ioptio, ios
512      !!
513      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
514         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
515         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
516      !!----------------------------------------------------------------------
517      !
518      REWIND( numnam_ref )              ! Namelist nam_vvl in reference namelist :
519      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
520901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist', lwp )
521      REWIND( numnam_cfg )              ! Namelist nam_vvl in configuration namelist : Parameters of the run
522      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
523902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist', lwp )
524      IF(lwm) WRITE ( numond, nam_vvl )
525      !
526      IF(lwp) THEN                    ! Namelist print
527         WRITE(numout,*)
528         WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate'
529         WRITE(numout,*) '~~~~~~~~~~~'
530         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
531         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
532         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
533         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
534         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
535         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
536         WRITE(numout,*) '      !'
537         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
538         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
539         IF( ln_vvl_ztilde_as_zstar ) THEN
540            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
541            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
542            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
543            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
544            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
545            WRITE(numout,*) '                         rn_lf_cutoff   = 1/rn_Dt'
546         ELSE
547            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
548            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
549         ENDIF
550         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
551      ENDIF
552      !
553
554!!gm
555      IF ( ln_vvl_ztilde .OR. ln_vvl_ztilde_as_zstar )   CALL ctl_stop( 'z-tilde NOT available in this branch' )
556!!gm
557
558      !
559      ioptio = 0                      ! Parameter control
560      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
561      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
562      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
563      IF( ln_vvl_layer           )   ioptio = ioptio + 1
564      !
565      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
566      IF( .NOT. ln_vvl_zstar .AND. ln_isf ) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' )
567      !
568      IF(lwp) THEN                   ! Print the choice
569         WRITE(numout,*)
570         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
571         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
572         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
573         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
574      ENDIF
575      !
576#if defined key_agrif
577      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
578#endif
579      !
580   END SUBROUTINE dom_vvl_ctl
581
582   
583   SUBROUTINE ssh2e3_now
584      !!----------------------------------------------------------------------
585      !!                  ***  ROUTINE ssh2e3_now  ***
586      !!----------------------------------------------------------------------
587      INTEGER ::   ji, jj, jk
588      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h, zsshf_h
589      !!----------------------------------------------------------------------
590      !
591      !                             !==  ssh at u- and v-points  ==!
592      !
593      DO jj = 1, jpjm1                    ! start from 1 due to f-point
594         DO ji = 1, jpim1
595            zsshu_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji+1,jj  ,Nnn) ) * ssumask(ji,jj)
596            zsshv_h(ji,jj) = 0.50_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn) ) * ssvmask(ji,jj)
597            zsshf_h(ji,jj) = 0.25_wp * ( ssh(ji  ,jj,Nnn) + ssh(ji  ,jj+1,Nnn)   & 
598               &                       + ssh(ji+1,jj,Nnn) + ssh(ji+1,jj+1,Nnn) ) * ssfmask(ji,jj)
599         END DO
600      END DO     
601      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp , zsshf_h(:,:),'F', 1._wp )
602      !
603      !                             !==  ht, hu and hv  == !   (and their inverse)
604      !
605      ht_n   (:,:) = ht_0(:,:) +  ssh   (:,:,Nnn)
606      hu_n   (:,:) = hu_0(:,:) + zsshu_h(:,:)
607      hv_n   (:,:) = hv_0(:,:) + zsshv_h(:,:)
608      r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )    ! ss mask mask due to ISF
609      r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )
610      !     
611      !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
612      !
613      zssht_h(:,:) =  ssh   (:,:,Nnn) * r1_ht_0(:,:)
614      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
615      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
616      zsshf_h(:,:) = zsshf_h(:,:)     * r1_hf_0(:,:)
617      !
618      !                             !==  e3t, e3w  ,  e3u, e3uw ,  e3v, e3vw  , and e3f  ==!
619      !     
620      DO jk = 1, jpkm1
621          e3t_n(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) *      tmask(:,:,jk)                     )
622          e3w_n(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
623          !
624          e3u_n(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) *  umask(:,:,jk) )
625         e3uw_n(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h(:,:) * wumask(:,:,jk) )
626         !
627          e3v_n(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) *  vmask(:,:,jk) )
628         e3vw_n(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h(:,:) * wvmask(:,:,jk) )
629          !
630          e3f_n(:,:,jk) =  e3f_0(:,:,jk) * ( 1._wp + zsshf_h(:,:) * fmask(:,:,jk) )
631      END DO
632      !     
633      !                             !== depth of t- and w-points  ==!
634      !
635      zssht_h(:,:) = 1._wp + zssht_h(:,:)     ! = 1 + ssh(Nnn) / ht_0
636      !
637      IF( ln_isfcav ) THEN    ! ISF cavities : ssh scaling not applied over the iceshelf thickness
638         DO jk = 1, jpkm1
639            gdept_n(:,:,jk) = ( gdept_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
640            gdepw_n(:,:,jk) = ( gdepw_0(:,:,jk) - risfdep(:,:) ) * zssht_h(:,:) + risfdep(:,:)
641            gde3w_n(:,:,jk) =   gdept_n(:,:,jk) - ssh(:,:,Nnn)
642         END DO
643      ELSE                    ! no ISF cavities
644!!gm BUG ???    gdept should be updated down to the ocean floor !  ===>> jpk NOT jpkm1 !!!
645         DO jk = 1, jpkm1
646            gdept_n(:,:,jk) = gdept_0(:,:,jk) * zssht_h(:,:)
647            gdepw_n(:,:,jk) = gdepw_0(:,:,jk) * zssht_h(:,:)
648            gde3w_n(:,:,jk) = gdept_n(:,:,jk) - ssh(:,:,Nnn)
649         END DO
650      ENDIF
651      !
652   END SUBROUTINE ssh2e3_now
653   
654   
655   SUBROUTINE ssh2e3_before
656      !!----------------------------------------------------------------------
657      !!                  ***  ROUTINE ssh2e3_before  ***
658      !!----------------------------------------------------------------------
659      INTEGER ::   ji, jj, jk
660      REAL(wp), DIMENSION(jpi,jpj) ::   zssht_h, zsshu_h, zsshv_h
661      !!----------------------------------------------------------------------
662      !
663      !                             !==  ssh at u- and v-points  ==!
664      DO jj = 2, jpjm1
665         DO ji = 2, jpim1
666            zsshu_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji+1,jj  ,Nbb) ) * ssumask(ji,jj)
667            zsshv_h(ji,jj) = 0.5_wp  * ( ssh(ji,jj,Nbb) + ssh(ji  ,jj+1,Nbb) ) * ssvmask(ji,jj)
668         END DO
669      END DO     
670      CALL lbc_lnk_multi( zsshu_h(:,:),'U', 1._wp , zsshv_h(:,:),'V', 1._wp )
671      !
672      !                             !==  ht, hu and hv  == !   (and their inverse)
673         hu_b(:,:) = hu_0(:,:) + zsshu_h(:,:)
674         hv_b(:,:) = hv_0(:,:) + zsshv_h(:,:)
675      r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )
676      r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )
677      !
678      !     
679      !                             !==  ssh / h  factor at t-, u- ,v- & f-points  ==!
680      zssht_h(:,:) = ssh    (:,:,Nbb) * r1_ht_0(:,:)
681      zsshu_h(:,:) = zsshu_h(:,:)     * r1_hu_0(:,:)
682      zsshv_h(:,:) = zsshv_h(:,:)     * r1_hv_0(:,:)
683      !
684      !                             !==  e3t, e3w  ,  e3u, e3uw , and  e3v, e3vw  ==!
685      DO jk = 1, jpkm1
686          e3t_b(:,:,jk) =  e3t_0(:,:,jk) * ( 1._wp + zssht_h(:,:) *      tmask(:,:,jk)                     )
687          e3w_b(:,:,jk) =  e3w_0(:,:,jk) * ( 1._wp + zssht_h(:,:) * MAX( tmask(:,:,jk) , tmask(:,:,jk+1) ) )
688          !
689          e3u_b(:,:,jk) =  e3u_0(:,:,jk) * ( 1._wp + zsshu_h  (:,:) *  umask(:,:,jk) )
690         e3uw_b(:,:,jk) = e3uw_0(:,:,jk) * ( 1._wp + zsshu_h  (:,:) * wumask(:,:,jk) )
691          !
692          e3v_b(:,:,jk) =  e3v_0(:,:,jk) * ( 1._wp + zsshv_h  (:,:) *  vmask(:,:,jk) )
693         e3vw_b(:,:,jk) = e3vw_0(:,:,jk) * ( 1._wp + zsshv_h  (:,:) * wvmask(:,:,jk) )
694      END DO
695      !   
696   END SUBROUTINE ssh2e3_before
697   
698   !!======================================================================
699END MODULE domvvl
Note: See TracBrowser for help on using the repository browser.