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/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DYN – NEMO

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DYN/sshwzv.F90 @ 10986

Last change on this file since 10986 was 10986, checked in by andmirek, 5 years ago

GMED 462 add flush

File size: 17.1 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_swp       : 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_swp    ! 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 )
57      !!----------------------------------------------------------------------
58      !!                ***  ROUTINE ssh_nxt  ***
59      !!                   
60      !! ** Purpose :   compute the after ssh (ssha)
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  :   ssha, after sea surface height
67      !!
68      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
69      !!----------------------------------------------------------------------
70      INTEGER, INTENT(in) ::   kt   ! time step
71      !
72      INTEGER  ::   jk            ! dummy loop indice
73      REAL(wp) ::   z2dt, zcoef   ! local scalars
74      REAL(wp), DIMENSION(jpi,jpj) ::   zhdiv   ! 2D workspace
75      !!----------------------------------------------------------------------
76      !
77      IF( ln_timing )   CALL timing_start('ssh_nxt')
78      !
79      IF( kt == nit000 .AND. lwp) THEN
80         WRITE(numout,*)
81         WRITE(numout,*) 'ssh_nxt : after sea surface height'
82         WRITE(numout,*) '~~~~~~~ '
83         IF(lflush) CALL FLUSH(numout)
84      ENDIF
85      !
86      z2dt = 2._wp * rdt                          ! set time step size (Euler/Leapfrog)
87      IF( neuler == 0 .AND. kt == nit000 )   z2dt = rdt
88      zcoef = 0.5_wp * r1_rau0
89
90      !                                           !------------------------------!
91      !                                           !   After Sea Surface Height   !
92      !                                           !------------------------------!
93      IF(ln_wd_il) THEN
94         CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)
95      ENDIF
96
97      CALL div_hor( kt )                               ! Horizontal divergence
98      !
99      zhdiv(:,:) = 0._wp
100      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
101        zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk)
102      END DO
103      !                                                ! Sea surface elevation time stepping
104      ! In time-split case we need a first guess of the ssh after (using the baroclinic timestep) in order to
105      ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
106      !
107      ssha(:,:) = (  sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * ssmask(:,:)
108      !
109#if defined key_agrif
110      CALL agrif_ssh( kt )
111#endif
112      !
113      IF ( .NOT.ln_dynspg_ts ) THEN
114         IF( ln_bdy ) THEN
115            CALL lbc_lnk( 'sshwzv', ssha, 'T', 1. )    ! Not sure that's necessary
116            CALL bdy_ssh( ssha )             ! Duplicate sea level across open boundaries
117         ENDIF
118      ENDIF
119      !                                           !------------------------------!
120      !                                           !           outputs            !
121      !                                           !------------------------------!
122      !
123      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha  - : ', mask1=tmask )
124      !
125      IF( ln_timing )   CALL timing_stop('ssh_nxt')
126      !
127   END SUBROUTINE ssh_nxt
128
129   
130   SUBROUTINE wzv( kt )
131      !!----------------------------------------------------------------------
132      !!                ***  ROUTINE wzv  ***
133      !!                   
134      !! ** Purpose :   compute the now vertical velocity
135      !!
136      !! ** Method  : - Using the incompressibility hypothesis, the vertical
137      !!      velocity is computed by integrating the horizontal divergence 
138      !!      from the bottom to the surface minus the scale factor evolution.
139      !!        The boundary conditions are w=0 at the bottom (no flux) and.
140      !!
141      !! ** action  :   wn      : now vertical velocity
142      !!
143      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
144      !!----------------------------------------------------------------------
145      INTEGER, INTENT(in) ::   kt   ! time step
146      !
147      INTEGER  ::   ji, jj, jk   ! dummy loop indices
148      REAL(wp) ::   z1_2dt       ! local scalars
149      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zhdiv
150      !!----------------------------------------------------------------------
151      !
152      IF( ln_timing )   CALL timing_start('wzv')
153      !
154      IF( kt == nit000 ) THEN
155         IF(lwp) THEN
156            WRITE(numout,*)
157            WRITE(numout,*) 'wzv : now vertical velocity '
158            WRITE(numout,*) '~~~~~ '
159            IF(lflush) CALL FLUSH(numout)
160         ENDIF
161         !
162         wn(:,:,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 hdivn. Perhaps just for security
185         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
186            ! computation of w
187            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)    &
188               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )     ) * tmask(:,:,jk)
189         END DO
190         !          IF( ln_vvl_layer ) wn(:,:,:) = 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            wn(:,:,jk) = wn(:,:,jk+1) - (  e3t_n(:,:,jk) * hdivn(:,:,jk)                 &
196               &                         + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) )  ) * tmask(:,:,jk)
197         END DO
198      ENDIF
199
200      IF( ln_bdy ) THEN
201         DO jk = 1, jpkm1
202            wn(:,:,jk) = wn(:,:,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)) wn(nlci-1 , :     ,:) = 0.e0      ! east
209         IF ((nbondi == -1).OR.(nbondi == 2)) wn(2      , :     ,:) = 0.e0      ! west
210         IF ((nbondj ==  1).OR.(nbondj == 2)) wn(:      ,nlcj-1 ,:) = 0.e0      ! north
211         IF ((nbondj == -1).OR.(nbondj == 2)) wn(:      ,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_swp( kt )
221      !!----------------------------------------------------------------------
222      !!                    ***  ROUTINE ssh_nxt  ***
223      !!
224      !! ** Purpose :   achieve the sea surface  height time stepping by
225      !!              applying Asselin time filter and swapping the arrays
226      !!              ssha  already computed in ssh_nxt 
227      !!
228      !! ** Method  : - apply Asselin time fiter to now ssh (excluding the forcing
229      !!              from the filter, see Leclair and Madec 2010) and swap :
230      !!                sshn = ssha + atfp * ( sshb -2 sshn + ssha )
231      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
232      !!                sshn = ssha
233      !!
234      !! ** action  : - sshb, sshn   : before & now sea surface height
235      !!                               ready for the next time step
236      !!
237      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
238      !!----------------------------------------------------------------------
239      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
240      !
241      REAL(wp) ::   zcoef   ! local scalar
242      !!----------------------------------------------------------------------
243      !
244      IF( ln_timing )   CALL timing_start('ssh_swp')
245      !
246      IF( kt == nit000 .AND. lwp) THEN
247         WRITE(numout,*)
248         WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
249         WRITE(numout,*) '~~~~~~~ '
250         IF(lflush) CALL FLUSH(numout)
251      ENDIF
252      !              !==  Euler time-stepping: no filter, just swap  ==!
253      IF ( neuler == 0 .AND. kt == nit000 ) THEN
254         sshn(:,:) = ssha(:,:)                              ! now    <-- after  (before already = now)
255         !
256      ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==!
257         !                                                  ! before <-- now filtered
258         sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )
259         IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed
260            zcoef = atfp * rdt * r1_rau0
261            sshb(:,:) = sshb(:,:) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
262               &                             -    rnf_b(:,:) + rnf   (:,:)   &
263               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
264         ENDIF
265         sshn(:,:) = ssha(:,:)                              ! now <-- after
266      ENDIF
267      !
268      IF(ln_ctl)   CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb  - : ', mask1=tmask )
269      !
270      IF( ln_timing )   CALL timing_stop('ssh_swp')
271      !
272   END SUBROUTINE ssh_swp
273
274   SUBROUTINE wAimp( kt )
275      !!----------------------------------------------------------------------
276      !!                ***  ROUTINE wAimp  ***
277      !!                   
278      !! ** Purpose :   compute the Courant number and partition vertical velocity
279      !!                if a proportion needs to be treated implicitly
280      !!
281      !! ** Method  : -
282      !!
283      !! ** action  :   wn      : now vertical velocity (to be handled explicitly)
284      !!            :   wi      : now vertical velocity (for implicit treatment)
285      !!
286      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
287      !!----------------------------------------------------------------------
288      INTEGER, INTENT(in) ::   kt   ! time step
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) THEN
302            WRITE(numout,*)
303            WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
304            WRITE(numout,*) '~~~~~ '
305            IF(lflush) CALL FLUSH(numout)
306         ENDIF
307         !
308         Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all)
309      ENDIF
310      !
311      DO jk = 1, jpkm1            ! calculate Courant numbers
312         DO jj = 2, jpjm1
313            DO ji = 2, fs_jpim1   ! vector opt.
314               z1_e3w = 1._wp / e3w_n(ji,jj,jk)
315               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) )    &
316                  &                      + ( MAX( e2u(ji  ,jj)*e3uw_n(ji  ,jj,jk)*un(ji  ,jj,jk), 0._wp ) -   &
317                  &                          MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) )   &
318                  &                        * r1_e1e2t(ji,jj)                                                  &
319                  &                      + ( MAX( e1v(ji,jj  )*e3vw_n(ji,jj  ,jk)*vn(ji,jj  ,jk), 0._wp ) -   &
320                  &                          MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) )   &
321                  &                        * r1_e1e2t(ji,jj)                                                  &
322                  &                      ) * z1_e3w
323            END DO
324         END DO
325      END DO
326      !
327      CALL iom_put("Courant",Cu_adv)
328      !
329      wi(:,:,:) = 0._wp                                 ! Includes top and bottom values set to zero
330      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
331         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition
332            DO jj = 2, jpjm1                            ! w where necessary
333               DO ji = 2, fs_jpim1   ! vector opt.
334                  !
335                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
336                  !
337                  IF( zCu < Cu_min ) THEN               !<-- Fully explicit
338                     zcff = 0._wp
339                  ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
340                     zcff = ( zCu - Cu_min )**2
341                     zcff = zcff / ( Fcu + zcff )
342                  ELSE                                  !<-- Mostly implicit
343                     zcff = ( zCu - Cu_max )/ zCu
344                  ENDIF
345                  zcff = MIN(1._wp, zcff)
346                  !
347                  wi(ji,jj,jk) =           zcff   * wn(ji,jj,jk)
348                  wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)
349                  !
350                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient
351               END DO
352            END DO
353         END DO
354      ELSE
355         ! Fully explicit everywhere
356         Cu_adv = 0.0_wp                                ! Reuse array to output coefficient
357      ENDIF
358      CALL iom_put("wimp",wi) 
359      CALL iom_put("wi_cff",Cu_adv)
360      CALL iom_put("wexp",wn)
361      !
362      IF( ln_timing )   CALL timing_stop('wAimp')
363      !
364   END SUBROUTINE wAimp
365   !!======================================================================
366END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.