source: flexpart.git/src/partpos_average.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: 5.2 KB
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4
5subroutine partpos_average(itime,j)
6
7
8!**********************************************************************
9! This subroutine averages particle quantities, to be used for particle
10! dump (in partoutput.f90). Averaging is done over output interval.
11!**********************************************************************
12
13  use par_mod
14  use com_mod
15
16  implicit none
17
18  integer :: itime,j,ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
19  real :: xlon,ylat,x,y,z
20  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
21  real :: topo,hm(2),hmixi,pv1(2),pvprof(2),pvi,qv1(2),qvprof(2),qvi
22  real :: tt1(2),ttprof(2),tti,rho1(2),rhoprof(2),rhoi
23  real :: uu1(2),uuprof(2),uui,vv1(2),vvprof(2),vvi
24  real :: tr(2),tri,energy
25
26
27
28 ! Some variables needed for temporal interpolation
29  !*************************************************
30
31  dt1=real(itime-memtime(1))
32  dt2=real(memtime(2)-itime)
33  dtt=1./(dt1+dt2)
34
35  xlon=xlon0+xtra1(j)*dx
36  ylat=ylat0+ytra1(j)*dy
37
38  !*****************************************************************************
39  ! Interpolate several variables (PV, specific humidity, etc.) to particle position
40  !*****************************************************************************
41
42  ix=xtra1(j)
43  jy=ytra1(j)
44  ixp=ix+1
45  jyp=jy+1
46  ddx=xtra1(j)-real(ix)
47  ddy=ytra1(j)-real(jy)
48  rddx=1.-ddx
49  rddy=1.-ddy
50  p1=rddx*rddy
51  p2=ddx*rddy
52  p3=rddx*ddy
53  p4=ddx*ddy
54
55! eso: Temporary fix for particle exactly at north pole
56  if (jyp >= nymax) then
57  !  write(*,*) 'WARNING: conccalc.f90 jyp >= nymax'
58    jyp=jyp-1
59  end if
60
61  ! Topography
62  !***********
63
64  topo=p1*oro(ix,jy)+p2*oro(ixp,jy)+p3*oro(ix,jyp)+p4*oro(ixp,jyp)
65
66 ! Potential vorticity, specific humidity, temperature, and density
67  !*****************************************************************
68
69  do il=2,nz
70    if (height(il).gt.ztra1(j)) then
71      indz=il-1
72      indzp=il
73      goto 6
74    endif
75  end do
766 continue
77
78  dz1=ztra1(j)-height(indz)
79  dz2=height(indzp)-ztra1(j)
80  dz=1./(dz1+dz2)
81
82
83  do ind=indz,indzp
84    do m=1,2
85      indexh=memind(m)
86
87  ! Potential vorticity
88      pv1(m)=p1*pv(ix ,jy ,ind,indexh) &
89            +p2*pv(ixp,jy ,ind,indexh) &
90            +p3*pv(ix ,jyp,ind,indexh) &
91            +p4*pv(ixp,jyp,ind,indexh)
92  ! Specific humidity
93      qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
94            +p2*qv(ixp,jy ,ind,indexh) &
95            +p3*qv(ix ,jyp,ind,indexh) &
96            +p4*qv(ixp,jyp,ind,indexh)
97  ! Temperature
98      tt1(m)=p1*tt(ix ,jy ,ind,indexh) &
99            +p2*tt(ixp,jy ,ind,indexh) &
100            +p3*tt(ix ,jyp,ind,indexh) &
101            +p4*tt(ixp,jyp,ind,indexh)
102  ! U wind
103      uu1(m)=p1*uu(ix ,jy ,ind,indexh) &
104            +p2*uu(ixp,jy ,ind,indexh) &
105            +p3*uu(ix ,jyp,ind,indexh) &
106            +p4*uu(ixp,jyp,ind,indexh)
107  ! V wind
108      vv1(m)=p1*vv(ix ,jy ,ind,indexh) &
109            +p2*vv(ixp,jy ,ind,indexh) &
110            +p3*vv(ix ,jyp,ind,indexh) &
111            +p4*vv(ixp,jyp,ind,indexh)
112  ! Density
113      rho1(m)=p1*rho(ix ,jy ,ind,indexh) &
114             +p2*rho(ixp,jy ,ind,indexh) &
115             +p3*rho(ix ,jyp,ind,indexh) &
116             +p4*rho(ixp,jyp,ind,indexh)
117    end do
118    pvprof(ind-indz+1)=(pv1(1)*dt2+pv1(2)*dt1)*dtt
119    qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt
120    ttprof(ind-indz+1)=(tt1(1)*dt2+tt1(2)*dt1)*dtt
121    uuprof(ind-indz+1)=(uu1(1)*dt2+uu1(2)*dt1)*dtt
122    vvprof(ind-indz+1)=(vv1(1)*dt2+vv1(2)*dt1)*dtt
123    rhoprof(ind-indz+1)=(rho1(1)*dt2+rho1(2)*dt1)*dtt
124  end do
125  pvi=(dz1*pvprof(2)+dz2*pvprof(1))*dz
126  qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz
127  tti=(dz1*ttprof(2)+dz2*ttprof(1))*dz
128  uui=(dz1*uuprof(2)+dz2*uuprof(1))*dz
129  vvi=(dz1*vvprof(2)+dz2*vvprof(1))*dz
130  rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
131
132  ! Tropopause and PBL height
133  !**************************
134
135  do m=1,2
136    indexh=memind(m)
137
138  ! Tropopause
139    tr(m)=p1*tropopause(ix ,jy ,1,indexh) &
140        + p2*tropopause(ixp,jy ,1,indexh) &
141        + p3*tropopause(ix ,jyp,1,indexh) &
142        + p4*tropopause(ixp,jyp,1,indexh)
143
144  ! PBL height
145    hm(m)=p1*hmix(ix ,jy ,1,indexh) &
146        + p2*hmix(ixp,jy ,1,indexh) &
147        + p3*hmix(ix ,jyp,1,indexh) &
148        + p4*hmix(ixp,jyp,1,indexh)
149  end do
150
151  hmixi=(hm(1)*dt2+hm(2)*dt1)*dtt
152  tri=(tr(1)*dt2+tr(2)*dt1)*dtt
153
154
155  energy=tti*cpa+(ztra1(j)+topo)*9.81+qvi*2501000.+(uui**2+vvi**2)/2.
156
157  ! Add new values to sum and increase counter by one
158  !**************************************************
159
160  npart_av(j)=npart_av(j)+1
161
162  ! Calculate Cartesian 3D coordinates suitable for averaging
163  !**********************************************************
164
165  xlon=xlon*pi180
166  ylat=ylat*pi180
167  x = cos(ylat)*sin(xlon)
168  y = -1.*cos(ylat)*cos(xlon)
169  z = sin(ylat)
170
171  part_av_cartx(j)=part_av_cartx(j)+x
172  part_av_carty(j)=part_av_carty(j)+y
173  part_av_cartz(j)=part_av_cartz(j)+z
174  part_av_z(j)=part_av_z(j)+ztra1(j)
175  part_av_topo(j)=part_av_topo(j)+topo
176  part_av_pv(j)=part_av_pv(j)+pvi
177  part_av_qv(j)=part_av_qv(j)+qvi
178  part_av_tt(j)=part_av_tt(j)+tti
179  part_av_uu(j)=part_av_uu(j)+uui
180  part_av_vv(j)=part_av_vv(j)+vvi
181  part_av_rho(j)=part_av_rho(j)+rhoi
182  part_av_tro(j)=part_av_tro(j)+tri
183  part_av_hmix(j)=part_av_hmix(j)+hmixi
184  part_av_energy(j)=part_av_energy(j)+energy
185
186
187return
188end subroutine partpos_average
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG