source: flexpart.git/src/partoutput_mpi.f90 @ 38b7917

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 38b7917 was 5f9d14a, checked in by Espen Sollum ATMOS <eso@…>, 9 years ago

Updated wet depo scheme

  • Property mode set to 100644
File size: 7.8 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(itime)
23  !                           i
24  !*****************************************************************************
25  !                                                                            *
26  !     Dump all particle positions                                            *
27  !                                                                            *
28  !     Author: A. Stohl                                                       *
29  !                                                                            *
30  !     12 March 1999                                                          *
31  !                                                                            *
32  !     12/2014 eso: Version for MPI                                           *
33  !                  processes sequentially access and append data to file     *
34  !                                                                            *
35  !*****************************************************************************
36  !                                                                            *
37  ! Variables:                                                                 *
38  !                                                                            *
39  !*****************************************************************************
40
41  use par_mod
42  use com_mod
43  use mpi_mod
44
45  implicit none
46
47  real(kind=dp) :: jul
48  integer :: itime,i,j,jjjjmmdd,ihmmss
49  integer :: ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
50  real :: xlon,ylat
51  real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
52  real :: topo,hm(2),hmixi,pv1(2),pvprof(2),pvi,qv1(2),qvprof(2),qvi
53  real :: tt1(2),ttprof(2),tti,rho1(2),rhoprof(2),rhoi
54  real :: tr(2),tri
55  character :: adate*8,atime*6
56  character(LEN=8) :: file_stat='OLD'
57
58  ! MPI root process creates the file, other processes append to it
59  if (lroot) file_stat='REPLACE'
60
61  ! Determine current calendar date, needed for the file name
62  !**********************************************************
63
64  jul=bdate+real(itime,kind=dp)/86400._dp
65  call caldate(jul,jjjjmmdd,ihmmss)
66  write(adate,'(i8.8)') jjjjmmdd
67  write(atime,'(i6.6)') ihmmss
68
69
70  ! Some variables needed for temporal interpolation
71  !*************************************************
72
73  dt1=real(itime-memtime(1))
74  dt2=real(memtime(2)-itime)
75  dtt=1./(dt1+dt2)
76
77  ! Open output file and write the output
78  !**************************************
79
80  if (ipout.eq.1) then
81    open(unitpartout,file=path(2)(1:length(2))//'partposit_'//adate// &
82         atime,form='unformatted',status=file_stat,position='append')
83  else
84    open(unitpartout,file=path(2)(1:length(2))//'partposit_end', &
85         form='unformatted',status=file_stat,position='append')
86  endif
87
88  ! Write current time to file
89  !***************************
90
91  if (lroot) write(unitpartout) itime ! MPI root process only
92  do i=1,numpart
93
94  ! Take only valid particles
95  !**************************
96
97    if (itra1(i).eq.itime) then
98      xlon=xlon0+xtra1(i)*dx
99      ylat=ylat0+ytra1(i)*dy
100
101  !*****************************************************************************
102  ! Interpolate several variables (PV, specific humidity, etc.) to particle position
103  !*****************************************************************************
104
105      ix=xtra1(i)
106      jy=ytra1(i)
107      ixp=ix+1
108      jyp=jy+1
109      ddx=xtra1(i)-real(ix)
110      ddy=ytra1(i)-real(jy)
111      rddx=1.-ddx
112      rddy=1.-ddy
113      p1=rddx*rddy
114      p2=ddx*rddy
115      p3=rddx*ddy
116      p4=ddx*ddy
117
118  ! Topography
119  !***********
120
121      topo=p1*oro(ix ,jy) &
122           + p2*oro(ixp,jy) &
123           + p3*oro(ix ,jyp) &
124           + p4*oro(ixp,jyp)
125
126  ! Potential vorticity, specific humidity, temperature, and density
127  !*****************************************************************
128
129      do il=2,nz
130        if (height(il).gt.ztra1(i)) then
131          indz=il-1
132          indzp=il
133          goto 6
134        endif
135      end do
1366     continue
137
138      dz1=ztra1(i)-height(indz)
139      dz2=height(indzp)-ztra1(i)
140      dz=1./(dz1+dz2)
141
142
143      do ind=indz,indzp
144        do m=1,2
145          indexh=memind(m)
146
147  ! Potential vorticity
148          pv1(m)=p1*pv(ix ,jy ,ind,indexh) &
149               +p2*pv(ixp,jy ,ind,indexh) &
150               +p3*pv(ix ,jyp,ind,indexh) &
151               +p4*pv(ixp,jyp,ind,indexh)
152  ! Specific humidity
153          qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
154               +p2*qv(ixp,jy ,ind,indexh) &
155               +p3*qv(ix ,jyp,ind,indexh) &
156               +p4*qv(ixp,jyp,ind,indexh)
157  ! Temperature
158          tt1(m)=p1*tt(ix ,jy ,ind,indexh) &
159               +p2*tt(ixp,jy ,ind,indexh) &
160               +p3*tt(ix ,jyp,ind,indexh) &
161               +p4*tt(ixp,jyp,ind,indexh)
162  ! Density
163          rho1(m)=p1*rho(ix ,jy ,ind,indexh) &
164               +p2*rho(ixp,jy ,ind,indexh) &
165               +p3*rho(ix ,jyp,ind,indexh) &
166               +p4*rho(ixp,jyp,ind,indexh)
167        end do
168        pvprof(ind-indz+1)=(pv1(1)*dt2+pv1(2)*dt1)*dtt
169        qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt
170        ttprof(ind-indz+1)=(tt1(1)*dt2+tt1(2)*dt1)*dtt
171        rhoprof(ind-indz+1)=(rho1(1)*dt2+rho1(2)*dt1)*dtt
172      end do
173      pvi=(dz1*pvprof(2)+dz2*pvprof(1))*dz
174      qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz
175      tti=(dz1*ttprof(2)+dz2*ttprof(1))*dz
176      rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
177
178  ! Tropopause and PBL height
179  !**************************
180
181      do m=1,2
182        indexh=memind(m)
183
184  ! Tropopause
185        tr(m)=p1*tropopause(ix ,jy ,1,indexh) &
186             + p2*tropopause(ixp,jy ,1,indexh) &
187             + p3*tropopause(ix ,jyp,1,indexh) &
188             + p4*tropopause(ixp,jyp,1,indexh)
189
190  ! PBL height
191        hm(m)=p1*hmix(ix ,jy ,1,indexh) &
192             + p2*hmix(ixp,jy ,1,indexh) &
193             + p3*hmix(ix ,jyp,1,indexh) &
194             + p4*hmix(ixp,jyp,1,indexh)
195      end do
196
197      hmixi=(hm(1)*dt2+hm(2)*dt1)*dtt
198      tri=(tr(1)*dt2+tr(2)*dt1)*dtt
199
200
201  ! Write the output
202  !*****************
203      if (mp_measure_time) call mpif_mtime('iotime',0)
204
205      write(unitpartout) npoint(i),xlon,ylat,ztra1(i), &
206           itramem(i),topo,pvi,qvi,rhoi,hmixi,tri,tti, &
207           (xmass1(i,j),j=1,nspec)
208
209      if (mp_measure_time) call mpif_mtime('iotime',1)
210    endif
211  end do
212
213  ! Only last MPI process writes EOF info
214  if (mp_partid.eq.mp_partgroup_np-1) then
215    write(unitpartout) -99999,-9999.9,-9999.9,-9999.9,-99999, &
216         -9999.9,-9999.9,-9999.9,-9999.9,-9999.9,-9999.9,-9999.9, &
217         (-9999.9,j=1,nspec)
218  end if
219
220
221  close(unitpartout)
222
223end subroutine partoutput
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG