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

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

sources for flexwrf v3.1

File size: 13.2 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      subroutine calcpv_nests(l,n,uuhn,vvhn,pvhn)
24!                             i i  i    i    o
25!*******************************************************************************
26!                                                                              *
27!     Note:  This is the FLEXPART_WRF version of subroutine calcpv_nests.      *
28!            The computational grid is the WRF x-y grid rather than lat-lon.   *
29!                                                                              *
30!  Calculation of potential vorticity on 3-d nested grid                       *
31!                                                                              *
32!     Author: P. James                                                         *
33!     22 February 2000                                                         *
34!                                                                              *
35!    11 Nov 2005, R. Easter - changes associated with WRF horizontal grid.     *
36!                             For pressure use pph instead of (akz + bkz*ps)   *
37!    *** Note -- see ??? comments below regarding the pvh calculation.         *
38!                                                                              *
39!*******************************************************************************
40!                                                                              *
41! Variables:                                                                   *
42! n                  temporal index for meteorological fields (1 to 2)         *
43! l                  index of current nest                                     *
44!                                                                              *
45! Constants:                                                                   *
46!                                                                              *
47!*******************************************************************************
48
49!     include 'includepar'
50!     include 'includecom'
51!
52!     integer n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
53!     integer jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
54!     integer nlck,l
55!     real vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
56!     real theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
57!     real altit(nuvzmax),ppml(nuvzmax)
58!     real thup,thdn,eps,p0
59!     parameter(eps=1.e-5,p0=101325)
60!     real thhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
61!     real uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
62!     real vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
63!     real pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
64
65  use par_mod
66  use com_mod
67
68  implicit none
69
70  integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
71  integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
72  integer :: nlck,l
73  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
74  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt,ppmk
75  real :: ppml(nuvzmax)
76  real :: thup,thdn
77  real,parameter :: eps=1.e-5,p0=101325
78  real(kind=4) :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
79  real(kind=4) :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
80  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
81
82  real :: dumlon,dumlat
83  real :: thhn(0:nxmax-1,0:nymax-1,nuvzmax,maxnests)
84  real :: altit(nuvzmax)
85
86! Set number of levels to check for adjacent theta
87      nlck=nuvz/3
88! FLEXPART_WRF -- altit is never used, so don't calculate it
89!      do 5 k=1,nuvz
90!        altit(k)=akz(k)/p0+bkz(k)
91!5     continue
92! *** Precalculate all theta levels for efficiency
93        do jy=0,nyn(l)-1
94        do kl=1,nuvz
95        do ix=0,nxn(l)-1
96! FLEXPART_WRF -- use pph here
97!         ppmk=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
98          ppmk=pphn(ix,jy,kl,n,l)
99          thhn(ix,jy,kl,l)=tthn(ix,jy,kl,n,l)*(100000./ppmk)**kappa
100      enddo
101      enddo
102      enddo
103!
104! Loop over entire grid
105!**********************
106      do jy=0,nyn(l)-1
107! for FLEXPART_WRF, x & y coords are in meters
108! and true latitude varies with both i and j
109!       phi = (ylat0n(l) + jy * dyn(l)) * pi / 180.
110!       f = 0.00014585 * sin(phi)
111!       tanphi = tan(phi)
112!       cosphi = cos(phi)
113! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
114          jyvp=jy+1
115          jyvm=jy-1
116          if (jy.eq.0) jyvm=0
117          if (jy.eq.nyn(l)-1) jyvp=nyn(l)-1
118! Define absolute gap length
119          jumpy=2
120          if (jy.eq.0.or.jy.eq.nyn(l)-1) jumpy=1
121          juy=jumpy
122!
123        do ix=0,nxn(l)-1
124
125! for FLEXPART_WRF, x & y coords are in meters
126! and true latitude varies with both i and j
127          call xyindex_to_ll_wrf(  &
128                  l, real(ix), real(jy), dumlon, dumlat )
129          phi = dumlat * pi / 180.
130          f = 0.00014585 * sin(phi)
131          tanphi = tan(phi)
132          cosphi = cos(phi)
133
134! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
135          ixvp=ix+1
136          ixvm=ix-1
137          jumpx=2
138          if (ix.eq.0) ixvm=0
139          if (ix.eq.nxn(l)-1) ixvp=nxn(l)-1
140          ivrp=ixvp
141          ivrm=ixvm
142! Define absolute gap length
143          if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1
144          jux=jumpx
145! Precalculate pressure values for efficiency
146          do kl=1,nuvz
147! FLEXPART_WRF -- use pph here
148!           ppml(kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
149            ppml(kl)=pphn(ix,jy,kl,n,l)
150          enddo
151!
152! Loop over the vertical
153!***********************
154
155          do kl=1,nuvz
156            theta=thhn(ix,jy,kl,l)
157            klvrp=kl+1
158            klvrm=kl-1
159            klpt=kl
160! If top or bottom level, dthetadp is evaluated between the current
161! level and the level inside, otherwise between level+1 and level-1
162!
163            if (klvrp.gt.nuvz) klvrp=nuvz
164            if (klvrm.lt.1) klvrm=1
165            thetap=thhn(ix,jy,klvrp,l)
166            thetam=thhn(ix,jy,klvrm,l)
167            dthetadp=(thetap-thetam)/(ppml(klvrp)-ppml(klvrm))
168           
169! Compute vertical position at pot. temperature surface on subgrid
170! and the wind at that position
171!*****************************************************************
172! a) in x direction
173            ii=0
174            do i=ixvm,ixvp,jumpx
175              ivr=i
176              ii=ii+1
177! Search adjacent levels for current theta value
178! Spiral out from current level for efficiency
179              kup=klpt-1
180              kdn=klpt
181              kch=0
18240            continue
183! Upward branch
184              kup=kup+1
185              if (kch.ge.nlck) goto 21     ! No more levels to check,
186!                                            ! and no values found
187              if (kup.ge.nuvz) goto 41
188              kch=kch+1
189              k=kup
190              thdn=thhn(ivr,jy,k,l)
191              thup=thhn(ivr,jy,k+1,l)
192          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
193          ((thdn.le.theta).and.(thup.ge.theta))) then
194                  dt1=abs(theta-thdn)
195                  dt2=abs(theta-thup)
196                  dt=dt1+dt2
197                  if (dt.lt.eps) then   ! Avoid division by zero error
198                    dt1=0.5             ! G.W., 10.4.1996
199                    dt2=0.5
200                    dt=1.0
201                  endif
202        vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt
203                  goto 20
204                endif
20541            continue
206! Downward branch
207              kdn=kdn-1
208              if (kdn.lt.1) goto 40
209              kch=kch+1
210              k=kdn
211              thdn=thhn(ivr,jy,k,l)
212              thup=thhn(ivr,jy,k+1,l)
213          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
214          ((thdn.le.theta).and.(thup.ge.theta))) then
215                  dt1=abs(theta-thdn)
216                  dt2=abs(theta-thup)
217                  dt=dt1+dt2
218                  if (dt.lt.eps) then   ! Avoid division by zero error
219                    dt1=0.5             ! G.W., 10.4.1996
220                    dt2=0.5
221                    dt=1.0
222                  endif
223        vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt
224                  goto 20
225                endif
226                goto 40
227! This section used when no values were found
22821          continue
229! Must use vv at current level and long. jux becomes smaller by 1
230            vx(ii)=vvhn(ix,jy,kl,l)
231            jux=jux-1
232! Otherwise OK
23320          continue
234        end do
235
236          if (jux.gt.0) then
237! for FLEXPART_WRF, dx & dy are in meters.
238!         dvdx=(vx(2)-vx(1))/real(jux)/(dxn(l)*pi/180.)
239          dvdx=(vx(2)-vx(1))/real(jux)/dxn(l)
240          else
241          dvdx=vvhn(ivrp,jy,kl,l)-vvhn(ivrm,jy,kl,l)
242!         dvdx=dvdx/real(jumpx)/(dxn(l)*pi/180.)
243          dvdx=dvdx/real(jumpx)/dxn(l)
244! Only happens if no equivalent theta value
245! can be found on either side, hence must use values
246! from either side, same pressure level.
247          end if
248
249! b) in y direction
250
251            jj=0
252            do j=jyvm,jyvp,jumpy
253              jj=jj+1
254! Search adjacent levels for current theta value
255! Spiral out from current level for efficiency
256              kup=klpt-1
257              kdn=klpt
258              kch=0
25970            continue
260! Upward branch
261              kup=kup+1
262              if (kch.ge.nlck) goto 51     ! No more levels to check,
263!                                          ! and no values found
264              if (kup.ge.nuvz) goto 71
265              kch=kch+1
266              k=kup
267              thdn=thhn(ix,j,k,l)
268              thup=thhn(ix,j,k+1,l)
269          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
270          ((thdn.le.theta).and.(thup.ge.theta))) then
271                  dt1=abs(theta-thdn)
272                  dt2=abs(theta-thup)
273                  dt=dt1+dt2
274                  if (dt.lt.eps) then   ! Avoid division by zero error
275                    dt1=0.5             ! G.W., 10.4.1996
276                    dt2=0.5
277                    dt=1.0
278                  endif
279            uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt
280                  goto 50
281                endif
28271            continue
283! Downward branch
284              kdn=kdn-1
285              if (kdn.lt.1) goto 70
286              kch=kch+1
287              k=kdn
288              thdn=thhn(ix,j,k,l)
289              thup=thhn(ix,j,k+1,l)
290          if (((thdn.ge.theta).and.(thup.le.theta)).or. &
291          ((thdn.le.theta).and.(thup.ge.theta))) then
292                  dt1=abs(theta-thdn)
293                  dt2=abs(theta-thup)
294                  dt=dt1+dt2
295                  if (dt.lt.eps) then   ! Avoid division by zero error
296                    dt1=0.5             ! G.W., 10.4.1996
297                    dt2=0.5
298                    dt=1.0
299                  endif
300            uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt
301                  goto 50
302                endif
303                goto 70
304! This section used when no values were found
30551          continue
306! Must use uu at current level and lat. juy becomes smaller by 1
307            uy(jj)=uuhn(ix,jy,kl,l)
308            juy=juy-1
309! Otherwise OK
31050          continue
311        end do
312
313          if (juy.gt.0) then
314! for FLEXPART_WRF, dx & dy are in meters.
315!         dudy=(uy(2)-uy(1))/real(juy)/(dyn(l)*pi/180.)
316          dudy=(uy(2)-uy(1))/real(juy)/dyn(l)
317          else
318          dudy=uuhn(ix,jyvp,kl,l)-uuhn(ix,jyvm,kl,l)
319!         dudy=dudy/real(jumpy)/(dyn(l)*pi/180.)
320          dudy=dudy/real(jumpy)/dyn(l)
321          end if
322!
323! for FLEXPART_WRF, dx & dy are in meters.
324!   don't need to divide by r_earth when doing d/dy
325!   don't need to divide by r_earth*cosphi when doing d/dx
326! ??? I don't understand the uuhn*tanphi term, but leave it in for now ???
327! ??? What is the "-1.e6" factor ???
328!
329!         pvhn(ix,jy,kl,l)=dthetadp*(f+(dvdx/cosphi-dudy
330!    +    +uuhn(ix,jy,kl,l)*tanphi)/r_earth)*(-1.e6)*9.81
331          pvhn(ix,jy,kl,l)=dthetadp*( f + dvdx - dudy &
332          + (uuhn(ix,jy,kl,l)*tanphi/r_earth) )*(-1.e6)*9.81
333!
334! Resest jux and juy
335          jux=jumpx
336          juy=jumpy
337      end do
338    end do
339  end do
340  !
341end subroutine calcpv_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG