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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/DYN/sshwzv.F90 @ 11466

Last change on this file since 11466 was 11466, checked in by davestorkey, 5 years ago

dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap: Bug fixes for ssh_atf routine in sshwzv.
  • Property svn:keywords set to Id
File size: 17.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   !!----------------------------------------------------------------------
12
13   !!----------------------------------------------------------------------
14   !!   ssh_nxt       : after ssh
15   !!   ssh_atf       : filter ans swap the ssh arrays
16   !!   wzv           : compute now vertical velocity
17   !!----------------------------------------------------------------------
18   USE oce            ! ocean dynamics and tracers variables
19   USE dom_oce        ! ocean space and time domain variables
20   USE sbc_oce        ! surface boundary condition: ocean
21   USE domvvl         ! Variable volume
22   USE divhor         ! horizontal divergence
23   USE phycst         ! physical constants
24   USE bdy_oce , ONLY : ln_bdy, bdytmask   ! Open BounDarY
25   USE bdydyn2d       ! bdy_ssh routine
26#if defined key_agrif
27   USE agrif_oce_interp
28#endif
29   !
30   USE iom 
31   USE in_out_manager ! I/O manager
32   USE restart        ! only for lrst_oce
33   USE prtctl         ! Print control
34   USE lbclnk         ! ocean lateral boundary condition (or mpp link)
35   USE lib_mpp        ! MPP library
36   USE timing         ! Timing
37   USE wet_dry        ! Wetting/Drying flux limiting
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   ssh_nxt    ! called by step.F90
43   PUBLIC   wzv        ! called by step.F90
44   PUBLIC   wAimp      ! called by step.F90
45   PUBLIC   ssh_atf    ! called by step.F90
46
47   !! * Substitutions
48#  include "vectopt_loop_substitute.h90"
49   !!----------------------------------------------------------------------
50   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
51   !! $Id$
52   !! Software governed by the CeCILL license (see ./LICENSE)
53   !!----------------------------------------------------------------------
54CONTAINS
55
56   SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa )
57      !!----------------------------------------------------------------------
58      !!                ***  ROUTINE ssh_nxt  ***
59      !!                   
60      !! ** Purpose :   compute the after ssh (ssh(Kaa))
61      !!
62      !! ** Method  : - Using the incompressibility hypothesis, the ssh increment
63      !!      is computed by integrating the horizontal divergence and multiply by
64      !!      by the time step.
65      !!
66      !! ** action  :   ssh(:,:,Kaa), after sea surface height
67      !!
68      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
69      !!----------------------------------------------------------------------
70      INTEGER                         , INTENT(in   ) ::   kt             ! time step
71      INTEGER                         , INTENT(in   ) ::   Kbb, Kmm, Kaa  ! time level index
72      REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) ::   pssh           ! sea-surface height
73      !
74      INTEGER  ::   jk            ! dummy loop indice
75      REAL(wp) ::   z2dt, zcoef   ! local scalars
76      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
77      !!----------------------------------------------------------------------
78      !
79      IF( ln_timing )   CALL timing_start('ssh_nxt')
80      !
81      IF( kt == nit000 ) THEN
82         IF(lwp) WRITE(numout,*)
83         IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height'
84         IF(lwp) WRITE(numout,*) '~~~~~~~ '
85      ENDIF
86      !
87      z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog)
88      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
89      zcoef = 0.5_wp * r1_rau0
90
91      !                                           !------------------------------!
92      !                                           !   After Sea Surface Height   !
93      !                                           !------------------------------!
94      IF(ln_wd_il) THEN
95         CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt, Kmm, uu, vv )
96      ENDIF
97
98      CALL div_hor( kt, Kbb, Kmm )                     ! Horizontal divergence
99      !
100      zhdiv(:,:) = 0._wp
101      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
102        zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk)
103      END DO
104      !                                                ! Sea surface elevation time stepping
105      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
106      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
107      !
108      pssh(:,:,Kaa) = (  pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
109      !
110#if defined key_agrif
111      Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt )
112#endif
113      !
114      IF ( .NOT.ln_dynspg_ts ) THEN
115         IF( ln_bdy ) THEN
116            CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. )    ! Not sure that's necessary
117            CALL bdy_ssh( pssh(:,:,Kaa) )             ! Duplicate sea level across open boundaries
118         ENDIF
119      ENDIF
120      !                                           !------------------------------!
121      !                                           !           outputs            !
122      !                                           !------------------------------!
123      !
124      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa)  - : ', mask1=tmask )
125      !
126      IF( ln_timing )   CALL timing_stop('ssh_nxt')
127      !
128   END SUBROUTINE ssh_nxt
129
130   
131   SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa )
132      !!----------------------------------------------------------------------
133      !!                ***  ROUTINE wzv  ***
134      !!                   
135      !! ** Purpose :   compute the now vertical velocity
136      !!
137      !! ** Method  : - Using the incompressibility hypothesis, the vertical
138      !!      velocity is computed by integrating the horizontal divergence 
139      !!      from the bottom to the surface minus the scale factor evolution.
140      !!        The boundary conditions are w=0 at the bottom (no flux) and.
141      !!
142      !! ** action  :   pww      : now vertical velocity
143      !!
144      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
145      !!----------------------------------------------------------------------
146      INTEGER                         , INTENT(in)    ::   kt             ! time step
147      INTEGER                         , INTENT(in)    ::   Kbb, Kmm, Kaa  ! time level indices
148      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pww            ! now vertical velocity
149      !
150      INTEGER  ::   ji, jj, jk   ! dummy loop indices
151      REAL(wp) ::   z1_2dt       ! local scalars
152      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
153      !!----------------------------------------------------------------------
154      !
155      IF( ln_timing )   CALL timing_start('wzv')
156      !
157      IF( kt == nit000 ) THEN
158         IF(lwp) WRITE(numout,*)
159         IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
160         IF(lwp) WRITE(numout,*) '~~~~~ '
161         !
162         pww(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
163      ENDIF
164      !                                           !------------------------------!
165      !                                           !     Now Vertical Velocity    !
166      !                                           !------------------------------!
167      z1_2dt = 1. / ( 2. * rdt )                         ! set time step size (Euler/Leapfrog)
168      IF( neuler == 0 .AND. kt == nit000 )   z1_2dt = 1. / rdt
169      !
170      IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN      ! z_tilde and layer cases
171         ALLOCATE( zhdiv(jpi,jpj,jpk) ) 
172         !
173         DO jk = 1, jpkm1
174            ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
175            ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
176            DO jj = 2, jpjm1
177               DO ji = fs_2, fs_jpim1   ! vector opt.
178                  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) )
179               END DO
180            END DO
181         END DO
182         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
183         !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
184         !                             ! Same question holds for hdiv. Perhaps just for security
185         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
186            ! computation of w
187            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk)    &
188               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )     ) * tmask(:,:,jk)
189         END DO
190         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
191         DEALLOCATE( zhdiv ) 
192      ELSE   ! z_star and linear free surface cases
193         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
194            ! computation of w
195            pww(:,:,jk) = pww(:,:,jk+1) - (  e3t(:,:,jk,Kmm) * hdiv(:,:,jk)                 &
196               &                         + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) )  ) * tmask(:,:,jk)
197         END DO
198      ENDIF
199
200      IF( ln_bdy ) THEN
201         DO jk = 1, jpkm1
202            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
203         END DO
204      ENDIF
205      !
206#if defined key_agrif 
207      IF( .NOT. AGRIF_Root() ) THEN
208         IF ((nbondi ==  1).OR.(nbondi == 2)) pww(nlci-1 , :     ,:) = 0.e0      ! east
209         IF ((nbondi == -1).OR.(nbondi == 2)) pww(2      , :     ,:) = 0.e0      ! west
210         IF ((nbondj ==  1).OR.(nbondj == 2)) pww(:      ,nlcj-1 ,:) = 0.e0      ! north
211         IF ((nbondj == -1).OR.(nbondj == 2)) pww(:      ,2      ,:) = 0.e0      ! south
212      ENDIF 
213#endif 
214      !
215      IF( ln_timing )   CALL timing_stop('wzv')
216      !
217   END SUBROUTINE wzv
218
219
220   SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
221      !!----------------------------------------------------------------------
222      !!                    ***  ROUTINE ssh_atf  ***
223      !!
224      !!      !!!!!!!!!!!! REWRITE HEADER COMMENTS !!!!!!!!!!!!
225      !!
226      !! ** Purpose :   achieve the sea surface  height time stepping by
227      !!              applying Asselin time filter and swapping the arrays
228      !!              pssh(:,:,Kaa)  already computed in ssh_nxt 
229      !!
230      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
231      !!              from the filter, see Leclair and Madec 2010) and swap :
232      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa) + atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) )
233      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
234      !!                pssh(:,:,Kmm) = pssh(:,:,Kaa)
235      !!
236      !! ** action  : - pssh(:,:,Kbb), pssh(:,:,Kmm)   : before & now sea surface height
237      !!                               ready for the next time step
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_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
264         ENDIF
265      ENDIF
266      !
267      IF(ln_ctl)   CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm)  - : ', mask1=tmask )
268      !
269      IF( ln_timing )   CALL timing_stop('ssh_atf')
270      !
271   END SUBROUTINE ssh_atf
272
273   SUBROUTINE wAimp( kt, Kmm )
274      !!----------------------------------------------------------------------
275      !!                ***  ROUTINE wAimp  ***
276      !!                   
277      !! ** Purpose :   compute the Courant number and partition vertical velocity
278      !!                if a proportion needs to be treated implicitly
279      !!
280      !! ** Method  : -
281      !!
282      !! ** action  :   ww      : now vertical velocity (to be handled explicitly)
283      !!            :   wi      : now vertical velocity (for implicit treatment)
284      !!
285      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
286      !!----------------------------------------------------------------------
287      INTEGER, INTENT(in) ::   kt   ! time step
288      INTEGER, INTENT(in) ::   Kmm  ! time level index
289      !
290      INTEGER  ::   ji, jj, jk   ! dummy loop indices
291      REAL(wp)             ::   zCu, zcff, z1_e3w                     ! local scalars
292      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
293      REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters
294      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
295      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
296      !!----------------------------------------------------------------------
297      !
298      IF( ln_timing )   CALL timing_start('wAimp')
299      !
300      IF( kt == nit000 ) THEN
301         IF(lwp) WRITE(numout,*)
302         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
303         IF(lwp) WRITE(numout,*) '~~~~~ '
304         !
305         Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all)
306      ENDIF
307      !
308      DO jk = 1, jpkm1            ! calculate Courant numbers
309         DO jj = 2, jpjm1
310            DO ji = 2, fs_jpim1   ! vector opt.
311               z1_e3w = 1._wp / e3w(ji,jj,jk,Kmm)
312               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )    &
313                  &                      + ( MAX( e2u(ji  ,jj)*e3uw(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
314                  &                          MIN( e2u(ji-1,jj)*e3uw(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
315                  &                        * r1_e1e2t(ji,jj)                                                  &
316                  &                      + ( MAX( e1v(ji,jj  )*e3vw(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   &
317                  &                          MIN( e1v(ji,jj-1)*e3vw(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   &
318                  &                        * r1_e1e2t(ji,jj)                                                  &
319                  &                      ) * z1_e3w
320            END DO
321         END DO
322      END DO
323      !
324      CALL iom_put("Courant",Cu_adv)
325      !
326      wi(:,:,:) = 0._wp                                 ! Includes top and bottom values set to zero
327      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
328         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition
329            DO jj = 2, jpjm1                            ! w where necessary
330               DO ji = 2, fs_jpim1   ! vector opt.
331                  !
332                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
333                  !
334                  IF( zCu < Cu_min ) THEN               !<-- Fully explicit
335                     zcff = 0._wp
336                  ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
337                     zcff = ( zCu - Cu_min )**2
338                     zcff = zcff / ( Fcu + zcff )
339                  ELSE                                  !<-- Mostly implicit
340                     zcff = ( zCu - Cu_max )/ zCu
341                  ENDIF
342                  zcff = MIN(1._wp, zcff)
343                  !
344                  wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
345                  ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
346                  !
347                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient
348               END DO
349            END DO
350         END DO
351      ELSE
352         ! Fully explicit everywhere
353         Cu_adv = 0.0_wp                                ! Reuse array to output coefficient
354      ENDIF
355      CALL iom_put("wimp",wi) 
356      CALL iom_put("wi_cff",Cu_adv)
357      CALL iom_put("wexp",ww)
358      !
359      IF( ln_timing )   CALL timing_stop('wAimp')
360      !
361   END SUBROUTINE wAimp
362   !!======================================================================
363END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.