source: flexpart.git/src/concoutput_nest_mpi.f90 @ 8a65cb0

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

Added code, makefile for dev branch

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