Context Navigation

api.F90 @ 1228

Last change on this file since 1228 was 763, checked in by acolaitis, 12 years ago

###################################################
# PYTHON / PLANETPLOT #
###################################################

# ------------------- XY plots ------------------ #

# Added a new category of plot to unidim, contour, etc... called "xy"
# - xy plots are plots that do not use time,vert,lat or lon as axis
# - variables to be plotted are stored in plot_x and plot_y, which is done in
# select_getfield. (there is no "what_I_plot" var for "wy" plots)
# - plot_x and plot_y are also subject to reduce field. A None value indicates
# the plot is not 'xy'
# - "xy" plots are used for a specific subset of plots : histograms, fourier
# transforms, hodographs
# - added option --analysis to perform certain kinds of analysis on the data
# (corresponding to xy plots).
# - One selects the fields he wants to plot (e.g. -v UV --lon 0 --lat 20) and
# chooses a kind of analysis :

--analysis fft # in the particular case given above, this mode corresponds

# to the mean of the ampitude spectrum of the vertical spatial
# fast fourier transform taken at each time index
# note that for now, fft are always done along the vertical
# axis. This could be made more flexible.
# not that the minimum wavelength you can plot depends on the
# vertical step of your simulation. THIS STEP HAS TO BE
# CONSTANT, hence you MUST use API or ZRECAST with constant
# spacing.

--analysis histo # histogram on the flattened data. If the user asks for --lon 0,20,

# the average is done before computing the histogram (contrary to the fft).
# However, if given a 2D array (with only --lon and --lat on a
# 4D field for example), the data is flattened before computation
# so that the result is still a 1D histogram.

--analysis histodensity # histogram with a kernel density estimate to a

# gaussian distribution, also giving the mean, variance,
# skewness and kurtosis. Other distributions are available in
# the scipy.stats package and could be implemented.

# - added variables "-v hodograph" and "-v hodograph_2". First one is a
# regular hodograph with "u" and "v" as axis, with labels of the local times
# (use --axtime lt). This is a "xy" plot (you must specify a vertical level
# as well, usually --vert 0 with an interpolation at 5m (-i 4 -l 0.005)).
# Second one is the variation with local time (for exemple) of the wind
# rotation (arctan(v/u)). This is not a "xy" plot but "unidim".
# For --ope plots in "xy" cases, only the operation plot is displayed.

# ------------------- Operations ------------------ #

# - For operations --ope -, the histogram (fourth plot) has been removed. To get it
# back, call --ope -_histo.
# - _only has been added to "+","-","-%" operations (eg "-_only")
# - For operations "-","+","-%","-_only","+_only","-%_only","-_histo", it is now
# possible to call multiple files (sill one variable) and compare each of them
# with the unique given reference file. Ex: -f file1,file2 --fref file3 --ope - will give 6 plots:

file1 file3 file1-file3
file2 file3 file2-file3

# In the case of "xy" plots, the multiple operation plots are regrouped on a
# single plot (using multiple lines). the title of this plot is not 'fig(2)-fig(1)'
# (default) but the argument of the --title command. Labels have to be given
# as following : --labels "dummy","dummy","file1-file3 label","dummy","dummy","file2-file3 label" (dummy can be anything. this is to be improved)
# To be able to run on multiple files and easily introduce the correct number and order of plots, these operations have been moved outside of the main loop
# on namefiles and variables. Operation of the type "add_var" or "cat" have
# been left inside the loop and are unchanged.

# ------------------- Localtime ------------------ #

# Changed the way localtime is computed. Reasons:
# - it was assuming one timestep per hour in computations mixing indices and
# actual times
# - to determine the starting date, it was using the name of the file (for Meso), instead of using the netcdf
# attribute (START_DATE)
# - it was using the computed mean longitude of the domain, which is not correct for
# hemispheric domain. => better to use the netcdf attribute CEN_LON
# - it was using this mean longitude even for profile plots at given
# longitude => better to get the local time at this given lon, especially
# important for large domains
# => new localtime() is in myplot (old one is commented). Interv is obsolete (but not removed yet). Case "Geo" has not been looked at.
# new localtime in myplot correctly account for starting date of the file in
# all cases
# accounts for local longitude of the plot
# accounts for files that do not have per hour outputs, but per timestep.
# specific cases can be added in myplot in localtime()

# ------------------ Misc ------------------ #

# - added option --xlog to get x logarithmic axis (--ylog already existed)
# - added the possibility to use 2 files of different gridding (-f
# file1,file2) although of the same type (meso for exemple). For that purpose,
# lon, lat, alt and vert arrays are now indexed with 'index_f', as for
# all_var. => all_lon, all_lat, all_lat, all_vert
# - added function teta_to_tk in myplot so that a call to pp.py on standard
# LMD mesoscale file with "-v tk" can be done without the need to call API.
# The temperature is computed from T and PTOT, knowing P0, T0 and R_CP.
# - bug correction in determineplot() that was causing wrong plot number and
# slices number when not calling averaged lon, lat, vert or times. (which is
# often !)
# - added new options to redope: --redope edge_x1, edge_x2, edge_y1, edge_y2
# which plots the boundaries of the domain. This is different compared to
# asking for a fixed longitude, because the domain boundary might not be at constant
# longitude (hemispheric domain for exemple). x1 is the western boundary, x2
# the eastern, y1 the southern and y2 the northern. x1 and x2 reduce the
# dimension along --lon, y1 and y2 reduce the dimension along --lat.
# - added control in windamplitude() to determine whether winds are staggered or
# unstaggered. This is usefull when dealing with non LMD_MMM files.
# - corrected a bug in reduce_field where the mean was computed on the wrong
# axis !!! (pretty serious bug)

# Exemple of plots you can do with these new options can be found in
$YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/Plots_MasterScript/bam.sh

# ------------------ API ------------------ #

# - changed maximum of levels from 299 to 1000 in API (interpolation on 1000 levels is
# usefull to get larger bandwidth in fourier transform)
# the following concerns users of MRAMS files.
# - API has not been modified for MRAMS files. Instead, a python script
# (ic.py) is run on MRAMS .ctl and .dat files, which automatically format those files
# to be API and pp.py compatible.
# - ic.py is in $YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/File_conversion

###################################################
# INTERCOMPARISON TOOLS: #
###################################################

#ic.py in
$YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/File_conversion
#CDO installer with import_binary in
$YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/CDO
#Plotting scripts in
$YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/Plots_MasterScript
#1D sensibility tool in
$YOUR_SVN/trunk/UTIL/PYTHON/Intercomparison/1D_sensibility_tool

# See README in each of these folders for details.

File size: 71.6 KB

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

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

Download in other formats:

Original Format