Context Navigation

process_domain_module_new.F90 @ 357

Last change on this file since 357 was 207, checked in by aslmd, 13 years ago
MESOSCALE: A GENERAL CLEAN-UP FOLLOWING UPDATING THE USER MANUAL. EVERYTHING ESSENTIAL IS IN MESOSCALE (much lighter than before). EVERYTHING FOR DEVELOPPERS OR EXPERTS IS IN MESOSCALE_DEV.
File size: 27.2 KB

Line
1	MODULE process_domain_module
2
3	USE input_module
4	USE output_module
5	USE module_model_basics
6	USE module_arrays
7
8	integer :: next_file
9	logical :: run_out_of_files
10	logical :: file_is_open
11	integer :: vertical !!! needed for v5dcreation
12
13
14	CONTAINS
15
16	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
17	! Name: process_domain
18	! Purpose: Process the input data
19	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20	SUBROUTINE process_domain ()
21
22	USE date_pack
23	USE gridinfo_module
24	USE misc_definitions_module
25	USE module_debug
26	USE v5d_module
27
28	IMPLICIT NONE
29
30	! ------ some vis5d stuff, contain MAXVARS,MAXTIMES,MAXROWS,MAXCOLUMNS,MAXLEVELS
31	include 'v5df.h'
32
33	! Local variables
34	integer :: idiff, n_times, i
35	character (len=19) :: valid_date, temp_date, datestr
36	character (len=128) :: cname, stagger, cunits, cdesc
37	integer :: istatus, reclength
38	integer :: ret !!! v5d stuff
39	integer, dimension(MAXLEVELS) :: v5d_nl
40	real, dimension(MAXLEVELS) :: v5d_levels
41	integer, dimension(MAXTIMES) :: timestamp, datestamp
42	integer :: projection
43	real, dimension(100) :: proj_args
44
45
46	! Compute number of times that we will process
47	CALL geth_idts(end_date, start_date, idiff)
48	CALL mprintf((idiff < 0),ERROR,'Ending date is earlier than starting date in namelist.')
49
50
51	! Check that the interval evenly divides the range of times to process
52	n_times = idiff / interval_seconds
53	CALL mprintf((mod(idiff, interval_seconds) /= 0),WARN, &
54	'In namelist, interval_seconds does not evenly divide '// &
55	'(end_date - start_date). Only %i time periods '// &
56	'will be processed.', i2=n_times)
57
58	!
59	! DO TIME_INDEPENDANT PROCESSING
60	!
61
62	! Initialize the input module to read input data
63	CALL input_init(1, istatus)
64	CALL mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input.')
65
66	! Read global attributes from the input file
67	CALL read_global_attrs ()
68	IF ( iprogram .le. 1 ) n_times = 0
69
70
71	IF (( output_type == 'grads' ) .OR. ( output_type == 'idl' )) THEN !!! We are making grads type files
72
73	! Open .dat and .ctl files
74	DO cunit=10,100
75	INQUIRE(unit=cunit, opened=is_used)
76	if (.not. is_used) exit
77	END DO
78	ctlfile = trim(output_root_name)//'.ctl'
79	OPEN(cunit,file=ctlfile)
80
81	DO dunit=10,100
82	INQUIRE(unit=dunit, opened=is_used)
83	IF (.not. is_used) EXIT
84	END DO
85
86	IF ( west_east_dim .gt. 2 .AND. south_north_dim .gt. 2 ) THEN
87	reclength = (west_east_dim)*(south_north_dim)
88	ELSE IF ( west_east_dim .gt. 2 .AND. south_north_dim .eq. 2 ) THEN
89	reclength = west_east_dim
90	ELSE IF ( west_east_dim .eq. 2 .AND. south_north_dim .gt. 2 ) THEN
91	reclength = south_north_dim
92	ELSE IF ( west_east_dim .eq. 2 .AND. south_north_dim .eq. 2 ) THEN
93	reclength = 1
94	END IF
95	#ifdef RECL4
96	reclength = reclength*4
97	#endif
98	datfile = trim(output_root_name)//'.dat'
99	! OPEN(dunit,file=datfile,form='unformatted',access="direct", recl=reclength)
100	OPEN(dunit,file=datfile)
101
102	ELSE IF ( output_type == 'v5d' ) THEN !!! We are making grads type files
103
104	v5dfile = trim(output_root_name)//'.v5d'
105
106	END IF
107
108
109	CALL arw_get_next_time(handle, datestr, istatus)
110
111	!!!
112	!!! MARS: au cas ou il n y a pas le geop init pour calculer les hauteurs
113	!!!
114	!!CALL arw_get_next_time(handle, datestr, istatus)
115
116
117
118	vertical = 0 !!! needed for v5d creation (equally spaced levels in generic units)
119	! If interp_method /= 0, then we need to know if this can be done, for now only doing it for wrfout data
120	IF ( iprogram .ge. 6 .AND. interp_method /= 0 ) THEN
121	CALL get_interp_info (datestr)
122	ELSE
123	vertical_type = 'n' !! In case user set this for other programs
124	number_of_zlevs = bottom_top_dim
125	END IF
126
127
128
129	#ifdef V5D
130	IF ( output_type == 'v5d' ) THEN !!! Create the vis5d file
131	v5d_nl = IMISSING
132	CALL v5d_varnames (v5d_nl)
133	CALL v5d_times (n_times+1, timestamp, datestamp)
134	CALL v5d_proj (datestr, projection, proj_args)
135	v5d_levels = MISSING
136	IF (vertical == 0 ) THEN
137	v5d_levels(1) = 0.0 !!! Same basics set up for sigma level data
138	v5d_levels(2) = 1.0
139	ELSE
140	DO i=1,number_of_zlevs
141	v5d_levels(i) = interp_levels(i)
142	END DO
143	END IF
144	ret = v5dCreate( v5dfile, n_times+1, numvars, &
145	!!west_east_dim, south_north_dim, v5d_nl, &
146	south_north_dim, west_east_dim, v5d_nl, &
147	varnames, timestamp, datestamp, 1, projection, &
148	proj_args, vertical, v5d_levels )
149	IF ( ret == 0 ) THEN
150	CALL mprintf(.true.,ERROR, 'v5dCreate')
151	ELSE
152	CALL mprintf(.true.,STDOUT, ' ')
153	CALL mprintf(.true.,STDOUT, ' SUCCESS: v5d file has been Created')
154	END IF
155
156	END IF
157	#endif
158
159	CALL input_close() !!! Will be opened again by next process if needed
160
161	!
162	! BEGIN TIME-DEPENDANT PROCESSING
163	next_file = 1
164	ivars = 0
165
166	! Loop over all times to be processed
167	could_not_find = plot_these_fields !!! Keep a list of requested fields
168	run_out_of_files = .FALSE.
169	file_is_open = .FALSE.
170	rec = 0
171	all_files : DO time=0,n_times
172
173	IF ( iprogram == 1 ) then
174	temp_date = "0000-00-00_00:00:00"
175	ELSE
176	CALL geth_newdate(valid_date, trim(start_date), time*interval_seconds)
177	temp_date = ' '
178	write(temp_date,'(a19)') valid_date(1:10)//'_'//valid_date(12:19)
179	ENDIF
180	IF ( time == 0 ) tdef_date = temp_date
181
182	CALL get_fields(temp_date)
183
184	IF (run_out_of_files) EXIT all_files
185
186	ENDDO all_files ! Loop over n_times
187
188	CALL mprintf(.true.,STDOUT, ' ')
189	CALL mprintf(.true.,STDOUT, 'DONE Processing Data')
190
191	IF ( output_type == 'grads' ) THEN !!! We are making grads type files
192	! Create the .ctl file
193	CALL mprintf(.true.,STDOUT, ' ')
194	CALL mprintf(.true.,STDOUT, 'CREATING .ctl file')
195	CALL create_ctl ()
196	END IF
197
198
199	#ifdef V5D
200	IF ( output_type == 'v5d' ) THEN !!! Create the vis5d file
201	ret = v5dclose()
202	END IF
203	#endif
204
205
206	END SUBROUTINE process_domain
207
208
209	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
210	! Name: get_interp_info
211	! Purpose: Get extra info if we need to interpolate to pressure / height
212	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
213	SUBROUTINE get_interp_info (valid_date)
214
215	IMPLICIT NONE
216
217	! Arguments
218	character (len=19) :: valid_date
219
220	! Local variables
221	integer :: istatus, wrftype
222	real, pointer, dimension(:,:,:) :: real_array
223	character (len=3) :: memorder
224	character (len=128) :: cname, stagger
225	character (len=128), dimension(3) :: dimnames
226	integer, dimension(3) :: domain_start, domain_end
227	integer, dimension(2) :: loc_of_min_z
228	real, allocatable, dimension(:,:,:) :: ph_tmp, phb_tmp
229	integer :: found
230
231
232	found = 0
233
234	wrftype = 104
235	memorder = "XYZ"
236	domain_start = 1
237	domain_end(1) = west_east_dim
238	domain_end(2) = south_north_dim
239
240	IF (iprogram == 6 .AND. vertical_type == 'p') THEN
241	domain_end(3) = bottom_top_dim
242	stagger = '-'
243	istatus = 0
244	CALL read_spec_field(domain_start, domain_end, 'QVAPOR', wrftype, &
245	memorder, stagger, dimnames, real_array, valid_date, istatus)
246	found = found + istatus
247
248	memorder = "XY"
249	domain_end(3) = 1
250	istatus = 0
251	CALL read_spec_field(domain_start, domain_end, 'MU', wrftype, &
252	memorder, stagger, dimnames, real_array, valid_date, istatus)
253	found = found + istatus
254
255	istatus = 0
256	CALL read_spec_field(domain_start, domain_end, 'MUB', wrftype, &
257	memorder, stagger, dimnames, real_array, valid_date, istatus)
258	found = found + istatus
259
260	memorder = "0"
261	domain_end(1) = 1
262	domain_end(2) = 1
263	istatus = 0
264	CALL read_spec_field(domain_start, domain_end, 'P_TOP', wrftype, &
265	memorder, stagger, dimnames, real_array, valid_date, istatus)
266	found = found + istatus
267
268	memorder = "Z"
269	domain_end(1) = bottom_top_dim
270	istatus = 0
271	CALL read_spec_field(domain_start, domain_end, 'ZNU', wrftype, &
272	memorder, stagger, dimnames, real_array, valid_date, istatus)
273	found = found + istatus
274
275	domain_end(1) = bottom_top_dim + 1
276	stagger = 'Z'
277	istatus = 0
278	CALL read_spec_field(domain_start, domain_end, 'ZNW', wrftype, &
279	memorder, stagger, dimnames, real_array, valid_date, istatus)
280	found = found + istatus
281
282	IF ( found == 0 ) THEN
283	vertical = 3 !!! (unequally spaced levels in mb - v5d)
284	CALL mprintf(.true.,STDOUT, ' Interpolating to PRESSURE levels ')
285	CALL mprintf(.true.,STDOUT, ' ')
286	ELSE
287	CALL mprintf(.true.,STDOUT, ' WARNING: Asked to interpolate to PRESSURE, but we do not have enough information ' )
288	CALL mprintf(.true.,STDOUT, ' Will output data on MODEL LEVELS ' )
289	vertical_type = 'n'
290	ENDIF
291
292	ENDIF
293
294	IF (iprogram == 8 .AND. vertical_type == 'p') THEN
295	domain_end(3) = bottom_top_dim
296	stagger = '-'
297
298	istatus = 0
299	CALL read_spec_field(domain_start, domain_end, 'P', wrftype, &
300	memorder, stagger, dimnames, real_array, valid_date, istatus)
301	found = found + istatus
302
303	istatus = 0
304	CALL read_spec_field(domain_start, domain_end, 'PB', wrftype, &
305	memorder, stagger, dimnames, real_array, valid_date, istatus)
306	found = found + istatus
307
308	IF ( found == 0 ) THEN
309	vertical = 3 !!! (unequally spaced levels in mb - v5d)
310	CALL mprintf(.true.,STDOUT, ' Interpolating to PRESSURE levels ')
311	CALL mprintf(.true.,STDOUT, ' ')
312	ELSE
313	CALL mprintf(.true.,STDOUT, ' WARNING: Asked to interpolate to PRESSURE, but we do not have enough information ' )
314	CALL mprintf(.true.,STDOUT, ' Will output data on MODEL LEVELS ' )
315	vertical_type = 'n'
316	ENDIF
317
318	ENDIF
319
320
321	IF (vertical_type == 'z' .AND. interp_method == 1 ) THEN
322
323	domain_end(3) = bottom_top_dim + 1
324	stagger = 'Z'
325
326	!domain_end(3) = bottom_top_dim
327	!stagger = ''
328
329	! istatus = 0
330	! CALL read_spec_field(domain_start, domain_end, 'PH', wrftype, &
331	! memorder, stagger, dimnames, real_array, valid_date, istatus)
332	! found = found + istatus
333	!
334	! istatus = 0
335	! CALL read_spec_field(domain_start, domain_end, 'PHB', wrftype, &
336	! memorder, stagger, dimnames, real_array, valid_date, istatus)
337	! found = found + istatus
338
339	!!
340	!! MARS MARS
341	!!
342	!!print *, 'get_interp_info'
343	istatus = 0
344	CALL read_spec_field(domain_start, domain_end, 'PHTOT', wrftype, &
345	memorder, stagger, dimnames, real_array, valid_date, istatus)
346	found = found + istatus
347
348	IF ( found == 0 ) THEN
349	vertical = 2 !!! (unequally spaced levels in km - v5d)
350	CALL mprintf(.true.,STDOUT, ' Interpolating to USER SPECIFIED HEIGHT levels ')
351	CALL mprintf(.true.,STDOUT, ' ')
352	ELSE
353	CALL mprintf(.true.,STDOUT, ' WARNING: Asked to interpolate&
354	& to USER SPECIFIED HEIGHT, but we do not have enough information ' )
355	CALL mprintf(.true.,STDOUT, ' Will output data on MODEL LEVELS ' )
356	vertical_type = 'n'
357	ENDIF
358
359	ENDIF
360
361	IF (vertical_type == 'z' .AND. interp_method == -1 ) THEN
362	domain_end(3) = bottom_top_dim + 1
363	stagger = 'Z'
364
365	istatus = 0
366	CALL read_spec_field(domain_start, domain_end, 'PH', wrftype, &
367	memorder, stagger, dimnames, real_array, valid_date, istatus)
368	found = found + istatus
369	IF ( istatus == 0 ) THEN
370	allocate(ph_tmp(domain_end(1),domain_end(2),domain_end(3)-1))
371	ph_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3)))
372	ENDIF
373
374	istatus = 0
375	CALL read_spec_field(domain_start, domain_end, 'PHB', wrftype, &
376	memorder, stagger, dimnames, real_array, valid_date, istatus)
377	found = found + istatus
378	IF ( istatus == 0 ) THEN
379	allocate(phb_tmp(domain_end(1),domain_end(2),domain_end(3)-1))
380	phb_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3)))
381	ENDIF
382
383	IF ( found == 0 ) THEN
384	vertical = 2 !!! (unequally spaced levels in km - v5d)
385	!! ph_tmp = ( (ph_tmp+phb_tmp) / 9.81)/1000. !!! convert to height in km
386
387	ph_tmp = ( (ph_tmp+phb_tmp) / 3.72)/1000. !!! convert to height in km
388
389	!!! Generate nice heights.
390	!!! Get location of lowest height on model grid
391	!!! Use column above this point as our heights to interpolate to
392	!!! Adjust lowerst and heighest heights a bit
393	number_of_zlevs = bottom_top_dim
394	loc_of_min_z = minloc(ph_tmp(:,:,1))
395	interp_levels(1:number_of_zlevs) = &
396	ph_tmp(loc_of_min_z(1),loc_of_min_z(2),1:number_of_zlevs)
397
398	!!--- .002 problematic on Mars with low pressure
399	!!interp_levels(1) = interp_levels(1) + 0.002
400	!!interp_levels(1) = MAX(interp_levels(1), interp_levels(2)/2.0) !! no neg values
401	!!interp_levels(number_of_zlevs) = interp_levels(number_of_zlevs) - 0.002
402	vertical = 2
403	CALL mprintf(.true.,STDOUT, ' Interpolating to GENERATED HEIGHT levels ')
404	CALL mprintf(.true.,STDOUT, ' ')
405	ELSE
406	CALL mprintf(.true.,STDOUT, ' WARNING: Asked to interpolate to GENERATED &
407	&HEIGHT, but we do not have enough information ' )
408	CALL mprintf(.true.,STDOUT, ' Will output data on MODEL LEVELS ' )
409	vertical_type = 'n'
410	ENDIF
411
412	IF (ALLOCATED(ph_tmp)) DEALLOCATE(ph_tmp)
413	IF (ALLOCATED(phb_tmp)) DEALLOCATE(phb_tmp)
414	ENDIF
415
416
417	END SUBROUTINE get_interp_info
418
419	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
420	! Name: get_interp_array
421	! Purpose: Get array we will use when interpolting
422	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
423	SUBROUTINE get_interp_array (valid_date)
424
425	USE module_pressure
426
427	IMPLICIT NONE
428
429	! Arguments
430	character (len=19) :: valid_date
431
432	! Local variables
433	integer :: istatus, wrftype
434	real, pointer, dimension(:,:,:) :: real_array
435	character (len=3) :: memorder
436	character (len=128) :: cname, stagger
437	character (len=128), dimension(3) :: dimnames
438	integer, dimension(3) :: domain_start, domain_end
439	integer, dimension(2) :: loc_of_min_z
440	real, allocatable, dimension(:,:,:) :: p_tmp, pb_tmp, ph_tmp, phb_tmp
441	real, pointer, dimension(:,:,:) :: SCR3
442
443
444	wrftype = 104
445	memorder = "XYZ"
446	domain_start = 1
447	domain_end(1) = west_east_dim
448	domain_end(2) = south_north_dim
449
450	IF (iprogram == 6 .AND. vertical_type == 'p') THEN
451	domain_end(3) = bottom_top_dim
452	stagger = '-'
453	CALL read_spec_field(domain_start, domain_end, 'QVAPOR', wrftype, &
454	memorder, stagger, dimnames, real_array, valid_date, istatus)
455	allocate(QV(west_east_dim,south_north_dim,bottom_top_dim))
456	QV = real_array
457
458	memorder = "XY"
459	domain_end(3) = 1
460	stagger = '-'
461	CALL read_spec_field(domain_start, domain_end, 'MU', wrftype, &
462	memorder, stagger, dimnames, real_array, valid_date, istatus)
463	allocate(MU(west_east_dim,south_north_dim))
464	MU = real_array(:,:,1)
465	CALL read_spec_field(domain_start, domain_end, 'MUB', wrftype, &
466	memorder, stagger, dimnames, real_array, valid_date, istatus)
467	allocate(MUB(west_east_dim,south_north_dim))
468	MUB = real_array(:,:,1)
469
470	memorder = "0"
471	domain_end(1) = 1
472	domain_end(2) = 1
473	CALL read_spec_field(domain_start, domain_end, 'P_TOP', wrftype, &
474	memorder, stagger, dimnames, real_array, valid_date, istatus)
475	PTOP = real_array(1,1,1)
476
477	memorder = "Z"
478	domain_end(1) = bottom_top_dim
479	CALL read_spec_field(domain_start, domain_end, 'ZNU', wrftype, &
480	memorder, stagger, dimnames, real_array, valid_date, istatus)
481	allocate(ZNU(bottom_top_dim))
482	ZNU = real_array(:,1,1)
483
484	domain_end(1) = bottom_top_dim + 1
485	stagger = 'Z'
486	CALL read_spec_field(domain_start, domain_end, 'ZNW', wrftype, &
487	memorder, stagger, dimnames, real_array, valid_date, istatus)
488	allocate(ZNW(bottom_top_dim+1))
489	ZNW = real_array(:,1,1)
490
491
492	CALL pressure(SCR3)
493	ALLOCATE(vert_array(west_east_dim,south_north_dim,bottom_top_dim))
494	vert_array = 0.01 * SCR3 !! pressure array in HPa
495
496	IF (ALLOCATED(QV)) DEALLOCATE(QV)
497	IF (ALLOCATED(MU)) DEALLOCATE(MU)
498	IF (ALLOCATED(MUB)) DEALLOCATE(MUB)
499	IF (ALLOCATED(ZNU)) DEALLOCATE(ZNU)
500	IF (ALLOCATED(ZNW)) DEALLOCATE(ZNW)
501	ENDIF
502
503
504
505	IF (iprogram == 8 .AND. vertical_type == 'p') THEN
506	domain_end(3) = bottom_top_dim
507	stagger = '-'
508	CALL read_spec_field(domain_start, domain_end, 'P', wrftype, &
509	memorder, stagger, dimnames, real_array, valid_date, istatus)
510	allocate(p_tmp(west_east_dim,south_north_dim,bottom_top_dim))
511	p_tmp = real_array
512	CALL read_spec_field(domain_start, domain_end, 'PB', wrftype, &
513	memorder, stagger, dimnames, real_array, valid_date, istatus)
514	allocate(pb_tmp(west_east_dim,south_north_dim,bottom_top_dim))
515	pb_tmp = real_array
516
517	IF (ALLOCATED(vert_array)) DEALLOCATE(vert_array)
518	ALLOCATE(vert_array(west_east_dim,south_north_dim,bottom_top_dim))
519	vert_array = (p_tmp+pb_tmp)/100. !! pressure array in HPa
520
521	IF (ALLOCATED(p_tmp)) DEALLOCATE(p_tmp)
522	IF (ALLOCATED(pb_tmp)) DEALLOCATE(pb_tmp)
523	ENDIF
524
525
526
527	IF (vertical_type == 'z' ) THEN
528
529	domain_end(3) = bottom_top_dim + 1
530	stagger = 'Z'
531
532	!domain_end(3) = bottom_top_dim
533	!stagger = ''
534
535	!!!
536	!!! MARS MARS
537	!!!
538
539	! CALL read_spec_field(domain_start, domain_end, 'PH', wrftype, &
540	! memorder, stagger, dimnames, real_array, valid_date, istatus)
541	! allocate(ph_tmp(west_east_dim,south_north_dim,bottom_top_dim))
542	! ph_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3)))
543	! CALL read_spec_field(domain_start, domain_end, 'PHB', wrftype, &
544	! memorder, stagger, dimnames, real_array, valid_date, istatus)
545	! allocate(phb_tmp(west_east_dim,south_north_dim,bottom_top_dim))
546	! phb_tmp = 0.5*(real_array(:,:,1:domain_end(3)-1)+real_array(:,:,2:domain_end(3)))
547
548
549	CALL read_spec_field(domain_start, domain_end, 'PHTOT', wrftype, &
550	memorder, stagger, dimnames, real_array, valid_date, istatus)
551
552	allocate(ph_tmp(west_east_dim,south_north_dim,bottom_top_dim))
553	ph_tmp = 0.5*real_array
554
555	allocate(phb_tmp(west_east_dim,south_north_dim,bottom_top_dim))
556	phb_tmp = 0.5*real_array
557
558
559	IF (ALLOCATED(vert_array)) DEALLOCATE(vert_array)
560	ALLOCATE(vert_array(west_east_dim,south_north_dim,bottom_top_dim))
561	!! vert_array = ( (ph_tmp+phb_tmp) / 9.81)/1000. !! height array in km
562	vert_array = ( (ph_tmp+phb_tmp) / 3.72)/1000. !! height array in km
563
564	IF (ALLOCATED(ph_tmp)) DEALLOCATE(ph_tmp)
565	IF (ALLOCATED(phb_tmp)) DEALLOCATE(phb_tmp)
566	ENDIF
567
568
569	END SUBROUTINE get_interp_array
570
571	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
572	! Name: get_fields
573	! Purpose: Read all fields in input file and process required output fields
574	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
575	SUBROUTINE get_fields (valid_date)
576
577	USE module_get_file_names
578	USE module_diagnostics
579	USE module_interp
580
581	IMPLICIT NONE
582
583	! Local variables
584	integer :: istatus, i, j, k, ii, jj, kk
585	real, pointer, dimension(:,:,:) :: real_array
586	character (len=3) :: memorder
587	character (len=19) :: datestr, valid_date
588	character (len=128) :: cname, stagger, cunits, cdesc
589	character (len=128), dimension(3) :: dimnames
590	integer, dimension(3) :: domain_start, domain_end
591	real, pointer, dimension(:,:,:) :: data_out
592	real, allocatable, dimension(:,:,:):: SCR
593	integer :: nxout, nyout, nzout
594	character (len=20) :: dummy
595	integer :: is_there
596
597
598	! Initialize the input module to read static input data for this domain
599	CALL mprintf(.true.,STDOUT, ' ')
600	CALL mprintf(.true.,STDOUT, ' Processing time --- %s', s1=trim(valid_date))
601	get_right_time : DO
602	IF ( next_file > number_of_input_files ) THEN
603	CALL mprintf(.true.,STDOUT, ' --- WE HAVE RUN OUT OF INPUT FILES ---')
604	run_out_of_files = .TRUE.
605	RETURN
606	END IF
607	IF ( .not. file_is_open) THEN
608	CALL input_init(next_file, istatus)
609	CALL mprintf((istatus /= 0),ERROR, 'input_init(): Error opening input.')
610	file_is_open = .TRUE.
611	END IF
612	CALL arw_get_next_time(handle, datestr, istatus)
613	IF ( istatus /= 0 ) THEN ! might be in a next file
614	CALL input_close()
615	file_is_open = .FALSE.
616	CALL mprintf(.true.,STDOUT, ' Date not in this file - see if there are more files ')
617	next_file = next_file + 1
618	CYCLE get_right_time
619	END IF
620	IF ( TRIM(datestr) .LT. TRIM(valid_date) ) THEN
621	CALL mprintf(.true.,STDOUT, ' Not looking for %s - skipping to next time' , s1=datestr )
622	CYCLE get_right_time
623	ELSEIF ( TRIM(datestr) .EQ. TRIM(valid_date) ) THEN
624	CALL mprintf(.true.,STDOUT, ' Found the right date - continue ' )
625	EXIT get_right_time
626	ELSEIF ( TRIM(datestr) .GT. TRIM(valid_date) ) THEN
627	CALL mprintf(.true.,STDOUT, ' Found %s before %s', s1=trim(datestr),s2=trim(valid_date))
628	run_out_of_files = .TRUE.
629	RETURN
630	ENDIF
631	ENDDO get_right_time
632
633	CALL clobber_arrays
634	#ifdef IO_GRIB1
635	IF (io_form_input == GRIB1) THEN
636	REWIND (13)
637	END IF
638	#endif
639
640
641	!! IF we are interpolting we need the pressure/height array to interpolate to
642	IF ( vertical_type /= 'n' ) THEN
643	CALL get_interp_array (valid_date)
644	ENDIF
645
646
647	! Read fields using the input module; we know that there are no more
648	! fields to be read when read_next_field() returns a non-zero status.
649	istatus = 0
650	process_all_fields : DO WHILE (istatus == 0)
651	CALL read_next_field(domain_start, domain_end, cname, cunits, cdesc, &
652	memorder, stagger, dimnames, real_array, valid_date, istatus)
653	#ifdef IO_GRIB1
654	IF (io_form_input == GRIB1) THEN
655	IF (istatus == -5) THEN
656	istatus = 0
657	CYCLE process_all_fields
658	END IF
659	END IF
660	#endif
661	IF (istatus == 0) THEN
662	IF (cname(len_trim(cname)-1:len_trim(cname)) == "_U" ) CYCLE process_all_fields
663	IF (cname(len_trim(cname)-1:len_trim(cname)) == "_V" ) CYCLE process_all_fields
664
665	CALL keep_arrays(cname, real_array) !!! Keep some fields around
666
667	!!
668	!! DECIDE IF WE WANT THE FIELD AND IF YES, INTERPOLATE TO CORRECT GRID IF NEEEDED
669	!! assume we got the grid info up at the top somewhere
670	!!
671
672	dummy = ","//trim(cname)//","
673	is_there = INDEX(plot_these_fields,trim(dummy))
674	IF ( ( INDEX(plot,'all') /= 0 .OR. is_there /= 0) .AND. domain_end(2) /= 1) THEN
675	IF (ALLOCATED(SCR)) DEALLOCATE(SCR)
676	ALLOCATE(SCR(domain_end(1),domain_end(2),domain_end(3)))
677	SCR = real_array
678	CALL interp( SCR, domain_end(1), domain_end(2), domain_end(3), &
679	data_out, nxout, nyout, nzout, &
680	vert_array, interp_levels, number_of_zlevs)
681
682	! Write the fields we want out to the .dat file, also keeps a list of what is written out
683	CALL write_dat (data_out, nxout, nyout, nzout, cname, cdesc, cunits)
684	ENDIF
685
686	ENDIF !! end "istatus==0"
687
688	END DO process_all_fields
689
690
691
692	IF ( INDEX(plot,'list') /= 0 .OR. INDEX(plot,'file') /= 0 ) THEN
693	!! Do we have any DIAGNOSTICS to process?
694	CALL process_diagnostics ()
695	END IF
696
697
698	!! We are DONE for this time
699
700
701
702
703	!! print a list of the requested fields we could not find
704	IF ( len_trim(could_not_find) > 1 ) THEN
705	CALL mprintf(.true.,STDOUT, ' ')
706	CALL mprintf(.true.,STDOUT, ' WARNING: The following requested fields could not be found ' )
707
708	could_not_find = could_not_find(2:len_trim(could_not_find))
709	is_there = INDEX(could_not_find,",")
710	#ifdef V5D
711	IF ( output_type == 'v5d' .AND. is_there > 1 ) THEN !!! Warning for possible bad data file
712	CALL mprintf(.true.,STDOUT, ' WARNING: This will create a problem for VIS5D, please ' )
713	CALL mprintf(.true.,STDOUT, ' remove these fields from your input list and ' )
714	CALL mprintf(.true.,STDOUT, ' run the code again. ' )
715	CALL mprintf(.true.,STDOUT, ' ' )
716	END IF
717	#endif
718	DO WHILE ( is_there > 1 )
719	PRINT*," ", could_not_find(1:is_there-1)
720	could_not_find = could_not_find(is_there+1:len_trim(could_not_find))
721	is_there = INDEX(could_not_find,",")
722	END DO
723
724	CALL mprintf(.true.,STDOUT, ' ')
725	ENDIF
726
727
728	END SUBROUTINE get_fields
729
730	END MODULE process_domain_module

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

Download in other formats:

Original Format