source: trunk/src/calcpv_nests.f90 @ 28

Last change on this file since 28 was 24, checked in by igpis, 10 years ago

version 9.2 beta. Changes from HH, AST, MC, NIK, IP. Changes in vert transform. New SPECIES input includes scavenging coefficients

File size: 9.9 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 calcpv_nests(l,n,uuhn,vvhn,pvhn)
23  !                     i i  i    i    o
24  !*****************************************************************************
25  !                                                                            *
26  !  Calculation of potential vorticity on 3-d nested grid                     *
27  !                                                                            *
28  !     Author: P. James                                                       *
29  !     22 February 2000                                                       *
30  !                                                                            *
31  !*****************************************************************************
32  !                                                                            *
33  ! Variables:                                                                 *
34  ! n                  temporal index for meteorological fields (1 to 2)       *
35  ! l                  index of current nest                                   *
36  !                                                                            *
37  ! Constants:                                                                 *
38  !                                                                            *
39  !*****************************************************************************
40
41  use par_mod
42  use com_mod
43
44  implicit none
45
46  integer :: n,ix,jy,i,j,k,kl,ii,jj,klvrp,klvrm,klpt,kup,kdn,kch
47  integer :: jyvp,jyvm,ixvp,ixvm,jumpx,jumpy,jux,juy,ivrm,ivrp,ivr
48  integer :: nlck,l
49  real :: vx(2),uy(2),phi,tanphi,cosphi,dvdx,dudy,f
50  real :: theta,thetap,thetam,dthetadp,dt1,dt2,dt
51  real :: ppml(0:nxmaxn-1,0:nymaxn-1,nuvzmax),ppmk(0:nxmaxn-1,0:nymaxn-1,nuvzmax)
52  real :: thup,thdn
53  real,parameter :: eps=1.e-5,p0=101325
54  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
55  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
56  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
57
58  ! Set number of levels to check for adjacent theta
59  nlck=nuvz/3
60  !
61  ! Loop over entire grid
62  !**********************
63  do kl=1,nuvz
64    do jy=0,nyn(l)-1
65      do ix=0,nxn(l)-1
66         ppml(ix,jy,kl)=akz(kl)+bkz(kl)*psn(ix,jy,1,n,l)
67      enddo
68    enddo
69  enddo
70  ppmk=(100000./ppml)**kappa
71
72  do jy=0,nyn(l)-1
73    phi = (ylat0n(l) + jy * dyn(l)) * pi / 180.
74    f = 0.00014585 * sin(phi)
75    tanphi = tan(phi)
76    cosphi = cos(phi)
77  ! Provide a virtual jy+1 and jy-1 in case we are on domain edge (Lat)
78      jyvp=jy+1
79      jyvm=jy-1
80      if (jy.eq.0) jyvm=0
81      if (jy.eq.nyn(l)-1) jyvp=nyn(l)-1
82  ! Define absolute gap length
83      jumpy=2
84      if (jy.eq.0.or.jy.eq.nyn(l)-1) jumpy=1
85      juy=jumpy
86  !
87    do ix=0,nxn(l)-1
88  ! Provide a virtual ix+1 and ix-1 in case we are on domain edge (Long)
89      ixvp=ix+1
90      ixvm=ix-1
91      jumpx=2
92      if (ix.eq.0) ixvm=0
93      if (ix.eq.nxn(l)-1) ixvp=nxn(l)-1
94      ivrp=ixvp
95      ivrm=ixvm
96  ! Define absolute gap length
97      if (ix.eq.0.or.ix.eq.nxn(l)-1) jumpx=1
98      jux=jumpx
99  !
100  ! Loop over the vertical
101  !***********************
102
103      do kl=1,nuvz
104        theta=tthn(ix,jy,kl,n,l)*ppmk(ix,jy,kl)
105        klvrp=kl+1
106        klvrm=kl-1
107        klpt=kl
108  ! If top or bottom level, dthetadp is evaluated between the current
109  ! level and the level inside, otherwise between level+1 and level-1
110  !
111        if (klvrp.gt.nuvz) klvrp=nuvz
112        if (klvrm.lt.1) klvrm=1
113        thetap=tthn(ix,jy,klvrp,n,l)*ppmk(ix,jy,klvrp)
114        thetam=tthn(ix,jy,klvrm,n,l)*ppmk(ix,jy,klvrm)
115        dthetadp=(thetap-thetam)/(ppml(ix,jy,klvrp)-ppml(ix,jy,klvrm))
116
117  ! Compute vertical position at pot. temperature surface on subgrid
118  ! and the wind at that position
119  !*****************************************************************
120  ! a) in x direction
121        ii=0
122        do i=ixvm,ixvp,jumpx
123          ivr=i
124          ii=ii+1
125  ! Search adjacent levels for current theta value
126  ! Spiral out from current level for efficiency
127          kup=klpt-1
128          kdn=klpt
129          kch=0
13040        continue
131  ! Upward branch
132          kup=kup+1
133          if (kch.ge.nlck) goto 21     ! No more levels to check,
134  !                                       ! and no values found
135          if (kup.ge.nuvz) goto 41
136          kch=kch+1
137          k=kup
138          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
139          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
140
141      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
142           ((thdn.le.theta).and.(thup.ge.theta))) then
143              dt1=abs(theta-thdn)
144              dt2=abs(theta-thup)
145              dt=dt1+dt2
146              if (dt.lt.eps) then   ! Avoid division by zero error
147                dt1=0.5             ! G.W., 10.4.1996
148                dt2=0.5
149                dt=1.0
150              endif
151    vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt
152              goto 20
153            endif
15441        continue
155  ! Downward branch
156          kdn=kdn-1
157          if (kdn.lt.1) goto 40
158          kch=kch+1
159          k=kdn
160          thdn=tthn(ivr,jy,k,n,l)*ppmk(ivr,jy,k)
161          thup=tthn(ivr,jy,k+1,n,l)*ppmk(ivr,jy,k+1)
162      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
163           ((thdn.le.theta).and.(thup.ge.theta))) then
164              dt1=abs(theta-thdn)
165              dt2=abs(theta-thup)
166              dt=dt1+dt2
167              if (dt.lt.eps) then   ! Avoid division by zero error
168                dt1=0.5             ! G.W., 10.4.1996
169                dt2=0.5
170                dt=1.0
171              endif
172    vx(ii)=(vvhn(ivr,jy,k,l)*dt2+vvhn(ivr,jy,k+1,l)*dt1)/dt
173              goto 20
174            endif
175            goto 40
176  ! This section used when no values were found
17721      continue
178  ! Must use vv at current level and long. jux becomes smaller by 1
179        vx(ii)=vvhn(ix,jy,kl,l)
180        jux=jux-1
181  ! Otherwise OK
18220        continue
183        end do
184      if (jux.gt.0) then
185      dvdx=(vx(2)-vx(1))/real(jux)/(dxn(l)*pi/180.)
186      else
187      dvdx=vvhn(ivrp,jy,kl,l)-vvhn(ivrm,jy,kl,l)
188      dvdx=dvdx/real(jumpx)/(dxn(l)*pi/180.)
189  ! Only happens if no equivalent theta value
190  ! can be found on either side, hence must use values
191  ! from either side, same pressure level.
192      end if
193
194  ! b) in y direction
195
196        jj=0
197        do j=jyvm,jyvp,jumpy
198          jj=jj+1
199  ! Search adjacent levels for current theta value
200  ! Spiral out from current level for efficiency
201          kup=klpt-1
202          kdn=klpt
203          kch=0
20470        continue
205  ! Upward branch
206          kup=kup+1
207          if (kch.ge.nlck) goto 51     ! No more levels to check,
208  !                                     ! and no values found
209          if (kup.ge.nuvz) goto 71
210          kch=kch+1
211          k=kup
212          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
213          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
214      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
215           ((thdn.le.theta).and.(thup.ge.theta))) then
216              dt1=abs(theta-thdn)
217              dt2=abs(theta-thup)
218              dt=dt1+dt2
219              if (dt.lt.eps) then   ! Avoid division by zero error
220                dt1=0.5             ! G.W., 10.4.1996
221                dt2=0.5
222                dt=1.0
223              endif
224        uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt
225              goto 50
226            endif
22771        continue
228  ! Downward branch
229          kdn=kdn-1
230          if (kdn.lt.1) goto 70
231          kch=kch+1
232          k=kdn
233          thdn=tthn(ix,j,k,n,l)*ppmk(ix,j,k)
234          thup=tthn(ix,j,k+1,n,l)*ppmk(ix,j,k+1)
235      if (((thdn.ge.theta).and.(thup.le.theta)).or. &
236           ((thdn.le.theta).and.(thup.ge.theta))) then
237              dt1=abs(theta-thdn)
238              dt2=abs(theta-thup)
239              dt=dt1+dt2
240              if (dt.lt.eps) then   ! Avoid division by zero error
241                dt1=0.5             ! G.W., 10.4.1996
242                dt2=0.5
243                dt=1.0
244              endif
245        uy(jj)=(uuhn(ix,j,k,l)*dt2+uuhn(ix,j,k+1,l)*dt1)/dt
246              goto 50
247            endif
248            goto 70
249  ! This section used when no values were found
25051      continue
251  ! Must use uu at current level and lat. juy becomes smaller by 1
252        uy(jj)=uuhn(ix,jy,kl,l)
253        juy=juy-1
254  ! Otherwise OK
25550        continue
256        end do
257      if (juy.gt.0) then
258      dudy=(uy(2)-uy(1))/real(juy)/(dyn(l)*pi/180.)
259      else
260      dudy=uuhn(ix,jyvp,kl,l)-uuhn(ix,jyvm,kl,l)
261      dudy=dudy/real(jumpy)/(dyn(l)*pi/180.)
262      end if
263  !
264      pvhn(ix,jy,kl,l)=dthetadp*(f+(dvdx/cosphi-dudy &
265           +uuhn(ix,jy,kl,l)*tanphi)/r_earth)*(-1.e6)*9.81
266
267  !
268  ! Resest jux and juy
269      jux=jumpx
270      juy=jumpy
271      end do
272    end do
273  end do
274  !
275end subroutine calcpv_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG