source: flexpart.git/src/partoutput_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: 6.8 KB
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4subroutine partoutput_average(itime)
5  !                             i
6  !*****************************************************************************
7  !                                                                            *
8  !     Dump all particle positions                                            *
9  !                                                                            *
10  !     Author: A. Stohl                                                       *
11  !                                                                            *
12  !     12 March 1999                                                          *
13  !                                                                            *
14  !*****************************************************************************
15  !                                                                            *
16  ! Variables:                                                                 *
17  !                                                                            *
18  !*****************************************************************************
19
20  use par_mod
21  use com_mod
22
23  implicit none
24
25  real(kind=dp) :: jul
26  integer :: itime,i,j,jjjjmmdd,ihmmss
27  integer :: ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
28  real :: xlon,ylat
29  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
30  real :: topo,hm(2),hmixi,pv1(2),pvprof(2),pvi,qv1(2),qvprof(2),qvi
31  real :: tt1(2),ttprof(2),tti,rho1(2),rhoprof(2),rhoi
32  real :: tr(2),tri,zlim
33  character :: adate*8,atime*6
34
35  integer(kind=2) :: ishort_xlon(maxpart),ishort_ylat(maxpart),ishort_z(maxpart)
36  integer(kind=2) :: ishort_topo(maxpart),ishort_tro(maxpart),ishort_hmix(maxpart)
37  integer(kind=2) :: ishort_pv(maxpart),ishort_rho(maxpart),ishort_qv(maxpart)
38  integer(kind=2) :: ishort_tt(maxpart),ishort_uu(maxpart),ishort_vv(maxpart)
39  integer(kind=2) :: ishort_energy(maxpart)
40
41
42  ! Determine current calendar date, needed for the file name
43  !**********************************************************
44
45  jul=bdate+real(itime,kind=dp)/86400._dp
46  call caldate(jul,jjjjmmdd,ihmmss)
47  write(adate,'(i8.8)') jjjjmmdd
48  write(atime,'(i6.6)') ihmmss
49
50
51  ! Some variables needed for temporal interpolation
52  !*************************************************
53
54  dt1=real(itime-memtime(1))
55  dt2=real(memtime(2)-itime)
56  dtt=1./(dt1+dt2)
57
58  ! Open output file and write the output
59  !**************************************
60
61    open(unitpartout_average,file=path(2)(1:length(2))//'partposit_average_'//adate// &
62         atime,form='unformatted',access='direct',status='replace',recl=24)
63
64
65  ! Write current time to file
66  !***************************
67
68!  write(unitpartout_average) itime,numpart
69  do i=1,numpart
70
71  ! Take only valid particles
72  !**************************
73
74    if (itra1(i).eq.itime) then
75      part_av_topo(i)=part_av_topo(i)/float(npart_av(i))
76      part_av_z(i)=part_av_z(i)/float(npart_av(i))
77      part_av_pv(i)=part_av_pv(i)/float(npart_av(i))
78      part_av_qv(i)=part_av_qv(i)/float(npart_av(i))
79      part_av_tt(i)=part_av_tt(i)/float(npart_av(i))
80      part_av_uu(i)=part_av_uu(i)/float(npart_av(i))
81      part_av_vv(i)=part_av_vv(i)/float(npart_av(i))
82      part_av_rho(i)=part_av_rho(i)/float(npart_av(i))
83      part_av_tro(i)=part_av_tro(i)/float(npart_av(i))
84      part_av_hmix(i)=part_av_hmix(i)/float(npart_av(i))
85      part_av_energy(i)=part_av_energy(i)/float(npart_av(i))
86
87      part_av_cartx(i)=part_av_cartx(i)/float(npart_av(i))
88      part_av_carty(i)=part_av_carty(i)/float(npart_av(i))
89      part_av_cartz(i)=part_av_cartz(i)/float(npart_av(i))
90
91! Project Cartesian coordinates back onto Earth's surface
92! #######################################################
93      xlon=atan2(part_av_cartx(i),-1.*part_av_carty(i))
94      ylat=atan2(part_av_cartz(i),sqrt(part_av_cartx(i)*part_av_cartx(i)+ &
95           part_av_carty(i)*part_av_carty(i)))
96      xlon=xlon/pi180
97      ylat=ylat/pi180
98
99
100  ! Convert all data to integer*2 variables (from -32768 to 32767) for compressed output
101  !*************************************************************************************
102
103      if (xlon.gt.180.) xlon=xlon-360.
104      if (xlon.lt.-180.) xlon=xlon+360.
105      ishort_xlon(i)=nint(xlon*180.)
106      ishort_ylat(i)=nint(ylat*360.)
107
108      zlim=(part_av_z(i)*2.-32000.)
109      zlim=min(zlim,32766.)
110      zlim=max(zlim,-32766.)
111      ishort_z(i)=nint(zlim)
112
113      zlim=(part_av_topo(i)*2.-32000.)
114      zlim=min(zlim,32766.)
115      zlim=max(zlim,-32766.)
116      ishort_topo(i)=nint(zlim)
117
118      zlim=(part_av_tro(i)*2.-32000.)
119      zlim=min(zlim,32766.)
120      zlim=max(zlim,-32766.)
121      ishort_tro(i)=nint(zlim)
122
123      zlim=(part_av_hmix(i)*2.-32000.)
124      zlim=min(zlim,32766.)
125      zlim=max(zlim,-32766.)
126      ishort_hmix(i)=nint(zlim)
127
128      zlim=(part_av_rho(i)*20000.-32000.)
129      zlim=min(zlim,32766.)
130      zlim=max(zlim,-32766.)
131      ishort_rho(i)=nint(zlim)
132
133      zlim=(part_av_qv(i)*1000000.-32000.)
134      zlim=min(zlim,32766.)
135      zlim=max(zlim,-32766.)
136      ishort_qv(i)=nint(zlim)
137
138      zlim=(part_av_pv(i)*100.)
139      zlim=min(zlim,32766.)
140      zlim=max(zlim,-32766.)
141      ishort_pv(i)=nint(zlim)
142
143      zlim=((part_av_tt(i)-273.15))*300.
144      zlim=min(zlim,32766.)
145      zlim=max(zlim,-32766.)
146      ishort_tt(i)=nint(zlim)
147
148      zlim=(part_av_uu(i)*200.)
149      zlim=min(zlim,32766.)
150      zlim=max(zlim,-32766.)
151      ishort_uu(i)=nint(zlim)
152
153      zlim=(part_av_vv(i)*200.)
154      zlim=min(zlim,32766.)
155      zlim=max(zlim,-32766.)
156      ishort_vv(i)=nint(zlim)
157
158      zlim=(part_av_energy(i)-300000.)/30.
159      zlim=min(zlim,32766.)
160      zlim=max(zlim,-32766.)
161      ishort_energy(i)=nint(zlim)
162
163  ! Turn on for readable test output
164  !*********************************
165
166!        write(119,*) itime,i,xlon,ylat,part_av_z(i),part_av_topo(i),part_av_tro(i), &
167!        part_av_hmix(i),part_av_rho(i),part_av_qv(i),part_av_pv(i),part_av_tt(i), &
168!        ishort_uu(i),ishort_vv(i)
169    endif
170
171! Re-initialize averages, and set number of steps to zero
172    npart_av(i)=0
173    part_av_topo(i)=0.
174    part_av_z(i)=0.
175    part_av_cartx(i)=0.
176    part_av_carty(i)=0.
177    part_av_cartz(i)=0.
178    part_av_pv(i)=0.
179    part_av_qv(i)=0.
180    part_av_tt(i)=0.
181    part_av_uu(i)=0.
182    part_av_vv(i)=0.
183    part_av_rho(i)=0.
184    part_av_tro(i)=0.
185    part_av_hmix(i)=0.
186    part_av_energy(i)=0.
187   
188  end do
189
190  ! Write the output
191  !*****************
192  do i=1,numpart
193    if (itra1(i).eq.itime) then
194      write(unitpartout_average,rec=i) ishort_xlon(i),ishort_ylat(i),ishort_z(i), &
195           ishort_topo(i),ishort_tro(i),ishort_hmix(i),ishort_rho(i),ishort_qv(i),ishort_pv(i), &
196           ishort_tt(i),ishort_uu(i),ishort_vv(i)  ! ,ishort_energy(i)
197    endif
198  enddo
199
200
201  close(unitpartout_average)
202
203end subroutine partoutput_average
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG