Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

readwind_nests.f90 @ 8

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

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

Download in other formats:

Original Format