source: trunk/src/calcpar_nests.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 9.1 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 calcpar_nests(n,uuhn,vvhn,pvhn)
23  !                         i  i    i    o
24  !*****************************************************************************
25  !                                                                            *
26  !     Computation of several boundary layer parameters needed for the        *
27  !     dispersion calculation and calculation of dry deposition velocities.   *
28  !     All parameters are calculated over the entire grid.                    *
29  !     This routine is similar to calcpar, but is used for the nested grids.  *
30  !                                                                            *
31  !     Author: A. Stohl                                                       *
32  !                                                                            *
33  !     8 February 1999                                                        *
34  !                                                                            *
35  ! ------------------------------------------------------------------         *
36  !     Petra Seibert, Feb 2000:                                               *
37  !     convection scheme:                                                     *
38  !     new variables in call to richardson                                    *
39  !                                                                            *
40  !*****************************************************************************
41  !  Changes, Bernd C. Krueger, Feb. 2001:
42  !   Variables tth and qvh (on eta coordinates) in common block
43  !*****************************************************************************
44  !                                                                            *
45  ! Variables:                                                                 *
46  ! n                  temporal index for meteorological fields (1 to 3)       *
47  !                                                                            *
48  ! Constants:                                                                 *
49  !                                                                            *
50  !                                                                            *
51  ! Functions:                                                                 *
52  ! scalev             computation of ustar                                    *
53  ! obukhov            computatio of Obukhov length                            *
54  !                                                                            *
55  !*****************************************************************************
56
57  use par_mod
58  use com_mod
59
60  implicit none
61
62  integer :: n,ix,jy,i,l,kz,lz,kzmin
63  real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
64  real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
65  real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax)
66  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
67  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
68  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
69  real,parameter :: const=r_air/ga
70
71
72  ! Loop over all nests
73  !********************
74
75  do l=1,numbnests
76
77  ! Loop over entire grid
78  !**********************
79
80  do jy=0,nyn(l)-1
81
82  ! Set minimum height for tropopause
83  !**********************************
84
85    ylat=ylat0n(l)+real(jy)*dyn(l)
86    if ((ylat.ge.-20.).and.(ylat.le.20.)) then
87      altmin = 5000.
88    else
89      if ((ylat.gt.20.).and.(ylat.lt.40.)) then
90        altmin=2500.+(40.-ylat)*125.
91      else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
92        altmin=2500.+(40.+ylat)*125.
93      else
94        altmin=2500.
95      endif
96    endif
97
98    do ix=0,nxn(l)-1
99
100  ! 1) Calculation of friction velocity
101  !************************************
102
103      ustarn(ix,jy,1,n,l)=scalev(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
104           td2n(ix,jy,1,n,l),surfstrn(ix,jy,1,n,l))
105
106  ! 2) Calculation of inverse Obukhov length scale
107  !***********************************************
108
109      ol=obukhov(psn(ix,jy,1,n,l),tt2n(ix,jy,1,n,l), &
110           td2n(ix,jy,1,n,l),tthn(ix,jy,2,n,l),ustarn(ix,jy,1,n,l), &
111           sshfn(ix,jy,1,n,l),akm,bkm)
112      if (ol.ne.0.) then
113        olin(ix,jy,1,n,l)=1./ol
114      else
115        olin(ix,jy,1,n,l)=99999.
116      endif
117
118
119  ! 3) Calculation of convective velocity scale and mixing height
120  !**************************************************************
121
122      do i=1,nuvz
123        ulev(i)=uuhn(ix,jy,i,l)
124        vlev(i)=vvhn(ix,jy,i,l)
125        ttlev(i)=tthn(ix,jy,i,n,l)
126        qvlev(i)=qvhn(ix,jy,i,n,l)
127      end do
128
129      call richardson(psn(ix,jy,1,n,l),ustarn(ix,jy,1,n,l),ttlev, &
130           qvlev,ulev,vlev,nuvz,akz,bkz,sshfn(ix,jy,1,n,l), &
131           tt2n(ix,jy,1,n,l),td2n(ix,jy,1,n,l),hmixn(ix,jy,1,n,l), &
132           wstarn(ix,jy,1,n,l),hmixplus)
133
134      if(lsubgrid.eq.1) then
135        subsceff=min(excessoron(ix,jy,l),hmixplus)
136      else
137        subsceff=0
138      endif
139  !
140  ! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
141  !
142      hmixn(ix,jy,1,n,l)=hmixn(ix,jy,1,n,l)+subsceff
143      hmixn(ix,jy,1,n,l)=max(hmixmin,hmixn(ix,jy,1,n,l)) ! minim PBL height
144      hmixn(ix,jy,1,n,l)=min(hmixmax,hmixn(ix,jy,1,n,l)) ! maxim PBL height
145
146
147  ! 4) Calculation of dry deposition velocities
148  !********************************************
149
150      if (DRYDEP) then
151        z0(4)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
152        z0(9)=0.016*ustarn(ix,jy,1,n,l)*ustarn(ix,jy,1,n,l)/ga
153
154  ! Calculate relative humidity at surface
155  !***************************************
156        rh=ew(td2n(ix,jy,1,n,l))/ew(tt2n(ix,jy,1,n,l))
157
158        call getvdep_nests(n,ix,jy,ustarn(ix,jy,1,n,l), &
159             tt2n(ix,jy,1,n,l),psn(ix,jy,1,n,l),1./olin(ix,jy,1,n,l), &
160             ssrn(ix,jy,1,n,l),rh,lsprecn(ix,jy,1,n,l)+ &
161             convprecn(ix,jy,1,n,l),sdn(ix,jy,1,n,l),vd,l)
162
163        do i=1,nspec
164          vdepn(ix,jy,i,n,l)=vd(i)
165        end do
166
167      endif
168
169  !******************************************************
170  ! Calculate height of thermal tropopause (Hoinka, 1997)
171  !******************************************************
172
173  ! 1) Calculate altitudes of ECMWF model levels
174  !*********************************************
175
176      tvold=tt2n(ix,jy,1,n,l)*(1.+0.378*ew(td2n(ix,jy,1,n,l))/ &
177           psn(ix,jy,1,n,l))
178      pold=psn(ix,jy,1,n,l)
179      zold=0.
180      do kz=2,nuvz
181        pint=akz(kz)+bkz(kz)*psn(ix,jy,1,n,l)  ! pressure on model layers
182        tv=tthn(ix,jy,kz,n,l)*(1.+0.608*qvhn(ix,jy,kz,n,l))
183
184        if (abs(tv-tvold).gt.0.2) then
185         zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
186        else
187          zlev(kz)=zold+const*log(pold/pint)*tv
188        endif
189        tvold=tv
190        pold=pint
191        zold=zlev(kz)
192      end do
193
194  ! 2) Define a minimum level kzmin, from which upward the tropopause is
195  !    searched for. This is to avoid inversions in the lower troposphere
196  !    to be identified as the tropopause
197  !************************************************************************
198
199      do kz=1,nuvz
200        if (zlev(kz).ge.altmin) then
201          kzmin=kz
202          goto 45
203        endif
204      end do
20545    continue
206
207  ! 3) Search for first stable layer above minimum height that fulfills the
208  !    thermal tropopause criterion
209  !************************************************************************
210
211      do kz=kzmin,nuvz
212        do lz=kz+1,nuvz
213          if ((zlev(lz)-zlev(kz)).gt.2000.) then
214            if (((tthn(ix,jy,kz,n,l)-tthn(ix,jy,lz,n,l))/ &
215                 (zlev(lz)-zlev(kz))).lt.0.002) then
216              tropopausen(ix,jy,1,n,l)=zlev(kz)
217              goto 51
218            endif
219            goto 50
220          endif
221        end do
22250      continue
223      end do
22451    continue
225
226
227    end do
228  end do
229
230
231    call calcpv_nests(l,n,uuhn,vvhn,pvhn)
232
233  end do
234
235
236end subroutine calcpar_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG