source: trunk/src/calcpv_nests.f90 @ 28

Last change on this file since 28 was 24, checked in by igpis, 8 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