Context Navigation

readwind.F90 @ 496c607

FPv9.3.1FPv9.3.1b_testingFPv9.3.2fp9.3.1-20161214-nc4grib2nc4_repair

Last change on this file since 496c607 was 496c607, checked in by Don Morton <Don.Morton@…>, 8 years ago

Initial commit of FPv9.3.1

Currently, this is a clone of snapshot FPv9.3.0

Property mode set to 100644

File size: 25.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
22	subroutine readwind_ecmwf(indj,n,uuh,vvh,wwh)
23
24	!**********************************************************************
25	! *
26	! TRAJECTORY MODEL SUBROUTINE READWIND *
27	! *
28	!**********************************************************************
29	! *
30	! AUTHOR: G. WOTAWA *
31	! DATE: 1997-08-05 *
32	! LAST UPDATE: 2000-10-17, Andreas Stohl *
33	! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
34	! ECMWF grib_api *
35	! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
36	! ECMWF grib_api *
37	! *
38	!**********************************************************************
39	! Changes, Bernd C. Krueger, Feb. 2001:
40	! Variables tth and qvh (on eta coordinates) in common block
41	!**********************************************************************
42	! *
43	! DESCRIPTION: *
44	! *
45	! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE *
46	! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE *
47	! *
48	!**********************************************************************
49	! Changes Arnold, D. and Morton, D. (2015): *
50	! -- description of local and common variables *
51	!**********************************************************************
52
53	use GRIB_API
54	use par_mod
55	use com_mod
56	use class_vtable
57
58	implicit none
59
60	!***********************************************************************
61	! Subroutine Parameters: *
62	! input *
63	! indj indicates number of the wind field to be read in *
64	! n temporal index for meteorological fields (1 to 3)*
65	! uuh,vvh, wwh wind components in ECMWF model levels *
66	integer :: indj,n
67	real(kind=4) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
68	real(kind=4) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
69	real(kind=4) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
70	!***********************************************************************
71
72	!***********************************************************************
73	! Local variables *
74	! *
75	! HSO variables for grib_api: *
76	! ifile grib file to be opened and read in *
77	! iret is a return code for successful or not open *
78	! igrib grib edition number (whether GRIB1 or GRIB2) *
79	! gribVer where info on igrib is kept, 1 for GRIB1 and *
80	! 2 for GRIB2 *
81	! parCat parameter category ( = number) , how FLEXPART *
82	! identifies fields *
83	! parNum parameter number by product discipline and *
84	! parameter category *
85	! typSurf type of first fixed surface *
86	! valSurf level *
87	! discipl discipline of processed data contained in a *
88	! GRIB message *
89	integer :: ifile
90	integer :: iret
91	integer :: igrib
92	integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl
93	#if defined CTBTO
94	integer :: parId
95	#endif
96	integer :: gotGrid
97
98
99	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
100	!!!! Vtable related variables
101
102	! Paths to Vtables (current implementation will assume they are in the cwd
103	! This is a crappy place for these parameters. Need to move them.
104	character(LEN=255), parameter :: VTABLE_ECMWF_GRIB1_PATH = &
105	"Vtable_ecmwf_grib1", &
106	VTABLE_ECMWF_GRIB2_PATH = &
107	"Vtable_ecmwf_grib2", &
108	VTABLE_ECMWF_GRIB1_2_PATH = &
109	"Vtable_ecmwf_grib1_2"
110
111	integer :: gribfile_type
112	integer :: current_grib_level ! This "was" isec1(8) in previous version
113	character(len=255) :: gribfile_name
114	character(len=255) :: vtable_path
115	character(len=15) :: fpname
116
117	type(Vtable) :: my_vtable ! unallocated
118	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
119
120
121
122	!
123	! Variables and arrays for grib decoding: *
124	! dimension of isec2 at least (22+n), where n is the number of *
125	! parallels or meridians in a quasi-regular (reduced) Gaussian or *
126	! lat/long grid *
127	! dimension of zsec2 at least (10+nn), where nn is the number of *
128	! vertical coordinate parameters *
129	! *
130	! isec1 grib definition section (version, center, ... *
131	! isec2 grid description section *
132	! zsec4 the binary data section *
133	! xaux
134	! yaux
135	! xauxin
136	! yauxin
137	! xaux0 auxiliary variable for xlon0 *
138	! yaux0 auxiliary variable for xlat0 *
139	! nsss NS surface stress *
140	! ewss EW surface stress *
141	! plev1 pressure of the first model layer *
142	! pmean mean pressure between ?? *
143	! tv virtual temperature *
144	! fu for conversion from pressure to meters *
145	! hlev1 height of the first model layer *
146	! ff10m wind speed at 10 m *
147	! fflev1 wind speed at the first model layere *
148	! hflswitch logical variable to check existence of flux data *
149	! strswitch logical variable to check existence of stress *
150	! data *
151
152	integer :: isec1(56),isec2(22+nxmax+nymax)
153	real(kind=4) :: zsec4(jpunp)
154	real(kind=4) :: xaux,yaux,xaux0,yaux0
155	#if defined CTBTO
156	real(kind=4) :: conversion_factor
157	#endif
158	real(kind=8) :: xauxin,yauxin
159	real(kind=4) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
160	real :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1
161	logical :: hflswitch,strswitch
162
163	! Other variables:
164	! i,j,k loop control indices in each direction *
165	! levdiff2 number of model levels - number of model levels *
166	! for the staggered vertical wind +1 *
167	! ifield index to control field read in *
168	! iumax
169	! iwmax
170	integer :: i,j,k,levdiff2,ifield,iumax,iwmax
171	!***********************************************************************
172
173	!***********************************************************************
174	! Local constants *
175	real,parameter :: eps=1.e-4
176	!HSO grib api error messages:
177	character(len=24) :: gribErrorMsg = 'Error reading grib file'
178	character(len=20) :: gribFunction = 'readwind'
179	!***********************************************************************
180
181	!***********************************************************************
182	! Global variables *
183	! from par_mod and com_mod *
184	! wfname File name of data to be read in *
185	! nx,ny,nuvz,nwz expected field dimensions *
186	! nxfield same as nx but for limited area fields *
187	! nxmin1,nymin1 nx-1, ny-1, respectively *
188	! nxmax,nymax maximum dimension of wind fields in x and y *
189	! direction, respectively *
190	! nuvzmax,nwzmax maximum dimension of (u,v) and (w) wind fields in z
191	! direction (for fields on eta levels) *
192	! nlev_ec number of vertical levels ecmwf model *
193	! u10,v10 wind fields at 10 m *
194	! tt2,td2 temperature and dew point temperature at 2 m *
195	! tth,qvh temperature and specific humidity on original *
196	! eta levels *
197	! ps surface pressure *
198	! sd snow depth (but not used afterwards!) *
199	! msl mean sea level pressure *
200	! ttc total cloud cover *
201	! lsprec large scale precipitation *
202	! convprec convective precipitation *
203	! sshf sensible heat flux *
204	! ssr solar radiation *
205	! surfstr surface stress *
206	! oro orography *
207	! excessoro standard deviation of orography *
208	! lsm land sea mask *
209	! r_air individual gas constant for dry air [J/kg/K] *
210	! ga gravity acceleration of earth [m/s*2]
211	!***********************************************************************
212
213	!-----------------------------------------------------------------------------
214
215
216	hflswitch=.false.
217	strswitch=.false.
218	levdiff2=nlev_ec-nwz+1
219	iumax=0
220	iwmax=0
221
222	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
223	!!!!!!!!!!!!!!!!!!! VTABLE code
224	!!!!!!! Vtable choice
225	gribfile_name = path(3)(1:length(3))//trim(wfname(indj))
226	print *, 'gribfile_name: ', gribfile_name
227
228	gribfile_type = vtable_detect_gribfile_type( gribfile_name )
229
230	print *, 'gribfile_type: ', gribfile_type
231
232	if (gribfile_type .eq. VTABLE_GRIBFILE_TYPE_ECMWF_GRIB1) then
233	vtable_path = VTABLE_ECMWF_GRIB1_PATH
234	else if (gribfile_type .eq. VTABLE_GRIBFILE_TYPE_ECMWF_GRIB2) then
235	vtable_path = VTABLE_ECMWF_GRIB2_PATH
236	else if (gribfile_type .eq. VTABLE_GRIBFILE_TYPE_ECMWF_GRIB1_2) then
237	vtable_path = VTABLE_ECMWF_GRIB1_2_PATH
238	else
239	print *, 'Unsupported gribfile_type: ', gribfile_type
240	stop
241	endif
242
243
244	! Load the Vtable into 'my_vtable'
245	print *, 'Loading Vtable: ', vtable_path
246	call vtable_load_by_name(vtable_path, my_vtable)
247	print *, 'Vtable Initialized: ', my_vtable%initialized
248	print *, 'Vtable num_entries: ', my_vtable%num_entries
249	!!!!!!!!!!!!!!!!!!! VTABLE code
250	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
251
252
253
254
255
256	!
257	! OPENING OF DATA FILE (GRIB CODE)
258	!
259	5 call grib_open_file(ifile,path(3)(1:length(3)) &
260	//trim(wfname(indj)),'r',iret)
261	if (iret.ne.GRIB_SUCCESS) then
262	goto 888 ! ERROR DETECTED
263	endif
264	!turn on support for multi fields messages */
265	!call grib_multi_support_on
266
267	gotGrid=0
268	ifield=0
269	10 ifield=ifield+1
270	!
271	! GET NEXT FIELDS
272	!
273	call grib_new_from_file(ifile,igrib,iret)
274	if (iret.eq.GRIB_END_OF_FILE) then
275	goto 50 ! EOF DETECTED
276	elseif (iret.ne.GRIB_SUCCESS) then
277	goto 888 ! ERROR DETECTED
278	endif
279
280	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
281	!!!!!!!!!!!!!!!!!!! VTABLE code
282	! Get the fpname
283	fpname = vtable_get_fpname(igrib, my_vtable)
284	!print *, 'fpname: ', trim(fpname)
285
286
287	!!!!!!!!!!!!!!!!!!! VTABLE code
288	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
289
290	#if defined CTBTO
291	conversion_factor=1.0
292	#endif
293
294	!first see if we read GRIB1 or GRIB2
295	call grib_get_int(igrib,'editionNumber',gribVer,iret)
296	call grib_check(iret,gribFunction,gribErrorMsg)
297
298	!!!!!!! DJM - candidate for consolidation
299	if (gribVer.eq.1) then ! GRIB Edition 1
300
301	call grib_get_int(igrib,'level', current_grib_level,iret)
302	call grib_check(iret,gribFunction,gribErrorMsg)
303
304
305	else ! GRIB Edition 2
306
307	#if defined CTBTO
308	!print *, 'fpname: ', fpname
309	call grib2check(igrib, fpname, conversion_factor)
310	#endif
311
312	call grib_get_int(igrib,'level', current_grib_level,iret)
313	call grib_check(iret,gribFunction,gribErrorMsg)
314
315
316	endif ! END IF Grib Version selection
317
318
319
320	!HSO get the size and data of the values array
321	if (trim(fpname) .ne. 'None') then
322	call grib_get_real4_array(igrib,'values',zsec4,iret)
323	#if defined CTBTO
324	zsec4=zsec4/conversion_factor
325	#endif
326	call grib_check(iret,gribFunction,gribErrorMsg)
327	endif
328
329	!HSO get the required fields from section 2 in a gribex compatible manner
330	if (ifield.eq.1) then
331	call grib_get_int(igrib,'numberOfPointsAlongAParallel', &
332	isec2(2),iret)
333	call grib_check(iret,gribFunction,gribErrorMsg)
334	call grib_get_int(igrib,'numberOfPointsAlongAMeridian', &
335	isec2(3),iret)
336	call grib_check(iret,gribFunction,gribErrorMsg)
337	call grib_get_int(igrib,'numberOfVerticalCoordinateValues', &
338	isec2(12))
339	call grib_check(iret,gribFunction,gribErrorMsg)
340	! CHECK GRID SPECIFICATIONS
341	if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
342	if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
343	if(isec2(12)/2-1.ne.nlev_ec) &
344	stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
345	endif ! ifield
346
347	!HSO get the second part of the grid dimensions only from GRiB1 messages
348	#if defined CTBTO
349	if (gotGrid.eq.0) then
350	#else
351	!!! DJM 2016-05-29 -- this check for gribVer doesn't make sense, and causes
352	!!! failure with pure GRIB2, so I'm commenting out
353	!if ((gribVer.eq.1).and.(gotGrid.eq.0)) then
354	if (gotGrid.eq.0) then
355	#endif
356	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
357	xauxin,iret)
358	call grib_check(iret,gribFunction,gribErrorMsg)
359	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
360	yauxin,iret)
361	call grib_check(iret,gribFunction,gribErrorMsg)
362	xaux=xauxin+real(nxshift)*dx
363	yaux=yauxin
364	xaux0=xlon0
365	yaux0=ylat0
366	if(xaux.lt.0.) xaux=xaux+360.
367	if(yaux.lt.0.) yaux=yaux+360.
368	if(xaux0.lt.0.) xaux0=xaux0+360.
369	if(yaux0.lt.0.) yaux0=yaux0+360.
370	if(abs(xaux-xaux0).gt.eps) &
371	stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
372	if(abs(yaux-yaux0).gt.eps) &
373	stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
374	gotGrid=1
375	endif ! gotGrid
376
377
378	k = current_grib_level
379
380	!! TEMPERATURE
381	if(trim(fpname) .eq. 'TT') then
382	do j=0,nymin1
383	do i=0,nxfield-1
384	tth(i,j,nlev_ec-k+2,n) = zsec4(nxfield*(ny-j-1)+i+1)
385	enddo
386	enddo
387	endif
388
389
390	!! U VELOCITY
391	if(trim(fpname) .eq. 'UU') then
392	iumax=max(iumax,nlev_ec-k+1)
393	do j=0,nymin1
394	do i=0,nxfield-1
395	uuh(i,j,nlev_ec-k+2) = zsec4(nxfield*(ny-j-1)+i+1)
396	enddo
397	enddo
398	endif
399
400
401	!! V VELOCITY
402	if(trim(fpname) .eq. 'VV') then
403	do j=0,nymin1
404	do i=0,nxfield-1
405	vvh(i,j,nlev_ec-k+2) = zsec4(nxfield*(ny-j-1)+i+1)
406	enddo
407	enddo
408	endif
409
410	!! W VELOCITY
411	if(trim(fpname) .eq. 'ETADOT') then
412	iwmax=max(iwmax,nlev_ec-k+1)
413	do j=0,nymin1
414	do i=0,nxfield-1
415	wwh(i,j,nlev_ec-k+1) = zsec4(nxfield*(ny-j-1)+i+1)
416	enddo
417	enddo
418	endif
419
420	!! SPEC HUMIDITY
421	if(trim(fpname) .eq. 'QV') then
422	do j=0,nymin1
423	do i=0,nxfield-1
424	qvh(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
425	if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) qvh(i,j,nlev_ec-k+2,n) = 0.
426	! this is necessary because the gridded data may contain
427	! spurious negative values
428	enddo
429	enddo
430	endif
431
432	!! SURF. PRESS.
433	if(trim(fpname) .eq. 'PS') then
434	do j=0,nymin1
435	do i=0,nxfield-1
436	ps(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
437	enddo
438	enddo
439	endif
440
441	!! SNOW DEPTH
442	if(trim(fpname) .eq. 'SD') then
443	do j=0,nymin1
444	do i=0,nxfield-1
445	sd(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
446	enddo
447	enddo
448	endif
449
450	!! SEA LEVEL PRESS.
451	if(trim(fpname) .eq. 'MSL') then
452	do j=0,nymin1
453	do i=0,nxfield-1
454	msl(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
455	enddo
456	enddo
457	endif
458
459	!! CLOUD COVER
460	if(trim(fpname) .eq. 'TCC') then
461	do j=0,nymin1
462	do i=0,nxfield-1
463	tcc(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
464	enddo
465	enddo
466	endif
467
468	!! 10 M U VELOCITY
469	if(trim(fpname) .eq. 'U10') then
470	do j=0,nymin1
471	do i=0,nxfield-1
472	u10(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
473	enddo
474	enddo
475	endif
476
477	!! 10 M V VELOCITY
478	if(trim(fpname) .eq. 'V10') then
479	do j=0,nymin1
480	do i=0,nxfield-1
481	v10(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
482	enddo
483	enddo
484	endif
485
486	!! 2 M TEMPERATURE
487	if(trim(fpname) .eq. 'T2') then
488	do j=0,nymin1
489	do i=0,nxfield-1
490	tt2(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
491	enddo
492	enddo
493	endif
494
495	!! 2 M DEW POINT
496	if(trim(fpname) .eq. 'TD2') then
497	do j=0,nymin1
498	do i=0,nxfield-1
499	td2(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
500	enddo
501	enddo
502	endif
503
504	!! LARGE SCALE PREC.
505	if(trim(fpname) .eq. 'LSPREC') then
506	do j=0,nymin1
507	do i=0,nxfield-1
508	lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
509	if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
510	enddo
511	enddo
512	endif
513
514	!! CONVECTIVE PREC.
515	if(trim(fpname) .eq. 'CONVPREC') then
516	do j=0,nymin1
517	do i=0,nxfield-1
518	convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
519	if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
520	enddo
521	enddo
522	endif
523
524	!! SENS. HEAT FLUX
525	if(trim(fpname) .eq. 'SHF') then
526	do j=0,nymin1
527	do i=0,nxfield-1
528	sshf(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
529	if (zsec4(nxfield*(ny-j-1)+i+1).ne.0.) hflswitch=.true. ! Heat flux available
530	enddo
531	enddo
532	endif
533
534	!! SOLAR RADIATION
535	if(trim(fpname) .eq. 'SR') then
536	do j=0,nymin1
537	do i=0,nxfield-1
538	ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
539	if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
540	enddo
541	enddo
542	endif
543
544	!! EW SURFACE STRESS
545	if(trim(fpname) .eq. 'EWSS') then
546	do j=0,nymin1
547	do i=0,nxfield-1
548	ewss(i,j)=zsec4(nxfield*(ny-j-1)+i+1)
549	if (zsec4(nxfield*(ny-j-1)+i+1).ne.0.) strswitch=.true. ! stress available
550	enddo
551	enddo
552	endif
553
554	!! NS SURFACE STRESS
555	if(trim(fpname) .eq. 'NSSS') then
556	do j=0,nymin1
557	do i=0,nxfield-1
558	nsss(i,j)=zsec4(nxfield*(ny-j-1)+i+1)
559	if (zsec4(nxfield*(ny-j-1)+i+1).ne.0.) strswitch=.true. ! stress available
560	enddo
561	enddo
562	endif
563
564	!! ECMWF OROGRAPHY
565	if(trim(fpname) .eq. 'ORO') then
566	do j=0,nymin1
567	do i=0,nxfield-1
568	oro(i,j)=zsec4(nxfield*(ny-j-1)+i+1)/ga
569	enddo
570	enddo
571	endif
572
573	!! ECMWF OROGRAPHY
574	if(trim(fpname) .eq. 'EXCESSORO') then
575	do j=0,nymin1
576	do i=0,nxfield-1
577	excessoro(i,j)=zsec4(nxfield*(ny-j-1)+i+1)
578	enddo
579	enddo
580	endif
581
582	!! ECMWF LAND SEA MASK
583	if(trim(fpname) .eq. 'LSM') then
584	do j=0,nymin1
585	do i=0,nxfield-1
586	lsm(i,j)=zsec4(nxfield*(ny-j-1)+i+1)
587	enddo
588	enddo
589	endif
590
591
592	call grib_release(igrib)
593	goto 10 !! READ NEXT LEVEL OR PARAMETER
594	!
595	! CLOSING OF INPUT DATA FILE
596	!
597
598	50 call grib_close_file(ifile)
599
600	!error message if no fields found with correct first longitude in it
601	if (gotGrid.eq.0) then
602	print,'**ERROR: no fields found with correct first longitude'// &
603	'messages'
604	stop
605	endif
606
607	if(levdiff2.eq.0) then
608	iwmax=nlev_ec+1
609	do i=0,nxmin1
610	do j=0,nymin1
611	wwh(i,j,nlev_ec+1)=0.
612	end do
613	end do
614	endif
615
616	! For global fields, assign the leftmost data column also to the rightmost
617	! data column; if required, shift whole grid by nxshift grid points
618	!*************************************************************************
619
620	if (xglobal) then
621	call shift_field_0(ewss,nxfield,ny)
622	call shift_field_0(nsss,nxfield,ny)
623	call shift_field_0(oro,nxfield,ny)
624	call shift_field_0(excessoro,nxfield,ny)
625	call shift_field_0(lsm,nxfield,ny)
626	call shift_field(ps,nxfield,ny,1,1,2,n)
627	call shift_field(sd,nxfield,ny,1,1,2,n)
628	call shift_field(msl,nxfield,ny,1,1,2,n)
629	call shift_field(tcc,nxfield,ny,1,1,2,n)
630	call shift_field(u10,nxfield,ny,1,1,2,n)
631	call shift_field(v10,nxfield,ny,1,1,2,n)
632	call shift_field(tt2,nxfield,ny,1,1,2,n)
633	call shift_field(td2,nxfield,ny,1,1,2,n)
634	call shift_field(lsprec,nxfield,ny,1,1,2,n)
635	call shift_field(convprec,nxfield,ny,1,1,2,n)
636	call shift_field(sshf,nxfield,ny,1,1,2,n)
637	call shift_field(ssr,nxfield,ny,1,1,2,n)
638	call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
639	call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
640	call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
641	call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
642	call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
643	endif
644
645	do i=0,nxmin1
646	do j=0,nymin1
647	surfstr(i,j,1,n)=sqrt(ewss(i,j)2+nsss(i,j)2)
648	end do
649	end do
650
651	if ((.not.hflswitch).or.(.not.strswitch)) then
652	write(,) 'WARNING: No flux data contained in GRIB file ', &
653	wfname(indj)
654
655	! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
656	! As ECMWF has increased the model resolution, such that now the first model
657	! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
658	! (3rd model level in FLEXPART) for the profile method
659	!***************************************************************************
660
661	do i=0,nxmin1
662	do j=0,nymin1
663	plev1=akz(3)+bkz(3)*ps(i,j,1,n)
664	pmean=0.5*(ps(i,j,1,n)+plev1)
665	tv=tth(i,j,3,n)(1.+0.61qvh(i,j,3,n))
666	fu=-r_air*tv/ga/pmean
667	hlev1=fu*(plev1-ps(i,j,1,n)) ! HEIGTH OF FIRST MODEL LAYER
668	ff10m= sqrt(u10(i,j,1,n)2+v10(i,j,1,n)2)
669	fflev1=sqrt(uuh(i,j,3)2+vvh(i,j,3)2)
670	call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
671	tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
672	surfstr(i,j,1,n),sshf(i,j,1,n))
673	if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
674	if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
675	end do
676	end do
677	endif
678
679
680	! Assign 10 m wind to model level at eta=1.0 to have one additional model
681	! level at the ground
682	! Specific humidity is taken the same as at one level above
683	! Temperature is taken as 2 m temperature
684	!**************************************************************************
685
686	do i=0,nxmin1
687	do j=0,nymin1
688	uuh(i,j,1)=u10(i,j,1,n)
689	vvh(i,j,1)=v10(i,j,1,n)
690	qvh(i,j,1,n)=qvh(i,j,2,n)
691	tth(i,j,1,n)=tt2(i,j,1,n)
692	end do
693	end do
694
695	if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT'
696	if(iwmax.ne.nwz) stop 'READWIND: NWZ NOT CONSISTENT'
697
698	return
699	888 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
700	write(,) ' #### ',wfname(indj),' #### '
701	write(,) ' #### IS NOT GRIB FORMAT !!! #### '
702	stop 'Execution terminated'
703	999 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
704	write(,) ' #### ',wfname(indj),' #### '
705	write(,) ' #### CANNOT BE OPENED !!! #### '
706	stop 'Execution terminated'
707
708	end subroutine readwind_ecmwf

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

Download in other formats:

Original Format