Context Navigation

api.F90 @ 87

Last change on this file since 87 was 87, checked in by aslmd, 14 years ago
MAJ 000-README-svn et 000-USERS. Version reference mars et mesoscale (avant nouvelle campagne de modification discutee en reunion).
File size: 69.5 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,*) "I used Cp, R : ",Cp,Rd
1131	write(6,*)
1132	write(6,*) "##########################################"
1133	write(6,*) " END of API PROGRAM - SUCCESS so far"
1134	write(6,*) "##########################################"
1135
1136	END PROGRAM api
1137	!---------------------------------------------------------------------
1138	!---------------------------------------------------------------------
1139	!---------------------------------------------------------------------
1140	SUBROUTINE handle_err(rcode)
1141	INTEGER rcode
1142	write(6,*) 'Error number ',rcode
1143	stop
1144	END SUBROUTINE
1145	!---------------------------------------------------------------------
1146	SUBROUTINE interp (data_out, data_in, pres_field, interp_levels, psfc, ter, tk, qv, ix, iy, iz, it, &
1147	num_metgrid_levels, LINLOG, extrapolate, GEOPT, MISSING)
1148
1149	INTEGER :: ix, iy, iz, it
1150	INTEGER :: num_metgrid_levels, LINLOG
1151	REAL, DIMENSION(ix, iy, num_metgrid_levels, it) :: data_out
1152	REAL, DIMENSION(ix, iy, iz, it) :: data_in, pres_field, tk, qv
1153	REAL, DIMENSION(ix, iy, it) :: psfc
1154	REAL, DIMENSION(ix, iy) :: ter
1155	REAL, DIMENSION(num_metgrid_levels) :: interp_levels
1156
1157	INTEGER :: i, j, itt, k, kk, kin
1158	REAL, DIMENSION(num_metgrid_levels) :: data_out1D
1159	REAL, DIMENSION(iz) :: data_in1D, pres_field1D
1160	INTEGER :: extrapolate
1161	REAL :: MISSING
1162	REAL, DIMENSION(ix, iy, num_metgrid_levels, it) :: N
1163	REAL :: sumA, sumN, AVE_geopt
1164	LOGICAL :: GEOPT
1165
1166	N = 1.0
1167
1168	do itt = 1, it
1169	!PRINT *, 'TIME... ', itt, ' in ', iz, ' out ', num_metgrid_levels
1170	do j = 1, iy
1171	do i = 1, ix
1172	data_in1D(:) = data_in(i,j,:,itt)
1173	pres_field1D(:) = pres_field(i,j,:,itt)
1174	IF ( LINLOG .le. 2 ) THEN
1175	CALL int1D (data_out1D, data_in1D, pres_field1D, interp_levels, iz, num_metgrid_levels, LINLOG, MISSING)
1176	ELSE
1177	CALL interp_1d (data_in1D, pres_field1D, iz, data_out1D, interp_levels, num_metgrid_levels, 'z', MISSING)
1178	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1179	!interp_1d (data_in1D, pres_field1D, num_metgrid_levels, data_out1D, interp_levels, iz, 'z')
1180	!CALL interp_1d( data_in_1d, z_data_1d, bottom_top_dim, data_out_1d, z_levs, number_of_zlevs, vertical_type)
1181	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1182	ENDIF
1183	data_out(i,j,:,itt) = data_out1D(:)
1184	end do
1185	end do
1186	end do
1187	!PRINT *, 'ok'
1188
1189	! Fill in missing values
1190	IF ( extrapolate == 0 ) RETURN !! no extrapolation - we are out of here
1191
1192
1193
1194	IF ( LINLOG .ge. 0 ) RETURN !! TEMPORARY: no extrapolation
1195	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1196	!! CAUTION BELOW: METHOD NOT ADAPTED TO MARS
1197	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1198
1199	!!! MARS MARS MARS
1200	!expon=287.04*.0065/9.81
1201	expon=192.*.0045/3.72
1202
1203	! First find where about 400 hPa is located
1204	kk = 0
1205	find_kk : do k = 1, num_metgrid_levels
1206	kk = k
1207	if ( interp_levels(k) <= 40000. ) exit find_kk
1208	end do find_kk
1209
1210
1211	IF ( GEOPT ) THEN !! geopt is treated different below ground
1212
1213	do itt = 1, it
1214	do k = 1, kk
1215	do j = 1, iy
1216	do i = 1, ix
1217	IF ( data_out(i,j,k,itt) == MISSING .AND. interp_levels(k) < psfc(i,j,itt) ) THEN
1218
1219	! We are below the first model level, but above the ground
1220
1221	!!! MARS MARS
1222	data_out(i,j,k,itt) = ((interp_levels(k) - pres_field(i,j,1,itt))ter(i,j)3.72 + &
1223	(psfc(i,j,itt) - interp_levels(k))*data_in(i,j,1,itt) ) / &
1224	(psfc(i,j,itt) - pres_field(i,j,1,itt))
1225
1226	ELSEIF ( data_out(i,j,k,itt) == MISSING ) THEN
1227
1228	! We are below both the ground and the lowest data level.
1229
1230	! First, find the model level that is closest to a "target" pressure
1231	! level, where the "target" pressure is delta-p less that the local
1232	! value of a horizontally smoothed surface pressure field. We use
1233	! delta-p = 150 hPa here. A standard lapse rate temperature profile
1234	! passing through the temperature at this model level will be used
1235	! to define the temperature profile below ground. This is similar
1236	! to the Benjamin and Miller (1990) method, except that for
1237	! simplicity, they used 700 hPa everywhere for the "target" pressure.
1238	! Code similar to what is implemented in RIP4
1239
1240
1241	!!! MARS MARS MARS
1242	ptarget = (psfc(i,j,itt)*.01) - 1.50 !150.
1243	dpmin=1.e2 !1.e4
1244	kupper = 0
1245	loop_kIN : do kin=iz,1,-1
1246	kupper = kin
1247	dp=abs( (pres_field(i,j,kin,itt)*.01) - ptarget )
1248	if (dp.gt.dpmin) exit loop_kIN
1249	dpmin=min(dpmin,dp)
1250	enddo loop_kIN
1251
1252	pbot=max(pres_field(i,j,1,itt),psfc(i,j,itt))
1253	zbot=min(data_in(i,j,1,itt)/3.72,ter(i,j)) !!MARS MARS
1254
1255	tbotextrap=tk(i,j,kupper,itt)(pbot/pres_field(i,j,kupper,itt))*expon
1256	tvbotextrap=virtual(tbotextrap,qv(i,j,1,itt))
1257
1258	data_out(i,j,k,itt) = (zbot+tvbotextrap/.0045(1.-(interp_levels(k)/pbot)expon))3.72 !!MARS MARS
1259
1260	ENDIF
1261	enddo
1262	enddo
1263	enddo
1264	enddo
1265
1266
1267	!!! Code for filling missing data with an average - we don't want to do this
1268	!!do itt = 1, it
1269	!!loop_levels : do k = 1, num_metgrid_levels
1270	!!sumA = SUM(data_out(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
1271	!!sumN = SUM(N(:,:,k,itt), MASK = data_out(:,:,k,itt) /= MISSING)
1272	!!IF ( sumN == 0. ) CYCLE loop_levels
1273	!!AVE_geopt = sumA/sumN
1274	!!WHERE ( data_out(:,:,k,itt) == MISSING )
1275	!!data_out(:,:,k,itt) = AVE_geopt
1276	!!END WHERE
1277	!!end do loop_levels
1278	!!end do
1279
1280	END IF
1281
1282	!!! All other fields and geopt at higher levels come here
1283	do itt = 1, it
1284	do j = 1, iy
1285	do i = 1, ix
1286	do k = 1, kk
1287	if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,1,itt)
1288	end do
1289	do k = kk+1, num_metgrid_levels
1290	if ( data_out(i,j,k,itt) == MISSING ) data_out(i,j,k,itt) = data_in(i,j,iz,itt)
1291	end do
1292	end do
1293	end do
1294	end do
1295
1296
1297	END SUBROUTINE interp
1298	!------------------------------------------------------------------------------
1299	!--------------------------------------------------------
1300	SUBROUTINE int1D(xxout, xxin, ppin, ppout, npin, npout, LINLOG, MISSING)
1301
1302	! Modified from int2p - NCL code
1303	! routine to interpolate from one set of pressure levels
1304	! . to another set using linear or ln(p) interpolation
1305	!
1306	! NCL: xout = int2p (pin,xin,pout,linlog)
1307	! This code was originally written for a specific purpose.
1308	! . Several features were added for incorporation into NCL's
1309	! . function suite including linear extrapolation.
1310	!
1311	! nomenclature:
1312	!
1313	! . ppin - input pressure levels. The pin can be
1314	! . be in ascending or descending order
1315	! . xxin - data at corresponding input pressure levels
1316	! . npin - number of input pressure levels >= 2
1317	! . ppout - output pressure levels (input by user)
1318	! . same (ascending or descending) order as pin
1319	! . xxout - data at corresponding output pressure levels
1320	! . npout - number of output pressure levels
1321	! . linlog - if abs(linlog)=1 use linear interp in pressure
1322	! . if abs(linlog)=2 linear interp in ln(pressure)
1323	! . missing- missing data code.
1324
1325	! ! input types
1326	INTEGER :: npin,npout,linlog,ier
1327	real :: ppin(npin),xxin(npin),ppout(npout)
1328	real :: MISSING
1329	logical :: AVERAGE
1330	! ! output
1331	real :: xxout(npout)
1332	INTEGER :: j1,np,nl,nin,nlmax,nplvl
1333	INTEGER :: nlsave,np1,no1,n1,n2,nlstrt
1334	real :: slope,pa,pb,pc
1335
1336	! automatic arrays
1337	real :: pin(npin),xin(npin),p(npin),x(npin)
1338	real :: pout(npout),xout(npout)
1339
1340
1341	xxout = MISSING
1342	pout = ppout
1343	p = ppin
1344	x = xxin
1345	nlmax = npin
1346
1347	! exact p-level matches
1348	nlstrt = 1
1349	nlsave = 1
1350	do np = 1,npout
1351	xout(np) = MISSING
1352	do nl = nlstrt,nlmax
1353	if (pout(np).eq.p(nl)) then
1354	xout(np) = x(nl)
1355	nlsave = nl + 1
1356	go to 10
1357	end if
1358	end do
1359	10 nlstrt = nlsave
1360	end do
1361
1362	if (LINLOG.eq.1) then
1363	do np = 1,npout
1364	do nl = 1,nlmax - 1
1365	if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
1366	slope = (x(nl)-x(nl+1))/ (p(nl)-p(nl+1))
1367	xout(np) = x(nl+1) + slope* (pout(np)-p(nl+1))
1368	end if
1369	end do
1370	end do
1371	elseif (LINLOG.eq.2) then
1372	do np = 1,npout
1373	do nl = 1,nlmax - 1
1374	if (pout(np).lt.p(nl) .and. pout(np).gt.p(nl+1)) then
1375	pa = log(p(nl))
1376	pb = log(pout(np))
1377	! special case: in case someone inadvertently enter p=0.
1378	if (p(nl+1).gt.0.d0) then
1379	pc = log(p(nl+1))
1380	else
1381	pc = log(1.d-4)
1382	end if
1383
1384	slope = (x(nl)-x(nl+1))/ (pa-pc)
1385	xout(np) = x(nl+1) + slope* (pb-pc)
1386	end if
1387	end do
1388	end do
1389	end if
1390
1391
1392	! place results in the return array;
1393	xxout = xout
1394
1395	END SUBROUTINE int1D
1396
1397	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1398	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1399	SUBROUTINE interp_1d( a, xa, na, b, xb, nb, vertical_type, MISSING)
1400
1401	implicit none
1402
1403	! Arguments
1404	integer, intent(in) :: na, nb
1405	real, intent(in), dimension(na) :: a, xa
1406	real, intent(in), dimension(nb) :: xb
1407	real, intent(out), dimension(nb) :: b
1408	character (len=1) :: vertical_type
1409
1410	!! Local variables
1411	real :: MISSING
1412	integer :: n_in, n_out
1413	real :: w1, w2
1414	logical :: interp
1415
1416	IF ( vertical_type == 'p' ) THEN
1417
1418	DO n_out = 1, nb
1419
1420	b(n_out) = MISSING
1421	interp = .false.
1422	n_in = 1
1423
1424	DO WHILE ( (.not.interp) .and. (n_in < na) )
1425	IF( (xa(n_in) >= xb(n_out)) .and. &
1426	(xa(n_in+1) <= xb(n_out)) ) THEN
1427	interp = .true.
1428	w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
1429	w2 = 1. - w1
1430	b(n_out) = w1a(n_in) + w2a(n_in+1)
1431	END IF
1432	n_in = n_in +1
1433	ENDDO
1434
1435	ENDDO
1436
1437	ELSE
1438
1439	DO n_out = 1, nb
1440
1441	b(n_out) = MISSING
1442	interp = .false.
1443	n_in = 1
1444
1445	DO WHILE ( (.not.interp) .and. (n_in < na) )
1446	IF( (xa(n_in) <= xb(n_out)) .and. &
1447	(xa(n_in+1) >= xb(n_out)) ) THEN
1448	interp = .true.
1449	w1 = (xa(n_in+1)-xb(n_out))/(xa(n_in+1)-xa(n_in))
1450	w2 = 1. - w1
1451	b(n_out) = w1a(n_in) + w2a(n_in+1)
1452	END IF
1453	n_in = n_in +1
1454	ENDDO
1455
1456	ENDDO
1457
1458	END IF
1459
1460	END SUBROUTINE interp_1d
1461	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1462	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1463
1464	!------------------------------------------------------------------------------
1465	FUNCTION virtual (tmp,rmix)
1466	! This function returns virtual temperature in K, given temperature
1467	! in K and mixing ratio in kg/kg.
1468
1469	real :: tmp, rmix, virtual
1470
1471	virtual=tmp(0.622+rmix)/(0.622(1.+rmix))
1472
1473	END FUNCTION virtual
1474
1475	!------------------------------------------------------------------------------
1476	SUBROUTINE all_spaces ( command , length_of_char )
1477
1478	IMPLICIT NONE
1479
1480	INTEGER :: length_of_char
1481	CHARACTER (LEN=length_of_char) :: command
1482	INTEGER :: loop
1483
1484	DO loop = 1 , length_of_char
1485	command(loop:loop) = ' '
1486	END DO
1487
1488	END SUBROUTINE all_spaces
1489
1490	!------------------------------------------------------------------------------
1491	SUBROUTINE def_var (mcid, jvar, cval, itype, idm, jshape, order, desc, units, stag, coord, missing )
1492
1493	IMPLICIT NONE
1494
1495	INCLUDE 'netcdf.inc'
1496
1497	INTEGER :: mcid, jvar
1498	CHARACTER (LEN = 4) :: cval
1499	INTEGER :: itype, idm
1500	REAL, DIMENSION(6) :: jshape
1501	CHARACTER (LEN = 3) :: order
1502	CHARACTER (LEN = 19) :: desc
1503	CHARACTER (LEN = 2) :: units
1504	CHARACTER (LEN = 1) :: stag
1505	CHARACTER (LEN = 10) :: coord
1506
1507	INTEGER :: rcode, ilen
1508	CHARACTER (LEN=30) :: att_text
1509
1510	REAL :: missing
1511
1512
1513	IF ( itype == 5 ) THEN
1514	rcode = nf_redef(mcid)
1515	rcode = nf_def_var(mcid, trim(cval), NF_REAL, idm, jshape, jvar)
1516	rcode = nf_put_att_int(mcid, jvar, "FieldType", NF_INT, 1, 104)
1517	ENDIF
1518
1519	att_text = order
1520	ilen = len_trim(att_text)
1521	rcode = nf_put_att_text(mcid, jvar, "MemoryOrder", ilen, att_text(1:ilen) )
1522
1523	att_text = desc
1524	ilen = len_trim(att_text)
1525	rcode = nf_put_att_text(mcid, jvar, "description", ilen, att_text(1:ilen) )
1526
1527	att_text = units
1528	ilen = len_trim(att_text)
1529	rcode = nf_put_att_text(mcid, jvar, "units", ilen, att_text(1:ilen) )
1530
1531	att_text = stag
1532	ilen = len_trim(att_text)
1533	rcode = nf_put_att_text(mcid, jvar, "stagger", ilen, att_text(1:ilen) )
1534
1535	att_text = coord
1536	ilen = len_trim(att_text)
1537	rcode = nf_put_att_text(mcid, jvar, "coordinates", ilen, att_text(1:ilen) )
1538
1539	rcode = nf_put_att_real(mcid, jvar, "missing_value", NF_FLOAT, 1, MISSING )
1540
1541	rcode = nf_enddef(mcid)
1542
1543	END SUBROUTINE def_var
1544	!------------------------------------------------------------------------------
1545	!------------------------------------------------------------------------------
1546	!! Diagnostics: U & V on earth coordinates
1547	SUBROUTINE calc_uvmet( UUU, VVV, &
1548	UUUmet, VVVmet, &
1549	truelat1, truelat2, &
1550	stand_lon, map_proj, &
1551	longi, lati, &
1552	west_east_dim, south_north_dim, bottom_top_dim)
1553
1554	IMPLICIT NONE
1555
1556	!Arguments
1557	integer :: west_east_dim, south_north_dim, bottom_top_dim
1558	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: UUU
1559	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: VVV
1560	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: UUUmet
1561	real, dimension(west_east_dim,south_north_dim,bottom_top_dim) :: VVVmet
1562	real, dimension(west_east_dim,south_north_dim) :: longi, lati
1563	real :: truelat1, truelat2, stand_lon
1564	integer :: map_proj
1565
1566	!Local
1567	integer :: i, j, k
1568	real :: cone
1569	real, dimension(west_east_dim,south_north_dim) :: diff, alpha
1570
1571	real, parameter :: PI = 3.141592653589793
1572	real, parameter :: DEG_PER_RAD = 180./PI
1573	real, parameter :: RAD_PER_DEG = PI/180.
1574
1575	! print *, 'map ', map_proj
1576	IF ( map_proj .ge. 3 ) THEN ! No need to rotate
1577	!PRINT *, 'NO NEED TO ROTATE !!!! equivalent to output U,V with unstagger_grid'
1578	UUUmet(:,:,:) = UUU
1579	VVVmet(:,:,:) = VVV
1580	ELSE
1581	!END IF
1582
1583	cone = 1. ! PS
1584	IF ( map_proj .eq. 1) THEN ! Lambert Conformal mapping
1585	IF (ABS(truelat1-truelat2) .GT. 0.1) THEN
1586	cone=(ALOG(COS(truelat1*RAD_PER_DEG))- &
1587	ALOG(COS(truelat2*RAD_PER_DEG))) / &
1588	(ALOG(TAN((90.-ABS(truelat1))RAD_PER_DEG0.5 ))- &
1589	ALOG(TAN((90.-ABS(truelat2))RAD_PER_DEG0.5 )) )
1590	ELSE
1591	cone = SIN(ABS(truelat1)*RAD_PER_DEG )
1592	ENDIF
1593	END IF
1594
1595
1596	diff = longi - stand_lon
1597	DO i = 1, west_east_dim
1598	DO j = 1, south_north_dim
1599	IF ( diff(i,j) .gt. 180. ) THEN
1600	diff(i,j) = diff(i,j) - 360.
1601	END IF
1602	IF ( diff(i,j) .lt. -180. ) THEN
1603	diff(i,j) = diff(i,j) + 360.
1604	END IF
1605	END DO
1606	END DO
1607
1608	!print *, longi(10,10)
1609	!print *, lati(10,10)
1610
1611
1612	DO i = 1, west_east_dim
1613	DO j = 1, south_north_dim
1614	IF ( lati(i,j) .lt. 0. ) THEN
1615	alpha(i,j) = - diff(i,j) * cone * RAD_PER_DEG
1616	ELSE
1617	alpha(i,j) = diff(i,j) * cone * RAD_PER_DEG
1618	END IF
1619	END DO
1620	END DO
1621
1622
1623	DO k = 1,bottom_top_dim
1624	UUUmet(:,:,k) = VVV(:,:,k)sin(alpha) + UUU(:,:,k)cos(alpha)
1625	VVVmet(:,:,k) = VVV(:,:,k)cos(alpha) - UUU(:,:,k)sin(alpha)
1626	!print *,UUU(10,10,k), UUUmet(10,10,k)
1627	!print *,VVV(10,10,k), VVVmet(10,10,k)
1628	END DO
1629	END IF
1630
1631	END SUBROUTINE calc_uvmet
1632

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

Download in other formats:

Original Format