source: flexpart.git/src/concoutput_nest_mpi.f90 @ 02095e3

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

Reworked the domain-filling option (MPI). Fixed a slow loop which had errors in loop counter (MPI)

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