Context Navigation

process_domain_module.F90 @ 777

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

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

Download in other formats:

Original Format