source: flexpart.git/src/concoutput_nest.f90 @ 7d02c2f

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

Added BW scavenging for nested domains

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