Forgot to add new files on previous commit

src/FLEXPART_MPI.f90

@@ -412,7 +412,7 @@ program flexpart
     if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
   end if ! (mpif_pid == 0) 
 
-  if (mp_measure_time) call mpif_mtime('iotime',0)
+  if (mp_measure_time) call mpif_mtime('iotime',1)
 
   if (verbosity.gt.0 .and. lroot) then
     print*,'call openreceptors'

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
+
+  implicit none
+
+  real(kind=dp) :: jul
+  integer :: itime,i,j,jjjjmmdd,ihmmss
+  integer :: ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
+  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)
+
+
+  ! 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)
+
+  ! Open output file and write the output
+  !**************************************
+
+    open(unitpartout_average,file=path(2)(1:length(2))//'partposit_average_'//adate// &
+         atime,form='unformatted',access='direct',status='replace',recl=24)
+
+
+  ! Write current time to file
+  !***************************
+
+!  write(unitpartout_average) itime,numpart
+  do i=1,numpart
+
+  ! Take only valid particles
+  !**************************
+
+    if (itra1(i).eq.itime) then
+      part_av_topo(i)=part_av_topo(i)/float(npart_av(i))
+      part_av_z(i)=part_av_z(i)/float(npart_av(i))
+      part_av_pv(i)=part_av_pv(i)/float(npart_av(i))
+      part_av_qv(i)=part_av_qv(i)/float(npart_av(i))
+      part_av_tt(i)=part_av_tt(i)/float(npart_av(i))
+      part_av_uu(i)=part_av_uu(i)/float(npart_av(i))
+      part_av_vv(i)=part_av_vv(i)/float(npart_av(i))
+      part_av_rho(i)=part_av_rho(i)/float(npart_av(i))
+      part_av_tro(i)=part_av_tro(i)/float(npart_av(i))
+      part_av_hmix(i)=part_av_hmix(i)/float(npart_av(i))
+      part_av_energy(i)=part_av_energy(i)/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))
+
+! 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
+
+
+  ! Convert all data to integer*2 variables (from -32768 to 32767) for compressed output
+  !*************************************************************************************
+
+      if (xlon.gt.180.) xlon=xlon-360.
+      if (xlon.lt.-180.) xlon=xlon+360.
+      ishort_xlon(i)=nint(xlon*180.)
+      ishort_ylat(i)=nint(ylat*360.)
+
+      zlim=(part_av_z(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_z(i)=nint(zlim)
+
+      zlim=(part_av_topo(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_topo(i)=nint(zlim)
+
+      zlim=(part_av_tro(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_tro(i)=nint(zlim)
+
+      zlim=(part_av_hmix(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_hmix(i)=nint(zlim)
+
+      zlim=(part_av_rho(i)*20000.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_rho(i)=nint(zlim)
+
+      zlim=(part_av_qv(i)*1000000.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_qv(i)=nint(zlim)
+
+      zlim=(part_av_pv(i)*100.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_pv(i)=nint(zlim)
+
+      zlim=((part_av_tt(i)-273.15))*300.
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_tt(i)=nint(zlim)
+
+      zlim=(part_av_uu(i)*200.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_uu(i)=nint(zlim)
+
+      zlim=(part_av_vv(i)*200.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_vv(i)=nint(zlim)
+
+      zlim=(part_av_energy(i)-300000.)/30.
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_energy(i)=nint(zlim)
+
+  ! Turn on for readable test output
+  !*********************************
+
+!        write(119,*) itime,i,xlon,ylat,part_av_z(i),part_av_topo(i),part_av_tro(i), &
+!        part_av_hmix(i),part_av_rho(i),part_av_qv(i),part_av_pv(i),part_av_tt(i), &
+!        ishort_uu(i),ishort_vv(i)
+    endif
+
+! Re-initialize averages, and set number of steps to zero
+    npart_av(i)=0
+    part_av_topo(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.
+    
+  end do
+
+  ! Write the output
+  !*****************
+  do i=1,numpart
+    if (itra1(i).eq.itime) then
+      write(unitpartout_average,rec=i) ishort_xlon(i),ishort_ylat(i),ishort_z(i), &
+           ishort_topo(i),ishort_tro(i),ishort_hmix(i),ishort_rho(i),ishort_qv(i),ishort_pv(i), &
+           ishort_tt(i),ishort_uu(i),ishort_vv(i)  ! ,ishort_energy(i)
+    endif
+  enddo
+
+
+  close(unitpartout_average)
+
+end subroutine partoutput_average

src/partoutput_average_mpi.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,irec)
+  !                             i
+  !*****************************************************************************
+  !                                                                            *
+  !     Dump all particle positions                                            *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     12 March 1999                                                          *
+  !     03/2019 AST: Version for MPI                                           *
+  !                  processes sequentially access and append data to file     *
+  !                                                                            *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  !                                                                            *
+  !*****************************************************************************
+
+  use par_mod
+  use com_mod
+  use mpi_mod
+
+
+  implicit none
+
+  real(kind=dp) :: jul
+  integer :: itime,i,j,jjjjmmdd,ihmmss,irec
+  integer :: ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
+  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
+  character(LEN=8) :: file_stat='OLD'
+
+  integer(kind=2) :: ishort_xlon(maxpart_mpi),ishort_ylat(maxpart_mpi),ishort_z(maxpart_mpi)
+  integer(kind=2) :: ishort_topo(maxpart_mpi),ishort_tro(maxpart_mpi),ishort_hmix(maxpart_mpi)
+  integer(kind=2) :: ishort_pv(maxpart_mpi),ishort_rho(maxpart_mpi),ishort_qv(maxpart_mpi)
+  integer(kind=2) :: ishort_tt(maxpart_mpi),ishort_uu(maxpart_mpi),ishort_vv(maxpart_mpi)
+  integer(kind=2) :: ishort_energy(maxpart_mpi)
+
+  ! MPI root process creates/overwrites the file, other processes append to it
+  if (lroot) file_stat='REPLACE'
+
+  ! 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)
+
+  ! Open output file and write the output
+  !**************************************
+
+  open(unitpartout_average,file=path(2)(1:length(2))//'partposit_average_'//adate// &
+       atime,form='unformatted',access='direct',status=file_stat,recl=24)
+
+
+  ! Write current time to file
+  !***************************
+
+!  write(unitpartout_average) itime,numpart
+  do i=1,numpart
+
+  ! Take only valid particles
+  !**************************
+
+    if (itra1(i).eq.itime) then
+      part_av_topo(i)=part_av_topo(i)/float(npart_av(i))
+      part_av_z(i)=part_av_z(i)/float(npart_av(i))
+      part_av_pv(i)=part_av_pv(i)/float(npart_av(i))
+      part_av_qv(i)=part_av_qv(i)/float(npart_av(i))
+      part_av_tt(i)=part_av_tt(i)/float(npart_av(i))
+      part_av_uu(i)=part_av_uu(i)/float(npart_av(i))
+      part_av_vv(i)=part_av_vv(i)/float(npart_av(i))
+      part_av_rho(i)=part_av_rho(i)/float(npart_av(i))
+      part_av_tro(i)=part_av_tro(i)/float(npart_av(i))
+      part_av_hmix(i)=part_av_hmix(i)/float(npart_av(i))
+      part_av_energy(i)=part_av_energy(i)/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))
+
+! 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
+
+
+  ! Convert all data to integer*2 variables (from -32768 to 32767) for compressed output
+  !*************************************************************************************
+
+      if (xlon.gt.180.) xlon=xlon-360.
+      if (xlon.lt.-180.) xlon=xlon+360.
+      ishort_xlon(i)=nint(xlon*180.)
+      ishort_ylat(i)=nint(ylat*360.)
+
+      zlim=(part_av_z(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_z(i)=nint(zlim)
+
+      zlim=(part_av_topo(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_topo(i)=nint(zlim)
+
+      zlim=(part_av_tro(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_tro(i)=nint(zlim)
+
+      zlim=(part_av_hmix(i)*2.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_hmix(i)=nint(zlim)
+
+      zlim=(part_av_rho(i)*20000.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_rho(i)=nint(zlim)
+
+      zlim=(part_av_qv(i)*1000000.-32000.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_qv(i)=nint(zlim)
+
+      zlim=(part_av_pv(i)*100.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_pv(i)=nint(zlim)
+
+      zlim=((part_av_tt(i)-273.15))*300.
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_tt(i)=nint(zlim)
+
+      zlim=(part_av_uu(i)*200.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_uu(i)=nint(zlim)
+
+      zlim=(part_av_vv(i)*200.)
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_vv(i)=nint(zlim)
+
+      zlim=(part_av_energy(i)-300000.)/30.
+      zlim=min(zlim,32766.)
+      zlim=max(zlim,-32766.)
+      ishort_energy(i)=nint(zlim)
+
+
+  ! Turn on for readable test output
+  !*********************************
+
+!        write(119,*) itime,i,xlon,ylat,part_av_z(i),part_av_topo(i),part_av_tro(i), &
+!        part_av_hmix(i),part_av_rho(i),part_av_qv(i),part_av_pv(i),part_av_tt(i), &
+!        ishort_uu(i),ishort_vv(i)
+    endif
+
+! Re-initialize averages, and set number of steps to zero
+    npart_av(i)=0
+    part_av_topo(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.
+    
+  end do
+
+  ! Write the output
+  !*****************
+
+  do i=1,numpart
+    if (itra1(i).eq.itime) then
+      write(unitpartout_average,rec=irec+i) ishort_xlon(i),ishort_ylat(i),ishort_z(i), &
+           ishort_topo(i),ishort_tro(i),ishort_hmix(i),ishort_rho(i),ishort_qv(i),ishort_pv(i), &
+           ishort_tt(i),ishort_uu(i),ishort_vv(i)  ! ,ishort_energy(i)
+    endif
+  enddo
+
+  close(unitpartout_average)
+
+end subroutine partoutput_average

src/partpos_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 partpos_average(itime,j)
+
+
+!**********************************************************************
+! This subroutine averages particle quantities, to be used for particle
+! dump (in partoutput.f90). Averaging is done over output interval.
+!**********************************************************************
+
+  use par_mod
+  use com_mod
+
+  implicit none
+
+  integer :: itime,j,ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
+  real :: xlon,ylat,x,y,z
+  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 :: uu1(2),uuprof(2),uui,vv1(2),vvprof(2),vvi
+  real :: tr(2),tri,energy
+
+
+
+ ! Some variables needed for temporal interpolation
+  !*************************************************
+
+  dt1=real(itime-memtime(1))
+  dt2=real(memtime(2)-itime)
+  dtt=1./(dt1+dt2)
+
+  xlon=xlon0+xtra1(j)*dx
+  ylat=ylat0+ytra1(j)*dy
+
+  !*****************************************************************************
+  ! Interpolate several variables (PV, specific humidity, etc.) to particle position
+  !*****************************************************************************
+
+  ix=xtra1(j)
+  jy=ytra1(j)
+  ixp=ix+1
+  jyp=jy+1
+  ddx=xtra1(j)-real(ix)
+  ddy=ytra1(j)-real(jy)
+  rddx=1.-ddx
+  rddy=1.-ddy
+  p1=rddx*rddy
+  p2=ddx*rddy
+  p3=rddx*ddy
+  p4=ddx*ddy
+
+! eso: Temporary fix for particle exactly at north pole
+  if (jyp >= nymax) then
+  !  write(*,*) 'WARNING: conccalc.f90 jyp >= nymax'
+    jyp=jyp-1
+  end if
+
+  ! Topography
+  !***********
+
+  topo=p1*oro(ix,jy)+p2*oro(ixp,jy)+p3*oro(ix,jyp)+p4*oro(ixp,jyp)
+
+ ! Potential vorticity, specific humidity, temperature, and density
+  !*****************************************************************
+
+  do il=2,nz
+    if (height(il).gt.ztra1(j)) then
+      indz=il-1
+      indzp=il
+      goto 6
+    endif
+  end do
+6 continue
+
+  dz1=ztra1(j)-height(indz)
+  dz2=height(indzp)-ztra1(j)
+  dz=1./(dz1+dz2)
+
+
+  do ind=indz,indzp
+    do m=1,2
+      indexh=memind(m)
+
+  ! Potential vorticity
+      pv1(m)=p1*pv(ix ,jy ,ind,indexh) &
+            +p2*pv(ixp,jy ,ind,indexh) &
+            +p3*pv(ix ,jyp,ind,indexh) &
+            +p4*pv(ixp,jyp,ind,indexh)
+  ! Specific humidity
+      qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
+            +p2*qv(ixp,jy ,ind,indexh) &
+            +p3*qv(ix ,jyp,ind,indexh) &
+            +p4*qv(ixp,jyp,ind,indexh)
+  ! Temperature
+      tt1(m)=p1*tt(ix ,jy ,ind,indexh) &
+            +p2*tt(ixp,jy ,ind,indexh) &
+            +p3*tt(ix ,jyp,ind,indexh) &
+            +p4*tt(ixp,jyp,ind,indexh)
+  ! U wind
+      uu1(m)=p1*uu(ix ,jy ,ind,indexh) &
+            +p2*uu(ixp,jy ,ind,indexh) &
+            +p3*uu(ix ,jyp,ind,indexh) &
+            +p4*uu(ixp,jyp,ind,indexh)
+  ! V wind
+      vv1(m)=p1*vv(ix ,jy ,ind,indexh) &
+            +p2*vv(ixp,jy ,ind,indexh) &
+            +p3*vv(ix ,jyp,ind,indexh) &
+            +p4*vv(ixp,jyp,ind,indexh)
+  ! Density
+      rho1(m)=p1*rho(ix ,jy ,ind,indexh) &
+             +p2*rho(ixp,jy ,ind,indexh) &
+             +p3*rho(ix ,jyp,ind,indexh) &
+             +p4*rho(ixp,jyp,ind,indexh)
+    end do
+    pvprof(ind-indz+1)=(pv1(1)*dt2+pv1(2)*dt1)*dtt
+    qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt
+    ttprof(ind-indz+1)=(tt1(1)*dt2+tt1(2)*dt1)*dtt
+    uuprof(ind-indz+1)=(uu1(1)*dt2+uu1(2)*dt1)*dtt
+    vvprof(ind-indz+1)=(vv1(1)*dt2+vv1(2)*dt1)*dtt
+    rhoprof(ind-indz+1)=(rho1(1)*dt2+rho1(2)*dt1)*dtt
+  end do
+  pvi=(dz1*pvprof(2)+dz2*pvprof(1))*dz
+  qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz
+  tti=(dz1*ttprof(2)+dz2*ttprof(1))*dz
+  uui=(dz1*uuprof(2)+dz2*uuprof(1))*dz
+  vvi=(dz1*vvprof(2)+dz2*vvprof(1))*dz
+  rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
+
+  ! Tropopause and PBL height
+  !**************************
+
+  do m=1,2
+    indexh=memind(m)
+
+  ! Tropopause
+    tr(m)=p1*tropopause(ix ,jy ,1,indexh) &
+        + p2*tropopause(ixp,jy ,1,indexh) &
+        + p3*tropopause(ix ,jyp,1,indexh) &
+        + p4*tropopause(ixp,jyp,1,indexh)
+
+  ! PBL height
+    hm(m)=p1*hmix(ix ,jy ,1,indexh) &
+        + p2*hmix(ixp,jy ,1,indexh) &
+        + p3*hmix(ix ,jyp,1,indexh) &
+        + p4*hmix(ixp,jyp,1,indexh)
+  end do
+
+  hmixi=(hm(1)*dt2+hm(2)*dt1)*dtt
+  tri=(tr(1)*dt2+tr(2)*dt1)*dtt
+
+
+  energy=tti*cpa+(ztra1(j)+topo)*9.81+qvi*2501000.+(uui**2+vvi**2)/2.
+
+  ! Add new values to sum and increase counter by one
+  !**************************************************
+
+  npart_av(j)=npart_av(j)+1
+
+  ! Calculate Cartesian 3D coordinates suitable for averaging
+  !**********************************************************
+
+  xlon=xlon*pi180
+  ylat=ylat*pi180
+  x = cos(ylat)*sin(xlon)
+  y = -1.*cos(ylat)*cos(xlon)