source: branches/jerome/src_flexwrf_v3.1/partoutput.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago

sources for flexwrf v3.1
File size: 12.0 KB
Line 
1!***********************************************************************
2!* Copyright 2012,2013                                                *
3!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
4!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
6!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
7!*                                                                     *
8!* This file is part of FLEXPART WRF                                   *
9!*                                                                     *
10!* FLEXPART is free software: you can redistribute it and/or modify    *
11!* it under the terms of the GNU General Public License as published by*
12!* the Free Software Foundation, either version 3 of the License, or   *
13!* (at your option) any later version.                                 *
14!*                                                                     *
15!* FLEXPART is distributed in the hope that it will be useful,         *
16!* but WITHOUT ANY WARRANTY; without even the implied warranty of      *
17!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
18!* GNU General Public License for more details.                        *
19!*                                                                     *
20!* You should have received a copy of the GNU General Public License   *
21!* along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
22!***********************************************************************
23
24      subroutine partoutput(itime)
25!                             i
26!*******************************************************************************
27!                                                                              *
28!     Note:  This is the FLEXPART_WRF version of subroutine partoutput.        *
29!                                                                              *
30!     Dump all particle positions                                              *
31!                                                                              *
32!     Author: A. Stohl                                                         *
33!     12 March 1999                                                            *
34!                                                                              *
35!     Dec 2005, J. Fast - Output files can be either binary or ascii.          *
36!                   Topo,pv,qv,... at particle positions are calculated        *
37!                   using nested fields when (partoutput_use_nested .gt. 0)    *
38!                   Particle xy coords can be either lat-lon or grid-meters.   *
39!                   Changed names of "*lon0*" & "*lat0*" variables             *
40!                                                                              *
41!*******************************************************************************
42!                                                                              *
43! Variables:                                                                   *
44!                                                                              *
45!*******************************************************************************
46
47  use par_mod
48  use com_mod
49
50  implicit none
51!      include 'includepar'
52!      include 'includecom'
53
54  real(kind=dp) :: jul
55      integer :: itime,i,j,jjjjmmdd,ihmmss
56      integer :: ix,jy,ixp,jyp,indexh,m,il,ind,indz,indzp
57      integer :: numpart_out
58      real :: xlon,ylat,xtmp,ytmp
59      real :: dt1,dt2,dtt,ddx,ddy,rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
60      real :: topo,hm(2),hmixi,pv1(2),pvprof(2),pvi,qv1(2),qvprof(2),qvi
61      real :: tt1(2),ttprof(2),tti,rho1(2),rhoprof(2),rhoi
62      real :: tr(2),tri
63      character :: adate*8,atime*6
64
65      integer :: k, ngrid
66      real :: xtn, ytn
67
68
69! Determine current calendar date, needed for the file name
70!**********************************************************
71
72  jul=bdate+real(itime,kind=dp)/86400._dp    ! this is the current day
73
74      call caldate(jul,jjjjmmdd,ihmmss)
75      write(adate,'(i8.8)') jjjjmmdd
76      write(atime,'(i6.6)') ihmmss
77
78
79! Some variables needed for temporal interpolation
80!*************************************************
81
82      dt1=float(itime-memtime(1))
83      dt2=float(memtime(2)-itime)
84      dtt=1./(dt1+dt2)
85
86! Open output file and write the output
87!**************************************
88
89      if (ipout.eq.1) then
90        if (iouttype.eq.0) &
91        open(unitpartout,file=path(1)(1:length(1))//'partposit_'//adate// &
92        atime,form='unformatted')
93        if (iouttype.eq.1) &
94        open(unitpartout,file=path(1)(1:length(1))//'partposit_'//adate// &
95        atime,form='formatted')
96      else
97        if (iouttype.eq.0) &
98        open(unitpartout,file=path(1)(1:length(1))//'partposit_end', &
99        form='unformatted')
100        if (iouttype.eq.1) &
101        open(unitpartout,file=path(1)(1:length(1))//'partposit_end', &
102        form='formatted')
103      endif
104
105! Write current time to file
106!***************************
107
108      numpart_out = 0
109      do i=1,numpart
110        if (itra1(i).eq.itime) numpart_out = numpart_out + 1
111      enddo
112
113      if (iouttype.eq.0) write(unitpartout)   itime, &
114          numpart_out, outgrid_option
115      if (iouttype.eq.1) write(unitpartout,*) itime, &
116          numpart_out, outgrid_option
117
118      do i=1,numpart
119
120! Take only valid particles
121!**************************
122
123        if (itra1(i).eq.itime) then
124          xlon=xmet0+xtra1(i)*dx
125          ylat=ymet0+ytra1(i)*dy
126
127!*********************************************************************************
128! Interpolate several variables (PV, specific humidity, etc.) to particle position
129!*********************************************************************************
130
131! If partoutput_use_nested=0, set ngrid=0, and use the outermost grid
132!    for calculating topo, pv, qv, ... at the particle position
133! Otherwise, determine the nest we are in
134          ngrid=0
135          if (partoutput_use_nested .gt. 0) then
136          do k=numbnests,1,-1
137             if ((xtra1(i).gt.xln(k)).and. &
138                 (xtra1(i).lt.xrn(k)).and. &
139                 (ytra1(i).gt.yln(k)).and. &
140                 (ytra1(i).lt.yrn(k))) then
141                    ngrid=k
142                    goto 26
143             endif
144          enddo
14526        continue
146          endif
147
148          if (ngrid .le. 0) then
149             ix=int(xtra1(i))
150             jy=int(ytra1(i))
151             ddy=ytra1(i)-float(jy)
152             ddx=xtra1(i)-float(ix)
153          else
154             xtn=(xtra1(i)-xln(ngrid))*xresoln(ngrid)
155             ytn=(ytra1(i)-yln(ngrid))*yresoln(ngrid)
156             ix=int(xtn)
157             jy=int(ytn)
158             ddy=ytn-float(jy)
159             ddx=xtn-float(ix)
160          endif
161          ixp=ix+1
162          jyp=jy+1
163          rddx=1.-ddx
164          rddy=1.-ddy
165          p1=rddx*rddy
166          p2=ddx*rddy
167          p3=rddx*ddy
168          p4=ddx*ddy
169
170! Topography
171!***********
172          if (ngrid .le. 0) then
173      topo=p1*oro(ix ,jy) &
174           + p2*oro(ixp,jy) &
175           + p3*oro(ix ,jyp) &
176           + p4*oro(ixp,jyp)
177          else
178             topo=p1*oron(ix ,jy ,ngrid) &
179                + p2*oron(ixp,jy ,ngrid) &
180                + p3*oron(ix ,jyp,ngrid) &
181                + p4*oron(ixp,jyp,ngrid)
182          endif
183
184! Potential vorticity, specific humidity, temperature, and density
185!*****************************************************************
186
187          do il=2,nz
188            if (height(il).gt.ztra1(i)) then
189              indz=il-1
190              indzp=il
191              goto 56
192            endif
193          enddo
19456        continue
195
196          dz1=ztra1(i)-height(indz)
197          dz2=height(indzp)-ztra1(i)
198          dz=1./(dz1+dz2)
199
200
201          do ind=indz,indzp
202            do m=1,2
203              indexh=memind(m)
204
205              if (ngrid .le. 0) then
206  ! Potential vorticity
207          pv1(m)=p1*pv(ix ,jy ,ind,indexh) &
208               +p2*pv(ixp,jy ,ind,indexh) &
209               +p3*pv(ix ,jyp,ind,indexh) &
210               +p4*pv(ixp,jyp,ind,indexh)
211  ! Specific humidity
212          qv1(m)=p1*qv(ix ,jy ,ind,indexh) &
213               +p2*qv(ixp,jy ,ind,indexh) &
214               +p3*qv(ix ,jyp,ind,indexh) &
215               +p4*qv(ixp,jyp,ind,indexh)
216  ! Temperature
217          tt1(m)=p1*tt(ix ,jy ,ind,indexh) &
218               +p2*tt(ixp,jy ,ind,indexh) &
219               +p3*tt(ix ,jyp,ind,indexh) &
220               +p4*tt(ixp,jyp,ind,indexh)
221  ! Density
222          rho1(m)=p1*rho(ix ,jy ,ind,indexh) &
223               +p2*rho(ixp,jy ,ind,indexh) &
224               +p3*rho(ix ,jyp,ind,indexh) &
225               +p4*rho(ixp,jyp,ind,indexh)
226              else
227              pv1(m)=p1*pvn(ix ,jy ,ind,indexh,ngrid) &
228                    +p2*pvn(ixp,jy ,ind,indexh,ngrid) &
229                    +p3*pvn(ix ,jyp,ind,indexh,ngrid) &
230                    +p4*pvn(ixp,jyp,ind,indexh,ngrid)
231              qv1(m)=p1*qvn(ix ,jy ,ind,indexh,ngrid) &
232                    +p2*qvn(ixp,jy ,ind,indexh,ngrid) &
233                    +p3*qvn(ix ,jyp,ind,indexh,ngrid) &
234                    +p4*qvn(ixp,jyp,ind,indexh,ngrid)
235              tt1(m)=p1*ttn(ix ,jy ,ind,indexh,ngrid) &
236                    +p2*ttn(ixp,jy ,ind,indexh,ngrid) &
237                    +p3*ttn(ix ,jyp,ind,indexh,ngrid) &
238                    +p4*ttn(ixp,jyp,ind,indexh,ngrid)
239              rho1(m)=p1*rhon(ix ,jy ,ind,indexh,ngrid) &
240                     +p2*rhon(ixp,jy ,ind,indexh,ngrid) &
241                     +p3*rhon(ix ,jyp,ind,indexh,ngrid) &
242                     +p4*rhon(ixp,jyp,ind,indexh,ngrid)
243
244              endif
245
246         enddo
247            pvprof(ind-indz+1)=(pv1(1)*dt2+pv1(2)*dt1)*dtt
248            qvprof(ind-indz+1)=(qv1(1)*dt2+qv1(2)*dt1)*dtt
249            ttprof(ind-indz+1)=(tt1(1)*dt2+tt1(2)*dt1)*dtt
250            rhoprof(ind-indz+1)=(rho1(1)*dt2+rho1(2)*dt1)*dtt
251         enddo
252          pvi=(dz1*pvprof(2)+dz2*pvprof(1))*dz
253          qvi=(dz1*qvprof(2)+dz2*qvprof(1))*dz
254          tti=(dz1*ttprof(2)+dz2*ttprof(1))*dz
255          rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
256
257! Tropopause and PBL height
258!**************************
259
260          do m=1,2
261            indexh=memind(m)
262
263            if (ngrid .le. 0) then
264! Tropopause
265            tr(m)=p1*tropopause(ix ,jy ,1,indexh) &
266                + p2*tropopause(ixp,jy ,1,indexh) &
267                + p3*tropopause(ix ,jyp,1,indexh) &
268                + p4*tropopause(ixp,jyp,1,indexh)
269! PBL height
270            hm(m)=p1*hmix(ix ,jy ,1,indexh) &
271                + p2*hmix(ixp,jy ,1,indexh) &
272                + p3*hmix(ix ,jyp,1,indexh) &
273                + p4*hmix(ixp,jyp,1,indexh)
274 
275            else
276            tr(m)=p1*tropopausen(ix ,jy ,1,indexh,ngrid) &
277                + p2*tropopausen(ixp,jy ,1,indexh,ngrid) &
278                + p3*tropopausen(ix ,jyp,1,indexh,ngrid) &
279                + p4*tropopausen(ixp,jyp,1,indexh,ngrid)
280            hm(m)=p1*hmixn(ix ,jy ,1,indexh,ngrid) &
281                + p2*hmixn(ixp,jy ,1,indexh,ngrid) &
282                + p3*hmixn(ix ,jyp,1,indexh,ngrid) &
283                + p4*hmixn(ixp,jyp,1,indexh,ngrid)
284
285            endif
286
287        enddo
288
289          hmixi=(hm(1)*dt2+hm(2)*dt1)*dtt
290          tri=(tr(1)*dt2+tr(2)*dt1)*dtt
291
292
293! Write the output
294!*****************
295
296          if (outgrid_option .eq. 1) then
297             xtmp = xlon
298             ytmp = ylat
299             call xymeter_to_ll_wrf( xtmp, ytmp, xlon, ylat )
300          endif
301          if(iouttype.eq.0)  &   
302          write(unitpartout) npoint(i),xlon,ylat,ztra1(i), &
303          itramem(i),topo,pvi,qvi,rhoi,hmixi,tri,tti, &
304          (xmass1(i,j),j=1,nspec)
305          if(iouttype.eq.1)    &
306          write(unitpartout,101) npoint(i),itramem(i),xlon,ylat, &
307          ztra1(i),topo,pvi,qvi,rhoi,hmixi,tri,tti, &
308          (xmass1(i,j),j=1,nspec)
309        endif
310
311       enddo
312
313      if(iouttype.eq.0)  &
314      write(unitpartout) -99999,-9999.9,-9999.9,-9999.9,-99999, &
315      -9999.9,-9999.9,-9999.9,-9999.9,-9999.9,-9999.9,-9999.9, &
316      (-9999.9,j=1,nspec)
317      if(iouttype.eq.1)  &
318      write(unitpartout,101) -99999,-99999,-9999.9,-9999.9,-9999.9, &
319      -9999.9,-9999.9,-9999.9,-9999.9,-9999.9,-9999.9,-9999.9, &
320      (-9999.9,j=1,nspec)
321
322101   format( 2i10, 1p, 21e14.6 )
323
324
325      close(unitpartout)
326
327end subroutine partoutput
328
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG