Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_si_io_nmm.F @ 1

Last change on this file since 1 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 112.6 KB

Line
1	MODULE module_si_io_nmm
2
3	USE module_optional_input
4
5	IMPLICIT NONE
6	!! FROM MODULE_KINDS
7
8	! The numerical data types defined in this module are:
9	! i_byte - specification kind for byte (1-byte) integer variable
10	! i_short - specification kind for short (2-byte) integer variable
11	! i_long - specification kind for long (4-byte) integer variable
12	! i_llong - specification kind for double long (8-byte) integer variable
13	! r_single - specification kind for single precision (4-byte) real variable
14	! r_double - specification kind for double precision (8-byte) real variable
15	! r_quad - specification kind for quad precision (16-byte) real variable
16	!
17	! i_kind - generic specification kind for default integer
18	! r_kind - generic specification kind for default floating point
19	!
20	!
21	! Integer type definitions below
22
23	! Integer types
24	integer, parameter, public :: i_byte = selected_int_kind(1) ! byte integer
25	integer, parameter, public :: i_short = selected_int_kind(4) ! short integer
26	integer, parameter, public :: i_long = selected_int_kind(8) ! long integer
27	integer, parameter, private :: llong_t = selected_int_kind(16) ! llong integer
28	integer, parameter, public :: i_llong = max( llong_t, i_long )
29
30	! Expected 8-bit byte sizes of the integer kinds
31	integer, parameter, public :: num_bytes_for_i_byte = 1
32	integer, parameter, public :: num_bytes_for_i_short = 2
33	integer, parameter, public :: num_bytes_for_i_long = 4
34	integer, parameter, public :: num_bytes_for_i_llong = 8
35
36	! Define arrays for default definition
37	integer, parameter, private :: num_i_kinds = 4
38	integer, parameter, dimension( num_i_kinds ), private :: integer_types = (/ &
39	i_byte, i_short, i_long, i_llong /)
40	integer, parameter, dimension( num_i_kinds ), private :: integer_byte_sizes = (/ &
41	num_bytes_for_i_byte, num_bytes_for_i_short, &
42	num_bytes_for_i_long, num_bytes_for_i_llong /)
43
44	! Default values
45	! ** CHANGE THE FOLLOWING TO CHANGE THE DEFAULT INTEGER TYPE KIND *
46	integer, parameter, private :: default_integer = 2 ! 1=byte,
47	! 2=short,
48	! 3=long,
49	! 4=llong
50	integer, parameter, public :: i_kind = integer_types( default_integer )
51	integer, parameter, public :: num_bytes_for_i_kind = &
52	integer_byte_sizes( default_integer )
53
54
55	! Real definitions below
56
57	! Real types
58	integer, parameter, public :: r_single = selected_real_kind(6) ! single precision
59	integer, parameter, public :: r_double = selected_real_kind(15) ! double precision
60	integer, parameter, private :: quad_t = selected_real_kind(20) ! quad precision
61	integer, parameter, public :: r_quad = max( quad_t, r_double )
62
63	! Expected 8-bit byte sizes of the real kinds
64	integer, parameter, public :: num_bytes_for_r_single = 4
65	integer, parameter, public :: num_bytes_for_r_double = 8
66	integer, parameter, public :: num_bytes_for_r_quad = 16
67
68	! Define arrays for default definition
69	integer, parameter, private :: num_r_kinds = 3
70	integer, parameter, dimension( num_r_kinds ), private :: real_kinds = (/ &
71	r_single, r_double, r_quad /)
72	integer, parameter, dimension( num_r_kinds ), private :: real_byte_sizes = (/ &
73	num_bytes_for_r_single, num_bytes_for_r_double, &
74	num_bytes_for_r_quad /)
75
76	! Default values
77	! ** CHANGE THE FOLLOWING TO CHANGE THE DEFAULT REAL TYPE KIND *
78	integer, parameter, private :: default_real = 2 ! 1=single,
79	! 2=double,
80	!! END FROM MODULE_KINDS
81
82	! Input 3D meteorological fields.
83
84	REAL , DIMENSION(:,:,:) , ALLOCATABLE :: u_input , v_input , &
85	q_input , t_input
86
87	! Input 3D LSM fields.
88
89	REAL , DIMENSION(:,:,:) , ALLOCATABLE :: landuse_frac_input , &
90	soil_top_cat_input , &
91	soil_bot_cat_input
92
93	REAL, ALLOCATABLE:: htm_in(:,:,:),vtm_in(:,:,:)
94
95	! Input 2D surface fields.
96
97	REAL , DIMENSION(:,:) , ALLOCATABLE,save :: soilt010_input , soilt040_input , &
98	soilt100_input , soilt200_input , &
99	soilm010_input , soilm040_input , &
100	soilm100_input , soilm200_input , &
101	psfc_in,pmsl
102
103	REAL , DIMENSION(:,:) , ALLOCATABLE :: lat_wind, lon_wind
104
105	REAL , DIMENSION(:) , ALLOCATABLE :: DETA_in, AETA_in, ETAX_in
106	REAL , DIMENSION(:) , ALLOCATABLE :: DETA1_in, AETA1_in, ETA1_in
107	REAL , DIMENSION(:) , ALLOCATABLE :: DETA2_in, AETA2_in, ETA2_in, DFL_in
108
109	REAL , DIMENSION(:,:,:), ALLOCATABLE,save :: st_inputx , sm_inputx, sw_inputx
110
111	! Local input arrays
112
113	REAL,DIMENSION(:,:),ALLOCATABLE :: dum2d
114	INTEGER,DIMENSION(:,:),ALLOCATABLE :: idum2d
115	REAL,DIMENSION(:,:,:),ALLOCATABLE :: dum3d
116
117	LOGICAL , SAVE :: first_time_in = .TRUE.
118
119	INTEGER :: flag_soilt010 , flag_soilt100 , flag_soilt200 , &
120	flag_soilm010 , flag_soilm100 , flag_soilm200
121
122	! Some constants to allow simple dimensions in the defined types
123	! given below.
124
125	INTEGER, PARAMETER :: var_maxdims = 5
126	INTEGER, PARAMETER :: max_staggers_xy_new = 4
127	INTEGER, PARAMETER :: max_staggers_xy_old = 3
128	INTEGER, PARAMETER :: max_staggers_z = 2
129	INTEGER, PARAMETER :: max_standard_lats = 4
130	INTEGER, PARAMETER :: max_standard_lons = 4
131	INTEGER, PARAMETER :: max_fg_variables = 200
132	INTEGER, PARAMETER :: max_vertical_levels = 2000
133
134	! This module defines the items needed for the WRF metadata
135	! which is broken up into three levels:
136	! Global metadata: Those things which apply to the
137	! entire simulation that are
138	! independent of time, domain, or
139	! variable
140	!
141	! Domain metadata: Those things which apply to
142	! a single domain (this may
143	! or may not be time dependent)
144	!
145	! Variable metadata: Those things which apply to
146	! a specific variable at a
147	! specific time
148	!
149	! The variable names and definitions can be
150	! found in the wrf_metadata spec, which is still
151	! a living document as coding goes on. The names
152	! may not match exactly, but you should be able
153	! to figure things out.
154	!
155
156	TYPE wrf_var_metadata
157	CHARACTER (LEN=8) :: name
158	CHARACTER (LEN=16) :: units
159	CHARACTER (LEN=80) :: description
160	INTEGER :: domain_id
161	INTEGER :: ndim
162	INTEGER :: dim_val (var_maxdims)
163	CHARACTER(LEN=4) :: dim_desc (var_maxdims)
164	INTEGER :: start_index(var_maxdims)
165	INTEGER :: stop_index(var_maxdims)
166	INTEGER :: h_stagger_index
167	INTEGER :: v_stagger_index
168	CHARACTER(LEN=8) :: array_order
169	CHARACTER(LEN=4) :: field_type
170	CHARACTER(LEN=8) :: field_source_prog
171	CHARACTER(LEN=80) :: source_desc
172	CHARACTER(LEN=8) :: field_time_type
173	INTEGER :: vt_date_start
174	REAL :: vt_time_start
175	INTEGER :: vt_date_stop
176	REAL :: vt_time_stop
177	END TYPE wrf_var_metadata
178
179	TYPE(wrf_var_metadata) :: var_meta , var_info
180
181	TYPE wrf_domain_metadata
182	INTEGER :: id
183	INTEGER :: parent_id
184	CHARACTER(LEN=8) :: dyn_init_src
185	CHARACTER(LEN=8) :: static_init_src
186	INTEGER :: vt_date
187	REAL :: vt_time
188	INTEGER :: origin_parent_x
189	INTEGER :: origin_parent_y
190	INTEGER :: ratio_to_parent
191	REAL :: delta_x
192	REAL :: delta_y
193	REAL :: top_level
194	INTEGER :: origin_parent_z
195	REAL :: corner_lats_new(4,max_staggers_xy_new)
196	REAL :: corner_lons_new(4,max_staggers_xy_new)
197	REAL :: corner_lats_old(4,max_staggers_xy_old)
198	REAL :: corner_lons_old(4,max_staggers_xy_old)
199	INTEGER :: xdim
200	INTEGER :: ydim
201	INTEGER :: zdim
202	END TYPE wrf_domain_metadata
203	TYPE(wrf_domain_metadata) :: dom_meta
204
205	TYPE wrf_global_metadata
206	CHARACTER(LEN=80) :: simulation_name
207	CHARACTER(LEN=80) :: user_desc
208	INTEGER :: si_version
209	INTEGER :: analysis_version
210	INTEGER :: wrf_version
211	INTEGER :: post_version
212	CHARACTER(LEN=32) :: map_projection
213	REAL :: moad_known_lat
214	REAL :: moad_known_lon
215	CHARACTER(LEN=8) :: moad_known_loc
216	REAL :: moad_stand_lats(max_standard_lats)
217	REAL :: moad_stand_lons(max_standard_lons)
218	REAL :: moad_delta_x
219	REAL :: moad_delta_y
220	CHARACTER(LEN=4) :: horiz_stagger_type
221	INTEGER :: num_stagger_xy
222	REAL :: stagger_dir_x_new(max_staggers_xy_new)
223	REAL :: stagger_dir_y_new(max_staggers_xy_new)
224	REAL :: stagger_dir_x_old(max_staggers_xy_old)
225	REAL :: stagger_dir_y_old(max_staggers_xy_old)
226	INTEGER :: num_stagger_z
227	REAL :: stagger_dir_z(max_staggers_z)
228	CHARACTER(LEN=8) :: vertical_coord
229	INTEGER :: num_domains
230	INTEGER :: init_date
231	REAL :: init_time
232	INTEGER :: end_date
233	REAL :: end_time
234	CHARACTER(LEN=4) :: lu_source
235	INTEGER :: lu_water
236	INTEGER :: lu_ice
237	END TYPE wrf_global_metadata
238	TYPE(wrf_global_metadata) :: global_meta
239
240	CONTAINS
241
242	SUBROUTINE read_si ( grid, file_date_string )
243
244	USE module_soil_pre
245	USE module_domain
246
247	IMPLICIT NONE
248
249	TYPE(domain) , INTENT(INOUT) :: grid
250	CHARACTER (LEN=19) , INTENT(IN) :: file_date_string
251
252	INTEGER :: ids,ide,jds,jde,kds,kde &
253	,ims,ime,jms,jme,kms,kme &
254	,its,ite,jts,jte,kts,kte
255
256	INTEGER :: i , j , k , loop, IMAX, JMAX
257
258	REAL :: dummy
259
260	CHARACTER (LEN= 8) :: dummy_char
261	CHARACTER (LEN=256) :: dummy_char_256
262
263	INTEGER :: ok , map_proj , ok_open
264	REAL :: pt
265	INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
266
267	SELECT CASE ( model_data_order )
268	CASE ( DATA_ORDER_ZXY )
269	kds = grid%sd31 ; kde = grid%ed31 ;
270	ids = grid%sd32 ; ide = grid%ed32 ;
271	jds = grid%sd33 ; jde = grid%ed33 ;
272
273	kms = grid%sm31 ; kme = grid%em31 ;
274	ims = grid%sm32 ; ime = grid%em32 ;
275	jms = grid%sm33 ; jme = grid%em33 ;
276
277	kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch
278	its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch
279	jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
280
281	CASE ( DATA_ORDER_XYZ )
282	ids = grid%sd31 ; ide = grid%ed31 ;
283	jds = grid%sd32 ; jde = grid%ed32 ;
284	kds = grid%sd33 ; kde = grid%ed33 ;
285
286	ims = grid%sm31 ; ime = grid%em31 ;
287	jms = grid%sm32 ; jme = grid%em32 ;
288	kms = grid%sm33 ; kme = grid%em33 ;
289
290	its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
291	jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch
292	kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch
293
294	CASE ( DATA_ORDER_XZY )
295	ids = grid%sd31 ; ide = grid%ed31 ;
296	kds = grid%sd32 ; kde = grid%ed32 ;
297	jds = grid%sd33 ; jde = grid%ed33 ;
298
299	ims = grid%sm31 ; ime = grid%em31 ;
300	kms = grid%sm32 ; kme = grid%em32 ;
301	jms = grid%sm33 ; jme = grid%em33 ;
302
303	its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
304	kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch
305	jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
306
307	END SELECT
308
309	! Initialize what soil temperature and moisture is available.
310
311	! write(0,*) 'dum3d I allocs: ', ids,ide-1
312	! write(0,*) 'dum3d J allocs: ', jds,jde-1
313	! write(0,*) 'dum3d K allocs: ', kds,kde-1
314
315	flag_st000010 = 0
316	flag_st010040 = 0
317	flag_st040100 = 0
318	flag_st100200 = 0
319	flag_sm000010 = 0
320	flag_sm010040 = 0
321	flag_sm040100 = 0
322	flag_sm100200 = 0
323	flag_st010200 = 0
324	flag_sm010200 = 0
325
326	flag_soilt010 = 0
327	flag_soilt040 = 0
328	flag_soilt100 = 0
329	flag_soilt200 = 0
330	flag_soilm010 = 0
331	flag_soilm040 = 0
332	flag_soilm100 = 0
333	flag_soilm200 = 0
334
335	flag_sst = 0
336	flag_toposoil = 0
337
338	! How many soil levels have we found? Well, right now, none.
339
340	num_st_levels_input = 0
341	num_sm_levels_input = 0
342	st_levels_input = -1
343	sm_levels_input = -1
344
345	! Get the space for the data if this is the first time here.
346
347	write(6,*) 'enter read_si...first_time_in:: ', first_time_in
348
349	IF ( first_time_in ) THEN
350
351	CLOSE(12)
352	OPEN ( FILE = 'real_input_nm.global.metadata' , &
353	UNIT = 12 , &
354	STATUS = 'OLD' , &
355	ACCESS = 'SEQUENTIAL' , &
356	FORM = 'UNFORMATTED' , &
357	IOSTAT = ok_open )
358
359	IF ( ok_open .NE. 0 ) THEN
360	PRINT '(A)','You asked for WRF SI data, but no real_input_nm.global.metadata file exists.'
361	STOP 'No_real_input_nm.global.metadata_exists'
362	END IF
363
364	READ(12) global_meta%simulation_name, global_meta%user_desc, &
365	global_meta%si_version, global_meta%analysis_version, &
366	global_meta%wrf_version, global_meta%post_version
367
368	REWIND (12)
369
370	IF ( global_meta%si_version .EQ. 1 ) THEN
371	READ(12) global_meta%simulation_name, global_meta%user_desc, &
372	global_meta%si_version, global_meta%analysis_version, &
373	global_meta%wrf_version, global_meta%post_version, &
374	global_meta%map_projection, global_meta%moad_known_lat, &
375	global_meta%moad_known_lon, global_meta%moad_known_loc, &
376	global_meta%moad_stand_lats, global_meta%moad_stand_lons, &
377	global_meta%moad_delta_x, global_meta%moad_delta_y, &
378	global_meta%horiz_stagger_type, global_meta%num_stagger_xy, &
379	global_meta%stagger_dir_x_old, global_meta%stagger_dir_y_old, &
380	global_meta%num_stagger_z, global_meta%stagger_dir_z, &
381	global_meta%vertical_coord, global_meta%num_domains, &
382	global_meta%init_date, global_meta%init_time, &
383	global_meta%end_date, global_meta%end_time
384	ELSE IF ( global_meta%si_version .EQ. 2 ) THEN
385	READ(12) global_meta%simulation_name, global_meta%user_desc, &
386	global_meta%si_version, global_meta%analysis_version, &
387	global_meta%wrf_version, global_meta%post_version, &
388	global_meta%map_projection, global_meta%moad_known_lat, &
389	global_meta%moad_known_lon, global_meta%moad_known_loc, &
390	global_meta%moad_stand_lats, global_meta%moad_stand_lons, &
391	global_meta%moad_delta_x, global_meta%moad_delta_y, &
392	global_meta%horiz_stagger_type, global_meta%num_stagger_xy, &
393	global_meta%stagger_dir_x_new, global_meta%stagger_dir_y_new, &
394	global_meta%num_stagger_z, global_meta%stagger_dir_z, &
395	global_meta%vertical_coord, global_meta%num_domains, &
396	global_meta%init_date, global_meta%init_time, &
397	global_meta%end_date, global_meta%end_time , &
398	global_meta%lu_source, global_meta%lu_water, global_meta%lu_ice
399	END IF
400	CLOSE (12)
401
402	print *,'GLOBAL METADATA'
403	print *,'global_meta%simulation_name', global_meta%simulation_name
404	print *,'global_meta%user_desc', global_meta%user_desc
405	print *,'global_meta%user_desc', global_meta%user_desc
406	print *,'global_meta%si_version', global_meta%si_version
407	print *,'global_meta%analysis_version', global_meta%analysis_version
408	print *,'global_meta%wrf_version', global_meta%wrf_version
409	print *,'global_meta%post_version', global_meta%post_version
410	print *,'global_meta%map_projection', global_meta%map_projection
411	print *,'global_meta%moad_known_lat', global_meta%moad_known_lat
412	print *,'global_meta%moad_known_lon', global_meta%moad_known_lon
413	print *,'global_meta%moad_known_loc', global_meta%moad_known_loc
414	print *,'global_meta%moad_stand_lats', global_meta%moad_stand_lats
415	print *,'global_meta%moad_stand_lons', global_meta%moad_stand_lons
416	print *,'global_meta%moad_delta_x', global_meta%moad_delta_x
417	print *,'global_meta%moad_delta_y', global_meta%moad_delta_y
418	print *,'global_meta%horiz_stagger_type', global_meta%horiz_stagger_type
419	print *,'global_meta%num_stagger_xy', global_meta%num_stagger_xy
420	IF ( global_meta%si_version .EQ. 1 ) THEN
421	print *,'global_meta%stagger_dir_x', global_meta%stagger_dir_x_old
422	print *,'global_meta%stagger_dir_y', global_meta%stagger_dir_y_old
423	ELSE IF ( global_meta%si_version .EQ. 2 ) THEN
424	print *,'global_meta%stagger_dir_x', global_meta%stagger_dir_x_new
425	print *,'global_meta%stagger_dir_y', global_meta%stagger_dir_y_new
426	END IF
427	print *,'global_meta%num_stagger_z', global_meta%num_stagger_z
428	print *,'global_meta%stagger_dir_z', global_meta%stagger_dir_z
429	print *,'global_meta%vertical_coord', global_meta%vertical_coord
430	print *,'global_meta%num_domains', global_meta%num_domains
431	print *,'global_meta%init_date', global_meta%init_date
432	print *,'global_meta%init_time', global_meta%init_time
433	print *,'global_meta%end_date', global_meta%end_date
434	print *,'global_meta%end_time', global_meta%end_time
435	IF ( global_meta%si_version .EQ. 2 ) THEN
436	print *,'global_meta%lu_source', global_meta%lu_source
437	print *,'global_meta%lu_water', global_meta%lu_water
438	print *,'global_meta%lu_ice', global_meta%lu_ice
439	END IF
440	print *,' '
441
442	! 1D - this is the definition of the vertical coordinate.
443
444	IF (.NOT. ALLOCATED (DETA_in)) ALLOCATE(DETA_in(kds:kde-1))
445	IF (.NOT. ALLOCATED (AETA_in)) ALLOCATE(AETA_in(kds:kde-1))
446	IF (.NOT. ALLOCATED (ETAX_in)) ALLOCATE(ETAX_in(kds:kde))
447
448	IF (.NOT. ALLOCATED (DETA1_in)) ALLOCATE(DETA1_in(kds:kde-1))
449	IF (.NOT. ALLOCATED (AETA1_in)) ALLOCATE(AETA1_in(kds:kde-1))
450	IF (.NOT. ALLOCATED (ETA1_in)) ALLOCATE(ETA1_in(kds:kde))
451
452	IF (.NOT. ALLOCATED (DETA2_in)) ALLOCATE(DETA2_in(kds:kde-1))
453	IF (.NOT. ALLOCATED (AETA2_in)) ALLOCATE(AETA2_in(kds:kde-1))
454	IF (.NOT. ALLOCATED (ETA2_in)) ALLOCATE(ETA2_in(kds:kde))
455
456	IF (.NOT. ALLOCATED (DFL_in)) ALLOCATE(DFL_in(kds:kde))
457
458	! 3D met
459
460	IF (.NOT. ALLOCATED (u_input) ) ALLOCATE ( u_input(its:ite,jts:jte,kts:kte) )
461	IF (.NOT. ALLOCATED (v_input) ) ALLOCATE ( v_input(its:ite,jts:jte,kts:kte) )
462	IF (.NOT. ALLOCATED (q_input) ) ALLOCATE ( q_input(its:ite,jts:jte,kts:kte) )
463	IF (.NOT. ALLOCATED (t_input) ) ALLOCATE ( t_input(its:ite,jts:jte,kts:kte) )
464	IF (.NOT. ALLOCATED (htm_in) ) ALLOCATE ( htm_in(its:ite,jts:jte,kts:kte) )
465	IF (.NOT. ALLOCATED (vtm_in) ) ALLOCATE ( vtm_in(its:ite,jts:jte,kts:kte) )
466
467	! 2D pressure fields
468
469	IF (.NOT. ALLOCATED (pmsl) ) ALLOCATE ( pmsl(its:ite,jts:jte) )
470	IF (.NOT. ALLOCATED (psfc_in) ) ALLOCATE ( psfc_in(its:ite,jts:jte) )
471
472	! 2D - for LSM, these are computed from the categorical precentage values.
473
474	! 2D - for LSM, the various soil temperature and moisture levels that are available.
475
476	IF (.NOT. ALLOCATED (st_inputx)) ALLOCATE (st_inputx(its:ite,jts:jte,num_st_levels_alloc))
477	IF (.NOT. ALLOCATED (sm_inputx)) ALLOCATE (sm_inputx(its:ite,jts:jte,num_st_levels_alloc))
478	IF (.NOT. ALLOCATED (sw_inputx)) ALLOCATE (sw_inputx(its:ite,jts:jte,num_st_levels_alloc))
479
480	IF (.NOT. ALLOCATED (soilt010_input) ) ALLOCATE ( soilt010_input(its:ite,jts:jte) )
481	IF (.NOT. ALLOCATED (soilt040_input) ) ALLOCATE ( soilt040_input(its:ite,jts:jte) )
482	IF (.NOT. ALLOCATED (soilt100_input) ) ALLOCATE ( soilt100_input(its:ite,jts:jte) )
483	IF (.NOT. ALLOCATED (soilt200_input) ) ALLOCATE ( soilt200_input(its:ite,jts:jte) )
484	IF (.NOT. ALLOCATED (soilm010_input) ) ALLOCATE ( soilm010_input(its:ite,jts:jte) )
485	IF (.NOT. ALLOCATED (soilm040_input) ) ALLOCATE ( soilm040_input(its:ite,jts:jte) )
486	IF (.NOT. ALLOCATED (soilm100_input) ) ALLOCATE ( soilm100_input(its:ite,jts:jte) )
487	IF (.NOT. ALLOCATED (soilm200_input) ) ALLOCATE ( soilm200_input(its:ite,jts:jte) )
488
489	IF (.NOT. ALLOCATED (lat_wind) ) ALLOCATE (lat_wind(its:ite,jts:jte))
490	IF (.NOT. ALLOCATED (lon_wind) ) ALLOCATE (lon_wind(its:ite,jts:jte))
491
492	! Local arrays
493	IF (.NOT. ALLOCATED (dum2d) ) ALLOCATE (dum2d(IDS:IDE-1,JDS:JDE-1))
494	IF (.NOT. ALLOCATED (idum2d) ) ALLOCATE (idum2d(IDS:IDE-1,JDS:JDE-1))
495	IF (.NOT. ALLOCATED (dum3d) ) ALLOCATE (dum3d(IDS:IDE-1,JDS:JDE-1,KDS:KDE-1))
496
497
498	END IF
499
500	CLOSE(13)
501
502	write(6,*) 'file_date_string: ', file_date_string
503	write(6,*) 'opening real_input_nm.d01.'//file_date_string//' as unit 13'
504	OPEN ( FILE = 'real_input_nm.d01.'//file_date_string , &
505	UNIT = 13 , &
506	STATUS = 'OLD' , &
507	ACCESS = 'SEQUENTIAL' , &
508	FORM = 'UNFORMATTED' )
509
510	IF ( global_meta%si_version .EQ. 1 ) THEN
511	READ (13) dom_meta%id,dom_meta%parent_id,dom_meta%dyn_init_src,&
512	dom_meta%static_init_src, dom_meta%vt_date, dom_meta%vt_time, &
513	dom_meta%origin_parent_x, dom_meta%origin_parent_y, &
514	dom_meta%ratio_to_parent, dom_meta%delta_x, dom_meta%delta_y, &
515	dom_meta%top_level, dom_meta%origin_parent_z, &
516	dom_meta%corner_lats_old, dom_meta%corner_lons_old, dom_meta%xdim, &
517	dom_meta%ydim, dom_meta%zdim
518	ELSE IF ( global_meta%si_version .EQ. 2 ) THEN
519	READ (13) dom_meta%id,dom_meta%parent_id,dom_meta%dyn_init_src,&
520	dom_meta%static_init_src, dom_meta%vt_date, dom_meta%vt_time, &
521	dom_meta%origin_parent_x, dom_meta%origin_parent_y, &
522	dom_meta%ratio_to_parent, dom_meta%delta_x, dom_meta%delta_y, &
523	dom_meta%top_level, dom_meta%origin_parent_z, &
524	dom_meta%corner_lats_new, dom_meta%corner_lons_new, dom_meta%xdim, &
525	dom_meta%ydim, dom_meta%zdim
526	END IF
527
528	print *,'DOMAIN METADATA'
529	print *,'dom_meta%id=', dom_meta%id
530	print *,'dom_meta%parent_id=', dom_meta%parent_id
531	print *,'dom_meta%dyn_init_src=', dom_meta%dyn_init_src
532	print *,'dom_meta%static_init_src=', dom_meta%static_init_src
533	print *,'dom_meta%vt_date=', dom_meta%vt_date
534	print *,'dom_meta%vt_time=', dom_meta%vt_time
535	print *,'dom_meta%origin_parent_x=', dom_meta%origin_parent_x
536	print *,'dom_meta%origin_parent_y=', dom_meta%origin_parent_y
537	print *,'dom_meta%ratio_to_parent=', dom_meta%ratio_to_parent
538	print *,'dom_meta%delta_x=', dom_meta%delta_x
539	print *,'dom_meta%delta_y=', dom_meta%delta_y
540	print *,'dom_meta%top_level=', dom_meta%top_level
541	print *,'dom_meta%origin_parent_z=', dom_meta%origin_parent_z
542	IF ( global_meta%si_version .EQ. 1 ) THEN
543	print *,'dom_meta%corner_lats=', dom_meta%corner_lats_old
544	print *,'dom_meta%corner_lons=', dom_meta%corner_lons_old
545	ELSE IF ( global_meta%si_version .EQ. 2 ) THEN
546	print *,'dom_meta%corner_lats=', dom_meta%corner_lats_new
547	print *,'dom_meta%corner_lons=', dom_meta%corner_lons_new
548	END IF
549	print *,'dom_meta%xdim=', dom_meta%xdim
550	print *,'dom_meta%ydim=', dom_meta%ydim
551	print *,'dom_meta%zdim=', dom_meta%zdim
552	print *,' '
553
554	! A simple domain size test.
555
556
557	!! relax constraint, as model namelist has +1 for i and j, while
558	!! si data has true dimensions
559
560	IF ( abs(dom_meta%xdim - (ide-1)) .gt. 1 &
561	.OR. abs(dom_meta%ydim - (jde-1)) .gt. 1 &
562	.OR. abs(dom_meta%zdim - (kde-1)) .gt. 1) THEN
563	PRINT '(A)','Namelist does not match the input data.'
564	PRINT '(A,3I5,A)','Namelist dimensions =',ide-1,jde-1,kde-1,'.'
565	PRINT '(A,3I5,A)','Input data dimensions =',dom_meta%xdim,dom_meta%ydim,dom_meta%zdim,'.'
566	STOP 'Wrong_data_size'
567	END IF
568
569	! How about the grid distance? Is it the same as in the namelist?
570
571	IF ( global_meta%si_version .EQ. 1 ) THEN
572	CALL nl_set_cen_lat ( grid%id , ( dom_meta%corner_lats_old(1,1) + dom_meta%corner_lats_old(2,1) + &
573	dom_meta%corner_lats_old(3,1) + dom_meta%corner_lats_old(4,1) ) * 0.25 )
574	ELSE IF ( ( global_meta%si_version .EQ. 2 ) .AND. ( global_meta%moad_known_loc(1:6) .EQ. 'CENTER' ) ) THEN
575	CALL nl_set_cen_lat ( grid%id , global_meta%moad_known_lat )
576	ELSE IF ( global_meta%si_version .EQ. 2 ) THEN
577	CALL nl_set_cen_lat ( grid%id , ( dom_meta%corner_lats_new(1,1) + dom_meta%corner_lats_new(2,1) + &
578	dom_meta%corner_lats_new(3,1) + dom_meta%corner_lats_new(4,1) ) * 0.25 )
579	END IF
580
581
582	!!! might be trouble here
583
584	CALL nl_set_cen_lon ( grid%id , global_meta%moad_stand_lons(1) )
585	!!!!!
586	write(6,*) 'set_cen_lat... global_meta%moad_stand_lats(1): ', global_meta%moad_stand_lats(1)
587	CALL nl_set_cen_lat ( grid%id , global_meta%moad_stand_lats(1) )
588	!!!!!
589	CALL nl_set_truelat1 ( grid%id , global_meta%moad_stand_lats(1) )
590	CALL nl_set_truelat2 ( grid%id , global_meta%moad_stand_lats(2) )
591
592	pt = dom_meta%top_level
593
594	IF ( global_meta%map_projection(1:17) .EQ. 'LAMBERT CONFORMAL' ) THEN
595	map_proj = 1
596	ELSE IF ( global_meta%map_projection(1:19) .EQ. 'POLAR STEREOGRAPHIC' ) THEN
597	map_proj = 2
598	ELSE IF ( global_meta%map_projection(1: 8) .EQ. 'MERCATOR' ) THEN
599	map_proj = 3
600	ELSE IF ( global_meta%map_projection(1:14) .EQ. 'ROTATED LATLON' ) THEN
601	map_proj = 203 !?
602	ELSE
603	PRINT '(A,A,A)','Undefined map projection: ',TRIM(global_meta%map_projection(1:20)),'.'
604	STOP 'Undefined_map_proj_si'
605	END IF
606	CALL nl_set_map_proj ( grid%id , map_proj )
607
608	! write(0,*) 'global_meta%si_version: ', global_meta%si_version
609	! write(0,*) 'global_meta%lu_source: ', global_meta%lu_source
610	! write(0,*) 'global_meta%lu_water: ', global_meta%lu_water
611	IF ( global_meta%si_version .EQ. 1 ) THEN
612	CALL nl_set_mminlu (grid%id, 'USGS' )
613	CALL nl_set_iswater (grid%id, 16 )
614	ELSE IF ( global_meta%si_version .EQ. 2 ) THEN
615	CALL nl_set_mminlu ( grid%id, global_meta%lu_source )
616	CALL nl_set_iswater (grid%id, global_meta%lu_water )
617	CALL nl_set_isice (grid%id, global_meta%lu_ice )
618	END IF
619
620	CALL nl_set_gmt (grid%id, dom_meta%vt_time / 3600. )
621	CALL nl_set_julyr (grid%id, dom_meta%vt_date / 1000 )
622	CALL nl_set_julday (grid%id, dom_meta%vt_date - ( dom_meta%vt_date / 1000 ) * 1000 )
623
624	write(6,*) 'start reading from unit 13'
625	read_all_the_data : DO
626
627	READ (13,IOSTAT=OK) var_info%name, var_info%units, &
628	var_info%description, var_info%domain_id, var_info%ndim, &
629	var_info%dim_val, var_info%dim_desc, var_info%start_index, &
630	var_info%stop_index, var_info%h_stagger_index, var_info%v_stagger_index,&
631	var_info%array_order, var_info%field_type, var_info%field_source_prog, &
632	var_info%source_desc, var_info%field_time_type, var_info%vt_date_start, &
633	var_info%vt_time_start, var_info%vt_date_stop, var_info%vt_time_stop
634
635	IF ( OK .NE. 0 ) THEN
636	PRINT '(A,A,A)','End of file found for real_input_nm.d01.',file_date_string,'.'
637	EXIT read_all_the_data
638	END IF
639
640	! print *,'VARIABLE METADATA'
641	PRINT '(A,A)','var_info%name=', var_info%name
642	! print *,'var_info%units=', var_info%units
643	! print *,'var_info%description=', var_info%description
644	! print *,'var_info%domain_id=', var_info%domain_id
645	! print *,'var_info%ndim=', var_info%ndim
646	! print *,'var_info%dim_val=', var_info%dim_val
647	! print *,'var_info%dim_desc=', var_info%dim_desc
648	! print *,'var_info%start_index=', var_info%start_index
649	! print *,'var_info%stop_index=', var_info%stop_index
650	! print *,'var_info%h_stagger_index=', var_info%h_stagger_index
651	! print *,'var_info%v_stagger_index=', var_info%v_stagger_index
652	! print *,'var_info%array_order=', var_info%array_order
653	! print *,'var_info%field_type=', var_info%field_type
654	! print *,'var_info%field_source_prog=', var_info%field_source_prog
655	! print *,'var_info%source_desc=', var_info%source_desc
656	! print *,'var_info%field_time_type=', var_info%field_time_type
657	! print *,'var_info%vt_date_start=', var_info%vt_date_start
658	! print *,'var_info%vt_time_start=', var_info%vt_time_start
659	! print *,'var_info%vt_date_stop=', var_info%vt_date_stop
660	! print *,'var_info%vt_time_stop=', var_info%vt_time_stop
661
662	JMAX=min(JDE-1,JTE)
663	IMAX=min(IDE-1,ITE)
664	! 3D meteorological fields.
665
666	! write(0,*)' read_si var_info%name=',var_info%name(1:8)
667
668	IF ( var_info%name(1:8) .EQ. 'T ' ) THEN
669	READ (13) dum3d
670	do k=kts,kte-1
671	do j=jts,JMAX
672	do i=its,IMAX
673	t_input(i,j,k)=dum3d(i,j,k)
674	enddo
675	enddo
676	enddo
677
678	ELSE IF ( var_info%name(1:8) .EQ. 'U ' ) THEN
679	READ (13) dum3d
680	do k=kts,kte-1
681	do j=jts,JMAX
682	do i=its,IMAX
683	u_input(i,j,k)=dum3d(i,j,k)
684	enddo
685	enddo
686	enddo
687
688	ELSE IF ( var_info%name(1:8) .EQ. 'V ' ) THEN
689	READ (13) dum3d
690	do k=kts,kte-1
691	do j=jts,JMAX
692	do i=its,IMAX
693	v_input(i,j,k)=dum3d(i,j,k)
694	enddo
695	enddo
696	enddo
697
698	ELSE IF ( var_info%name(1:8) .EQ. 'Q ' ) THEN
699	READ (13) dum3d
700	do k=kts,kte-1
701	do j=jts,JMAX
702	do i=its,IMAX
703	q_input(i,j,k)=dum3d(i,j,k)
704	enddo
705	enddo
706	enddo
707
708	! 3D LSM fields. Don't know the 3rd dimension until we read it in.
709
710	ELSE IF ( var_info%name(1:8) .EQ. 'LANDUSEF' ) THEN
711	IF ( ( first_time_in ) .AND. ( .NOT. ALLOCATED ( landuse_frac_input) ) ) THEN
712	ALLOCATE (landuse_frac_input(its:ite,jts:jte,var_info%dim_val(3)) )
713	END IF
714	READ (13) (((dum3d(i,j,k),i=ids,ide-1),j=jds,jde-1),k=1,var_info%dim_val(3))
715	do k=1,var_info%dim_val(3)
716	do j=jts,JMAX
717	do i=its,IMAX
718	landuse_frac_input(i,j,k)=dum3d(i,j,k)
719	enddo
720	enddo
721	enddo
722	ELSE IF ( var_info%name(1:8) .EQ. 'SOILCTOP' ) THEN
723	IF ( ( first_time_in ) .AND. ( .NOT. ALLOCATED ( soil_top_cat_input) ) ) THEN
724	ALLOCATE (soil_top_cat_input(its:ite,jts:jte,var_info%dim_val(3)) )
725	END IF
726	READ (13) (((dum3d(i,j,k),i=ids,ide-1),j=jds,jde-1),k=1,var_info%dim_val(3))
727	do k=1,var_info%dim_val(3)
728	do j=jts,JMAX
729	do i=its,IMAX
730	soil_top_cat_input(i,j,k)=dum3d(i,j,k)
731	enddo
732	enddo
733	enddo
734	ELSE IF ( var_info%name(1:8) .EQ. 'SOILCBOT' ) THEN
735	IF ( ( first_time_in ) .AND. ( .NOT. ALLOCATED ( soil_bot_cat_input) ) ) THEN
736	ALLOCATE (soil_bot_cat_input(its:ite,jts:jte,var_info%dim_val(3)) )
737	END IF
738	READ (13) (((dum3d(i,j,k),i=ids,ide-1),j=jds,jde-1),k=1,var_info%dim_val(3))
739	do k=1,var_info%dim_val(3)
740	do j=jts,JMAX
741	do i=its,IMAX
742	soil_bot_cat_input(i,j,k)=dum3d(i,j,k)
743	enddo
744	enddo
745	enddo
746
747	! 2D dry pressure minus ptop.
748
749	ELSE IF ( var_info%name(1:8) .EQ. 'PD ' ) THEN
750	READ (13) dum2d
751	do j=jts,JMAX
752	do i=its,IMAX
753	grid%pd(i,j)=dum2d(i,j)
754	enddo
755	enddo
756	ELSE IF ( var_info%name(1:8) .EQ. 'PSFC ' ) THEN
757	READ (13) dum2d
758	do j=jts,JMAX
759	do i=its,IMAX
760	psfc_in(i,j)=dum2d(i,j)
761	enddo
762	enddo
763	ELSE IF ( var_info%name(1:8) .EQ. 'PMSL ' ) THEN
764	READ (13) dum2d
765	do j=jts,JMAX
766	do i=its,IMAX
767	pmsl(i,j)=dum2d(i,j)
768	enddo
769	enddo
770	ELSE IF ( var_info%name(1:8) .EQ. 'PDTOP ' ) THEN
771	READ (13) grid%pdtop
772
773	ELSE IF ( var_info%name(1:8) .EQ. 'PT ' ) THEN
774	READ (13) grid%pt
775
776	! 2D surface fields.
777
778	ELSE IF ( var_info%name(1:8) .eq. 'GLAT ' ) THEN
779	READ (13) dum2d
780	do j=jts,JMAX
781	do i=its,IMAX
782	grid%glat(i,j)=dum2d(i,j)
783	enddo
784	enddo
785	ELSE IF ( var_info%name(1:8) .eq. 'GLON ' ) THEN
786	READ (13) dum2d
787	do j=jts,JMAX
788	do i=its,IMAX
789	grid%glon(i,j)=dum2d(i,j)
790	enddo
791	enddo
792	ELSE IF ( var_info%name(1:8) .eq. 'LAT_V ' ) THEN
793	READ (13) dum2d
794	do j=jts,JMAX
795	do i=its,IMAX
796	lat_wind(i,j)=dum2d(i,j)
797	enddo
798	enddo
799	ELSE IF ( var_info%name(1:8) .eq. 'LON_V ' ) THEN
800	READ (13) dum2d
801	do j=jts,JMAX
802	do i=its,IMAX
803	lon_wind(i,j)=dum2d(i,j)
804	enddo
805	enddo
806
807	ELSE IF ( var_info%name(1:8) .EQ. 'ST000010' ) THEN
808	READ (13) dum2d
809	do j=jts,JMAX
810	do i=its,IMAX
811	grid%st000010(i,j)=dum2d(i,j)
812	enddo
813	enddo
814	flag_st000010 = 1
815	num_st_levels_input = num_st_levels_input + 1
816	st_levels_input(num_st_levels_input) = char2int2(var_info%name(3:8))
817	do j=jts,JMAX
818	do i=its,IMAX
819	st_inputx(I,J,num_st_levels_input + 1) = grid%st000010(i,j)
820	enddo
821	enddo
822
823	ELSE IF ( var_info%name(1:8) .EQ. 'ST010040' ) THEN
824	READ (13) dum2d
825	do j=jts,JMAX
826	do i=its,IMAX
827	grid%st010040(i,j)=dum2d(i,j)
828	enddo
829	enddo
830	flag_st010040 = 1
831	num_st_levels_input = num_st_levels_input + 1
832	st_levels_input(num_st_levels_input) = char2int2(var_info%name(3:8))
833	do j=jts,JMAX
834	do i=its,IMAX
835	st_inputx(I,J,num_st_levels_input + 1) = grid%st010040(i,j)
836	enddo
837	enddo
838
839	ELSE IF ( var_info%name(1:8) .EQ. 'ST040100' ) THEN
840	READ (13) dum2d
841	do j=jts,JMAX
842	do i=its,IMAX
843	grid%st040100(i,j)=dum2d(i,j)
844	enddo
845	enddo
846	flag_st040100 = 1
847	num_st_levels_input = num_st_levels_input + 1
848	st_levels_input(num_st_levels_input) = char2int2(var_info%name(3:8))
849	do j=jts,JMAX
850	do i=its,IMAX
851	st_inputx(I,J,num_st_levels_input + 1) = grid%st040100(i,j)
852	enddo
853	enddo
854
855	ELSE IF ( var_info%name(1:8) .EQ. 'ST100200' ) THEN
856	READ (13) dum2d
857	do j=jts,JMAX
858	do i=its,IMAX
859	grid%st100200(i,j)=dum2d(i,j)
860	enddo
861	enddo
862	flag_st100200 = 1
863	num_st_levels_input = num_st_levels_input + 1
864	st_levels_input(num_st_levels_input) = char2int2(var_info%name(3:8))
865	do j=jts,JMAX
866	do i=its,IMAX
867	st_inputx(I,J,num_st_levels_input + 1) = grid%st100200(i,j)
868	enddo
869	enddo
870
871	ELSE IF ( var_info%name(1:8) .EQ. 'ST010200' ) THEN
872	READ (13) dum2d
873	do j=jts,JMAX
874	do i=its,IMAX
875	grid%st010200(i,j)=dum2d(i,j)
876	enddo
877	enddo
878	flag_st010200 = 1
879	num_st_levels_input = num_st_levels_input + 1
880	st_levels_input(num_st_levels_input) = char2int2(var_info%name(3:8))
881	do j=jts,JMAX
882	do i=its,IMAX
883	st_inputx(I,J,num_st_levels_input + 1) = grid%st010200(i,j)
884	enddo
885	enddo
886
887	ELSE IF ( var_info%name(1:8) .EQ. 'SM000010' ) THEN
888	READ (13) dum2d
889	do j=jts,JMAX
890	do i=its,IMAX
891	grid%sm000010(i,j)=dum2d(i,j)
892	enddo
893	enddo
894	flag_sm000010 = 1
895	num_sm_levels_input = num_sm_levels_input + 1
896	sm_levels_input(num_sm_levels_input) = char2int2(var_info%name(3:8))
897	do j=jts,JMAX
898	do i=its,IMAX
899	sm_inputx(I,J,num_sm_levels_input + 1) = grid%sm000010(i,j)
900	enddo
901	enddo
902
903	ELSE IF ( var_info%name(1:8) .EQ. 'SM010040' ) THEN
904	READ (13) dum2d
905	do j=jts,JMAX
906	do i=its,IMAX
907	grid%sm010040(i,j)=dum2d(i,j)
908	enddo
909	enddo
910	flag_sm010040 = 1
911	num_sm_levels_input = num_sm_levels_input + 1
912	sm_levels_input(num_sm_levels_input) = char2int2(var_info%name(3:8))
913	do j=jts,JMAX
914	do i=its,IMAX
915	sm_inputx(I,J,num_sm_levels_input + 1) = grid%sm010040(i,j)
916	enddo
917	enddo
918
919	ELSE IF ( var_info%name(1:8) .EQ. 'SM040100' ) THEN
920	READ (13) dum2d
921	do j=jts,JMAX
922	do i=its,IMAX
923	grid%sm040100(i,j)=dum2d(i,j)
924	enddo
925	enddo
926	flag_sm040100 = 1
927	num_sm_levels_input = num_sm_levels_input + 1
928	sm_levels_input(num_sm_levels_input) = char2int2(var_info%name(3:8))
929	do j=jts,JMAX
930	do i=its,IMAX
931	sm_inputx(I,J,num_sm_levels_input + 1) = grid%sm040100(i,j)
932	enddo
933	enddo
934
935	ELSE IF ( var_info%name(1:8) .EQ. 'SM100200' ) THEN
936	READ (13) dum2d
937	do j=jts,JMAX
938	do i=its,IMAX
939	grid%sm100200(i,j)=dum2d(i,j)
940	enddo
941	enddo
942	flag_sm100200 = 1
943	num_sm_levels_input = num_sm_levels_input + 1
944	sm_levels_input(num_sm_levels_input) = char2int2(var_info%name(3:8))
945	do j=jts,JMAX
946	do i=its,IMAX
947	sm_inputx(I,J,num_sm_levels_input + 1) = grid%sm100200(i,j)
948	enddo
949	enddo
950
951	ELSE IF ( var_info%name(1:8) .EQ. 'SM010200' ) THEN
952	READ (13) dum2d
953	do j=jts,JMAX
954	do i=its,IMAX
955	grid%sm010200(i,j)=dum2d(i,j)
956	enddo
957	enddo
958	flag_sm010200 = 1
959	num_sm_levels_input = num_sm_levels_input + 1
960	sm_levels_input(num_sm_levels_input) = char2int2(var_info%name(3:8))
961	do j=jts,JMAX
962	do i=its,IMAX
963	sm_inputx(I,J,num_sm_levels_input + 1) = grid%sm010200(i,j)
964	enddo
965	enddo
966
967	ELSE IF ( var_info%name(1:8) .EQ. 'SOILT010' ) THEN
968	READ (13) dum2d
969	do j=jts,JMAX
970	do i=its,IMAX
971	soilt010_input(i,j)=dum2d(i,j)
972	enddo
973	enddo
974	flag_soilt010 = 1
975	num_st_levels_input = num_st_levels_input + 1
976	st_levels_input(num_st_levels_input) = char2int1(var_info%name(6:8))
977	!mp st_inputx(:,:,num_st_levels_input + 1) = soilt010_input
978	do j=jts,JMAX
979	do i=its,IMAX
980	st_inputx(I,J,num_st_levels_input + 1) = soilt010_input(I,J)
981	enddo
982	enddo
983	write(6,*) 'num_st_levels_input=',num_st_levels_input
984	ELSE IF ( var_info%name(1:8) .EQ. 'SOILT040' ) THEN
985	READ (13) dum2d
986	do j=jts,JMAX
987	do i=its,IMAX
988	soilt040_input(i,j)=dum2d(i,j)
989	enddo
990	enddo
991	flag_soilt040 = 1
992	num_st_levels_input = num_st_levels_input + 1
993	st_levels_input(num_st_levels_input) = char2int1(var_info%name(6:8))
994	!mp st_inputx(:,:,num_st_levels_input + 1) = soilt040_input
995	do j=jts,JMAX
996	do i=its,IMAX
997	st_inputx(I,J,num_st_levels_input + 1) = soilt040_input(I,J)
998	enddo
999	enddo
1000	write(6,*) 'num_st_levels_input=',num_st_levels_input
1001	ELSE IF ( var_info%name(1:8) .EQ. 'SOILT100' ) THEN
1002	READ (13) dum2d
1003	do j=jts,JMAX
1004	do i=its,IMAX
1005	soilt100_input(i,j)=dum2d(i,j)
1006	enddo
1007	enddo
1008	flag_soilt100 = 1
1009	num_st_levels_input = num_st_levels_input + 1
1010	st_levels_input(num_st_levels_input) = char2int1(var_info%name(6:8))
1011	!mp st_inputx(:,:,num_st_levels_input + 1) = soilt100_input
1012	do j=jts,JMAX
1013	do i=its,IMAX
1014	st_inputx(I,J,num_st_levels_input + 1) = soilt100_input(I,J)
1015	enddo
1016	enddo
1017	write(6,*) 'num_st_levels_input=',num_st_levels_input
1018	ELSE IF ( var_info%name(1:8) .EQ. 'SOILT200' ) THEN
1019	READ (13) dum2d
1020	do j=jts,JMAX
1021	do i=its,IMAX
1022	soilt200_input(i,j)=dum2d(i,j)
1023	enddo
1024	enddo
1025	flag_soilt200 = 1
1026	num_st_levels_input = num_st_levels_input + 1
1027	st_levels_input(num_st_levels_input) = char2int1(var_info%name(6:8))
1028	!mp st_inputx(:,:,num_st_levels_input + 1) = soilt200_input
1029	do j=jts,JMAX
1030	do i=its,IMAX
1031	st_inputx(I,J,num_st_levels_input + 1) = soilt200_input(I,J)
1032	enddo
1033	enddo
1034	write(6,*) 'num_st_levels_input=',num_st_levels_input
1035	ELSE IF ( var_info%name(1:8) .EQ. 'SOILM010' ) THEN
1036	READ (13) dum2d
1037	do j=jts,JMAX
1038	do i=its,IMAX
1039	soilm010_input(i,j)=dum2d(i,j)
1040	enddo
1041	enddo
1042	flag_soilm010 = 1
1043	num_sm_levels_input = num_sm_levels_input + 1
1044	sm_levels_input(num_sm_levels_input) = char2int1(var_info%name(6:8))
1045	!mp sm_inputx(:,:,num_sm_levels_input + 1) = soilm010_input
1046	do j=jts,JMAX
1047	do i=its,IMAX
1048	sm_inputx(I,J,num_sm_levels_input + 1) = soilm010_input(I,J)
1049	enddo
1050	enddo
1051
1052	ELSE IF ( var_info%name(1:8) .EQ. 'SOILM040' ) THEN
1053	READ (13) dum2d
1054	do j=jts,JMAX
1055	do i=its,IMAX
1056	soilm040_input(i,j)=dum2d(i,j)
1057	enddo
1058	enddo
1059	flag_soilm040 = 1
1060	num_sm_levels_input = num_sm_levels_input + 1
1061	sm_levels_input(num_sm_levels_input) = char2int1(var_info%name(6:8))
1062	!mp sm_inputx(:,:,num_sm_levels_input + 1) = soilm040_input
1063	do j=jts,JMAX
1064	do i=its,IMAX
1065	sm_inputx(I,J,num_sm_levels_input + 1) = soilm040_input(I,J)
1066	enddo
1067	enddo
1068	ELSE IF ( var_info%name(1:8) .EQ. 'SOILM100' ) THEN
1069	READ (13) dum2d
1070	do j=jts,JMAX
1071	do i=its,IMAX
1072	soilm100_input(i,j)=dum2d(i,j)
1073	enddo
1074	enddo
1075	flag_soilm100 = 1
1076	num_sm_levels_input = num_sm_levels_input + 1
1077	sm_levels_input(num_sm_levels_input) = char2int1(var_info%name(6:8))
1078	!mp sm_inputx(:,:,num_sm_levels_input + 1) = soilm100_input
1079	do j=jts,JMAX
1080	do i=its,IMAX
1081	sm_inputx(I,J,num_sm_levels_input + 1) = soilm100_input(I,J)
1082	enddo
1083	enddo
1084
1085	ELSE IF ( var_info%name(1:8) .EQ. 'SOILM200' ) THEN
1086	READ (13) dum2d
1087	do j=jts,JMAX
1088	do i=its,IMAX
1089	soilm200_input(i,j)=dum2d(i,j)
1090	enddo
1091	enddo
1092	flag_soilm200 = 1
1093	num_sm_levels_input = num_sm_levels_input + 1
1094	sm_levels_input(num_sm_levels_input) = char2int1(var_info%name(6:8))
1095	!mp sm_inputx(:,:,num_sm_levels_input + 1) = soilm200_input
1096	do j=jts,JMAX
1097	do i=its,IMAX
1098	sm_inputx(I,J,num_sm_levels_input + 1) = soilm200_input(I,J)
1099	enddo
1100	enddo
1101
1102	ELSE IF ( var_info%name(1:8) .EQ. 'SEAICE ' ) THEN
1103	READ (13) dum2d
1104	do j=jts,JMAX
1105	do i=its,IMAX
1106	grid%xice(i,j)=dum2d(i,j)
1107	enddo
1108	enddo
1109	ELSE IF ( var_info%name(1:8) .EQ. 'WEASD ' ) THEN
1110	READ (13) dum2d
1111	do j=jts,JMAX
1112	do i=its,IMAX
1113	grid%weasd(i,j)=dum2d(i,j)
1114	enddo
1115	enddo
1116	ELSE IF ( var_info%name(1:8) .EQ. 'CANWAT ' ) THEN
1117	READ (13) dum2d
1118	do j=jts,JMAX
1119	do i=its,IMAX
1120	grid%canwat(i,j)=dum2d(i,j)
1121	enddo
1122	enddo
1123	ELSE IF ( var_info%name(1:8) .EQ. 'LANDMASK' ) THEN
1124	READ (13) dum2d
1125	do j=jts,JMAX
1126	do i=its,IMAX
1127	grid%landmask(i,j)=dum2d(i,j)
1128	enddo
1129	enddo
1130	ELSE IF ( var_info%name(1:8) .EQ. 'SKINTEMP' ) THEN
1131	READ (13) dum2d
1132	do j=jts,JMAX
1133	do i=its,IMAX
1134	grid%nmm_tsk(i,j)=dum2d(i,j)
1135	enddo
1136	enddo
1137	ELSE IF ( var_info%name(1:8) .EQ. 'TGROUND ' ) THEN
1138	READ (13) dum2d
1139	do j=jts,JMAX
1140	do i=its,IMAX
1141	grid%tg(i,j)=dum2d(i,j)
1142	enddo
1143	enddo
1144	ELSE IF ( var_info%name(1:8) .EQ. 'SOILTB ' ) THEN
1145	READ (13) dum2d
1146	do j=jts,JMAX
1147	do i=its,IMAX
1148	grid%soiltb(i,j)=dum2d(i,j)
1149	enddo
1150	enddo
1151	ELSE IF ( var_info%name(1:8) .EQ. 'SST ' ) THEN
1152	READ (13) dum2d
1153	do j=jts,JMAX
1154	do i=its,IMAX
1155	grid%sst(i,j)=dum2d(i,j)
1156	enddo
1157	enddo
1158	flag_sst = 1
1159	ELSE IF ( var_info%name(1:8) .EQ. 'GREENFRC' ) THEN
1160	READ (13) dum2d
1161	do j=jts,JMAX
1162	do i=its,IMAX
1163	grid%vegfrc(i,j)=dum2d(i,j)
1164	enddo
1165	enddo
1166	ELSE IF ( var_info%name(1:8) .EQ. 'ISLOPE ' ) THEN
1167	READ (13) dum2d
1168	do j=jts,JMAX
1169	do i=its,IMAX
1170	grid%islope(i,j)=nint(dum2d(i,j))
1171	enddo
1172	enddo
1173	ELSE IF ( var_info%name(1:8) .EQ. 'GREENMAX' ) THEN
1174	READ (13) dum2d
1175	do j=jts,JMAX
1176	do i=its,IMAX
1177	grid%greenmax(i,j)=dum2d(i,j)
1178	enddo
1179	enddo
1180	ELSE IF ( var_info%name(1:8) .EQ. 'GREENMIN' ) THEN
1181	READ (13) dum2d
1182	do j=jts,JMAX
1183	do i=its,IMAX
1184	grid%greenmin(i,j)=dum2d(i,j)
1185	enddo
1186	enddo
1187	ELSE IF ( var_info%name(1:8) .EQ. 'FIS ' ) THEN
1188	READ (13) dum2d
1189	do j=jts,JMAX
1190	do i=its,IMAX
1191	grid%fis(i,j)=dum2d(i,j)
1192	enddo
1193	enddo
1194	ELSE IF ( var_info%name(1:8) .EQ. 'Z0 ' ) THEN
1195	! ELSE IF ( var_info%name(1:8) .EQ. 'STDEV ' ) THEN
1196	READ (13) dum2d
1197	do j=jts,JMAX
1198	do i=its,IMAX
1199	grid%z0(i,j)=dum2d(i,j)
1200	enddo
1201	enddo
1202	ELSE IF ( var_info%name(1:8) .EQ. 'CMC ' ) THEN
1203	READ (13) dum2d
1204	do j=jts,JMAX
1205	do i=its,IMAX
1206	grid%cmc(i,j)=dum2d(i,j)
1207	enddo
1208	enddo
1209	ELSE IF ( var_info%name(1:8) .EQ. 'HTM ' ) THEN
1210	READ (13) dum3d
1211	do k=kts,kte-1
1212	do j=jts,JMAX
1213	do i=its,IMAX
1214	htm_in(i,j,k)=dum3d(i,j,k)
1215	enddo
1216	enddo
1217	enddo
1218	ELSE IF ( var_info%name(1:8) .EQ. 'VTM ' ) THEN
1219	READ (13) dum3d
1220	do k=kts,kte-1
1221	do j=jts,JMAX
1222	do i=its,IMAX
1223	vtm_in(i,j,k)=dum3d(i,j,k)
1224	enddo
1225	enddo
1226	enddo
1227	ELSE IF ( var_info%name(1:8) .EQ. 'SM ' ) THEN
1228	READ (13) dum2d
1229	do j=jts,JMAX
1230	do i=its,IMAX
1231	grid%sm(i,j)=dum2d(i,j)
1232	enddo
1233	enddo
1234	ELSE IF ( var_info%name(1:8) .EQ. 'ALBASE ' ) THEN
1235	READ (13) dum2d
1236	do j=jts,JMAX
1237	do i=its,IMAX
1238	grid%albase(i,j)=dum2d(i,j)
1239	enddo
1240	enddo
1241	ELSE IF ( var_info%name(1:8) .EQ. 'MXSNAL ' ) THEN
1242	READ (13) dum2d
1243	do j=jts,JMAX
1244	do i=its,IMAX
1245	grid%mxsnal(i,j)=dum2d(i,j)
1246	enddo
1247	enddo
1248
1249	! 1D vertical coordinate.
1250
1251	ELSE IF ( var_info%name(1:8) .EQ. 'DETA ' ) THEN
1252	READ(13) DETA_in
1253	ELSE IF ( var_info%name(1:8) .EQ. 'DETA1 ' ) THEN
1254	READ(13) DETA1_in
1255	ELSE IF ( var_info%name(1:8) .EQ. 'DETA2 ' ) THEN
1256	READ(13) DETA2_in
1257	ELSE IF ( var_info%name(1:8) .EQ. 'ETAX ' ) THEN
1258	READ(13) ETAX_in
1259	ELSE IF ( var_info%name(1:8) .EQ. 'ETA1 ' ) THEN
1260	READ(13) ETA1_in
1261	ELSE IF ( var_info%name(1:8) .EQ. 'ETA2 ' ) THEN
1262	READ(13) ETA2_in
1263	ELSE IF ( var_info%name(1:8) .EQ. 'AETA ' ) THEN
1264	READ(13) AETA_in
1265	ELSE IF ( var_info%name(1:8) .EQ. 'AETA1 ' ) THEN
1266	READ(13) AETA1_in
1267	ELSE IF ( var_info%name(1:8) .EQ. 'AETA2 ' ) THEN
1268	READ(13) AETA2_in
1269	ELSE IF ( var_info%name(1:8) .EQ. 'DFL ' ) THEN
1270	READ(13) DFL_in
1271
1272	! ELSE IF ( var_info%name(1:8) .EQ. 'ETAPHALF' ) THEN
1273	! READ (13) etahalf
1274	! ELSE IF ( var_info%name(1:8) .EQ. 'ETAPFULL' ) THEN
1275	! READ (13) etafull
1276
1277	! wrong input data.
1278
1279	ELSE IF ( var_info%name(1:8) .EQ. 'ZETAFULL' ) THEN
1280	PRINT '(A)','Oops, you put in the height data.'
1281	STOP 'this_is_mass_not_height'
1282
1283
1284	! Stuff that we do not want or need is just skipped over.
1285
1286	ELSE
1287	print *,'------------------> skipping ', var_info%name(1:8)
1288	READ (13) dummy
1289	END IF
1290
1291	END DO read_all_the_data
1292
1293	CLOSE (13)
1294
1295	first_time_in = .FALSE.
1296
1297	!new
1298	sw_inputx=0.
1299	!new
1300
1301	do k=kts,kte-1
1302	do j=jts,JMAX
1303	do i=its,IMAX
1304	grid%U(I,J,K)=U_input(I,J,K)
1305	grid%V(I,J,K)=V_input(I,J,K)
1306	grid%T(I,J,K)=T_input(I,J,K)
1307	grid%Q(I,J,K)=Q_input(I,J,K)
1308	enddo
1309	enddo
1310	enddo
1311
1312	! write(0,*) 'size sw_input: ', size(sw_input,dim=1),size(sw_input,dim=2),size(sw_input,dim=3)
1313	! write(0,*) 'size sw_inputx: ', size(sw_inputx,dim=1),size(sw_inputx,dim=2),size(sw_inputx,dim=3)
1314	sw_input=0.
1315
1316	! write(0,*) 'maxval st_inputx(1): ', maxval(st_input(:,:,1))
1317	! write(0,*) 'maxval st_inputx(2): ', maxval(st_input(:,:,2))
1318	! write(0,*) 'maxval st_inputx(3): ', maxval(st_input(:,:,3))
1319	! write(0,*) 'maxval st_inputx(4): ', maxval(st_input(:,:,4))
1320
1321
1322	do J=JTS,min(JDE-1,JTE)
1323	do K=1,num_st_levels_alloc
1324	do I=ITS,min(IDE-1,ITE)
1325	st_input(I,K,J)=st_inputx(I,J,K)
1326	sm_input(I,K,J)=sm_inputx(I,J,K)
1327	sw_input(I,K,J)=sw_inputx(I,J,K)
1328	enddo
1329	enddo
1330	enddo
1331
1332	! write(0,*) 'maxval st_input(1): ', maxval(st_input(:,1,:))
1333	! write(0,*) 'maxval st_input(2): ', maxval(st_input(:,2,:))
1334	! write(0,*) 'maxval st_input(3): ', maxval(st_input(:,3,:))
1335	! write(0,*) 'maxval st_input(4): ', maxval(st_input(:,4,:))
1336
1337
1338	num_veg_cat = SIZE ( grid%landusef , DIM=2 )
1339	num_soil_top_cat = SIZE ( grid%soilctop , DIM=2 )
1340	num_soil_bot_cat = SIZE ( grid%soilcbot , DIM=2 )
1341
1342	do J=JTS,min(JDE-1,JTE)
1343	do K=1,num_soil_top_cat
1344	do I=ITS,min(IDE-1,ITE)
1345	grid%SOILCTOP(I,K,J)=soil_top_cat_input(I,J,K)
1346	grid%SOILCTOP_gc(I,J,K)=soil_top_cat_input(I,J,K)
1347	enddo
1348	enddo
1349	enddo
1350
1351	do J=JTS,min(JDE-1,JTE)
1352	do K=1,num_soil_bot_cat
1353	do I=ITS,min(IDE-1,ITE)
1354	grid%SOILCBOT(I,K,J)=soil_bot_cat_input(I,J,K)
1355	grid%SOILCBOT_gc(I,J,K)=soil_bot_cat_input(I,J,K)
1356	enddo
1357	enddo
1358	enddo
1359
1360	do J=JTS,min(JDE-1,JTE)
1361	do K=1,num_veg_cat
1362	do I=ITS,min(IDE-1,ITE)
1363	grid%LANDUSEF(I,K,J)=landuse_frac_input(I,J,K)
1364	grid%LANDUSEF_gc(I,J,K)=landuse_frac_input(I,J,K)
1365	enddo
1366	enddo
1367	enddo
1368
1369
1370	do K=KDS,KDE
1371	grid%ETAX(K)=ETAX_in(KDE-K+1)
1372	grid%ETA1(K)=ETA1_in(KDE-K+1)
1373	grid%ETA2(K)=ETA2_in(KDE-K+1)
1374	grid%DFL(K)=DFL_in(KDE-K+1)
1375	enddo
1376
1377	do K=KDS,KDE-1
1378	grid%DETA(K)=DETA_in(KDE-K)
1379	grid%DETA1(K)=DETA1_in(KDE-K)
1380	grid%DETA2(K)=DETA2_in(KDE-K)
1381	grid%AETA(K)=AETA_in(KDE-K)
1382	grid%AETA1(K)=AETA1_in(KDE-K)
1383	grid%AETA2(K)=AETA2_in(KDE-K)
1384	enddo
1385
1386	END SUBROUTINE read_si
1387
1388	! ------------------------------------------------------------
1389	! ------------------------------------------------------------
1390	! ------------------------------------------------------------
1391	! ------------------------------------------------------------
1392	! ------------------------------------------------------------
1393	! ------------------------------------------------------------
1394
1395	SUBROUTINE read_wps ( grid, filename, file_date_string, num_metgrid_levels )
1396
1397	USE module_soil_pre
1398	USE module_domain
1399
1400	IMPLICIT NONE
1401
1402	TYPE(domain) , INTENT(INOUT) :: grid
1403	CHARACTER (LEN=19) , INTENT(IN) :: file_date_string
1404	CHARACTER (LEN=19) :: VarName
1405	CHARACTER (LEN=150) :: chartemp
1406	CHARACTER (LEN=256) :: mminlu_loc, message
1407	CHARACTER (*) , INTENT(IN) :: filename
1408
1409	INTEGER :: ids,ide,jds,jde,kds,kde &
1410	,ims,ime,jms,jme,kms,kme &
1411	,its,ite,jts,jte,kts,kte
1412
1413	INTEGER :: i , j , k , loop, IMAX, JMAX
1414	INTEGER :: DATAHANDLE, num_metgrid_levels
1415	INTEGER :: Sysdepinfo, Status, N
1416	INTEGER :: istatus,ioutcount,iret,index,ierr
1417
1418	integer :: nrecs,iunit, L,hor_size
1419
1420	!!
1421	character*132, allocatable,save :: datestr_all(:)
1422	character*132, allocatable,save :: varname_all(:)
1423	integer, allocatable,save :: domainend_all(:,:)
1424	integer, allocatable,save :: start_block(:)
1425	integer, allocatable,save :: end_block(:)
1426	integer, allocatable,save :: start_byte(:)
1427	integer, allocatable,save :: end_byte(:)
1428	integer(kind=i_llong), allocatable,save :: file_offset(:)
1429	!!
1430
1431	REAL :: dummy,tmp,garb
1432	REAL, ALLOCATABLE,save :: dumdata(:,:,:)
1433
1434	CHARACTER (LEN= 8) :: dummy_char
1435
1436	INTEGER :: ok , map_proj , ok_open, igarb
1437	REAL :: pt
1438	INTEGER :: num_veg_cat , num_soil_top_cat , num_soil_bot_cat
1439
1440	#if defined(DM_PARALLEL) && !defined(STUBMPI)
1441	INCLUDE "mpif.h"
1442
1443	SELECT CASE ( model_data_order )
1444	CASE ( DATA_ORDER_ZXY )
1445	kds = grid%sd31 ; kde = grid%ed31 ;
1446	ids = grid%sd32 ; ide = grid%ed32 ;
1447	jds = grid%sd33 ; jde = grid%ed33 ;
1448
1449	kms = grid%sm31 ; kme = grid%em31 ;
1450	ims = grid%sm32 ; ime = grid%em32 ;
1451	jms = grid%sm33 ; jme = grid%em33 ;
1452
1453	kts = grid%sp31 ; kte = grid%ep31 ; ! tile is entire patch
1454	its = grid%sp32 ; ite = grid%ep32 ; ! tile is entire patch
1455	jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
1456
1457	CASE ( DATA_ORDER_XYZ )
1458	ids = grid%sd31 ; ide = grid%ed31 ;
1459	jds = grid%sd32 ; jde = grid%ed32 ;
1460	kds = grid%sd33 ; kde = grid%ed33 ;
1461
1462	ims = grid%sm31 ; ime = grid%em31 ;
1463	jms = grid%sm32 ; jme = grid%em32 ;
1464	kms = grid%sm33 ; kme = grid%em33 ;
1465
1466	its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
1467	jts = grid%sp32 ; jte = grid%ep32 ; ! tile is entire patch
1468	kts = grid%sp33 ; kte = grid%ep33 ; ! tile is entire patch
1469
1470	CASE ( DATA_ORDER_XZY )
1471	ids = grid%sd31 ; ide = grid%ed31 ;
1472	kds = grid%sd32 ; kde = grid%ed32 ;
1473	jds = grid%sd33 ; jde = grid%ed33 ;
1474
1475	ims = grid%sm31 ; ime = grid%em31 ;
1476	kms = grid%sm32 ; kme = grid%em32 ;
1477	jms = grid%sm33 ; jme = grid%em33 ;
1478
1479	its = grid%sp31 ; ite = grid%ep31 ; ! tile is entire patch
1480	kts = grid%sp32 ; kte = grid%ep32 ; ! tile is entire patch
1481	jts = grid%sp33 ; jte = grid%ep33 ; ! tile is entire patch
1482
1483	END SELECT
1484
1485	! Initialize what soil temperature and moisture is available.
1486
1487
1488	flag_st000010 = 0
1489	flag_st010040 = 0
1490	flag_st040100 = 0
1491	flag_st100200 = 0
1492	flag_sm000010 = 0
1493	flag_sm010040 = 0
1494	flag_sm040100 = 0
1495	flag_sm100200 = 0
1496	flag_st010200 = 0
1497	flag_sm010200 = 0
1498
1499	flag_soilt010 = 0
1500	flag_soilt040 = 0
1501	flag_soilt100 = 0
1502	flag_soilt200 = 0
1503	flag_soilm010 = 0
1504	flag_soilm040 = 0
1505	flag_soilm100 = 0
1506	flag_soilm200 = 0
1507
1508	flag_sst = 0
1509	flag_toposoil = 0
1510
1511	! How many soil levels have we found? Well, right now, none.
1512
1513	num_st_levels_input = 0
1514	num_sm_levels_input = 0
1515	st_levels_input = -1
1516	sm_levels_input = -1
1517
1518	mminlu_loc=''
1519	mminlu_loc(1:4)='USGS'
1520
1521	CALL nl_set_mminlu ( grid%id, mminlu_loc)
1522	CALL nl_set_iswater (grid%id, 16 )
1523	CALL nl_set_isice (grid%id, 24 )
1524	CALL nl_set_islake (grid%id, -1)
1525
1526
1527	! Get the space for the data if this is the first time here.
1528
1529	IF (.NOT. ALLOCATED (pmsl) ) ALLOCATE ( pmsl(its:ite,jts:jte) )
1530	IF (.NOT. ALLOCATED (psfc_in) ) ALLOCATE ( psfc_in(its:ite,jts:jte) )
1531	IF (.NOT. ALLOCATED (st_inputx)) ALLOCATE (st_inputx(its:ite,jts:jte,num_st_levels_alloc))
1532	IF (.NOT. ALLOCATED (sm_inputx)) ALLOCATE (sm_inputx(its:ite,jts:jte,num_st_levels_alloc))
1533	IF (.NOT. ALLOCATED (sw_inputx)) ALLOCATE (sw_inputx(its:ite,jts:jte,num_st_levels_alloc))
1534	IF (.NOT. ALLOCATED (soilt010_input) ) ALLOCATE ( soilt010_input(its:ite,jts:jte) )
1535	IF (.NOT. ALLOCATED (soilt040_input) ) ALLOCATE ( soilt040_input(its:ite,jts:jte) )
1536	IF (.NOT. ALLOCATED (soilt100_input) ) ALLOCATE ( soilt100_input(its:ite,jts:jte) )
1537	IF (.NOT. ALLOCATED (soilt200_input) ) ALLOCATE ( soilt200_input(its:ite,jts:jte) )
1538	IF (.NOT. ALLOCATED (soilm010_input) ) ALLOCATE ( soilm010_input(its:ite,jts:jte) )
1539	IF (.NOT. ALLOCATED (soilm040_input) ) ALLOCATE ( soilm040_input(its:ite,jts:jte) )
1540	IF (.NOT. ALLOCATED (soilm100_input) ) ALLOCATE ( soilm100_input(its:ite,jts:jte) )
1541	IF (.NOT. ALLOCATED (soilm200_input) ) ALLOCATE ( soilm200_input(its:ite,jts:jte) )
1542
1543	! Local arrays
1544
1545	!!! MPI IO
1546
1547	iunit=33
1548	call count_recs_wrf_binary_file(iunit, trim(fileName), nrecs)
1549
1550	if (.NOT. ALLOCATED (datestr_all)) allocate (datestr_all(nrecs))
1551	if (.NOT. ALLOCATED (varname_all)) allocate (varname_all(nrecs))
1552	if (.NOT. ALLOCATED (domainend_all)) allocate (domainend_all(3,nrecs))
1553	if (.NOT. ALLOCATED (start_block)) allocate (start_block(nrecs))
1554	if (.NOT. ALLOCATED (end_block)) allocate (end_block(nrecs))
1555	if (.NOT. ALLOCATED (start_byte)) allocate (start_byte(nrecs))
1556	if (.NOT. ALLOCATED (end_byte)) allocate (end_byte(nrecs))
1557	if (.NOT. ALLOCATED (file_offset)) allocate (file_offset(nrecs))
1558
1559	call inventory_wrf_binary_file(iunit, trim(filename), nrecs, &
1560	datestr_all,varname_all,domainend_all, &
1561	start_block,end_block,start_byte,end_byte,file_offset)
1562
1563	do N=1,NRECS
1564	write(message,*) 'N,varname_all(N): ',N, varname_all(N)
1565	call wrf_debug(100,message)
1566	enddo
1567
1568	call mpi_file_open(mpi_comm_world, trim(filename), &
1569	mpi_mode_rdonly,mpi_info_null, iunit, ierr)
1570	if (ierr /= 0) then
1571	call wrf_error_fatal("Error opening file with mpi io")
1572	end if
1573
1574	VarName='CEN_LAT'
1575	call retrieve_index(index,VarName,varname_all,nrecs,iret)
1576	if (iret /= 0) then
1577	write(message,*) VarName, ' not found in file'
1578	call wrf_message(message)
1579	else
1580
1581	call mpi_file_read_at(iunit,file_offset(index)+5*4, &
1582	garb,1,mpi_real4, &
1583	mpi_status_ignore, ierr)
1584
1585	if (ierr /= 0) then
1586	write(message,*) 'Error reading ', VarName,' using MPIIO'
1587	call wrf_message(message)
1588	else
1589	write(message,*) VarName, ' from MPIIO READ= ', garb
1590	call wrf_message(message)
1591	CALL nl_set_cen_lat ( grid%id , garb )
1592	end if
1593	end if
1594
1595	VarName='CEN_LON'
1596	call retrieve_index(index,VarName,varname_all,nrecs,iret)
1597	call mpi_file_read_at(iunit,file_offset(index)+5*4, &
1598	garb,1,mpi_real4, &
1599	mpi_status_ignore, ierr)
1600	CALL nl_set_cen_lon ( grid%id , garb )
1601	CALL nl_set_stand_lon ( grid%id , garb )
1602
1603	VarName='TRUELAT1'
1604	call retrieve_index(index,VarName,varname_all,nrecs,iret)
1605	call mpi_file_read_at(iunit,file_offset(index)+5*4, &
1606	garb,1,mpi_real4, &
1607	mpi_status_ignore, ierr)
1608	CALL nl_set_truelat1 ( grid%id , garb )
1609
1610	VarName='TRUELAT2'
1611	call retrieve_index(index,VarName,varname_all,nrecs,iret)
1612	call mpi_file_read_at(iunit,file_offset(index)+5*4, &
1613	garb,1,mpi_real4, &
1614	mpi_status_ignore, ierr)
1615	CALL nl_set_truelat2 ( grid%id , garb )
1616
1617	VarName='MAP_PROJ'
1618	call retrieve_index(index,VarName,varname_all,nrecs,iret)
1619	call mpi_file_read_at(iunit,file_offset(index)+5*4, &
1620	igarb,1,mpi_integer4, &
1621	mpi_status_ignore, ierr)
1622
1623	CALL nl_set_map_proj( grid%id, igarb)
1624
1625	hor_size=(IDE-IDS)*(JDE-JDS)
1626
1627	varName='PRES'
1628	if (.not. allocated(dumdata)) then
1629	allocate(dumdata(IDS:IDE-1,JDS:JDE-1,num_metgrid_levels))
1630	endif
1631
1632	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1633	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1634	dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1635	mpi_status_ignore, ierr)
1636
1637	DO K=1,num_metgrid_levels
1638	DO J=JTS,min(JTE,JDE-1)
1639	DO I=ITS,min(ITE,IDE-1)
1640	grid%p_gc(I,J,K)=dumdata(I,J,K)
1641	ENDDO
1642	ENDDO
1643	ENDDO
1644
1645	varName='GHT'
1646
1647	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1648	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1649	dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1650	mpi_status_ignore, ierr)
1651
1652
1653	DO K=1,num_metgrid_levels
1654	DO J=JTS,min(JTE,JDE-1)
1655	DO I=ITS,min(ITE,IDE-1)
1656	grid%ght_gc(I,J,K)=dumdata(I,J,K)
1657	ENDDO
1658	ENDDO
1659	ENDDO
1660
1661
1662	varName='VEGCAT'
1663
1664	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1665	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1666	dumdata,hor_size,mpi_real4, &
1667	mpi_status_ignore, ierr)
1668
1669	DO J=JTS,min(JTE,JDE-1)
1670	DO I=ITS,min(ITE,IDE-1)
1671	grid%vegcat(I,J)=dumdata(I,J,1)
1672	ENDDO
1673	ENDDO
1674
1675	varName='SOIL_CAT'
1676	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1677	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1678	dumdata,hor_size,mpi_real4, &
1679	mpi_status_ignore, ierr)
1680	DO J=JTS,min(JTE,JDE-1)
1681	DO I=ITS,min(ITE,IDE-1)
1682	grid%input_soil_cat(I,J)=dumdata(I,J,1)
1683	ENDDO
1684	ENDDO
1685
1686	varName='CANWAT'
1687	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1688	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1689	dumdata,hor_size,mpi_real4, &
1690	mpi_status_ignore, ierr)
1691
1692	DO J=JTS,min(JTE,JDE-1)
1693	DO I=ITS,min(ITE,IDE-1)
1694	grid%canwat(I,J)=dumdata(I,J,1)
1695	ENDDO
1696	ENDDO
1697
1698	varName='SNOW'
1699	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1700	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1701	dumdata,hor_size,mpi_real4, &
1702	mpi_status_ignore, ierr)
1703	DO J=JTS,min(JTE,JDE-1)
1704	DO I=ITS,min(ITE,IDE-1)
1705	grid%snow(I,J)=dumdata(I,J,1)
1706	ENDDO
1707	ENDDO
1708
1709	varName='SKINTEMP'
1710	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1711	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1712	dumdata,hor_size,mpi_real4, &
1713	mpi_status_ignore, ierr)
1714	DO J=JTS,min(JTE,JDE-1)
1715	DO I=ITS,min(ITE,IDE-1)
1716	grid%tsk_gc(I,J)=dumdata(I,J,1)
1717	ENDDO
1718	ENDDO
1719
1720	varName='SOILHGT'
1721	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1722	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1723	dumdata,hor_size,mpi_real4, &
1724	mpi_status_ignore, ierr)
1725	DO J=JTS,min(JTE,JDE-1)
1726	DO I=ITS,min(ITE,IDE-1)
1727	grid%toposoil(I,J)=dumdata(I,J,1)
1728	ENDDO
1729	ENDDO
1730
1731	varName='SEAICE'
1732	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1733	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1734	dumdata,hor_size,mpi_real4, &
1735	mpi_status_ignore, ierr)
1736	DO J=JTS,min(JTE,JDE-1)
1737	DO I=ITS,min(ITE,IDE-1)
1738	grid%xice_gc(I,J)=dumdata(I,J,1)
1739	ENDDO
1740	ENDDO
1741
1742	varName='ST100200'
1743	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1744	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1745	dumdata,hor_size,mpi_real4, &
1746	mpi_status_ignore, ierr)
1747	DO J=JTS,min(JTE,JDE-1)
1748	DO I=ITS,min(ITE,IDE-1)
1749	grid%st100200(I,J)=dumdata(I,J,1)
1750	ENDDO
1751	ENDDO
1752
1753	varName='ST040100'
1754	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1755	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1756	dumdata,hor_size,mpi_real4, &
1757	mpi_status_ignore, ierr)
1758	DO J=JTS,min(JTE,JDE-1)
1759	DO I=ITS,min(ITE,IDE-1)
1760	grid%st040100(I,J)=dumdata(I,J,1)
1761	ENDDO
1762	ENDDO
1763
1764	varName='ST010040'
1765	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1766	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1767	dumdata,hor_size,mpi_real4, &
1768	mpi_status_ignore, ierr)
1769
1770	DO J=JTS,min(JTE,JDE-1)
1771	DO I=ITS,min(ITE,IDE-1)
1772	grid%st010040(I,J)=dumdata(I,J,1)
1773	ENDDO
1774	ENDDO
1775
1776	varName='ST000010'
1777	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1778	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1779	dumdata,hor_size,mpi_real4, &
1780	mpi_status_ignore, ierr)
1781
1782	DO J=JTS,min(JTE,JDE-1)
1783	DO I=ITS,min(ITE,IDE-1)
1784	grid%st000010(I,J)=dumdata(I,J,1)
1785	ENDDO
1786	ENDDO
1787
1788	varName='SM100200'
1789	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1790	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1791	dumdata,hor_size,mpi_real4, &
1792	mpi_status_ignore, ierr)
1793
1794	DO J=JTS,min(JTE,JDE-1)
1795	DO I=ITS,min(ITE,IDE-1)
1796	grid%sm100200(I,J)=dumdata(I,J,1)
1797	ENDDO
1798	ENDDO
1799
1800	varName='SM040100'
1801	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1802	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1803	dumdata,hor_size,mpi_real4, &
1804	mpi_status_ignore, ierr)
1805
1806	DO J=JTS,min(JTE,JDE-1)
1807	DO I=ITS,min(ITE,IDE-1)
1808	grid%sm040100(I,J)=dumdata(I,J,1)
1809	ENDDO
1810	ENDDO
1811
1812	varName='SM010040'
1813	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1814	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1815	dumdata,hor_size,mpi_real4, &
1816	mpi_status_ignore, ierr)
1817
1818	DO J=JTS,min(JTE,JDE-1)
1819	DO I=ITS,min(ITE,IDE-1)
1820	grid%sm010040(I,J)=dumdata(I,J,1)
1821	ENDDO
1822	ENDDO
1823
1824	varName='SM000010'
1825	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1826	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1827	dumdata,hor_size,mpi_real4, &
1828	mpi_status_ignore, ierr)
1829
1830	DO J=JTS,min(JTE,JDE-1)
1831	DO I=ITS,min(ITE,IDE-1)
1832	grid%sm000010(I,J)=dumdata(I,J,1)
1833	ENDDO
1834	ENDDO
1835
1836	varName='PSFC'
1837	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1838	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1839	dumdata,hor_size,mpi_real4, &
1840	mpi_status_ignore, ierr)
1841
1842	DO J=JTS,min(JTE,JDE-1)
1843	DO I=ITS,min(ITE,IDE-1)
1844	grid%psfc(I,J)=dumdata(I,J,1)
1845	ENDDO
1846	ENDDO
1847
1848	varName='RH'
1849	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1850	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1851	dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1852	mpi_status_ignore, ierr)
1853
1854	DO K=1,num_metgrid_levels
1855	DO J=JTS,min(JTE,JDE-1)
1856	DO I=ITS,min(ITE,IDE-1)
1857	grid%rh_gc(I,J,K)=dumdata(I,J,K)
1858	ENDDO
1859	ENDDO
1860	ENDDO
1861
1862	varName='VV'
1863	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1864	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1865	dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1866	mpi_status_ignore, ierr)
1867	DO K=1,num_metgrid_levels
1868	DO J=JTS,min(JTE,JDE-1)
1869	DO I=ITS,min(ITE,IDE-1)
1870	grid%v_gc(I,J,K)=dumdata(I,J,K)
1871	ENDDO
1872	ENDDO
1873	ENDDO
1874
1875	varName='UU'
1876	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1877	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1878	dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1879	mpi_status_ignore, ierr)
1880	DO K=1,num_metgrid_levels
1881	DO J=JTS,min(JTE,JDE-1)
1882	DO I=ITS,min(ITE,IDE-1)
1883	grid%u_gc(I,J,K)=dumdata(I,J,K)
1884	ENDDO
1885	ENDDO
1886	ENDDO
1887
1888	varName='TT'
1889	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1890
1891	if (IRET .eq. 0) then
1892	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1893	dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1894	mpi_status_ignore, ierr)
1895	DO K=1,num_metgrid_levels
1896	DO J=JTS,min(JTE,JDE-1)
1897	DO I=ITS,min(ITE,IDE-1)
1898	grid%t_gc(I,J,K)=dumdata(I,J,K)
1899	ENDDO
1900	ENDDO
1901	ENDDO
1902	endif
1903
1904	! varName='RWMR'
1905	! CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1906	! if (iret .ne. 0) then
1907	! grid%rwmr_gc=0.
1908	! grid%snmr_gc=0.
1909	! grid%clwmr_gc=0.
1910	! grid%cice_gc=0.
1911	! grid%rimef_gc=0.
1912	! write(0,*) 'skipping cloud reads'
1913	! goto 979
1914	! endif
1915
1916	! CALL mpi_file_read_at(iunit,file_offset(index+1), &
1917	! dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1918	! mpi_status_ignore, ierr)
1919	! DO K=1,num_metgrid_levels
1920	! DO J=JTS,min(JTE,JDE-1)
1921	! DO I=ITS,min(ITE,IDE-1)
1922	! grid%rwmr_gc(I,J,K)=dumdata(I,J,K)
1923	! ENDDO
1924	! ENDDO
1925	! ENDDO
1926
1927	! varName='SNMR'
1928	! CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1929	! CALL mpi_file_read_at(iunit,file_offset(index+1), &
1930	! dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1931	! mpi_status_ignore, ierr)
1932	! DO K=1,num_metgrid_levels
1933	! DO J=JTS,min(JTE,JDE-1)
1934	! DO I=ITS,min(ITE,IDE-1)
1935	! grid%snmr_gc(I,J,K)=dumdata(I,J,K)
1936	! ENDDO
1937	! ENDDO
1938	! ENDDO
1939
1940	! varName='CLWMR'
1941	! CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1942	! CALL mpi_file_read_at(iunit,file_offset(index+1), &
1943	! dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1944	! mpi_status_ignore, ierr)
1945	! DO K=1,num_metgrid_levels
1946	! DO J=JTS,min(JTE,JDE-1)
1947	! DO I=ITS,min(ITE,IDE-1)
1948	! grid%clwmr_gc(I,J,K)=dumdata(I,J,K)
1949	! ENDDO
1950	! ENDDO
1951	! ENDDO
1952
1953	! varName='CICE'
1954	! CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1955	! CALL mpi_file_read_at(iunit,file_offset(index+1), &
1956	! dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1957	! mpi_status_ignore, ierr)
1958	! DO K=1,num_metgrid_levels
1959	! DO J=JTS,min(JTE,JDE-1)
1960	! DO I=ITS,min(ITE,IDE-1)
1961	! grid%cice_gc(I,J,K)=dumdata(I,J,K)
1962	! ENDDO
1963	! ENDDO
1964	! ENDDO
1965
1966	! varName='FRIMEF'
1967	! CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1968	! CALL mpi_file_read_at(iunit,file_offset(index+1), &
1969	! dumdata,hor_size*num_metgrid_levels,mpi_real4, &
1970	! mpi_status_ignore, ierr)
1971	! DO K=1,num_metgrid_levels
1972	! DO J=JTS,min(JTE,JDE-1)
1973	! DO I=ITS,min(ITE,IDE-1)
1974	! grid%rimef_gc(I,J,K)=dumdata(I,J,K)
1975	! ENDDO
1976	! ENDDO
1977	! ENDDO
1978
1979	979 continue
1980
1981	! varName='PMSL'
1982
1983	varName='SLOPECAT'
1984	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1985	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1986	dumdata,hor_size,mpi_real4, &
1987	mpi_status_ignore, ierr)
1988	DO J=JTS,min(JTE,JDE-1)
1989	DO I=ITS,min(ITE,IDE-1)
1990	grid%slopecat(I,J)=dumdata(I,J,1)
1991	ENDDO
1992	ENDDO
1993
1994	varName='SNOALB'
1995	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
1996	CALL mpi_file_read_at(iunit,file_offset(index+1), &
1997	dumdata,hor_size,mpi_real4, &
1998	mpi_status_ignore, ierr)
1999
2000	DO J=JTS,min(JTE,JDE-1)
2001	DO I=ITS,min(ITE,IDE-1)
2002	grid%snoalb(I,J)=dumdata(I,J,1)
2003	ENDDO
2004	ENDDO
2005
2006	num_veg_cat = SIZE ( grid%landusef_gc , DIM=3 )
2007	num_soil_top_cat = SIZE ( grid%soilctop_gc , DIM=3 )
2008	num_soil_bot_cat = SIZE ( grid%soilcbot_gc , DIM=3 )
2009
2010	varName='GREENFRAC'
2011	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2012	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2013	dumdata,hor_size*12,mpi_real4, &
2014	mpi_status_ignore, ierr)
2015
2016	DO K=1,12
2017	DO J=JTS,min(JTE,JDE-1)
2018	DO I=ITS,min(ITE,IDE-1)
2019	grid%greenfrac_gc(I,J,K)=dumdata(I,J,K)
2020	ENDDO
2021	ENDDO
2022	ENDDO
2023
2024	varName='ALBEDO12M'
2025	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2026	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2027	dumdata,hor_size*12,mpi_real4, &
2028	mpi_status_ignore, ierr)
2029
2030	DO K=1,12
2031	DO J=JTS,min(JTE,JDE-1)
2032	DO I=ITS,min(ITE,IDE-1)
2033	grid%albedo12m_gc(I,J,K)=dumdata(I,J,K)
2034	ENDDO
2035	ENDDO
2036	ENDDO
2037
2038	varName='SOILCBOT'
2039	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2040	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2041	dumdata,hor_size*num_soil_bot_cat,mpi_real4, &
2042	mpi_status_ignore, ierr)
2043
2044	DO K=1,num_soil_bot_cat
2045	DO J=JTS,min(JTE,JDE-1)
2046	DO I=ITS,min(ITE,IDE-1)
2047	grid%soilcbot_gc(I,J,K)=dumdata(I,J,K)
2048	ENDDO
2049	ENDDO
2050	ENDDO
2051
2052	varName='SOILCAT'
2053	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2054	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2055	dumdata,hor_size,mpi_real4, &
2056	mpi_status_ignore, ierr)
2057
2058	DO J=JTS,min(JTE,JDE-1)
2059	DO I=ITS,min(ITE,IDE-1)
2060	grid%soilcat(I,J)=dumdata(I,J,1)
2061	ENDDO
2062	ENDDO
2063
2064
2065	varName='SOILCTOP'
2066	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2067	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2068	dumdata,hor_size*num_soil_top_cat,mpi_real4, &
2069	mpi_status_ignore, ierr)
2070	DO K=1,num_soil_top_cat
2071	DO J=JTS,min(JTE,JDE-1)
2072	DO I=ITS,min(ITE,IDE-1)
2073	grid%soilctop_gc(I,J,K)=dumdata(I,J,K)
2074	ENDDO
2075	ENDDO
2076	ENDDO
2077
2078	varName='SOILTEMP'
2079	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2080	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2081	dumdata,hor_size,mpi_real4, &
2082	mpi_status_ignore, ierr)
2083
2084	DO J=JTS,min(JTE,JDE-1)
2085	DO I=ITS,min(ITE,IDE-1)
2086	grid%tmn_gc(I,J)=dumdata(I,J,1)
2087	ENDDO
2088	ENDDO
2089
2090	varName='HGT_V'
2091	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2092	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2093	dumdata,hor_size,mpi_real4, &
2094	mpi_status_ignore, ierr)
2095
2096	DO J=JTS,min(JTE,JDE-1)
2097	DO I=ITS,min(ITE,IDE-1)
2098	grid%htv_gc(I,J)=dumdata(I,J,1)
2099	ENDDO
2100	ENDDO
2101
2102	varName='HGT_M'
2103	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2104	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2105	dumdata,hor_size,mpi_real4, &
2106	mpi_status_ignore, ierr)
2107
2108	DO J=JTS,min(JTE,JDE-1)
2109	DO I=ITS,min(ITE,IDE-1)
2110	grid%ht_gc(I,J)=dumdata(I,J,1)
2111	ENDDO
2112	ENDDO
2113
2114	varName='LU_INDEX'
2115	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2116	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2117	dumdata,hor_size,mpi_real4, &
2118	mpi_status_ignore, ierr)
2119
2120	DO J=JTS,min(JTE,JDE-1)
2121	DO I=ITS,min(ITE,IDE-1)
2122	grid%lu_index(I,J)=dumdata(I,J,1)
2123	ENDDO
2124	ENDDO
2125
2126	varName='LANDUSEF'
2127	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2128	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2129	dumdata,hor_size*num_veg_cat,mpi_real4, &
2130	mpi_status_ignore, ierr)
2131
2132	DO K=1,num_veg_cat
2133	DO J=JTS,min(JTE,JDE-1)
2134	DO I=ITS,min(ITE,IDE-1)
2135	grid%landusef_gc(I,J,K)=dumdata(I,J,K)
2136	ENDDO
2137	ENDDO
2138	ENDDO
2139
2140	varName='LANDMASK'
2141	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2142	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2143	dumdata,hor_size,mpi_real4, &
2144	mpi_status_ignore, ierr)
2145
2146	DO J=JTS,min(JTE,JDE-1)
2147	DO I=ITS,min(ITE,IDE-1)
2148	grid%landmask(I,J)=dumdata(I,J,1)
2149	ENDDO
2150	ENDDO
2151
2152
2153	varName='XLONG_V'
2154	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2155	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2156	dumdata,hor_size,mpi_real4, &
2157	mpi_status_ignore, ierr)
2158
2159	DO J=JTS,min(JTE,JDE-1)
2160	DO I=ITS,min(ITE,IDE-1)
2161	grid%vlon_gc(I,J)=dumdata(I,J,1)
2162	ENDDO
2163	ENDDO
2164
2165	varName='XLAT_V'
2166	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2167	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2168	dumdata,hor_size,mpi_real4, &
2169	mpi_status_ignore, ierr)
2170
2171	DO J=JTS,min(JTE,JDE-1)
2172	DO I=ITS,min(ITE,IDE-1)
2173	grid%vlat_gc(I,J)=dumdata(I,J,1)
2174	ENDDO
2175	ENDDO
2176
2177	varName='XLONG_M'
2178	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2179	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2180	dumdata,hor_size,mpi_real4, &
2181	mpi_status_ignore, ierr)
2182
2183	DO J=JTS,min(JTE,JDE-1)
2184	DO I=ITS,min(ITE,IDE-1)
2185	grid%hlon_gc(I,J)=dumdata(I,J,1)
2186	ENDDO
2187	ENDDO
2188
2189	varName='XLAT_M'
2190	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2191	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2192	dumdata,hor_size,mpi_real4, &
2193	mpi_status_ignore, ierr)
2194
2195	DO J=JTS,min(JTE,JDE-1)
2196	DO I=ITS,min(ITE,IDE-1)
2197	grid%hlat_gc(I,J)=dumdata(I,J,1)
2198	ENDDO
2199	ENDDO
2200
2201
2202	!!! GWD fields
2203
2204	varName='HSTDV'
2205	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2206	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2207	dumdata,hor_size,mpi_real4, &
2208	mpi_status_ignore, ierr)
2209
2210	DO J=JTS,min(JTE,JDE-1)
2211	DO I=ITS,min(ITE,IDE-1)
2212	grid%hstdv(I,J)=dumdata(I,J,1)
2213	ENDDO
2214	ENDDO
2215
2216	varName='HCNVX'
2217	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2218	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2219	dumdata,hor_size,mpi_real4, &
2220	mpi_status_ignore, ierr)
2221
2222	DO J=JTS,min(JTE,JDE-1)
2223	DO I=ITS,min(ITE,IDE-1)
2224	grid%hcnvx(I,J)=dumdata(I,J,1)
2225	ENDDO
2226	ENDDO
2227
2228	varName='HASYW'
2229	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2230	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2231	dumdata,hor_size,mpi_real4, &
2232	mpi_status_ignore, ierr)
2233
2234	DO J=JTS,min(JTE,JDE-1)
2235	DO I=ITS,min(ITE,IDE-1)
2236	grid%hasyw(I,J)=dumdata(I,J,1)
2237	ENDDO
2238	ENDDO
2239
2240	varName='HASYS'
2241	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2242	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2243	dumdata,hor_size,mpi_real4, &
2244	mpi_status_ignore, ierr)
2245
2246	DO J=JTS,min(JTE,JDE-1)
2247	DO I=ITS,min(ITE,IDE-1)
2248	grid%hasys(I,J)=dumdata(I,J,1)
2249	ENDDO
2250	ENDDO
2251
2252	varName='HASYSW'
2253	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2254	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2255	dumdata,hor_size,mpi_real4, &
2256	mpi_status_ignore, ierr)
2257
2258	DO J=JTS,min(JTE,JDE-1)
2259	DO I=ITS,min(ITE,IDE-1)
2260	grid%hasysw(I,J)=dumdata(I,J,1)
2261	ENDDO
2262	ENDDO
2263
2264	varName='HASYNW'
2265	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2266	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2267	dumdata,hor_size,mpi_real4, &
2268	mpi_status_ignore, ierr)
2269
2270	DO J=JTS,min(JTE,JDE-1)
2271	DO I=ITS,min(ITE,IDE-1)
2272	grid%hasynw(I,J)=dumdata(I,J,1)
2273	ENDDO
2274	ENDDO
2275
2276	varName='HLENW'
2277	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2278	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2279	dumdata,hor_size,mpi_real4, &
2280	mpi_status_ignore, ierr)
2281
2282	DO J=JTS,min(JTE,JDE-1)
2283	DO I=ITS,min(ITE,IDE-1)
2284	grid%hlenw(I,J)=dumdata(I,J,1)
2285	ENDDO
2286	ENDDO
2287
2288	varName='HLENS'
2289	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2290	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2291	dumdata,hor_size,mpi_real4, &
2292	mpi_status_ignore, ierr)
2293
2294	DO J=JTS,min(JTE,JDE-1)
2295	DO I=ITS,min(ITE,IDE-1)
2296	grid%hlens(I,J)=dumdata(I,J,1)
2297	ENDDO
2298	ENDDO
2299
2300	varName='HLENSW'
2301	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2302	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2303	dumdata,hor_size,mpi_real4, &
2304	mpi_status_ignore, ierr)
2305
2306	DO J=JTS,min(JTE,JDE-1)
2307	DO I=ITS,min(ITE,IDE-1)
2308	grid%hlensw(I,J)=dumdata(I,J,1)
2309	ENDDO
2310	ENDDO
2311
2312	varName='HLENNW'
2313	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2314	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2315	dumdata,hor_size,mpi_real4, &
2316	mpi_status_ignore, ierr)
2317
2318	DO J=JTS,min(JTE,JDE-1)
2319	DO I=ITS,min(ITE,IDE-1)
2320	grid%hlennw(I,J)=dumdata(I,J,1)
2321	ENDDO
2322	ENDDO
2323
2324	varName='HANGL'
2325	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2326	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2327	dumdata,hor_size,mpi_real4, &
2328	mpi_status_ignore, ierr)
2329
2330	DO J=JTS,min(JTE,JDE-1)
2331	DO I=ITS,min(ITE,IDE-1)
2332	grid%hangl(I,J)=dumdata(I,J,1)
2333	ENDDO
2334	ENDDO
2335
2336	varName='HANIS'
2337	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2338	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2339	dumdata,hor_size,mpi_real4, &
2340	mpi_status_ignore, ierr)
2341
2342	DO J=JTS,min(JTE,JDE-1)
2343	DO I=ITS,min(ITE,IDE-1)
2344	grid%hanis(I,J)=dumdata(I,J,1)
2345	ENDDO
2346	ENDDO
2347
2348	varName='HSLOP'
2349	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2350	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2351	dumdata,hor_size,mpi_real4, &
2352	mpi_status_ignore, ierr)
2353
2354	DO J=JTS,min(JTE,JDE-1)
2355	DO I=ITS,min(ITE,IDE-1)
2356	grid%hslop(I,J)=dumdata(I,J,1)
2357	ENDDO
2358	ENDDO
2359
2360	varName='HZMAX'
2361	CALL retrieve_index(index,VarName,varname_all,nrecs,iret)
2362	CALL mpi_file_read_at(iunit,file_offset(index+1), &
2363	dumdata,hor_size,mpi_real4, &
2364	mpi_status_ignore, ierr)
2365
2366	DO J=JTS,min(JTE,JDE-1)
2367	DO I=ITS,min(ITE,IDE-1)
2368	grid%hzmax(I,J)=dumdata(I,J,1)
2369	ENDDO
2370	ENDDO
2371	write(0,*) 'maxval grid%hzmax: ', maxval(grid%hzmax)
2372
2373	call mpi_file_close(mpi_comm_world, ierr)
2374
2375	varName='ST000010'
2376	flag_st000010 = 1
2377	num_st_levels_input = num_st_levels_input + 1
2378	st_levels_input(num_st_levels_input) = char2int2(varName(3:8))
2379	DO J=JTS,min(JTE,JDE-1)
2380	DO I=ITS,min(ITE,IDE-1)
2381	st_inputx(I,J,num_st_levels_input + 1) = grid%st000010(i,j)
2382	ENDDO
2383	ENDDO
2384
2385	varName='ST010040'
2386	flag_st010040 = 1
2387	num_st_levels_input = num_st_levels_input + 1
2388	st_levels_input(num_st_levels_input) = char2int2(varName(3:8))
2389	DO J=JTS,min(JTE,JDE-1)
2390	DO I=ITS,min(ITE,IDE-1)
2391	st_inputx(I,J,num_st_levels_input + 1) = grid%st010040(i,j)
2392	ENDDO
2393	ENDDO
2394
2395	varName='ST040100'
2396	flag_st040100 = 1
2397	num_st_levels_input = num_st_levels_input + 1
2398	st_levels_input(num_st_levels_input) = char2int2(varName(3:8))
2399	DO J=JTS,min(JTE,JDE-1)
2400	DO I=ITS,min(ITE,IDE-1)
2401	st_inputx(I,J,num_st_levels_input + 1) = grid%st040100(i,j)
2402	ENDDO
2403	ENDDO
2404
2405	varName='ST100200'
2406	flag_st100200 = 1
2407	num_st_levels_input = num_st_levels_input + 1
2408	st_levels_input(num_st_levels_input) = char2int2(varName(3:8))
2409	DO J=JTS,min(JTE,JDE-1)
2410	DO I=ITS,min(ITE,IDE-1)
2411	st_inputx(I,J,num_st_levels_input + 1) = grid%st100200(i,j)
2412	ENDDO
2413	ENDDO
2414
2415	varName='SM000010'
2416	flag_sm000010 = 1
2417	num_sm_levels_input = num_sm_levels_input + 1
2418	sm_levels_input(num_sm_levels_input) = char2int2(varName(3:8))
2419	DO J=JTS,min(JTE,JDE-1)
2420	DO I=ITS,min(ITE,IDE-1)
2421	sm_inputx(I,J,num_sm_levels_input + 1) = grid%sm000010(i,j)
2422	ENDDO
2423	ENDDO
2424
2425	varName='SM010040'
2426	flag_sm010040 = 1
2427	num_sm_levels_input = num_sm_levels_input + 1
2428	sm_levels_input(num_sm_levels_input) = char2int2(varName(3:8))
2429	DO J=JTS,min(JTE,JDE-1)
2430	DO I=ITS,min(ITE,IDE-1)
2431	sm_inputx(I,J,num_sm_levels_input + 1) = grid%sm010040(i,j)
2432	ENDDO
2433	ENDDO
2434
2435	varName='SM040100'
2436	flag_sm040100 = 1
2437	num_sm_levels_input = num_sm_levels_input + 1
2438	sm_levels_input(num_sm_levels_input) = char2int2(varName(3:8))
2439	DO J=JTS,min(JTE,JDE-1)
2440	DO I=ITS,min(ITE,IDE-1)
2441	sm_inputx(I,J,num_sm_levels_input + 1) = grid%sm040100(i,j)
2442	ENDDO
2443	ENDDO
2444
2445	varName='SM100200'
2446	flag_sm100200 = 1
2447	num_sm_levels_input = num_sm_levels_input + 1
2448	sm_levels_input(num_sm_levels_input) = char2int2(varName(3:8))
2449	DO J=JTS,min(JTE,JDE-1)
2450	DO I=ITS,min(ITE,IDE-1)
2451	sm_inputx(I,J,num_sm_levels_input + 1) = grid%sm100200(i,j)
2452	ENDDO
2453	ENDDO
2454
2455	flag_sst = 1
2456
2457	!new
2458	sw_inputx=0.
2459	!new
2460
2461	sw_input=0.
2462
2463
2464	do J=JTS,min(JDE-1,JTE)
2465	do K=1,num_st_levels_alloc
2466	do I=ITS,min(IDE-1,ITE)
2467	st_input(I,K,J)=st_inputx(I,J,K)
2468	sm_input(I,K,J)=sm_inputx(I,J,K)
2469	sw_input(I,K,J)=sw_inputx(I,J,K)
2470	enddo
2471	enddo
2472	enddo
2473
2474	! DEALLOCATE(pmsl)
2475	! DEALLOCATE(psfc_in)
2476	! DEALLOCATE(st_inputx)
2477	! DEALLOCATE(sm_inputx)
2478	! DEALLOCATE(sw_inputx)
2479	! DEALLOCATE(soilt010_input)
2480	! DEALLOCATE(soilt040_input)
2481	! DEALLOCATE(soilt100_input)
2482	! DEALLOCATE(soilt200_input)
2483	! DEALLOCATE(soilm010_input)
2484	! DEALLOCATE(soilm040_input)
2485	! DEALLOCATE(soilm100_input)
2486	! DEALLOCATE(soilm200_input)
2487	! DEALLOCATE(dumdata)
2488
2489	#endif
2490	end subroutine read_wps
2491
2492	! -------------------------------------------------------------------------
2493	! -------------------------------------------------------------------------
2494	! -------------------------------------------------------------------------
2495	! -------------------------------------------------------------------------
2496	!!!! MPI-IO pieces
2497
2498	subroutine retrieve_index(index,string,varname_all,nrecs,iret)
2499	!$$$ subprogram documentation block
2500	! . . . .
2501	! subprogram: retrieve_index get record number of desired variable
2502	! prgmmr: parrish org: np22 date: 2004-11-29
2503	!
2504	! abstract: by examining previously generated inventory of wrf binary restart file,
2505	! find record number that contains the header record for variable
2506	! identified by input character variable "string".
2507	!
2508	! program history log:
2509	! 2004-11-29 parrish
2510	!
2511	! input argument list:
2512	! string - mnemonic for variable desired
2513	! varname_all - list of all mnemonics obtained from inventory of file
2514	! nrecs - total number of sequential records counted in wrf
2515	! binary restart file
2516	!
2517	! output argument list:
2518	! index - desired record number
2519	! iret - return status, set to 0 if variable was found,
2520	! non-zero if not.
2521	!
2522	! attributes:
2523	! language: f90
2524	! machine: ibm RS/6000 SP
2525	!
2526	!$$$
2527
2528	implicit none
2529
2530	integer,intent(out)::iret
2531	integer,intent(in)::nrecs
2532	integer,intent(out):: index
2533	character(*), intent(in):: string
2534	character(132),intent(in)::varname_all(nrecs)
2535
2536	integer i
2537
2538	iret=0
2539
2540	do i=1,nrecs
2541	if(trim(string) == trim(varname_all(i))) then
2542	index=i
2543	return
2544	end if
2545	end do
2546
2547	write(6,*)' problem reading wrf nmm binary file, rec id "',trim(string),'" not found'
2548
2549	iret=-1
2550
2551	end subroutine retrieve_index
2552	subroutine next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2553	!$$$ subprogram documentation block
2554	! . . . .
2555	! subprogram: next_buf bring in next direct access block
2556	! prgmmr: parrish org: np22 date: 2004-11-29
2557	!
2558	! abstract: bring in next direct access block when needed, as the file is scanned
2559	! from beginning to end during counting and inventory of records.
2560	! (subroutines count_recs_wrf_binary_file and inventory_wrf_binary_file)
2561	!
2562	! program history log:
2563	! 2004-11-29 parrish
2564	!
2565	! input argument list:
2566	! in_unit - fortran unit number where input file is opened through.
2567	! nextbyte - byte number from beginning of file that is desired
2568	! locbyte - byte number from beginning of last block read for desired byt
2569	! lrecl - direct access block length
2570	! nreads - number of blocks read before now (for diagnostic information
2571	! lastbuf - logical, if true, then no more blocks, so return
2572	!
2573	! output argument list:
2574	! buf - output array containing contents of next block
2575	! locbyte - byte number from beginning of new block read for desired byte
2576	! thisblock - number of new block being read by this routine
2577	! nreads - number of blocks read now (for diagnostic information only)
2578	! lastbuf - logical, if true, then at end of file.
2579	!
2580	! attributes:
2581	! language: f90
2582	! machine: ibm RS/6000 SP
2583	!
2584	!$$$
2585	! use kinds, only: i_byte,i_llong
2586	implicit none
2587
2588	integer(i_llong) lrecl
2589	integer in_unit,nreads
2590	integer(i_byte) buf(lrecl)
2591	integer(i_llong) nextbyte,locbyte,thisblock
2592	logical lastbuf
2593
2594	integer ierr
2595
2596	if(lastbuf) return
2597
2598	ierr=0
2599	nreads=nreads+1
2600
2601	! compute thisblock:
2602
2603	thisblock = 1_i_llong + (nextbyte-1_i_llong)/lrecl
2604
2605	locbyte = 1_i_llong+mod(locbyte-1_i_llong,lrecl)
2606
2607	read(in_unit,rec=thisblock,iostat=ierr)buf
2608	lastbuf = ierr /= 0
2609
2610	end subroutine next_buf
2611
2612	subroutine inventory_wrf_binary_file(in_unit,wrf_ges_filename,nrecs, &
2613	datestr_all,varname_all,domainend_all, &
2614	start_block,end_block,start_byte,end_byte,file_offset)
2615	!$$$ subprogram documentation block
2616	! . . . .
2617	! subprogram: inventory_wrf_binary_file get contents of wrf binary file
2618	! prgmmr: parrish org: np22 date: 2004-11-29
2619	!
2620	! abstract: generate list of contents and map of wrf binary file which can be
2621	! used for reading and writing with mpi io routines.
2622	! same basic routine as count_recs_wrf_binary_file, except
2623	! now wrf unpacking routines are used to decode wrf header
2624	! records, and send back lists of variable mnemonics, dates,
2625	! grid dimensions, and byte addresses relative to start of
2626	! file for each field (this is used by mpi io routines).
2627	!
2628	! program history log:
2629	! 2004-11-29 parrish
2630	!
2631	!
2632	! input argument list:
2633	! in_unit - fortran unit number where input file is opened through.
2634	! wrf_ges_filename - filename of input wrf binary restart file
2635	! nrecs - number of sequential records found on input wrf binary restart file.
2636	! (obtained by a previous call to count_recs_wrf_binary_file)
2637	!
2638	! output argument list: (all following dimensioned nrecs)
2639	! datestr_all - date character string for each field, where applicable (or else blanks)
2640	! varname_all - wrf mnemonic for each variable, where applicable (or blank)
2641	! domainend_all - dimensions of each field, where applicable (up to 3 dimensions)
2642	! start_block - direct access block number containing 1st byte of record
2643	! (after 4 byte record mark)
2644	! end_block - direct access block number containing last byte of record
2645	! (before 4 byte record mark)
2646	! start_byte - relative byte address in direct access block of 1st byte of record
2647	! end_byte - relative byte address in direct access block of last byte of record
2648	! file_offset - absolute address of byte before 1st byte of record (used by mpi io)
2649	!
2650	! attributes:
2651	! language: f90
2652	! machine: ibm RS/6000 SP
2653	!
2654	!$$$
2655	! use kinds, only: r_single,i_byte,i_long,i_llong
2656	use module_internal_header_util
2657	implicit none
2658
2659	integer,intent(in)::in_unit,nrecs
2660	character(*),intent(in)::wrf_ges_filename
2661	character(132),intent(out)::datestr_all(nrecs),varname_all(nrecs)
2662	integer,intent(out)::domainend_all(3,nrecs)
2663	integer,intent(out)::start_block(nrecs),end_block(nrecs)
2664	integer,intent(out)::start_byte(nrecs),end_byte(nrecs)
2665	integer(i_llong),intent(out)::file_offset(nrecs)
2666
2667	integer irecs
2668	integer(i_llong) nextbyte,locbyte,thisblock
2669	integer(i_byte) lenrec4(4)
2670	integer(i_long) lenrec,lensave
2671	equivalence (lenrec4(1),lenrec)
2672	integer(i_byte) missing4(4)
2673	integer(i_long) missing
2674	equivalence (missing,missing4(1))
2675	integer(i_llong),parameter:: lrecl=2**16-1
2676	integer(i_byte) buf(lrecl)
2677	integer i,loc_count,nreads
2678	logical lastbuf
2679	integer(i_byte) hdrbuf4(2048)
2680	integer(i_long) hdrbuf(512)
2681	equivalence(hdrbuf(1),hdrbuf4(1))
2682	integer,parameter:: int_field = 530
2683	integer,parameter:: int_dom_ti_char = 220
2684	integer,parameter:: int_dom_ti_real = 140
2685	integer,parameter:: int_dom_ti_integer = 180
2686	integer hdrbufsize
2687	integer inttypesize
2688	integer datahandle,count
2689	character(128) element,dumstr,strdata
2690	integer loccode
2691	character(132) blanks
2692	integer typesize
2693	integer fieldtype,comm,iocomm
2694	integer domaindesc
2695	character(132) memoryorder,stagger,dimnames(3)
2696	integer domainstart(3),domainend(3)
2697	integer memorystart(3),memoryend(3)
2698	integer patchstart(3),patchend(3)
2699	character(132) datestr,varname
2700	real(r_single) dummy_field(1)
2701	! integer dummy_field
2702	! real dummy_field
2703	integer itypesize
2704	integer idata(1)
2705	real rdata(16)
2706
2707	call wrf_sizeof_integer(itypesize)
2708	inttypesize=itypesize
2709
2710	blanks=trim(' ')
2711
2712	open(in_unit,file=trim(wrf_ges_filename),access='direct',recl=lrecl)
2713	irecs=0
2714	missing=-9999
2715	nextbyte=0_i_llong
2716	locbyte=lrecl
2717	nreads=0
2718	lastbuf=.false.
2719	do
2720
2721	! get length of next record
2722
2723	do i=1,4
2724	nextbyte=nextbyte+1_i_llong
2725	locbyte=locbyte+1_i_llong
2726	if(locbyte > lrecl .and. lastbuf) go to 900
2727	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2728	lenrec4(i)=buf(locbyte)
2729	end do
2730	if(lenrec <= 0 .and. lastbuf) go to 900
2731	if(lenrec <= 0 .and. .not. lastbuf) go to 885
2732	nextbyte=nextbyte+1_i_llong
2733	locbyte=locbyte+1_i_llong
2734	if(locbyte > lrecl .and. lastbuf) go to 900
2735	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2736
2737	irecs=irecs+1
2738	start_block(irecs)=thisblock
2739	start_byte(irecs)=locbyte
2740	file_offset(irecs)=nextbyte-1_i_llong
2741	hdrbuf4(1)=buf(locbyte)
2742	hdrbuf4(2:4)=missing4(2:4)
2743	hdrbuf4(5:8)=missing4(1:4)
2744	datestr_all(irecs)=blanks
2745	varname_all(irecs)=blanks
2746	domainend_all(1:3,irecs)=0
2747
2748	loc_count=1
2749	do i=2,8
2750	if(loc_count.ge.lenrec) exit
2751	loc_count=loc_count+1
2752	nextbyte=nextbyte+1_i_llong
2753	locbyte=locbyte+1_i_llong
2754	if(locbyte > lrecl .and. lastbuf) go to 900
2755	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2756	hdrbuf4(i)=buf(locbyte)
2757	end do
2758
2759	! if(lenrec==2048) write(6,*)' irecs,hdrbuf(2),int_dom_ti_char,int_field=', &
2760	! irecs,hdrbuf(2),int_dom_ti_char,int_field
2761
2762	if(lenrec==2048.and.(hdrbuf(2) == int_dom_ti_char .or. hdrbuf(2) == int_field &
2763	.or. hdrbuf(2) == int_dom_ti_real .or. hdrbuf(2) == int_dom_ti_integer)) then
2764
2765	! bring in next full record, so we can unpack datestr, varname, and domainend
2766	do i=9,lenrec
2767	loc_count=loc_count+1
2768	nextbyte=nextbyte+1_i_llong
2769	locbyte=locbyte+1_i_llong
2770	if(locbyte > lrecl .and. lastbuf) go to 900
2771	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2772	hdrbuf4(i)=buf(locbyte)
2773	end do
2774
2775	if(hdrbuf(2) == int_dom_ti_char) then
2776
2777	call int_get_ti_header_char(hdrbuf,hdrbufsize,inttypesize, &
2778	datahandle,element,dumstr,strdata,loccode)
2779	varname_all(irecs)=trim(element)
2780	datestr_all(irecs)=trim(strdata)
2781	! write(6,*)' irecs,varname,datestr = ',irecs,trim(varname_all(irecs)),trim(datestr_all(irecs))
2782
2783	else if(hdrbuf(2) == int_dom_ti_real) then
2784
2785	call int_get_ti_header_real(hdrbuf,hdrbufsize,inttypesize,typesize, &
2786	datahandle,element,rdata,count,loccode)
2787	varname_all(irecs)=trim(element)
2788	! datestr_all(irecs)=trim(strdata)
2789	! write(6,*) 'count: ', count
2790	! write(6,*)' irecs,varname,datestr = ',irecs,trim(varname_all(irecs)),rdata(1:count)
2791
2792	else if(hdrbuf(2) == int_dom_ti_integer) then
2793
2794	call int_get_ti_header_integer(hdrbuf,hdrbufsize,inttypesize,typesize, &
2795	datahandle,element,idata,count,loccode)
2796	varname_all(irecs)=trim(element)
2797	! datestr_all(irecs)=trim(strdata)
2798	! write(6,*)' irecs,varname,datestr = ',irecs,trim(varname_all(irecs)),idata(1:count)
2799
2800	else
2801
2802	call int_get_write_field_header(hdrbuf,hdrbufsize,inttypesize,typesize, &
2803	datahandle,datestr,varname,dummy_field,fieldtype,comm,iocomm, &
2804	domaindesc,memoryorder,stagger,dimnames, &
2805	domainstart,domainend,memorystart,memoryend,patchstart,patchend)
2806	varname_all(irecs)=trim(varname)
2807	datestr_all(irecs)=trim(datestr)
2808	domainend_all(1:3,irecs)=domainend(1:3)
2809	! write(6,*)' irecs,datestr,domend,varname = ', &
2810	! irecs,trim(datestr_all(irecs)),domainend_all(1:3,irecs),trim(varname_all(irecs))
2811
2812	end if
2813	end if
2814
2815	nextbyte=nextbyte-loc_count+lenrec
2816	locbyte=locbyte-loc_count+lenrec
2817	if(locbyte > lrecl .and. lastbuf) go to 900
2818	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2819	end_block(irecs)=thisblock
2820	end_byte(irecs)=locbyte
2821	lensave=lenrec
2822	do i=1,4
2823	nextbyte=nextbyte+1_i_llong
2824	locbyte=locbyte+1_i_llong
2825	if(locbyte > lrecl .and. lastbuf) go to 900
2826	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2827	lenrec4(i)=buf(locbyte)
2828	end do
2829	if(lenrec /= lensave) go to 890
2830
2831	end do
2832
2833	880 continue
2834	write(6,*)' reached impossible place in inventory_wrf_binary_file'
2835	close(in_unit)
2836	return
2837
2838	885 continue
2839	write(6,*)' problem in inventory_wrf_binary_file, lenrec has bad value before end of file'
2840	write(6,*)' lenrec =',lenrec
2841	close(in_unit)
2842	return
2843
2844	890 continue
2845	write(6,*)' problem in inventory_wrf_binary_file, beginning and ending rec len words unequal'
2846	write(6,*)' begining reclen =',lensave
2847	write(6,*)' ending reclen =',lenrec
2848	write(6,*)' irecs =',irecs
2849	write(6,*)' nrecs =',nrecs
2850	close(in_unit)
2851	return
2852
2853	900 continue
2854	write(6,*)' normal end of file reached in inventory_wrf_binary_file'
2855	write(6,*)' nblocks=',thisblock
2856	write(6,*)' irecs,nrecs=',irecs,nrecs
2857	write(6,*)' nreads=',nreads
2858	close(in_unit)
2859
2860	end subroutine inventory_wrf_binary_file
2861
2862	SUBROUTINE wrf_sizeof_integer( retval )
2863	IMPLICIT NONE
2864	INTEGER retval
2865	! 4 is defined by CPP
2866	retval = 4
2867	RETURN
2868	END SUBROUTINE wrf_sizeof_integer
2869
2870	SUBROUTINE wrf_sizeof_real( retval )
2871	IMPLICIT NONE
2872	INTEGER retval
2873	! 4 is defined by CPP
2874	retval = 4
2875	RETURN
2876	END SUBROUTINE wrf_sizeof_real
2877	subroutine count_recs_wrf_binary_file(in_unit,wrf_ges_filename,nrecs)
2878	!$$$ subprogram documentation block
2879	! . . . .
2880	! subprogram: count_recs_binary_file count # recs on wrf binary file
2881	! prgmmr: parrish org: np22 date: 2004-11-29
2882	!
2883	! abstract: count number of sequential records contained in wrf binary
2884	! file. this is done by opening the file in direct access
2885	! mode with block length of 2**20, the size of the physical
2886	! blocks on ibm "blue" and "white" machines. for optimal
2887	! performance, change block length to correspond to the
2888	! physical block length of host machine disk space.
2889	! records are counted by looking for the 4 byte starting
2890	! and ending sequential record markers, which contain the
2891	! record size in bytes. only blocks are read which are known
2892	! by simple calculation to contain these record markers.
2893	! even though this is done on one processor, it is still
2894	! very fast, and the time will always scale by the number of
2895	! sequential records, not their size. this step and the
2896	! following inventory step consistently take less than 0.1 seconds
2897	! to complete.
2898	!
2899	! program history log:
2900	! 2004-11-29 parrish
2901	!
2902	! input argument list:
2903	! in_unit - fortran unit number where input file is opened through.
2904	! wrf_ges_filename - filename of input wrf binary restart file
2905	!
2906	! output argument list:
2907	! nrecs - number of sequential records found on input wrf binary restart fil
2908	!
2909	! attributes:
2910	! language: f90
2911	! machine: ibm RS/6000 SP
2912	!
2913	!$$$
2914
2915	! do an initial read through of a wrf binary file, and get total number of sequential fil
2916
2917	! use kinds, only: r_single,i_byte,i_long,i_llong
2918	implicit none
2919
2920	integer,intent(in)::in_unit
2921	character(*),intent(in)::wrf_ges_filename
2922	integer,intent(out)::nrecs
2923
2924	integer(i_llong) nextbyte,locbyte,thisblock
2925	integer(i_byte) lenrec4(4)
2926	integer(i_long) lenrec,lensave
2927	equivalence (lenrec4(1),lenrec)
2928	integer(i_byte) missing4(4)
2929	integer(i_long) missing
2930	equivalence (missing,missing4(1))
2931	integer(i_llong),parameter:: lrecl=2**16-1
2932	integer(i_byte) buf(lrecl)
2933	integer i,loc_count,nreads
2934	logical lastbuf
2935
2936	open(in_unit,file=trim(wrf_ges_filename),access='direct',recl=lrecl)
2937	nrecs=0
2938	missing=-9999
2939	nextbyte=0_i_llong
2940	locbyte=lrecl
2941	nreads=0
2942	lastbuf=.false.
2943	do
2944
2945	! get length of next record
2946
2947	do i=1,4
2948	nextbyte=nextbyte+1_i_llong
2949	locbyte=locbyte+1_i_llong
2950	if(locbyte > lrecl .and. lastbuf) go to 900
2951	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2952	lenrec4(i)=buf(locbyte)
2953	end do
2954	if(lenrec <= 0 .and. lastbuf) go to 900
2955	if(lenrec <= 0 .and. .not.lastbuf) go to 885
2956	nextbyte=nextbyte+1_i_llong
2957	locbyte=locbyte+1_i_llong
2958	if(locbyte > lrecl .and. lastbuf) go to 900
2959	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2960
2961	nrecs=nrecs+1
2962
2963	loc_count=1
2964	do i=2,4
2965	if(loc_count.ge.lenrec) exit
2966	loc_count=loc_count+1
2967	nextbyte=nextbyte+1_i_llong
2968	locbyte=locbyte+1_i_llong
2969	if(locbyte > lrecl .and. lastbuf) go to 900
2970	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2971	end do
2972	do i=1,4
2973	if(loc_count.ge.lenrec) exit
2974	loc_count=loc_count+1
2975	nextbyte=nextbyte+1_i_llong
2976	locbyte=locbyte+1_i_llong
2977	if(locbyte > lrecl .and. lastbuf) go to 900
2978	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2979	end do
2980	nextbyte=nextbyte-loc_count+lenrec
2981	locbyte=locbyte-loc_count+lenrec
2982	if(locbyte > lrecl .and. lastbuf) go to 900
2983	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2984	lensave=lenrec
2985	do i=1,4
2986	nextbyte=nextbyte+1_i_llong
2987	locbyte=locbyte+1_i_llong
2988	if(locbyte > lrecl .and. lastbuf) go to 900
2989	if(locbyte > lrecl) call next_buf(in_unit,buf,nextbyte,locbyte,thisblock,lrecl,nreads,lastbuf)
2990	lenrec4(i)=buf(locbyte)
2991	end do
2992	if(lenrec /= lensave) go to 890
2993
2994	end do
2995
2996	880 continue
2997	write(6,*)' reached impossible place in count_recs_wrf_binary_file'
2998	close(in_unit)
2999	return
3000
3001	885 continue
3002	write(6,*)' problem in count_recs_wrf_binary_file, lenrec has bad value before end of file'
3003	write(6,*)' lenrec =',lenrec
3004	close(in_unit)
3005	return
3006
3007	890 continue
3008	write(6,*)' problem in count_recs_wrf_binary_file, beginning and ending rec len words unequal'
3009	write(6,*)' begining reclen =',lensave
3010	write(6,*)' ending reclen =',lenrec
3011	close(in_unit)
3012	return
3013
3014	900 continue
3015	write(6,*)' normal end of file reached in count_recs_wrf_binary_file'
3016	write(6,*)' nblocks=',thisblock
3017	write(6,*)' nrecs=',nrecs
3018	write(6,*)' nreads=',nreads
3019	close(in_unit)
3020
3021	end subroutine count_recs_wrf_binary_file
3022
3023	subroutine retrieve_field(in_unit,wrfges,out,start_block,end_block,start_byte,end_byte)
3024	!$$$ subprogram documentation block
3025	! . . . .
3026	! subprogram: retrieve_field retrieve field from wrf binary file
3027	! prgmmr: parrish org: np22 date: 2004-11-29
3028	!
3029	! abstract: still using direct access, retrieve a field from the wrf binary restart file.
3030	!
3031	! program history log:
3032	! 2004-11-29 parrish
3033	!
3034	! input argument list:
3035	! in_unit - fortran unit number where input file is opened through.
3036	! wrfges - filename of input wrf binary restart file
3037	! start_block - direct access block number containing 1st byte of record
3038	! (after 4 byte record mark)
3039	! end_block - direct access block number containing last byte of record
3040	! (before 4 byte record mark)
3041	! start_byte - relative byte address in direct access block of 1st byte of record
3042	! end_byte - relative byte address in direct access block of last byte of record
3043	!
3044	! output argument list:
3045	! out - output buffer where desired field is deposited
3046	!
3047	! attributes:
3048	! language: f90
3049	! machine: ibm RS/6000 SP
3050	!
3051	!$$$
3052
3053	! use kinds, only: i_byte,i_kind
3054	implicit none
3055
3056	integer(i_kind),intent(in)::in_unit
3057	character(50),intent(in)::wrfges
3058	integer(i_kind),intent(in)::start_block,end_block,start_byte,end_byte
3059	integer(i_byte),intent(out)::out(*)
3060
3061	integer(i_llong),parameter:: lrecl=2**16-1
3062	integer(i_byte) buf(lrecl)
3063	integer(i_kind) i,ii,k
3064	integer ierr, ibegin, iend
3065
3066	open(in_unit,file=trim(wrfges),access='direct',recl=lrecl)
3067
3068	write(6,*)' in retrieve_field, start_block,end_block=',start_block,end_block
3069	write(6,*)' in retrieve_field, start_byte,end_byte=',start_byte,end_byte
3070	ii=0
3071	ierr=0
3072	do k=start_block,end_block
3073	read(in_unit,rec=k,iostat=ierr)buf
3074	ibegin=1 ; iend=lrecl
3075	if(k == start_block) ibegin=start_byte
3076	if(k == end_block) iend=end_byte
3077	do i=ibegin,iend
3078	ii=ii+1
3079	out(ii)=buf(i)
3080	end do
3081	end do
3082	close(in_unit)
3083
3084	end subroutine retrieve_field
3085
3086
3087	END MODULE module_si_io_nmm

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

Download in other formats:

Original Format