source: flexpart.git/src/concoutput_nest_mpi.f90 @ 3481cc1

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

move license from headers to a different file

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