source: flexpart.git/src/concoutput_nest_mpi.f90

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

add SPDX-License-Identifier to all .f90 files

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