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_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN – NEMO

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/sshwzv.F90 @ 13247

Last change on this file since 13247 was 13247, checked in by francesca, 4 years ago

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