source: flexpart.git/src/fluxoutput.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

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