source: trunk/src/concoutput_surf.f90 @ 28

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

move version 9.1.8 form branches to trunk. Contributions from HSO, saeck, pesei, NIK, RT, XKF, IP and others

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', ACCESS='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