*
L. Bakels 2022, This module contains routines that are relevant for the * plume trajectory computations * *
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt ! SPDX-License-Identifier: GPL-3.0-or-later !***************************************************************************** ! * ! L. Bakels 2022, This module contains routines that are relevant for the * ! plume trajectory computations * ! * !***************************************************************************** module plume_mod implicit none private :: centerofmass,clustering,distance,distance2 public :: plumetraj,openouttraj contains subroutine plumetraj(itime) ! i !***************************************************************************** ! * ! Determines a plume centroid trajectory for each release site, and manages * ! clustering of particle locations. Certain parameters (average PV, * ! tropopause height, etc., are provided along the plume trajectories. * ! At the end, output is written to file 'trajectories.txt'. * ! * ! Author: A. Stohl * ! * ! 24 January 2002 * ! * ! Variables: * ! fclust fraction of particles belonging to each cluster * ! hmixcenter mean mixing height for all particles * ! ncluster number of clusters to be used * ! pvcenter mean PV for all particles * ! pvfract fraction of particles with PV<2pvu * ! rms total horizontal rms distance after clustering * ! rmsdist total horizontal rms distance before clustering * ! rmsclust horizontal rms distance for each individual cluster * ! topocenter mean topography underlying all particles * ! tropocenter mean tropopause height at the positions of particles * ! tropofract fraction of particles within the troposphere * ! zrms total vertical rms distance after clustering * ! zrmsdist total vertical rms distance before clustering * ! xclust,yclust, Cluster centroid positions * ! zclust * ! * !***************************************************************************** use point_mod use par_mod use com_mod use mean_mod use particle_mod #ifdef ETA use coord_ecmwf_mod #endif use windfields_mod implicit none integer :: ii,itime,ix,jy,ixp,jyp,indexh,i,j,k,m,n,il,ind,indz,indzp ! real :: xl(maxpart),yl(maxpart),zl(maxpart) ! moved to particle_mod and now xplum,yplum,zplum real :: xcenter,ycenter,zcenter,dist,rmsdist,zrmsdist real :: xclust(ncluster),yclust(ncluster),zclust(ncluster) real :: fclust(ncluster),rms,rmsclust(ncluster),zrms real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz real :: topo,topocenter,hm(2),hmixi,hmixfract,hmixcenter real :: pv1(2),pvprof(2),pvi,pvcenter,pvfract,tr(2),tri,tropofract real :: tropocenter dt1=real(itime-memtime(1)) dt2=real(memtime(2)-itime) dtt=1./(dt1+dt2) ! Loop about all release points !****************************** do j=1,numpoint if (abs(ireleasestart(j)-itime).gt.lage(nageclass)) cycle topocenter=0. hmixcenter=0. hmixfract=0. tropocenter=0. tropofract=0. pvfract=0. pvcenter=0. rmsdist=0. zrmsdist=0. n=0 do ii=1,count%alive i=count%ialive(ii) if (part(i)%npoint.ne.j) cycle n=n+1 xplum(n)=xlon0+real(part(i)%xlon)*dx yplum(n)=ylat0+real(part(i)%ylat)*dy #ifdef ETA call update_zeta_to_z(itime,i) #endif zplum(n)=real(part(i)%z) ! Interpolate PBL height, PV, and tropopause height to each ! particle position in order to determine fraction of particles ! within the PBL, above tropopause height, and average PV. ! Interpolate topography, too, and convert to altitude asl !************************************************************** ix=int(part(i)%xlon) jy=int(part(i)%ylat) ixp=ix+1 jyp=jy+1 ! eso: Temporary fix for particle exactly at north pole if (jyp >= nymax) then write(*,*) 'WARNING: plume_mod.f90 jyp >= nymax. xt,yt:',part(i)%xlon,part(i)%ylat jyp=jyp-1 end if if (ixp >= nxmax) then write(*,*) 'WARNING: plume_mod.f90 ixp >= nxmax. xt,yt:',part(i)%xlon,part(i)%ylat ixp=ixp-nxmax end if ddx=real(part(i)%xlon)-real(ix) ddy=real(part(i)%ylat)-real(jy) rddx=1.-ddx rddy=1.-ddy p1=rddx*rddy p2=ddx*rddy p3=rddx*ddy p4=ddx*ddy ! Topography !*********** topo=p1*oro(ix ,jy) & + p2*oro(ixp,jy) & + p3*oro(ix ,jyp) & + p4*oro(ixp,jyp) topocenter=topocenter+topo ! Potential vorticity !******************** indz=nz-1 indzp=nz do il=2,nz if (height(il).gt.zplum(n)) then indz=il-1 indzp=il exit endif end do dz1=zplum(n)-height(indz) dz2=height(indzp)-zplum(n) dz=1./(dz1+dz2) do ind=indz,indzp do m=1,2 indexh=memind(m) 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) end do pvprof(ind-indz+1)=(pv1(1)*dt2+pv1(2)*dt1)*dtt end do pvi=(dz1*pvprof(2)+dz2*pvprof(1))*dz pvcenter=pvcenter+pvi if (yplum(n).gt.0.) then if (pvi.lt.2.) pvfract=pvfract+1. else if (pvi.gt.-2.) pvfract=pvfract+1. endif ! Tropopause and PBL height !************************** do m=1,2 indexh=memind(m) 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) 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 if (zplum(n).lt.tri) tropofract=tropofract+1. tropocenter=tropocenter+tri+topo if (zplum(n).lt.hmixi) hmixfract=hmixfract+1. zplum(n)=zplum(n)+topo ! convert to height asl hmixcenter=hmixcenter+hmixi end do ! Make statistics for all plumes with n>0 particles !************************************************** if (n.gt.0) then topocenter=topocenter/real(n) hmixcenter=hmixcenter/real(n) pvcenter=pvcenter/real(n) tropocenter=tropocenter/real(n) hmixfract=100.*hmixfract/real(n) pvfract=100.*pvfract/real(n) tropofract=100.*tropofract/real(n) ! Cluster the particle positions !******************************* call clustering(n,xclust,yclust,zclust,fclust,rms, & rmsclust,zrms) ! Determine center of mass position on earth and average height !************************************************************** call centerofmass(xplum,yplum,n,xcenter,ycenter) call mean(zplum,zcenter,zrmsdist,n) ! Root mean square distance from center of mass !********************************************** do k=1,n dist=distance(yplum(k),xplum(k),ycenter,xcenter) rmsdist=rmsdist+dist*dist end do if (rmsdist.gt.0.) rmsdist=sqrt(rmsdist/real(n)) rmsdist=max(rmsdist,0.) ! Write out results in trajectory data file !****************************************** write(unitouttraj,'(i5,i8,2f9.4,4f8.1,f8.2,4f8.1,3f6.1,& &5(2f8.3,f7.0,f6.1,f8.1))')& &j,itime-(ireleasestart(j)+ireleaseend(j))/2, & xcenter,ycenter,zcenter,topocenter,hmixcenter,tropocenter, & pvcenter,rmsdist,rms,zrmsdist,zrms,hmixfract,pvfract, & tropofract, & (xclust(k),yclust(k),zclust(k),fclust(k),rmsclust(k), & k=1,ncluster) endif end do end subroutine plumetraj subroutine centerofmass(xl,yl,n,xcenter,ycenter) ! i i i o o !***************************************************************************** ! * ! This routine calculates the center of mass of n points on the Earth. * ! Input are the longitudes (xl) and latitudes (yl) of the individual * ! points, output is the longitude and latitude of the centre of mass. * ! * ! Author: A. Stohl * ! * ! 24 January 2002 * ! * !***************************************************************************** use par_mod implicit none integer :: n,l real :: xl(n),yl(n),xll,yll,xav,yav,zav,x,y,z,xcenter,ycenter xav=0. yav=0. zav=0. do l=1,n ! Convert longitude and latitude from degrees to radians !******************************************************* xll=xl(l)*pi180 yll=yl(l)*pi180 ! Calculate 3D coordinates from longitude and latitude !***************************************************** x = cos(yll)*sin(xll) y = -1.*cos(yll)*cos(xll) z = sin(yll) ! Find the mean location in Cartesian coordinates !************************************************ xav=xav+x yav=yav+y zav=zav+z end do xav=xav/real(n) yav=yav/real(n) zav=zav/real(n) ! Project the point back onto Earth's surface !******************************************** xcenter=atan2(xav,-1.*yav) ycenter=atan2(zav,sqrt(xav*xav+yav*yav)) ! Convert back to degrees !************************ xcenter=xcenter/pi180 ycenter=ycenter/pi180 end subroutine centerofmass subroutine clustering(n,xclust,yclust,zclust,fclust,rms, & rmsclust,zrms) ! i i i i o o o o o ! o o !***************************************************************************** ! * ! This routine clusters the particle position into ncluster custers. * ! Input are the longitudes (xl) and latitudes (yl) of the individual * ! points, output are the cluster mean positions (xclust,yclust). * ! Vertical positions are not directly used for the clustering. * ! * ! For clustering, the procedure described in Dorling et al. (1992) is used.* ! * ! Dorling, S.R., Davies, T.D. and Pierce, C.E. (1992): * ! Cluster analysis: a technique for estimating the synoptic meteorological * ! controls on air and precipitation chemistry - method and applications. * ! Atmospheric Environment 26A, 2575-2581. * ! * ! * ! Author: A. Stohl * ! * ! 1 February 2002 * ! * ! Variables: * ! fclust fraction of particles belonging to each cluster * ! ncluster number of clusters to be used * ! rms total horizontal rms distance after clustering * ! rmsclust horizontal rms distance for each individual cluster * ! zrms total vertical rms distance after clustering * ! xclust,yclust, Cluster centroid positions * ! zclust * ! xl,yl,zl particle positions * ! * !***************************************************************************** use par_mod use particle_mod implicit none integer :: n,i,j,l,numb(ncluster),ncl real :: xclust(ncluster),yclust(ncluster),x,y,z real :: zclust(ncluster),distances,distancemin,rms,rmsold real :: xav(ncluster),yav(ncluster),zav(ncluster),fclust(ncluster) real :: rmsclust(ncluster) real :: zdist,zrms if (n.lt.ncluster) return rmsold=-5. ! Convert longitude and latitude from degrees to radians !******************************************************* do i=1,n nclust(i)=i xplum(i)=xplum(i)*pi180 yplum(i)=yplum(i)*pi180 end do ! Generate a seed for each cluster !********************************* do j=1,ncluster zclust(j)=0. xclust(j)=xplum(j*n/ncluster) yclust(j)=yplum(j*n/ncluster) end do ! Iterative loop to compute the cluster means !******************************************** do l=1,100 ! Assign each particle to a cluster: criterion minimum distance to the ! cluster mean position !********************************************************************* do i=1,n distancemin=10.**10. ncl=1 do j=1,ncluster distances=distance2(yplum(i),xplum(i),yclust(j),xclust(j)) if (distances.lt.distancemin) then distancemin=distances ncl=j endif end do nclust(i)=ncl end do ! Recalculate the cluster centroid position: convert to 3D Cartesian coordinates, ! calculate mean position, and re-project this point onto the Earth's surface !***************************************************************************** do j=1,ncluster xav(j)=0. yav(j)=0. zav(j)=0. rmsclust(j)=0. numb(j)=0 end do rms=0. do i=1,n numb(nclust(i))=numb(nclust(i))+1 distances=distance2(yplum(i),xplum(i), & yclust(nclust(i)),xclust(nclust(i))) ! rms is the total rms of all particles ! rmsclust is the rms for a particular cluster !********************************************* rms=rms+distances*distances rmsclust(nclust(i))=rmsclust(nclust(i))+distances*distances ! Calculate Cartesian 3D coordinates from longitude and latitude !*************************************************************** x = cos(yplum(i))*sin(xplum(i)) y = -1.*cos(yplum(i))*cos(xplum(i)) z = sin(yplum(i)) xav(nclust(i))=xav(nclust(i))+x yav(nclust(i))=yav(nclust(i))+y zav(nclust(i))=zav(nclust(i))+z end do rms=sqrt(rms/real(n)) ! Find the mean location in Cartesian coordinates !************************************************ do j=1,ncluster if (numb(j).gt.0) then rmsclust(j)=sqrt(rmsclust(j)/real(numb(j))) xav(j)=xav(j)/real(numb(j)) yav(j)=yav(j)/real(numb(j)) zav(j)=zav(j)/real(numb(j)) ! Project the point back onto Earth's surface !******************************************** xclust(j)=atan2(xav(j),-1.*yav(j)) yclust(j)=atan2(zav(j),sqrt(xav(j)*xav(j)+yav(j)*yav(j))) endif end do ! Leave the loop if the RMS distance decreases only slightly between 2 iterations !***************************************************************************** if ((l.gt.1).and.(abs(rms-rmsold)/rmsold.lt.0.005)) exit rmsold=rms end do ! Convert longitude and latitude from radians to degrees !******************************************************* do i=1,n xplum(i)=xplum(i)/pi180 yplum(i)=yplum(i)/pi180 zclust(nclust(i))=zclust(nclust(i))+zplum(i) end do do j=1,ncluster xclust(j)=xclust(j)/pi180 yclust(j)=yclust(j)/pi180 if (numb(j).gt.0) zclust(j)=zclust(j)/real(numb(j)) fclust(j)=100.*real(numb(j))/real(n) end do ! Determine total vertical RMS deviation !*************************************** zrms=0. do i=1,n zdist=zplum(i)-zclust(nclust(i)) zrms=zrms+zdist*zdist end do if (zrms.gt.0.) zrms=sqrt(zrms/real(n)) end subroutine clustering real function distance(rlat1,rlon1,rlat2,rlon2) !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: GCDIST COMPUTE GREAT CIRCLE DISTANCE ! PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 ! ! ABSTRACT: THIS SUBPROGRAM COMPUTES GREAT CIRCLE DISTANCE ! BETWEEN TWO POINTS ON THE EARTH. ! ! PROGRAM HISTORY LOG: ! 96-04-10 IREDELL ! ! USAGE: ...GCDIST(RLAT1,RLON1,RLAT2,RLON2) ! ! INPUT ARGUMENT LIST: !rlat1 - REAL LATITUDE OF POINT 1 IN DEGREES !rlon1 - REAL LONGITUDE OF POINT 1 IN DEGREES !rlat2 - REAL LATITUDE OF POINT 2 IN DEGREES !rlon2 - REAL LONGITUDE OF POINT 2 IN DEGREES ! ! OUTPUT ARGUMENT LIST: !distance - REAL GREAT CIRCLE DISTANCE IN KILOMETERS ! ! ATTRIBUTES: ! LANGUAGE: Fortran 90 ! !$$$ use par_mod, only: dp implicit none real :: rlat1,rlon1,rlat2,rlon2 real(kind=dp) :: clat1,clat2,slat1,slat2,cdlon,crd real(kind=dp),parameter :: rerth=6.3712e6_dp real(kind=dp),parameter :: pi=3.14159265358979_dp, dpr=180.0_dp/pi ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - if ((abs(rlat1-rlat2).lt.0.03).and. & (abs(rlon1-rlon2).lt.0.03)) then distance=0. else clat1=cos(real(rlat1,kind=dp)/dpr) slat1=sin(real(rlat1,kind=dp)/dpr) clat2=cos(real(rlat2,kind=dp)/dpr) slat2=sin(real(rlat2,kind=dp)/dpr) cdlon=cos(real((rlon1-rlon2),kind=dp)/dpr) crd=slat1*slat2+clat1*clat2*cdlon distance=real(rerth*acos(crd)/1000.0_dp) endif ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function distance real function distance2(rlat1,rlon1,rlat2,rlon2) !$$$ SUBPROGRAM DOCUMENTATION BLOCK ! ! SUBPROGRAM: GCDIST COMPUTE GREAT CIRCLE DISTANCE ! PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 ! ! ABSTRACT: THIS SUBPROGRAM COMPUTES GREAT CIRCLE DISTANCE ! BETWEEN TWO POINTS ON THE EARTH. COORDINATES ARE GIVEN IN RADIANS! ! ! PROGRAM HISTORY LOG: ! 96-04-10 IREDELL ! ! USAGE: ...GCDIST(RLAT1,RLON1,RLAT2,RLON2) ! ! INPUT ARGUMENT LIST: !rlat1 - REAL LATITUDE OF POINT 1 IN RADIANS !rlon1 - REAL LONGITUDE OF POINT 1 IN RADIANS !rlat2 - REAL LATITUDE OF POINT 2 IN RADIANS !rlon2 - REAL LONGITUDE OF POINT 2 IN RADIANS ! ! OUTPUT ARGUMENT LIST: !distance2 - REAL GREAT CIRCLE DISTANCE IN KM ! ! ATTRIBUTES: ! LANGUAGE: Fortran 90 ! !$$$ use par_mod, only: dp implicit none real :: rlat1,rlon1,rlat2,rlon2 real(kind=dp) :: clat1,clat2,slat1,slat2,cdlon,crd real(kind=dp),parameter :: rerth=6.3712e6_dp real(kind=dp),parameter :: pi=3.14159265358979_dp if ((abs(rlat1-rlat2).lt.0.0003).and. & (abs(rlon1-rlon2).lt.0.0003)) then distance2=0.0_dp else ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - clat1=cos(real(rlat1,kind=dp)) slat1=sin(real(rlat1,kind=dp)) clat2=cos(real(rlat2,kind=dp)) slat2=sin(real(rlat2,kind=dp)) cdlon=cos(real(rlon1-rlon2,kind=dp)) crd=slat1*slat2+clat1*clat2*cdlon distance2=real(rerth*acos(crd)/1000.0_dp) endif ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - end function distance2 subroutine openouttraj !***************************************************************************** ! * ! This routine opens the output file for the plume trajectory output * ! produced by the cluster analysis. * ! * ! Author: A. Stohl * ! * ! 27 January 2001 * ! * !***************************************************************************** ! * ! Variables: * ! * !***************************************************************************** use point_mod use par_mod use com_mod implicit none integer :: i real :: xp1,yp1,xp2,yp2 ! Open output file for trajectory output !*************************************** open(unitouttraj,file=path(2)(1:length(2))//'trajectories.txt', & form='formatted',err=998) if (ldirect.eq.1) then write(unitouttraj,'(i8,1x,i6,1x,a)') ibdate,ibtime, trim(flexversion) else write(unitouttraj,'(i8,1x,i6,1x,a)') iedate,ietime, trim(flexversion) endif write(unitouttraj,*) method,lsubgrid,lconvection write(unitouttraj,*) numpoint do i=1,numpoint xp1=xpoint1(i)*dx+xlon0 yp1=ypoint1(i)*dy+ylat0 xp2=xpoint2(i)*dx+xlon0 yp2=ypoint2(i)*dy+ylat0 write(unitouttraj,*) ireleasestart(i),ireleaseend(i), & xp1,yp1,xp2,yp2,zpoint1(i),zpoint2(i),kindz(i),npart(i) if (numpoint.le.1000) then write(unitouttraj,'(a)') compoint(i)(1:40) else write(unitouttraj,'(a)') compoint(1001)(1:40) endif end do return 998 write(*,*) ' #### FLEXPART MODEL ERROR! THE FILE #### ' write(*,*) ' #### trajectories.txt #### ' write(*,*) ' #### CANNOT BE OPENED. IF A FILE WITH THIS #### ' write(*,*) ' #### NAME ALREADY EXISTS, DELETE IT AND START #### ' write(*,*) ' #### THE PROGRAM AGAIN. #### ' error stop end subroutine openouttraj end module plume_mod