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