Context Navigation

process_domain_module_old.F90 @ 1068

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

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

Download in other formats:

Original Format