Context Navigation

api.F90 @ 189

Last change on this file since 189 was 186, checked in by aslmd, 14 years ago
MESOSCALE: python post-processing. wrapper with my Fortran interpolator (small modifications done to api.F90). added the corresponding options to winds.py which is now a pretty complete script
File size: 71.3 KB

Line
1	PROGRAM api
2
3	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4	! API !
5	! Altitude and Pressure Interpolator !
6	! !
7	! This program is based on !
8	! ------------------------ !
9	! p_interp v1.0 !
10	! http://www.mmm.ucar.edu/wrf/src/p_interp.tar.gz !
11	! September 2008 - Cindy Bruyere [NCAR, USA] !
12	! !
13	! Modifications !
14	! ------------- !
15	! - Presentation of the program (cosmetics) !
16	! - Additional routine and arguments in existing routines !
17	! - Improve memory management !
18	! - Martian constants + Addition of 'grav' variable !
19	! - Interpolation to height above areoid (MOLA zero datum) !
20	! October 2009 - Aymeric Spiga [The Open University, UK] !
21	! - Interpolation to height above surface !
22	! - Rotated winds (e.g. for polar projections) !
23	! November 2009 - AS !
24	! !
25	! Purpose !
26	! ------- !
27	! Program to read wrfout data (possibly several files) !
28	! and interpolate to pressure or height levels !
29	! A new NETCDF file is generated with appropriate suffix !
30	! The program reads namelist.api !
31	! !
32	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
33
34	IMPLICIT NONE
35	INCLUDE 'netcdf.inc'
36
37	!
38	! VARIABLES
39	!
40	CHARACTER (LEN=500) :: path_to_input
41	CHARACTER (LEN=500) :: path_to_output
42	CHARACTER (LEN=500) :: input_name
43	CHARACTER (LEN=500) :: output_name
44	CHARACTER (LEN=20) :: process
45	CHARACTER (LEN=2000) :: fields
46	REAL, DIMENSION(299) :: interp_levels
47	INTEGER :: interp_method=1
48	INTEGER :: extrapolate=0
49	LOGICAL :: debug=.FALSE.
50	LOGICAL :: unstagger_grid=.FALSE.
51	LOGICAL :: bit64=.FALSE.
52	LOGICAL :: oldvar=.TRUE.
53
54	INTEGER :: funit,ios
55	LOGICAL :: is_used
56
57	!
58	! NAMELISTS
59	!
60	NAMELIST /io/ path_to_input, input_name, path_to_output, output_name, &
61	process, fields, debug, bit64, oldvar
62	NAMELIST /interp_in/ interp_levels, interp_method, extrapolate, unstagger_grid
63
64	!
65	! DEFAULT VALUES for VARIABLES
66	!
67	path_to_input = './'
68	path_to_output = './'
69	output_name = ' '
70	interp_levels = -99999.
71	process = 'all'
72
73
74	!
75	! READ NAMELIST
76	!
77	DO funit=10,100
78	INQUIRE(unit=funit, opened=is_used)
79	IF (.not. is_used) EXIT
80	END DO
81	OPEN(funit,file='namelist.api',status='old',form='formatted',iostat=ios)
82	IF ( ios /= 0 ) STOP "ERROR opening namelist.api"
83	READ(funit,io)
84	READ(funit,interp_in)
85	CLOSE(funit)
86
87	!!! MAIN CALL
88	CALL api_main ( path_to_input, input_name, path_to_output, output_name, &
89	process, fields, debug, bit64, oldvar, &
90	interp_levels, interp_method, extrapolate, unstagger_grid, -99999. )
91
92	END PROGRAM api
93
94	SUBROUTINE api_main ( path_to_input, input_name, path_to_output, output_name, &
95	process, fields, debug, bit64, oldvar, &
96	interp_levels, interp_method, extrapolate, unstagger_grid, onelevel )
97
98	IMPLICIT NONE
99	INCLUDE 'netcdf.inc'
100
101	!!
102	!! EARTH CONSTANTS
103	!!
104	!REAL, PARAMETER :: Rd = 287.04 ! gas constant m2 s-2 K-1
105	!REAL, PARAMETER :: Cp = 7.Rd/2. ! r=8.314511E+01000.E+0/mugaz
106	!REAL, PARAMETER :: RCP = Rd/Cp
107	!REAL, PARAMETER :: p0 = 100000.
108	!REAL, PARAMETER :: grav = 9.81
109	!REAL, PARAMETER :: tpot0 = 300.
110
111	!
112	! MARS CONSTANTS
113	!
114	!REAL, PARAMETER :: Rd = 192. ! gas constant m2 s-2 K-1
115	!REAL, PARAMETER :: Cp = 844.6 ! r= 8.314511E+0*1000.E+0/mugaz
116	REAL, PARAMETER :: Rd = 191.0
117	REAL, PARAMETER :: Cp = 744.5
118	REAL, PARAMETER :: RCP = Rd/Cp
119	REAL, PARAMETER :: p0 = 610.
120	REAL, PARAMETER :: grav = 3.72
121	REAL, PARAMETER :: tpot0 = 220. ! reference potential temperature
122
123	!
124	! VARIABLES
125	!
126	CHARACTER, ALLOCATABLE, DIMENSION(:,:,:,:) :: text
127	CHARACTER (LEN=31),ALLOCATABLE, DIMENSION(:) :: dnamei, dnamej
128	CHARACTER(LEN=500),ALLOCATABLE, DIMENSION(:) :: input_file_names
129	CHARACTER(LEN=500),ALLOCATABLE, DIMENSION(:) :: output_file_names
130	DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:,:) :: ddata1, ddata2
131	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: data1, data2, data3
132	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: pres_field, pres_out
133	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: pres_stagU, pres_stagV
134	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: ght, phb, qv, tk, rh, tpot, u, v, umet, vmet
135	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: psfc
136	REAL, ALLOCATABLE, DIMENSION(:,:) :: ter
137	REAL, ALLOCATABLE, DIMENSION(:,:) :: longi, lati
138	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: longit, latit
139	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: interm1, interm2
140	INTEGER, ALLOCATABLE, DIMENSION(:) :: dvali, dvalj
141	INTEGER, ALLOCATABLE, DIMENSION(:,:,:,:) :: idata1, idata2
142	INTEGER, DIMENSION(4) :: start_dims = 1
143	INTEGER, DIMENSION(4) :: dims_in, dims_out
144	INTEGER, DIMENSION(6) :: ishape, jshape
145	CHARACTER (LEN=500) :: cval !80
146	CHARACTER (LEN=31) :: cname, test_dim_name
147	CHARACTER (LEN=500) :: input_file, output_file, att_text !80
148	CHARACTER (LEN=500) :: path_to_input
149	CHARACTER (LEN=500) :: path_to_output
150	CHARACTER (LEN=500) :: input_name
151	CHARACTER (LEN=500) :: output_name, tmp_name
152	CHARACTER (LEN=10) :: option
153	CHARACTER (LEN=132) :: command
154	CHARACTER (LEN=20) :: process, dummy
155	CHARACTER (LEN=2000) :: fields, process_these_fields
156	REAL, DIMENSION(299) :: interp_levels
157	REAL :: rval
158	REAL :: MISSING=1.e36
159	REAL :: truelat1, truelat2, stand_lon
160	INTEGER :: map_proj
161	INTEGER :: LINLOG = 1
162	INTEGER :: interp_method!=1
163	INTEGER :: extrapolate!=0
164	INTEGER :: ncid, mcid, rcode
165	INTEGER :: idm, ndims, nvars, natt, ngatts
166	INTEGER :: nunlimdimid
167	INTEGER :: i, ii, j, jj, ix, iy, iweg, isng, ibtg
168	INTEGER :: ivar, jvar
169	INTEGER :: times_in_file
170	INTEGER :: ilen, itype, ival, na, funit, ios
171	INTEGER :: num_metgrid_levels
172	INTEGER :: number_of_input_files
173	INTEGER :: ierr, loop, loslen, strlen, lent
174	INTEGER :: is_there
175	INTEGER :: kk
176	LOGICAL :: is_used
177	LOGICAL :: debug!=.FALSE.
178	LOGICAL :: interpolate=.FALSE.
179	LOGICAL :: unstagger_grid!=.FALSE.
180	LOGICAL :: fix_meta_stag=.FALSE.
181	LOGICAL :: bit64!=.FALSE.
182	LOGICAL :: first=.TRUE.
183	LOGICAL :: oldvar!=.FALSE.
184
185	REAL :: onelevel
186	if ( onelevel .ne. -99999. ) then
187	interp_levels(1) = onelevel
188	interp_levels(2:) = -99999.
189	endif
190
191	!
192	! INPUT FILE NAMES
193	!
194	lent = len_trim(path_to_input)
195	!IF ( path_to_input(lent:lent) /= "/" ) THEN
196	! path_to_input = TRIM(path_to_input)//"/"
197	!ENDIF
198	!lent = len_trim(path_to_output)
199	!IF ( path_to_output(lent:lent) /= "/" ) THEN
200	! path_to_output = TRIM(path_to_output)//"/"
201	!ENDIF
202	input_name = TRIM(path_to_input)//TRIM(input_name)
203	!
204	! BUILD a UNIX COMMAND based on "ls" of all of the input files
205	!
206	loslen = LEN ( command )
207	CALL all_spaces ( command , loslen )
208	WRITE ( command , FMT='("ls -1 ",A," > .foo")' ) TRIM ( input_name )
209
210	!
211	! NUMBER OF INPUTS FILES
212	!
213	! We stuck all of the matching files in the ".foo" file. Now we place the
214	! number of the those file (i.e. how many there are) in ".foo1".
215	CALL SYSTEM ( TRIM ( command ) )
216	CALL SYSTEM ( '( cat .foo \| wc -l > .foo1 )' )
217	! Read the number of files.
218	OPEN (FILE = '.foo1' , &
219	UNIT = 112 , &
220	STATUS = 'OLD' , &
221	ACCESS = 'SEQUENTIAL' , &
222	FORM = 'FORMATTED' )
223	READ ( 112 , * ) number_of_input_files
224	CLOSE ( 112 )
225	! If there are zero files, we are toast.
226	IF ( number_of_input_files .LE. 0 ) THEN
227	print*, ' Oops, we need at least ONE input file for the program to read.'
228	print*, ' Make sure you have the path, and name file(s) correct,'
229	print*, ' including wild characters if needed.'
230	STOP
231	END IF
232
233	!
234	! GET FILENAMES IN THE PROGRAM
235	!
236	! Allocate space for this many files.
237	ALLOCATE ( input_file_names(number_of_input_files) , STAT=ierr )
238	ALLOCATE ( output_file_names(number_of_input_files) , STAT=ierr )
239	! Did the allocate work OK?
240	IF ( ierr .NE. 0 ) THEN
241	print*, ' tried to allocate ', number_of_input_files, ' input files, (look at ./foo)'
242	STOP
243	END IF
244	! Initialize all of the file names to blank.
245	input_file_names = ' ' // &
246	' ' // &
247	' '
248	output_file_names = ' ' // &
249	' ' // &
250	' '
251	! Open the file that has the list of filenames.
252	OPEN (FILE = '.foo' , &
253	UNIT = 111 , &
254	STATUS = 'OLD' , &
255	ACCESS = 'SEQUENTIAL' , &
256	FORM = 'FORMATTED' )
257	! Read all of the file names and store them - define output file name.
258	LINLOG = interp_method
259	DO loop = 1 , number_of_input_files
260	READ ( 111 , FMT='(A)' ) input_file_names(loop)
261	IF ( output_name == ' ' ) THEN
262	ilen = INDEX(TRIM(input_file_names(loop)),'/',.TRUE.)
263	output_file_names(loop) = TRIM(path_to_output)//input_file_names(loop)(ilen+1:)
264	IF ( LINLOG .le. 2 ) output_file_names(loop) = TRIM(output_file_names(loop))//"_p"
265	IF ( LINLOG == 3 ) output_file_names(loop) = TRIM(output_file_names(loop))//"_z"
266	IF ( LINLOG == 4 ) output_file_names(loop) = TRIM(output_file_names(loop))//"_zabg"
267	ELSE
268	IF ( number_of_input_files == 1 ) THEN
269	output_file_names(loop) = TRIM(path_to_output)//TRIM(output_name)
270	ELSE
271	write(tmp_name,'(A,A,"_",I4.4)') TRIM(path_to_output), TRIM(output_name), loop
272	output_file_names(loop) = tmp_name
273	ENDIF
274	ENDIF
275	END DO
276	CLOSE ( 112 )
277	print*, " "
278	! We clean up our own messes.
279	CALL SYSTEM ( '/bin/rm -f .foo' )
280	CALL SYSTEM ( '/bin/rm -f .foo1' )
281
282	!
283	! LIST OF FIELD WANTED for OUTPUT
284	!
285	! Do we have a list of field that we want on output?
286	! process_these_fields = ','
287	!!!! ICI CHAMPS SORTIS PAR DEFAULT
288	!!!! ICI CHAMPS SORTIS PAR DEFAULT
289	!!!! ICI CHAMPS SORTIS PAR DEFAULT
290	process_these_fields = ',Times,XLAT,XLONG,' !! les premiere et derniere virgules sont importantes
291	!!!! ICI CHAMPS SORTIS PAR DEFAULT
292	!!!! ICI CHAMPS SORTIS PAR DEFAULT
293	!!!! ICI CHAMPS SORTIS PAR DEFAULT
294	IF ( INDEX(process,'list') /= 0) THEN
295	DO i = 1 , len(fields)
296	IF (fields(i:i) /= ' ' ) THEN
297	process_these_fields = trim(process_these_fields)//fields(i:i)
298	ENDIF
299	END DO
300	process_these_fields = trim(process_these_fields)//","
301	END IF
302
303	!
304	! HELLO to WORLD
305	!
306	write(6,*)
307	write(6,*) "##############################################################"
308	write(6,'(A,i4,A)') " RUNNING api on ", number_of_input_files, " file(s)."
309	write(6,*)
310	write(6,'(A,$)') " A) INTERPOLATION METHOD: "
311	IF ( LINLOG == 1 ) write(6,*) " Pressure - linear in p"
312	IF ( LINLOG == 2 ) write(6,*) " Pressure - linear in log p"
313	IF ( LINLOG == 3 ) write(6,*) " Height above areoid (MOLA zero datum)"
314	IF ( LINLOG == 4 ) write(6,*) " Height above surface"
315	write(6,'(A,$)') " B) VERTICAL EXTRAPOLATION: "
316	IF (extrapolate == 0) write(6,*) " BELOW GROUND and ABOVE model top will be set to missing values "
317	IF (extrapolate == 1) write(6,*) " BELOW GROUND will be extrapolated and ABOVE model top will be set to values at model top"
318	write(6,'(A,$)') " C) HORIZONTAL GRID: "
319	IF (.not. unstagger_grid) write(6,*) " Data will be output on C-grid "
320	IF (unstagger_grid) write(6,*) " Data will be output on unstaggered grid "
321
322	!
323	! GET LEVELS (pressure or altitude) TO INTERPOLATE TO
324	!
325	write(6,*)
326	IF ( LINLOG .le. 2 ) write(6,*) "INTERPOLATING TO PRESSURE LEVELS: "
327	IF ( LINLOG == 3 ) write(6,*) "INTERPOLATING TO ALTITUDE LEVELS: "
328	IF ( LINLOG == 4 ) write(6,*) "INTERPOLATING TO ABG LEVELS: "
329	num_metgrid_levels = 0
330	DO
331	IF (interp_levels(num_metgrid_levels+1) == -99999.) EXIT
332	num_metgrid_levels = num_metgrid_levels + 1
333	if (mod(num_metgrid_levels,8) /= 0 )write(6,'(f8.3,$)') interp_levels(num_metgrid_levels)
334	if (mod(num_metgrid_levels,8) == 0 )write(6,'(f8.3)') interp_levels(num_metgrid_levels)
335	IF ( LINLOG .le. 2 ) interp_levels(num_metgrid_levels) = interp_levels(num_metgrid_levels) * 100.0 !!! Pa
336	IF ( LINLOG .gt. 2 ) interp_levels(num_metgrid_levels) = interp_levels(num_metgrid_levels) * 1000.0 !!! km
337	END DO
338	write(6,*)
339	write(6,*)
340
341	!
342	! LOOP on FILES
343	!
344	DO loop = 1, number_of_input_files
345
346	IF (debug) write(6,*) "##############################################################"
347	input_file = input_file_names(loop)
348	output_file = output_file_names(loop)
349
350	IF ( .not. debug ) write(6,*) " Output will be written to: ",trim(output_file)
351
352	IF (debug) THEN
353	write(6,*) " INPUT FILE: ",trim(input_file)
354	write(6,*) " OUTPUT FILE: ",trim(output_file)
355	write(6,*) " "
356	ENDIF
357
358	!
359	! OPEN INPUT AND OUTPUT FILE
360	!
361	rcode = nf_open(input_file, 0, ncid)
362	if (rcode .ne. nf_noerr) call handle_err(rcode)
363	if (bit64) then
364	rcode = nf_create(output_file, NF_64BIT_OFFSET, mcid)
365	else
366	rcode = nf_create(output_file, 0, mcid)
367	endif
368	if (rcode .ne. nf_noerr) call handle_err(rcode)
369
370	!
371	! GET BASIC INFORMATION ABOUT THE FILE
372	! most important
373	! ndims: number of dimensions
374	! nvars: number of variables
375	! ngatts: number of global attributes
376	rcode = nf_inq(ncid, ndims, nvars, ngatts, nunlimdimid)
377	if (rcode .ne. nf_noerr) call handle_err(rcode)
378	IF (debug) THEN
379	write(6,*) ' INPUT file has = ',ndims, ' dimensions, '
380	write(6,*) ' ',nvars, ' variables, and '
381	write(6,*) ' ',ngatts,' global attributes '
382	write(6,*) " "
383	ENDIF
384	rcode = nf_get_att_int (ncid, nf_global, 'WEST-EAST_GRID_DIMENSION', iweg)
385	rcode = nf_get_att_int (ncid, nf_global, 'SOUTH-NORTH_GRID_DIMENSION', isng)
386	rcode = nf_get_att_int (ncid, nf_global, 'BOTTOM-TOP_GRID_DIMENSION', ibtg)
387
388	!
389	! ALLOCATE SOME VARIABLES
390	!
391	IF (ALLOCATED(dnamei)) deallocate(dnamei)
392	ALLOCATE (dnamei(20))
393	IF (ALLOCATED(dnamej)) deallocate(dnamej)
394	ALLOCATE (dnamej(20))
395	IF (ALLOCATED(dvali)) deallocate(dvali)
396	ALLOCATE (dvali(20))
397	IF (ALLOCATED(dvalj)) deallocate(dvalj)
398	ALLOCATE (dvalj(20))
399
400	!
401	! READ ALL DIMS FROM INPUT FILE AND CREATE SOME DIMS FOR OUTPUT FILE
402	!
403	j = 0
404	DO i = 1, ndims
405	rcode = nf_inq_dim(ncid, i, dnamei(i), dvali(i))
406	IF ( dnamei(i) == "Time" ) THEN
407	j = j + 1
408	dnamej(j) = dnamei(i)
409	dvalj(j) = dvali(i)
410	rcode = nf_def_dim(mcid, dnamej(j), NF_UNLIMITED, j)
411	times_in_file = dvali(i)
412	ENDIF
413	ENDDO
414	!!! Create the new height/pressure dims
415	j = j + 1
416	IF ( LINLOG .le. 2 ) THEN
417	dnamej(j) = 'pressure'
418	dvalj(j) = num_metgrid_levels
419	ELSE IF ( LINLOG .eq. 3 ) THEN
420	dnamej(j) = 'bottom_top' !'altitude'
421	dvalj(j) = num_metgrid_levels
422	ELSE
423	dnamej(j) = 'bottom_top' !'altitude_abg'
424	dvalj(j) = num_metgrid_levels
425	ENDIF
426	rcode = nf_def_dim(mcid, dnamej(j), dvalj(j), j)
427
428	!
429	! DEALING WITH THE GLOBAL ATTRIBUTES
430	!
431	IF (debug) THEN
432	write(6,*)
433	write(6,*) " OUTPUT FILE attributes:"
434	ENDIF
435	do i = 1, ngatts
436	rcode = nf_inq_attname(ncid, nf_global, i, cname)
437	rcode = nf_inq_atttype(ncid, nf_global, cname, itype)
438	rcode = nf_inq_attlen (ncid, nf_global, cname, ilen)
439	if ( itype .eq. 2 ) then ! characters
440	rcode = nf_get_att_text (ncid, nf_global, cname, cval)
441	if(cname(1:5) .eq. 'TITLE') then
442	IF ( LINLOG .le. 2 ) THEN
443	cval = cval(1:ilen)//" - ON PRES LEVELS"
444	ELSE
445	cval = cval(1:ilen)//" - ON z LEVELS"
446	ENDIF
447	ilen = len_trim(cval)
448	endif
449	IF (debug) &
450	write(6,'(" i = ",i2," : ",A," = ",A)') &
451	i,cname,cval(1:ilen)
452	rcode = nf_put_att_text(mcid, nf_global, cname, ilen,&
453	cval(1:ilen))
454	elseif ( itype .eq. 4 ) then ! integers
455	rcode = nf_get_att_int (ncid, nf_global, cname, ival)
456	IF ( INDEX(cname,'BOTTOM-TOP_PATCH') == 0 ) THEN
457	IF (cname .eq. 'BOTTOM-TOP_GRID_DIMENSION') ival = num_metgrid_levels
458	IF (debug) &
459	write(6,'(" i = ",i2," : ",A," = ",i7)') &
460	i,cname,ival
461	rcode = nf_put_att_int(mcid, nf_global, cname, itype,&
462	ilen, ival)
463	ENDIF
464	IF (cname .eq. 'MAP_PROJ') map_proj = ival
465	elseif ( itype .eq. 5 ) then ! real
466	rcode = nf_get_att_real (ncid, nf_global, cname, rval)
467	IF (debug) &
468	write(6,'(" i = ",i2," : ",A," = ",G18.10E2)') &
469	i,cname,rval
470	rcode = nf_put_att_real(mcid, nf_global, cname, itype,&
471	ilen, rval)
472	IF (cname .eq. 'TRUELAT1') truelat1 = rval
473	IF (cname .eq. 'TRUELAT2') truelat2 = rval
474	IF (cname .eq. 'STAND_LON') stand_lon = rval
475	end if
476	enddo
477	rcode = nf_enddef(mcid)
478
479	!
480	! WE NEED SOME BASIC FIELDS
481	!
482	! --> P_TOP
483	!
484	IF ( LINLOG .le. 2 ) THEN
485	IF ( DEBUG ) PRINT *, 'P_TOP'
486	IF (ALLOCATED(data1)) deallocate(data1)
487	allocate (data1(times_in_file,1,1,1))
488	rcode = nf_inq_varid ( ncid, "P_TOP", i )
489	rcode = nf_get_var_real ( ncid, i, data1 )
490	IF ( first ) THEN
491	IF ( extrapolate == 1 .AND. &
492	(data1(1,1,1,1)-interp_levels(num_metgrid_levels)) > 0.0 ) THEN
493	write(6,*)
494	write(6,*) "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
495	write(6,*) " WARNING: Highest requested pressure level is above PTOP."
496	write(6,'(A,F7.2,A)') " Use all pressure level data above", data1(1,1,1,1)*.01, " mb"
497	write(6,*) " with caution."
498	write(6,*) "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
499	ENDIF
500	first = .FALSE.
501	ENDIF
502	deallocate (data1)
503	ENDIF
504	!
505	! --> INTERPOLATION LEVELS
506	!
507	IF (ALLOCATED(pres_out)) deallocate(pres_out)
508	allocate (pres_out(iweg-1, isng-1, num_metgrid_levels, times_in_file))
509	do i = 1, num_metgrid_levels
510	pres_out (:,:,i,:) = interp_levels(i)
511	enddo
512	!
513	! --> LONGITUDES + everything for ROTATED WINDS
514	!
515	IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'uvmet') /= 0 ) THEN
516	IF (ALLOCATED(longi)) deallocate(longi)
517	IF (ALLOCATED(longit)) deallocate(longit)
518	allocate (longi(iweg-1, isng-1))
519	allocate (longit(iweg-1, isng-1, times_in_file))
520	IF ( DEBUG ) PRINT *, 'LONGITUDE'
521	rcode = nf_inq_varid ( ncid, "XLONG", i )
522	rcode = nf_get_var_real ( ncid, i, longit )
523	longi = longit(:,:,1)
524	deallocate(longit)
525	IF (ALLOCATED(lati)) deallocate(lati)
526	IF (ALLOCATED(latit)) deallocate(latit)
527	allocate (lati(iweg-1, isng-1))
528	allocate (latit(iweg-1, isng-1, times_in_file))
529	IF ( DEBUG ) PRINT *, 'LATITUDE'
530	rcode = nf_inq_varid ( ncid, "XLAT", i )
531	rcode = nf_get_var_real ( ncid, i, latit )
532	lati = latit(:,:,1)
533	deallocate(latit)
534	IF (ALLOCATED(u)) deallocate(u)
535	IF (oldvar) THEN
536	allocate (u (iweg, isng-1, ibtg-1, times_in_file ))
537	IF ( DEBUG ) PRINT *, 'U COMPONENT'
538	rcode = nf_inq_varid ( ncid, "U", i )
539	rcode = nf_get_var_real ( ncid, i, u )
540	ELSE
541	allocate (u (iweg-1, isng-1, ibtg-1, times_in_file ))
542	IF ( DEBUG ) PRINT *, 'UAVE COMPONENT'
543	rcode = nf_inq_varid ( ncid, "UAVE", i )
544	rcode = nf_get_var_real ( ncid, i, u )
545	ENDIF
546	IF (ALLOCATED(v)) deallocate(v)
547	IF (oldvar) THEN
548	allocate (v (iweg-1, isng, ibtg-1, times_in_file ))
549	IF ( DEBUG ) PRINT *, 'V COMPONENT'
550	rcode = nf_inq_varid ( ncid, "V", i )
551	rcode = nf_get_var_real ( ncid, i, v )
552	ELSE
553	allocate (v (iweg-1, isng-1, ibtg-1, times_in_file ))
554	IF ( DEBUG ) PRINT *, 'VAVE COMPONENT'
555	rcode = nf_inq_varid ( ncid, "VAVE", i )
556	rcode = nf_get_var_real ( ncid, i, v )
557	ENDIF
558	IF (ALLOCATED(umet)) deallocate(umet)
559	allocate (umet (iweg-1, isng-1, ibtg-1, times_in_file ))
560	IF (ALLOCATED(vmet)) deallocate(vmet)
561	allocate (vmet (iweg-1, isng-1, ibtg-1, times_in_file ))
562	IF (ALLOCATED(interm1)) deallocate(interm1)
563	allocate (interm1(iweg-1, isng-1, ibtg-1 ))
564	IF (ALLOCATED(interm2)) deallocate(interm2)
565	allocate (interm2(iweg-1, isng-1, ibtg-1 ))
566	IF ( DEBUG ) PRINT *, 'ROTATE WINDS'
567	write(6,*) " Data will be output on unstaggered grid "
568	do kk = 1, times_in_file
569	!IF ( DEBUG ) print *, kk
570	IF (oldvar) THEN
571	interm1(1:iweg-1,:,:) = ( u(1:iweg-1,:,:,kk) + u(2:iweg,:,:,kk) ) * .5
572	interm2(:,1:isng-1,:) = ( v(:,1:isng-1,:,kk) + v(:,2:isng,:,kk) ) * .5
573	ELSE
574	interm1(1:iweg-1,:,:) = u(1:iweg-1,:,:,kk)
575	interm2(:,1:isng-1,:) = v(:,1:isng-1,:,kk)
576	ENDIF
577	if (unstagger_grid) map_proj=3 !!! AS: unstagger_grid=T makes the program to put unstaggered U and V in uvmet
578	CALL calc_uvmet( interm1, interm2, umet(:,:,:,kk), vmet(:,:,:,kk), &
579	truelat1, truelat2, stand_lon, map_proj, longi, lati, iweg-1, isng-1, ibtg-1 )
580	enddo
581	deallocate(lati)
582	deallocate(longi)
583	deallocate(u)
584	deallocate(v)
585	deallocate(interm1)
586	deallocate(interm2)
587	ENDIF
588
589	!
590	! --> GEOPOTENTIAL HEIGHT (in mode 1-2, no need if geopotential is not requested)
591	!
592	IF ( LINLOG .gt. 2 .OR. INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'GHT') /= 0 ) THEN
593	IF ( DEBUG ) PRINT *, 'GEOPOTENTIAL'
594	!IF (ALLOCATED(ght)) deallocate(ght)
595	!allocate (ght(iweg-1, isng-1, ibtg-1, times_in_file))
596	IF (ALLOCATED(phb)) deallocate(phb)
597	allocate (phb(iweg-1, isng-1, ibtg, times_in_file ))
598	! IF (ALLOCATED(data1)) deallocate(data1)
599	! allocate (data1(iweg-1, isng-1, ibtg, times_in_file ))
600	! rcode = nf_inq_varid ( ncid, "PH", i )
601	! rcode = nf_get_var_real ( ncid, i, data1 )
602	! rcode = nf_inq_varid ( ncid, "PHB", i )
603	! rcode = nf_get_var_real ( ncid, i, phb )
604	! data1 = (data1 + phb)
605	! ght(:,:,1:ibtg-1,:) = ( data1(:,:,1:ibtg-1,:) + data1(:,:,2:ibtg,:) )*.5
606	! deallocate (data1)
607	rcode = nf_inq_varid ( ncid, "PHTOT", i )
608	rcode = nf_get_var_real ( ncid, i, phb )
609	!ght(:,:,1:ibtg-1,:) = phb(:,:,1:ibtg-1,:) !! costly in memory !!
610	!deallocate (phb)
611	ENDIF
612	!
613	! --> PRESSURE (in mode 3-4, no need if regular temperature is not requested)
614	!
615	IF ( LINLOG .le. 2 .OR. INDEX(process,'all') /= 0 &
616	.OR. INDEX(process_these_fields,'tk') /= 0 &
617	.OR. INDEX(process_these_fields,'tpot') /= 0 ) THEN
618	IF ( DEBUG ) PRINT *, 'PRESSURE'
619	IF (ALLOCATED(pres_field)) deallocate(pres_field)
620	allocate (pres_field(iweg-1, isng-1, ibtg-1, times_in_file ))
621	! IF (ALLOCATED(data1)) deallocate(data1)
622	! allocate (data1(iweg-1, isng-1, ibtg-1, times_in_file ))
623	! rcode = nf_inq_varid ( ncid, "P", i )
624	! rcode = nf_get_var_real ( ncid, i, pres_field )
625	! rcode = nf_inq_varid ( ncid, "PB", i )
626	! rcode = nf_get_var_real ( ncid, i, data1 )
627	! pres_field = pres_field + data1
628	rcode = nf_inq_varid ( ncid, "PTOT", i )
629	rcode = nf_get_var_real ( ncid, i, pres_field )
630	! deallocate (data1)
631	!
632	! --> REGULAR and POTENTIAL TEMPERATURE
633	!
634	IF ( DEBUG ) PRINT *, 'TEMP'
635	IF (ALLOCATED(tpot)) deallocate(tpot)
636	allocate (tpot(iweg-1, isng-1, ibtg-1, times_in_file ))
637	IF (oldvar) THEN
638	rcode = nf_inq_varid ( ncid, "T", i )
639	rcode = nf_get_var_real ( ncid, i, tpot )
640	tpot = tpot + tpot0
641	ELSE
642	rcode = nf_inq_varid ( ncid, "TAVE", i )
643	rcode = nf_get_var_real ( ncid, i, tpot )
644	ENDIF
645	IF ( LINLOG .le. 2 .OR. INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'tk') /= 0 ) THEN
646	IF (ALLOCATED(tk)) deallocate(tk)
647	allocate (tk(iweg-1, isng-1, ibtg-1, times_in_file ))
648	tk = tpot * ( pres_field / p0 )**RCP
649	ENDIF
650	IF (INDEX(process_these_fields,'tpot') == 0 .AND. INDEX(process,'all') == 0) deallocate(tpot)
651	ENDIF
652	!
653	! --> SURFACE PRESSURE
654	!
655	IF ( LINLOG .le. 2 ) THEN
656	IF ( DEBUG ) PRINT *, 'PSFC'
657	IF (ALLOCATED(psfc)) deallocate(psfc)
658	allocate (psfc(iweg-1, isng-1, times_in_file ))
659	IF (ALLOCATED(data1)) deallocate(data1)
660	allocate (data1(iweg-1, isng-1, 1, times_in_file ))
661	rcode = nf_inq_varid ( ncid, "PSFC", i )
662	rcode = nf_get_var_real ( ncid, i, data1 )
663	psfc(:,:,:) = data1(:,:,1,:)
664	deallocate (data1)
665	ENDIF
666	!
667	! --> TOPOGRAPHY
668	!
669	IF ( LINLOG .ne. 3 ) THEN
670	IF ( DEBUG ) PRINT *, 'TOPO'
671	IF (ALLOCATED(ter)) deallocate(ter)
672	allocate (ter(iweg-1, isng-1))
673	IF (ALLOCATED(data1)) deallocate(data1)
674	allocate (data1(iweg-1, isng-1, 1, times_in_file ))
675	rcode = nf_inq_varid ( ncid, "HGT", i )
676	rcode = nf_get_var_real ( ncid, i, data1 )
677	ter(:,:) = data1(:,:,1,1)
678	IF ( LINLOG .ne. 4 ) deallocate (data1)
679	ENDIF
680
681	! IF (ALLOCATED(qv)) deallocate(qv)
682	! allocate (qv(iweg-1, isng-1, ibtg-1, times_in_file ))
683	! rcode = nf_inq_varid ( ncid, "QVAPOR", i )
684	! rcode = nf_get_var_real ( ncid, i, qv )
685	!
686	! NOT SUITABLE for MARS
687	!
688	! IF (ALLOCATED(rh)) deallocate(rh)
689	! allocate (rh(iweg-1, isng-1, ibtg-1, times_in_file ))
690	! IF (ALLOCATED(data1)) deallocate(data1)
691	! IF (ALLOCATED(data2)) deallocate(data2)
692	! allocate (data1(iweg-1, isng-1, ibtg-1, times_in_file ))
693	! allocate (data2(iweg-1, isng-1, ibtg-1, times_in_file ))
694	! data1 = 10.0.6112exp(17.67*(tk-273.16)/(TK-29.65))
695	! data2 = 0.622data1/(0.01 pres_field - (1.-0.622)*data1)
696	! rh = 100.*AMAX1(AMIN1(qv/data2,1.0),0.0)
697	! deallocate (data1)
698	! deallocate (data2)
699	!
700	! --> USEFUL for INTERPOLATION
701	!
702	!!! at that point pres_field is not used for calculation anymore
703	!!! it is used to set the initial coordinate for interpolation
704	IF ( LINLOG == 3) THEN
705	IF (ALLOCATED(pres_field)) deallocate(pres_field)
706	allocate ( pres_field(iweg-1, isng-1, ibtg-1, times_in_file) )
707	pres_field = phb(:,:,1:ibtg-1,:) / grav
708	ENDIF
709	IF ( LINLOG == 4) THEN
710	IF (ALLOCATED(pres_field)) deallocate(pres_field)
711	allocate ( pres_field(iweg-1, isng-1, ibtg-1, times_in_file))
712	DO kk = 1, ibtg-1
713	pres_field(:,:,kk,:) = - data1(:,:,1,:) + phb(:,:,kk,:) / grav
714	ENDDO
715	deallocate (data1)
716	!PRINT *, pres_field(10,10,:,1)
717	ENDIF
718
719
720	IF (ALLOCATED(pres_stagU)) deallocate(pres_stagU)
721	IF (ALLOCATED(pres_stagV)) deallocate(pres_stagV)
722	IF ( DEBUG ) PRINT *, 'STAGGERED COORDINATES'
723	allocate (pres_stagU(iweg, isng-1, ibtg-1, times_in_file ))
724	allocate (pres_stagV(iweg-1, isng, ibtg-1, times_in_file ))
725	pres_stagU(1,:,:,:) = pres_field(1,:,:,:)
726	pres_stagU(iweg,:,:,:) = pres_field(iweg-1,:,:,:)
727	pres_stagU(2:iweg-1,:,:,:) = (pres_field(1:iweg-2,:,:,:) + pres_field(2:iweg-1,:,:,:))*.5
728	pres_stagV(:,1,:,:) = pres_field(:,1,:,:)
729	pres_stagV(:,isng,:,:) = pres_field(:,isng-1,:,:)
730	pres_stagV(:,2:isng-1,:,:) = (pres_field(:,1:isng-2,:,:) + pres_field(:,2:isng-1,:,:))*.5
731
732	!-----------------------------
733	!-----------------------------
734	! TRAIN FILE
735	!-----------------------------
736	!-----------------------------
737	IF (debug) THEN
738	write(6,*)
739	write(6,*)
740	write(6,*) "FILE variables:"
741	ENDIF
742	jvar = 0
743	loop_variables : DO ivar = 1, nvars
744
745	rcode = nf_inq_var(ncid, ivar, cval, itype, idm, ishape, natt)
746
747	!!! Do we want this variable ?
748	! IF ( trim(cval) == 'P' .OR. trim(cval) == 'PB' ) CYCLE loop_variables
749	! IF ( trim(cval) == 'PH' .OR. trim(cval) == 'PHB' ) CYCLE loop_variables
750	!IF ( trim(cval) == 'PTOT' ) CYCLE loop_variables ! PTOT: no
751	IF ( trim(cval) == 'PHTOT' ) CYCLE loop_variables ! PHTOT: no
752	IF ( trim(cval) == 'T' ) CYCLE loop_variables ! T: no (tk and tpot calculated above)
753	IF ( unstagger_grid .AND. (INDEX(cval,'_U') /= 0) ) CYCLE loop_variables !!! no sense in keeping these
754	IF ( unstagger_grid .AND. (INDEX(cval,'_V') /= 0) ) CYCLE loop_variables !!! no sense in keeping these
755	IF ( INDEX(process,'all') == 0 ) THEN
756	!!! Only want some variables - see which
757	dummy = ","//trim(cval)//","
758	is_there = INDEX(process_these_fields,trim(dummy))
759	IF ( is_there == 0 ) THEN
760	IF ( debug ) print*,"NOTE: ", trim(cval), " - Not requested"
761	CYCLE loop_variables !!! don't want this one
762	ENDIF
763	ENDIF
764
765	IF ( idm >= 4 .AND. itype == 4 ) THEN
766	print*,"NOTE: We cannot deal with 3D integers - maybe later"
767	CYCLE loop_variables
768	ENDIF
769	IF ( itype == 6 ) THEN
770	print*,"NOTE: We cannot deal with double precision data - maybe later"
771	CYCLE loop_variables
772	ENDIF
773	IF ( itype == 2 .OR. itype == 4 .OR. itype == 5 ) THEN
774	!!! OK I know what to do this this
775	ELSE
776	print*,"NOTE: Do not understand this data type ", itype, " skip field."
777	CYCLE loop_variables
778	ENDIF
779
780	!IF ( trim(cval) == 'U' ) cval = 'UU'
781	!IF ( trim(cval) == 'V' ) cval = 'VV'
782	!IF ( trim(cval) == 'T' ) cval = 'TT'
783
784	!!! OK - we want this - lets continue
785	jvar = jvar + 1
786	jshape = 0
787	interpolate = .FALSE.
788	fix_meta_stag = .FALSE.
789	rcode = nf_redef(mcid)
790	DO ii = 1, idm
791	test_dim_name = dnamei(ishape(ii))
792	IF ( test_dim_name == 'bottom_top' .OR. test_dim_name == 'bottom_top_stag' ) THEN
793	IF ( test_dim_name == 'bottom_top_stag' ) fix_meta_stag = .TRUE.
794	IF (LINLOG .le. 2) test_dim_name = 'pressure'
795	IF (LINLOG .eq. 3) test_dim_name = 'bottom_top' !'altitude'
796	IF (LINLOG .eq. 4) test_dim_name = 'bottom_top' !'altitude_abg'
797	interpolate = .TRUE.
798	ENDIF
799	IF ( unstagger_grid .AND. test_dim_name == 'west_east_stag' ) THEN
800	test_dim_name = 'west_east'
801	fix_meta_stag = .TRUE.
802	ENDIF
803	IF ( unstagger_grid .AND. test_dim_name == 'south_north_stag' ) THEN
804	test_dim_name = 'south_north'
805	fix_meta_stag = .TRUE.
806	ENDIF
807	DO jj = 1,j
808	IF ( test_dim_name == dnamej(jj) ) THEN
809	jshape(ii) = jj
810	ENDIF
811	ENDDO
812
813	IF ( jshape(ii) == 0 ) THEN
814	j = j + 1
815	jshape(ii) = j
816	dnamej(j) = dnamei(ishape(ii))
817	dvalj(j) = dvali(ishape(ii))
818	rcode = nf_def_dim(mcid, dnamej(j), dvalj(j), j)
819	ENDIF
820	ENDDO
821	rcode = nf_def_var(mcid, cval, itype, idm, jshape, jvar)
822	!print *, 'mcid', mcid
823	!print *, 'cval', cval
824	!print *, 'itype', itype
825	!print *, 'idm', idm
826	!print *, 'jshape', jshape
827	!print *, 'jvar', jvar
828
829	DO na = 1, natt
830	rcode = nf_inq_attname(ncid, ivar, na, cname)
831	IF ( fix_meta_stag .AND. trim(cname) == 'stagger' ) THEN
832	att_text = "-"
833	ilen = len_trim(att_text)
834	rcode = nf_put_att_text(mcid, jvar, cname, ilen, att_text(1:ilen) )
835	ELSEIF ( fix_meta_stag .AND. trim(cname) == 'coordinates' ) THEN
836	att_text = "XLONG XLAT"
837	ilen = len_trim(att_text)
838	rcode = nf_put_att_text(mcid, jvar, cname, ilen, att_text(1:ilen) )
839	ELSE
840	rcode = nf_copy_att(ncid, ivar, cname, mcid, jvar)
841	ENDIF
842	ENDDO
843
844	IF ( extrapolate == 0 ) THEN
845	rcode = nf_put_att_real(mcid, jvar, 'missing_value', NF_FLOAT, 1, MISSING )
846	ENDIF
847
848	rcode = nf_enddef(mcid)
849
850	!
851	! GET THE DIMS FOR INPUT AND OUTPUT FROM THE SHAPE
852	!
853	dims_in = 1
854	dims_out = 1
855	DO ii = 1,idm
856	dims_in(ii) = dvali(ishape(ii))
857	dims_out(ii) = dvalj(jshape(ii))
858	ENDDO
859	IF (debug) write(6,*) 'VAR: ',trim(cval)
860	IF (debug) THEN
861	write(6,*) ' DIMS IN: ',dims_in
862	write(6,*) ' DIMS OUT: ',dims_out
863	ENDIF
864
865	!
866	! ALLOCATE THE INPUT AND OUTPUT ARRAYS
867	! READ THE DATA FROM INPUT FILE
868	!
869	IF (itype == 2) THEN ! character
870	allocate (text(dims_in(1), dims_in(2), dims_in(3), dims_in(4)))
871	rcode = nf_get_var_text(ncid, ivar, text)
872	rcode = nf_put_vara_text (mcid, jvar, start_dims, dims_in, text)
873	IF (debug) write(6,*) ' SAMPLE VALUE = ',text(:,1,1,1)
874	deallocate (text)
875
876	ELSEIF (itype == 4) THEN ! integer
877	allocate (idata1(dims_in(1), dims_in(2), dims_in(3), dims_in(4)))
878	rcode = nf_get_var_int(ncid, ivar, idata1)
879	rcode = nf_put_vara_int (mcid, jvar, start_dims, dims_in, idata1)
880	IF (debug) write(6,*) ' SAMPLE VALUE = ',idata1(dims_in(1)/2,dims_in(2)/2,1,1)
881	deallocate (idata1)
882
883	ELSEIF (itype == 5) THEN ! real
884	allocate (data1(dims_in(1), dims_in(2), dims_in(3), dims_in(4)))
885	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
886	rcode = nf_get_var_real(ncid, ivar, data1)
887
888	IF (idm >= 4 .AND. interpolate) THEN
889	IF (debug) write(6,*) ' THIS IS A FIELD WE NEED TO INTERPOLATE'
890
891	IF ( dims_in(1) == iweg .AND. .not. unstagger_grid ) THEN
892	CALL interp (data2, data1, pres_stagU, interp_levels, psfc, ter, tk, qv, &
893	iweg, isng-1, ibtg-1, dims_in(4), &
894	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
895	ELSEIF ( dims_in(2) == isng .AND. .not. unstagger_grid ) THEN
896	CALL interp (data2, data1, pres_stagV, interp_levels, psfc, ter, tk, qv, &
897	iweg-1, isng, ibtg-1, dims_in(4), &
898	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
899	ELSEIF ( dims_in(1) == iweg .AND. unstagger_grid ) THEN
900	allocate (data3(iweg-1, isng-1, ibtg-1, dims_in(4)))
901	data3(1:iweg-1,:,:,:) = (data1(1:iweg-1,:,:,:) + data1(2:iweg,:,:,:)) * .5
902	CALL interp (data2, data3, pres_field, interp_levels, psfc, ter, tk, qv, &
903	iweg-1, isng-1, ibtg-1, dims_in(4), &
904	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
905	deallocate(data3)
906	ELSEIF ( dims_in(2) == isng .AND. unstagger_grid ) THEN
907	allocate (data3(iweg-1, isng-1, ibtg-1, dims_in(4)))
908	data3(:,1:isng-1,:,:) = (data1(:,1:isng-1,:,:) + data1(:,2:isng,:,:)) * .5
909	CALL interp (data2, data3, pres_field, interp_levels, psfc, ter, tk, qv, &
910	iweg-1, isng-1, ibtg-1, dims_in(4), &
911	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
912	deallocate(data3)
913	ELSEIF ( dims_in(3) == ibtg ) THEN
914	allocate (data3(iweg-1, isng-1, ibtg-1, dims_in(4)))
915	IF (debug) PRINT *, 'VERTICAL UNSTAGGERING'
916	data3(:,:,1:ibtg-1,:) = (data1(:,:,1:ibtg-1,:) + data1(:,:,2:ibtg,:)) * .5
917	CALL interp (data2, data3, pres_field, interp_levels, psfc, ter, tk, qv, &
918	dims_in(1), dims_in(2), ibtg-1, dims_in(4), &
919	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
920	deallocate(data3)
921	ELSE
922	CALL interp (data2, data1, pres_field, interp_levels, psfc, ter, tk, qv, &
923	dims_in(1), dims_in(2), dims_in(3), dims_in(4), &
924	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
925	END IF
926
927	IF (debug) write(6,*) ' SAMPLE VALUE IN = ',data1(dims_in(1)/2,dims_in(2)/2,1,1)
928	IF (debug) write(6,*) ' SAMPLE VALUE OUT = ',data2(dims_out(1)/2,dims_out(2)/2,1,1)
929
930	ELSEIF (idm == 3 .AND. unstagger_grid ) THEN
931	IF ( dims_in(1) == iweg ) THEN
932	data2(1:iweg-1,:,:,:) = (data1(1:iweg-1,:,:,:) + data1(2:iweg,:,:,:)) * .5
933	ELSEIF ( dims_in(2) == isng ) THEN
934	data2(:,1:isng-1,:,:) = (data1(:,1:isng-1,:,:) + data1(:,2:isng,:,:)) * .5
935	ELSE
936	data2 = data1
937	ENDIF
938	IF (debug) write(6,*) ' SAMPLE VALUE = ',data1(dims_in(1)/2,dims_in(2)/2,1,1)
939
940	ELSE
941	data2 = data1
942	IF (debug) write(6,*) ' SAMPLE VALUE = ',data1(dims_in(1)/2,dims_in(2)/2,1,1)
943
944	ENDIF
945
946	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
947	!print *, 'mcid ', mcid
948	!print *, 'jvar ', jvar
949	!print *, 'start_dims ', start_dims
950	!print *, 'dims_out ', dims_out
951
952	deallocate (data1)
953	deallocate (data2)
954
955	ENDIF
956
957	ENDDO loop_variables
958
959	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
960	!!!!!! ADDITIONAL VARIABLES !!!!!!
961	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
962	IF ( debug ) print*," "
963	IF ( debug ) print*,"Calculating some diagnostics"
964	interpolate = .TRUE.
965
966	!!! NB: what follows is not useful because we'd like diagnostics for each history timestep
967	jshape = 0
968	DO ii = 1, 4
969	IF ( ii == 1 ) test_dim_name = 'west_east'
970	IF ( ii == 2 ) test_dim_name = 'south_north'
971	IF (( ii == 3 ) .and. (LINLOG .le. 2)) test_dim_name = 'pressure'
972	IF (( ii == 3 ) .and. (LINLOG .eq. 3)) test_dim_name = 'bottom_top' !'altitude'
973	IF (( ii == 3 ) .and. (LINLOG .eq. 4)) test_dim_name = 'bottom_top' !'altitude_abg'
974	IF ( ii == 4 ) test_dim_name = 'Time'
975	DO jj = 1,j
976	IF ( test_dim_name == dnamej(jj) ) THEN
977	jshape(ii) = jj
978	ENDIF
979	ENDDO
980	IF ( jshape(ii) == 0 ) THEN
981	j = j + 1
982	jshape(ii) = j
983	dnamej(j) = dnamei(ishape(ii))
984	dvalj(j) = dvali(ishape(ii))
985	rcode = nf_def_dim(mcid, dnamej(j), dvalj(j), j)
986	ENDIF
987	ENDDO
988	dims_in = 1
989	dims_out = 1
990	DO ii = 1,4
991	dims_out(ii) = dvalj(jshape(ii))
992	!print *, dims_out(ii)
993	ENDDO
994	!!! NB: what follows is useful because we'd like diagnostics for each history timestep
995	dims_in = dims_out
996
997	!IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'PRES') /= 0 ) THEN
998	! !!! PRES
999	! jvar = jvar + 1
1000	! CALL def_var (mcid, jvar, "PRES", 5, 4, jshape, "XZY", "Pressure ", "Pa", "-", "XLONG XLAT")
1001	! IF (debug) THEN
1002	! write(6,*) 'VAR: PRES'
1003	! write(6,*) ' DIMS OUT: ',dims_out
1004	! ENDIF
1005	! rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, pres_out)
1006	! IF (debug) write(6,*) ' SAMPLE VALUE OUT = ',pres_out(dims_out(1)/2,dims_out(2)/2,1,1)
1007	!ENDIF
1008
1009	!
1010	! EN FAIT IL FAUDRAIT QUE LE jshape SOIT BIEN DEFINI POUR CE QUI SUIT
1011	!
1012
1013	!
1014	! OUTPUT REGULAR TEMPERATURE
1015	!
1016	IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'tk') /= 0 ) THEN
1017	jvar = jvar + 1
1018	CALL def_var (mcid, jvar, "tk ", 5, 4, jshape, "XZY", "Temperature ", "K ", "-", "XLONG XLAT", MISSING)
1019	IF (debug) THEN
1020	write(6,*) 'VAR: tk'
1021	write(6,*) ' DIMS OUT: ',dims_out
1022	ENDIF
1023	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
1024	CALL interp (data2, tk, pres_field, interp_levels, psfc, ter, tk, qv, &
1025	iweg-1, isng-1, ibtg-1, dims_in(4), &
1026	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
1027	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
1028	IF (debug) write(6,*) ' SAMPLE VALUE OUT = ',data2(dims_out(1)/2,dims_out(2)/2,1,1)
1029	deallocate(data2)
1030	ENDIF
1031
1032	!
1033	! OUTPUT ROTATED WINDS
1034	!
1035	IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'uvmet') /= 0 ) THEN
1036	jvar = jvar + 1
1037	CALL def_var (mcid, jvar, "Um ", 5, 4, jshape, "XZY", "U rotated wind ", "K ", "-", "XLONG XLAT", MISSING)
1038	IF (debug) THEN
1039	write(6,*) 'VAR: u (rotated)'
1040	write(6,*) ' DIMS OUT: ',dims_out
1041	ENDIF
1042	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
1043	CALL interp (data2, umet, pres_field, interp_levels, psfc, ter, tk, qv, &
1044	iweg-1, isng-1, ibtg-1, dims_in(4), &
1045	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
1046	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
1047	IF (debug) write(6,*) ' SAMPLE VALUE OUT=',data2(dims_out(1)/2,dims_out(2)/2,1,1)
1048	deallocate(data2)
1049
1050	jvar = jvar + 1
1051	CALL def_var (mcid, jvar, "Vm ", 5, 4, jshape, "XZY", "V rotated wind ", "K ", "-", "XLONG XLAT", MISSING)
1052	IF (debug) THEN
1053	write(6,*) 'VAR: v (rotated)'
1054	write(6,*) ' DIMS OUT: ',dims_out
1055	ENDIF
1056	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
1057	CALL interp (data2, vmet, pres_field, interp_levels, psfc, ter, tk, qv, &
1058	iweg-1, isng-1, ibtg-1, dims_in(4), &
1059	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
1060	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
1061	IF (debug) write(6,*) ' SAMPLE VALUE OUT=',data2(dims_out(1)/2,dims_out(2)/2,1,1)
1062	deallocate(data2)
1063	ENDIF
1064
1065	!
1066	! OUTPUT POTENTIAL TEMPERATURE
1067	!
1068	IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'tpot') /= 0 ) THEN
1069	jvar = jvar + 1
1070	CALL def_var (mcid, jvar, "tpot", 5, 4, jshape, "XZY", "Potential Temperature ", "K ", "-", "XLONG XLAT", MISSING)
1071	IF (debug) THEN
1072	write(6,*) 'VAR: tpot'
1073	write(6,*) ' DIMS OUT: ',dims_out
1074	ENDIF
1075	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
1076	CALL interp (data2, tpot, pres_field, interp_levels, psfc, ter, tk, qv, &
1077	iweg-1, isng-1, ibtg-1, dims_in(4), &
1078	num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
1079	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
1080	IF (debug) write(6,*) ' SAMPLE VALUE OUT =',data2(dims_out(1)/2,dims_out(2)/2,1,1)
1081	deallocate(data2)
1082	ENDIF
1083
1084	!
1085	! OUTPUT GEOPOTENTIAL HEIGHT
1086	!
1087	IF (LINLOG .le. 2) THEN
1088	IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'GHT') /= 0 ) THEN
1089	jvar = jvar + 1
1090	CALL def_var (mcid, jvar, "GHT ", 5, 4, jshape, "XZY", "Geopotential Height", "m ", "-", "XLONG XLAT", MISSING)
1091	IF (debug) THEN
1092	write(6,*) 'VAR: GHT'
1093	write(6,*) ' DIMS OUT: ',dims_out
1094	ENDIF
1095	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
1096	CALL interp (data2, phb(:,:,1:ibtg-1,:), pres_field, interp_levels, psfc, ter, tk, qv, &
1097	iweg-1, isng-1, ibtg-1, dims_in(4), &
1098	num_metgrid_levels, LINLOG, extrapolate, .TRUE., MISSING)
1099	data2 = data2/grav
1100	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
1101	IF (debug) write(6,*) ' SAMPLE VALUE OUT = ',data2(dims_out(1)/2,dims_out(2)/2,1,1)
1102	deallocate(data2)
1103	ENDIF
1104	ELSE
1105	PRINT *, 'Geopotential height is actually vertical coordinate'
1106	ENDIF
1107
1108	! IF ( INDEX(process,'all') /= 0 .OR. INDEX(process_these_fields,'RH') /= 0 ) THEN
1109	! !!! RH
1110	! jvar = jvar + 1
1111	! CALL def_var (mcid, jvar, "RH ", 5, 4, jshape, "XZY", "Relative Humidity ", "% ", "-", "XLONG XLAT")
1112	! IF (debug) THEN
1113	! write(6,*) 'VAR: RH'
1114	! write(6,*) ' DIMS OUT: ',dims_out
1115	! ENDIF
1116	! allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
1117	! CALL interp (data2, rh, pres_field, interp_levels, psfc, ter, tk, qv, &
1118	! iweg-1, isng-1, ibtg-1, dims_in(4), &
1119	! num_metgrid_levels, LINLOG, extrapolate, .FALSE., MISSING)
1120	! WHERE ( rh < 0.0 )
1121	! rh = 0.0
1122	! ENDWHERE
1123	! WHERE ( rh > 100.0 )
1124	! rh = 100.0
1125	! ENDWHERE
1126	! rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
1127	! IF (debug) write(6,*) ' SAMPLE VALUE OUT = ',data2(dims_out(1)/2,dims_out(2)/2,1,1)
1128	! deallocate(data2)
1129	! ENDIF
1130
1131
1132	!=====================================================================================
1133	!=====================================================================================
1134	!
1135	! VERTICAL COORDINATES
1136	!
1137	jvar = jvar + 1
1138	jshape(:)=0
1139	jshape(1)=2
1140	CALL def_var (mcid, jvar, "vert", 5, 1, jshape, " Z ", "Vert. coord. ","m ", "-", "XLONG XLAT", MISSING)
1141	start_dims(1)=1
1142	start_dims(2)=1
1143	start_dims(3)=1
1144	start_dims(4)=1
1145	dims_out(1)=dims_out(3)
1146	dims_out(2)=1
1147	dims_out(3)=1
1148	dims_out(4)=1
1149	allocate (data2(dims_out(1),dims_out(2),dims_out(3),dims_out(4)))
1150	do kk=1,dims_out(1)
1151	data2(kk,1,1,1) = interp_levels(kk)
1152	enddo
1153	rcode = nf_put_vara_real (mcid, jvar, start_dims, dims_out, data2)
1154	deallocate(data2)
1155	!=====================================================================================
1156	!=====================================================================================
1157
1158
1159
1160	!
1161	! CLOSE FILES ON EXIT
1162	!
1163	rcode = nf_close(ncid)
1164	rcode = nf_close(mcid)
1165	write(6,*)
1166	ENDDO
1167
1168	!
1169	! WELL...
1170	!
1171	write(6,*) "I used Cp, R : ",Cp,Rd
1172	write(6,*)
1173	write(6,*) "##########################################"
1174	write(6,*) " END of API PROGRAM - SUCCESS so far"
1175	write(6,*) "##########################################"
1176
1177	END SUBROUTINE
1178	! END PROGRAM api
1179	!---------------------------------------------------------------------
1180	!---------------------------------------------------------------------
1181	!---------------------------------------------------------------------
1182	SUBROUTINE handle_err(rcode)
1183	INTEGER rcode
1184	write(6,*) 'Error number ',rcode
1185	stop
1186	END SUBROUTINE
1187	!---------------------------------------------------------------------
1188	SUBROUTINE interp (data_out, data_in, pres_field, interp_levels, psfc, ter, tk, qv, ix, iy, iz, it, &
1189	num_metgrid_levels, LINLOG, extrapolate, GEOPT, MISSING)
1190
1191	INTEGER :: ix, iy, iz, it
1192	INTEGER :: num_metgrid_levels, LINLOG
1193	REAL, DIMENSION(ix, iy, num_metgrid_levels, it) :: data_out
1194	REAL, DIMENSION(ix, iy, iz, it) :: data_in, pres_field, tk, qv
1195	REAL, DIMENSION(ix, iy, it) :: psfc
1196	REAL, DIMENSION(ix, iy) :: ter
1197	REAL, DIMENSION(num_metgrid_levels) :: interp_levels
1198
1199	INTEGER :: i, j, itt, k, kk, kin
1200	REAL, DIMENSION(num_metgrid_levels) :: data_out1D
1201	REAL, DIMENSION(iz) :: data_in1D, pres_field1D
1202	INTEGER :: extrapolate
1203	REAL :: MISSING
1204	REAL, DIMENSION(ix, iy, num_metgrid_levels, it) :: N
1205	REAL :: sumA, sumN, AVE_geopt
1206	LOGICAL :: GEOPT
1207
1208	N = 1.0
1209
1210	do itt = 1, it
1211	!PRINT *, 'TIME... ', itt, ' in ', iz, ' out ', num_metgrid_levels
1212	do j = 1, iy
1213	do i = 1, ix
1214	data_in1D(:) = data_in(i,j,:,itt)
1215	pres_field1D(:) = pres_field(i,j,:,itt)
1216	IF ( LINLOG .le. 2 ) THEN
1217	CALL int1D (data_out1D, data_in1D, pres_field1D, interp_levels, iz, num_metgrid_levels, LINLOG, MISSING)
1218	ELSE
1219	CALL interp_1d (data_in1D, pres_field1D, iz, data_out1D, interp_levels, num_metgrid_levels, 'z', MISSING)
1220	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1221	!interp_1d (data_in1D, pres_field1D, num_metgrid_levels, data_out1D, interp_levels, iz, 'z')
1222	!CALL interp_1d( data_in_1d, z_data_1d, bottom_top_dim, data_out_1d, z_levs, number_of_zlevs, vertical_type)
1223	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1224	ENDIF
1225	data_out(i,j,:,itt) = data_out1D(:)
1226	end do
1227	end do
1228	end do
1229	!PRINT *, 'ok'
1230
1231	! Fill in missing values
1232	IF ( extrapolate == 0 ) RETURN !! no extrapolation - we are out of here
1233
1234
1235
1236	IF ( LINLOG .ge. 0 ) RETURN !! TEMPORARY: no extrapolation
1237	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1238	!! CAUTION BELOW: METHOD NOT ADAPTED TO MARS
1239	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1240
1241	!!! MARS MARS MARS
1242	!expon=287.04*.0065/9.81
1243	expon=192.*.0045/3.72
1244
1245	! First find where about 400 hPa is located
1246	kk = 0
1247	find_kk : do k = 1, num_metgrid_levels
1248	kk = k
1249	if ( interp_levels(k) <= 40000. ) exit find_kk
1250	end do find_kk
1251
1252
1253	IF ( GEOPT ) THEN !! geopt is treated different below ground
1254
1255	do itt = 1, it
1256	do k = 1, kk
1257	do j = 1, iy
1258	do i = 1, ix
1259	IF ( data_out(i,j,k,itt) == MISSING .AND. interp_levels(k) < psfc(i,j,itt) ) THEN
1260
1261	! We are below the first model level, but above the ground
1262
1263	!!! MARS MARS
1264	data_out(i,j,k,itt) = ((interp_levels(k) - pres_field(i,j,1,itt))ter(i,j)3.72 + &
1265	(psfc(i,j,itt) - interp_levels(k))*data_in(i,j,1,itt) ) / &
1266	(psfc(i,j,itt) - pres_field(i,j,1,itt))
1267
1268	ELSEIF ( data_out(i,j,k,itt) == MISSING ) THEN
1269
1270	! We are below both the ground and the lowest data level.
1271
1272	! First, find the model level that is closest to a "target" pressure
1273	! level, where the "target" pressure is delta-p less that the local
1274	! value of a horizontally smoothed surface pressure field. We use
1275	! delta-p = 150 hPa here. A standard lapse rate temperature profile
1276	! passing through the temperature at this model level will be used
1277	! to define the temperature profile below ground. This is similar
1278	! to the Benjamin and Miller (1990) method, except that for
1279	! simplicity, they used 700 hPa everywhere for the "target" pressure.
1280	! Code similar to what is implemented in RIP4
1281
1282
1283	!!! MARS MARS MARS
1284	ptarget = (psfc(i,j,itt)*.01) - 1.50 !150.
1285	dpmin=1.e2 !1.e4
1286	kupper = 0
1287	loop_kIN : do kin=iz,1,-1
1288	kupper = kin
1289	dp=abs( (pres_field(i,j,kin,itt)*.01) - ptarget )
1290	if (dp.gt.dpmin) exit loop_kIN
1291	dpmin=min(dpmin,dp)
1292	enddo loop_kIN
1293
1294	pbot=max(pres_field(i,j,1,itt),psfc(i,j,itt))
1295	zbot=min(data_in(i,j,1,itt)/3.72,ter(i,j)) !!MARS MARS
1296
1297	tbotextrap=tk(i,j,kupper,itt)(pbot/pres_field(i,j,kupper,itt))*expon
1298	tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
1299
1300	data_out(i,j,k,itt) = (zbot+tvbotextrap/.0045(1.-(interp_levels(k)/pbot)expon))3.72 !!MARS MARS
1301
1302	ENDIF
1303	enddo
1304	enddo
1305	enddo
1306	enddo
1307
1308
1309	!!! Code for filling missing data with an average - we don't want to do this
1310	!!do itt = 1, it
1311	!!loop_levels : do k = 1, num_metgrid_levels
1312	!!sumA = SUM(data_out(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
1313	!!sumN = SUM(N(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
1314	!!IF ( sumN == 0. ) CYCLE loop_levels
1315	!!AVE_geopt = sumA/sumN
1316	!!WHERE ( data_out(:,:,k,itt) == MISSING )
1317	!!data_out(:,:,k,itt) = AVE_geopt
1318	!!END WHERE
1319	!!end do loop_levels
1320	!!end do
1321
1322	END IF
1323
1324	!!! All other fields and geopt at higher levels come here
1325	do itt = 1, it
1326	do j = 1, iy
1327	do i = 1, ix
1328	do k = 1, kk
1329	if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,1,itt)
1330	end do
1331	do k = kk+1, num_metgrid_levels
1332	if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,iz,itt)
1333	end do
1334	end do
1335	end do
1336	end do
1337
1338
1339	END SUBROUTINE interp
1340	!------------------------------------------------------------------------------
1341	!--------------------------------------------------------
1342	SUBROUTINE int1D(xxout, xxin, ppin, ppout, npin, npout, LINLOG, MISSING)
1343
1344	! Modified from int2p - NCL code
1345	! routine to interpolate from one set of pressure levels
1346	! . to another set using linear or ln(p) interpolation
1347	!
1348	! NCL: xout = int2p (pin,xin,pout,linlog)
1349	! This code was originally written for a specific purpose.
1350	! . Several features were added for incorporation into NCL's
1351	! . function suite including linear extrapolation.
1352	!
1353	! nomenclature:
1354	!
1355	! . ppin - input pressure levels. The pin can be
1356	! . be in ascending or descending order
1357	! . xxin - data at corresponding input pressure levels
1358	! . npin - number of input pressure levels >= 2
1359	! . ppout - output pressure levels (input by user)
1360	! . same (ascending or descending) order as pin
1361	! . xxout - data at corresponding output pressure levels
1362	! . npout - number of output pressure levels
1363	! . linlog - if abs(linlog)=1 use linear interp in pressure
1364	! . if abs(linlog)=2 linear interp in ln(pressure)
1365	! . missing- missing data code.
1366
1367	! ! input types
1368	INTEGER :: npin,npout,linlog,ier
1369	real :: ppin(npin),xxin(npin),ppout(npout)
1370	real :: MISSING
1371	logical :: AVERAGE
1372	! ! output
1373	real :: xxout(npout)
1374	INTEGER :: j1,np,nl,nin,nlmax,nplvl
1375	INTEGER :: nlsave,np1,no1,n1,n2,nlstrt
1376	real :: slope,pa,pb,pc
1377
1378	! automatic arrays
1379	real :: pin(npin),xin(npin),p(npin),x(npin)
1380	real :: pout(npout),xout(npout)
1381
1382
1383	xxout = MISSING
1384	pout = ppout
1385	p = ppin
1386	x = xxin
1387	nlmax = npin
1388
1389	! exact p-level matches
1390	nlstrt = 1
1391	nlsave = 1
1392	do np = 1,npout
1393	xout(np) = MISSING
1394	do nl = nlstrt,nlmax
1395	if (pout(np).eq.p(nl)) then
1396	xout(np) = x(nl)
1397	nlsave = nl + 1
1398	go to 10
1399	end if
1400	end do
1401	10 nlstrt = nlsave
1402	end do
1403
1404	if (LINLOG.eq.1) then
1405	do np = 1,npout
1406	do nl = 1,nlmax - 1
1407	if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
1408	slope = (x(nl)-x(nl+1))/ (p(nl)-p(nl+1))
1409	xout(np) = x(nl+1) + slope* (pout(np)-p(nl+1))
1410	end if
1411	end do
1412	end do
1413	elseif (LINLOG.eq.2) then
1414	do np = 1,npout
1415	do nl = 1,nlmax - 1
1416	if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
1417	pa = log(p(nl))
1418	pb = log(pout(np))
1419	! special case: in case someone inadvertently enter p=0.
1420	if (p(nl+1).gt.0.d0) then
1421	pc = log(p(nl+1))
1422	else
1423	pc = log(1.d-4)
1424	end if
1425
1426	slope = (x(nl)-x(nl+1))/ (pa-pc)
1427	xout(np) = x(nl+1) + slope* (pb-pc)
1428	end if
1429	end do
1430	end do
1431	end if
1432
1433
1434	! place results in the return array;
1435	xxout = xout
1436
1437	END SUBROUTINE int1D
1438
1439	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1440	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1441	SUBROUTINE interp_1d( a, xa, na, b, xb, nb, vertical_type, MISSING)
1442
1443	implicit none
1444
1445	! Arguments
1446	integer, intent(in) :: na, nb
1447	real, intent(in), dimension(na) :: a, xa
1448	real, intent(in), dimension(nb) :: xb
1449	real, intent(out), dimension(nb) :: b
1450	character (len=1) :: vertical_type
1451
1452	!! Local variables
1453	real :: MISSING
1454	integer :: n_in, n_out
1455	real :: w1, w2
1456	logical :: interp
1457
1458	IF ( vertical_type == 'p' ) THEN
1459
1460	DO n_out = 1, nb
1461
1462	b(n_out) = MISSING
1463	interp = .false.
1464	n_in = 1
1465
1466	DO WHILE ( (.not.interp) .and. (n_in < na) )
1467	IF( (xa(n_in) >= xb(n_out)) .and. &
1468	(xa(n_in+1) <= xb(n_out)) ) THEN
1469	interp = .true.
1470	w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
1471	w2 = 1. - w1
1472	b(n_out) = w1a(n_in) + w2a(n_in+1)
1473	END IF
1474	n_in = n_in +1
1475	ENDDO
1476
1477	ENDDO
1478
1479	ELSE
1480
1481	DO n_out = 1, nb
1482
1483	b(n_out) = MISSING
1484	interp = .false.
1485	n_in = 1
1486
1487	DO WHILE ( (.not.interp) .and. (n_in < na) )
1488	IF( (xa(n_in) <= xb(n_out)) .and. &
1489	(xa(n_in+1) >= xb(n_out)) ) THEN
1490	interp = .true.
1491	w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
1492	w2 = 1. - w1
1493	b(n_out) = w1a(n_in) + w2a(n_in+1)
1494	END IF
1495	n_in = n_in +1
1496	ENDDO
1497
1498	ENDDO
1499
1500	END IF
1501
1502	END SUBROUTINE interp_1d
1503	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1504	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1505
1506	!------------------------------------------------------------------------------
1507	FUNCTION virtual (tmp,rmix)
1508	! This function returns virtual temperature in K, given temperature
1509	! in K and mixing ratio in kg/kg.
1510
1511	real :: tmp, rmix, virtual
1512
1513	virtual=tmp(0.622+rmix)/(0.622(1.+rmix))
1514
1515	END FUNCTION virtual
1516
1517	!------------------------------------------------------------------------------
1518	SUBROUTINE all_spaces ( command , length_of_char )
1519
1520	IMPLICIT NONE
1521
1522	INTEGER :: length_of_char
1523	CHARACTER (LEN=length_of_char) :: command
1524	INTEGER :: loop
1525
1526	DO loop = 1 , length_of_char
1527	command(loop:loop) = ' '
1528	END DO
1529
1530	END SUBROUTINE all_spaces
1531
1532	!------------------------------------------------------------------------------
1533	SUBROUTINE def_var (mcid, jvar, cval, itype, idm, jshape, order, desc, units, stag, coord, missing )
1534
1535	IMPLICIT NONE
1536
1537	INCLUDE 'netcdf.inc'
1538
1539	INTEGER :: mcid, jvar
1540	CHARACTER (LEN = 4) :: cval
1541	INTEGER :: itype, idm
1542	REAL, DIMENSION(6) :: jshape
1543	CHARACTER (LEN = 3) :: order
1544	CHARACTER (LEN = 19) :: desc
1545	CHARACTER (LEN = 2) :: units
1546	CHARACTER (LEN = 1) :: stag
1547	CHARACTER (LEN = 10) :: coord
1548
1549	INTEGER :: rcode, ilen
1550	CHARACTER (LEN=30) :: att_text
1551
1552	REAL :: missing
1553
1554
1555	IF ( itype == 5 ) THEN
1556	rcode = nf_redef(mcid)
1557	rcode = nf_def_var(mcid, trim(cval), NF_REAL, idm, jshape, jvar)
1558	rcode = nf_put_att_int(mcid, jvar, "FieldType", NF_INT, 1, 104)
1559	ENDIF
1560
1561	att_text = order
1562	ilen = len_trim(att_text)
1563	rcode = nf_put_att_text(mcid, jvar, "MemoryOrder", ilen, att_text(1:ilen) )
1564
1565	att_text = desc
1566	ilen = len_trim(att_text)
1567	rcode = nf_put_att_text(mcid, jvar, "description", ilen, att_text(1:ilen) )
1568
1569	att_text = units
1570	ilen = len_trim(att_text)
1571	rcode = nf_put_att_text(mcid, jvar, "units", ilen, att_text(1:ilen) )
1572
1573	att_text = stag
1574	ilen = len_trim(att_text)
1575	rcode = nf_put_att_text(mcid, jvar, "stagger", ilen, att_text(1:ilen) )
1576
1577	att_text = coord
1578	ilen = len_trim(att_text)
1579	rcode = nf_put_att_text(mcid, jvar, "coordinates", ilen, att_text(1:ilen) )
1580
1581	rcode = nf_put_att_real(mcid, jvar, "missing_value", NF_FLOAT, 1, MISSING )
1582
1583	rcode = nf_enddef(mcid)
1584
1585	END SUBROUTINE def_var
1586	!------------------------------------------------------------------------------
1587	!------------------------------------------------------------------------------
1588	!! Diagnostics: U & V on earth coordinates
1589	SUBROUTINE calc_uvmet( UUU, VVV, &
1590	UUUmet, VVVmet, &
1591	truelat1, truelat2, &
1592	stand_lon, map_proj, &
1593	longi, lati, &
1594	west_east_dim, south_north_dim, bottom_top_dim)
1595
1596	IMPLICIT NONE
1597
1598	!Arguments
1599	integer :: west_east_dim, south_north_dim, bottom_top_dim
1600	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: UUU
1601	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: VVV
1602	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: UUUmet
1603	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: VVVmet
1604	real, dimension(west_east_dim,south_north_dim) :: longi, lati
1605	real :: truelat1, truelat2, stand_lon
1606	integer :: map_proj
1607
1608	!Local
1609	integer :: i, j, k
1610	real :: cone
1611	real, dimension(west_east_dim,south_north_dim) :: diff, alpha
1612
1613	real, parameter :: PI = 3.141592653589793
1614	real, parameter :: DEG_PER_RAD = 180./PI
1615	real, parameter :: RAD_PER_DEG = PI/180.
1616
1617	IF ( map_proj .ge. 3 ) THEN ! No need to rotate
1618	!PRINT *, 'NO NEED TO ROTATE !!!! equivalent to output U,V with unstagger_grid'
1619	UUUmet(:,:,:) = UUU
1620	VVVmet(:,:,:) = VVV
1621	ELSE
1622
1623	cone = 1. ! PS
1624	IF ( map_proj .eq. 1) THEN ! Lambert Conformal mapping
1625	IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
1626	cone=(ALOG(COS(truelat1*RAD_PER_DEG))- &
1627	ALOG(COS(truelat2*RAD_PER_DEG))) / &
1628	(ALOG(TAN((90.-ABS(truelat1))RAD_PER_DEG0.5 ))- &
1629	ALOG(TAN((90.-ABS(truelat2))RAD_PER_DEG0.5 )) )
1630	ELSE
1631	cone = SIN(ABS(truelat1)*RAD_PER_DEG )
1632	ENDIF
1633	END IF
1634
1635
1636	diff = longi - stand_lon
1637	DO i = 1, west_east_dim
1638	DO j = 1, south_north_dim
1639	IF ( diff(i,j) .gt. 180. ) THEN
1640	diff(i,j) = diff(i,j) - 360.
1641	END IF
1642	IF ( diff(i,j) .lt. -180. ) THEN
1643	diff(i,j) = diff(i,j) + 360.
1644	END IF
1645	END DO
1646	END DO
1647
1648
1649	DO i = 1, west_east_dim
1650	DO j = 1, south_north_dim
1651	IF ( lati(i,j) .lt. 0. ) THEN
1652	alpha(i,j) = - diff(i,j) * cone * RAD_PER_DEG
1653	ELSE
1654	alpha(i,j) = diff(i,j) * cone * RAD_PER_DEG
1655	END IF
1656	END DO
1657	END DO
1658
1659
1660	DO k = 1,bottom_top_dim
1661	UUUmet(:,:,k) = VVV(:,:,k)sin(alpha) + UUU(:,:,k)cos(alpha)
1662	VVVmet(:,:,k) = VVV(:,:,k)cos(alpha) - UUU(:,:,k)sin(alpha)
1663	END DO
1664	END IF
1665
1666	END SUBROUTINE calc_uvmet
1667

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

Download in other formats:

Original Format