source: branches/petra/src/concoutput_surf.f90 @ 36

Last change on this file since 36 was 36, checked in by pesei, 9 years ago

Implement switch for incremental deposition, see ticket:113 and many small changes, see changelog.txt

File size: 26.6 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_surf(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  if (verbosity.eq.1) then
104     print*,'inside concoutput_surf '
105     CALL SYSTEM_CLOCK(count_clock)
106     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
107  endif
108
109  ! Determine current calendar date, needed for the file name
110  !**********************************************************
111
112  jul=bdate+real(itime,kind=dp)/86400._dp
113  call caldate(jul,jjjjmmdd,ihmmss)
114  write(adate,'(i8.8)') jjjjmmdd
115  write(atime,'(i6.6)') ihmmss
116  !write(unitdates,'(a)') adate//atime
117
118    open(unitdates,file=path(2)(1:length(2))//'dates', position='append')
119      write(unitdates,'(a)') adate//atime
120    close(unitdates)
121
122  ! For forward simulations, output fields have dimension MAXSPEC,
123  ! for backward simulations, output fields have dimension MAXPOINT.
124  ! Thus, make loops either about nspec, or about numpoint
125  !*****************************************************************
126
127
128  if (ldirect.eq.1) then
129    do ks=1,nspec
130      do kp=1,maxpointspec_act
131        tot_mu(ks,kp)=1
132      end do
133    end do
134  else
135    do ks=1,nspec
136      do kp=1,maxpointspec_act
137        tot_mu(ks,kp)=xmass(kp,ks)
138      end do
139    end do
140  endif
141
142
143  if (verbosity.eq.1) then
144     print*,'concoutput_surf 2'
145     CALL SYSTEM_CLOCK(count_clock)
146     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
147  endif
148
149  !*******************************************************************
150  ! Compute air density: sufficiently accurate to take it
151  ! from coarse grid at some time
152  ! Determine center altitude of output layer, and interpolate density
153  ! data to that altitude
154  !*******************************************************************
155
156  do kz=1,numzgrid
157    if (kz.eq.1) then
158      halfheight=outheight(1)/2.
159    else
160      halfheight=(outheight(kz)+outheight(kz-1))/2.
161    endif
162    do kzz=2,nz
163      if ((height(kzz-1).lt.halfheight).and. &
164           (height(kzz).gt.halfheight)) goto 46
165    end do
16646   kzz=max(min(kzz,nz),2)
167    dz1=halfheight-height(kzz-1)
168    dz2=height(kzz)-halfheight
169    dz=dz1+dz2
170    do jy=0,numygrid-1
171      do ix=0,numxgrid-1
172        xl=outlon0+real(ix)*dxout
173        yl=outlat0+real(jy)*dyout
174        xl=(xl-xlon0)/dx
175        yl=(yl-ylat0)/dx
176        iix=max(min(nint(xl),nxmin1),0)
177        jjy=max(min(nint(yl),nymin1),0)
178        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
179             rho(iix,jjy,kzz-1,2)*dz2)/dz
180      end do
181    end do
182  end do
183
184    do i=1,numreceptor
185      xl=xreceptor(i)
186      yl=yreceptor(i)
187      iix=max(min(nint(xl),nxmin1),0)
188      jjy=max(min(nint(yl),nymin1),0)
189      densityoutrecept(i)=rho(iix,jjy,1,2)
190    end do
191
192
193  ! Output is different for forward and backward simulations
194    do kz=1,numzgrid
195      do jy=0,numygrid-1
196        do ix=0,numxgrid-1
197          if (ldirect.eq.1) then
198            factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum
199          else
200            factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
201          endif
202        end do
203      end do
204    end do
205
206  !*********************************************************************
207  ! Determine the standard deviation of the mean concentration or mixing
208  ! ratio (uncertainty of the output) and the dry and wet deposition
209  !*********************************************************************
210
211  if (verbosity.eq.1) then
212     print*,'concoutput_surf 3 (sd)'
213     CALL SYSTEM_CLOCK(count_clock)
214     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
215  endif
216  gridtotal=0.
217  gridsigmatotal=0.
218  gridtotalunc=0.
219  wetgridtotal=0.
220  wetgridsigmatotal=0.
221  wetgridtotalunc=0.
222  drygridtotal=0.
223  drygridsigmatotal=0.
224  drygridtotalunc=0.
225
226  do ks=1,nspec
227
228  write(anspec,'(i3.3)') ks
229  if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
230    if (ldirect.eq.1) then
231      open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_'//adate// &
232           atime//'_'//anspec,form='unformatted')
233    else
234      open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_'//adate// &
235           atime//'_'//anspec,form='unformatted')
236    endif
237     write(unitoutgrid) itime
238   endif
239
240  if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
241   open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_'//adate// &
242        atime//'_'//anspec,form='unformatted')
243
244    write(unitoutgridppt) itime
245  endif
246
247  do kp=1,maxpointspec_act
248  do nage=1,nageclass
249
250    do jy=0,numygrid-1
251      do ix=0,numxgrid-1
252
253  ! WET DEPOSITION
254        if ((WETDEP).and.(ldirect.gt.0)) then
255            do l=1,nclassunc
256              auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
257            end do
258            call mean(auxgrid,wetgrid(ix,jy), &
259                 wetgridsigma(ix,jy),nclassunc)
260  ! Multiply by number of classes to get total concentration
261            wetgrid(ix,jy)=wetgrid(ix,jy) &
262                 *nclassunc
263            wetgridtotal=wetgridtotal+wetgrid(ix,jy)
264  ! Calculate standard deviation of the mean
265            wetgridsigma(ix,jy)= &
266                 wetgridsigma(ix,jy)* &
267                 sqrt(real(nclassunc))
268            wetgridsigmatotal=wetgridsigmatotal+ &
269                 wetgridsigma(ix,jy)
270        endif
271
272  ! DRY DEPOSITION
273        if ((DRYDEP).and.(ldirect.gt.0)) then
274            do l=1,nclassunc
275              auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
276            end do
277            call mean(auxgrid,drygrid(ix,jy), &
278                 drygridsigma(ix,jy),nclassunc)
279  ! Multiply by number of classes to get total concentration
280            drygrid(ix,jy)=drygrid(ix,jy)* &
281                 nclassunc
282            drygridtotal=drygridtotal+drygrid(ix,jy)
283  ! Calculate standard deviation of the mean
284            drygridsigma(ix,jy)= &
285                 drygridsigma(ix,jy)* &
286                 sqrt(real(nclassunc))
287125         drygridsigmatotal=drygridsigmatotal+ &
288                 drygridsigma(ix,jy)
289        endif
290
291  ! CONCENTRATION OR MIXING RATIO
292        do kz=1,numzgrid
293            do l=1,nclassunc
294              auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
295            end do
296            call mean(auxgrid,grid(ix,jy,kz), &
297                 gridsigma(ix,jy,kz),nclassunc)
298  ! Multiply by number of classes to get total concentration
299            grid(ix,jy,kz)= &
300                 grid(ix,jy,kz)*nclassunc
301            gridtotal=gridtotal+grid(ix,jy,kz)
302  ! Calculate standard deviation of the mean
303            gridsigma(ix,jy,kz)= &
304                 gridsigma(ix,jy,kz)* &
305                 sqrt(real(nclassunc))
306            gridsigmatotal=gridsigmatotal+ &
307                 gridsigma(ix,jy,kz)
308        end do
309      end do
310    end do
311
312
313
314
315  !*******************************************************************
316  ! Generate output: may be in concentration (ng/m3) or in mixing
317  ! ratio (ppt) or both
318  ! Output the position and the values alternated multiplied by
319  ! 1 or -1, first line is number of values, number of positions
320  ! For backward simulations, the unit is seconds, stored in grid_time
321  !*******************************************************************
322  if (verbosity.eq.1) then
323     print*,'concoutput_surf 4 (output)'
324     CALL SYSTEM_CLOCK(count_clock)
325     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
326  endif
327
328  ! Concentration output
329  !*********************
330  if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
331  if (verbosity.eq.1) then
332     print*,'concoutput_surf (Wet deposition)'
333     CALL SYSTEM_CLOCK(count_clock)
334     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
335  endif
336
337  ! Wet deposition
338         sp_count_i=0
339         sp_count_r=0
340         sp_fact=-1.
341         sp_zer=.true.
342         if ((ldirect.eq.1).and.(WETDEP)) then
343         do jy=0,numygrid-1
344            do ix=0,numxgrid-1
345  !oncentraion greater zero
346              if (wetgrid(ix,jy).gt.smallnum) then
347                 if (sp_zer.eqv..true.) then ! first non zero value
348                    sp_count_i=sp_count_i+1
349                    sparse_dump_i(sp_count_i)=ix+jy*numxgrid
350                    sp_zer=.false.
351                    sp_fact=sp_fact*(-1.)
352                 endif
353                 sp_count_r=sp_count_r+1
354                 sparse_dump_r(sp_count_r)= &
355                      sp_fact*1.e12*wetgrid(ix,jy)/area(ix,jy)
356                 sparse_dump_u(sp_count_r)= &
357                      1.e12*wetgridsigma(ix,jy)/area(ix,jy)
358              else ! concentration is zero
359                  sp_zer=.true.
360              endif
361            end do
362         end do
363         else
364            sp_count_i=0
365            sp_count_r=0
366         endif
367         write(unitoutgrid) sp_count_i
368         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
369         write(unitoutgrid) sp_count_r
370         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
371!         write(unitoutgrid) sp_count_r
372         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
373
374  if (verbosity.eq.1) then
375     print*,'concoutput_surf (Dry deposition)'
376     CALL SYSTEM_CLOCK(count_clock)
377     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
378  endif
379  ! Dry deposition
380         sp_count_i=0
381         sp_count_r=0
382         sp_fact=-1.
383         sp_zer=.true.
384         if ((ldirect.eq.1).and.(DRYDEP)) then
385          do jy=0,numygrid-1
386            do ix=0,numxgrid-1
387              if (drygrid(ix,jy).gt.smallnum) then
388                 if (sp_zer.eqv..true.) then ! first non zero value
389                    sp_count_i=sp_count_i+1
390                    sparse_dump_i(sp_count_i)=ix+jy*numxgrid
391                    sp_zer=.false.
392                    sp_fact=sp_fact*(-1.)
393                 endif
394                 sp_count_r=sp_count_r+1
395                 sparse_dump_r(sp_count_r)= &
396                      sp_fact* &
397                      1.e12*drygrid(ix,jy)/area(ix,jy)
398                  sparse_dump_u(sp_count_r)= &
399                      1.e12*drygridsigma(ix,jy)/area(ix,jy)
400              else ! concentration is zero
401                  sp_zer=.true.
402              endif
403            end do
404          end do
405         else
406            sp_count_i=0
407            sp_count_r=0
408         endif
409         write(unitoutgrid) sp_count_i
410         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
411         write(unitoutgrid) sp_count_r
412         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
413!         write(unitoutgrid) sp_count_r
414         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
415
416
417
418  if (verbosity.eq.1) then
419     print*,'concoutput_surf (Concentrations)'
420     CALL SYSTEM_CLOCK(count_clock)
421     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
422  endif
423  ! Concentrations
424
425  ! if surf_only write only 1st layer
426
427         if(surf_only.eq.1) then
428         sp_count_i=0
429         sp_count_r=0
430         sp_fact=-1.
431         sp_zer=.true.
432          do kz=1,1
433            do jy=0,numygrid-1
434              do ix=0,numxgrid-1
435                if (grid(ix,jy,kz).gt.smallnum) then
436                  if (sp_zer.eqv..true.) then ! first non zero value
437                    sp_count_i=sp_count_i+1
438                    sparse_dump_i(sp_count_i)= &
439                         ix+jy*numxgrid+kz*numxgrid*numygrid
440                    sp_zer=.false.
441                    sp_fact=sp_fact*(-1.)
442                   endif
443                   sp_count_r=sp_count_r+1
444                   sparse_dump_r(sp_count_r)= &
445                        sp_fact* &
446                        grid(ix,jy,kz)* &
447                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
448  !                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
449  !    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
450                   sparse_dump_u(sp_count_r)= &
451                        gridsigma(ix,jy,kz)* &
452                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
453              else ! concentration is zero
454                  sp_zer=.true.
455              endif
456              end do
457            end do
458          end do
459         write(unitoutgrid) sp_count_i
460         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
461         write(unitoutgrid) sp_count_r
462         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
463!         write(unitoutgrid) sp_count_r
464         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
465         else
466
467  ! write full vertical resolution
468  if (verbosity.eq.1) then
469     print*,'concoutput_surf (write full vertical resolution)'
470     CALL SYSTEM_CLOCK(count_clock)
471     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
472  endif
473
474         sp_count_i=0
475         sp_count_r=0
476         sp_fact=-1.
477         sp_zer=.true.
478          do kz=1,numzgrid
479            do jy=0,numygrid-1
480              do ix=0,numxgrid-1
481                if (grid(ix,jy,kz).gt.smallnum) then
482                  if (sp_zer.eqv..true.) then ! first non zero value
483                    sp_count_i=sp_count_i+1
484                    sparse_dump_i(sp_count_i)= &
485                         ix+jy*numxgrid+kz*numxgrid*numygrid
486                    sp_zer=.false.
487                    sp_fact=sp_fact*(-1.)
488                   endif
489                   sp_count_r=sp_count_r+1
490                   sparse_dump_r(sp_count_r)= &
491                        sp_fact* &
492                        grid(ix,jy,kz)* &
493                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
494  !                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
495  !    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
496                   sparse_dump_u(sp_count_r)= &
497                        gridsigma(ix,jy,kz)* &
498                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
499              else ! concentration is zero
500                  sp_zer=.true.
501              endif
502              end do
503            end do
504          end do
505         write(unitoutgrid) sp_count_i
506         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
507         write(unitoutgrid) sp_count_r
508         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
509!         write(unitoutgrid) sp_count_r
510         write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)         
511         endif ! surf_only
512
513    endif !  concentration output
514
515  ! Mixing ratio output
516  !********************
517
518  if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
519
520  ! Wet deposition
521         sp_count_i=0
522         sp_count_r=0
523         sp_fact=-1.
524         sp_zer=.true.
525         if ((ldirect.eq.1).and.(WETDEP)) then
526          do jy=0,numygrid-1
527            do ix=0,numxgrid-1
528                if (wetgrid(ix,jy).gt.smallnum) then
529                  if (sp_zer.eqv..true.) then ! first non zero value
530                    sp_count_i=sp_count_i+1
531                    sparse_dump_i(sp_count_i)= &
532                         ix+jy*numxgrid
533                    sp_zer=.false.
534                    sp_fact=sp_fact*(-1.)
535                 endif
536                 sp_count_r=sp_count_r+1
537                 sparse_dump_r(sp_count_r)= &
538                      sp_fact* &
539                      1.e12*wetgrid(ix,jy)/area(ix,jy)
540                 sparse_dump_u(sp_count_r)= &
541                      1.e12*wetgridsigma(ix,jy)/area(ix,jy)
542              else ! concentration is zero
543                  sp_zer=.true.
544              endif
545            end do
546          end do
547         else
548           sp_count_i=0
549           sp_count_r=0
550         endif
551         write(unitoutgridppt) sp_count_i
552         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
553         write(unitoutgridppt) sp_count_r
554         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
555!         write(unitoutgridppt) sp_count_r
556         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
557
558
559  ! Dry deposition
560         sp_count_i=0
561         sp_count_r=0
562         sp_fact=-1.
563         sp_zer=.true.
564         if ((ldirect.eq.1).and.(DRYDEP)) then
565          do jy=0,numygrid-1
566            do ix=0,numxgrid-1
567                if (drygrid(ix,jy).gt.smallnum) then
568                  if (sp_zer.eqv..true.) then ! first non zero value
569                    sp_count_i=sp_count_i+1
570                    sparse_dump_i(sp_count_i)= &
571                         ix+jy*numxgrid
572                    sp_zer=.false.
573                    sp_fact=sp_fact*(-1)
574                 endif
575                 sp_count_r=sp_count_r+1
576                 sparse_dump_r(sp_count_r)= &
577                      sp_fact* &
578                      1.e12*drygrid(ix,jy)/area(ix,jy)
579                 sparse_dump_u(sp_count_r)= &
580                      1.e12*drygridsigma(ix,jy)/area(ix,jy)
581              else ! concentration is zero
582                  sp_zer=.true.
583              endif
584            end do
585          end do
586         else
587           sp_count_i=0
588           sp_count_r=0
589         endif
590         write(unitoutgridppt) sp_count_i
591         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
592         write(unitoutgridppt) sp_count_r
593         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
594!         write(unitoutgridppt) sp_count_r
595         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
596
597
598  ! Mixing ratios
599
600  ! if surf_only write only 1st layer
601
602         if(surf_only.eq.1) then
603         sp_count_i=0
604         sp_count_r=0
605         sp_fact=-1.
606         sp_zer=.true.
607          do kz=1,1
608            do jy=0,numygrid-1
609              do ix=0,numxgrid-1
610                if (grid(ix,jy,kz).gt.smallnum) then
611                  if (sp_zer.eqv..true.) then ! first non zero value
612                    sp_count_i=sp_count_i+1
613                    sparse_dump_i(sp_count_i)= &
614                         ix+jy*numxgrid+kz*numxgrid*numygrid
615                    sp_zer=.false.
616                    sp_fact=sp_fact*(-1.)
617                 endif
618                 sp_count_r=sp_count_r+1
619                 sparse_dump_r(sp_count_r)= &
620                      sp_fact* &
621                      1.e12*grid(ix,jy,kz) &
622                      /volume(ix,jy,kz)/outnum* &
623                      weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
624                 sparse_dump_u(sp_count_r)= &
625                      1.e12*gridsigma(ix,jy,kz)/volume(ix,jy,kz)/ &
626                      outnum*weightair/weightmolar(ks)/ &
627                      densityoutgrid(ix,jy,kz)
628              else ! concentration is zero
629                  sp_zer=.true.
630              endif
631              end do
632            end do
633          end do
634         write(unitoutgridppt) sp_count_i
635         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
636         write(unitoutgridppt) sp_count_r
637         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
638!         write(unitoutgridppt) sp_count_r
639         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
640         else
641
642  ! write full vertical resolution
643
644         sp_count_i=0
645         sp_count_r=0
646         sp_fact=-1.
647         sp_zer=.true.
648          do kz=1,numzgrid
649            do jy=0,numygrid-1
650              do ix=0,numxgrid-1
651                if (grid(ix,jy,kz).gt.smallnum) then
652                  if (sp_zer.eqv..true.) then ! first non zero value
653                    sp_count_i=sp_count_i+1
654                    sparse_dump_i(sp_count_i)= &
655                         ix+jy*numxgrid+kz*numxgrid*numygrid
656                    sp_zer=.false.
657                    sp_fact=sp_fact*(-1.)
658                 endif
659                 sp_count_r=sp_count_r+1
660                 sparse_dump_r(sp_count_r)= &
661                      sp_fact* &
662                      1.e12*grid(ix,jy,kz) &
663                      /volume(ix,jy,kz)/outnum* &
664                      weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
665                 sparse_dump_u(sp_count_r)= &
666                      1.e12*gridsigma(ix,jy,kz)/volume(ix,jy,kz)/ &
667                      outnum*weightair/weightmolar(ks)/ &
668                      densityoutgrid(ix,jy,kz)
669              else ! concentration is zero
670                  sp_zer=.true.
671              endif
672              end do
673            end do
674          end do
675         write(unitoutgridppt) sp_count_i
676         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
677         write(unitoutgridppt) sp_count_r
678         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
679!         write(unitoutgridppt) sp_count_r
680         write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
681         endif ! surf_only
682
683      endif ! output for ppt
684
685  end do
686  end do
687
688    close(unitoutgridppt)
689    close(unitoutgrid)
690
691  end do
692
693  if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
694  if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
695       wetgridtotal
696  if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
697       drygridtotal
698
699  ! Dump of receptor concentrations
700
701    if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3)  ) then
702      write(unitoutreceptppt) itime
703      do ks=1,nspec
704        write(unitoutreceptppt) (1.e12*creceptor(i,ks)/outnum* &
705             weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)
706      end do
707    endif
708
709  ! Dump of receptor concentrations
710
711    if (numreceptor.gt.0) then
712      write(unitoutrecept) itime
713      do ks=1,nspec
714        write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum, &
715             i=1,numreceptor)
716      end do
717    endif
718
719
720
721  ! Reinitialization of grid
722  !*************************
723
724  do ks=1,nspec
725  do kp=1,maxpointspec_act
726    do i=1,numreceptor
727      creceptor(i,ks)=0.
728    end do
729    do jy=0,numygrid-1
730      do ix=0,numxgrid-1
731        do l=1,nclassunc
732          do nage=1,nageclass
733            do kz=1,numzgrid
734              gridunc(ix,jy,kz,ks,kp,l,nage)=0.
735            end do
736          end do
737        end do
738      end do
739    end do
740  end do
741  end do
742
743
744end subroutine concoutput_surf
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG