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

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

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

Last change on this file since 10001 was 10001, checked in by gm, 6 years ago

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