source: flexpart.git/src/concoutput_nest_mpi.f90 @ 38b7917

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

Parallelization of domain fill option (save/restart not implemented yet)

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