source: flexpart.git/src/concoutput_nest_mpi.f90 @ f9ce123

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

Updated wet depo scheme

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