source: flexpart.git/src/partoutput_average.f90 @ 332fbbd

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since 332fbbd was 332fbbd, checked in by Ignacio Pisso <ip@…>, 4 years ago

Revert "move license from headers to a different file"
Adding the full GPL license it was apparent that the FSF recommends to
have the statements as file headers

This reverts commit 3481cc180a6b19f71d9d2d04240a91d23f0a2800.
  • Property mode set to 100644
File size: 8.1 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine partoutput_average(itime)
23  !                             i
24  !*****************************************************************************
25  !                                                                            *
26  !     Dump all particle positions                                            *
27  !                                                                            *
28  !     Author: A. Stohl                                                       *
29  !                                                                            *
30  !     12 March 1999                                                          *
31  !                                                                            *
32  !*****************************************************************************
33  !                                                                            *
34  ! Variables:                                                                 *
35  !                                                                            *
36  !*****************************************************************************
37
38  use par_mod
39  use com_mod
40
41  implicit none
42
43  real(kind=dp) :: jul
44  integer :: itime,i,j,jjjjmmdd,ihmmss
45  integer :: ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
46  real :: xlon,ylat
47  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
48  real :: topo,hm(2),hmixi,pv1(2),pvprof(2),pvi,qv1(2),qvprof(2),qvi
49  real :: tt1(2),ttprof(2),tti,rho1(2),rhoprof(2),rhoi
50  real :: tr(2),tri,zlim
51  character :: adate*8,atime*6
52
53  integer(kind=2) :: ishort_xlon(maxpart),ishort_ylat(maxpart),ishort_z(maxpart)
54  integer(kind=2) :: ishort_topo(maxpart),ishort_tro(maxpart),ishort_hmix(maxpart)
55  integer(kind=2) :: ishort_pv(maxpart),ishort_rho(maxpart),ishort_qv(maxpart)
56  integer(kind=2) :: ishort_tt(maxpart),ishort_uu(maxpart),ishort_vv(maxpart)
57  integer(kind=2) :: ishort_energy(maxpart)
58
59
60  ! Determine current calendar date, needed for the file name
61  !**********************************************************
62
63  jul=bdate+real(itime,kind=dp)/86400._dp
64  call caldate(jul,jjjjmmdd,ihmmss)
65  write(adate,'(i8.8)') jjjjmmdd
66  write(atime,'(i6.6)') ihmmss
67
68
69  ! Some variables needed for temporal interpolation
70  !*************************************************
71
72  dt1=real(itime-memtime(1))
73  dt2=real(memtime(2)-itime)
74  dtt=1./(dt1+dt2)
75
76  ! Open output file and write the output
77  !**************************************
78
79    open(unitpartout_average,file=path(2)(1:length(2))//'partposit_average_'//adate// &
80         atime,form='unformatted',access='direct',status='replace',recl=24)
81
82
83  ! Write current time to file
84  !***************************
85
86!  write(unitpartout_average) itime,numpart
87  do i=1,numpart
88
89  ! Take only valid particles
90  !**************************
91
92    if (itra1(i).eq.itime) then
93      part_av_topo(i)=part_av_topo(i)/float(npart_av(i))
94      part_av_z(i)=part_av_z(i)/float(npart_av(i))
95      part_av_pv(i)=part_av_pv(i)/float(npart_av(i))
96      part_av_qv(i)=part_av_qv(i)/float(npart_av(i))
97      part_av_tt(i)=part_av_tt(i)/float(npart_av(i))
98      part_av_uu(i)=part_av_uu(i)/float(npart_av(i))
99      part_av_vv(i)=part_av_vv(i)/float(npart_av(i))
100      part_av_rho(i)=part_av_rho(i)/float(npart_av(i))
101      part_av_tro(i)=part_av_tro(i)/float(npart_av(i))
102      part_av_hmix(i)=part_av_hmix(i)/float(npart_av(i))
103      part_av_energy(i)=part_av_energy(i)/float(npart_av(i))
104
105      part_av_cartx(i)=part_av_cartx(i)/float(npart_av(i))
106      part_av_carty(i)=part_av_carty(i)/float(npart_av(i))
107      part_av_cartz(i)=part_av_cartz(i)/float(npart_av(i))
108
109! Project Cartesian coordinates back onto Earth's surface
110! #######################################################
111      xlon=atan2(part_av_cartx(i),-1.*part_av_carty(i))
112      ylat=atan2(part_av_cartz(i),sqrt(part_av_cartx(i)*part_av_cartx(i)+ &
113           part_av_carty(i)*part_av_carty(i)))
114      xlon=xlon/pi180
115      ylat=ylat/pi180
116
117
118  ! Convert all data to integer*2 variables (from -32768 to 32767) for compressed output
119  !*************************************************************************************
120
121      if (xlon.gt.180.) xlon=xlon-360.
122      if (xlon.lt.-180.) xlon=xlon+360.
123      ishort_xlon(i)=nint(xlon*180.)
124      ishort_ylat(i)=nint(ylat*360.)
125
126      zlim=(part_av_z(i)*2.-32000.)
127      zlim=min(zlim,32766.)
128      zlim=max(zlim,-32766.)
129      ishort_z(i)=nint(zlim)
130
131      zlim=(part_av_topo(i)*2.-32000.)
132      zlim=min(zlim,32766.)
133      zlim=max(zlim,-32766.)
134      ishort_topo(i)=nint(zlim)
135
136      zlim=(part_av_tro(i)*2.-32000.)
137      zlim=min(zlim,32766.)
138      zlim=max(zlim,-32766.)
139      ishort_tro(i)=nint(zlim)
140
141      zlim=(part_av_hmix(i)*2.-32000.)
142      zlim=min(zlim,32766.)
143      zlim=max(zlim,-32766.)
144      ishort_hmix(i)=nint(zlim)
145
146      zlim=(part_av_rho(i)*20000.-32000.)
147      zlim=min(zlim,32766.)
148      zlim=max(zlim,-32766.)
149      ishort_rho(i)=nint(zlim)
150
151      zlim=(part_av_qv(i)*1000000.-32000.)
152      zlim=min(zlim,32766.)
153      zlim=max(zlim,-32766.)
154      ishort_qv(i)=nint(zlim)
155
156      zlim=(part_av_pv(i)*100.)
157      zlim=min(zlim,32766.)
158      zlim=max(zlim,-32766.)
159      ishort_pv(i)=nint(zlim)
160
161      zlim=((part_av_tt(i)-273.15))*300.
162      zlim=min(zlim,32766.)
163      zlim=max(zlim,-32766.)
164      ishort_tt(i)=nint(zlim)
165
166      zlim=(part_av_uu(i)*200.)
167      zlim=min(zlim,32766.)
168      zlim=max(zlim,-32766.)
169      ishort_uu(i)=nint(zlim)
170
171      zlim=(part_av_vv(i)*200.)
172      zlim=min(zlim,32766.)
173      zlim=max(zlim,-32766.)
174      ishort_vv(i)=nint(zlim)
175
176      zlim=(part_av_energy(i)-300000.)/30.
177      zlim=min(zlim,32766.)
178      zlim=max(zlim,-32766.)
179      ishort_energy(i)=nint(zlim)
180
181  ! Turn on for readable test output
182  !*********************************
183
184!        write(119,*) itime,i,xlon,ylat,part_av_z(i),part_av_topo(i),part_av_tro(i), &
185!        part_av_hmix(i),part_av_rho(i),part_av_qv(i),part_av_pv(i),part_av_tt(i), &
186!        ishort_uu(i),ishort_vv(i)
187    endif
188
189! Re-initialize averages, and set number of steps to zero
190    npart_av(i)=0
191    part_av_topo(i)=0.
192    part_av_z(i)=0.
193    part_av_cartx(i)=0.
194    part_av_carty(i)=0.
195    part_av_cartz(i)=0.
196    part_av_pv(i)=0.
197    part_av_qv(i)=0.
198    part_av_tt(i)=0.
199    part_av_uu(i)=0.
200    part_av_vv(i)=0.
201    part_av_rho(i)=0.
202    part_av_tro(i)=0.
203    part_av_hmix(i)=0.
204    part_av_energy(i)=0.
205   
206  end do
207
208  ! Write the output
209  !*****************
210  do i=1,numpart
211    if (itra1(i).eq.itime) then
212      write(unitpartout_average,rec=i) ishort_xlon(i),ishort_ylat(i),ishort_z(i), &
213           ishort_topo(i),ishort_tro(i),ishort_hmix(i),ishort_rho(i),ishort_qv(i),ishort_pv(i), &
214           ishort_tt(i),ishort_uu(i),ishort_vv(i)  ! ,ishort_energy(i)
215    endif
216  enddo
217
218
219  close(unitpartout_average)
220
221end subroutine partoutput_average
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG