source: branches/UKMO/dev_r7750_GO6_package_oasis_timers/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_flt.F90 @ 9252

Last change on this file since 9252 was 9252, checked in by andmirek, 3 years ago

#1978 timers for nn_timing 3

File size: 20.3 KB
Line 
1MODULE dynspg_flt
2   !!======================================================================
3   !!                   ***  MODULE  dynspg_flt  ***
4   !! Ocean dynamics:  surface pressure gradient trend
5   !!======================================================================
6   !! History    OPA  !  1998-05  (G. Roullet)  free surface
7   !!                 !  1998-10  (G. Madec, M. Imbard)  release 8.2
8   !!   NEMO     O.1  !  2002-08  (G. Madec)  F90: Free form and module
9   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
10   !!            1.0  !  2004-08  (C. Talandier) New trends organization
11   !!             -   !  2005-11  (V. Garnier) Surface pressure gradient organization
12   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
13   !!             -   !  2006-08  (J.Chanut, A.Sellar) Calls to BDY routines.
14   !!            3.2  !  2009-03  (G. Madec, M. Leclair, R. Benshila) introduce sshwzv module
15   !!            3.7  !  2014-04  (F. Roquet, G. Madec)  add some trends diag
16   !!----------------------------------------------------------------------
17#if defined key_dynspg_flt   ||   defined key_esopa 
18   !!----------------------------------------------------------------------
19   !!   'key_dynspg_flt'                              filtered free surface
20   !!----------------------------------------------------------------------
21   !!   dyn_spg_flt  : update the momentum trend with the surface pressure gradient in the filtered free surface case
22   !!   flt_rst      : read/write the time-splitting restart fields in the ocean restart file
23   !!----------------------------------------------------------------------
24   USE oce             ! ocean dynamics and tracers
25   USE dom_oce         ! ocean space and time domain
26   USE zdf_oce         ! ocean vertical physics
27   USE sbc_oce         ! surface boundary condition: ocean
28   USE bdy_oce         ! Lateral open boundary condition
29   USE sol_oce         ! ocean elliptic solver
30   USE phycst          ! physical constants
31   USE domvvl          ! variable volume
32   USE dynadv          ! advection
33   USE solmat          ! matrix construction for elliptic solvers
34   USE solpcg          ! preconditionned conjugate gradient solver
35   USE solsor          ! Successive Over-relaxation solver
36   USE bdydyn          ! ocean open boundary condition on dynamics
37   USE bdyvol          ! ocean open boundary condition (bdy_vol routine)
38   USE cla             ! cross land advection
39   USE trd_oce         ! trends: ocean variables
40   USE trddyn          ! trend manager: dynamics
41   !
42   USE in_out_manager  ! I/O manager
43   USE lib_mpp         ! distributed memory computing library
44   USE wrk_nemo        ! Memory Allocation
45   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
46   USE prtctl          ! Print control
47   USE iom
48   USE lib_fortran
49   USE timing          ! Timing
50#if defined key_agrif
51   USE agrif_opa_interp
52#endif
53
54   IMPLICIT NONE
55   PRIVATE
56
57   PUBLIC   dyn_spg_flt  ! routine called by step.F90
58   PUBLIC   flt_rst      ! routine called by istate.F90
59
60   !! * Substitutions
61#  include "domzgr_substitute.h90"
62#  include "vectopt_loop_substitute.h90"
63   !!----------------------------------------------------------------------
64   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
65   !! $Id$
66   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
67   !!----------------------------------------------------------------------
68CONTAINS
69
70   SUBROUTINE dyn_spg_flt( kt, kindic )
71      !!----------------------------------------------------------------------
72      !!                  ***  routine dyn_spg_flt  ***
73      !!
74      !! ** Purpose :   Compute the now trend due to the surface pressure
75      !!      gradient in case of filtered free surface formulation  and add
76      !!      it to the general trend of momentum equation.
77      !!
78      !! ** Method  :   Filtered free surface formulation. The surface
79      !!      pressure gradient is given by:
80      !!         spgu = 1/rau0 d/dx(ps) =  1/e1u di( sshn + btda )
81      !!         spgv = 1/rau0 d/dy(ps) =  1/e2v dj( sshn + btda )
82      !!      where sshn is the free surface elevation and btda is the after
83      !!      time derivative of the free surface elevation
84      !!       -1- evaluate the surface presure trend (including the addi-
85      !!      tional force) in three steps:
86      !!        a- compute the right hand side of the elliptic equation:
87      !!            gcb = 1/(e1t e2t) [ di(e2u spgu) + dj(e1v spgv) ]
88      !!         where (spgu,spgv) are given by:
89      !!            spgu = vertical sum[ e3u (ub+ 2 rdt ua ) ]
90      !!                 - grav 2 rdt hu /e1u di[sshn + (emp-rnf)]
91      !!            spgv = vertical sum[ e3v (vb+ 2 rdt va) ]
92      !!                 - grav 2 rdt hv /e2v dj[sshn + (emp-rnf)]
93      !!         and define the first guess from previous computation :
94      !!            zbtd = btda
95      !!            btda = 2 zbtd - btdb
96      !!            btdb = zbtd
97      !!        b- compute the relative accuracy to be reached by the
98      !!         iterative solver
99      !!        c- apply the solver by a call to sol... routine
100      !!       -2- compute and add the free surface pressure gradient inclu-
101      !!      ding the additional force used to stabilize the equation.
102      !!
103      !! ** Action : - Update (ua,va) with the surf. pressure gradient trend
104      !!
105      !! References : Roullet and Madec, JGR, 2000.
106      !!---------------------------------------------------------------------
107      INTEGER, INTENT(in   ) ::   kt       ! ocean time-step index
108      INTEGER, INTENT(  out) ::   kindic   ! solver convergence flag (<0 if not converge)
109      !
110      INTEGER  ::   ji, jj, jk   ! dummy loop indices
111      REAL(wp) ::   z2dt, z2dtg, zgcb, zbtd, ztdgu, ztdgv   ! local scalars
112      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdu, ztrdv
113      REAL(wp), POINTER, DIMENSION(:,:)   ::  zpw
114      !!----------------------------------------------------------------------
115      !
116      IF( nn_timing == 1 )  CALL timing_start('dyn_spg_flt')
117      !
118      IF( kt == nit000 ) THEN
119         IF(lwp) WRITE(numout,*)
120         IF(lwp) WRITE(numout,*) 'dyn_spg_flt : surface pressure gradient trend'
121         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   (free surface constant volume case)'
122       
123         ! set to zero free surface specific arrays
124         spgu(:,:) = 0._wp                     ! surface pressure gradient (i-direction)
125         spgv(:,:) = 0._wp                     ! surface pressure gradient (j-direction)
126
127         ! read filtered free surface arrays in restart file
128         ! when using agrif, sshn, gcx have to be read in istate
129         IF(.NOT. lk_agrif)   CALL flt_rst( nit000, 'READ' )      ! read or initialize the following fields:
130         !                                                        ! gcx, gcxb
131      ENDIF
132
133      ! Local constant initialization
134      z2dt = 2. * rdt                                             ! time step: leap-frog
135      IF( neuler == 0 .AND. kt == nit000   )   z2dt = rdt         ! time step: Euler if restart from rest
136      IF( neuler == 0 .AND. kt == nit000+1 )   CALL sol_mat( kt )
137      z2dtg  = grav * z2dt
138
139      ! Evaluate the masked next velocity (effect of the additional force not included)
140      ! --------------------------------- 
141      IF( lk_vvl ) THEN          ! variable volume  (surface pressure gradient already included in dyn_hpg)
142         !
143         IF( ln_dynadv_vec ) THEN      ! vector form : applied on velocity
144            DO jk = 1, jpkm1
145               DO jj = 2, jpjm1
146                  DO ji = fs_2, fs_jpim1   ! vector opt.
147                     ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ua(ji,jj,jk)  ) * umask(ji,jj,jk)
148                     va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * va(ji,jj,jk)  ) * vmask(ji,jj,jk)
149                  END DO
150               END DO
151            END DO
152            !
153         ELSE                          ! flux form : applied on thickness weighted velocity
154            DO jk = 1, jpkm1
155               DO jj = 2, jpjm1
156                  DO ji = fs_2, fs_jpim1   ! vector opt.
157                     ua(ji,jj,jk) = (        ub(ji,jj,jk) * fse3u_b(ji,jj,jk)      &
158                        &           + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk)  )   &
159                        &         / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
160                     va(ji,jj,jk) = (        vb(ji,jj,jk) * fse3v_b(ji,jj,jk)      &
161                        &           + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk)  )   &
162                        &         / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
163                 END DO
164               END DO
165            END DO
166            !
167         ENDIF
168        IF( l_trddyn )   THEN                      ! Put here so code doesn't crash when doing KE trend but needs to be done properly
169            CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )
170         ENDIF
171         !
172      ELSE                       ! fixed volume  (add the surface pressure gradient + unweighted time stepping)
173         !
174         DO jj = 2, jpjm1              ! Surface pressure gradient (now)
175            DO ji = fs_2, fs_jpim1   ! vector opt.
176               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
177               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
178            END DO
179         END DO
180         DO jk = 1, jpkm1              ! unweighted time stepping
181            DO jj = 2, jpjm1
182               DO ji = fs_2, fs_jpim1   ! vector opt.
183                  ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) )  ) * umask(ji,jj,jk)
184                  va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) )  ) * vmask(ji,jj,jk)
185               END DO
186            END DO
187         END DO
188         !
189         IF( l_trddyn )   THEN                      ! temporary save of spg trends
190            CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv )
191            DO jk = 1, jpkm1              ! unweighted time stepping
192               DO jj = 2, jpjm1
193                  DO ji = fs_2, fs_jpim1   ! vector opt.
194                     ztrdu(ji,jj,jk) = spgu(ji,jj) * umask(ji,jj,jk)
195                     ztrdv(ji,jj,jk) = spgv(ji,jj) * vmask(ji,jj,jk)
196                  END DO
197               END DO
198            END DO
199            CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgexp, kt )
200         ENDIF
201         !
202      ENDIF
203
204#if defined key_bdy
205      IF( lk_bdy ) CALL bdy_dyn( kt )   ! Update velocities on each open boundary
206      IF( lk_bdy ) CALL bdy_vol( kt )   ! Correction of the barotropic component velocity to control the volume of the system
207#endif
208#if defined key_agrif
209      CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces
210#endif
211      IF( nn_cla == 1 .AND. cp_cfg == 'orca' .AND. jp_cfg == 2 )   CALL cla_dynspg( kt )      ! Cross Land Advection (update (ua,va))
212
213      ! compute the next vertically averaged velocity (effect of the additional force not included)
214      ! ---------------------------------------------
215      DO jj = 2, jpjm1
216         DO ji = fs_2, fs_jpim1   ! vector opt.
217            spgu(ji,jj) = fse3u_a(ji,jj,1) * ua(ji,jj,1)
218            spgv(ji,jj) = fse3v_a(ji,jj,1) * va(ji,jj,1)
219         END DO
220      END DO
221      DO jk = 2, jpkm1                     ! vertical sum
222         DO jj = 2, jpjm1
223            DO ji = fs_2, fs_jpim1   ! vector opt.
224               spgu(ji,jj) = spgu(ji,jj) + fse3u_a(ji,jj,jk) * ua(ji,jj,jk)
225               spgv(ji,jj) = spgv(ji,jj) + fse3v_a(ji,jj,jk) * va(ji,jj,jk)
226            END DO
227         END DO
228      END DO
229
230      DO jj = 2, jpjm1                     ! transport: multiplied by the horizontal scale factor
231         DO ji = fs_2, fs_jpim1   ! vector opt.
232            spgu(ji,jj) = spgu(ji,jj) * e2u(ji,jj)
233            spgv(ji,jj) = spgv(ji,jj) * e1v(ji,jj)
234         END DO
235      END DO
236      CALL lbc_lnk( spgu, 'U', -1. )       ! lateral boundary conditions
237      CALL lbc_lnk( spgv, 'V', -1. )
238
239      IF( lk_vvl ) CALL sol_mat( kt )      ! build the matrix at kt (vvl case only)
240
241      ! Right hand side of the elliptic equation and first guess
242      ! --------------------------------------------------------
243      DO jj = 2, jpjm1
244         DO ji = fs_2, fs_jpim1   ! vector opt.
245            ! Divergence of the after vertically averaged velocity
246            zgcb =  spgu(ji,jj) - spgu(ji-1,jj)   &
247                  + spgv(ji,jj) - spgv(ji,jj-1)
248            gcb(ji,jj) = gcdprc(ji,jj) * zgcb
249            ! First guess of the after barotropic transport divergence
250            zbtd = gcx(ji,jj)
251            gcx (ji,jj) = 2. * zbtd   - gcxb(ji,jj)
252            gcxb(ji,jj) =      zbtd
253         END DO
254      END DO
255      ! applied the lateral boundary conditions
256      IF( nn_solv == 2 .AND. MAX( jpr2di, jpr2dj ) > 0 )   CALL lbc_lnk_e( gcb, c_solver_pt, 1., jpr2di, jpr2dj )   
257
258#if defined key_agrif
259      IF( .NOT. AGRIF_ROOT() ) THEN
260         ! add contribution of gradient of after barotropic transport divergence
261         IF( nbondi == -1 .OR. nbondi == 2 )   gcb(3     ,:) =   &
262            &    gcb(3     ,:) - z2dtg * z2dt * laplacu(2     ,:) * gcdprc(3     ,:) * hu(2     ,:) * e2u(2     ,:)
263         IF( nbondi ==  1 .OR. nbondi == 2 )   gcb(nlci-2,:) =   &
264            &    gcb(nlci-2,:) + z2dtg * z2dt * laplacu(nlci-2,:) * gcdprc(nlci-2,:) * hu(nlci-2,:) * e2u(nlci-2,:)
265         IF( nbondj == -1 .OR. nbondj == 2 )   gcb(:     ,3) =   &
266            &    gcb(:,3     ) - z2dtg * z2dt * laplacv(:,2     ) * gcdprc(:,3     ) * hv(:,2     ) * e1v(:,2     )
267         IF( nbondj ==  1 .OR. nbondj == 2 )   gcb(:,nlcj-2) =   &
268            &    gcb(:,nlcj-2) + z2dtg * z2dt * laplacv(:,nlcj-2) * gcdprc(:,nlcj-2) * hv(:,nlcj-2) * e1v(:,nlcj-2)
269      ENDIF
270#endif
271
272
273      ! Relative precision (computation on one processor)
274      ! ------------------
275      rnorme =0.e0
276      rnorme = GLOB_SUM( gcb(1:jpi,1:jpj) * gcdmat(1:jpi,1:jpj) * gcb(1:jpi,1:jpj) * bmask(:,:) )
277
278      epsr = eps * eps * rnorme
279      ncut = 0
280      ! if rnorme is 0, the solution is 0, the solver is not called
281      IF( rnorme == 0._wp ) THEN
282         gcx(:,:) = 0._wp
283         res   = 0._wp
284         niter = 0
285         ncut  = 999
286      ENDIF
287
288      ! Evaluate the next transport divergence
289      ! --------------------------------------
290      !    Iterarive solver for the elliptic equation (except IF sol.=0)
291      !    (output in gcx with boundary conditions applied)
292      kindic = 0
293      IF( ncut == 0 ) THEN
294         IF    ( nn_solv == 1 ) THEN   ;   CALL sol_pcg( kindic )      ! diagonal preconditioned conjuguate gradient
295         ELSEIF( nn_solv == 2 ) THEN   ;   CALL sol_sor( kindic )      ! successive-over-relaxation
296         ENDIF
297      ENDIF
298
299      ! Transport divergence gradient multiplied by z2dt
300      ! --------------------------------------------====
301      DO jj = 2, jpjm1
302         DO ji = fs_2, fs_jpim1   ! vector opt.
303            ! trend of Transport divergence gradient
304            ztdgu = z2dtg * (gcx(ji+1,jj  ) - gcx(ji,jj) ) / e1u(ji,jj)
305            ztdgv = z2dtg * (gcx(ji  ,jj+1) - gcx(ji,jj) ) / e2v(ji,jj)
306            ! multiplied by z2dt
307#if defined key_bdy
308            IF(lk_bdy) THEN
309            ! caution : grad D = 0 along open boundaries
310               spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
311               spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)
312            ELSE
313               spgu(ji,jj) = z2dt * ztdgu
314               spgv(ji,jj) = z2dt * ztdgv
315            ENDIF
316#else
317            spgu(ji,jj) = z2dt * ztdgu
318            spgv(ji,jj) = z2dt * ztdgv
319#endif
320         END DO
321      END DO
322
323#if defined key_agrif     
324      IF( .NOT. Agrif_Root() ) THEN
325         ! caution : grad D (fine) = grad D (coarse) at coarse/fine interface
326         IF( nbondi == -1 .OR. nbondi == 2 ) spgu(2     ,:) = z2dtg * z2dt * laplacu(2     ,:) * umask(2     ,:,1)
327         IF( nbondi ==  1 .OR. nbondi == 2 ) spgu(nlci-2,:) = z2dtg * z2dt * laplacu(nlci-2,:) * umask(nlci-2,:,1)
328         IF( nbondj == -1 .OR. nbondj == 2 ) spgv(:,2     ) = z2dtg * z2dt * laplacv(:,2     ) * vmask(:     ,2,1)
329         IF( nbondj ==  1 .OR. nbondj == 2 ) spgv(:,nlcj-2) = z2dtg * z2dt * laplacv(:,nlcj-2) * vmask(:,nlcj-2,1)
330      ENDIF
331#endif     
332
333      IF( l_trddyn )   THEN                     
334         ztrdu(:,:,:) = ua(:,:,:)                 ! save the after velocity before the filtered SPG
335         ztrdv(:,:,:) = va(:,:,:)
336         !
337         CALL wrk_alloc( jpi, jpj, zpw )
338         !
339         zpw(:,:) = - z2dt * gcx(:,:)
340         CALL iom_put( "ssh_flt" , zpw )          ! output equivalent ssh modification due to implicit filter
341         !
342         !                                        ! save surface pressure flux: -pw at z=0
343         zpw(:,:) = - rau0 * grav * sshn(:,:) * wn(:,:,1) * tmask(:,:,1)
344         CALL iom_put( "pw0_exp" , zpw )
345         zpw(:,:) = wn(:,:,1)
346         CALL iom_put( "w0" , zpw )
347         zpw(:,:) =  rau0 * z2dtg * gcx(:,:) * wn(:,:,1) * tmask(:,:,1)
348         CALL iom_put( "pw0_flt" , zpw )
349         !
350         CALL wrk_dealloc( jpi, jpj, zpw ) 
351         !                                   
352      ENDIF
353     
354      ! Add the trends multiplied by z2dt to the after velocity
355      ! -------------------------------------------------------
356      !     ( c a u t i o n : (ua,va) here are the after velocity not the
357      !                       trend, the leap-frog time stepping will not
358      !                       be done in dynnxt.F90 routine)
359      DO jk = 1, jpkm1
360         DO jj = 2, jpjm1
361            DO ji = fs_2, fs_jpim1   ! vector opt.
362               ua(ji,jj,jk) = ( ua(ji,jj,jk) + spgu(ji,jj) ) * umask(ji,jj,jk)
363               va(ji,jj,jk) = ( va(ji,jj,jk) + spgv(ji,jj) ) * vmask(ji,jj,jk)
364            END DO
365         END DO
366      END DO
367
368      IF( l_trddyn )   THEN                      ! save the explicit SPG trends for further diagnostics
369         ztrdu(:,:,:) = ( ua(:,:,:) - ztrdu(:,:,:) ) / z2dt
370         ztrdv(:,:,:) = ( va(:,:,:) - ztrdv(:,:,:) ) / z2dt
371         CALL trd_dyn( ztrdu, ztrdv, jpdyn_spgflt, kt )
372         !
373         CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 
374      ENDIF
375
376      IF( lrst_oce )   CALL flt_rst( kt, 'WRITE' )      ! write filtered free surface arrays in restart file
377      !
378      IF( nn_timing == 1 )   CALL timing_stop('dyn_spg_flt')
379      !
380   END SUBROUTINE dyn_spg_flt
381
382
383   SUBROUTINE flt_rst( kt, cdrw )
384      !!---------------------------------------------------------------------
385      !!                   ***  ROUTINE ts_rst  ***
386      !!
387      !! ** Purpose : Read or write filtered free surface arrays in restart file
388      !!----------------------------------------------------------------------
389      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
390      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
391      !!----------------------------------------------------------------------
392      !
393      IF( TRIM(cdrw) == 'READ' ) THEN
394         IF( iom_varid( numror, 'gcx', ldstop = .FALSE. ) > 0 ) THEN
395! Caution : extra-hallow
396! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
397            CALL iom_get( numror, jpdom_autoglo, 'gcx' , gcx (1:jpi,1:jpj) )
398            CALL iom_get( numror, jpdom_autoglo, 'gcxb', gcxb(1:jpi,1:jpj) )
399            IF( neuler == 0 )   gcxb(:,:) = gcx (:,:)
400         ELSE
401            gcx (:,:) = 0.e0
402            gcxb(:,:) = 0.e0
403         ENDIF
404      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
405! Caution : extra-hallow
406! gcx and gcxb are defined as: DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj)
407         IF(nn_timing == 3)  CALL timing_start('rst_put')
408         CALL iom_rstput( kt, nitrst, numrow, 'gcx' , gcx (1:jpi,1:jpj) )
409         CALL iom_rstput( kt, nitrst, numrow, 'gcxb', gcxb(1:jpi,1:jpj) )
410         IF(nn_timing == 3)  CALL timing_stop('rst_put')
411      ENDIF
412      !
413   END SUBROUTINE flt_rst
414
415#else
416   !!----------------------------------------------------------------------
417   !!   Default case :   Empty module   No standart free surface cst volume
418   !!----------------------------------------------------------------------
419CONTAINS
420   SUBROUTINE dyn_spg_flt( kt, kindic )       ! Empty routine
421      WRITE(*,*) 'dyn_spg_flt: You should not have seen this print! error?', kt, kindic
422   END SUBROUTINE dyn_spg_flt
423   SUBROUTINE flt_rst    ( kt, cdrw )         ! Empty routine
424      INTEGER         , INTENT(in) ::   kt         ! ocean time-step
425      CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
426      WRITE(*,*) 'flt_rst: You should not have seen this print! error?', kt, cdrw
427   END SUBROUTINE flt_rst
428#endif
429   
430   !!======================================================================
431END MODULE dynspg_flt
Note: See TracBrowser for help on using the repository browser.