Context Navigation

readwind_ecmwf_mpi.f90

bugfixes+enhancements

Last change on this file was dba4221, checked in by Pirmin Kaufmann <pirmin.kaufmann@…>, 18 months ago

Bugfixes:

options/SPECIES/SPECIES_009: corrected wrong number format, replaced comma by decimal point
options/SPECIES/SPECIES_028: corrected wrong number format, moved sign of exponent to after the E
options/SPECIES/specoverview.f90: added namelist parameters that appear in SPECIES files but were missing here
src/FLEXPART.f90: replaced compiler-specific command line argument routines by standard Fortran intrinsic routines
src/FLEXPART_MPI.f90: ditto
src/gridcheck_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
src/gridcheck_nests.f90: ditto
src/readwind_ecmwf.f90: corrected handling of vertical levels when input files do not contain uppermost layers
readwind_ecmwf_mpi.f90: ditto

Code enhancements:

options/OUTGRID: added comments describing contents
options/SPECIES/SPECIES_*: aligned comments
options/SPECIES/specoverview.f90: removed commented lines, rectified lines indenting
src/FLEXPART.f90: rectified lines indenting, updated date in version string
src/FLEXPART_MPI.f90: ditto, and realigned code with src/FLEXPART.f90
src/gridcheck_*.f90: added code to write out name of file before it is opened (helps a lot when an input file causes troubles)
src/par_mod.f90: added comment explaining relevance of nuvzmax for GRIB input
src/readreleases.f90: write out warning if too few particles are used to randomize release
src/readspecies.f90: write out name of SPECIES file before it is read
src/readwind_*.f90: write out name of input file before opening it
src/writeheader_txt.f90: removed wrong comment

Property mode set to 100644

File size: 20.7 KB

