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

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

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

dynatfQCO.F90, stepLF.F90 : fixed (remove pe3. from dyn_atf_qco input arguments), all : remove e3. tables and include gurvan's feedbacks
  • Property svn:keywords set to Id
File size: 19.9 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
117      CALL agrif_ssh( kt )
118#endif
119      !
120      IF ( .NOT.ln_dynspg_ts ) THEN
121         IF( ln_bdy ) THEN
122            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. )    ! Not sure that's necessary
123            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries
124         ENDIF
125      ENDIF
126      !                                           !------------------------------!
127      !                                           !           outputs            !
128      !                                           !------------------------------!
129      !
130      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask )
131      !
132      IF( ln_timing )   CALL timing_stop('ssh_nxt')
133      !
134   END SUBROUTINE ssh_nxt
135
136
137   SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
138      !!----------------------------------------------------------------------
139      !!                ***  ROUTINE wzv  ***
140      !!
141      !! ** Purpose :   compute the now vertical velocity
142      !!
143      !! ** Method  : - Using the incompressibility hypothesis, the vertical
144      !!      velocity is computed by integrating the horizontal divergence
145      !!      from the bottom to the surface minus the scale factor evolution.
146      !!        The boundary conditions are w=0 at the bottom (no flux) and.
147      !!
148      !! ** action  :   pww      : now vertical velocity
149      !!
150      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
151      !!----------------------------------------------------------------------
152      INTEGER                         , INTENT(in)    ::   kt             ! time step
153      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices
154      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! vertical velocity at Kmm
155      !
156      INTEGER  ::   ji, jj, jk   ! dummy loop indices
157      REAL(wp) ::   z1_2dt       ! local scalars
158      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
159      !!----------------------------------------------------------------------
160      !
161      IF( ln_timing )   CALL timing_start('wzv')
162      !
163      IF( kt == nit000 ) THEN
164         IF(lwp) WRITE(numout,*)
165         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
166         IF(lwp) WRITE(numout,*) '~~~~~ '
167         !
168         pww(:,:,jpk) = 0._wp           ! bottom boundary condition: w=0 (set once for all)
169      ENDIF
170      !
171      z1_2dt = 1. / ( 2. * rdt )        ! set time step size (Euler/Leapfrog)
172      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
173      !
174      !                                               !===============================!
175      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      !==  z_tilde and layer cases  ==!
176         !                                            !===============================!
177         ALLOCATE( zhdiv(jpi,jpj,jpk) )
178         !
179         DO jk = 1, jpkm1
180            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
181            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
182            DO_2D_00_00
183               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) )
184            END_2D
185         END DO
186         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
187         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
188         !                             ! Same question holds for hdiv. Perhaps just for security
189         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
190            ! computation of w
191            pww(:,:,jk) = pww(:,:,jk+1) - (   e3t(:,:,jk,Kmm) * hdiv(:,:,jk)   &
192               &                            +                  zhdiv(:,:,jk)   &
193               &                            + z1_2dt * ( e3t(:,:,jk,Kaa)       &
194               &                                       - e3t(:,:,jk,Kbb) )   ) * tmask(:,:,jk)
195         END DO
196         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
197         DEALLOCATE( zhdiv )
198         !                                            !=================================!
199      ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==!
200         !                                            !=================================!
201         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
202            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)  ) * tmask(:,:,jk)
203         END DO
204         !                                            !==========================================!
205      ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
206         !                                            !==========================================!
207         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
208            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)    &
209               &                         + z1_2dt * ( e3t(:,:,jk,Kaa)          &
210               &                                    - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
211         END DO
212      ENDIF
213
214      IF( ln_bdy ) THEN
215         DO jk = 1, jpkm1
216            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
217         END DO
218      ENDIF
219      !
220#if defined key_agrif
221      IF( .NOT. AGRIF_Root() ) THEN
222         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east
223         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west
224         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north
225         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south
226      ENDIF
227#endif
228      !
229      IF( ln_timing )   CALL timing_stop('wzv')
230      !
231   END SUBROUTINE wzv
232
233
234   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
235      !!----------------------------------------------------------------------
236      !!                    ***  ROUTINE ssh_atf  ***
237      !!
238      !! ** Purpose :   Apply Asselin time filter to now SSH.
239      !!
240      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
241      !!              from the filter, see Leclair and Madec 2010) and swap :
242      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa) + atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
243      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
244      !!
245      !! ** action  : - pssh(:,:,Kmm) time filtered
246      !!
247      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
248      !!----------------------------------------------------------------------
249      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index
250      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! ocean time level indices
251      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! SSH field
252      !
253      REAL(wp) ::   zcoef   ! local scalar
254      !!----------------------------------------------------------------------
255      !
256      IF( ln_timing )   CALL timing_start('ssh_atf')
257      !
258      IF( kt == nit000 ) THEN
259         IF(lwp) WRITE(numout,*)
260         IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height'
261         IF(lwp) WRITE(numout,*) '~~~~~~~ '
262      ENDIF
263      !              !==  Euler time-stepping: no filter, just swap  ==!
264      IF ( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN   ! Only do time filtering for leapfrog timesteps
265         !                                                  ! filtered "now" field
266         pssh(:,:,Kmm) = pssh(:,:,Kmm) + atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) )
267         IF( .NOT.ln_linssh ) THEN                          ! "now" <-- with forcing removed
268            zcoef = atfp * rdt * r1_rau0
269            pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
270               &                             - rnf_b(:,:)        + rnf   (:,:)       &
271               &                             + fwfisf_cav_b(:,:) - fwfisf_cav(:,:)   &
272               &                             + fwfisf_par_b(:,:) - fwfisf_par(:,:)   ) * ssmask(:,:)
273
274            ! ice sheet coupling
275            IF ( ln_isf .AND. ln_isfcpl .AND. kt == nit000+1) pssh(:,:,Kbb) = pssh(:,:,Kbb) - atfp * rdt * ( risfcpl_ssh(:,:) - 0.0 ) * ssmask(:,:)
276
277         ENDIF
278      ENDIF
279      !
280      IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm)  - : ', mask1=tmask )
281      !
282      IF( ln_timing )   CALL timing_stop('ssh_atf')
283      !
284   END SUBROUTINE ssh_atf
285
286   SUBROUTINE wAimp( kt, Kmm )
287      !!----------------------------------------------------------------------
288      !!                ***  ROUTINE wAimp  ***
289      !!
290      !! ** Purpose :   compute the Courant number and partition vertical velocity
291      !!                if a proportion needs to be treated implicitly
292      !!
293      !! ** Method  : -
294      !!
295      !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
296      !!            :   wi      : now vertical velocity (for implicit treatment)
297      !!
298      !! Reference  : Shchepetkin, A. F. (2015): An adaptive, Courant-number-dependent
299      !!              implicit scheme for vertical advection in oceanic modeling.
300      !!              Ocean Modelling, 91, 38-69.
301      !!----------------------------------------------------------------------
302      INTEGER, INTENT(in) ::   kt   ! time step
303      INTEGER, INTENT(in) ::   Kmm  ! time level index
304      !
305      INTEGER  ::   ji, jj, jk   ! dummy loop indices
306      REAL(wp)             ::   zCu, zcff, z1_e3t                     ! local scalars
307      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
308      REAL(wp) , PARAMETER ::   Cu_max = 0.30_wp                      ! local parameters
309      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
310      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
311      !!----------------------------------------------------------------------
312      !
313      IF( ln_timing )   CALL timing_start('wAimp')
314      !
315      IF( kt == nit000 ) THEN
316         IF(lwp) WRITE(numout,*)
317         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
318         IF(lwp) WRITE(numout,*) '~~~~~ '
319         wi(:,:,:) = 0._wp
320      ENDIF
321      !
322      ! Calculate Courant numbers
323      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN
324         DO_3D_00_00( 1, jpkm1 )
325            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
326            ! 2*rdt and not r2dt (for restartability)
327            Cu_adv(ji,jj,jk) =  2._wp * rdt *   &
328               &  ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )            &
329               &  + ( MAX( e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm)                                  &
330               &                        * uu (ji  ,jj,jk,Kmm) + un_td(ji  ,jj,jk), 0._wp ) -   &
331               &      MIN( e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm)                                  &
332               &                        * uu (ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) )   &
333               &    * r1_e1e2t(ji  ,jj)                                                        &
334               &  + ( MAX( e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm)                                  &
335               &                        * vv (ji,jj  ,jk,Kmm) + vn_td(ji,jj  ,jk), 0._wp ) -   &
336               &      MIN( e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm)                                  &
337               &                        * vv (ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) )   &
338               &    * r1_e1e2t(ji,jj  )                                                        &
339               &  ) * z1_e3t
340         END_3D
341      ELSE
342         DO_3D_00_00( 1, jpkm1 )
343            z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
344            ! 2*rdt and not r2dt (for restartability)
345            Cu_adv(ji,jj,jk) =   2._wp * rdt *   &
346               &  ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )         &
347               &  + ( MAX( e2u(ji  ,jj)*e3u(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
348               &      MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
349               &    * r1_e1e2t(ji,jj)                                                       &
350               &  + ( MAX( e1v(ji,jj  )*e3v(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   &
351               &      MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   &
352               &    * r1_e1e2t(ji,jj)                                                       &
353               &  ) * z1_e3t
354         END_3D
355      ENDIF
356      CALL lbc_lnk( 'sshwzv', Cu_adv, 'T', 1. )
357      !
358      CALL iom_put("Courant",Cu_adv)
359      !
360      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
361         DO_3DS_11_11( jpkm1, 2, -1 )
362            !
363            zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk-1) )
364! alt:
365!                  IF ( ww(ji,jj,jk) > 0._wp ) THEN
366!                     zCu =  Cu_adv(ji,jj,jk)
367!                  ELSE
368!                     zCu =  Cu_adv(ji,jj,jk-1)
369!                  ENDIF
370            !
371            IF( zCu <= Cu_min ) THEN              !<-- Fully explicit
372               zcff = 0._wp
373            ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
374               zcff = ( zCu - Cu_min )**2
375               zcff = zcff / ( Fcu + zcff )
376            ELSE                                  !<-- Mostly implicit
377               zcff = ( zCu - Cu_max )/ zCu
378            ENDIF
379            zcff = MIN(1._wp, zcff)
380            !
381            wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
382            ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
383            !
384            Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient below and in stp_ctl
385         END_3D
386         Cu_adv(:,:,1) = 0._wp
387      ELSE
388         ! Fully explicit everywhere
389         Cu_adv(:,:,:) = 0._wp                          ! Reuse array to output coefficient below and in stp_ctl
390         wi    (:,:,:) = 0._wp
391      ENDIF
392      CALL iom_put("wimp",wi)
393      CALL iom_put("wi_cff",Cu_adv)
394      CALL iom_put("wexp",ww)
395      !
396      IF( ln_timing )   CALL timing_stop('wAimp')
397      !
398   END SUBROUTINE wAimp
399   !!======================================================================
400END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.