source: NEMO/branches/2020/dev_r12377_KERNEL-06_techene_e3/src/OCE/DYN/sshwzv.F90 @ 12622

Last change on this file since 12622 was 12622, checked in by techene, 8 months ago

all: add e3 substitute (sometimes it requires to add ze3t/u/v/w) and limit precompiled files lines to about 130 character

  • Property svn:keywords set to Id
File size: 19.2 KB
Line 
1MODULE sshwzv
2   !!==============================================================================
3   !!                       ***  MODULE  sshwzv  ***
4   !! Ocean dynamics : sea surface height and vertical velocity
5   !!==============================================================================
6   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
7   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
8   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
9   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
10   !!            3.3  !  2011-10  (M. Leclair) split former ssh_wzv routine and remove all vvl related work
11   !!            4.0  !  2018-12  (A. Coward) add mixed implicit/explicit advection
12   !!            4.1  !  2019-08  (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering.
13   !!----------------------------------------------------------------------
14
15   !!----------------------------------------------------------------------
16   !!   ssh_nxt       : after ssh
17   !!   ssh_atf       : time filter the ssh arrays
18   !!   wzv           : compute now vertical velocity
19   !!----------------------------------------------------------------------
20   USE oce            ! ocean dynamics and tracers variables
21   USE isf_oce        ! ice shelf
22   USE dom_oce        ! ocean space and time domain variables
23   USE sbc_oce        ! surface boundary condition: ocean
24   USE domvvl         ! Variable volume
25   USE divhor         ! horizontal divergence
26   USE phycst         ! physical constants
27   USE bdy_oce , ONLY : ln_bdy, bdytmask   ! Open BounDarY
28   USE bdydyn2d       ! bdy_ssh routine
29#if defined key_agrif
30   USE agrif_oce_interp
31#endif
32   !
33   USE iom
34   USE in_out_manager ! I/O manager
35   USE restart        ! only for lrst_oce
36   USE prtctl         ! Print control
37   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
38   USE lib_mpp        ! MPP library
39   USE timing         ! Timing
40   USE wet_dry        ! Wetting/Drying flux limiting
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   ssh_nxt    ! called by step.F90
46   PUBLIC   wzv        ! called by step.F90
47   PUBLIC   wAimp      ! called by step.F90
48   PUBLIC   ssh_atf    ! called by step.F90
49
50   !! * Substitutions
51#  include "do_loop_substitute.h90"
52#  include "domzgr_substitute.h90"
53
54   !!----------------------------------------------------------------------
55   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
56   !! $Id$
57   !! Software governed by the CeCILL license (see ./LICENSE)
58   !!----------------------------------------------------------------------
59CONTAINS
60
61   SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
62      !!----------------------------------------------------------------------
63      !!                ***  ROUTINE ssh_nxt  ***
64      !!
65      !! ** Purpose :   compute the after ssh (ssh(Kaa))
66      !!
67      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
68      !!      is computed by integrating the horizontal divergence and multiply by
69      !!      by the time step.
70      !!
71      !! ** action  :   ssh(:,:,Kaa), after sea surface height
72      !!
73      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
74      !!----------------------------------------------------------------------
75      INTEGER                         , INTENT(in   ) ::   kt             ! time step
76      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index
77      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height
78      !
79      INTEGER  ::   jk            ! dummy loop indice
80      REAL(wp) ::   z2dt, zcoef   ! local scalars
81      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
82      !!----------------------------------------------------------------------
83      !
84      IF( ln_timing )   CALL timing_start('ssh_nxt')
85      !
86      IF( kt == nit000 ) THEN
87         IF(lwp) WRITE(numout,*)
88         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
89         IF(lwp) WRITE(numout,*) '~~~~~~~ '
90      ENDIF
91      !
92      z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog)
93      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
94      zcoef = 0.5_wp * r1_rau0
95
96      !                                           !------------------------------!
97      !                                           !   After Sea Surface Height   !
98      !                                           !------------------------------!
99      IF(ln_wd_il) THEN
100         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt, Kmm, uu, vv )
101      ENDIF
102
103      CALL div_hor( kt, Kbb, Kmm )                     ! Horizontal divergence
104      !
105      zhdiv(:,:) = 0._wp
106      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
107        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
108      END DO
109      !                                                ! Sea surface elevation time stepping
110      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
111      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
112      !
113      pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
114      !
115#if defined key_agrif
116      Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt )
117#endif
118      !
119      IF ( .NOT.ln_dynspg_ts ) THEN
120         IF( ln_bdy ) THEN
121            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. )    ! Not sure that's necessary
122            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries
123         ENDIF
124      ENDIF
125      !                                           !------------------------------!
126      !                                           !           outputs            !
127      !                                           !------------------------------!
128      !
129      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask )
130      !
131      IF( ln_timing )   CALL timing_stop('ssh_nxt')
132      !
133   END SUBROUTINE ssh_nxt
134
135
136   SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa )
137      !!----------------------------------------------------------------------
138      !!                ***  ROUTINE wzv  ***
139      !!
140      !! ** Purpose :   compute the now vertical velocity
141      !!
142      !! ** Method  : - Using the incompressibility hypothesis, the vertical
143      !!      velocity is computed by integrating the horizontal divergence
144      !!      from the bottom to the surface minus the scale factor evolution.
145      !!        The boundary conditions are w=0 at the bottom (no flux) and.
146      !!
147      !! ** action  :   pww      : now vertical velocity
148      !!
149      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
150      !!----------------------------------------------------------------------
151      INTEGER                         , INTENT(in)    ::   kt             ! time step
152      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices
153      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! now vertical velocity
154      !
155      INTEGER  ::   ji, jj, jk   ! dummy loop indices
156      REAL(wp) ::   z1_2dt       ! local scalars
157      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
158      !!----------------------------------------------------------------------
159      !
160      IF( ln_timing )   CALL timing_start('wzv')
161      !
162      IF( kt == nit000 ) THEN
163         IF(lwp) WRITE(numout,*)
164         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
165         IF(lwp) WRITE(numout,*) '~~~~~ '
166         !
167         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
168      ENDIF
169      !                                           !------------------------------!
170      !                                           !     Now Vertical Velocity    !
171      !                                           !------------------------------!
172      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
173      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
174      !
175      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
176         ALLOCATE( zhdiv(jpi,jpj,jpk) )
177         !
178         DO jk = 1, jpkm1
179            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
180            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
181            DO_2D_00_00
182               zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
183            END_2D
184         END DO
185         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
186         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
187         !                             ! Same question holds for hdiv. Perhaps just for security
188         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
189            ! computation of w
190            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)    &
191               &                         + zhdiv(:,:,jk)                       &
192               &                         + z1_2dt * ( e3t(:,:,jk,Kaa)          &
193               &                                    - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
194         END DO
195         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
196         DEALLOCATE( zhdiv )
197      ELSE   ! z_star and linear free surface cases
198         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
199            ! computation of w
200            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)    &
201               &                         + z1_2dt * ( e3t(:,:,jk,Kaa)          &
202               &                                    - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
203         END DO
204      ENDIF
205
206      IF( ln_bdy ) THEN
207         DO jk = 1, jpkm1
208            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
209         END DO
210      ENDIF
211      !
212#if defined key_agrif
213      IF( .NOT. AGRIF_Root() ) THEN
214         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east
215         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west
216         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north
217         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south
218      ENDIF
219#endif
220      !
221      IF( ln_timing )   CALL timing_stop('wzv')
222      !
223   END SUBROUTINE wzv
224
225
226   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
227      !!----------------------------------------------------------------------
228      !!                    ***  ROUTINE ssh_atf  ***
229      !!
230      !! ** Purpose :   Apply Asselin time filter to now SSH.
231      !!
232      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
233      !!              from the filter, see Leclair and Madec 2010) and swap :
234      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa) + atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
235      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
236      !!
237      !! ** action  : - pssh(:,:,Kmm) time filtered
238      !!
239      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
240      !!----------------------------------------------------------------------
241      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
242      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices
243      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! SSH field
244      !
245      REAL(wp) ::   zcoef   ! local scalar
246      !!----------------------------------------------------------------------
247      !
248      IF( ln_timing )   CALL timing_start('ssh_atf')
249      !
250      IF( kt == nit000 ) THEN
251         IF(lwp) WRITE(numout,*)
252         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
253         IF(lwp) WRITE(numout,*) '~~~~~~~ '
254      ENDIF
255      !              !==  Euler time-stepping: no filter, just swap  ==!
256      IF ( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN   ! Only do time filtering for leapfrog timesteps
257         !                                                  ! filtered "now" field
258         pssh(:,:,Kmm) = pssh(:,:,Kmm) + atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
259         IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed
260            zcoef = atfp * rdt * r1_rau0
261            pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
262               &                             - rnf_b(:,:)        + rnf   (:,:)       &
263               &                             + fwfisf_cav_b(:,:) - fwfisf_cav(:,:)   &
264               &                             + fwfisf_par_b(:,:) - fwfisf_par(:,:)   ) * ssmask(:,:)
265
266            ! ice sheet coupling
267            IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - atfp * rdt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
268
269         ENDIF
270      ENDIF
271      !
272      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm)  - : ', mask1=tmask )
273      !
274      IF( ln_timing )   CALL timing_stop('ssh_atf')
275      !
276   END SUBROUTINE ssh_atf
277
278   SUBROUTINE wAimp( kt, Kmm )
279      !!----------------------------------------------------------------------
280      !!                ***  ROUTINE wAimp  ***
281      !!
282      !! ** Purpose :   compute the Courant number and partition vertical velocity
283      !!                if a proportion needs to be treated implicitly
284      !!
285      !! ** Method  : -
286      !!
287      !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
288      !!            :   wi      : now vertical velocity (for implicit treatment)
289      !!
290      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
291      !!              implicit scheme for vertical advection in oceanic modeling.
292      !!              Ocean Modelling, 91, 38-69.
293      !!----------------------------------------------------------------------
294      INTEGER, INTENT(in) ::   kt   ! time step
295      INTEGER, INTENT(in) ::   Kmm  ! time level index
296      !
297      INTEGER  ::   ji, jj, jk   ! dummy loop indices
298      REAL(wp)             ::   zCu, zcff, z1_e3t                     ! local scalars
299      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
300      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters
301      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
302      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
303      !!----------------------------------------------------------------------
304      !
305      IF( ln_timing )   CALL timing_start('wAimp')
306      !
307      IF( kt == nit000 ) THEN
308         IF(lwp) WRITE(numout,*)
309         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
310         IF(lwp) WRITE(numout,*) '~~~~~ '
311         wi(:,:,:) = 0._wp
312      ENDIF
313      !
314      ! Calculate Courant numbers
315      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
316         DO_3D_00_00( 1, jpkm1 )
317            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
318            ! 2*rdt and not r2dt (for restartability)
319            Cu_adv(ji,jj,jk) =  2._wp * rdt *   &
320               &  ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )            &
321               &  + ( MAX( e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)                                  &
322               &                        * uu (ji  ,jj,jk,Kmm) + un_td(ji  ,jj,jk), 0._wp ) -   &
323               &      MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)                                  &
324               &                        * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) )   &
325               &    * r1_e1e2t(ji  ,jj)                                                        &
326               &  + ( MAX( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)                                  &
327               &                        * vv (ji,jj  ,jk,Kmm) + vn_td(ji,jj  ,jk), 0._wp ) -   &
328               &      MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)                                  &
329               &                        * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) )   &
330               &    * r1_e1e2t(ji,jj  )                                                        &
331               &  ) * z1_e3t
332         END_3D
333      ELSE
334         DO_3D_00_00( 1, jpkm1 )
335            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
336            ! 2*rdt and not r2dt (for restartability)
337            Cu_adv(ji,jj,jk) =   2._wp * rdt *   &
338               &  ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )         &
339               &  + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
340               &      MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
341               &    * r1_e1e2t(ji,jj)                                                       &
342               &  + ( MAX( e1v(ji,jj  )*e3v(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   &
343               &      MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   &
344               &    * r1_e1e2t(ji,jj)                                                       &
345               &  ) * z1_e3t
346         END_3D
347      ENDIF
348      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )
349      !
350      CALL iom_put("Courant",Cu_adv)
351      !
352      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
353         DO_3DS_11_11( jpkm1, 2, -1 )
354            !
355            zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
356! alt:
357!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN
358!                     zCu =  Cu_adv(ji,jj,jk)
359!                  ELSE
360!                     zCu =  Cu_adv(ji,jj,jk-1)
361!                  ENDIF
362            !
363            IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
364               zcff = 0._wp
365            ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
366               zcff = ( zCu - Cu_min )**2
367               zcff = zcff / ( Fcu + zcff )
368            ELSE                                  !<-- Mostly implicit
369               zcff = ( zCu - Cu_max )/ zCu
370            ENDIF
371            zcff = MIN(1._wp, zcff)
372            !
373            wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
374            ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
375            !
376            Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl
377         END_3D
378         Cu_adv(:,:,1) = 0._wp
379      ELSE
380         ! Fully explicit everywhere
381         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient below and in stp_ctl
382         wi    (:,:,:) = 0._wp
383      ENDIF
384      CALL iom_put("wimp",wi)
385      CALL iom_put("wi_cff",Cu_adv)
386      CALL iom_put("wexp",ww)
387      !
388      IF( ln_timing )   CALL timing_stop('wAimp')
389      !
390   END SUBROUTINE wAimp
391   !!======================================================================
392END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.