Line
1	! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2	! SPDX-License-Identifier: GPL-3.0-or-later
3
4	subroutine readwind_ecmwf(indj,n,uuh,vvh,wwh)
5
6	!**********************************************************************
7	! *
8	! TRAJECTORY MODEL SUBROUTINE READWIND *
9	! *
10	!**********************************************************************
11	! *
12	! AUTHOR: G. WOTAWA *
13	! DATE: 1997-08-05 *
14	! LAST UPDATE: 2000-10-17, Andreas Stohl *
15	! CHANGE: 11/01/2008, Harald Sodemann, GRIB1/2 input with *
16	! ECMWF grib_api *
17	! CHANGE: 03/12/2008, Harald Sodemann, update to f90 with *
18	! ECMWF grib_api *
19	! *
20	!**********************************************************************
21	! Changes, Bernd C. Krueger, Feb. 2001:
22	! Variables tth and qvh (on eta coordinates) in common block
23	! eso 2014:
24	! This version for reading windfields in advance by a separate
25	! MPI process.
26	! TODO: can be merged with readwind.f90 if using numwfmem in
27	! shift_field
28	!
29	! Unified ECMWF and GFS builds
30	! Marian Harustak, 12.5.2017
31	! - Renamed from reawind to readwind_ecmwf
32	!
33	!**********************************************************************
34	! *
35	! DESCRIPTION: *
36	! *
37	! READING OF ECMWF METEOROLOGICAL FIELDS FROM INPUT DATA FILES. THE *
38	! INPUT DATA FILES ARE EXPECTED TO BE AVAILABLE IN GRIB CODE *
39	! *
40	! INPUT: *
41	! indj indicates number of the wind field to be read in *
42	! n temporal index for meteorological fields (1 to 3)*
43	! *
44	! IMPORTANT VARIABLES FROM COMMON BLOCK: *
45	! *
46	! wfname File name of data to be read in *
47	! nx,ny,nuvz,nwz expected field dimensions *
48	! nlev_ec number of vertical levels ecmwf model *
49	! uu,vv,ww wind fields *
50	! tt,qv temperature and specific humidity *
51	! ps surface pressure *
52	! *
53	!**********************************************************************
54
55	use GRIB_API
56	use par_mod
57	use com_mod
58	use mpi_mod
59
60	implicit none
61
62	! include 'grib_api.h'
63
64	!HSO parameters for grib_api
65	integer :: ifile
66	integer :: iret
67	integer :: igrib
68	integer :: gribVer,parCat,parNum,typSurf,valSurf,discipl,parId
69	integer :: gotGrid
70	!HSO end
71
72	real(sp), intent(inout) :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
73	real(sp), intent(inout) :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
74	real(sp), intent(inout) :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
75	integer, intent(in) :: indj,n
76	integer :: i,j,k,levdiff2,ifield,iumax,iwmax
77
78	! VARIABLES AND ARRAYS NEEDED FOR GRIB DECODING
79
80	! dimension of isec2 at least (22+n), where n is the number of parallels or
81	! meridians in a quasi-regular (reduced) Gaussian or lat/long grid
82
83	! dimension of zsec2 at least (10+nn), where nn is the number of vertical
84	! coordinate parameters
85
86	integer :: isec1(56),isec2(22+nxmax+nymax)
87	real(sp) :: zsec4(jpunp)
88	real(sp) :: xaux,yaux
89	real(dp) :: xauxin,yauxin
90	real(sp),parameter :: eps=1.e-4
91	real(sp) :: nsss(0:nxmax-1,0:nymax-1),ewss(0:nxmax-1,0:nymax-1)
92	real(sp) :: plev1,pmean,tv,fu,hlev1,ff10m,fflev1,conversion_factor
93
94	logical :: hflswitch,strswitch !,readcloud
95
96	!HSO grib api error messages
97	character(len=24) :: gribErrorMsg = 'Error reading grib file'
98	character(len=20) :: gribFunction = 'readwind'
99
100	! measure time
101	if (mp_measure_time) call mpif_mtime('readwind',0)
102
103	hflswitch=.false.
104	strswitch=.false.
105	!ZHG test the grib fields that have lcwc without using them
106	! readcloud=.false.
107
108	levdiff2=nlev_ec-nwz+1
109	iumax=0
110	iwmax=0
111
112	!
113	! OPENING OF DATA FILE (GRIB CODE)
114	!
115	write(,) 'Reading: '//path(3)(1:length(3)) &
116	//trim(wfname(indj))
117	5 call grib_open_file(ifile,path(3)(1:length(3)) &
118	//trim(wfname(indj)),'r',iret)
119	if (iret.ne.GRIB_SUCCESS) then
120	goto 888 ! ERROR DETECTED
121	endif
122	!turn on support for multi fields messages */
123	!call grib_multi_support_on
124
125	gotGrid=0
126	ifield=0
127	10 ifield=ifield+1
128	!
129	! GET NEXT FIELDS
130	!
131	call grib_new_from_file(ifile,igrib,iret)
132	if (iret.eq.GRIB_END_OF_FILE) then
133	goto 50 ! EOF DETECTED
134	elseif (iret.ne.GRIB_SUCCESS) then
135	goto 888 ! ERROR DETECTED
136	endif
137
138	!first see if we read GRIB1 or GRIB2
139	call grib_get_int(igrib,'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,'indicatorOfParameter',isec1(6),iret)
147	call grib_check(iret,gribFunction,gribErrorMsg)
148	call grib_get_int(igrib,'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	conversion_factor=1.
157
158	else
159
160	!print*,'GRiB Edition 2'
161	!read the grib2 identifiers
162	call grib_get_int(igrib,'discipline',discipl,iret)
163	call grib_check(iret,gribFunction,gribErrorMsg)
164	call grib_get_int(igrib,'parameterCategory',parCat,iret)
165	call grib_check(iret,gribFunction,gribErrorMsg)
166	call grib_get_int(igrib,'parameterNumber',parNum,iret)
167	call grib_check(iret,gribFunction,gribErrorMsg)
168	call grib_get_int(igrib,'typeOfFirstFixedSurface',typSurf,iret)
169	call grib_check(iret,gribFunction,gribErrorMsg)
170	call grib_get_int(igrib,'level',valSurf,iret)
171	call grib_check(iret,gribFunction,gribErrorMsg)
172	call grib_get_int(igrib,'paramId',parId,iret)
173	call grib_check(iret,gribFunction,gribErrorMsg)
174
175	!print*,discipl,parCat,parNum,typSurf,valSurf
176
177	!convert to grib1 identifiers
178	isec1(6)=-1
179	isec1(7)=-1
180	isec1(8)=-1
181	isec1(8)=valSurf ! level
182	conversion_factor=1.
183	if ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! T
184	isec1(6)=130 ! indicatorOfParameter
185	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.105)) then ! U
186	isec1(6)=131 ! indicatorOfParameter
187	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.105)) then ! V
188	isec1(6)=132 ! indicatorOfParameter
189	elseif ((parCat.eq.1).and.(parNum.eq.0).and.(typSurf.eq.105)) then ! Q
190	isec1(6)=133 ! indicatorOfParameter
191	! ESO Cloud water is in a) fields CLWC and CIWC, or b) field qc
192	elseif ((parCat.eq.1).and.(parNum.eq.83).and.(typSurf.eq.105)) then ! clwc
193	isec1(6)=246 ! indicatorOfParameter
194	elseif ((parCat.eq.1).and.(parNum.eq.84).and.(typSurf.eq.105)) then ! ciwc
195	isec1(6)=247 ! indicatorOfParameter
196	! ESO qc(=clwc+ciwc):
197	elseif ((parCat.eq.201).and.(parNum.eq.31).and.(typSurf.eq.105)) then ! qc
198	isec1(6)=201031 ! indicatorOfParameter
199	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.1)) then !SP
200	isec1(6)=134 ! indicatorOfParameter
201	elseif ((parCat.eq.2).and.(parNum.eq.32)) then ! W, actually eta dot
202	isec1(6)=135 ! indicatorOfParameter
203	elseif ((parCat.eq.128).and.(parNum.eq.77)) then ! W, actually eta dot
204	isec1(6)=135 ! indicatorOfParameter
205	elseif ((parCat.eq.3).and.(parNum.eq.0).and.(typSurf.eq.101)) then !SLP
206	isec1(6)=151 ! indicatorOfParameter
207	elseif ((parCat.eq.2).and.(parNum.eq.2).and.(typSurf.eq.103)) then ! 10U
208	isec1(6)=165 ! indicatorOfParameter
209	elseif ((parCat.eq.2).and.(parNum.eq.3).and.(typSurf.eq.103)) then ! 10V
210	isec1(6)=166 ! indicatorOfParameter
211	elseif ((parCat.eq.0).and.(parNum.eq.0).and.(typSurf.eq.103)) then ! 2T
212	isec1(6)=167 ! indicatorOfParameter
213	elseif ((parCat.eq.0).and.(parNum.eq.6).and.(typSurf.eq.103)) then ! 2D
214	isec1(6)=168 ! indicatorOfParameter
215	elseif ((parCat.eq.1).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SD
216	isec1(6)=141 ! indicatorOfParameter
217	conversion_factor=1000.
218	elseif ((parCat.eq.6).and.(parNum.eq.1) .or. parId .eq. 164) then ! CC
219	isec1(6)=164 ! indicatorOfParameter
220	elseif ((parCat.eq.1).and.(parNum.eq.9) .or. parId .eq. 142) then ! LSP
221	isec1(6)=142 ! indicatorOfParameter
222	elseif ((parCat.eq.1).and.(parNum.eq.10)) then ! CP
223	isec1(6)=143 ! indicatorOfParameter
224	conversion_factor=1000.
225	elseif ((parCat.eq.0).and.(parNum.eq.11).and.(typSurf.eq.1)) then ! SHF
226	isec1(6)=146 ! indicatorOfParameter
227	elseif ((parCat.eq.4).and.(parNum.eq.9).and.(typSurf.eq.1)) then ! SR
228	isec1(6)=176 ! indicatorOfParameter
229	! elseif ((parCat.eq.2).and.(parNum.eq.17) .or. parId .eq. 180) then ! EWSS --wrong
230	elseif ((parCat.eq.2).and.(parNum.eq.38) .or. parId .eq. 180) then ! EWSS --correct
231	isec1(6)=180 ! indicatorOfParameter
232	! elseif ((parCat.eq.2).and.(parNum.eq.18) .or. parId .eq. 181) then ! NSSS --wrong
233	elseif ((parCat.eq.2).and.(parNum.eq.37) .or. parId .eq. 181) then ! NSSS --correct
234	isec1(6)=181 ! indicatorOfParameter
235	elseif ((parCat.eq.3).and.(parNum.eq.4)) then ! ORO
236	isec1(6)=129 ! indicatorOfParameter
237	elseif ((parCat.eq.3).and.(parNum.eq.7) .or. parId .eq. 160) then ! SDO
238	isec1(6)=160 ! indicatorOfParameter
239	elseif ((discipl.eq.2).and.(parCat.eq.0).and.(parNum.eq.0).and. &
240	(typSurf.eq.1)) then ! LSM
241	isec1(6)=172 ! indicatorOfParameter
242	elseif (parNum.eq.152) then
243	isec1(6)=152 ! avoid warning for lnsp
244	else
245	print,'**WARNING: undefined GRiB2 message found!',discipl, &
246	parCat,parNum,typSurf
247	endif
248	if(parId .ne. isec1(6) .and. parId .ne. 77) then
249	write(,) 'parId',parId, 'isec1(6)',isec1(6)
250	! stop
251	endif
252
253	endif
254
255	!HSO get the size and data of the values array
256	if (isec1(6).ne.-1) then
257	call grib_get_real4_array(igrib,'values',zsec4,iret)
258	call grib_check(iret,gribFunction,gribErrorMsg)
259	endif
260
261	!HSO get the required fields from section 2 in a gribex compatible manner
262	if (ifield.eq.1) then
263	call grib_get_int(igrib,'numberOfPointsAlongAParallel',isec2(2),iret)
264	call grib_check(iret,gribFunction,gribErrorMsg)
265	call grib_get_int(igrib,'numberOfPointsAlongAMeridian',isec2(3),iret)
266	call grib_check(iret,gribFunction,gribErrorMsg)
267	call grib_get_int(igrib,'numberOfVerticalCoordinateValues',isec2(12))
268	call grib_check(iret,gribFunction,gribErrorMsg)
269	! CHECK GRID SPECIFICATIONS
270	if(isec2(2).ne.nxfield) stop 'READWIND: NX NOT CONSISTENT'
271	if(isec2(3).ne.ny) stop 'READWIND: NY NOT CONSISTENT'
272	if(isec2(12)/2-1.ne.nlev_ec) &
273	stop 'READWIND: VERTICAL DISCRETIZATION NOT CONSISTENT'
274	endif ! ifield
275
276	!HSO get the second part of the grid dimensions only from GRiB1 messages
277	if (isec1(6) .eq. 167 .and. (gotGrid.eq.0)) then
278	call grib_get_real8(igrib,'longitudeOfFirstGridPointInDegrees', &
279	xauxin,iret)
280	call grib_check(iret,gribFunction,gribErrorMsg)
281	call grib_get_real8(igrib,'latitudeOfLastGridPointInDegrees', &
282	yauxin,iret)
283	call grib_check(iret,gribFunction,gribErrorMsg)
284	if (xauxin.gt.180.) xauxin=xauxin-360.0
285	if (xauxin.lt.-180.) xauxin=xauxin+360.0
286
287	xaux=xauxin+real(nxshift)*dx
288	yaux=yauxin
289	if (xaux.gt.180.) xaux=xaux-360.0
290	if(abs(xaux-xlon0).gt.eps) &
291	stop 'READWIND: LOWER LEFT LONGITUDE NOT CONSISTENT'
292	if(abs(yaux-ylat0).gt.eps) &
293	stop 'READWIND: LOWER LEFT LATITUDE NOT CONSISTENT'
294	gotGrid=1
295	endif ! gotGrid
296
297	do j=0,nymin1
298	do i=0,nxfield-1
299	k=isec1(8)
300	if(isec1(6).eq.130) tth(i,j,nlev_ec-k+2,n)= &!! TEMPERATURE
301	zsec4(nxfield*(ny-j-1)+i+1)
302	if(isec1(6).eq.131) uuh(i,j,nlev_ec-k+2)= &!! U VELOCITY
303	zsec4(nxfield*(ny-j-1)+i+1)
304	if(isec1(6).eq.132) vvh(i,j,nlev_ec-k+2)= &!! V VELOCITY
305	zsec4(nxfield*(ny-j-1)+i+1)
306	if(isec1(6).eq.133) then !! SPEC. HUMIDITY
307	qvh(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
308	if (qvh(i,j,nlev_ec-k+2,n) .lt. 0.) &
309	qvh(i,j,nlev_ec-k+2,n) = 0.
310	! this is necessary because the gridded data may contain
311	! spurious negative values
312	endif
313	if(isec1(6).eq.134) ps(i,j,1,n)= &!! SURF. PRESS.
314	zsec4(nxfield*(ny-j-1)+i+1)
315
316	if(isec1(6).eq.135) wwh(i,j,nlev_ec-k+1)= &!! W VELOCITY
317	zsec4(nxfield*(ny-j-1)+i+1)
318	if(isec1(6).eq.141) sd(i,j,1,n)= &!! SNOW DEPTH
319	zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
320	if(isec1(6).eq.151) msl(i,j,1,n)= &!! SEA LEVEL PRESS.
321	zsec4(nxfield*(ny-j-1)+i+1)
322	if(isec1(6).eq.164) tcc(i,j,1,n)= &!! CLOUD COVER
323	zsec4(nxfield*(ny-j-1)+i+1)
324	if(isec1(6).eq.165) u10(i,j,1,n)= &!! 10 M U VELOCITY
325	zsec4(nxfield*(ny-j-1)+i+1)
326	if(isec1(6).eq.166) v10(i,j,1,n)= &!! 10 M V VELOCITY
327	zsec4(nxfield*(ny-j-1)+i+1)
328	if(isec1(6).eq.167) tt2(i,j,1,n)= &!! 2 M TEMPERATURE
329	zsec4(nxfield*(ny-j-1)+i+1)
330	if(isec1(6).eq.168) td2(i,j,1,n)= &!! 2 M DEW POINT
331	zsec4(nxfield*(ny-j-1)+i+1)
332	if(isec1(6).eq.142) then !! LARGE SCALE PREC.
333	lsprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
334	if (lsprec(i,j,1,n).lt.0.) lsprec(i,j,1,n)=0.
335	endif
336	if(isec1(6).eq.143) then !! CONVECTIVE PREC.
337	convprec(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)/conversion_factor
338	if (convprec(i,j,1,n).lt.0.) convprec(i,j,1,n)=0.
339	endif
340	if(isec1(6).eq.146) sshf(i,j,1,n)= &!! SENS. HEAT FLUX
341	zsec4(nxfield*(ny-j-1)+i+1)
342	if((isec1(6).eq.146).and.(zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) &
343	hflswitch=.true. ! Heat flux available
344	if(isec1(6).eq.176) then !! SOLAR RADIATION
345	ssr(i,j,1,n)=zsec4(nxfield*(ny-j-1)+i+1)
346	if (ssr(i,j,1,n).lt.0.) ssr(i,j,1,n)=0.
347	endif
348	if(isec1(6).eq.180) ewss(i,j)= &!! EW SURFACE STRESS
349	zsec4(nxfield*(ny-j-1)+i+1)
350	if(isec1(6).eq.181) nsss(i,j)= &!! NS SURFACE STRESS
351	zsec4(nxfield*(ny-j-1)+i+1)
352	if(((isec1(6).eq.180).or.(isec1(6).eq.181)).and. &
353	(zsec4(nxfield*(ny-j-1)+i+1).ne.0.)) strswitch=.true. ! stress available
354	!sec strswitch=.true.
355	if(isec1(6).eq.129) oro(i,j)= &!! ECMWF OROGRAPHY
356	zsec4(nxfield*(ny-j-1)+i+1)/ga
357	if(isec1(6).eq.160) excessoro(i,j)= &!! STANDARD DEVIATION OF OROGRAPHY
358	zsec4(nxfield*(ny-j-1)+i+1)
359	if(isec1(6).eq.172) lsm(i,j)= &!! ECMWF LAND SEA MASK
360	zsec4(nxfield*(ny-j-1)+i+1)
361	if(isec1(6).eq.131) iumax=max(iumax,nlev_ec-k+1)
362	if(isec1(6).eq.135) iwmax=max(iwmax,nlev_ec-k+1)
363	!ZHG READING CLOUD FIELDS ASWELL
364	if(isec1(6).eq.246) then !! CLWC Cloud liquid water content [kg/kg]
365	clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
366	readclouds = .true.
367	sumclouds=.false.
368	endif
369	if(isec1(6).eq.247) then !! CIWC Cloud ice water content
370	ciwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
371	endif
372	!ZHG end
373	!ESO read qc (=clwc+ciwc)
374	if(isec1(6).eq.201031) then !! QC Cloud liquid water content [kg/kg]
375	clwch(i,j,nlev_ec-k+2,n)=zsec4(nxfield*(ny-j-1)+i+1)
376	readclouds=.true.
377	sumclouds=.true.
378	endif
379
380	end do
381	end do
382
383	call grib_release(igrib)
384	goto 10 !! READ NEXT LEVEL OR PARAMETER
385	!
386	! CLOSING OF INPUT DATA FILE
387	!
388
389	50 call grib_close_file(ifile)
390
391	!error message if no fields found with correct first longitude in it
392	if (gotGrid.eq.0) then
393	print,'**ERROR: input file needs to contain GRiB1 formatted'// &
394	'messages'
395	stop
396	endif
397
398	! if(levdiff2.eq.0) then
399	iwmax=iwmax+1
400	do i=0,nxmin1
401	do j=0,nymin1
402	wwh(i,j,iwmax)=0.
403	end do
404	end do
405	! endif
406
407	! For global fields, assign the leftmost data column also to the rightmost
408	! data column; if required, shift whole grid by nxshift grid points
409	!*************************************************************************
410
411	if (xglobal) then
412	call shift_field_0(ewss,nxfield,ny)
413	call shift_field_0(nsss,nxfield,ny)
414	call shift_field_0(oro,nxfield,ny)
415	call shift_field_0(excessoro,nxfield,ny)
416	call shift_field_0(lsm,nxfield,ny)
417	! call shift_field(ps,nxfield,ny,1,1,2,n)
418	! call shift_field(sd,nxfield,ny,1,1,2,n)
419	! call shift_field(msl,nxfield,ny,1,1,2,n)
420	! call shift_field(tcc,nxfield,ny,1,1,2,n)
421	! call shift_field(u10,nxfield,ny,1,1,2,n)
422	! call shift_field(v10,nxfield,ny,1,1,2,n)
423	! call shift_field(tt2,nxfield,ny,1,1,2,n)
424	! call shift_field(td2,nxfield,ny,1,1,2,n)
425	! call shift_field(lsprec,nxfield,ny,1,1,2,n)
426	! call shift_field(convprec,nxfield,ny,1,1,2,n)
427	! call shift_field(sshf,nxfield,ny,1,1,2,n)
428	! call shift_field(ssr,nxfield,ny,1,1,2,n)
429	! call shift_field(tth,nxfield,ny,nuvzmax,nuvz,2,n)
430	! call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,2,n)
431	call shift_field(uuh,nxfield,ny,nuvzmax,nuvz,1,1)
432	call shift_field(vvh,nxfield,ny,nuvzmax,nuvz,1,1)
433	call shift_field(wwh,nxfield,ny,nwzmax,nwz,1,1)
434	! numwfmem
435	call shift_field(ps,nxfield,ny,1,1,numwfmem,n)
436	call shift_field(sd,nxfield,ny,1,1,numwfmem,n)
437	call shift_field(msl,nxfield,ny,1,1,numwfmem,n)
438	call shift_field(tcc,nxfield,ny,1,1,numwfmem,n)
439	call shift_field(u10,nxfield,ny,1,1,numwfmem,n)
440	call shift_field(v10,nxfield,ny,1,1,numwfmem,n)
441	call shift_field(tt2,nxfield,ny,1,1,numwfmem,n)
442	call shift_field(td2,nxfield,ny,1,1,numwfmem,n)
443	call shift_field(lsprec,nxfield,ny,1,1,numwfmem,n)
444	call shift_field(convprec,nxfield,ny,1,1,numwfmem,n)
445	call shift_field(sshf,nxfield,ny,1,1,numwfmem,n)
446	call shift_field(ssr,nxfield,ny,1,1,numwfmem,n)
447	call shift_field(tth,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
448	call shift_field(qvh,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
449	!hg
450	call shift_field(clwch,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
451	if (.not.sumclouds) call shift_field(ciwch,nxfield,ny,nuvzmax,nuvz,numwfmem,n)
452	!hg end
453
454	endif
455
456	do i=0,nxmin1
457	do j=0,nymin1
458	surfstr(i,j,1,n)=sqrt(ewss(i,j)2+nsss(i,j)2)
459	end do
460	end do
461
462	if (mp_measure_time) call mpif_mtime('readwind',1)
463
464	if ((.not.hflswitch).or.(.not.strswitch)) then
465	write(,) 'WARNING: No flux data contained in GRIB file ', &
466	wfname(indj)
467
468	! CALCULATE USTAR AND SSHF USING THE PROFILE METHOD
469	! As ECMWF has increased the model resolution, such that now the first model
470	! level is at about 10 m (where 10-m wind is given), use the 2nd ECMWF level
471	! (3rd model level in FLEXPART) for the profile method
472	!***************************************************************************
473
474	do i=0,nxmin1
475	do j=0,nymin1
476	plev1=akz(3)+bkz(3)*ps(i,j,1,n)
477	pmean=0.5*(ps(i,j,1,n)+plev1)
478	tv=tth(i,j,3,n)(1.+0.61qvh(i,j,3,n))
479	fu=-r_air*tv/ga/pmean
480	hlev1=fu*(plev1-ps(i,j,1,n)) ! HEIGTH OF FIRST MODEL LAYER
481	ff10m= sqrt(u10(i,j,1,n)2+v10(i,j,1,n)2)
482	fflev1=sqrt(uuh(i,j,3)2+vvh(i,j,3)2)
483	call pbl_profile(ps(i,j,1,n),td2(i,j,1,n),hlev1, &
484	tt2(i,j,1,n),tth(i,j,3,n),ff10m,fflev1, &
485	surfstr(i,j,1,n),sshf(i,j,1,n))
486	if(sshf(i,j,1,n).gt.200.) sshf(i,j,1,n)=200.
487	if(sshf(i,j,1,n).lt.-400.) sshf(i,j,1,n)=-400.
488	end do
489	end do
490	endif
491
492
493	! Assign 10 m wind to model level at eta=1.0 to have one additional model
494	! level at the ground
495	! Specific humidity is taken the same as at one level above
496	! Temperature is taken as 2 m temperature
497	!**************************************************************************
498
499	do i=0,nxmin1
500	do j=0,nymin1
501	uuh(i,j,1)=u10(i,j,1,n)
502	vvh(i,j,1)=v10(i,j,1,n)
503	qvh(i,j,1,n)=qvh(i,j,2,n)
504	tth(i,j,1,n)=tt2(i,j,1,n)
505	end do
506	end do
507
508	if(iumax.ne.nuvz-1) stop 'READWIND: NUVZ NOT CONSISTENT'
509	if(iwmax.ne.nwz) stop 'READWIND: NWZ NOT CONSISTENT'
510
511	return
512
513	888 write(,) ' #### FLEXPART MODEL ERROR! WINDFIELD #### '
514	write(,) ' #### ',wfname(indj),' #### '
515	write(,) ' #### IS NOT GRIB FORMAT !!! #### '
516	stop 'Execution terminated'
517
518	end subroutine readwind_ecmwf
519

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

Download in other formats:

Original Format