source: NEMO/trunk/src/OCE/DYN/sshwzv.F90

Last change on this file was 13497, checked in by techene, 7 days ago

re-introduce comments that have been erased by loops transformation see #2525

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