files for calculating averages in partpos are added

src/partoutput_average.f90

+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+subroutine partoutput_average(itime)
+  !                             i
+  !*****************************************************************************
+  !                                                                            *
+  !     Dump all particle positions                                            *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     12 March 1999                                                          *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  !                                                                            *
+  !*****************************************************************************
+
+  use par_mod
+  use com_mod
+  use netcdf_output_mod 
+
+  implicit none
+
+  real(kind=dp) :: jul
+  integer :: itime,i,j,jjjjmmdd,ihmmss
+  integer :: ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp,indpart,indts
+  real :: xlon,ylat
+  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
+  real :: topo,hm(2),hmixi,pv1(2),pvprof(2),pvi,qv1(2),qvprof(2),qvi
+  real :: tt1(2),ttprof(2),tti,rho1(2),rhoprof(2),rhoi
+  real :: tr(2),tri,zlim
+  character :: adate*8,atime*6
+
+  integer(kind=2) :: ishort_xlon(maxpart),ishort_ylat(maxpart),ishort_z(maxpart)
+  integer(kind=2) :: ishort_topo(maxpart),ishort_tro(maxpart),ishort_hmix(maxpart)
+  integer(kind=2) :: ishort_pv(maxpart),ishort_rho(maxpart),ishort_qv(maxpart)
+  integer(kind=2) :: ishort_tt(maxpart),ishort_uu(maxpart),ishort_vv(maxpart)
+  integer(kind=2) :: ishort_energy(maxpart)
+
+! character(len=255)  :: fname
+! integer             :: cache_size, ncid
+! integer npointDimID,topoID,zID,xID,yID,pvID,qvID,ttID,uuID,vvID,rhoID,troID,hmixID,itimeID
+
+  ! Determine current calendar date, needed for the file name
+  !**********************************************************
+
+  jul=bdate+real(itime,kind=dp)/86400._dp
+  call caldate(jul,jjjjmmdd,ihmmss)
+  write(adate,'(i8.8)') jjjjmmdd
+  write(atime,'(i6.6)') ihmmss
+
+
+  ! Some variables needed for temporal interpolation
+  !*************************************************
+
+  dt1=real(itime-memtime(1))
+  dt2=real(memtime(2)-itime)
+  dtt=1./(dt1+dt2)
+
+  ! Write current time to file
+  !***************************
+
+  do indpart=1,numpart
+
+  ! Take only valid particles
+  !**************************
+
+    if (itra1(indpart).eq.itime) then
+      i=npoint(indpart)
+      if (npart_av(i).gt.0) then
+        indts=part_av_filled(i)+1
+        if (i.lt.10) then
+            write (*,*) indts,i,indpart,npart_av(i),part_av_tt(i,indts),part_av_tt(i,indts)
+        endif
+        part_av_topo(i,indts)=part_av_topo(i,indts)/float(npart_av(i))
+        part_av_z(i,indts)=part_av_z(i,indts)/float(npart_av(i))
+        part_av_pv(i,indts)=part_av_pv(i,indts)/float(npart_av(i))
+        part_av_qv(i,indts)=part_av_qv(i,indts)/float(npart_av(i))
+        part_av_tt(i,indts)=part_av_tt(i,indts)/float(npart_av(i))
+        part_av_uu(i,indts)=part_av_uu(i,indts)/float(npart_av(i))
+        part_av_vv(i,indts)=part_av_vv(i,indts)/float(npart_av(i))
+        part_av_rho(i,indts)=part_av_rho(i,indts)/float(npart_av(i))
+        part_av_tro(i,indts)=part_av_tro(i,indts)/float(npart_av(i))
+        part_av_hmix(i,indts)=part_av_hmix(i,indts)/float(npart_av(i))
+!        part_av_energy(i,indts)=part_av_energy(i,indts)/float(npart_av(i))
+
+        part_av_cartx(i)=part_av_cartx(i)/float(npart_av(i))
+        part_av_carty(i)=part_av_carty(i)/float(npart_av(i))
+        part_av_cartz(i)=part_av_cartz(i)/float(npart_av(i))
+      
+        part_av_itime(i,indts)=itramem(i) !maybe needs to be averaged!
+
+! Project Cartesian coordinates back onto Earth's surface
+! #######################################################
+        xlon=atan2(part_av_cartx(i),-1.*part_av_carty(i))
+        ylat=atan2(part_av_cartz(i),sqrt(part_av_cartx(i)*part_av_cartx(i)+ &
+           part_av_carty(i)*part_av_carty(i)))
+        xlon=xlon/pi180
+        ylat=ylat/pi180
+
+
+        if (xlon.gt.180.) xlon=xlon-360.
+        if (xlon.lt.-180.) xlon=xlon+360.
+        part_av_x(i,indts)=xlon*180.
+        part_av_y(i,indts)=ylat*360.
+
+        part_av_filled(i)=part_av_filled(i)+1
+      endif
+    endif
+
+  end do
+
+  ! Write the output
+  !*****************
+!   write(*,*) 'partoutput_average',numpart
+!   call writeaveragepart_netcdf(adate,atime)
+
+! Re-initialize averages, and set number of steps to zero
+   do i=1,maxpart
+     npart_av(i)=0
+!    part_av_topo(i)=0.
+!    part_av_x(i)=0.
+!    part_av_y(i)=0.
+!    part_av_z(i)=0.
+!    part_av_cartx(i)=0.
+!    part_av_carty(i)=0.
+!    part_av_cartz(i)=0.
+!    part_av_pv(i)=0.
+!    part_av_qv(i)=0.
+!    part_av_tt(i)=0.
+!    part_av_uu(i)=0.
+!    part_av_vv(i)=0.
+!    part_av_rho(i)=0.
+!    part_av_tro(i)=0.
+!    part_av_hmix(i)=0.
+!    part_av_energy(i)=0.
+!    part_av_itime(i)=0.
+   enddo
+  
+end subroutine partoutput_average