source: trunk/src/readwind_nests_emos.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
  • Property svn:executable set to *
File size: 11.7 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 readwind_nests(indj,n,uuhn,vvhn,wwhn)
23  !                           i   i  o    o    o
24  !*****************************************************************************
25  !                                                                            *
26  !     This routine reads the wind fields for the nested model domains.       *
27  !     It is similar to subroutine readwind, which reads the mother domain.   *
28  !                                                                            *
29  !     Authors: A. Stohl, G. Wotawa                                           *
30  !                                                                            *
31  !     8 February 1999                                                        *
32  !                                                                            *
33  !     Last update: 17 October 2000, A. Stohl                                 *
34  !                                                                            *
35  !*****************************************************************************
36  !  Changes, Bernd C. Krueger, Feb. 2001:
37  !   Variables tthn and qvhn (on eta coordinates) in common block
38  !*****************************************************************************
39
40  use par_mod
41  use com_mod
42
43  implicit none
44
45  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
46  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
47  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
48  integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax,lunit,l
49
50  ! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
51
52  ! dimension of isec2 at least (22+n), where n is the number of parallels or
53  ! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
54
55  ! dimension of zsec2 at least (10+nn), where nn is the number of vertical
56  ! coordinate parameters
57
58  integer :: isec0(2),isec1(56),isec2(22+nxmaxn+nymaxn),isec3(2)
59  integer :: isec4(64),inbuff(jpack),ilen,ierr,iword
60  !integer iswap
61  real :: zsec2(60+2*nuvzmax),zsec3(2),zsec4(jpunp)
62  real :: xaux,yaux,xaux0,yaux0
63  real :: ewss(0:nxmaxn-1,0:nymaxn-1),nsss(0:nxmaxn-1,0:nymaxn-1)
64  real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
65
66  character(len=1) :: yoper = 'D'
67  logical :: hflswitch,strswitch
68
69
70  do l=1,numbnests
71    hflswitch=.false.
72    strswitch=.false.
73    levdiff2=nlev_ec-nwz+1
74    iumax=0
75    iwmax=0
76
77  !
78  ! OPENING OF DATA FILE (GRIB CODE)
79  !
805   call pbopen(lunit,path(numpath+2*(l-1)+1) &
81         (1:length(numpath+2*(l-1)+1))//wfnamen(l,indj),'r',ierr)
82    if(ierr.lt.0) goto 999
83
84    ifield=0
8510    ifield=ifield+1
86  !
87  ! GET NEXT FIELDS
88  !
89    call pbgrib(lunit,inbuff,jpack,ilen,ierr)
90    if(ierr.eq.-1) goto 50    ! EOF DETECTED
91    if(ierr.lt.-1) goto 888   ! ERROR DETECTED
92
93    ierr=1
94
95  ! Check whether we are on a little endian or on a big endian computer
96  !********************************************************************
97
98  !  if (inbuff(1).eq.1112101447) then         ! little endian, swap bytes
99  !    iswap=1+ilen/4
100  !    call swap32(inbuff,iswap)
101  !  else if (inbuff(1).ne.1196575042) then    ! big endian
102  !    stop 'subroutine gridcheck: corrupt GRIB data'
103  !  endif
104
105    call gribex(isec0,isec1,isec2,zsec2,isec3,zsec3,isec4, &
106         zsec4,jpunp,inbuff,jpack,iword,yoper,ierr)
107    if (ierr.ne.0) goto 888   ! ERROR DETECTED
108
109    if(ifield.eq.1) then
110
111  ! CHECK GRID SPECIFICATIONS
112
113      if(isec2(2).ne.nxn(l)) stop &
114           'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
115      if(isec2(3).ne.nyn(l)) stop &
116           'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
117      if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&
118           &IZATION NOT CONSISTENT FOR A NESTING LEVEL'
119      xaux=real(isec2(5))/1000.
120      yaux=real(isec2(7))/1000.
121      xaux0=xlon0n(l)
122      yaux0=ylat0n(l)
123      if(xaux.lt.0.) xaux=xaux+360.
124      if(yaux.lt.0.) yaux=yaux+360.
125      if(xaux0.lt.0.) xaux0=xaux0+360.
126      if(yaux0.lt.0.) yaux0=yaux0+360.
127      if(xaux.ne.xaux0) &
128           stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT FOR A NES&
129           &TING LEVEL'
130      if(yaux.ne.yaux0) &
131           stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT FOR A NEST&
132           &ING LEVEL'
133    endif
134
135
136    do j=0,nyn(l)-1
137      do i=0,nxn(l)-1
138        k=isec1(8)
139        if(isec1(6).eq.130) tthn(i,j,nlev_ec-k+2,n,l)= &!! TEMPERATURE
140             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
141        if(isec1(6).eq.131) uuhn(i,j,nlev_ec-k+2,l)= &!! U VELOCITY
142             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
143        if(isec1(6).eq.132) vvhn(i,j,nlev_ec-k+2,l)= &!! V VELOCITY
144             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
145        if(isec1(6).eq.133) then                         !! SPEC. HUMIDITY
146          qvhn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
147          if (qvhn(i,j,nlev_ec-k+2,n,l) .lt. 0.) &
148               qvhn(i,j,nlev_ec-k+2,n,l) = 0.
149  !          this is necessary because the gridded data may contain
150  !          spurious negative values
151        endif
152        if(isec1(6).eq.134) psn(i,j,1,n,l)= &!! SURF. PRESS.
153             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
154        if(isec1(6).eq.135) wwhn(i,j,nlev_ec-k+1,l)= &!! W VELOCITY
155             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
156        if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH
157             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
158        if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
159             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
160        if(isec1(6).eq.164) tccn(i,j,1,n,l)= &!! CLOUD COVER
161             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
162        if(isec1(6).eq.165) u10n(i,j,1,n,l)= &!! 10 M U VELOCITY
163             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
164        if(isec1(6).eq.166) v10n(i,j,1,n,l)= &!! 10 M V VELOCITY
165             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
166        if(isec1(6).eq.167) tt2n(i,j,1,n,l)= &!! 2 M TEMPERATURE
167             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
168        if(isec1(6).eq.168) td2n(i,j,1,n,l)= &!! 2 M DEW POINT
169             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
170        if(isec1(6).eq.142) then                         !! LARGE SCALE PREC.
171          lsprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
172          if (lsprecn(i,j,1,n,l).lt.0.) lsprecn(i,j,1,n,l)=0.
173        endif
174        if(isec1(6).eq.143) then                         !! CONVECTIVE PREC.
175          convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
176          if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
177        endif
178        if(isec1(6).eq.146) sshfn(i,j,1,n,l)= &!! SENS. HEAT FLUX
179             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
180        if((isec1(6).eq.146).and. &
181             (zsec4(nxn(l)*(nyn(l)-j-1)+i+1).ne.0.)) hflswitch=.true.    ! Heat flux available
182        if(isec1(6).eq.176) then                         !! SOLAR RADIATION
183          ssrn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
184          if (ssrn(i,j,1,n,l).lt.0.) ssrn(i,j,1,n,l)=0.
185        endif
186        if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
187             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
188        if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
189             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
190        if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
191             (zsec4(nxn(l)*(nyn(l)-j-1)+i+1).ne.0.)) strswitch=.true.    ! stress available
192        if(isec1(6).eq.129) oron(i,j,l)= &!! ECMWF OROGRAPHY
193             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
194        if(isec1(6).eq.160) excessoron(i,j,l)= &!! STANDARD DEVIATION OF OROGRAPHY
195             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
196        if(isec1(6).eq.172) lsmn(i,j,l)= &!! ECMWF LAND SEA MASK
197             zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
198        if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
199        if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
200      end do
201    end do
202
203    goto 10                      !! READ NEXT LEVEL OR PARAMETER
204  !
205  ! CLOSING OF INPUT DATA FILE
206  !
20750   call pbclose(lunit,ierr)     !! FINISHED READING / CLOSING GRIB FILE
208
209    if(levdiff2.eq.0) then
210      iwmax=nlev_ec+1
211      do i=0,nxn(l)-1
212        do j=0,nyn(l)-1
213          wwhn(i,j,nlev_ec+1,l)=0.
214        end do
215      end do
216    endif
217
218    do i=0,nxn(l)-1
219      do j=0,nyn(l)-1
220        surfstrn(i,j,1,n,l)=sqrt(ewss(i,j)**2+nsss(i,j)**2)
221      end do
222    end do
223
224    if ((.not.hflswitch).or.(.not.strswitch)) then
225      write(*,*) 'WARNING: No flux data contained in GRIB file ', &
226           wfnamen(l,indj)
227
228  ! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
229  ! As ECMWF has increased the model resolution, such that now the first model
230  ! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
231  ! (3rd model level in FLEXPART) for the profile method
232  !***************************************************************************
233
234      do i=0,nxn(l)-1
235        do j=0,nyn(l)-1
236        plev1=akz(3)+bkz(3)*psn(i,j,1,n,l)
237        pmean=0.5*(psn(i,j,1,n,l)+plev1)
238        tv=tthn(i,j,3,n,l)*(1.+0.61*qvhn(i,j,3,n,l))
239        fu=-r_air*tv/ga/pmean
240        hlev1=fu*(plev1-psn(i,j,1,n,l))   ! HEIGTH OF FIRST MODEL LAYER
241        ff10m= sqrt(u10n(i,j,1,n,l)**2+v10n(i,j,1,n,l)**2)
242        fflev1=sqrt(uuhn(i,j,3,l)**2+vvhn(i,j,3,l)**2)
243        call pbl_profile(psn(i,j,1,n,l),td2n(i,j,1,n,l),hlev1, &
244             tt2n(i,j,1,n,l),tthn(i,j,3,n,l),ff10m,fflev1, &
245             surfstrn(i,j,1,n,l),sshfn(i,j,1,n,l))
246        if(sshfn(i,j,1,n,l).gt.200.) sshfn(i,j,1,n,l)=200.
247        if(sshfn(i,j,1,n,l).lt.-400.) sshfn(i,j,1,n,l)=-400.
248        end do
249      end do
250    endif
251
252
253  ! Assign 10 m wind to model level at eta=1.0 to have one additional model
254  ! level at the ground
255  ! Specific humidity is taken the same as at one level above
256  ! Temperature is taken as 2 m temperature
257  !**************************************************************************
258
259    do i=0,nxn(l)-1
260      do j=0,nyn(l)-1
261        uuhn(i,j,1,l)=u10n(i,j,1,n,l)
262        vvhn(i,j,1,l)=v10n(i,j,1,n,l)
263        qvhn(i,j,1,n,l)=qvhn(i,j,2,n,l)
264        tthn(i,j,1,n,l)=tt2n(i,j,1,n,l)
265      end do
266    end do
267
268    if(iumax.ne.nuvz-1) stop &
269         'READWIND: NUVZ NOT CONSISTENT FOR A NESTING LEVEL'
270    if(iwmax.ne.nwz) stop &
271         'READWIND: NWZ NOT CONSISTENT FOR A NESTING LEVEL'
272
273  end do
274
275  return
276888   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
277  write(*,*) ' #### ',wfnamen(l,indj),' FOR NESTING LEVEL  #### '
278  write(*,*) ' #### ',l,' IS NOT GRIB FORMAT !!!           #### '
279  stop 'Execution terminated'
280
281
282999   write(*,*) ' #### FLEXPART MODEL ERROR! WINDFIELD         #### '
283  write(*,*) ' #### ',wfnamen(l,indj),'                    #### '
284  write(*,*) ' #### CANNOT BE OPENED FOR NESTING LEVEL ',l,'####'
285
286end subroutine readwind_nests
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG