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

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

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

File size: 19.4 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_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  !     PS, 2/2015: option to produce incremental deposition output
48  !                 access= -> position=
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)
93  real :: halfheight,dz,dz1,dz2,tot_mu(maxspec,maxpointspec_act)
94  real,parameter :: smallnum = tiny(0.0) ! smallest number that can be handled
95  real,parameter :: weightair=28.97
96  logical :: sp_zer
97  character :: adate*8,atime*6
98  character(len=3) :: anspec
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  do kz=1,numzgrid
139    if (kz.eq.1) then
140      halfheight=outheight(1)/2.
141    else
142      halfheight=(outheight(kz)+outheight(kz-1))/2.
143    endif
144    do kzz=2,nz
145      if ((height(kzz-1).lt.halfheight).and. &
146           (height(kzz).gt.halfheight)) goto 46
147    end do
14846   kzz=max(min(kzz,nz),2)
149    dz1=halfheight-height(kzz-1)
150    dz2=height(kzz)-halfheight
151    dz=dz1+dz2
152    do jy=0,numygridn-1
153      do ix=0,numxgridn-1
154        xl=outlon0n+real(ix)*dxoutn
155        yl=outlat0n+real(jy)*dyoutn
156        xl=(xl-xlon0)/dx
157        yl=(yl-ylat0)/dy
158        iix=max(min(nint(xl),nxmin1),0)
159        jjy=max(min(nint(yl),nymin1),0)
160        densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,2)*dz1+ &
161             rho(iix,jjy,kzz-1,2)*dz2)/dz
162      end do
163    end do
164  end do
165
166  do i=1,numreceptor
167    xl=xreceptor(i)
168    yl=yreceptor(i)
169    iix=max(min(nint(xl),nxmin1),0)
170    jjy=max(min(nint(yl),nymin1),0)
171    densityoutrecept(i)=rho(iix,jjy,1,2)
172  end do
173
174
175! Output is different for forward and backward simulations
176  do kz=1,numzgrid
177    do jy=0,numygridn-1
178      do ix=0,numxgridn-1
179        if (ldirect.eq.1) then
180          factor3d(ix,jy,kz)=1.e12/volumen(ix,jy,kz)/outnum
181        else
182          factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
183        endif
184      end do
185    end do
186  end do
187
188  !*********************************************************************
189  ! Determine the standard deviation of the mean concentration or mixing
190  ! ratio (uncertainty of the output) and the dry and wet deposition
191  !*********************************************************************
192
193  do ks=1,nspec
194
195    write(anspec,'(i3.3)') ks
196    if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
197      if (ldirect.eq.1) then
198        open(unitoutgrid,file=path(2)(1:length(2))//'grid_conc_nest_' &
199             //adate// atime//'_'//anspec,form='unformatted')
200      else
201        open(unitoutgrid,file=path(2)(1:length(2))//'grid_time_nest_' &
202             //adate// 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// atime//'_'//anspec,form='unformatted')
210
211      write(unitoutgridppt) itime
212    endif
213
214    do kp=1,maxpointspec_act
215    do nage=1,nageclass
216
217      do jy=0,numygridn-1
218        do ix=0,numxgridn-1
219
220    ! WET DEPOSITION
221          if ((WETDEP).and.(ldirect.gt.0)) then
222            do l=1,nclassunc
223              auxgrid(l)=wetgriduncn(ix,jy,ks,kp,l,nage)
224            end do
225            call mean(auxgrid,wetgrid(ix,jy), wetgridsigma(ix,jy),nclassunc)
226    ! Multiply by number of classes to get total concentration
227            wetgrid(ix,jy)=wetgrid(ix,jy)*nclassunc
228    ! Calculate standard deviation of the mean
229            wetgridsigma(ix,jy)= wetgridsigma(ix,jy)* sqrt(real(nclassunc))
230          endif
231
232    ! DRY DEPOSITION
233          if ((DRYDEP).and.(ldirect.gt.0)) then
234            do l=1,nclassunc
235              auxgrid(l)=drygriduncn(ix,jy,ks,kp,l,nage)
236            end do
237            call mean(auxgrid,drygrid(ix,jy), drygridsigma(ix,jy),nclassunc)
238    ! Multiply by number of classes to get total concentration
239            drygrid(ix,jy)=drygrid(ix,jy)* nclassunc
240    ! Calculate standard deviation of the mean
241            drygridsigma(ix,jy)= drygridsigma(ix,jy)* sqrt(real(nclassunc))
242          endif
243
244    ! CONCENTRATION OR MIXING RATIO
245          do kz=1,numzgrid
246            do l=1,nclassunc
247              auxgrid(l)=griduncn(ix,jy,kz,ks,kp,l,nage)
248            end do
249            call mean(auxgrid,grid(ix,jy,kz), gridsigma(ix,jy,kz),nclassunc)
250    ! Multiply by number of classes to get total concentration
251            grid(ix,jy,kz)= grid(ix,jy,kz)*nclassunc
252    ! Calculate standard deviation of the mean
253            gridsigma(ix,jy,kz)= gridsigma(ix,jy,kz)* sqrt(real(nclassunc))
254          end do
255        end do
256      end do
257
258
259    !*******************************************************************
260    ! Generate output: may be in concentration (ng/m3) or in mixing
261    ! ratio (ppt) or both
262    ! Output the position and the values alternated multiplied by
263    ! 1 or -1, first line is number of values, number of positions
264    ! For backward simulations, the unit is seconds, stored in grid_time
265    !*******************************************************************
266
267    ! Concentration output
268    !*********************
269    if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
270
271    ! Wet deposition
272      sp_count_i=0
273      sp_count_r=0
274      sp_fact=-1.
275      sp_zer=.true.
276      if ((ldirect.eq.1).and.(WETDEP)) then
277      do jy=0,numygridn-1
278        do ix=0,numxgridn-1
279    ! concentration greater zero
280          if (wetgrid(ix,jy).gt.smallnum) then
281            if (sp_zer.eqv..true.) then ! first non zero value
282              sp_count_i=sp_count_i+1
283              sparse_dump_i(sp_count_i)=ix+jy*numxgridn
284              sp_zer=.false.
285              sp_fact=sp_fact*(-1.)
286            endif
287            sp_count_r=sp_count_r+1
288            sparse_dump_r(sp_count_r)= sp_fact*1.e12*wetgrid(ix,jy)/arean(ix,jy)
289  !                sparse_dump_u(sp_count_r)=
290  !+                1.e12*wetgridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
291          else ! concentration is zero
292            sp_zer=.true.
293          endif
294       end do
295      end do
296      else
297        sp_count_i=0
298        sp_count_r=0
299      endif
300      write(unitoutgrid) sp_count_i
301      write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
302      write(unitoutgrid) sp_count_r
303      write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
304    !       write(unitoutgrid) sp_count_u
305    !       write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
306
307    ! Dry deposition
308      sp_count_i=0
309      sp_count_r=0
310      sp_fact=-1.
311      sp_zer=.true.
312      if ((ldirect.eq.1).and.(DRYDEP)) then
313       do jy=0,numygridn-1
314         do ix=0,numxgridn-1
315           if (drygrid(ix,jy).gt.smallnum) then
316             if (sp_zer.eqv..true.) then ! first non zero value
317                sp_count_i=sp_count_i+1
318                sparse_dump_i(sp_count_i)=ix+jy*numxgridn
319                sp_zer=.false.
320                sp_fact=sp_fact*(-1.)
321             endif
322             sp_count_r=sp_count_r+1
323             sparse_dump_r(sp_count_r)= sp_fact* &
324                  1.e12*drygrid(ix,jy)/arean(ix,jy)
325  !                sparse_dump_u(sp_count_r)=
326  !+                1.e12*drygridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
327           else ! concentration is zero
328             sp_zer=.true.
329           endif
330         end do
331       end do
332      else
333        sp_count_i=0
334        sp_count_r=0
335      endif
336      write(unitoutgrid) sp_count_i
337      write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
338      write(unitoutgrid) sp_count_r
339      write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
340    !       write(*,*) sp_count_u
341    !       write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
342
343
344
345    ! Concentrations
346      sp_count_i=0
347      sp_count_r=0
348      sp_fact=-1.
349      sp_zer=.true.
350      do kz=1,numzgrid
351        do jy=0,numygridn-1
352          do ix=0,numxgridn-1
353            if (grid(ix,jy,kz).gt.smallnum) then
354              if (sp_zer.eqv..true.) then ! first non zero value
355                sp_count_i=sp_count_i+1
356                sparse_dump_i(sp_count_i)= &
357                      ix+jy*numxgridn+kz*numxgridn*numygridn
358                sp_zer=.false.
359                sp_fact=sp_fact*(-1.)
360              endif
361              sp_count_r=sp_count_r+1
362              sparse_dump_r(sp_count_r)= sp_fact* grid(ix,jy,kz)* &
363                     factor3d(ix,jy,kz)/tot_mu(ks,kp)
364  !                 if ((factor(ix,jy,kz)/tot_mu(ks,kp)).eq.0)
365  !    +              write (*,*) factor(ix,jy,kz),tot_mu(ks,kp),ks,kp
366  !                sparse_dump_u(sp_count_r)=
367  !+               ,gridsigma(ix,jy,kz,ks,kp,nage)*
368  !+               factor(ix,jy,kz)/tot_mu(ks,kp)
369            else ! concentration is zero
370              sp_zer=.true.
371            endif
372          end do
373        end do
374      end do
375      write(unitoutgrid) sp_count_i
376      write(unitoutgrid) (sparse_dump_i(i),i=1,sp_count_i)
377      write(unitoutgrid) sp_count_r
378      write(unitoutgrid) (sparse_dump_r(i),i=1,sp_count_r)
379    !       write(unitoutgrid) sp_count_u
380    !       write(unitoutgrid) (sparse_dump_u(i),i=1,sp_count_r)
381
382
383
384      endif !  concentration output
385
386    ! Mixing ratio output
387    !********************
388
389      if ((iout.eq.2).or.(iout.eq.3)) then      ! mixing ratio
390
391      ! Wet deposition
392        sp_count_i=0
393        sp_count_r=0
394        sp_fact=-1.
395        sp_zer=.true.
396        if ((ldirect.eq.1).and.(WETDEP)) then
397          do jy=0,numygridn-1
398            do ix=0,numxgridn-1
399              if (wetgrid(ix,jy).gt.smallnum) then
400                if (sp_zer.eqv..true.) then ! first non zero value
401                  sp_count_i=sp_count_i+1
402                  sparse_dump_i(sp_count_i)= ix+jy*numxgridn
403                  sp_zer=.false.
404                  sp_fact=sp_fact*(-1.)
405                endif
406                sp_count_r=sp_count_r+1
407                sparse_dump_r(sp_count_r)= sp_fact* &
408                   1.e12*wetgrid(ix,jy)/arean(ix,jy)
409    !                sparse_dump_u(sp_count_r)=
410    !    +            ,1.e12*wetgridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
411              else ! concentration is zero
412                sp_zer=.true.
413              endif
414            end do
415          end do
416        else
417          sp_count_i=0
418          sp_count_r=0
419        endif
420        write(unitoutgridppt) sp_count_i
421        write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
422        write(unitoutgridppt) sp_count_r
423        write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
424    !       write(unitoutgridppt) sp_count_u
425    !       write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
426
427
428      ! Dry deposition
429        sp_count_i=0
430        sp_count_r=0
431        sp_fact=-1.
432        sp_zer=.true.
433        if ((ldirect.eq.1).and.(DRYDEP)) then
434          do jy=0,numygridn-1
435            do ix=0,numxgridn-1
436              if (drygrid(ix,jy).gt.smallnum) then
437                if (sp_zer.eqv..true.) then ! first non zero value
438                  sp_count_i=sp_count_i+1
439                  sparse_dump_i(sp_count_i)= ix+jy*numxgridn
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)= sp_fact* &
445                     1.e12*drygrid(ix,jy)/arean(ix,jy)
446    !                sparse_dump_u(sp_count_r)=
447    !    +            ,1.e12*drygridsigma(ix,jy,ks,kp,nage)/area(ix,jy)
448              else ! concentration is zero
449                sp_zer=.true.
450              endif
451            end do
452          end do
453        else
454          sp_count_i=0
455          sp_count_r=0
456        endif
457        write(unitoutgridppt) sp_count_i
458        write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
459        write(unitoutgridppt) sp_count_r
460        write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
461    !       write(unitoutgridppt) sp_count_u
462    !       write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
463
464
465      ! Mixing ratios
466        sp_count_i=0
467        sp_count_r=0
468        sp_fact=-1.
469        sp_zer=.true.
470        do kz=1,numzgrid
471           do jy=0,numygridn-1
472              do ix=0,numxgridn-1
473                if (grid(ix,jy,kz).gt.smallnum) then
474                  if (sp_zer.eqv..true.) then ! first non zero value
475                    sp_count_i=sp_count_i+1
476                    sparse_dump_i(sp_count_i)= &
477                        ix+jy*numxgridn+kz*numxgridn*numygridn
478                    sp_zer=.false.
479                    sp_fact=sp_fact*(-1.)
480                 endif
481                 sp_count_r=sp_count_r+1
482                 sparse_dump_r(sp_count_r)= sp_fact* 1.e12*grid(ix,jy,kz) &
483                     /volumen(ix,jy,kz)/outnum* &
484                     weightair/weightmolar(ks)/densityoutgrid(ix,jy,kz)
485    !                sparse_dump_u(sp_count_r)=
486    !+              ,1.e12*gridsigma(ix,jy,kz,ks,kp,nage)/volume(ix,jy,kz)/
487    !+              outnum*weightair/weightmolar(ks)/
488    !+              densityoutgrid(ix,jy,kz)
489                else ! concentration is zero
490                sp_zer=.true.
491              endif
492            end do
493          end do
494        end do
495        write(unitoutgridppt) sp_count_i
496        write(unitoutgridppt) (sparse_dump_i(i),i=1,sp_count_i)
497        write(unitoutgridppt) sp_count_r
498        write(unitoutgridppt) (sparse_dump_r(i),i=1,sp_count_r)
499    !       write(unitoutgridppt) sp_count_u
500    !       write(unitoutgridppt) (sparse_dump_u(i),i=1,sp_count_r)
501
502      endif ! output for ppt
503
504    end do
505    end do
506
507    close(unitoutgridppt)
508    close(unitoutgrid)
509
510  end do
511
512
513
514  ! Reinitialization of grid
515  !*************************
516
517  griduncn=0.
518  if (ldep_incr) then ! incremental deposition output
519    wetgriduncn=0.
520    drygriduncn=0.
521  endif
522
523end subroutine concoutput_nest
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG