Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

readwind_nests.f90 @ 23

Last change on this file since 23 was 23, checked in by igpis, 10 years ago
start tracking test environment directory FP_AI
File size: 18.3 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
22	subroutine 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	! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with ECMWF grib_api *
39	! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with ECMWF grib_api *
40	!*****************************************************************************
41
42	use grib_api
43	use par_mod
44	use com_mod
45
46	implicit none
47
48	!HSO parameters for grib_api
49	integer :: ifile
50	integer :: iret
51	integer :: igrib
52	integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
53	integer :: parId !!added by mc for making it consistent with new readwind.f90
54	integer :: gotGrid
55	!HSO end
56
57	real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
58	real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
59	real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
60	integer :: indj,i,j,k,n,levdiff2,ifield,iumax,iwmax,l
61
62	! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
63
64	! dimension of isec2 at least (22+n), where n is the number of parallels or
65	! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
66
67	! dimension of zsec2 at least (10+nn), where nn is the number of vertical
68	! coordinate parameters
69
70	integer :: isec1(56),isec2(22+nxmaxn+nymaxn)
71	real(kind=4) :: zsec4(jpunp)
72	real(kind=4) :: xaux,yaux
73	real(kind=8) :: xauxin,yauxin
74	real,parameter :: eps=1.e-4
75	real :: ewss(0:nxmaxn-1,0:nymaxn-1),nsss(0:nxmaxn-1,0:nymaxn-1)
76	real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
77	real :: conversion_factor !added by mc to make it consistent with new gridchek.f90
78
79	logical :: hflswitch,strswitch
80
81	!HSO grib api error messages
82	character(len=24) :: gribErrorMsg = 'Error reading grib file'
83	character(len=20) :: gribFunction = 'readwind_nests'
84
85	do l=1,numbnests
86	hflswitch=.false.
87	strswitch=.false.
88	levdiff2=nlev_ec-nwz+1
89	iumax=0
90	iwmax=0
91
92	ifile=0
93	igrib=0
94	iret=0
95
96	!
97	! OPENING OF DATA FILE (GRIB CODE)
98	!
99
100	5 call grib_open_file(ifile,path(numpath+2*(l-1)+1) &
101	(1:length(numpath+2*(l-1)+1))//trim(wfnamen(l,indj)),'r')
102	if (iret.ne.GRIB_SUCCESS) then
103	goto 888 ! ERROR DETECTED
104	endif
105	!turn on support for multi fields messages */
106	!call grib_multi_support_on
107
108	gotGrid=0
109	ifield=0
110	10 ifield=ifield+1
111	!
112	! GET NEXT FIELDS
113	!
114	call grib_new_from_file(ifile,igrib,iret)
115	if (iret.eq.GRIB_END_OF_FILE) then
116	goto 50 ! EOF DETECTED
117	elseif (iret.ne.GRIB_SUCCESS) then
118	goto 888 ! ERROR DETECTED
119	endif
120
121	!first see if we read GRIB1 or GRIB2
122	call grib_get_int(igrib,'editionNumber',gribVer,iret)
123	call grib_check(iret,gribFunction,gribErrorMsg)
124
125	if (gribVer.eq.1) then ! GRIB Edition 1
126
127	!print*,'GRiB Edition 1'
128	!read the grib2 identifiers
129	call grib_get_int(igrib,'indicatorOfParameter',isec1(6),iret)
130	call grib_check(iret,gribFunction,gribErrorMsg)
131	call grib_get_int(igrib,'level',isec1(8),iret)
132	call grib_check(iret,gribFunction,gribErrorMsg)
133
134	!change code for etadot to code for omega
135	if (isec1(6).eq.77) then
136	isec1(6)=135
137	endif
138
139	conversion_factor=1.
140
141
142	else
143
144	!print*,'GRiB Edition 2'
145	!read the grib2 identifiers
146	call grib_get_int(igrib,'discipline',discipl,iret)
147	call grib_check(iret,gribFunction,gribErrorMsg)
148	call grib_get_int(igrib,'parameterCategory',parCat,iret)
149	call grib_check(iret,gribFunction,gribErrorMsg)
150	call grib_get_int(igrib,'parameterNumber',parNum,iret)
151	call grib_check(iret,gribFunction,gribErrorMsg)
152	call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
153	call grib_check(iret,gribFunction,gribErrorMsg)
154	call grib_get_int(igrib,'level',valSurf,iret)
155	call grib_check(iret,gribFunction,gribErrorMsg)
156	call grib_get_int(igrib,'paramId',parId,iret) !added by mc to make it consisitent with new readwind.f90
157	call grib_check(iret,gribFunction,gribErrorMsg) !added by mc to make it consisitent with new readwind.f90
158
159	!print*,discipl,parCat,parNum,typSurf,valSurf
160
161	!convert to grib1 identifiers
162	isec1(6)=-1
163	isec1(7)=-1
164	isec1(8)=-1
165	isec1(8)=valSurf ! level
166	conversion_factor=1.
167	if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
168	isec1(6)=130 ! indicatorOfParameter
169	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
170	isec1(6)=131 ! indicatorOfParameter
171	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
172	isec1(6)=132 ! indicatorOfParameter
173	elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
174	isec1(6)=133 ! indicatorOfParameter
175	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
176	isec1(6)=134 ! indicatorOfParameter
177	elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot !
178	isec1(6)=135 ! indicatorOfParameter
179	elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot !added by mc to make it consisitent with new readwind.f90
180	isec1(6)=135 ! indicatorOfParameter !added by mc to make it consisitent with new readwind.f90
181	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
182	isec1(6)=151 ! indicatorOfParameter
183	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
184	isec1(6)=165 ! indicatorOfParameter
185	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
186	isec1(6)=166 ! indicatorOfParameter
187	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
188	isec1(6)=167 ! indicatorOfParameter
189	elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
190	isec1(6)=168 ! indicatorOfParameter
191	elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
192	isec1(6)=141 ! indicatorOfParameter
193	conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
194	elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC !added by mc to make it consisitent with new readwind.f90
195	isec1(6)=164 ! indicatorOfParameter
196	elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP !added by mc to make it consisitent with new readwind.f90
197	isec1(6)=142 ! indicatorOfParameter
198	elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
199	isec1(6)=143 ! indicatorOfParameter
200	conversion_factor=1000. !added by mc to make it consisitent with new readwind.f90
201	elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
202	isec1(6)=146 ! indicatorOfParameter
203	elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
204	isec1(6)=176 ! indicatorOfParameter
205	elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS !added by mc to make it consisitent with new readwind.f90
206	isec1(6)=180 ! indicatorOfParameter
207	elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS !added by mc to make it consisitent with new readwind.f90
208	isec1(6)=181 ! indicatorOfParameter
209	elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
210	isec1(6)=129 ! indicatorOfParameter
211	elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO !added by mc to make it consisitent with new readwind.f90
212	isec1(6)=160 ! indicatorOfParameter
213	elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
214	(typSurf.eq.1)) then ! LSM
215	isec1(6)=172 ! indicatorOfParameter
216	else
217	print,'**ERROR: undefined GRiB2 message found!',discipl, &
218	parCat,parNum,typSurf
219	endif
220	if(parId .ne. isec1(6) .and. parId .ne. 77) then !added by mc to make it consisitent with new readwind.f90
221	write(,) 'parId',parId, 'isec1(6)',isec1(6) !
222	! stop
223	endif
224
225	endif
226
227	!HSO get the size and data of the values array
228	if (isec1(6).ne.-1) then
229	call grib_get_real4_array(igrib,'values',zsec4,iret)
230	call grib_check(iret,gribFunction,gribErrorMsg)
231	endif
232
233	!HSO get the required fields from section 2 in a gribex compatible manner
234	if(ifield.eq.1) then
235	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
236	isec2(2),iret)
237	call grib_check(iret,gribFunction,gribErrorMsg)
238	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
239	isec2(3),iret)
240	call grib_check(iret,gribFunction,gribErrorMsg)
241	call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
242	isec2(12))
243	call grib_check(iret,gribFunction,gribErrorMsg)
244	! CHECK GRID SPECIFICATIONS
245	if(isec2(2).ne.nxn(l)) stop &
246	'READWIND: NX NOT CONSISTENT FOR A NESTING LEVEL'
247	if(isec2(3).ne.nyn(l)) stop &
248	'READWIND: NY NOT CONSISTENT FOR A NESTING LEVEL'
249	if(isec2(12)/2-1.ne.nlev_ec) stop 'READWIND: VERTICAL DISCRET&
250	IZATION NOT CONSISTENT FOR A NESTING LEVEL'
251	endif ! ifield
252
253	!HSO get the second part of the grid dimensions only from GRiB1 messages
254	if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then ! !added by mc to make it consisitent with new readwind.f90
255	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
256	xauxin,iret)
257	call grib_check(iret,gribFunction,gribErrorMsg)
258	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
259	yauxin,iret)
260	call grib_check(iret,gribFunction,gribErrorMsg)
261	if (xauxin.gt.180.) xauxin=xauxin-360.0
262	if (xauxin.lt.-180.) xauxin=xauxin+360.0
263
264	xaux=xauxin
265	yaux=yauxin
266	if (abs(xaux-xlon0n(l)).gt.eps) &
267	stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT FOR A NESTING LEVEL'
268	if (abs(yaux-ylat0n(l)).gt.eps) &
269	stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT FOR A NESTING LEVEL'
270	gotGrid=1
271	endif
272
273	do j=0,nyn(l)-1
274	do i=0,nxn(l)-1
275	k=isec1(8)
276	if(isec1(6).eq.130) tthn(i,j,nlev_ec-k+2,n,l)= &!! TEMPERATURE
277	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
278	if(isec1(6).eq.131) uuhn(i,j,nlev_ec-k+2,l)= &!! U VELOCITY
279	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
280	if(isec1(6).eq.132) vvhn(i,j,nlev_ec-k+2,l)= &!! V VELOCITY
281	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
282	if(isec1(6).eq.133) then !! SPEC. HUMIDITY
283	qvhn(i,j,nlev_ec-k+2,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
284	if (qvhn(i,j,nlev_ec-k+2,n,l) .lt. 0.) &
285	qvhn(i,j,nlev_ec-k+2,n,l) = 0.
286	! this is necessary because the gridded data may contain
287	! spurious negative values
288	endif
289	if(isec1(6).eq.134) psn(i,j,1,n,l)= &!! SURF. PRESS.
290	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
291	if(isec1(6).eq.135) wwhn(i,j,nlev_ec-k+1,l)= &!! W VELOCITY
292	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
293	if(isec1(6).eq.141) sdn(i,j,1,n,l)= &!! SNOW DEPTH
294	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90!
295	if(isec1(6).eq.151) msln(i,j,1,n,l)= &!! SEA LEVEL PRESS.
296	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
297	if(isec1(6).eq.164) tccn(i,j,1,n,l)= &!! CLOUD COVER
298	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
299	if(isec1(6).eq.165) u10n(i,j,1,n,l)= &!! 10 M U VELOCITY
300	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
301	if(isec1(6).eq.166) v10n(i,j,1,n,l)= &!! 10 M V VELOCITY
302	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
303	if(isec1(6).eq.167) tt2n(i,j,1,n,l)= &!! 2 M TEMPERATURE
304	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
305	if(isec1(6).eq.168) td2n(i,j,1,n,l)= &!! 2 M DEW POINT
306	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
307	if(isec1(6).eq.142) then !! LARGE SCALE PREC.
308	lsprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
309	if (lsprecn(i,j,1,n,l).lt.0.) lsprecn(i,j,1,n,l)=0.
310	endif
311	if(isec1(6).eq.143) then !! CONVECTIVE PREC.
312	convprecn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/conversion_factor !added by mc to make it consisitent with new readwind.f90
313	if (convprecn(i,j,1,n,l).lt.0.) convprecn(i,j,1,n,l)=0.
314	endif
315	if(isec1(6).eq.146) sshfn(i,j,1,n,l)= &!! SENS. HEAT FLUX
316	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
317	if((isec1(6).eq.146).and. &
318	(zsec4(nxn(l)*(nyn(l)-j-1)+i+1).ne.0.)) hflswitch=.true. ! Heat flux available
319	if(isec1(6).eq.176) then !! SOLAR RADIATION
320	ssrn(i,j,1,n,l)=zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
321	if (ssrn(i,j,1,n,l).lt.0.) ssrn(i,j,1,n,l)=0.
322	endif
323	if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
324	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
325	if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
326	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
327	if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
328	(zsec4(nxn(l)*(nyn(l)-j-1)+i+1).ne.0.)) strswitch=.true. ! stress available
329	if(isec1(6).eq.129) oron(i,j,l)= &!! ECMWF OROGRAPHY
330	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)/ga
331	if(isec1(6).eq.160) excessoron(i,j,l)= &!! STANDARD DEVIATION OF OROGRAPHY
332	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
333	if(isec1(6).eq.172) lsmn(i,j,l)= &!! ECMWF LAND SEA MASK
334	zsec4(nxn(l)*(nyn(l)-j-1)+i+1)
335	if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
336	if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
337
338	end do
339	end do
340
341	call grib_release(igrib)
342	goto 10 !! READ NEXT LEVEL OR PARAMETER
343	!
344	! CLOSING OF INPUT DATA FILE
345	!
346	50 call grib_close_file(ifile)
347
348	!error message if no fields found with correct first longitude in it
349	if (gotGrid.eq.0) then
350	print,'**ERROR: input file needs to contain GRiB1 formatted'// &
351	'messages'
352	stop
353	endif
354
355	if(levdiff2.eq.0) then
356	iwmax=nlev_ec+1
357	do i=0,nxn(l)-1
358	do j=0,nyn(l)-1
359	wwhn(i,j,nlev_ec+1,l)=0.
360	end do
361	end do
362	endif
363
364	do i=0,nxn(l)-1
365	do j=0,nyn(l)-1
366	surfstrn(i,j,1,n,l)=sqrt(ewss(i,j)2+nsss(i,j)2)
367	end do
368	end do
369
370	if ((.not.hflswitch).or.(.not.strswitch)) then
371	write(,) 'WARNING: No flux data contained in GRIB file ', &
372	wfnamen(l,indj)
373
374	! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
375	! As ECMWF has increased the model resolution, such that now the first model
376	! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
377	! (3rd model level in FLEXPART) for the profile method
378	!***************************************************************************
379
380	do i=0,nxn(l)-1
381	do j=0,nyn(l)-1
382	plev1=akz(3)+bkz(3)*psn(i,j,1,n,l)
383	pmean=0.5*(psn(i,j,1,n,l)+plev1)
384	tv=tthn(i,j,3,n,l)(1.+0.61qvhn(i,j,3,n,l))
385	fu=-r_air*tv/ga/pmean
386	hlev1=fu*(plev1-psn(i,j,1,n,l)) ! HEIGTH OF FIRST MODEL LAYER
387	ff10m= sqrt(u10n(i,j,1,n,l)2+v10n(i,j,1,n,l)2)
388	fflev1=sqrt(uuhn(i,j,3,l)2+vvhn(i,j,3,l)2)
389	call pbl_profile(psn(i,j,1,n,l),td2n(i,j,1,n,l),hlev1, &
390	tt2n(i,j,1,n,l),tthn(i,j,3,n,l),ff10m,fflev1, &
391	surfstrn(i,j,1,n,l),sshfn(i,j,1,n,l))
392	if(sshfn(i,j,1,n,l).gt.200.) sshfn(i,j,1,n,l)=200.
393	if(sshfn(i,j,1,n,l).lt.-400.) sshfn(i,j,1,n,l)=-400.
394	end do
395	end do
396	endif
397
398
399	! Assign 10 m wind to model level at eta=1.0 to have one additional model
400	! level at the ground
401	! Specific humidity is taken the same as at one level above
402	! Temperature is taken as 2 m temperature
403	!**************************************************************************
404
405	do i=0,nxn(l)-1
406	do j=0,nyn(l)-1
407	uuhn(i,j,1,l)=u10n(i,j,1,n,l)
408	vvhn(i,j,1,l)=v10n(i,j,1,n,l)
409	qvhn(i,j,1,n,l)=qvhn(i,j,2,n,l)
410	tthn(i,j,1,n,l)=tt2n(i,j,1,n,l)
411	end do
412	end do
413
414	if(iumax.ne.nuvz-1) stop &
415	'READWIND: NUVZ NOT CONSISTENT FOR A NESTING LEVEL'
416	if(iwmax.ne.nwz) stop &
417	'READWIND: NWZ NOT CONSISTENT FOR A NESTING LEVEL'
418
419	end do
420
421	return
422	888 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
423	write(,) ' #### ',wfnamen(l,indj),' FOR NESTING LEVEL #### '
424	write(,) ' #### ',l,' IS NOT GRIB FORMAT !!! #### '
425	stop 'Execution terminated'
426
427
428	999 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
429	write(,) ' #### ',wfnamen(l,indj),' #### '
430	write(,) ' #### CANNOT BE OPENED FOR NESTING LEVEL ',l,'####'
431
432	end subroutine readwind_nests

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format