source: flexpart.git/flexpart_code/concoutput_nest.F90 @ b398fb6

FPv9.3.1FPv9.3.1b_testingFPv9.3.2fp9.3.1-20161214-nc4grib2nc4_repair
Last change on this file since b398fb6 was 496c607, checked in by Don Morton <Don.Morton@…>, 8 years ago

Initial commit of FPv9.3.1

Currently, this is a clone of snapshot FPv9.3.0

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