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