source: branches/ignacio/FLEXPART_9.1.8/src/concoutput.f90 @ 18

Last change on this file since 18 was 18, checked in by igpis, 10 years ago

add verbose mode to version 9.1.7.1

File size: 21.9 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine concoutput(itime,outnum,gridtotalunc,wetgridtotalunc, &
23       drygridtotalunc)
24  !                        i     i          o             o
25  !       o
26  !*****************************************************************************
27  !                                                                            *
28  !     Output of the concentration grid and the receptor concentrations.      *
29  !                                                                            *
30  !     Author: A. Stohl                                                       *
31  !                                                                            *
32  !     24 May 1995                                                            *
33  !                                                                            *
34  !     13 April 1999, Major update: if output size is smaller, dump output    *
35  !                    in sparse matrix format; additional output of           *
36  !                    uncertainty                                             *
37  !                                                                            *
38  !     05 April 2000, Major update: output of age classes; output for backward*
39  !                    runs is time spent in grid cell times total mass of     *
40  !                    species.                                                *
41  !                                                                            *
42  !     17 February 2002, Appropriate dimensions for backward and forward runs *
43  !                       are now specified in file par_mod                    *
44  !                                                                            *
45  !     June 2006, write grid in sparse matrix with a single write command     *
46  !                in order to save disk space                                 *
47  !                                                                            *
48  !     2008 new sparse matrix format                                          *
49  !                                                                            *
50  !*****************************************************************************
51  !                                                                            *
52  ! Variables:                                                                 *
53  ! outnum          number of samples                                          *
54  ! ncells          number of cells with non-zero concentrations               *
55  ! sparse          .true. if in sparse matrix format, else .false.            *
56  ! tot_mu          1 for forward, initial mass mixing ration for backw. runs  *
57  !                                                                            *
58  !*****************************************************************************
59
60  use unc_mod
61  use point_mod
62  use outg_mod
63  use par_mod
64  use com_mod
65
66  implicit none
67
68  real(kind=dp) :: jul
69  integer :: itime,i,ix,jy,kz,ks,kp,l,iix,jjy,kzz,nage,jjjjmmdd,ihmmss
70  integer :: sp_count_i,sp_count_r
71  real :: sp_fact
72  real :: outnum,densityoutrecept(maxreceptor),xl,yl
73
74  !real densityoutgrid(0:numxgrid-1,0:numygrid-1,numzgrid),
75  !    +grid(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,maxpointspec_act,
76  !    +    maxageclass)
77  !real wetgrid(0:numxgrid-1,0:numygrid-1,maxspec,maxpointspec_act,
78  !    +       maxageclass)
79  !real drygrid(0:numxgrid-1,0:numygrid-1,maxspec,
80  !    +       maxpointspec_act,maxageclass)
81  !real gridsigma(0:numxgrid-1,0:numygrid-1,numzgrid,maxspec,
82  !    +       maxpointspec_act,maxageclass),
83  !    +     drygridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
84  !    +     maxpointspec_act,maxageclass),
85  !    +     wetgridsigma(0:numxgrid-1,0:numygrid-1,maxspec,
86  !    +     maxpointspec_act,maxageclass)
87  !real factor(0:numxgrid-1,0:numygrid-1,numzgrid)
88  !real sparse_dump_r(numxgrid*numygrid*numzgrid)
89  !integer sparse_dump_i(numxgrid*numygrid*numzgrid)
90
91  !real sparse_dump_u(numxgrid*numygrid*numzgrid)
92  real :: auxgrid(nclassunc),gridtotal,gridsigmatotal,gridtotalunc
93  real :: wetgridtotal,wetgridsigmatotal,wetgridtotalunc
94  real :: drygridtotal,drygridsigmatotal,drygridtotalunc
95  real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
96  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
97  real,parameter :: weightair=28.97
98  logical :: sp_zer
99  character :: adate*8,atime*6
100  character(len=3) :: anspec
101
102
103  ! Determine current calendar date, needed for the file name
104  !**********************************************************
105
106  jul=bdate+real(itime,kind=dp)/86400._dp
107  call caldate(jul,jjjjmmdd,ihmmss)
108  write(adate,'(i8.8)') jjjjmmdd
109  write(atime,'(i6.6)') ihmmss
110  open(unitdates,file=path(2)(1:length(2))//'dates', ACCESS='APPEND')
111  write(unitdates,'(a)') adate//atime
112  close(unitdates) 
113
114  ! For forward simulations, output fields have dimension MAXSPEC,
115  ! for backward simulations, output fields have dimension MAXPOINT.
116  ! Thus, make loops either about nspec, or about numpoint
117  !*****************************************************************
118
119
120  if (ldirect.eq.1) then
121    do ks=1,nspec
122      do kp=1,maxpointspec_act
123        tot_mu(ks,kp)=1
124      end do
125    end do
126  else
127    do ks=1,nspec
128      do kp=1,maxpointspec_act
129        tot_mu(ks,kp)=xmass(kp,ks)
130      end do
131    end do
132  endif
133
134
135  !*******************************************************************
136  ! Compute air density: sufficiently accurate to take it
137  ! from coarse grid at some time
138  ! Determine center altitude of output layer, and interpolate density
139  ! data to that altitude
140  !*******************************************************************
141
142  do kz=1,numzgrid
143    if (kz.eq.1) then
144      halfheight=outheight(1)/2.
145    else
146      halfheight=(outheight(kz)+outheight(kz-1))/2.
147    endif
148    do kzz=2,nz
149      if ((height(kzz-1).lt.halfheight).and. &
150           (height(kzz).gt.halfheight)) goto 46
151    end do
15246   kzz=max(min(kzz,nz),2)
153    dz1=halfheight-height(kzz-1)
154    dz2=height(kzz)-halfheight
155    dz=dz1+dz2
156    do jy=0,numygrid-1
157      do ix=0,numxgrid-1
158        xl=outlon0+real(ix)*dxout
159        yl=outlat0+real(jy)*dyout
160        xl=(xl-xlon0)/dx
161        yl=(yl-ylat0)/dy !v9.1.1
162        iix=max(min(nint(xl),nxmin1),0)
163        jjy=max(min(nint(yl),nymin1),0)
164        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
165             rho(iix,jjy,kzz-1,2)*dz2)/dz
166      end do
167    end do
168  end do
169
170    do i=1,numreceptor
171      xl=xreceptor(i)
172      yl=yreceptor(i)
173      iix=max(min(nint(xl),nxmin1),0)
174      jjy=max(min(nint(yl),nymin1),0)
175      densityoutrecept(i)=rho(iix,jjy,1,2)
176    end do
177
178
179  ! Output is different for forward and backward simulations
180    do kz=1,numzgrid
181      do jy=0,numygrid-1
182        do ix=0,numxgrid-1
183          if (ldirect.eq.1) then
184            factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum
185          else
186            factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
187          endif
188        end do
189      end do
190    end do
191
192  !*********************************************************************
193  ! Determine the standard deviation of the mean concentration or mixing
194  ! ratio (uncertainty of the output) and the dry and wet deposition
195  !*********************************************************************
196
197  gridtotal=0.
198  gridsigmatotal=0.
199  gridtotalunc=0.
200  wetgridtotal=0.
201  wetgridsigmatotal=0.
202  wetgridtotalunc=0.
203  drygridtotal=0.
204  drygridsigmatotal=0.
205  drygridtotalunc=0.
206
207  do ks=1,nspec
208
209  write(anspec,'(i3.3)') ks
210  if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
211    if (ldirect.eq.1) then
212      open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
213           atime//'_'//anspec,form='unformatted')
214    else
215      open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
216           atime//'_'//anspec,form='unformatted')
217    endif
218     write(unitoutgrid) itime
219   endif
220
221  if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
222   open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
223        atime//'_'//anspec,form='unformatted')
224
225    write(unitoutgridppt) itime
226  endif
227
228  do kp=1,maxpointspec_act
229  do nage=1,nageclass
230
231    do jy=0,numygrid-1
232      do ix=0,numxgrid-1
233
234  ! WET DEPOSITION
235        if ((WETDEP).and.(ldirect.gt.0)) then
236            do l=1,nclassunc
237              auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
238            end do
239            call mean(auxgrid,wetgrid(ix,jy), &
240                 wetgridsigma(ix,jy),nclassunc)
241  ! Multiply by number of classes to get total concentration
242            wetgrid(ix,jy)=wetgrid(ix,jy) &
243                 *nclassunc
244            wetgridtotal=wetgridtotal+wetgrid(ix,jy)
245  ! Calculate standard deviation of the mean
246            wetgridsigma(ix,jy)= &
247                 wetgridsigma(ix,jy)* &
248                 sqrt(real(nclassunc))
249            wetgridsigmatotal=wetgridsigmatotal+ &
250                 wetgridsigma(ix,jy)
251        endif
252
253  ! DRY DEPOSITION
254        if ((DRYDEP).and.(ldirect.gt.0)) then
255            do l=1,nclassunc
256              auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
257            end do
258            call mean(auxgrid,drygrid(ix,jy), &
259                 drygridsigma(ix,jy),nclassunc)
260  ! Multiply by number of classes to get total concentration
261            drygrid(ix,jy)=drygrid(ix,jy)* &
262                 nclassunc
263            drygridtotal=drygridtotal+drygrid(ix,jy)
264  ! Calculate standard deviation of the mean
265            drygridsigma(ix,jy)= &
266                 drygridsigma(ix,jy)* &
267                 sqrt(real(nclassunc))
268125         drygridsigmatotal=drygridsigmatotal+ &
269                 drygridsigma(ix,jy)
270        endif
271
272  ! CONCENTRATION OR MIXING RATIO
273        do kz=1,numzgrid
274            do l=1,nclassunc
275              auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
276            end do
277            call mean(auxgrid,grid(ix,jy,kz), &
278                 gridsigma(ix,jy,kz),nclassunc)
279  ! Multiply by number of classes to get total concentration
280            grid(ix,jy,kz)= &
281                 grid(ix,jy,kz)*nclassunc
282            gridtotal=gridtotal+grid(ix,jy,kz)
283  ! Calculate standard deviation of the mean
284            gridsigma(ix,jy,kz)= &
285                 gridsigma(ix,jy,kz)* &
286                 sqrt(real(nclassunc))
287            gridsigmatotal=gridsigmatotal+ &
288                 gridsigma(ix,jy,kz)
289        end do
290      end do
291    end do
292
293
294
295
296  !*******************************************************************
297  ! Generate output: may be in concentration (ng/m3) or in mixing
298  ! ratio (ppt) or both
299  ! Output the position and the values alternated multiplied by
300  ! 1 or -1, first line is number of values, number of positions
301  ! For backward simulations, the unit is seconds, stored in grid_time
302  !*******************************************************************
303
304  ! Concentration output
305  !*********************
306  if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
307
308  ! Wet deposition
309         sp_count_i=0
310         sp_count_r=0
311         sp_fact=-1.
312         sp_zer=.true.
313         if ((ldirect.eq.1).and.(WETDEP)) then
314         do jy=0,numygrid-1
315            do ix=0,numxgrid-1
316  !oncentraion greater zero
317              if (wetgrid(ix,jy).gt.smallnum) then
318                 if (sp_zer.eqv..true.) then ! first non zero value
319                    sp_count_i=sp_count_i+1
320                    sparse_dump_i(sp_count_i)=ix+jy*numxgrid
321                    sp_zer=.false.
322                    sp_fact=sp_fact*(-1.)
323                 endif
324                 sp_count_r=sp_count_r+1
325                 sparse_dump_r(sp_count_r)= &
326                      sp_fact*1.e12*wetgrid(ix,jy)/area(ix,jy)
327  !                sparse_dump_u(sp_count_r)=
328  !+                1.e12*wetgridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
329              else ! concentration is zero
330                  sp_zer=.true.
331              endif
332            end do
333         end do
334         else
335            sp_count_i=0
336            sp_count_r=0
337         endif
338         write(unitoutgrid) sp_count_i
339         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
340         write(unitoutgrid) sp_count_r
341         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
342  !       write(unitoutgrid) sp_count_u
343  !       write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
344
345  ! Dry deposition
346         sp_count_i=0
347         sp_count_r=0
348         sp_fact=-1.
349         sp_zer=.true.
350         if ((ldirect.eq.1).and.(DRYDEP)) then
351          do jy=0,numygrid-1
352            do ix=0,numxgrid-1
353              if (drygrid(ix,jy).gt.smallnum) then
354                 if (sp_zer.eqv..true.) then ! first non zero value
355                    sp_count_i=sp_count_i+1
356                    sparse_dump_i(sp_count_i)=ix+jy*numxgrid
357                    sp_zer=.false.
358                    sp_fact=sp_fact*(-1.)
359                 endif
360                 sp_count_r=sp_count_r+1
361                 sparse_dump_r(sp_count_r)= &
362                      sp_fact* &
363                      1.e12*drygrid(ix,jy)/area(ix,jy)
364  !                sparse_dump_u(sp_count_r)=
365  !+                1.e12*drygridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
366              else ! concentration is zero
367                  sp_zer=.true.
368              endif
369            end do
370          end do
371         else
372            sp_count_i=0
373            sp_count_r=0
374         endif
375         write(unitoutgrid) sp_count_i
376         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
377         write(unitoutgrid) sp_count_r
378         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
379  !       write(*,*) sp_count_u
380  !       write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
381
382
383
384  ! Concentrations
385         sp_count_i=0
386         sp_count_r=0
387         sp_fact=-1.
388         sp_zer=.true.
389          do kz=1,numzgrid
390            do jy=0,numygrid-1
391              do ix=0,numxgrid-1
392                if (grid(ix,jy,kz).gt.smallnum) then
393                  if (sp_zer.eqv..true.) then ! first non zero value
394                    sp_count_i=sp_count_i+1
395                    sparse_dump_i(sp_count_i)= &
396                         ix+jy*numxgrid+kz*numxgrid*numygrid
397                    sp_zer=.false.
398                    sp_fact=sp_fact*(-1.)
399                   endif
400                   sp_count_r=sp_count_r+1
401                   sparse_dump_r(sp_count_r)= &
402                        sp_fact* &
403                        grid(ix,jy,kz)* &
404                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
405  !                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
406  !    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
407  !                sparse_dump_u(sp_count_r)=
408  !+               ,gridsigma(ix,jy,kz,ks,kp,nage)*
409  !+               factor(ix,jy,kz)/tot_mu(ks,kp)
410              else ! concentration is zero
411                  sp_zer=.true.
412              endif
413              end do
414            end do
415          end do
416         write(unitoutgrid) sp_count_i
417         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
418         write(unitoutgrid) sp_count_r
419         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
420  !       write(unitoutgrid) sp_count_u
421  !       write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
422
423
424
425    endif !  concentration output
426
427  ! Mixing ratio output
428  !********************
429
430  if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
431
432  ! Wet deposition
433         sp_count_i=0
434         sp_count_r=0
435         sp_fact=-1.
436         sp_zer=.true.
437         if ((ldirect.eq.1).and.(WETDEP)) then
438          do jy=0,numygrid-1
439            do ix=0,numxgrid-1
440                if (wetgrid(ix,jy).gt.smallnum) then
441                  if (sp_zer.eqv..true.) then ! first non zero value
442                    sp_count_i=sp_count_i+1
443                    sparse_dump_i(sp_count_i)= &
444                         ix+jy*numxgrid
445                    sp_zer=.false.
446                    sp_fact=sp_fact*(-1.)
447                 endif
448                 sp_count_r=sp_count_r+1
449                 sparse_dump_r(sp_count_r)= &
450                      sp_fact* &
451                      1.e12*wetgrid(ix,jy)/area(ix,jy)
452  !                sparse_dump_u(sp_count_r)=
453  !    +            ,1.e12*wetgridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
454              else ! concentration is zero
455                  sp_zer=.true.
456              endif
457            end do
458          end do
459         else
460           sp_count_i=0
461           sp_count_r=0
462         endif
463         write(unitoutgridppt) sp_count_i
464         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
465         write(unitoutgridppt) sp_count_r
466         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
467  !       write(unitoutgridppt) sp_count_u
468  !       write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
469
470
471  ! Dry deposition
472         sp_count_i=0
473         sp_count_r=0
474         sp_fact=-1.
475         sp_zer=.true.
476         if ((ldirect.eq.1).and.(DRYDEP)) then
477          do jy=0,numygrid-1
478            do ix=0,numxgrid-1
479                if (drygrid(ix,jy).gt.smallnum) then
480                  if (sp_zer.eqv..true.) then ! first non zero value
481                    sp_count_i=sp_count_i+1
482                    sparse_dump_i(sp_count_i)= &
483                         ix+jy*numxgrid
484                    sp_zer=.false.
485                    sp_fact=sp_fact*(-1)
486                 endif
487                 sp_count_r=sp_count_r+1
488                 sparse_dump_r(sp_count_r)= &
489                      sp_fact* &
490                      1.e12*drygrid(ix,jy)/area(ix,jy)
491  !                sparse_dump_u(sp_count_r)=
492  !    +            ,1.e12*drygridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
493              else ! concentration is zero
494                  sp_zer=.true.
495              endif
496            end do
497          end do
498         else
499           sp_count_i=0
500           sp_count_r=0
501         endif
502         write(unitoutgridppt) sp_count_i
503         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
504         write(unitoutgridppt) sp_count_r
505         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
506  !       write(unitoutgridppt) sp_count_u
507  !       write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
508
509
510  ! Mixing ratios
511         sp_count_i=0
512         sp_count_r=0
513         sp_fact=-1.
514         sp_zer=.true.
515          do kz=1,numzgrid
516            do jy=0,numygrid-1
517              do ix=0,numxgrid-1
518                if (grid(ix,jy,kz).gt.smallnum) then
519                  if (sp_zer.eqv..true.) then ! first non zero value
520                    sp_count_i=sp_count_i+1
521                    sparse_dump_i(sp_count_i)= &
522                         ix+jy*numxgrid+kz*numxgrid*numygrid
523                    sp_zer=.false.
524                    sp_fact=sp_fact*(-1.)
525                 endif
526                 sp_count_r=sp_count_r+1
527                 sparse_dump_r(sp_count_r)= &
528                      sp_fact* &
529                      1.e12*grid(ix,jy,kz) &
530                      /volume(ix,jy,kz)/outnum* &
531                      weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
532  !                sparse_dump_u(sp_count_r)=
533  !+              ,1.e12*gridsigma(ix,jy,kz,ks,kp,nage)/volume(ix,jy,kz)/
534  !+              outnum*weightair/weightmolar(ks)/
535  !+              densityoutgrid(ix,jy,kz)
536              else ! concentration is zero
537                  sp_zer=.true.
538              endif
539              end do
540            end do
541          end do
542         write(unitoutgridppt) sp_count_i
543         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
544         write(unitoutgridppt) sp_count_r
545         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
546  !       write(unitoutgridppt) sp_count_u
547  !       write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
548
549      endif ! output for ppt
550
551  end do
552  end do
553
554    close(unitoutgridppt)
555    close(unitoutgrid)
556
557  end do
558
559  if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
560  if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
561       wetgridtotal
562  if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
563       drygridtotal
564
565  ! Dump of receptor concentrations
566
567    if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3)  ) then
568      write(unitoutreceptppt) itime
569      do ks=1,nspec
570        write(unitoutreceptppt) (1.e12*creceptor(i,ks)/outnum* &
571             weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)
572      end do
573    endif
574
575  ! Dump of receptor concentrations
576
577    if (numreceptor.gt.0) then
578      write(unitoutrecept) itime
579      do ks=1,nspec
580        write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum, &
581             i=1,numreceptor)
582      end do
583    endif
584
585
586
587  ! Reinitialization of grid
588  !*************************
589
590  do ks=1,nspec
591  do kp=1,maxpointspec_act
592    do i=1,numreceptor
593      creceptor(i,ks)=0.
594    end do
595    do jy=0,numygrid-1
596      do ix=0,numxgrid-1
597        do l=1,nclassunc
598          do nage=1,nageclass
599            do kz=1,numzgrid
600              gridunc(ix,jy,kz,ks,kp,l,nage)=0.
601            end do
602          end do
603        end do
604      end do
605    end do
606  end do
607  end do
608
609
610end subroutine concoutput
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG