Context Navigation

readwind_nests.F90 @ b398fb6

fp9.3.1-20161214-nc4

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

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

Download in other formats:

Original Format