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.
domqe.F90 in NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM – NEMO

source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DOM/domqe.F90 @ 12482

Last change on this file since 12482 was 12482, checked in by techene, 4 years ago

new reference without ztilde, duplicated modules and routines to be modified from zstar MLF to zstar LF
File size: 36.2 KB
Line 
1MODULE domqe
2   !!======================================================================
3   !!                       ***  MODULE domqe   ***
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   !!            4.x  ! 2020-02  (G. Madec, S. Techene) pure z* (quasi-eulerian) coordinate
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_update   : Swap vertical scale factors and update the vertical grid
18   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
19   !!   dom_vvl_rst      : read/write restart file
20   !!   dom_vvl_ctl      : Check the vvl options
21   !!----------------------------------------------------------------------
22   USE oce             ! ocean dynamics and tracers
23   USE phycst          ! physical constant
24   USE dom_oce         ! ocean space and time domain
25   USE sbc_oce         ! ocean surface boundary condition
26   USE wet_dry         ! wetting and drying
27   USE usrdef_istate   ! user defined initial state (wad only)
28   USE restart         ! ocean restart
29   !
30   USE in_out_manager  ! I/O manager
31   USE iom             ! I/O manager library
32   USE lib_mpp         ! distributed memory computing library
33   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
34   USE timing          ! Timing
35
36   IMPLICIT NONE
37   PRIVATE
38
39   PUBLIC  dom_qe_init       ! called by domain.F90
40   PUBLIC  dom_qe_zgr        ! called by isfcpl.F90
41   PUBLIC  dom_qe_sf_nxt     ! called by step.F90
42   PUBLIC  dom_qe_sf_update  ! called by step.F90
43   PUBLIC  dom_qe_interpol   ! called by dynnxt.F90
44
45   !                                                      !!* Namelist nam_vvl
46   LOGICAL , PUBLIC :: ln_vvl_zstar           = .FALSE.    ! zstar  vertical coordinate
47   LOGICAL , PUBLIC :: ln_vvl_ztilde          = .FALSE.    ! ztilde vertical coordinate
48   LOGICAL , PUBLIC :: ln_vvl_layer           = .FALSE.    ! level  vertical coordinate
49   LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE.    ! ztilde vertical coordinate
50   LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor  = .FALSE.    ! ztilde vertical coordinate
51   LOGICAL , PUBLIC :: ln_vvl_kepe            = .FALSE.    ! kinetic/potential energy transfer
52   !                                                       ! conservation: not used yet
53   REAL(wp)         :: rn_ahe3                             ! thickness diffusion coefficient
54   REAL(wp)         :: rn_rst_e3t                          ! ztilde to zstar restoration timescale [days]
55   REAL(wp)         :: rn_lf_cutoff                        ! cutoff frequency for low-pass filter  [days]
56   REAL(wp)         :: rn_zdef_max                         ! maximum fractional e3t deformation
57   LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE.                ! debug control prints
58
59   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td                ! thickness diffusion transport
60
61   !! * Substitutions
62#  include "do_loop_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
65   !! $Id: domvvl.F90 12377 2020-02-12 14:39:06Z acc $
66   !! Software governed by the CeCILL license (see ./LICENSE)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE dom_qe_init( Kbb, Kmm, Kaa )
71      !!----------------------------------------------------------------------
72      !!                ***  ROUTINE dom_qe_init  ***
73      !!
74      !! ** Purpose :  Initialization of all scale factors, depths
75      !!               and water column heights
76      !!
77      !! ** Method  :  - use restart file and/or initialize
78      !!               - interpolate scale factors
79      !!
80      !! ** Action  : - e3t_(n/b)
81      !!              - Regrid: e3[u/v](:,:,:,Kmm)
82      !!                        e3[u/v](:,:,:,Kmm)
83      !!                        e3w(:,:,:,Kmm)
84      !!                        e3[u/v]w_b
85      !!                        e3[u/v]w_n
86      !!                        gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm) and gde3w
87      !!              - h(t/u/v)_0
88      !!
89      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
90      !!----------------------------------------------------------------------
91      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
92      !
93      IF(lwp) WRITE(numout,*)
94      IF(lwp) WRITE(numout,*) 'dom_qe_init : Variable volume activated'
95      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
96      !
97      CALL dom_qe_ctl     ! choose vertical coordinate (z_star, z_tilde or layer)
98      !
99      !                    ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf
100      CALL dom_qe_rst( nit000, Kbb, Kmm, 'READ' )
101      e3t(:,:,jpk,Kaa) = e3t_0(:,:,jpk)  ! last level always inside the sea floor set one for all
102      !
103      CALL dom_qe_zgr(Kbb, Kmm, Kaa) ! interpolation scale factor, depth and water column
104      !
105      IF(lwxios) THEN   ! define variables in restart file when writing with XIOS
106         CALL iom_set_rstw_var_active('e3t_b')
107         CALL iom_set_rstw_var_active('e3t_n')
108      ENDIF
109      !
110   END SUBROUTINE dom_qe_init
111   !
112   SUBROUTINE dom_qe_zgr(Kbb, Kmm, Kaa)
113      !!----------------------------------------------------------------------
114      !!                ***  ROUTINE dom_qe_init  ***
115      !!
116      !! ** Purpose :  Interpolation of all scale factors,
117      !!               depths and water column heights
118      !!
119      !! ** Method  :  - interpolate scale factors
120      !!
121      !! ** Action  : - e3t_(n/b)
122      !!              - Regrid: e3(u/v)_n
123      !!                        e3(u/v)_b
124      !!                        e3w_n
125      !!                        e3(u/v)w_b
126      !!                        e3(u/v)w_n
127      !!                        gdept_n, gdepw_n and gde3w_n
128      !!              - h(t/u/v)_0
129      !!
130      !! Reference  : Leclair, M., and G. Madec, 2011, Ocean Modelling.
131      !!----------------------------------------------------------------------
132      INTEGER, INTENT(in) :: Kbb, Kmm, Kaa
133      !!----------------------------------------------------------------------
134      INTEGER ::   ji, jj, jk
135      INTEGER ::   ii0, ii1, ij0, ij1
136      REAL(wp)::   zcoef
137      !!----------------------------------------------------------------------
138      !
139      !                    !== Set of all other vertical scale factors  ==!  (now and before)
140      !                                ! Horizontal interpolation of e3t
141      CALL dom_qe_interpol( e3t(:,:,:,Kbb), e3u(:,:,:,Kbb), 'U' )    ! from T to U
142      CALL dom_qe_interpol( e3t(:,:,:,Kmm), e3u(:,:,:,Kmm), 'U' )
143      CALL dom_qe_interpol( e3t(:,:,:,Kbb), e3v(:,:,:,Kbb), 'V' )    ! from T to V
144      CALL dom_qe_interpol( e3t(:,:,:,Kmm), e3v(:,:,:,Kmm), 'V' )
145      CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' )    ! from U to F
146      !                                ! Vertical interpolation of e3t,u,v
147      CALL dom_qe_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W'  )  ! from T to W
148      CALL dom_qe_interpol( e3t(:,:,:,Kbb), e3w (:,:,:,Kbb), 'W'  )
149      CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )  ! from U to UW
150      CALL dom_qe_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
151      CALL dom_qe_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )  ! from V to UW
152      CALL dom_qe_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
153
154      ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...)
155      e3t(:,:,:,Kaa) = e3t(:,:,:,Kmm)
156      e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm)
157      e3v(:,:,:,Kaa) = e3v(:,:,:,Kmm)
158      !
159      !                    !==  depth of t and w-point  ==!   (set the isf depth as it is in the initial timestep)
160      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)       ! reference to the ocean surface (used for MLD and light penetration)
161      gdepw(:,:,1,Kmm) = 0.0_wp
162      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)  ! reference to a common level z=0 for hpg
163      gdept(:,:,1,Kbb) = 0.5_wp * e3w(:,:,1,Kbb)
164      gdepw(:,:,1,Kbb) = 0.0_wp
165      DO_3D_11_11( 2, jpk )
166         !    zcoef = tmask - wmask    ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
167         !                             ! 1 everywhere from mbkt to mikt + 1 or 1 (if no isf)
168         !                             ! 0.5 where jk = mikt
169!!gm ???????   BUG ?  gdept(:,:,:,Kmm) as well as gde3w  does not include the thickness of ISF ??
170         zcoef = ( tmask(ji,jj,jk) - wmask(ji,jj,jk) )
171         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
172         gdept(ji,jj,jk,Kmm) =      zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm))  &
173            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm))
174         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
175         gdepw(ji,jj,jk,Kbb) = gdepw(ji,jj,jk-1,Kbb) + e3t(ji,jj,jk-1,Kbb)
176         gdept(ji,jj,jk,Kbb) =      zcoef  * ( gdepw(ji,jj,jk  ,Kbb) + 0.5 * e3w(ji,jj,jk,Kbb))  &
177            &                + (1-zcoef) * ( gdept(ji,jj,jk-1,Kbb) +       e3w(ji,jj,jk,Kbb))
178      END_3D
179      !
180      !                    !==  thickness of the water column  !!   (ocean portion only)
181      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)   !!gm  BUG  :  this should be 1/2 * e3w(k=1) ....
182      hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1)
183      hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1)
184      hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1)
185      hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1)
186      DO jk = 2, jpkm1
187         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
188         hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk)
189         hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk)
190         hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)
191         hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)
192      END DO
193      !
194      !                    !==  inverse of water column thickness   ==!   (u- and v- points)
195      r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) )    ! _i mask due to ISF
196      r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) )
197      r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
198      r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) )
199      !
200   END SUBROUTINE dom_qe_zgr
201
202
203   SUBROUTINE dom_qe_sf_nxt( kt, Kbb, Kmm, Kaa, kcall )
204      !!----------------------------------------------------------------------
205      !!                ***  ROUTINE dom_qe_sf_nxt  ***
206      !!
207      !! ** Purpose :  - compute the after scale factors used in tra_zdf, dynnxt,
208      !!                 tranxt and dynspg routines
209      !!
210      !! ** Method  :  - z_star case:  Repartition of ssh INCREMENT proportionnaly to the level thickness.
211      !!
212      !! ** Action  :  - hdiv_lf    : restoring towards full baroclinic divergence in z_tilde case
213      !!               - tilde_e3t_a: after increment of vertical scale factor
214      !!                              in z_tilde case
215      !!               - e3(t/u/v)_a
216      !!
217      !! Reference  : Leclair, M., and Madec, G. 2011, Ocean Modelling.
218      !!----------------------------------------------------------------------
219      INTEGER, INTENT( in )           ::   kt             ! time step
220      INTEGER, INTENT( in )           ::   Kbb, Kmm, Kaa  ! time step
221      INTEGER, INTENT( in ), OPTIONAL ::   kcall          ! optional argument indicating call sequence
222      !
223      INTEGER                ::   ji, jj, jk            ! dummy loop indices
224      INTEGER , DIMENSION(3) ::   ijk_max, ijk_min      ! temporary integers
225      REAL(wp)               ::   z2dt, z_tmin, z_tmax  ! local scalars
226      LOGICAL                ::   ll_do_bclinic         ! local logical
227      REAL(wp), DIMENSION(jpi,jpj)     ::   zht, z_scale, zwu, zwv, zhdiv
228      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ze3t
229      !!----------------------------------------------------------------------
230      !
231      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
232      !
233      IF( ln_timing )   CALL timing_start('dom_qe_sf_nxt')
234      !
235      IF( kt == nit000 ) THEN
236         IF(lwp) WRITE(numout,*)
237         IF(lwp) WRITE(numout,*) 'dom_qe_sf_nxt : compute after scale factors'
238         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
239      ENDIF
240
241      ! ll_do_bclinic = .TRUE.
242      ! IF( PRESENT(kcall) ) THEN
243      !    IF( kcall == 2 .AND. ln_vvl_ztilde )   ll_do_bclinic = .FALSE.
244      ! ENDIF
245
246      ! ******************************* !
247      ! After acale factors at t-points !
248      ! ******************************* !
249      !                                                ! --------------------------------------------- !
250      !                                                ! z_star coordinate and barotropic z-tilde part !
251      !                                                ! --------------------------------------------- !
252      !
253      z_scale(:,:) = ( ssh(:,:,Kaa) - ssh(:,:,Kbb) ) * ssmask(:,:) / ( ht_0(:,:) + ssh(:,:,Kmm) + 1. - ssmask(:,:) )
254      DO jk = 1, jpkm1
255         ! formally this is the same as e3t(:,:,:,Kaa) = e3t_0*(1+ssha/ht_0)
256         e3t(:,:,jk,Kaa) = e3t(:,:,jk,Kbb) + e3t(:,:,jk,Kmm) * z_scale(:,:) * tmask(:,:,jk)
257      END DO
258      !
259      ! *********************************** !
260      ! After scale factors at u- v- points !
261      ! *********************************** !
262
263      CALL dom_qe_interpol( e3t(:,:,:,Kaa), e3u(:,:,:,Kaa), 'U' )
264      CALL dom_qe_interpol( e3t(:,:,:,Kaa), e3v(:,:,:,Kaa), 'V' )
265
266      ! *********************************** !
267      ! After depths at u- v points         !
268      ! *********************************** !
269
270      hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1)
271      hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1)
272      DO jk = 2, jpkm1
273         hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk)
274         hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)
275      END DO
276      !                                        ! Inverse of the local depth
277!!gm BUG ?  don't understand the use of umask_i here .....
278      r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) )
279      r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) )
280      !
281      IF( ln_timing )   CALL timing_stop('dom_qe_sf_nxt')
282      !
283   END SUBROUTINE dom_qe_sf_nxt
284
285
286   SUBROUTINE dom_qe_sf_update( kt, Kbb, Kmm, Kaa )
287      !!----------------------------------------------------------------------
288      !!                ***  ROUTINE dom_qe_sf_update  ***
289      !!
290      !! ** Purpose :  for z tilde case: compute time filter and swap of scale factors
291      !!               compute all depths and related variables for next time step
292      !!               write outputs and restart file
293      !!
294      !! ** Method  :  - reconstruct scale factor at other grid points (interpolate)
295      !!               - recompute depths and water height fields
296      !!
297      !! ** Action  :  - Recompute:
298      !!                    e3(u/v)_b
299      !!                    e3w(:,:,:,Kmm)
300      !!                    e3(u/v)w_b
301      !!                    e3(u/v)w_n
302      !!                    gdept(:,:,:,Kmm), gdepw(:,:,:,Kmm)  and gde3w
303      !!                    h(u/v) and h(u/v)r
304      !!
305      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
306      !!              Leclair, M., and G. Madec, 2011, Ocean Modelling.
307      !!----------------------------------------------------------------------
308      INTEGER, INTENT( in ) ::   kt              ! time step
309      INTEGER, INTENT( in ) ::   Kbb, Kmm, Kaa   ! time level indices
310      !
311      INTEGER  ::   ji, jj, jk   ! dummy loop indices
312      REAL(wp) ::   zcoef        ! local scalar
313      !!----------------------------------------------------------------------
314      !
315      IF( ln_linssh )   RETURN      ! No calculation in linear free surface
316      !
317      IF( ln_timing )   CALL timing_start('dom_qe_sf_update')
318      !
319      IF( kt == nit000 )   THEN
320         IF(lwp) WRITE(numout,*)
321         IF(lwp) WRITE(numout,*) 'dom_qe_sf_update : - interpolate scale factors and compute depths for next time step'
322         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~'
323      ENDIF
324      !
325      ! Compute all missing vertical scale factor and depths
326      ! ====================================================
327      ! Horizontal scale factor interpolations
328      ! --------------------------------------
329      ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt
330      ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also
331
332      CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F'  )
333
334      ! Vertical scale factor interpolations
335      CALL dom_qe_interpol( e3t(:,:,:,Kmm),  e3w(:,:,:,Kmm), 'W'  )
336      CALL dom_qe_interpol( e3u(:,:,:,Kmm), e3uw(:,:,:,Kmm), 'UW' )
337      CALL dom_qe_interpol( e3v(:,:,:,Kmm), e3vw(:,:,:,Kmm), 'VW' )
338      CALL dom_qe_interpol( e3t(:,:,:,Kbb),  e3w(:,:,:,Kbb), 'W'  )
339      CALL dom_qe_interpol( e3u(:,:,:,Kbb), e3uw(:,:,:,Kbb), 'UW' )
340      CALL dom_qe_interpol( e3v(:,:,:,Kbb), e3vw(:,:,:,Kbb), 'VW' )
341
342      ! t- and w- points depth (set the isf depth as it is in the initial step)
343      gdept(:,:,1,Kmm) = 0.5_wp * e3w(:,:,1,Kmm)
344      gdepw(:,:,1,Kmm) = 0.0_wp
345      gde3w(:,:,1) = gdept(:,:,1,Kmm) - ssh(:,:,Kmm)
346      DO_3D_11_11( 2, jpk )
347        !    zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))   ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt
348                                                           ! 1 for jk = mikt
349         zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk))
350         gdepw(ji,jj,jk,Kmm) = gdepw(ji,jj,jk-1,Kmm) + e3t(ji,jj,jk-1,Kmm)
351         gdept(ji,jj,jk,Kmm) =    zcoef  * ( gdepw(ji,jj,jk  ,Kmm) + 0.5 * e3w(ji,jj,jk,Kmm) )  &
352             &             + (1-zcoef) * ( gdept(ji,jj,jk-1,Kmm) +       e3w(ji,jj,jk,Kmm) )
353         gde3w(ji,jj,jk) = gdept(ji,jj,jk,Kmm) - ssh(ji,jj,Kmm)
354      END_3D
355
356      ! Local depth and Inverse of the local depth of the water
357      ! -------------------------------------------------------
358      !
359      ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1)
360      DO jk = 2, jpkm1
361         ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)
362      END DO
363
364      ! write restart file
365      ! ==================
366      IF( lrst_oce  )   CALL dom_qe_rst( kt, Kbb, Kmm, 'WRITE' )
367      !
368      IF( ln_timing )   CALL timing_stop('dom_qe_sf_update')
369      !
370   END SUBROUTINE dom_qe_sf_update
371
372
373   SUBROUTINE dom_qe_interpol( pe3_in, pe3_out, pout )
374      !!---------------------------------------------------------------------
375      !!                  ***  ROUTINE dom_qe_interpol  ***
376      !!
377      !! ** Purpose :   interpolate scale factors from one grid point to another
378      !!
379      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
380      !!                - horizontal interpolation: grid cell surface averaging
381      !!                - vertical interpolation: simple averaging
382      !!----------------------------------------------------------------------
383      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in   ) ::  pe3_in    ! input e3 to be interpolated
384      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  pe3_out   ! output interpolated e3
385      CHARACTER(LEN=*)                , INTENT(in   ) ::  pout      ! grid point of out scale factors
386      !                                                             !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
387      !
388      INTEGER ::   ji, jj, jk                                       ! dummy loop indices
389      REAL(wp) ::  zlnwd                                            ! =1./0. when ln_wd_il = T/F
390      !!----------------------------------------------------------------------
391      !
392      IF(ln_wd_il) THEN
393        zlnwd = 1.0_wp
394      ELSE
395        zlnwd = 0.0_wp
396      END IF
397      !
398      SELECT CASE ( pout )    !==  type of interpolation  ==!
399         !
400      CASE( 'U' )                   !* from T- to U-point : hor. surface weighted mean
401         DO_3D_10_10( 1, jpk )
402            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2u(ji,jj)   &
403               &                       * (   e1e2t(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0(ji  ,jj,jk) )     &
404               &                           + e1e2t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) )
405         END_3D
406         CALL lbc_lnk( 'domqe', pe3_out(:,:,:), 'U', 1._wp )
407         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:)
408         !
409      CASE( 'V' )                   !* from T- to V-point : hor. surface weighted mean
410         DO_3D_10_10( 1, jpk )
411            pe3_out(ji,jj,jk) = 0.5_wp * ( vmask(ji,jj,jk)  * (1.0_wp - zlnwd) + zlnwd ) * r1_e1e2v(ji,jj)   &
412               &                       * (   e1e2t(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
413               &                           + e1e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
414         END_3D
415         CALL lbc_lnk( 'domqe', pe3_out(:,:,:), 'V', 1._wp )
416         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:)
417         !
418      CASE( 'F' )                   !* from U-point to F-point : hor. surface weighted mean
419         DO_3D_10_10( 1, jpk )
420            pe3_out(ji,jj,jk) = 0.5_wp * (  umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) &
421               &                       *    r1_e1e2f(ji,jj)                                                  &
422               &                       * (   e1e2u(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3u_0(ji,jj  ,jk) )     &
423               &                           + e1e2u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) )
424         END_3D
425         CALL lbc_lnk( 'domqe', pe3_out(:,:,:), 'F', 1._wp )
426         pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:)
427         !
428      CASE( 'W' )                   !* from T- to W-point : vertical simple mean
429         !
430         pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1)
431         ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing
432!!gm BUG? use here wmask in case of ISF ?  to be checked
433         DO jk = 2, jpk
434            pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )   &
435               &                            * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) )                               &
436               &                            +            0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )     &
437               &                            * ( pe3_in(:,:,jk  ) - e3t_0(:,:,jk  ) )
438         END DO
439         !
440      CASE( 'UW' )                  !* from U- to UW-point : vertical simple mean
441         !
442         pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1)
443         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
444!!gm BUG? use here wumask in case of ISF ?  to be checked
445         DO jk = 2, jpk
446            pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
447               &                             * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) )                              &
448               &                             +            0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
449               &                             * ( pe3_in(:,:,jk  ) - e3u_0(:,:,jk  ) )
450         END DO
451         !
452      CASE( 'VW' )                  !* from V- to VW-point : vertical simple mean
453         !
454         pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1)
455         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
456!!gm BUG? use here wvmask in case of ISF ?  to be checked
457         DO jk = 2, jpk
458            pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) )  &
459               &                             * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) )                              &
460               &                             +            0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd )    &
461               &                             * ( pe3_in(:,:,jk  ) - e3v_0(:,:,jk  ) )
462         END DO
463      END SELECT
464      !
465   END SUBROUTINE dom_qe_interpol
466
467
468   SUBROUTINE dom_qe_rst( kt, Kbb, Kmm, cdrw )
469      !!---------------------------------------------------------------------
470      !!                   ***  ROUTINE dom_qe_rst  ***
471      !!
472      !! ** Purpose :   Read or write VVL file in restart file
473      !!
474      !! ** Method  :   use of IOM library
475      !!                if the restart does not contain vertical scale factors,
476      !!                they are set to the _0 values
477      !!                if the restart does not contain vertical scale factors increments (z_tilde),
478      !!                they are set to 0.
479      !!----------------------------------------------------------------------
480      INTEGER         , INTENT(in) ::   kt        ! ocean time-step
481      INTEGER         , INTENT(in) ::   Kbb, Kmm  ! ocean time level indices
482      CHARACTER(len=*), INTENT(in) ::   cdrw      ! "READ"/"WRITE" flag
483      !
484      INTEGER ::   ji, jj, jk
485      INTEGER ::   id1, id2     ! local integers
486      !!----------------------------------------------------------------------
487      !
488      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
489         !                                   ! ===============
490         IF( ln_rstart ) THEN                   !* Read the restart file
491            CALL rst_read_open                  !  open the restart file if necessary
492            CALL iom_get( numror, jpdom_autoglo, 'sshn'   , ssh(:,:,Kmm), ldxios = lrxios    )
493            !
494            id1 = iom_varid( numror, 'e3t_b', ldstop = .FALSE. )
495            id2 = iom_varid( numror, 'e3t_n', ldstop = .FALSE. )
496            !
497            !                             ! --------- !
498            !                             ! all cases !
499            !                             ! --------- !
500            !
501            IF( MIN( id1, id2 ) > 0 ) THEN       ! all required arrays exist
502               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
503               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
504               ! needed to restart if land processor not computed
505               IF(lwp) write(numout,*) 'dom_qe_rst : e3t(:,:,:,Kbb) and e3t(:,:,:,Kmm) found in restart files'
506               WHERE ( tmask(:,:,:) == 0.0_wp )
507                  e3t(:,:,:,Kmm) = e3t_0(:,:,:)
508                  e3t(:,:,:,Kbb) = e3t_0(:,:,:)
509               END WHERE
510               IF( neuler == 0 ) THEN
511                  e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
512               ENDIF
513            ELSE IF( id1 > 0 ) THEN
514               IF(lwp) write(numout,*) 'dom_qe_rst WARNING : e3t(:,:,:,Kmm) not found in restart files'
515               IF(lwp) write(numout,*) 'e3t_n set equal to e3t_b.'
516               IF(lwp) write(numout,*) 'neuler is forced to 0'
517               CALL iom_get( numror, jpdom_autoglo, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lrxios )
518               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
519               neuler = 0
520            ELSE IF( id2 > 0 ) THEN
521               IF(lwp) write(numout,*) 'dom_qe_rst WARNING : e3t(:,:,:,Kbb) not found in restart files'
522               IF(lwp) write(numout,*) 'e3t_b set equal to e3t_n.'
523               IF(lwp) write(numout,*) 'neuler is forced to 0'
524               CALL iom_get( numror, jpdom_autoglo, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lrxios )
525               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
526               neuler = 0
527            ELSE
528               IF(lwp) write(numout,*) 'dom_qe_rst WARNING : e3t(:,:,:,Kmm) not found in restart file'
529               IF(lwp) write(numout,*) 'Compute scale factor from sshn'
530               IF(lwp) write(numout,*) 'neuler is forced to 0'
531               DO jk = 1, jpk
532                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm) ) &
533                      &                          / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
534                      &          + e3t_0(:,:,jk)                               * (1._wp -tmask(:,:,jk))
535               END DO
536               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
537               neuler = 0
538            ENDIF
539            !
540         ELSE                                   !* Initialize at "rest"
541            !
542            IF( ll_wd ) THEN   ! MJB ll_wd edits start here - these are essential
543               !
544               IF( cn_cfg == 'wad' ) THEN
545                  ! Wetting and drying test case
546                  CALL usr_def_istate( gdept(:,:,:,Kbb), tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
547                  ts  (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb)       ! set now values from to before ones
548                  ssh (:,:,Kmm)     = ssh(:,:,Kbb)
549                  uu   (:,:,:,Kmm)   = uu  (:,:,:,Kbb)
550                  vv   (:,:,:,Kmm)   = vv  (:,:,:,Kbb)
551               ELSE
552                  ! if not test case
553                  ssh(:,:,Kmm) = -ssh_ref
554                  ssh(:,:,Kbb) = -ssh_ref
555
556                  DO_2D_11_11
557                     IF( ht_0(ji,jj)-ssh_ref <  rn_wdmin1 ) THEN ! if total depth is less than min depth
558                        ssh(ji,jj,Kbb) = rn_wdmin1 - (ht_0(ji,jj) )
559                        ssh(ji,jj,Kmm) = rn_wdmin1 - (ht_0(ji,jj) )
560                     ENDIF
561                  END_2D
562               ENDIF !If test case else
563
564               ! Adjust vertical metrics for all wad
565               DO jk = 1, jpk
566                  e3t(:,:,jk,Kmm) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kmm)  ) &
567                    &                            / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk)   &
568                    &            + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) )
569               END DO
570               e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)
571
572               DO ji = 1, jpi
573                  DO jj = 1, jpj
574                     IF ( ht_0(ji,jj) .LE. 0.0 .AND. NINT( ssmask(ji,jj) ) .EQ. 1) THEN
575                       CALL ctl_stop( 'dom_qe_rst: ht_0 must be positive at potentially wet points' )
576                     ENDIF
577                  END DO
578               END DO
579               !
580            ELSE
581               !
582               ! Just to read set ssh in fact, called latter once vertical grid
583               ! is set up:
584!               CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb)  )
585!               !
586!               DO jk=1,jpk
587!                  e3t(:,:,jk,Kbb) =  e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) &
588!                     &            / ( ht_0(:,:) + 1._wp -ssmask(:,:) ) * tmask(:,:,jk)
589!               END DO
590!               e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb)
591                ssh(:,:,Kmm)=0._wp
592                e3t(:,:,:,Kmm)=e3t_0(:,:,:)
593                e3t(:,:,:,Kbb)=e3t_0(:,:,:)
594               !
595            ENDIF           ! end of ll_wd edits
596            !
597         ENDIF
598         !
599      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
600         !                                   ! ===================
601         IF(lwp) WRITE(numout,*) '---- dom_qe_rst ----'
602         IF( lwxios ) CALL iom_swap(      cwxios_context          )
603         !                                           ! --------- !
604         !                                           ! all cases !
605         !                                           ! --------- !
606         CALL iom_rstput( kt, nitrst, numrow, 'e3t_b', e3t(:,:,:,Kbb), ldxios = lwxios )
607         CALL iom_rstput( kt, nitrst, numrow, 'e3t_n', e3t(:,:,:,Kmm), ldxios = lwxios )
608         !
609         IF( lwxios ) CALL iom_swap(      cxios_context          )
610      ENDIF
611      !
612   END SUBROUTINE dom_qe_rst
613
614
615   SUBROUTINE dom_qe_ctl
616      !!---------------------------------------------------------------------
617      !!                  ***  ROUTINE dom_qe_ctl  ***
618      !!
619      !! ** Purpose :   Control the consistency between namelist options
620      !!                for vertical coordinate
621      !!----------------------------------------------------------------------
622      INTEGER ::   ioptio, ios
623      !!
624      NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, &
625         &              ln_vvl_zstar_at_eqtor      , rn_ahe3     , rn_rst_e3t            , &
626         &              rn_lf_cutoff               , rn_zdef_max , ln_vvl_dbg                ! not yet implemented: ln_vvl_kepe
627      !!----------------------------------------------------------------------
628      !
629      READ  ( numnam_ref, nam_vvl, IOSTAT = ios, ERR = 901)
630901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_vvl in reference namelist' )
631      READ  ( numnam_cfg, nam_vvl, IOSTAT = ios, ERR = 902 )
632902   IF( ios >  0 ) CALL ctl_nam ( ios , 'nam_vvl in configuration namelist' )
633      IF(lwm) WRITE ( numond, nam_vvl )
634      !
635      IF(lwp) THEN                    ! Namelist print
636         WRITE(numout,*)
637         WRITE(numout,*) 'dom_qe_ctl : choice/control of the variable vertical coordinate'
638         WRITE(numout,*) '~~~~~~~~~~~'
639         WRITE(numout,*) '   Namelist nam_vvl : chose a vertical coordinate'
640         WRITE(numout,*) '      zstar                      ln_vvl_zstar   = ', ln_vvl_zstar
641         WRITE(numout,*) '      ztilde                     ln_vvl_ztilde  = ', ln_vvl_ztilde
642         WRITE(numout,*) '      layer                      ln_vvl_layer   = ', ln_vvl_layer
643         WRITE(numout,*) '      ztilde as zstar   ln_vvl_ztilde_as_zstar  = ', ln_vvl_ztilde_as_zstar
644         WRITE(numout,*) '      ztilde near the equator    ln_vvl_zstar_at_eqtor  = ', ln_vvl_zstar_at_eqtor
645         WRITE(numout,*) '      !'
646         WRITE(numout,*) '      thickness diffusion coefficient                      rn_ahe3      = ', rn_ahe3
647         WRITE(numout,*) '      maximum e3t deformation fractional change            rn_zdef_max  = ', rn_zdef_max
648         IF( ln_vvl_ztilde_as_zstar ) THEN
649            WRITE(numout,*) '      ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T) '
650            WRITE(numout,*) '         ignoring namelist timescale parameters and using:'
651            WRITE(numout,*) '            hard-wired : z-tilde to zstar restoration timescale (days)'
652            WRITE(numout,*) '                         rn_rst_e3t     = 0.e0'
653            WRITE(numout,*) '            hard-wired : z-tilde cutoff frequency of low-pass filter (days)'
654            WRITE(numout,*) '                         rn_lf_cutoff   = 1.0/rdt'
655         ELSE
656            WRITE(numout,*) '      z-tilde to zstar restoration timescale (days)        rn_rst_e3t   = ', rn_rst_e3t
657            WRITE(numout,*) '      z-tilde cutoff frequency of low-pass filter (days)   rn_lf_cutoff = ', rn_lf_cutoff
658         ENDIF
659         WRITE(numout,*) '         debug prints flag                                 ln_vvl_dbg   = ', ln_vvl_dbg
660      ENDIF
661      !
662      ioptio = 0                      ! Parameter control
663      IF( ln_vvl_ztilde_as_zstar )   ln_vvl_ztilde = .true.
664      IF( ln_vvl_zstar           )   ioptio = ioptio + 1
665      IF( ln_vvl_ztilde          )   ioptio = ioptio + 1
666      IF( ln_vvl_layer           )   ioptio = ioptio + 1
667      !
668      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' )
669      !
670      IF(lwp) THEN                   ! Print the choice
671         WRITE(numout,*)
672         IF( ln_vvl_zstar           ) WRITE(numout,*) '      ==>>>   zstar vertical coordinate is used'
673         IF( ln_vvl_ztilde          ) WRITE(numout,*) '      ==>>>   ztilde vertical coordinate is used'
674         IF( ln_vvl_layer           ) WRITE(numout,*) '      ==>>>   layer vertical coordinate is used'
675         IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) '      ==>>>   to emulate a zstar coordinate'
676      ENDIF
677      !
678#if defined key_agrif
679      IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) )   CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' )
680#endif
681      !
682   END SUBROUTINE dom_qe_ctl
683
684   !!======================================================================
685END MODULE domqe
Note: See TracBrowser for help on using the repository browser.