source: flexpart.git/src/concoutput_nest.f90 @ f9ce123

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since f9ce123 was 5f9d14a, checked in by Espen Sollum ATMOS <eso@…>, 9 years ago

Updated wet depo scheme

  • Property mode set to 100644
File size: 20.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_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  integer :: mind
98! mind        eso:added to ensure identical results between 2&3-fields versions
99
100
101  ! Determine current calendar date, needed for the file name
102  !**********************************************************
103
104  jul=bdate+real(itime,kind=dp)/86400._dp
105  call caldate(jul,jjjjmmdd,ihmmss)
106  write(adate,'(i8.8)') jjjjmmdd
107  write(atime,'(i6.6)') ihmmss
108
109
110  ! For forward simulations, output fields have dimension MAXSPEC,
111  ! for backward simulations, output fields have dimension MAXPOINT.
112  ! Thus, make loops either about nspec, or about numpoint
113  !*****************************************************************
114
115
116    if (ldirect.eq.1) then
117       do ks=1,nspec
118         do kp=1,maxpointspec_act
119           tot_mu(ks,kp)=1
120         end do
121       end do
122   else
123      do ks=1,nspec
124             do kp=1,maxpointspec_act
125               tot_mu(ks,kp)=xmass(kp,ks)
126             end do
127      end do
128    endif
129
130
131  !*******************************************************************
132  ! Compute air density: sufficiently accurate to take it
133  ! from coarse grid at some time
134  ! Determine center altitude of output layer, and interpolate density
135  ! data to that altitude
136  !*******************************************************************
137
138  mind=memind(2)
139  do kz=1,numzgrid
140    if (kz.eq.1) then
141      halfheight=outheight(1)/2.
142    else
143      halfheight=(outheight(kz)+outheight(kz-1))/2.
144    endif
145    do kzz=2,nz
146      if ((height(kzz-1).lt.halfheight).and. &
147           (height(kzz).gt.halfheight)) goto 46
148    end do
14946   kzz=max(min(kzz,nz),2)
150    dz1=halfheight-height(kzz-1)
151    dz2=height(kzz)-halfheight
152    dz=dz1+dz2
153    do jy=0,numygridn-1
154      do ix=0,numxgridn-1
155        xl=outlon0n+real(ix)*dxoutn
156        yl=outlat0n+real(jy)*dyoutn
157        xl=(xl-xlon0)/dx
158        yl=(yl-ylat0)/dy
159        iix=max(min(nint(xl),nxmin1),0)
160        jjy=max(min(nint(yl),nymin1),0)
161        ! densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
162        !      rho(iix,jjy,kzz-1,2)*dz2)/dz
163        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,mind)*dz1+ &
164             rho(iix,jjy,kzz-1,mind)*dz2)/dz
165      end do
166    end do
167  end do
168
169  do i=1,numreceptor
170    xl=xreceptor(i)
171    yl=yreceptor(i)
172    iix=max(min(nint(xl),nxmin1),0)
173    jjy=max(min(nint(yl),nymin1),0)
174    !densityoutrecept(i)=rho(iix,jjy,1,2)
175    densityoutrecept(i)=rho(iix,jjy,1,mind)
176  end do
177
178
179  ! Output is different for forward and backward simulations
180    do kz=1,numzgrid
181      do jy=0,numygridn-1
182        do ix=0,numxgridn-1
183          if (ldirect.eq.1) then
184            factor3d(ix,jy,kz)=1.e12/volumen(ix,jy,kz)/outnum
185          else
186            factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
187          endif
188        end do
189      end do
190    end do
191
192  !*********************************************************************
193  ! Determine the standard deviation of the mean concentration or mixing
194  ! ratio (uncertainty of the output) and the dry and wet deposition
195  !*********************************************************************
196
197  do ks=1,nspec
198
199  write(anspec,'(i3.3)') ks
200  if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
201    if (ldirect.eq.1) then
202      open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_' &
203           //adate// &
204           atime//'_'//anspec,form='unformatted')
205    else
206      open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_nest_' &
207           //adate// &
208           atime//'_'//anspec,form='unformatted')
209    endif
210     write(unitoutgrid) itime
211   endif
212
213  if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
214   open(unitoutgridppt,file=path(2)(1:length(2))//'grid_pptv_nest_' &
215        //adate// &
216        atime//'_'//anspec,form='unformatted')
217
218    write(unitoutgridppt) itime
219  endif
220
221  do kp=1,maxpointspec_act
222  do nage=1,nageclass
223
224    do jy=0,numygridn-1
225      do ix=0,numxgridn-1
226
227  ! WET DEPOSITION
228        if ((WETDEP).and.(ldirect.gt.0)) then
229            do l=1,nclassunc
230              auxgrid(l)=wetgriduncn(ix,jy,ks,kp,l,nage)
231            end do
232            call mean(auxgrid,wetgrid(ix,jy), &
233                 wetgridsigma(ix,jy),nclassunc)
234  ! Multiply by number of classes to get total concentration
235            wetgrid(ix,jy)=wetgrid(ix,jy) &
236                 *nclassunc
237  ! Calculate standard deviation of the mean
238            wetgridsigma(ix,jy)= &
239                 wetgridsigma(ix,jy)* &
240                 sqrt(real(nclassunc))
241        endif
242
243  ! DRY DEPOSITION
244        if ((DRYDEP).and.(ldirect.gt.0)) then
245            do l=1,nclassunc
246              auxgrid(l)=drygriduncn(ix,jy,ks,kp,l,nage)
247            end do
248            call mean(auxgrid,drygrid(ix,jy), &
249                 drygridsigma(ix,jy),nclassunc)
250  ! Multiply by number of classes to get total concentration
251            drygrid(ix,jy)=drygrid(ix,jy)* &
252                 nclassunc
253  ! Calculate standard deviation of the mean
254            drygridsigma(ix,jy)= &
255                 drygridsigma(ix,jy)* &
256                 sqrt(real(nclassunc))
257        endif
258
259  ! CONCENTRATION OR MIXING RATIO
260        do kz=1,numzgrid
261            do l=1,nclassunc
262              auxgrid(l)=griduncn(ix,jy,kz,ks,kp,l,nage)
263            end do
264            call mean(auxgrid,grid(ix,jy,kz), &
265                 gridsigma(ix,jy,kz),nclassunc)
266  ! Multiply by number of classes to get total concentration
267            grid(ix,jy,kz)= &
268                 grid(ix,jy,kz)*nclassunc
269  ! Calculate standard deviation of the mean
270            gridsigma(ix,jy,kz)= &
271                 gridsigma(ix,jy,kz)* &
272                 sqrt(real(nclassunc))
273        end do
274      end do
275    end do
276
277
278  !*******************************************************************
279  ! Generate output: may be in concentration (ng/m3) or in mixing
280  ! ratio (ppt) or both
281  ! Output the position and the values alternated multiplied by
282  ! 1 or -1, first line is number of values, number of positions
283  ! For backward simulations, the unit is seconds, stored in grid_time
284  !*******************************************************************
285
286  ! Concentration output
287  !*********************
288  if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
289
290  ! Wet deposition
291         sp_count_i=0
292         sp_count_r=0
293         sp_fact=-1.
294         sp_zer=.true.
295         if ((ldirect.eq.1).and.(WETDEP)) then
296         do jy=0,numygridn-1
297            do ix=0,numxgridn-1
298  !oncentraion greater zero
299              if (wetgrid(ix,jy).gt.smallnum) then
300                 if (sp_zer.eqv..true.) then ! first non zero value
301                    sp_count_i=sp_count_i+1
302                    sparse_dump_i(sp_count_i)=ix+jy*numxgridn
303                    sp_zer=.false.
304                    sp_fact=sp_fact*(-1.)
305                 endif
306                 sp_count_r=sp_count_r+1
307                 sparse_dump_r(sp_count_r)= &
308                      sp_fact*1.e12*wetgrid(ix,jy)/arean(ix,jy)
309  !                sparse_dump_u(sp_count_r)=
310  !+                1.e12*wetgridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
311              else ! concentration is zero
312                  sp_zer=.true.
313              endif
314            end do
315         end do
316         else
317            sp_count_i=0
318            sp_count_r=0
319         endif
320         write(unitoutgrid) sp_count_i
321         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
322         write(unitoutgrid) sp_count_r
323         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
324  !       write(unitoutgrid) sp_count_u
325  !       write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
326
327  ! Dry deposition
328         sp_count_i=0
329         sp_count_r=0
330         sp_fact=-1.
331         sp_zer=.true.
332         if ((ldirect.eq.1).and.(DRYDEP)) then
333          do jy=0,numygridn-1
334            do ix=0,numxgridn-1
335              if (drygrid(ix,jy).gt.smallnum) then
336                 if (sp_zer.eqv..true.) then ! first non zero value
337                    sp_count_i=sp_count_i+1
338                    sparse_dump_i(sp_count_i)=ix+jy*numxgridn
339                    sp_zer=.false.
340                    sp_fact=sp_fact*(-1.)
341                 endif
342                 sp_count_r=sp_count_r+1
343                 sparse_dump_r(sp_count_r)= &
344                      sp_fact* &
345                      1.e12*drygrid(ix,jy)/arean(ix,jy)
346  !                sparse_dump_u(sp_count_r)=
347  !+                1.e12*drygridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
348              else ! concentration is zero
349                  sp_zer=.true.
350              endif
351            end do
352          end do
353         else
354            sp_count_i=0
355            sp_count_r=0
356         endif
357         write(unitoutgrid) sp_count_i
358         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
359         write(unitoutgrid) sp_count_r
360         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
361  !       write(*,*) sp_count_u
362  !       write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
363
364
365
366  ! Concentrations
367         sp_count_i=0
368         sp_count_r=0
369         sp_fact=-1.
370         sp_zer=.true.
371          do kz=1,numzgrid
372            do jy=0,numygridn-1
373              do ix=0,numxgridn-1
374                if (grid(ix,jy,kz).gt.smallnum) then
375                  if (sp_zer.eqv..true.) then ! first non zero value
376                    sp_count_i=sp_count_i+1
377                    sparse_dump_i(sp_count_i)= &
378                         ix+jy*numxgridn+kz*numxgridn*numygridn
379                    sp_zer=.false.
380                    sp_fact=sp_fact*(-1.)
381                   endif
382                   sp_count_r=sp_count_r+1
383                   sparse_dump_r(sp_count_r)= &
384                        sp_fact* &
385                        grid(ix,jy,kz)* &
386                        factor3d(ix,jy,kz)/tot_mu(ks,kp)
387  !                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
388  !    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
389  !                sparse_dump_u(sp_count_r)=
390  !+               ,gridsigma(ix,jy,kz,ks,kp,nage)*
391  !+               factor(ix,jy,kz)/tot_mu(ks,kp)
392              else ! concentration is zero
393                  sp_zer=.true.
394              endif
395              end do
396            end do
397          end do
398         write(unitoutgrid) sp_count_i
399         write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
400         write(unitoutgrid) sp_count_r
401         write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
402  !       write(unitoutgrid) sp_count_u
403  !       write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
404
405
406
407    endif !  concentration output
408
409  ! Mixing ratio output
410  !********************
411
412  if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
413
414  ! Wet deposition
415         sp_count_i=0
416         sp_count_r=0
417         sp_fact=-1.
418         sp_zer=.true.
419         if ((ldirect.eq.1).and.(WETDEP)) then
420          do jy=0,numygridn-1
421            do ix=0,numxgridn-1
422                if (wetgrid(ix,jy).gt.smallnum) then
423                  if (sp_zer.eqv..true.) then ! first non zero value
424                    sp_count_i=sp_count_i+1
425                    sparse_dump_i(sp_count_i)= &
426                         ix+jy*numxgridn
427                    sp_zer=.false.
428                    sp_fact=sp_fact*(-1.)
429                 endif
430                 sp_count_r=sp_count_r+1
431                 sparse_dump_r(sp_count_r)= &
432                      sp_fact* &
433                      1.e12*wetgrid(ix,jy)/arean(ix,jy)
434  !                sparse_dump_u(sp_count_r)=
435  !    +            ,1.e12*wetgridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
436              else ! concentration is zero
437                  sp_zer=.true.
438              endif
439            end do
440          end do
441         else
442           sp_count_i=0
443           sp_count_r=0
444         endif
445         write(unitoutgridppt) sp_count_i
446         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
447         write(unitoutgridppt) sp_count_r
448         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
449  !       write(unitoutgridppt) sp_count_u
450  !       write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
451
452
453  ! Dry deposition
454         sp_count_i=0
455         sp_count_r=0
456         sp_fact=-1.
457         sp_zer=.true.
458         if ((ldirect.eq.1).and.(DRYDEP)) then
459          do jy=0,numygridn-1
460            do ix=0,numxgridn-1
461                if (drygrid(ix,jy).gt.smallnum) then
462                  if (sp_zer.eqv..true.) then ! first non zero value
463                    sp_count_i=sp_count_i+1
464                    sparse_dump_i(sp_count_i)= &
465                         ix+jy*numxgridn
466                    sp_zer=.false.
467                    sp_fact=sp_fact*(-1)
468                 endif
469                 sp_count_r=sp_count_r+1
470                 sparse_dump_r(sp_count_r)= &
471                      sp_fact* &
472                      1.e12*drygrid(ix,jy)/arean(ix,jy)
473  !                sparse_dump_u(sp_count_r)=
474  !    +            ,1.e12*drygridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
475              else ! concentration is zero
476                  sp_zer=.true.
477              endif
478            end do
479          end do
480         else
481           sp_count_i=0
482           sp_count_r=0
483         endif
484         write(unitoutgridppt) sp_count_i
485         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
486         write(unitoutgridppt) sp_count_r
487         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
488  !       write(unitoutgridppt) sp_count_u
489  !       write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
490
491
492  ! Mixing ratios
493         sp_count_i=0
494         sp_count_r=0
495         sp_fact=-1.
496         sp_zer=.true.
497          do kz=1,numzgrid
498            do jy=0,numygridn-1
499              do ix=0,numxgridn-1
500                if (grid(ix,jy,kz).gt.smallnum) then
501                  if (sp_zer.eqv..true.) then ! first non zero value
502                    sp_count_i=sp_count_i+1
503                    sparse_dump_i(sp_count_i)= &
504                         ix+jy*numxgridn+kz*numxgridn*numygridn
505                    sp_zer=.false.
506                    sp_fact=sp_fact*(-1.)
507                 endif
508                 sp_count_r=sp_count_r+1
509                 sparse_dump_r(sp_count_r)= &
510                      sp_fact* &
511                      1.e12*grid(ix,jy,kz) &
512                      /volumen(ix,jy,kz)/outnum* &
513                      weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
514  !                sparse_dump_u(sp_count_r)=
515  !+              ,1.e12*gridsigma(ix,jy,kz,ks,kp,nage)/volume(ix,jy,kz)/
516  !+              outnum*weightair/weightmolar(ks)/
517  !+              densityoutgrid(ix,jy,kz)
518              else ! concentration is zero
519                  sp_zer=.true.
520              endif
521              end do
522            end do
523          end do
524         write(unitoutgridppt) sp_count_i
525         write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
526         write(unitoutgridppt) sp_count_r
527         write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
528  !       write(unitoutgridppt) sp_count_u
529  !       write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
530
531      endif ! output for ppt
532
533  end do
534  end do
535
536    close(unitoutgridppt)
537    close(unitoutgrid)
538
539  end do
540
541
542
543  ! Reinitialization of grid
544  !*************************
545
546  do ks=1,nspec
547  do kp=1,maxpointspec_act
548    do i=1,numreceptor
549      creceptor(i,ks)=0.
550    end do
551    do jy=0,numygridn-1
552      do ix=0,numxgridn-1
553        do l=1,nclassunc
554          do nage=1,nageclass
555            do kz=1,numzgrid
556              griduncn(ix,jy,kz,ks,kp,l,nage)=0.
557            end do
558          end do
559        end do
560      end do
561    end do
562  end do
563  end do
564
565
566end subroutine concoutput_nest
567
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG