source: trunk/src/calcpar_gfs.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 9.4 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(n,uuh,vvh,pvh)
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  !                                                                            *
30  !     Author: A. Stohl                                                       *
31  !                                                                            *
32  !     21 May 1995                                                            *
33  !                                                                            *
34  ! ------------------------------------------------------------------         *
35  !     Petra Seibert, Feb 2000:                                               *
36  !     convection scheme:                                                     *
37  !     new variables in call to richardson                                    *
38  !                                                                            *
39  !*****************************************************************************
40  !  Changes, Bernd C. Krueger, Feb. 2001:
41  !   Variables tth and qvh (on eta coordinates) in common block
42  !*****************************************************************************
43  !                                                                            *
44  !  CHANGE 17/11/2005 Caroline Forster NCEP GFS version                       *
45  !                                                                            *
46  ! Variables:                                                                 *
47  ! n                  temporal index for meteorological fields (1 to 3)       *
48  !                                                                            *
49  ! Constants:                                                                 *
50  !                                                                            *
51  !                                                                            *
52  ! Functions:                                                                 *
53  ! scalev             computation of ustar                                    *
54  ! obukhov            computatio of Obukhov length                            *
55  !                                                                            *
56  !*****************************************************************************
57
58  use par_mod
59  use com_mod
60
61  implicit none
62
63  integer :: n,ix,jy,i,kz,lz,kzmin,llev
64  real :: ttlev(nuvzmax),qvlev(nuvzmax),obukhov,scalev,ol,hmixplus
65  real :: ulev(nuvzmax),vlev(nuvzmax),ew,rh,vd(maxspec),subsceff,ylat
66  real :: altmin,tvold,pold,zold,pint,tv,zlev(nuvzmax),hmixdummy
67  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
68  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
69  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
70  real,parameter :: const=r_air/ga
71
72
73  ! Loop over entire grid
74  !**********************
75
76  do jy=0,nymin1
77
78  ! Set minimum height for tropopause
79  !**********************************
80
81    ylat=ylat0+real(jy)*dy
82    if ((ylat.ge.-20.).and.(ylat.le.20.)) then
83      altmin = 5000.
84    else
85      if ((ylat.gt.20.).and.(ylat.lt.40.)) then
86        altmin=2500.+(40.-ylat)*125.
87      else if ((ylat.gt.-40.).and.(ylat.lt.-20.)) then
88        altmin=2500.+(40.+ylat)*125.
89      else
90        altmin=2500.
91      endif
92    endif
93
94    do ix=0,nxmin1
95
96  ! 1) Calculation of friction velocity
97  !************************************
98
99      ustar(ix,jy,1,n)=scalev(ps(ix,jy,1,n),tt2(ix,jy,1,n), &
100           td2(ix,jy,1,n),surfstr(ix,jy,1,n))
101      if (ustar(ix,jy,1,n).le.1.e-8) ustar(ix,jy,1,n)=1.e-8
102
103  ! 2) Calculation of inverse Obukhov length scale
104  !***********************************************
105
106  ! NCEP version: find first level above ground
107      llev = 0
108      do i=1,nuvz
109       if (ps(ix,jy,1,n).lt.akz(i)) llev=i
110      end do
111       llev = llev+1
112       if (llev.gt.nuvz) llev = nuvz-1
113  ! NCEP version
114
115  ! calculate inverse Obukhov length scale with tth(llev)
116      ol=obukhov(ps(ix,jy,1,n),tt2(ix,jy,1,n),td2(ix,jy,1,n), &
117           tth(ix,jy,llev,n),ustar(ix,jy,1,n),sshf(ix,jy,1,n),akz(llev))
118      if (ol.ne.0.) then
119        oli(ix,jy,1,n)=1./ol
120      else
121        oli(ix,jy,1,n)=99999.
122      endif
123
124
125  ! 3) Calculation of convective velocity scale and mixing height
126  !**************************************************************
127
128      do i=1,nuvz
129        ulev(i)=uuh(ix,jy,i)
130        vlev(i)=vvh(ix,jy,i)
131        ttlev(i)=tth(ix,jy,i,n)
132        qvlev(i)=qvh(ix,jy,i,n)
133      end do
134
135  ! NCEP version hmix has been read in in readwind.f, is therefore not calculated here
136      call richardson(ps(ix,jy,1,n),ustar(ix,jy,1,n),ttlev,qvlev, &
137           ulev,vlev,nuvz,akz,bkz,sshf(ix,jy,1,n),tt2(ix,jy,1,n), &
138           td2(ix,jy,1,n),hmixdummy,wstar(ix,jy,1,n),hmixplus)
139
140      if(lsubgrid.eq.1) then
141        subsceff=min(excessoro(ix,jy),hmixplus)
142      else
143        subsceff=0
144      endif
145  !
146  ! CALCULATE HMIX EXCESS ACCORDING TO SUBGRIDSCALE VARIABILITY AND STABILITY
147  !
148      hmix(ix,jy,1,n)=hmix(ix,jy,1,n)+subsceff
149      hmix(ix,jy,1,n)=max(hmixmin,hmix(ix,jy,1,n)) ! set minimum PBL height
150      hmix(ix,jy,1,n)=min(hmixmax,hmix(ix,jy,1,n)) ! set maximum PBL height
151
152  ! 4) Calculation of dry deposition velocities
153  !********************************************
154
155      if (DRYDEP) then
156  ! Sabine Eckhardt, Dec 06: use new index for z0 for water depending on
157  ! windspeed
158        z0(7)=0.016*ustar(ix,jy,1,n)*ustar(ix,jy,1,n)/ga
159
160  ! Calculate relative humidity at surface
161  !***************************************
162        rh=ew(td2(ix,jy,1,n))/ew(tt2(ix,jy,1,n))
163
164        call getvdep(n,ix,jy,ustar(ix,jy,1,n), &
165             tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), &
166             ssr(ix,jy,1,n),rh,lsprec(ix,jy,1,n)+convprec(ix,jy,1,n), &
167             sd(ix,jy,1,n),vd)
168
169        do i=1,nspec
170          vdep(ix,jy,i,n)=vd(i)
171        end do
172
173      endif
174
175  !******************************************************
176  ! Calculate height of thermal tropopause (Hoinka, 1997)
177  !******************************************************
178
179  ! 1) Calculate altitudes of NCEP model levels
180  !*********************************************
181
182      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
183           ps(ix,jy,1,n))
184      pold=ps(ix,jy,1,n)
185      zold=0.
186      do kz=llev,nuvz
187        pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)  ! pressure on model layers
188        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
189
190        if (abs(tv-tvold).gt.0.2) then
191         zlev(kz)=zold+const*log(pold/pint)*(tv-tvold)/log(tv/tvold)
192        else
193          zlev(kz)=zold+const*log(pold/pint)*tv
194        endif
195        tvold=tv
196        pold=pint
197        zold=zlev(kz)
198      end do
199
200  ! 2) Define a minimum level kzmin, from which upward the tropopause is
201  !    searched for. This is to avoid inversions in the lower troposphere
202  !    to be identified as the tropopause
203  !************************************************************************
204
205      do kz=llev,nuvz
206        if (zlev(kz).ge.altmin) then
207          kzmin=kz
208          goto 45
209        endif
210      end do
21145    continue
212
213  ! 3) Search for first stable layer above minimum height that fulfills the
214  !    thermal tropopause criterion
215  !************************************************************************
216
217      do kz=kzmin,nuvz
218        do lz=kz+1,nuvz
219          if ((zlev(lz)-zlev(kz)).gt.2000.) then
220            if (((tth(ix,jy,kz,n)-tth(ix,jy,lz,n))/ &
221                 (zlev(lz)-zlev(kz))).lt.0.002) then
222              tropopause(ix,jy,1,n)=zlev(kz)
223              goto 51
224            endif
225            goto 50
226          endif
227        end do
22850      continue
229      end do
23051    continue
231
232
233    end do
234  end do
235
236
237  ! Calculation of potential vorticity on 3-d grid
238  !***********************************************
239
240  call calcpv(n,uuh,vvh,pvh)
241
242
243end subroutine calcpar
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG