Context Navigation

api.F90 @ 73

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

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

Download in other formats:

Original Format