source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/OCE/TRA/traatf.F90 @ 11057

Last change on this file since 11057 was 11057, checked in by davestorkey, 18 months ago

2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap : Rewrite time filtering of tracers and SSH. This version now bit compares with the base code for ORCA1 starting from a restart.
To do:

  1. Check that it bit compares with an Euler timestep.
  2. Full SETTE tests.
  3. Check that it compiles with C1D.
  4. Check that TOP works (see 1).
  • Property svn:keywords set to Id
File size: 19.7 KB
Line 
1MODULE traatf
2   !!======================================================================
3   !!                       ***  MODULE  traatf  ***
4   !! Ocean active tracers:  time stepping on temperature and salinity
5   !!======================================================================
6   !! History :  OPA  !  1991-11  (G. Madec)  Original code
7   !!            7.0  !  1993-03  (M. Guyon)  symetrical conditions
8   !!            8.0  !  1996-02  (G. Madec & M. Imbard)  opa release 8.0
9   !!             -   !  1996-04  (A. Weaver)  Euler forward step
10   !!            8.2  !  1999-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
11   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
12   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
13   !!             -   !  2005-04  (C. Deltel) Add Asselin trend in the ML budget
14   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
15   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazdf
16   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
17   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  semi-implicit hpg with asselin filter + modified LF-RA
18   !!             -   !  2010-05  (C. Ethe, G. Madec)  merge TRC-TRA
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   tra_atf       : time stepping on tracers
23   !!   tra_atf_fix   : time stepping on tracers : fixed    volume case
24   !!   tra_atf_vvl   : time stepping on tracers : variable volume case
25   !!----------------------------------------------------------------------
26   USE oce             ! ocean dynamics and tracers variables
27   USE dom_oce         ! ocean space and time domain variables
28   USE sbc_oce         ! surface boundary condition: ocean
29   USE sbcrnf          ! river runoffs
30   USE sbcisf          ! ice shelf melting
31   USE zdf_oce         ! ocean vertical mixing
32   USE domvvl          ! variable volume
33   USE trd_oce         ! trends: ocean variables
34   USE trdtra          ! trends manager: tracers
35   USE traqsr          ! penetrative solar radiation (needed for nksr)
36   USE phycst          ! physical constant
37   USE ldftra          ! lateral physics : tracers
38   USE ldfslp          ! lateral physics : slopes
39   USE bdy_oce  , ONLY : ln_bdy
40   USE bdytra          ! open boundary condition (bdy_tra routine)
41   !
42   USE in_out_manager  ! I/O manager
43   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
44   USE prtctl          ! Print control
45   USE timing          ! Timing
46#if defined key_agrif
47   USE agrif_oce_interp
48#endif
49
50   IMPLICIT NONE
51   PRIVATE
52
53   PUBLIC   tra_atf       ! routine called by step.F90
54   PUBLIC   tra_atf_fix   ! to be used in trcnxt
55   PUBLIC   tra_atf_vvl   ! to be used in trcnxt
56
57   !! * Substitutions
58#  include "vectopt_loop_substitute.h90"
59   !!----------------------------------------------------------------------
60   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
61   !! $Id$
62   !! Software governed by the CeCILL license (see ./LICENSE)
63   !!----------------------------------------------------------------------
64CONTAINS
65
66   SUBROUTINE tra_atf( kt, Kbb, Kmm, Kaa, pts )
67      !!----------------------------------------------------------------------
68      !!                   ***  ROUTINE traatf  ***
69      !!
70      !!       !!!!!!!!!!!!!!!!!  REWRITE HEADER COMMENTS !!!!!!!!!!!!!!!
71      !!
72      !! ** Purpose :   Apply the boundary condition on the after temperature 
73      !!             and salinity fields, achieved the time stepping by adding
74      !!             the Asselin filter on now fields and swapping the fields.
75      !!
76      !! ** Method  :   At this stage of the computation, ta and sa are the
77      !!             after temperature and salinity as the time stepping has
78      !!             been performed in trazdf_imp or trazdf_exp module.
79      !!
80      !!              - Apply lateral boundary conditions on (ta,sa)
81      !!             at the local domain   boundaries through lbc_lnk call,
82      !!             at the one-way open boundaries (ln_bdy=T),
83      !!             at the AGRIF zoom   boundaries (lk_agrif=T)
84      !!
85      !!              - Update lateral boundary conditions on AGRIF children
86      !!             domains (lk_agrif=T)
87      !!
88      !! ** Action  : - ts(Kbb) & ts(Kmm) ready for the next time step
89      !!----------------------------------------------------------------------
90      INTEGER                                  , INTENT(in   ) :: kt             ! ocean time-step index
91      INTEGER                                  , INTENT(in   ) :: Kbb, Kmm, Kaa  ! time level indices
92      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts            ! active tracers
93      !!
94      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
95      REAL(wp) ::   zfact            ! local scalars
96      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
97      !!----------------------------------------------------------------------
98      !
99      IF( ln_timing )   CALL timing_start( 'tra_atf')
100      !
101      IF( kt == nit000 ) THEN
102         IF(lwp) WRITE(numout,*)
103         IF(lwp) WRITE(numout,*) 'tra_atf : achieve the time stepping by Asselin filter and array swap'
104         IF(lwp) WRITE(numout,*) '~~~~~~~'
105      ENDIF
106
107      ! Update after tracer on domain lateral boundaries
108      !
109#if defined key_agrif
110      CALL Agrif_tra                     ! AGRIF zoom boundaries
111#endif
112      !                                              ! local domain boundaries  (T-point, unchanged sign)
113      CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. )
114      !
115      IF( ln_bdy )   CALL bdy_tra( kt, Kbb, pts, Kaa )  ! BDY open boundaries
116 
117      ! set time step size (Euler/Leapfrog)
118      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dt =        rdt   ! at nit000             (Euler)
119      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dt = 2._wp* rdt   ! at nit000 or nit000+1 (Leapfrog)
120      ENDIF
121
122      ! trends computation initialisation
123      IF( l_trdtra )   THEN                   
124         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
125         ztrdt(:,:,jpk) = 0._wp
126         ztrds(:,:,jpk) = 0._wp
127         IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend
128            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt )
129            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_zdfp, ztrds )
130         ENDIF
131         ! total trend for the non-time-filtered variables.
132         zfact = 1.0 / rdt
133         ! G Nurser 23 Mar 2017. Recalculate trend as Delta(e3t*T)/e3tn; e3tn cancel from pts(Kmm) terms
134         DO jk = 1, jpkm1
135            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_tem,Kmm)) * zfact
136            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) / e3t(:,:,jk,Kmm) - pts(:,:,jk,jp_sal,Kmm)) * zfact
137         END DO
138         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_tot, ztrdt )
139         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_tot, ztrds )
140         IF( ln_linssh ) THEN       ! linear sea surface height only
141            ! Store now fields before applying the Asselin filter
142            ! in order to calculate Asselin filter trend later.
143            ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kmm) 
144            ztrds(:,:,:) = pts(:,:,:,jp_sal,Kmm)
145         ENDIF
146      ENDIF
147
148      IF( neuler == 0 .AND. kt == nit000 ) THEN       ! Euler time-stepping at first time-step (only swap)
149         DO jn = 1, jpts
150            DO jk = 1, jpkm1
151               pts(:,:,jk,jn,Kmm) = pts(:,:,jk,jn,Kaa)   
152            END DO
153         END DO
154         IF (l_trdtra .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
155            !                                        ! Asselin filter is output by tra_atf_vvl that is not called on this time step
156            ztrdt(:,:,:) = 0._wp
157            ztrds(:,:,:) = 0._wp
158            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt )
159            CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_atf, ztrds )
160         END IF
161         !
162      ELSE                                            ! Leap-Frog + Asselin filter time stepping
163         !
164         IF( ln_linssh ) THEN   ;   CALL tra_atf_fix( kt, Kbb, Kmm, Kaa, nit000,      'TRA', pts, jpts )  ! linear free surface
165         ELSE                   ;   CALL tra_atf_vvl( kt, Kbb, Kmm, Kaa, nit000, rdt, 'TRA', pts, sbc_tsc, sbc_tsc_b, jpts )  ! non-linear free surface
166         ENDIF
167         !
168         CALL lbc_lnk_multi( 'traatf', pts(:,:,:,jp_tem,Kbb) , 'T', 1., pts(:,:,:,jp_sal,Kbb) , 'T', 1., &
169                  &                    pts(:,:,:,jp_tem,Kmm) , 'T', 1., pts(:,:,:,jp_sal,Kmm) , 'T', 1., &
170                  &                    pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1.  )
171         !
172      ENDIF     
173      !
174      IF( l_trdtra .AND. ln_linssh ) THEN      ! trend of the Asselin filter (tb filtered - tb)/dt     
175         zfact = 1._wp / r2dt             
176         DO jk = 1, jpkm1
177            ztrdt(:,:,jk) = ( pts(:,:,jk,jp_tem,Kbb) - ztrdt(:,:,jk) ) * zfact
178            ztrds(:,:,jk) = ( pts(:,:,jk,jp_sal,Kbb) - ztrds(:,:,jk) ) * zfact
179         END DO
180         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_atf, ztrdt )
181         CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_sal, jptra_atf, ztrds )
182      END IF
183      IF( l_trdtra )   DEALLOCATE( ztrdt , ztrds )
184      !
185      !                        ! control print
186      IF(ln_ctl)   CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Kmm), clinfo1=' nxt  - Tn: ', mask1=tmask,   &
187         &                       tab3d_2=pts(:,:,:,jp_sal,Kmm), clinfo2=       ' Sn: ', mask2=tmask )
188      !
189      IF( ln_timing )   CALL timing_stop('tra_atf')
190      !
191   END SUBROUTINE tra_atf
192
193
194   SUBROUTINE tra_atf_fix( kt, kit000, Kbb, Kmm, Kaa, cdtype, pt, kjpt )
195      !!----------------------------------------------------------------------
196      !!                   ***  ROUTINE tra_atf_fix  ***
197      !!
198      !! ** Purpose :   fixed volume: apply the Asselin time filter and
199      !!                swap the tracer fields.
200      !!
201      !! ** Method  : - Apply a Asselin time filter on now fields.
202      !!              - swap tracer fields to prepare the next time_step.
203      !!
204      !! ** Action  : - ptb & ptn ready for the next time step
205      !!----------------------------------------------------------------------
206      INTEGER                                  , INTENT(in   ) ::  kt            ! ocean time-step index
207      INTEGER                                  , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
208      INTEGER                                  , INTENT(in   ) ::  kit000        ! first time step index
209      CHARACTER(len=3)                         , INTENT(in   ) ::  cdtype        ! =TRA or TRC (tracer indicator)
210      INTEGER                                  , INTENT(in   ) ::  kjpt          ! number of tracers
211      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::  pt            ! tracer fields
212      !
213      INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
214      REAL(wp) ::   ztn, ztd         ! local scalars
215      !!----------------------------------------------------------------------
216      !
217      IF( kt == kit000 )  THEN
218         IF(lwp) WRITE(numout,*)
219         IF(lwp) WRITE(numout,*) 'tra_atf_fix : time stepping', cdtype
220         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
221      ENDIF
222      !
223      DO jn = 1, kjpt
224         !
225         DO jk = 1, jpkm1
226            DO jj = 2, jpjm1
227               DO ji = fs_2, fs_jpim1
228                  ztn = pt(ji,jj,jk,jn,Kmm)                                   
229                  ztd = pt(ji,jj,jk,jn,Kaa) - 2._wp * ztn + pt(ji,jj,jk,jn,Kbb)  ! time laplacian on tracers
230                  !
231                  pt(ji,jj,jk,jn,Kmm) = ztn + atfp * ztd                      ! ptb <-- filtered ptn
232               END DO
233           END DO
234         END DO
235         !
236      END DO
237      !
238   END SUBROUTINE tra_atf_fix
239
240
241   SUBROUTINE tra_atf_vvl( kt, Kbb, Kmm, Kaa, kit000, p2dt, cdtype, pt, psbc_tc, psbc_tc_b, kjpt )
242      !!----------------------------------------------------------------------
243      !!                   ***  ROUTINE tra_atf_vvl  ***
244      !!
245      !! ** Purpose :   Time varying volume: apply the Asselin time filter 
246      !!                and swap the tracer fields.
247      !!
248      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
249      !!              - swap tracer fields to prepare the next time_step.
250      !!             tb  = ( e3t(Kmm)*tn + atfp*[ e3t(Kbb)*tb - 2 e3t(Kmm)*tn + e3t_a*ta ] )
251      !!                  /( e3t(Kmm)    + atfp*[ e3t(Kbb)    - 2 e3t(Kmm)    + e3t(Kaa)    ] )
252      !!             tn  = ta
253      !!
254      !! ** Action  : - ptb & ptn ready for the next time step
255      !!----------------------------------------------------------------------
256      INTEGER                                  , INTENT(in   ) ::  kt        ! ocean time-step index
257      INTEGER                                  , INTENT(in   ) ::  Kbb, Kmm, Kaa ! time level indices
258      INTEGER                                  , INTENT(in   ) ::  kit000    ! first time step index
259      REAL(wp)                                 , INTENT(in   ) ::  p2dt      ! time-step
260      CHARACTER(len=3)                         , INTENT(in   ) ::  cdtype    ! =TRA or TRC (tracer indicator)
261      INTEGER                                  , INTENT(in   ) ::  kjpt      ! number of tracers
262      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) ::  pt        ! tracer fields
263      REAL(wp), DIMENSION(jpi,jpj    ,kjpt)    , INTENT(in   ) ::  psbc_tc   ! surface tracer content
264      REAL(wp), DIMENSION(jpi,jpj    ,kjpt)    , INTENT(in   ) ::  psbc_tc_b ! before surface tracer content
265      !
266      LOGICAL  ::   ll_traqsr, ll_rnf, ll_isf   ! local logical
267      INTEGER  ::   ji, jj, jk, jn              ! dummy loop indices
268      REAL(wp) ::   zfact, zfact1, ztc_a , ztc_n , ztc_b , ztc_f , ztc_d    ! local scalar
269      REAL(wp) ::   zfact2, ze3t_b, ze3t_n, ze3t_a, ze3t_f, ze3t_d   !   -      -
270      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd_atf
271      !!----------------------------------------------------------------------
272      !
273      IF( kt == kit000 )  THEN
274         IF(lwp) WRITE(numout,*)
275         IF(lwp) WRITE(numout,*) 'tra_atf_vvl : time stepping', cdtype
276         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
277      ENDIF
278      !
279      IF( cdtype == 'TRA' )  THEN   
280         ll_traqsr  = ln_traqsr        ! active  tracers case  and  solar penetration
281         ll_rnf     = ln_rnf           ! active  tracers case  and  river runoffs
282         ll_isf     = ln_isf           ! active  tracers case  and  ice shelf melting
283      ELSE                          ! passive tracers case
284         ll_traqsr  = .FALSE.          ! NO solar penetration
285         ll_rnf     = .FALSE.          ! NO river runoffs ????          !!gm BUG ? 
286         ll_isf     = .FALSE.          ! NO ice shelf melting/freezing  !!gm BUG ??
287      ENDIF
288      !
289      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
290         ALLOCATE( ztrd_atf(jpi,jpj,jpk,kjpt) )
291         ztrd_atf(:,:,:,:) = 0.0_wp
292      ENDIF
293      zfact = 1._wp / p2dt
294      zfact1 = atfp * p2dt
295      zfact2 = zfact1 * r1_rau0
296      DO jn = 1, kjpt     
297         DO jk = 1, jpkm1
298            DO jj = 2, jpjm1
299               DO ji = fs_2, fs_jpim1
300                  ze3t_b = e3t(ji,jj,jk,Kbb)
301                  ze3t_n = e3t(ji,jj,jk,Kmm)
302                  ze3t_a = e3t(ji,jj,jk,Kaa)
303                  !                                         ! tracer content at Before, now and after
304                  ztc_b  = pt(ji,jj,jk,jn,Kbb) * ze3t_b
305                  ztc_n  = pt(ji,jj,jk,jn,Kmm) * ze3t_n
306                  ztc_a  = pt(ji,jj,jk,jn,Kaa) * ze3t_a
307                  !
308                  ze3t_d = ze3t_a - 2. * ze3t_n + ze3t_b
309                  ztc_d  = ztc_a  - 2. * ztc_n  + ztc_b
310                  !
311                  ze3t_f = ze3t_n + atfp * ze3t_d
312                  ztc_f  = ztc_n  + atfp * ztc_d
313                  !
314                  IF( jk == mikt(ji,jj) ) THEN           ! first level
315                     ze3t_f = ze3t_f - zfact2 * ( (emp_b(ji,jj)    - emp(ji,jj)   )  &
316                            &                   + (fwfisf_b(ji,jj) - fwfisf(ji,jj))  )
317                     ztc_f  = ztc_f  - zfact1 * ( psbc_tc(ji,jj,jn) - psbc_tc_b(ji,jj,jn) )
318                  ENDIF
319                  IF( ln_rnf_depth ) THEN
320                     ! Rivers are not just at the surface must go down to nk_rnf(ji,jj)
321                     IF( mikt(ji,jj) <=jk .and. jk <= nk_rnf(ji,jj)  ) THEN
322                        ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj) - rnf(ji,jj)   )  ) &
323                    &                            * ( e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj) ) 
324                     ENDIF
325                  ELSE
326                     IF( jk == mikt(ji,jj) ) THEN           ! first level
327                        ze3t_f = ze3t_f - zfact2 * ( - (rnf_b(ji,jj)    - rnf(ji,jj)   ) ) 
328                     ENDIF
329                  ENDIF
330
331                  !
332                  ! solar penetration (temperature only)
333                  IF( ll_traqsr .AND. jn == jp_tem .AND. jk <= nksr )                            & 
334                     &     ztc_f  = ztc_f  - zfact1 * ( qsr_hc(ji,jj,jk) - qsr_hc_b(ji,jj,jk) ) 
335                     !
336                  ! river runoff
337                  IF( ll_rnf .AND. jk <= nk_rnf(ji,jj) )                                          &
338                     &     ztc_f  = ztc_f  - zfact1 * ( rnf_tsc(ji,jj,jn) - rnf_tsc_b(ji,jj,jn) ) & 
339                     &                              * e3t(ji,jj,jk,Kmm) / h_rnf(ji,jj)
340                     !
341                  ! ice shelf
342                  IF( ll_isf ) THEN
343                     ! level fully include in the Losch_2008 ice shelf boundary layer
344                     IF ( jk >= misfkt(ji,jj) .AND. jk < misfkb(ji,jj) )                          &
345                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
346                               &                 * e3t(ji,jj,jk,Kmm) * r1_hisf_tbl (ji,jj)
347                     ! level partially include in Losch_2008 ice shelf boundary layer
348                     IF ( jk == misfkb(ji,jj) )                                                   &
349                        ztc_f  = ztc_f  - zfact1 * ( risf_tsc(ji,jj,jn) - risf_tsc_b(ji,jj,jn) )  &
350                               &                 * e3t(ji,jj,jk,Kmm) * r1_hisf_tbl (ji,jj) * ralpha(ji,jj)
351                  END IF
352                  !
353                  ze3t_f = 1.e0 / ze3t_f
354                  pt(ji,jj,jk,jn,Kmm) = ztc_f * ze3t_f    ! time filtered "now" field
355                  !
356                  IF( ( l_trdtra .and. cdtype == 'TRA' ) .OR. ( l_trdtrc .and. cdtype == 'TRC' ) ) THEN
357                     ztrd_atf(ji,jj,jk,jn) = (ztc_f - ztc_n) * zfact/ze3t_n
358                  ENDIF
359                  !
360               END DO
361            END DO
362         END DO
363         !
364      END DO
365      !
366      IF( ( l_trdtra .AND. cdtype == 'TRA' ) .OR. ( l_trdtrc .AND. cdtype == 'TRC' ) )   THEN
367         IF( l_trdtra .AND. cdtype == 'TRA' ) THEN
368            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_tem, jptra_atf, ztrd_atf(:,:,:,jp_tem) )
369            CALL trd_tra( kt, Kmm, Kaa, cdtype, jp_sal, jptra_atf, ztrd_atf(:,:,:,jp_sal) )
370         ENDIF
371         IF( l_trdtrc .AND. cdtype == 'TRC' ) THEN
372            DO jn = 1, kjpt
373               CALL trd_tra( kt, Kmm, Kaa, cdtype, jn, jptra_atf, ztrd_atf(:,:,:,jn) )
374            END DO
375         ENDIF
376         DEALLOCATE( ztrd_atf )
377      ENDIF
378      !
379   END SUBROUTINE tra_atf_vvl
380
381   !!======================================================================
382END MODULE traatf
Note: See TracBrowser for help on using the repository browser.