source: trunk/src/fluxoutput.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 11.0 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 fluxoutput(itime)
23  !                        i
24  !*****************************************************************************
25  !                                                                            *
26  !     Output of the gridded fluxes.                                          *
27  !     Eastward, westward, northward, southward, upward and downward gross    *
28  !     fluxes are written to output file in either sparse matrix or grid dump *
29  !     format, whichever is more efficient.                                   *
30  !                                                                            *
31  !     Author: A. Stohl                                                       *
32  !                                                                            *
33  !     04 April 2000                                                          *
34  !                                                                            *
35  !*****************************************************************************
36  !                                                                            *
37  ! Variables:                                                                 *
38  ! ncellse         number of cells with non-zero values for eastward fluxes   *
39  ! sparsee         .true. if in sparse matrix format, else .false.            *
40  !                                                                            *
41  !*****************************************************************************
42
43  use flux_mod
44  use outg_mod
45  use par_mod
46  use com_mod
47
48  implicit none
49
50  real(kind=dp) :: jul
51  integer :: itime,ix,jy,kz,k,nage,jjjjmmdd,ihmmss,kp,i
52  integer :: ncellse(maxspec,maxageclass),ncellsw(maxspec,maxageclass)
53  integer :: ncellss(maxspec,maxageclass),ncellsn(maxspec,maxageclass)
54  integer :: ncellsu(maxspec,maxageclass),ncellsd(maxspec,maxageclass)
55  logical :: sparsee(maxspec,maxageclass),sparsew(maxspec,maxageclass)
56  logical :: sparses(maxspec,maxageclass),sparsen(maxspec,maxageclass)
57  logical :: sparseu(maxspec,maxageclass),sparsed(maxspec,maxageclass)
58  character :: adate*8,atime*6
59
60
61  ! Determine current calendar date, needed for the file name
62  !**********************************************************
63
64  jul=bdate+real(itime,kind=dp)/86400._dp
65  call caldate(jul,jjjjmmdd,ihmmss)
66  write(adate,'(i8.8)') jjjjmmdd
67  write(atime,'(i6.6)') ihmmss
68
69
70  open(unitflux,file=path(2)(1:length(2))//'grid_flux_'//adate// &
71       atime,form='unformatted')
72
73  !**************************************************************
74  ! Check, whether output of full grid or sparse matrix format is
75  ! more efficient in terms of storage space. This is checked for
76  ! every species and for every age class
77  !**************************************************************
78
79  do k=1,nspec
80    do nage=1,nageclass
81      ncellse(k,nage)=0
82      ncellsw(k,nage)=0
83      ncellsn(k,nage)=0
84      ncellss(k,nage)=0
85      ncellsu(k,nage)=0
86      ncellsd(k,nage)=0
87    end do
88  end do
89
90  do k=1,nspec
91  do kp=1,maxpointspec_act
92    do nage=1,nageclass
93      do jy=0,numygrid-1
94        do ix=0,numxgrid-1
95          do kz=1,numzgrid
96            if (flux(2,ix,jy,kz,k,kp,nage).gt.0) ncellse(k,nage)= &
97                 ncellse(k,nage)+1
98            if (flux(1,ix,jy,kz,k,kp,nage).gt.0) ncellsw(k,nage)= &
99                 ncellsw(k,nage)+1
100            if (flux(4,ix,jy,kz,k,kp,nage).gt.0) ncellsn(k,nage)= &
101                 ncellsn(k,nage)+1
102            if (flux(3,ix,jy,kz,k,kp,nage).gt.0) ncellss(k,nage)= &
103                 ncellss(k,nage)+1
104            if (flux(5,ix,jy,kz,k,kp,nage).gt.0) ncellsu(k,nage)= &
105                 ncellsu(k,nage)+1
106            if (flux(6,ix,jy,kz,k,kp,nage).gt.0) ncellsd(k,nage)= &
107                 ncellsd(k,nage)+1
108          end do
109        end do
110      end do
111    end do
112  end do
113  end do
114
115  ! Output in sparse matrix format more efficient, if less than
116  ! 2/5 of all cells contains concentrations>0
117  !************************************************************
118
119  do k=1,nspec
120    do nage=1,nageclass
121      if (4*ncellse(k,nage).lt.numxgrid*numygrid*numzgrid) then
122        sparsee(k,nage)=.true.
123      else
124        sparsee(k,nage)=.false.
125      endif
126      if (4*ncellsw(k,nage).lt.numxgrid*numygrid*numzgrid) then
127        sparsew(k,nage)=.true.
128      else
129        sparsew(k,nage)=.false.
130      endif
131      if (4*ncellsn(k,nage).lt.numxgrid*numygrid*numzgrid) then
132        sparsen(k,nage)=.true.
133      else
134        sparsen(k,nage)=.false.
135      endif
136      if (4*ncellss(k,nage).lt.numxgrid*numygrid*numzgrid) then
137        sparses(k,nage)=.true.
138      else
139        sparses(k,nage)=.false.
140      endif
141      if (4*ncellsu(k,nage).lt.numxgrid*numygrid*numzgrid) then
142        sparseu(k,nage)=.true.
143      else
144        sparseu(k,nage)=.false.
145      endif
146      if (4*ncellsd(k,nage).lt.numxgrid*numygrid*numzgrid) then
147        sparsed(k,nage)=.true.
148      else
149        sparsed(k,nage)=.false.
150      endif
151    end do
152  end do
153
154
155
156  ! Flux output: divide by area and time to get flux in ng/m2/s
157  !************************************************************
158
159  write(unitflux) itime
160  do k=1,nspec
161  do kp=1,maxpointspec_act
162    do nage=1,nageclass
163
164      if (sparsee(k,nage)) then
165        write(unitflux) 1
166        do kz=1,numzgrid
167          do jy=0,numygrid-1
168            do ix=0,numxgrid-1
169              if (flux(2,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
170                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
171                   flux(2,ix,jy,kz,k,kp,nage)/areaeast(ix,jy,kz)/outstep
172            end do
173          end do
174        end do
175        write(unitflux) -999,999.
176      else
177        write(unitflux) 2
178        do kz=1,numzgrid
179          do ix=0,numxgrid-1
180            write(unitflux) (1.e12*flux(2,ix,jy,kz,k,kp,nage)/ &
181                 areaeast(ix,jy,kz)/outstep,jy=0,numygrid-1)
182          end do
183        end do
184      endif
185
186      if (sparsew(k,nage)) then
187        write(unitflux) 1
188        do kz=1,numzgrid
189          do jy=0,numygrid-1
190            do ix=0,numxgrid-1
191              if (flux(1,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
192                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
193                   flux(1,ix,jy,kz,k,kp,nage)/areaeast(ix,jy,kz)/outstep
194            end do
195          end do
196        end do
197        write(unitflux) -999,999.
198      else
199        write(unitflux) 2
200        do kz=1,numzgrid
201          do ix=0,numxgrid-1
202            write(unitflux) (1.e12*flux(1,ix,jy,kz,k,kp,nage)/ &
203                 areaeast(ix,jy,kz)/outstep,jy=0,numygrid-1)
204          end do
205        end do
206      endif
207
208      if (sparses(k,nage)) then
209        write(unitflux) 1
210        do kz=1,numzgrid
211          do jy=0,numygrid-1
212            do ix=0,numxgrid-1
213              if (flux(3,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
214                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
215                   flux(3,ix,jy,kz,k,kp,nage)/areanorth(ix,jy,kz)/outstep
216            end do
217          end do
218        end do
219        write(unitflux) -999,999.
220      else
221        write(unitflux) 2
222        do kz=1,numzgrid
223          do ix=0,numxgrid-1
224            write(unitflux) (1.e12*flux(3,ix,jy,kz,k,kp,nage)/ &
225                 areanorth(ix,jy,kz)/outstep,jy=0,numygrid-1)
226          end do
227        end do
228      endif
229
230      if (sparsen(k,nage)) then
231        write(unitflux) 1
232        do kz=1,numzgrid
233          do jy=0,numygrid-1
234            do ix=0,numxgrid-1 ! north
235              if (flux(4,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
236                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
237                   flux(4,ix,jy,kz,k,kp,nage)/areanorth(ix,jy,kz)/outstep
238            end do
239          end do
240        end do
241        write(unitflux) -999,999.
242      else
243        write(unitflux) 2
244        do kz=1,numzgrid
245          do ix=0,numxgrid-1
246            write(unitflux) (1.e12*flux(4,ix,jy,kz,k,kp,nage)/ &
247                 areanorth(ix,jy,kz)/outstep,jy=0,numygrid-1)
248          end do
249        end do
250      endif
251
252      if (sparseu(k,nage)) then
253        write(unitflux) 1
254        do kz=1,numzgrid
255          do jy=0,numygrid-1
256            do ix=0,numxgrid-1
257              if (flux(5,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
258                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
259                   flux(5,ix,jy,kz,k,kp,nage)/area(ix,jy)/outstep
260            end do
261          end do
262        end do
263        write(unitflux) -999,999.
264      else
265        write(unitflux) 2
266        do kz=1,numzgrid
267          do ix=0,numxgrid-1
268            write(unitflux) (1.e12*flux(5,ix,jy,kz,k,kp,nage)/ &
269                 area(ix,jy)/outstep,jy=0,numygrid-1)
270          end do
271        end do
272      endif
273
274      if (sparsed(k,nage)) then
275        write(unitflux) 1
276        do kz=1,numzgrid
277          do jy=0,numygrid-1
278            do ix=0,numxgrid-1
279              if (flux(6,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
280                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
281                   flux(6,ix,jy,kz,k,kp,nage)/area(ix,jy)/outstep
282            end do
283          end do
284        end do
285        write(unitflux) -999,999.
286      else
287        write(unitflux) 2
288        do kz=1,numzgrid
289          do ix=0,numxgrid-1
290            write(unitflux) (1.e12*flux(6,ix,jy,kz,k,kp,nage)/ &
291                 area(ix,jy)/outstep,jy=0,numygrid-1)
292          end do
293        end do
294      endif
295
296    end do
297  end do
298  end do
299
300
301  close(unitflux)
302
303
304  ! Reinitialization of grid
305  !*************************
306
307  do k=1,nspec
308  do kp=1,maxpointspec_act
309    do jy=0,numygrid-1
310      do ix=0,numxgrid-1
311          do kz=1,numzgrid
312            do nage=1,nageclass
313              do i=1,6
314                flux(i,ix,jy,kz,k,kp,nage)=0.
315              end do
316            end do
317          end do
318      end do
319    end do
320  end do
321  end do
322
323
324end subroutine fluxoutput
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG