source: trunk/src/concoutput_nest.f90 @ 4

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