Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

read_ncwrfout.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 10 years ago
sources for flexwrf v3.1
File size: 31.3 KB

Line
1	!***********************************************************************
2	!* Copyright 2012,2013 *
3	!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine, *
4	!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5	!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso, *
6	!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann, *
7	!* *
8	!* This file is part of FLEXPART WRF *
9	!* *
10	!* FLEXPART is free software: you can redistribute it and/or modify *
11	!* it under the terms of the GNU General Public License as published by*
12	!* the Free Software Foundation, either version 3 of the License, or *
13	!* (at your option) any later version. *
14	!* *
15	!* FLEXPART is distributed in the hope that it will be useful, *
16	!* but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	!* GNU General Public License for more details. *
19	!* *
20	!* You should have received a copy of the GNU General Public License *
21	!* along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
22	!***********************************************************************
23
24	!**********************************************************************
25	! FLEXPART SOURCE FILE READ_NCWRFOUT - CONTAINS *
26	! *
27	! SUBROUTINE READ_NCWRFOUT_GRIDINFO *
28	! SUBROUTINE READ_NCWRFOUT_1REALFIELD *
29	! SUBROUTINE READ_NCWRFOUT_1DATETIME *
30	! *
31	!**********************************************************************
32	! *
33	! AUTHOR: R. Easter & J. Fast, PNNL *
34	! DATE: 2005-autumn-?? *
35	! LAST UPDATE: same *
36	! *
37	!**********************************************************************
38	! *
39	! DESCRIPTION: *
40	! *
41	! These routines read the netcdf wrf output files. *
42	! *
43	! 13 JUNE 2007, add more arguments-- ext_scalar,pbl_physcis
44	!**********************************************************************
45
46
47	!-----------------------------------------------------------------------
48	subroutine read_ncwrfout_gridinfo( ierr, idiagaa, fnamenc, &
49	n_west_east, n_south_north, n_bottom_top, &
50	dx_met, dy_met, &
51	m_grid_id, m_parent_grid_id, m_parent_grid_ratio, &
52	i_parent_start, j_parent_start, &
53	map_proj_id, map_stdlon, map_truelat1, map_truelat2, &
54	ext_scalar,pbl_physics,microphysics)
55	!
56	! reads grid definition information from a netcdf wrf output file
57	!
58	! arguments
59	! ierr - output - if non-zero, an error occurred
60	! while opening or reading from the file
61	! idiagaa - input - if positive, testing diagnostics are printed
62	! fnamenc - input - path+filename of the wrf output file
63	! n_west_east - output - east_west dimension of the "T-grid"
64	! n_south_north - output - south_north dimension of the "T-grid"
65	! n_bottom_top - output - bottom_top dimension of the "T-grid"
66	! dx_met, dy_met - output - horizontal grid spacing (m)
67	! m_grid_id - output - grid id number
68	! m_parent_grid_id - output - grid id number of parent grid
69	! m_parent_grid_ratio - output - ratio of parent grid dxy to current grid dxy
70	! i_parent_start, j_parent_start - output - location of lower left corner
71	! of current grid relative to the parent grid.
72	! (if there is no parent grid, then the above 4 "...parent..." variables
73	! area set to -987.)
74	! map_proj_id - WRF map projection id (2=polar stereographic)
75	! map_stdlon - map projection standard longitude (deg)
76	! map_truelat1, truelat2 - map projection true latitudes (deg)
77	!
78	! ext_scalar -- dimension of ex_scalar
79	! pbl_physics -- type of PBL scheme
80	! microphysics -- micorphysice scheme used
81
82	include 'netcdf.inc'
83	! use netcdf
84	! implicit none
85
86	! arguments
87	integer :: ierr, idiagaa, &
88	n_west_east, n_south_north, n_bottom_top, &
89	m_grid_id, m_parent_grid_id, m_parent_grid_ratio, &
90	i_parent_start, j_parent_start, map_proj_id, &
91	ext_scalar,pbl_physics,microphysics
92
93	real :: dx_met, dy_met, map_stdlon, map_truelat1, map_truelat2
94
95	character() fnamenc
96
97	! local variables
98	integer,parameter :: maxdim=20
99	integer,parameter :: ibadaa= -987
100	integer,parameter :: xbadaa= -987
101
102	integer :: i, iatt, idimid_unlim, idum, iret, ivtype
103	integer :: l, lenatt, lendim(maxdim)
104	integer :: natts_tot, ncid, ndims_tot, nvars_tot
105	integer :: n_west_east_stag, n_south_north_stag, n_bottom_top_stag
106
107	real :: duma
108	real, allocatable, dimension(:) :: duma_alloc
109
110	character(len=80) :: dimname(maxdim)
111	character(len=80) :: attname
112	character(len=1000) :: dumch1000
113
114	! externals
115	! integer nf_close
116	! integer nf_inq
117	! integer nf_inq_dim
118	! integer nf_open
119
120
121
122	! initialize with "missing values"
123	n_west_east = ibadaa
124	n_south_north = ibadaa
125	n_bottom_top = ibadaa
126	dx_met = xbadaa
127	dy_met = xbadaa
128	m_grid_id = ibadaa
129	m_parent_grid_id = ibadaa
130	m_parent_grid_ratio = ibadaa
131	i_parent_start = ibadaa
132	j_parent_start = ibadaa
133	ext_scalar = ibadaa
134	pbl_physics = ibadaa
135	microphysics = ibadaa
136	!
137	! open the netcdf file
138	!
139	ncid = 10
140	!C write(,)'xxx inside read_ncwrfout.f fnamenc=',fnamenc
141	!C write(,)'mp_physics=',microphysics
142
143	! print*,'filename ',fnamenc
144	iret = nf_open( fnamenc, NF_NOWRITE, ncid )
145	if (iret .ne. nf_noerr) then
146	write(*,9100) 'error doing open 123',fnamenc
147	! print*, NF_NOWRITE, ncid,iret
148	ierr = -1
149	return
150	end if
151
152	9100 format( / '*** read_ncwrfout_gridinfo -- ', a / &
153	'file = ', a )
154	9110 format( / '*** read_ncwrfout_gridinfo -- ', a, 1x, i8 / &
155	'file = ', a )
156	9120 format( / '*** read_ncwrfout_gridinfo -- ', a, 2(1x,i8) / &
157	'file = ', a )
158
159	90030 format( a, 2i6, 2(2x,a) )
160
161	!
162	! get information on dimensions
163	!
164	iret = nf_inq( ncid, &
165	ndims_tot, nvars_tot, natts_tot, idimid_unlim )
166	if (iret .ne. nf_noerr) then
167	write(*,9100) 'error inquiring dimensions', fnamenc
168	ierr = -2
169	return
170	end if
171
172	n_west_east_stag = ibadaa
173	n_south_north_stag = ibadaa
174	n_bottom_top_stag = ibadaa
175
176	do i = 1, min(ndims_tot,maxdim)
177	iret = nf_inq_dim( ncid, i, dimname(i), lendim(i) )
178	if (iret .ne. nf_noerr) then
179	write(*,9110) 'error inquiring dimensions for dim#', i, fnamenc
180	ierr = -2
181	return
182	end if
183	end do
184
185	do i = 1, min(ndims_tot,maxdim)
186	if (dimname(i) .eq. 'west_east') &
187	n_west_east = lendim(i)
188	if (dimname(i) .eq. 'south_north') &
189	n_south_north = lendim(i)
190	if (dimname(i) .eq. 'bottom_top') &
191	n_bottom_top = lendim(i)
192	if (dimname(i) .eq. 'west_east_stag') &
193	n_west_east_stag = lendim(i)
194	if (dimname(i) .eq. 'south_north_stag') &
195	n_south_north_stag = lendim(i)
196	if (dimname(i) .eq. 'bottom_top_stag') &
197	n_bottom_top_stag = lendim(i)
198	if (dimname(i) .eq. 'ext_scalar') &
199	ext_scalar = lendim(i)
200
201	end do
202
203	if (idiagaa .gt. 0) then
204	write(*,9100) 'diagnostics', fnamenc
205	do i = 1, min(ndims_tot,maxdim)
206	write(*,90030) 'dim #, len, name =', &
207	i, lendim(i), dimname(i)
208	end do
209	end if
210
211	if ((n_west_east .le. 0) .or. &
212	(n_west_east+1 .ne. n_west_east_stag)) then
213	write(*,9120) 'bad n_west_east, n_west_east_stag = ', &
214	n_west_east, n_west_east_stag, fnamenc
215	ierr = -3
216	return
217	end if
218
219	if ((n_south_north .le. 0) .or. &
220	(n_south_north+1 .ne. n_south_north_stag)) then
221	write(*,9120) 'bad n_south_north, n_south_north_stag = ', &
222	n_south_north, n_south_north_stag, fnamenc
223	ierr = -3
224	return
225	end if
226
227	if ((n_bottom_top .le. 0) .or. &
228	(n_bottom_top+1 .ne. n_bottom_top_stag)) then
229	write(*,9120) 'bad n_bottom_top, n_bottom_top_stag = ', &
230	n_bottom_top, n_bottom_top_stag, fnamenc
231	ierr = -3
232	return
233	end if
234
235	!
236	! get information on global attributes
237	!
238
239	! first just do diagnostics
240	if (idiagaa .gt. 0) then
241	write(,)
242	write(,) 'attribute #, name, type, value'
243	end if
244	do iatt = 1, natts_tot
245	iret = nf_inq_attname( ncid, nf_global, iatt, attname)
246	if (iret .ne. nf_noerr) goto 3600
247	iret = nf_inq_att( ncid, nf_global, attname, ivtype, lenatt )
248	if (iret .ne. nf_noerr) goto 3600
249	if (ivtype .eq. 2) then
250	iret = nf_get_att_text( ncid, nf_global, attname, dumch1000 )
251	if (iret .ne. nf_noerr) goto 3600
252	i = max(1,min(1000,lenatt))
253	if (idiagaa .gt. 0) write(*,91010) &
254	iatt, attname(1:40), ivtype, lenatt, dumch1000(1:i)
255	else if (ivtype .eq. 4) then
256	iret = nf_get_att_int( ncid, nf_global, attname, idum )
257	if (iret .ne. nf_noerr) goto 3600
258	if (idiagaa .gt. 0) write(*,91020) &
259	iatt, attname(1:40), ivtype, lenatt, idum
260	else if ((ivtype .eq. 5) .and. (lenatt .eq. 1)) then
261	iret = nf_get_att_real( ncid, nf_global, attname, duma )
262	if (iret .ne. nf_noerr) goto 3600
263	if (idiagaa .gt. 0) write(*,91030) &
264	iatt, attname(1:40), ivtype, lenatt, duma
265	else if ((ivtype .eq. 5) .and. (lenatt .gt. 1)) then
266	allocate( duma_alloc(lenatt) )
267	iret = nf_get_att_real( ncid, nf_global, attname, duma_alloc )
268	if (iret .ne. nf_noerr) goto 3600
269	if (idiagaa .gt. 0) then
270	write(*,91010) iatt, attname(1:40), ivtype, lenatt
271	write(*,91040) (duma_alloc(i), i=1,lenatt)
272	end if
273	deallocate( duma_alloc )
274	else
275	if (idiagaa .gt. 0) write(*,'(i4,1x,a,2(1x,i6))') &
276	iatt, attname(1:40), ivtype, lenatt
277	goto 3400
278	endif
279
280	if (attname .eq. 'GRID_ID') then
281	m_grid_id = idum
282	else if (attname .eq. 'PARENT_ID') then
283	m_parent_grid_id = idum
284	else if (attname .eq. 'PARENT_GRID_RATIO') then
285	m_parent_grid_ratio = idum
286	else if (attname .eq. 'I_PARENT_START') then
287	i_parent_start = idum
288	else if (attname .eq. 'J_PARENT_START') then
289	j_parent_start = idum
290	else if (attname .eq. 'DX') then
291	dx_met = duma
292	else if (attname .eq. 'DY') then
293	dy_met = duma
294	else if (attname .eq. 'MAP_PROJ') then
295	map_proj_id = idum
296	else if (attname .eq. 'STAND_LON') then
297	map_stdlon = duma
298	else if (attname .eq. 'TRUELAT1') then
299	map_truelat1 = duma
300	else if (attname .eq. 'TRUELAT2') then
301	map_truelat2 = duma
302	else if (attname .eq. 'BL_PBL_PHYSICS') then
303	pbl_physics = idum
304	else if (attname .eq. 'MP_PHYSICS') then
305	microphysics = idum
306	end if
307	enddo
308	3400 continue
309	91010 format( i4, 1x, a, 2(1x,i6), 1x, a )
310	91020 format( i4, 1x, a, 2(1x,i6), 1x, i10 )
311	91030 format( i4, 1x, a, 2(1x,i6), 1x, 1pe12.4 )
312	91040 format(( 12x, 5(1pe12.4) ))
313
314	goto 3900
315
316	3600 write(*,9110) 'error inquiring attribute', iatt, fnamenc
317	stop
318
319	3900 continue
320
321	!C write(,)'mp_physics=',microphysics,pbl_physics
322
323	!
324	! close and return
325	!
326	iret = nf_close( ncid )
327	ierr = 0
328
329	return
330	end subroutine read_ncwrfout_gridinfo
331
332
333
334	!-----------------------------------------------------------------------
335	subroutine read_ncwrfout_1datetime( ierr, fnamenc, &
336	itime, jyyyymmdd, jhhmmss )
337	!
338	! a wrf output file may contain data at multiple time. This routine returns
339	! the date & time of the "itime" data group in the file.
340	!
341	! arguments
342	! ierr - output - if non-zero, an error occurred
343	! while opening or reading from the file,
344	! or itime < 0, or itime > number of times in the file.
345	! fnamenc - input - path+filename of the wrf output file
346	! itime - input - specifies which data/time to return.
347	! 1 for first, 2 for second, ...
348	! jyyyymmdd - output - date as 8 decimal digits (yyyymmdd).
349	! yyyy=year, mm=month, dd=day of month.
350	! jhhmmss - output - time of day as 6 decimal digits (hhmmss).
351	! hh=hour, mm=minute, ss=second
352	! if (jyyyymmdd=jhhmmss=-1, then ierr is non-zero, and vice-versa)
353	!
354
355	! use netcdf
356	include 'netcdf.inc'
357	!implicit none
358
359
360	! arguments
361	integer :: ierr, itime, jyyyymmdd, jhhmmss
362	character() fnamenc
363
364	! local variables
365	integer,parameter :: ndims_maxbb=4 ! max number of dimensions for a variable
366
367	integer :: i, id_var, iret, itype_var
368	integer :: iduma, idumb, idumc
369	integer :: id_dim(ndims_maxbb)
370	integer :: istart(ndims_maxbb), icount(ndims_maxbb)
371	integer :: lendim(ndims_maxbb)
372	integer :: natts_tot, ncid, ndims
373
374	character(len=32) timetext
375	character(len=80) varname, varnamenc
376
377	! externals
378	! integer nf_close
379	! integer nf_inq
380	! integer nf_inq_dim
381	! integer nf_open
382
383
384
385	jyyyymmdd = -1
386	jhhmmss = -1
387
388	!
389	! open the netcdf file
390	!
391	ncid = 10
392	iret = nf_open( fnamenc, NF_NOWRITE, ncid )
393	if (iret .ne. nf_noerr) then
394	write(*,9100) 'error doing open 370', fnamenc
395	ierr = -1
396	goto 8100
397	end if
398
399	9100 format( / '*** read_ncwrfout_1datetime -- ', a / &
400	'file = ', a )
401	9110 format( / '*** read_ncwrfout_1datetime -- ', a, 1x, i8 / &
402	'file = ', a )
403	9120 format( / '*** read_ncwrfout_1datetime -- ', a, 2(1x,i8) / &
404	'file = ', a )
405	9130 format( / '*** read_ncwrfout_1datetime -- ', a, 3(1x,i8) / &
406	'file = ', a )
407	9115 format( / '*** read_ncwrfout_1datetime -- ', a / a, 1x, i8 / &
408	'file = ', a )
409	9125 format( / '*** read_ncwrfout_1datetime -- ', a / a, 2(1x,i8) / &
410	'file = ', a )
411	9135 format( / '*** read_ncwrfout_1datetime -- ', a / a, 3(1x,i8) / &
412	'file = ', a )
413
414	90030 format( a, 2i6, 2(2x,a) )
415
416	!
417	! get information on the variable
418	!
419	varname = 'Times'
420	iret = nf_inq_varid( ncid, varname, id_var )
421	if (iret .ne. nf_noerr) then
422	write(*,9100) 'error inquiring var id for ' // varname, fnamenc
423	ierr = -1
424	goto 8100
425	end if
426
427	iret = nf_inq_var( ncid, id_var, &
428	varnamenc, itype_var, ndims, id_dim, natts_tot )
429	if (iret .ne. nf_noerr) then
430	write(*,9100) 'error inquiring var info for ' // varname, fnamenc
431	ierr = -1
432	goto 8100
433	end if
434
435	! check variable type
436	if (itype_var .ne. nf_char) then
437	write(*,9110) 'var type wrong for ' // varname, &
438	itype_var, fnamenc
439	ierr = -1
440	goto 8100
441	end if
442
443
444	! check number of dimensions
445	if (ndims .ne. 2) then
446	write(*,9115) 'var ndims is wrong for ' // varname, &
447	'ndims =', ndims, fnamenc
448	ierr = -1
449	goto 8100
450	end if
451
452	! get sizes of dimensions
453	! dimension 1 = # of characters in date/time string
454	! dimension 2 = # of times in the file
455	do i = 1, ndims
456	iret = nf_inq_dimlen( ncid, id_dim(i), lendim(i) )
457	if (iret .ne. nf_noerr) then
458	write(*,9115) 'error inquiring var dim len for ' // varname, &
459	'idim =', i, fnamenc
460	ierr = -1
461	goto 8100
462	end if
463	end do
464
465	if (itime .lt. 1) then
466	ierr = -11
467	goto 8100
468	else if (itime .gt. lendim(2)) then
469	ierr = -12
470	goto 8100
471	end if
472
473	! get the data and extract the data & time
474	do i = 1, ndims_maxbb
475	istart(i) = 1
476	icount(i) = 1
477	end do
478	istart(1) = 1
479	icount(1) = lendim(1)
480	istart(2) = itime
481	icount(2) = 1
482	iret = nf_get_vara_text( ncid, id_var, istart, icount, timetext )
483	if (iret .ne. nf_noerr) then
484	write(*,9100) 'error reading var data for ' // varname, &
485	fnamenc
486	ierr = -1
487	goto 8100
488	end if
489
490	read( timetext, '(i4,1x,i2,1x,i2)', iostat=iret ) &
491	iduma, idumb, idumc
492	if (iret .ne. 0) then
493	write(*,9125) &
494	'error reading from timetext = "' // timetext // '"', &
495	'itime, lendim(1) =', itime, lendim(1), fnamenc
496	ierr = -1
497	goto 8100
498	end if
499	jyyyymmdd = iduma10000 + idumb100 + idumc
500
501	read( timetext, '(11x,i2,1x,i2,1x,i2)', iostat=iret ) &
502	iduma, idumb, idumc
503	if (iret .ne. 0) then
504	write(*,9125) &
505	'error reading from timetext = "' // timetext // '"', &
506	'itime, lendim(1) =', itime, lendim(1), fnamenc
507	ierr = -1
508	goto 8100
509	end if
510	jhhmmss = iduma10000 + idumb100 + idumc
511
512	!
513	! success - close and return
514	!
515	iret = nf_close( ncid )
516	ierr = 0
517	return
518
519	!
520	! failure - close and return
521	!
522	8100 iret = nf_close( ncid )
523	return
524
525	end subroutine read_ncwrfout_1datetime
526
527
528
529	!-----------------------------------------------------------------------
530	subroutine read_ncwrfout_1realfield( ierr, idiagaa, fnamenc, &
531	varname, vardata, &
532	itime, &
533	ndims, ndims_exp, ndims_max, &
534	lendim, lendim_exp, lendim_max )
535	!
536	! reads of real (single precision) field at one time from a netcdf wrf output file
537	!
538	! arguments
539	! ierr - output - if non-zero, an error occurred
540	! while opening or reading from the file
541	! -1 = error opening file
542	! -2 = requested variable is not in the file
543	! -3 = error while inquiring about the variable
544	! -4 = variable type is other than real (single precision)
545	! ... = check below, in the code, for explanation of other ierr values.
546	! idiagaa - input - if positive, testing diagnostics are printed
547	! fnamenc - input - path+filename of the wrf output file
548	! varname - input - field name
549	! vardata - output - the data for the field
550	! itime - input - specifies which time to read.
551	! (1 for first time in the file, 2 for second, ...)
552	! ndims - output - number of (netcdf) dimensions of the field.
553	! This includes the time dimension.
554	! ndims_exp - input - expected number of dimensions of the field.
555	! An error occurs if ndims .ne. ndims_exp.
556	! ndims_max - input - The dimension/size of the lendim_... arrays.
557	! lendim - output - The size of each of the "ndims" dimensions.
558	! lendim_exp - input - The expected size of each of the "ndims" dimensions.
559	! If lendim_exp .gt. 0, then an error occurs if lendim .ne. lendim_exp.
560	! lendim_max - input - The maximum size of each dimension. These are
561	! The dimensions of the vardata array.
562	!
563	! use netcdf
564	include 'netcdf.inc'
565	!implicit none
566
567
568	! arguments
569	integer :: ierr, idiagaa, itime, &
570	ndims, ndims_exp, ndims_max, &
571	lendim(ndims_max), lendim_exp(ndims_max), &
572	lendim_max(ndims_max)
573
574	real :: vardata( lendim_max(1), lendim_max(2), lendim_max(3) )
575
576	character() fnamenc, varname
577
578	! local variables
579	integer,parameter :: ndims_maxbb=4 ! max number of dimensions for a variable
580	integer,parameter :: ibadaa=-987
581	integer,parameter :: xbadaa=-987
582	integer,parameter :: ijktestmax=500
583
584	integer :: i, iatt, id_var, iret, itype_var
585	integer :: itot, ijktot, ii
586	integer :: id_dim(ndims_maxbb)
587	integer :: istart(ndims_maxbb), icount(ndims_maxbb)
588	integer :: j, jtot, jj
589	integer :: k, ktot, kk
590	integer :: l, lenatt
591	integer :: lendim_use(ndims_maxbb)
592	integer :: m
593	integer :: natts_tot, ncid, ndiffa, ndiffb
594
595	real :: duma, dumb
596	real :: testavg(ijktestmax,3,2), &
597	testmin(ijktestmax,3,2), testmax(ijktestmax,3,2)
598
599	character(len=80) varnamenc
600	character(len=80) dimname(ndims_maxbb)
601
602	! externals
603	! integer nf_close
604	! integer nf_inq
605	! integer nf_inq_dim
606	! integer nf_open
607
608	! print*,'dim array',lendim_max(1), lendim_max(2), lendim_max(3)
609
610	!
611	! open the netcdf file
612	!
613	ncid = 10
614	iret = nf_open( fnamenc, NF_NOWRITE, ncid )
615	if (iret .ne. nf_noerr) then
616	write(*,9100) 'error doing open 592', fnamenc
617	ierr = -1
618	return
619	end if
620
621	9100 format( / '*** read_ncwrfout_1realfield -- ', a / &
622	'file = ', a )
623	9110 format( / '*** read_ncwrfout_1realfield -- ', a, 1x, i8 / &
624	'file = ', a )
625	9120 format( / '*** read_ncwrfout_1realfield -- ', a, 2(1x,i8) / &
626	'file = ', a )
627	9130 format( / '*** read_ncwrfout_1realfield -- ', a, 3(1x,i8) / &
628	'file = ', a )
629	9115 format( / '*** read_ncwrfout_1realfield -- ', a / a, 1x, i8 / &
630	'file = ', a )
631	9125 format( / '*** read_ncwrfout_1realfield -- ', a / a, 2(1x,i8) / &
632	'file = ', a )
633	9135 format( / '*** read_ncwrfout_1realfield -- ', a / a, 3(1x,i8) / &
634	'file = ', a )
635
636	90030 format( a, 2i6, 2(2x,a) )
637
638	!
639	! get information on the variable
640	!
641	iret = nf_inq_varid( ncid, varname, id_var )
642	if (iret .ne. nf_noerr) then
643	! write(*,9100) 'error inquiring var id for ' // varname, fnamenc
644	ierr = -2
645	goto 8100
646	end if
647
648	iret = nf_inq_var( ncid, id_var, &
649	varnamenc, itype_var, ndims, id_dim, natts_tot )
650	if (iret .ne. nf_noerr) then
651	write(*,9100) 'error inquiring var info for ' // varname, fnamenc
652	ierr = -3
653	goto 8100
654	end if
655
656	! check variable type
657	if (itype_var .ne. nf_real) then
658	write(*,9110) 'var type wrong for ' // varname, &
659	itype_var, fnamenc
660	ierr = -4
661	goto 8100
662	end if
663
664
665	! check number of dimensions
666	if (ndims_exp .le. 0) then
667	write(*,9115) &
668	'bad ndims_exp for ' // varname, &
669	'ndims_exp =', ndims_exp, fnamenc
670	ierr = -11
671	goto 8100
672	end if
673	if (ndims .ne. ndims_exp) then
674	write(*,9125) 'var ndims mismatch for ' // varname, &
675	'ndims_exp, ndims =', ndims_exp, ndims, fnamenc
676	ierr = -12
677	goto 8100
678	end if
679	if (ndims .gt. ndims_max) then
680	write(*,9125) 'var ndims > ndims_max for ' // varname, &
681	'ndims, ndims_max =', ndims, ndims_max, fnamenc
682	ierr = -13
683	goto 8100
684	end if
685	if (ndims .gt. ndims_maxbb) then
686	write(*,9125) 'var ndims > ndims_maxbb for ' // varname, &
687	'ndims, ndims_maxbb =', ndims, ndims_maxbb, fnamenc
688	ierr = -14
689	goto 8100
690	end if
691
692	! check size of each dimension
693	do i = 1, ndims_exp
694	iret = nf_inq_dimlen( ncid, id_dim(i), lendim(i) )
695	if (iret .ne. nf_noerr) then
696	write(*,9110) 'error inquiring var dim len for ' // varname, &
697	i, fnamenc
698	ierr = -15
699	goto 8100
700	end if
701	if ((i .lt. ndims_exp) .and. (lendim_exp(i) .gt. 0) .and. &
702	(lendim(i) .ne. lendim_exp(i))) then
703	! print*,i,ndims_exp,lendim_exp(i),lendim(i),lendim_exp(i)
704	write(*,9130) 'var lendim mismatch for ' // varname, &
705	i, lendim_exp(i), lendim(i), fnamenc
706	ierr = -16
707	goto 8100
708	end if
709	if ((i .lt. ndims_exp) .and. &
710	(lendim(i) .gt. lendim_max(i))) then
711	write(*,9130) 'var lendim too big for ' // varname, &
712	i, lendim_max(i), lendim(i), fnamenc
713	ierr = -17
714	goto 8100
715	end if
716	if ((i .eq. ndims_exp) .and. (lendim(i) .lt. itime)) then
717	write(*,9130) 'var itime < ntimes for ' // varname, &
718	i, itime, lendim(i), fnamenc
719	ierr = -18
720	goto 8100
721	end if
722	end do
723
724	! do diagnostics on the dimensions
725	if (idiagaa .gt. 0) then
726	write(*,'(/a)') &
727	'read_ncwrfout_1realfield - dim info for var = ' // &
728	varname(1:20)
729	do i = 1, ndims
730	iret = nf_inq_dim( ncid, id_dim(i), dimname(i), lendim(i) )
731	if (iret .ne. nf_noerr) then
732	write(*,9115) 'error inquiring var dim info for ' // varname, &
733	'idim =', i, fnamenc
734	ierr = -19
735	goto 8100
736	end if
737	write(*,'(a,3i5,2x,a)') ' i,id,len,name =', &
738	i, id_dim(i), lendim(i), dimname(i)(1:32)
739	end do
740	end if
741
742	!
743	! get the data
744	!
745	do i = 1, ndims_maxbb
746	istart(i) = 1
747	icount(i) = 1
748	end do
749	do i = 1, ndims_exp - 1
750	istart(i) = 1
751	icount(i) = lendim(i)
752	end do
753	! in wrfout files, the last dimension should always be time
754	istart(ndims_exp) = itime
755	icount(ndims_exp) = 1
756
757	iret = nf_get_vara_real( ncid, id_var, istart, icount, vardata )
758	if (iret .ne. nf_noerr) then
759	write(*,9120) 'error reading var data for ' // varname, &
760	fnamenc
761	ierr = -21
762	goto 8100
763	end if
764
765	!
766	! reorder the data
767	!
768	lendim_use(1) = lendim(1)
769	lendim_use(2) = 1
770	if (ndims_exp .ge. 3) lendim_use(2) = lendim(2)
771	lendim_use(3) = 1
772	if (ndims_exp .ge. 4) lendim_use(3) = lendim(3)
773
774	! print*,'dim to go',lendim_use(1),lendim_use(2),lendim_use(3)
775	! print*,'value in',lendim_use, lendim_max
776	call reorder_ncwrfout_1realfield( ierr, idiagaa, &
777	varname, vardata, vardata, &
778	lendim_use, lendim_max )
779
780	if (ierr .ne. 0) then
781	print*, 'error re-ordering var data for ', varname, &
782	fnamenc
783	ierr = -22
784	goto 8100
785	end if
786
787	!
788	! success - close and return
789	!
790	iret = nf_close( ncid )
791	ierr = 0
792	return
793
794	!
795	! error - close and return
796	!
797	8100 iret = nf_close( ncid )
798	return
799
800	end subroutine read_ncwrfout_1realfield
801
802
803
804	!-----------------------------------------------------------------------
805	subroutine reorder_ncwrfout_1realfield( ierr, idiagaa, &
806	varname, vardata_in, vardata_out, &
807	lendim_use, lendim_max )
808	!
809	! reorders a real (single precision) field
810	! the nf_get_vara_real loads the data for a field into
811	! a contiguous block of memory starting at vardata(1,1,1)
812	! it does not know if perhaps lendim() < lendim_max()
813	! this routine corrects for that, so that the data are
814	! loaded into non-contiguous blocks when lendim() < lendim_max()
815	!
816	! arguments
817	! ierr - output - if non-zero, an error occurred while re-ordering
818	! idiagaa - input - if positive, testing diagnostics are printed
819	! varname - input - field name
820	!
821	! vardata_in - input - the data for the field
822	! vardata_out - output - the data for the field
823	! In the calling program, vardata_in & vardata_out are usually
824	! the same array. This routine "pretends" that they are
825	! different, and specifies their dimensions differently,
826	! to facilitate the reordering.
827	!
828	! lendim_use - input - The actual size of the spatial dimensions of the field.
829	! lendim_max - input - The actual spatial dimensions of the vardata array.
830	! Most wrf fields are spatially either 1d (z), 2d (xy), or 3d (xyz).
831	! For a 1d spatial field (e.g., z only), set
832	! lendim_use(1) = nz, lendim_max(1) = nz_max
833	! lendim_use(2) = 1, lendim_max(2) = 1
834	! lendim_use(3) = 1, lendim_max(3) = 1
835	! For a 2d spatial field (e.g., xy only), set
836	! lendim_use(1) = nx, lendim_max(1) = nx_max
837	! lendim_use(2) = ny lendim_max(2) = ny_max
838	! lendim_use(3) = 1, lendim_max(3) = 1
839	! For a 3d spatial field (xyz), set
840	! lendim_use(1) = nx, lendim_max(1) = nx_max
841	! lendim_use(2) = ny lendim_max(2) = ny_max
842	! lendim_use(3) = nz lendim_max(3) = nz_max
843	!
844
845	! use netcdf
846	include 'netcdf.inc'
847	!implicit none
848
849
850	! arguments
851	integer :: ierr, idiagaa, &
852	lendim_use(3), lendim_max(3)
853
854	real :: vardata_in( lendim_use(1), lendim_use(2), lendim_use(3) )
855	real :: vardata_out( lendim_max(1), lendim_max(2), lendim_max(3) )
856
857	character() varname
858
859	! local variables
860	integer,parameter :: ijktestmax=500
861	integer,parameter :: check_reordering=1
862
863	integer :: i, j, k, m, n
864	integer :: itestend, jtestend, ktestend
865	integer :: ijk, ijktestend
866	integer :: ndiffa, ndiffb
867
868	real :: duma, dumb
869	real :: testavg(ijktestmax,3,2), &
870	testmin(ijktestmax,3,2), testmax(ijktestmax,3,2)
871
872	!
873	! the testavg/min/max are avg, min, and max values for
874	! a i (or j or k) fixed and j,k (or i,k or i,j) varying
875	! they are computed before and after the data has been reordered
876	! then compared at the end
877	! an error occurs if they do not match
878	!
879	! print*,'size in out',lendim_use(1), lendim_use(2), lendim_use(3)
880	! print*,'size in out',lendim_max(1), lendim_max(2), lendim_max(3)
881	if (check_reordering .gt. 0) then
882
883	do n = 1, 2
884	do m = 1, 3
885	do i = 1, ijktestmax
886	testavg(i,m,n) = 0.0
887	testmin(i,m,n) = +1.0e37
888	testmax(i,m,n) = -1.0e37
889	end do
890	end do
891	end do
892	! print*,varname
893	ktestend = min( ijktestmax, lendim_use(3) )
894	jtestend = min( ijktestmax, lendim_use(2) )
895	itestend = min( ijktestmax, lendim_use(1) )
896	! print*,'ktestend',ktestend
897	! pass 1 -- compute the test---(:,:,1) from vardata_in
898	do k = 1, ktestend
899	do j = 1, jtestend
900	do i = 1, itestend
901	duma = vardata_in(i,j,k)
902	testavg(i,1,1) = testavg(i,1,1) + duma
903	testmin(i,1,1) = min( testmin(i,1,1), duma )
904	testmax(i,1,1) = max( testmax(i,1,1), duma )
905	testavg(j,2,1) = testavg(j,2,1) + duma
906	testmin(j,2,1) = min( testmin(j,2,1), duma )
907	testmax(j,2,1) = max( testmax(j,2,1), duma )
908	testavg(k,3,1) = testavg(k,3,1) + duma
909	testmin(k,3,1) = min( testmin(k,3,1), duma )
910	testmax(k,3,1) = max( testmax(k,3,1), duma )
911	end do
912	end do
913	end do
914
915	end if ! if (check_reordering .gt. 0) then
916
917	! pass 2 -- shift the data values
918	! print*,'max',lendim_max
919	do k = lendim_use(3), 1, -1
920	do j = lendim_use(2), 1, -1
921	do i = lendim_use(1), 1, -1
922	! print*,i,j,k
923	vardata_out(i,j,k) = vardata_in(i,j,k)
924	end do
925	end do
926	end do
927
928	! pass 3 -- compute the test---(:,:,2) from vardata_out
929	if (check_reordering .gt. 0) then
930
931	do k = 1, ktestend
932	do j = 1, jtestend
933	do i = 1, itestend
934	duma = vardata_out(i,j,k)
935	testavg(i,1,2) = testavg(i,1,2) + duma
936	testmin(i,1,2) = min( testmin(i,1,2), duma )
937	testmax(i,1,2) = max( testmax(i,1,2), duma )
938	testavg(j,2,2) = testavg(j,2,2) + duma
939	testmin(j,2,2) = min( testmin(j,2,2), duma )
940	testmax(j,2,2) = max( testmax(j,2,2), duma )
941	testavg(k,3,2) = testavg(k,3,2) + duma
942	testmin(k,3,2) = min( testmin(k,3,2), duma )
943	testmax(k,3,2) = max( testmax(k,3,2), duma )
944	end do
945	end do
946	end do
947
948	! now compare the test---(:,:,1) & test---(:,:,2)
949	ndiffb = 0
950	do m = 1, 3
951	if (m .eq. 1) then
952	ijktestend = itestend
953	duma = 1.0/(jtestend*ktestend)
954	if (idiagaa .gt. 0) write(*,'(a,a)') varname(1:20), &
955	'i, testavg(i,1), testmin(i,1), testmax(i,1)'
956	else if (m .eq. 2) then
957	ijktestend = jtestend
958	duma = 1.0/(itestend*ktestend)
959	if (idiagaa .gt. 0) write(*,'(a,a)') varname(1:20), &
960	'j, testavg(j,2), testmin(j,2), testmax(j,2)'
961	else
962	ijktestend = ktestend
963	duma = 1.0/(itestend*jtestend)
964	if (idiagaa .gt. 0) write(*,'(a,a)') varname(1:20), &
965	'k, testavg(k,3), testmin(k,3), testmax(k,3)'
966	end if
967
968	ndiffa = 0
969	do ijk = 1, ijktestend
970	i = ijk
971	dumb = max( abs(testavg(i,m,1)), abs(testavg(i,m,2)) )*2.0e-7
972	if (abs(testavg(i,m,1)-testavg(i,m,2)) .gt. dumb) ndiffa = ndiffa + 1
973	dumb = max( abs(testmin(i,m,1)), abs(testmin(i,m,2)) )*2.0e-7
974	if (abs(testmin(i,m,1)-testmin(i,m,2)) .gt. dumb) ndiffa = ndiffa + 1
975	dumb = max( abs(testmax(i,m,1)), abs(testmax(i,m,2)) )*2.0e-7
976	if (abs(testmax(i,m,1)-testmax(i,m,2)) .gt. dumb) ndiffa = ndiffa + 1
977	end do
978
979	if (ndiffa .le. 0) then
980	if (idiagaa .gt. 0) write(,) ' *** no differences'
981	else
982	do ijk = 1, ijktestend
983	i = ijk
984	if (idiagaa .gt. 0) write(*,'(i3,1p,3(2x,2e11.3))') i, &
985	testavg(i,m,1)*duma, &
986	(testavg(i,m,1)-testavg(i,m,2))*duma, &
987	testmin(i,m,1), (testmin(i,m,1)-testmin(i,m,2)), &
988	testmax(i,m,1), (testmax(i,m,1)-testmax(i,m,2))
989	end do
990	end if
991
992	ndiffb = ndiffb + ndiffa
993	end do
994
995	if (ndiffb .gt. 0) then
996	ierr = -12
997	goto 8100
998	end if
999
1000	end if ! if (check_reordering .gt. 0) then
1001
1002	!
1003	! success
1004	!
1005	ierr = 0
1006	return
1007
1008	!
1009	! error
1010	!
1011	8100 return
1012
1013	end subroutine reorder_ncwrfout_1realfield
1014
1015
1016

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

Download in other formats:

Original Format