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

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

Last change on this file since 10978 was 10978, checked in by davestorkey, 4 years ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : DOM and sshwzv.F90. (sshwzv.F90 will be rewritten somewhat when we rewrite the time-level swapping but I've done it anyway). Passes all non-AGRIF SETTE tests.

  • Property svn:keywords set to Id
File size: 17.8 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, 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)
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      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_swp( kt, Kbb, Kmm, Kaa )
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      !!              ssh(:,:,Kaa)  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      !!                ssh(:,:,Kmm) = ssh(:,:,Kaa) + atfp * ( ssh(:,:,Kbb) -2 ssh(:,:,Kmm) + ssh(:,:,Kaa) )
231      !!                            - atfp * rdt * ( emp_b - emp ) / rau0
232      !!                ssh(:,:,Kmm) = ssh(:,:,Kaa)
233      !!
234      !! ** action  : - ssh(:,:,Kbb), ssh(:,:,Kmm)   : 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      INTEGER, INTENT(in) ::   Kbb, Kmm, Kaa   ! ocean time-step index
241      !
242      REAL(wp) ::   zcoef   ! local scalar
243      !!----------------------------------------------------------------------
244      !
245      IF( ln_timing )   CALL timing_start('ssh_swp')
246      !
247      IF( kt == nit000 ) THEN
248         IF(lwp) WRITE(numout,*)
249         IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'
250         IF(lwp) WRITE(numout,*) '~~~~~~~ '
251      ENDIF
252      !              !==  Euler time-stepping: no filter, just swap  ==!
253      IF ( neuler == 0 .AND. kt == nit000 ) THEN
254         ssh(:,:,Kmm) = ssh(:,:,Kaa)                              ! now    <-- after  (before already = now)
255         !
256      ELSE           !==  Leap-Frog time-stepping: Asselin filter + swap  ==!
257         !                                                  ! before <-- now filtered
258         ssh(:,:,Kbb) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa) )
259         IF( .NOT.ln_linssh ) THEN                          ! before <-- with forcing removed
260            zcoef = atfp * rdt * r1_rau0
261            ssh(:,:,Kbb) = ssh(:,:,Kbb) - zcoef * (     emp_b(:,:) - emp   (:,:)   &
262               &                             -    rnf_b(:,:) + rnf   (:,:)   &
263               &                             + fwfisf_b(:,:) - fwfisf(:,:)   ) * ssmask(:,:)
264         ENDIF
265         ssh(:,:,Kmm) = ssh(:,:,Kaa)                              ! now <-- after
266      ENDIF
267      !
268      IF(ln_ctl)   CALL prt_ctl( tab2d_1=ssh(:,:,Kbb), clinfo1=' ssh(:,:,Kbb)  - : ', mask1=tmask )
269      !
270      IF( ln_timing )   CALL timing_stop('ssh_swp')
271      !
272   END SUBROUTINE ssh_swp
273
274   SUBROUTINE wAimp( kt, Kmm )
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  :   ww      : 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      INTEGER, INTENT(in) ::   Kmm  ! time level index
290      !
291      INTEGER  ::   ji, jj, jk   ! dummy loop indices
292      REAL(wp)             ::   zCu, zcff, z1_e3w                     ! local scalars
293      REAL(wp) , PARAMETER ::   Cu_min = 0.15_wp                      ! local parameters
294      REAL(wp) , PARAMETER ::   Cu_max = 0.27                         ! local parameters
295      REAL(wp) , PARAMETER ::   Cu_cut = 2._wp*Cu_max - Cu_min        ! local parameters
296      REAL(wp) , PARAMETER ::   Fcu    = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters
297      !!----------------------------------------------------------------------
298      !
299      IF( ln_timing )   CALL timing_start('wAimp')
300      !
301      IF( kt == nit000 ) THEN
302         IF(lwp) WRITE(numout,*)
303         IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity '
304         IF(lwp) WRITE(numout,*) '~~~~~ '
305         !
306         Cu_adv(:,:,jpk) = 0._wp              ! bottom value : Cu_adv=0 (set once for all)
307      ENDIF
308      !
309      DO jk = 1, jpkm1            ! calculate Courant numbers
310         DO jj = 2, jpjm1
311            DO ji = 2, fs_jpim1   ! vector opt.
312               z1_e3w = 1._wp / e3w(ji,jj,jk,Kmm)
313               Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) )    &
314                  &                      + ( MAX( e2u(ji  ,jj)*e3uw(ji  ,jj,jk,Kmm)*uu(ji  ,jj,jk,Kmm), 0._wp ) -   &
315                  &                          MIN( e2u(ji-1,jj)*e3uw(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) )   &
316                  &                        * r1_e1e2t(ji,jj)                                                  &
317                  &                      + ( MAX( e1v(ji,jj  )*e3vw(ji,jj  ,jk,Kmm)*vv(ji,jj  ,jk,Kmm), 0._wp ) -   &
318                  &                          MIN( e1v(ji,jj-1)*e3vw(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) )   &
319                  &                        * r1_e1e2t(ji,jj)                                                  &
320                  &                      ) * z1_e3w
321            END DO
322         END DO
323      END DO
324      !
325      CALL iom_put("Courant",Cu_adv)
326      !
327      wi(:,:,:) = 0._wp                                 ! Includes top and bottom values set to zero
328      IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN       ! Quick check if any breaches anywhere
329         DO jk = 1, jpkm1                               ! or scan Courant criterion and partition
330            DO jj = 2, jpjm1                            ! w where necessary
331               DO ji = 2, fs_jpim1   ! vector opt.
332                  !
333                  zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) )
334                  !
335                  IF( zCu < Cu_min ) THEN               !<-- Fully explicit
336                     zcff = 0._wp
337                  ELSEIF( zCu < Cu_cut ) THEN           !<-- Mixed explicit
338                     zcff = ( zCu - Cu_min )**2
339                     zcff = zcff / ( Fcu + zcff )
340                  ELSE                                  !<-- Mostly implicit
341                     zcff = ( zCu - Cu_max )/ zCu
342                  ENDIF
343                  zcff = MIN(1._wp, zcff)
344                  !
345                  wi(ji,jj,jk) =           zcff   * ww(ji,jj,jk)
346                  ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk)
347                  !
348                  Cu_adv(ji,jj,jk) = zcff               ! Reuse array to output coefficient
349               END DO
350            END DO
351         END DO
352      ELSE
353         ! Fully explicit everywhere
354         Cu_adv = 0.0_wp                                ! Reuse array to output coefficient
355      ENDIF
356      CALL iom_put("wimp",wi) 
357      CALL iom_put("wi_cff",Cu_adv)
358      CALL iom_put("wexp",ww)
359      !
360      IF( ln_timing )   CALL timing_stop('wAimp')
361      !
362   END SUBROUTINE wAimp
363   !!======================================================================
364END MODULE sshwzv
Note: See TracBrowser for help on using the repository browser.