Context Navigation

zrecast.F90 @ 137

Last change on this file since 137 was 137, checked in by aslmd, 14 years ago

ATTENTION COMMIT MAJEUR NON-CONSERVATIF

CHANGEMENT ARBORESCENCE POUR SEPARATION CLAIRE DES COMPOSANTES et MODELES
UTILISATION DE LIENS SYMBOLIQUES POUR GARDER UNE BASE COMMUNE LMDZ.COMMON
VOIR LE FICHIER DOC/000-MODELS POUR EN SAVOIR PLUS !

EN CAS DE PROBLEME AVEC svn update IL EST CONSEILLE DE REVENIR A UN svn co

DERNIERE REVISION AVEC L'ANCIENNE ARBORESCENCE : 132

File size: 82.0 KB

Line
1
2
3	program zrecast
4
5	! This program reads 4D (lon-lat-alt-time) fields from GCM output files
6	! (ie: diagfi.nc time series or concat.nc or stats.nc files) and, by
7	! integrating the hydrostatic equation, recasts data along the vertical
8	! direction.
9	! The vertical coordinate can be either "above areoid altitudes" or
10	! "pressure". Some interpolation along the vertical direction is also
11	! done, following instructions given by user.
12	! For "above areoid altitudes" output, Atmospheric pressure is added to
13	! output dataset; for "pressure coordinate" outputs, the above areoid
14	! altitude of pressure is added to output dataset.
15	!
16	! Minimal requirements and dependencies:
17	! The dataset must include the following data:
18	! - surface pressure
19	! - atmospheric temperature
20	! - hybrid coordinates aps() and bps(), or sigma levels() (see section 1.3.2)
21	! - ground geopotential (in input file; if not found, it is sought
22	! in a 'diagfi.nc' file. If not found there, it is then sought in
23	! a 'phisinit.nc' file (see section 1.3.3 of program)
24	!
25	! - When integration the hydrostatic equation, we assume that R, the molecular
26	! Gas Constant, may not be constant, so it is computed as
27	! R=P/(rho*T) (P=Pressure, rho=density, T=temperature)
28	! If 'rho' is not available, then we use a constant R (see section 2.2)
29	!
30	! WARNING: Asking for many points along the vertical direction quickly
31	! leads to HUGE output files.
32	!
33	! EM 01/2006 : Corrected a bug in vertical (log) interpolation for pressure
34	! and density
35	! EM 10/2006 : Modified program so that it can now process 'stats.nc'
36	! files obtained from British GCM (ie: vertical coordinate
37	! given as sigma levels and geopotential read from file
38	! 'phisinit.nc')
39	! EM 02/2007 : Changed behavior for "altitude above surface" case
40	! (for MCD RMS computation). Number of levels is then set as
41	! number of levels in initial file,
42	! and the new set of above surface levels follow a more elaborate
43	! distribution (see build_zs.F90 routine).
44	! EM 08/2009 : User may now specify values of each vertical level,
45	! or just min,max and number of levels (as before)
46	! EM 01/2010 : Corrected bug in 'zs_coord_interp' to correctly handle the case
47	! when interpolating log-wise and the density field is not
48	! available.
49
50	implicit none
51
52	include "netcdf.inc" ! NetCDF definitions
53
54	character (len=128) :: infile ! input file name (diagfi.nc or stats.nc format)
55	character (len=128) :: infile2 ! second input file (may be needed for 'phisini')
56	character (len=128) :: outfile ! output file name
57
58	character (len=64) :: text ! to store some text
59	character (len=64) :: tmpvarname ! temporarily store a variable name
60	integer tmpvarid ! temporarily store a variable ID
61	integer tmpdimid ! temporarily store a dimension ID
62	integer tmpndims ! temporarily store # of dimensions of a variable
63	integer infid ! NetCDF input file ID (of diagfi.nc or stats.nc format)
64	integer infid2 ! NetCDF input file which contains 'phisini' dataset (diagfi.nc)
65	integer nbvarinfile ! # of variables in input file
66	integer nbattr ! # of attributes of a given variable in input file
67	integer nbvar4dinfile ! # of 4D (lon,lat,alt,time) variables in input file
68	integer outfid ! NetCDF output file ID
69	integer lon_dimid,lat_dimid,alt_dimid,time_dimid ! NetCDF dimension IDs
70	integer lon_varid,lat_varid,alt_varid,time_varid
71	integer gcm_layers_dimid ! NetCDF dimension ID for # of layers in GCM
72	integer sigma_varid,aps_varid,bps_varid
73	integer za_varid,p_varid ! above areoid and pressure data IDs
74
75
76
77	integer ps_varid ! surface pressure data ID
78	integer,dimension(4) :: datashape ! shape of 4D datasets
79	integer,dimension(3) :: surfdatashape ! shape of 3D (surface+time) datasets
80
81	real :: miss_val=-9.99e+33 ! special "missing value" to specify missing data
82	real,parameter :: miss_val_def=-9.99e+33 ! default value for "missing value"
83	character (len=64), dimension(:), allocatable :: var
84	! var(): names of variables that will be processed
85	integer nbvar ! # of variables to process
86	integer,dimension(:),allocatable :: var_id ! IDs of variables var() (in outfile)
87	real,dimension(:),allocatable :: lon ! longitude
88	integer lonlength ! # of grid points along longitude
89	real,dimension(:),allocatable :: lat ! latitude
90	integer latlength ! # of grid points along latitude
91	integer altlength ! # of grid point along altitude (of input datasets)
92	real,dimension(:),allocatable :: time ! time
93	integer timelength ! # of points along time
94	real,dimension(:),allocatable :: aps,bps ! hybrid vertical coordinates
95	real,dimension(:),allocatable :: sigma ! sigma levels
96	real,dimension(:,:),allocatable :: phisinit ! Ground geopotential
97	real,dimension(:,:,:),allocatable :: ps ! GCM surface pressure
98	real,dimension(:,:,:,:),allocatable :: press ! GCM atmospheric pressure
99	real,dimension(:,:,:,:),allocatable :: temp ! GCM atmospheric temperature
100	real,dimension(:,:,:,:),allocatable :: rho ! GCM atmospheric density
101	real,dimension(:,:,:,:),allocatable :: za_gcm ! GCM above areoid levels (m)
102	real,dimension(:,:,:,:),allocatable :: zs_gcm ! GCM above surface heights (m)
103	real,dimension(:,:,:,:),allocatable :: indata ! to store a GCM 4D dataset
104	real,dimension(:,:,:,:),allocatable :: outdata ! to store a 4D dataset
105	integer ierr ! NetCDF routines return code
106	integer i,j,ilon,ilat,ilev,itim ! for loops
107
108	integer ztype ! Flag for vertical coordinate of output
109	! ztype=1: pressure ztype=2: above areoid ztype=3: above local surface
110	integer nblev ! # of levels (along vertical coordinate) for output data
111	real pmin,pmax ! min and max values for output pressure coordinate
112	real,dimension(:),allocatable :: plevel ! Pressure levels for output data
113	real zamin,zamax ! min and max values for output above areoid coordinate
114	real,dimension(:),allocatable :: zareoid ! Above areoid heights for output data
115	real,dimension(:),allocatable :: zsurface ! Above surface heights for output
116	logical :: have_rho ! Flag: true if density 'rho' is available
117	logical :: have_sigma ! Flag: true if sigma levels are known (false if hybrid
118	! coordinates are used)
119	logical :: auto_vert_levels ! Flag: true if the positions of vertical levels
120	! has to be computed; false if these are given
121	! by the user
122	!===============================================================================
123	! 1. Input parameters
124	!===============================================================================
125
126	!===============================================================================
127	! 1.1 Input file
128	!===============================================================================
129
130	write(,) ""
131	write(,) " Program valid for diagfi.nc, concatnc.nc and stats.nc files"
132	write(,) "Enter input file name:"
133
134	read(*,'(a128)') infile
135	write(,) ""
136
137	! open input file
138
139	ierr = NF_OPEN(infile,NF_NOWRITE,infid)
140	if (ierr.ne.NF_NOERR) then
141	write(,) 'ERROR: Pb opening file ',trim(infile)
142	stop ""
143	endif
144
145	!===============================================================================
146	! 1.2 Get # and names of variables in input file
147	!===============================================================================
148
149	ierr=NF_INQ_NVARS(infid,nbvarinfile)
150	if (ierr.ne.NF_NOERR) then
151	write(,) 'ERROR: Failed geting number of variables from file'
152	stop
153	endif
154
155	write(,)" The following variables have been found:"
156	nbvar4dinfile=0
157	do i=1,nbvarinfile
158	! get name of variable # i
159	ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
160	! check if it is a 4D variable
161	ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
162	if (tmpndims.eq.4) then
163	nbvar4dinfile=nbvar4dinfile+1
164	write(,) trim(tmpvarname)
165	endif
166	enddo
167
168	allocate(var(nbvar4dinfile))
169
170	write(,) ""
171	write(,) "Which variable do you want to keep?"
172	write(,) "all or list of <variables> (separated by <Return>s)"
173	write(,) "(an empty line , i.e: just <Return>, implies end of list)"
174	nbvar=0
175	read(*,'(a64)') tmpvarname
176	do while ((tmpvarname.ne.' ').and.(trim(tmpvarname).ne.'all'))
177	! check if tmpvarname is valid
178	ierr=NF_INQ_VARID(infid,tmpvarname,tmpvarid)
179	if (ierr.eq.NF_NOERR) then ! valid name
180	! also check that it is indeed a 4D variable
181	ierr=NF_INQ_VARNDIMS(infid,tmpvarid,tmpndims)
182	if (tmpndims.eq.4) then
183	nbvar=nbvar+1
184	var(nbvar)=tmpvarname
185	else
186	write(,) 'Error: ',trim(tmpvarname),' is not a 4D variable'
187	write(,) ' (we''ll skip that one)'
188	endif
189	else ! invalid name
190	write(,) 'Error: ',trim(tmpvarname),' is not a valid name'
191	write(,) ' (we''ll skip that one)'
192	endif
193	read(*,'(a64)') tmpvarname
194	enddo
195
196	! handle "all" case
197	if (tmpvarname.eq.'all') then
198	nbvar=0
199	do i=1,nbvarinfile
200	! look for 4D variables
201	ierr=NF_INQ_VARNDIMS(infid,i,tmpndims)
202	if (tmpndims.eq.4) then
203	nbvar=nbvar+1
204	! get the corresponding name
205	ierr=NF_INQ_VARNAME(infid,i,tmpvarname)
206	var(nbvar)=tmpvarname
207	endif
208	enddo
209	endif
210
211	! Check that there is at least 1 variable to process
212	if (nbvar.eq.0) then
213	write(,) 'No variables to process !?'
214	write(,) 'Might as well stop here'
215	stop ""
216	else
217	write(,) ""
218	write(,) 'OK, the following variables will be processed:'
219	do i=1,nbvar
220	write(,) var(i)
221	enddo
222	endif
223
224	!===============================================================================
225	! 1.3 Get grids in lon,lat,alt,time,
226	! as well as hybrid coordinates aps() and bps() (or sigma levels sigma())
227	! and phisinit() from input file
228	!===============================================================================
229
230	! 1.3.1 longitude, latitude, altitude and time
231
232	! latitude
233	ierr=NF_INQ_DIMID(infid,"latitude",tmpdimid)
234	if (ierr.ne.NF_NOERR) then
235	write(,) "Could not get latitude dimension ID"
236	write(,) " looking for lat dimension instead... "
237	ierr=NF_INQ_DIMID(infid,"lat",tmpdimid)
238	if (ierr.ne.NF_NOERR) then
239	stop "Error: Failed to get lat dimension ID"
240	else
241	ierr=NF_INQ_VARID(infid,"lat",tmpvarid)
242	if (ierr.ne.NF_NOERR) then
243	stop "Error: Failed to get lat ID"
244	else
245	ierr=NF_INQ_DIMLEN(infid,tmpdimid,latlength)
246	if (ierr.ne.NF_NOERR) then
247	stop "Error: Failed to get lat length"
248	else
249	allocate(lat(latlength))
250	ierr=NF_GET_VAR_REAL(infid,tmpvarid,lat)
251	if (ierr.ne.NF_NOERR) then
252	stop "Error: Failed reading lat"
253	endif
254	endif
255	endif
256	endif
257	else
258	ierr=NF_INQ_VARID(infid,"latitude",tmpvarid)
259	if (ierr.ne.NF_NOERR) then
260	stop "Error: Failed to get latitude ID"
261	else
262	ierr=NF_INQ_DIMLEN(infid,tmpdimid,latlength)
263	if (ierr.ne.NF_NOERR) then
264	stop "Error: Failed to get latitude length"
265	else
266	allocate(lat(latlength))
267	ierr=NF_GET_VAR_REAL(infid,tmpvarid,lat)
268	if (ierr.ne.NF_NOERR) then
269	stop "Error: Failed reading latitude"
270	endif
271	endif
272	endif
273	endif
274
275	! longitude
276	ierr=NF_INQ_DIMID(infid,"longitude",tmpdimid)
277	if (ierr.ne.NF_NOERR) then
278	write(,) "Could not get longitude dimension ID"
279	write(,) " looking for lon dimension instead... "
280	ierr=NF_INQ_DIMID(infid,"lon",tmpdimid)
281	if (ierr.ne.NF_NOERR) then
282	stop "Error: Failed to get lon dimension ID"
283	else
284	ierr=NF_INQ_VARID(infid,"lon",tmpvarid)
285	if (ierr.ne.NF_NOERR) then
286	stop "Error: Failed to get lon ID"
287	else
288	ierr=NF_INQ_DIMLEN(infid,tmpdimid,lonlength)
289	if (ierr.ne.NF_NOERR) then
290	stop "Error: Failed to get lon length"
291	else
292	allocate(lon(lonlength))
293	ierr=NF_GET_VAR_REAL(infid,tmpvarid,lon)
294	if (ierr.ne.NF_NOERR) then
295	stop "Error: Failed reading longitude"
296	endif
297	endif
298	endif
299	endif
300	else
301	ierr=NF_INQ_VARID(infid,"longitude",tmpvarid)
302	if (ierr.ne.NF_NOERR) then
303	stop "Error: Failed to get longitude ID"
304	else
305	ierr=NF_INQ_DIMLEN(infid,tmpdimid,lonlength)
306	if (ierr.ne.NF_NOERR) then
307	stop "Error: Failed to get longitude length"
308	else
309	allocate(lon(lonlength))
310	ierr=NF_GET_VAR_REAL(infid,tmpvarid,lon)
311	if (ierr.ne.NF_NOERR) then
312	stop "Error: Failed reading longitude"
313	endif
314	endif
315	endif
316	endif
317
318	! time
319	ierr=NF_INQ_DIMID(infid,"Time",tmpdimid)
320	if (ierr.ne.NF_NOERR) then
321	write(,) "Could not get Time dimension ID"
322	write(,) " looking for time dimension instead... "
323	ierr=NF_INQ_DIMID(infid,"time",tmpdimid)
324	if (ierr.ne.NF_NOERR) then
325	stop "Error: Failed to get lon dimension ID"
326	else
327	ierr=NF_INQ_VARID(infid,"time",tmpvarid)
328	if (ierr.ne.NF_NOERR) then
329	stop "Error: Failed to get time ID"
330	else
331	ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength)
332	if (ierr.ne.NF_NOERR) then
333	stop "Error: Failed to get time length"
334	else
335	allocate(time(timelength))
336	ierr=NF_GET_VAR_REAL(infid,tmpvarid,time)
337	if (ierr.ne.NF_NOERR) then
338	stop "Error: Failed reading time"
339	endif
340	endif
341	endif
342	endif
343	else
344	ierr=NF_INQ_VARID(infid,"Time",tmpvarid)
345	if (ierr.ne.NF_NOERR) then
346	stop "Error: Failed to get Time ID"
347	else
348	ierr=NF_INQ_DIMLEN(infid,tmpdimid,timelength)
349	if (ierr.ne.NF_NOERR) then
350	stop "Error: Failed to get Time length"
351	else
352	allocate(time(timelength))
353	ierr=NF_GET_VAR_REAL(infid,tmpvarid,time)
354	if (ierr.ne.NF_NOERR) then
355	stop "Error: Failed reading Time"
356	endif
357	endif
358	endif
359	endif
360
361	! altlength
362	ierr=NF_INQ_DIMID(infid,"altitude",tmpdimid)
363	if (ierr.ne.NF_NOERR) then
364	write(,) "Could not get altitude dimension ID"
365	write(,) " looking for sigma dimension instead... "
366	ierr=NF_INQ_DIMID(infid,"sigma",tmpdimid)
367	if (ierr.ne.NF_NOERR) then
368	stop "Error: Failed to get sigma dimension ID"
369	else
370	ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
371	if (ierr.ne.NF_NOERR) then
372	stop "Error: Failed to get altitude length"
373	endif
374	endif
375	else
376	ierr=NF_INQ_DIMLEN(infid,tmpdimid,altlength)
377	if (ierr.ne.NF_NOERR) then
378	stop "Error: Failed to get altitude length"
379	endif
380	endif
381
382	! 1.3.2 Get hybrid coordinates (or sigma levels)
383
384	! start by looking for sigma levels
385	ierr=NF_INQ_VARID(infid,"sigma",tmpvarid)
386	if (ierr.ne.NF_NOERR) then
387	have_sigma=.false.
388	write(,) "Could not find sigma levels... will look for hybrid coordinates"
389	else
390	have_sigma=.true.
391	allocate(sigma(altlength))
392	ierr=NF_GET_VAR_REAL(infid,tmpvarid,sigma)
393	if (ierr.ne.NF_NOERR) then
394	stop "Error: Failed reading sigma"
395	endif
396	endif
397
398	! if no sigma levels, look for hybrid coordinates
399	if (.not.have_sigma) then
400	! hybrid coordinate aps
401	ierr=NF_INQ_VARID(infid,"aps",tmpvarid)
402	if (ierr.ne.NF_NOERR) then
403	stop "Error: Failed to get aps ID"
404	else
405	allocate(aps(altlength))
406	ierr=NF_GET_VAR_REAL(infid,tmpvarid,aps)
407	if (ierr.ne.NF_NOERR) then
408	stop "Error: Failed reading aps"
409	endif
410	endif
411
412	! hybrid coordinate bps
413	ierr=NF_INQ_VARID(infid,"bps",tmpvarid)
414	if (ierr.ne.NF_NOERR) then
415	stop "Error: Failed to get bps ID"
416	else
417	allocate(bps(altlength))
418	ierr=NF_GET_VAR_REAL(infid,tmpvarid,bps)
419	if (ierr.ne.NF_NOERR) then
420	stop "Error: Failed reading bps"
421	endif
422	endif
423	endif !of if (.not.have_sigma)
424
425	! 1.3.3 ground geopotential phisinit
426
427	allocate(phisinit(lonlength,latlength),stat=ierr)
428	if (ierr.ne.0) then
429	write(,) "Failed allocation of phisinit(lonlength,latlength) !!!"
430	write(,) "lonlength=",lonlength," latlength=",latlength
431	endif
432	! look for 'phisinit' in current file
433	ierr=NF_INQ_VARID(infid,"phisinit",tmpvarid)
434	if (ierr.ne.NF_NOERR) then
435	write(,) "Warning: Failed to get phisinit ID from file ",trim(infile)
436	infile2="diagfi.nc"
437	write(,) " Trying file ",trim(infile2)
438	ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
439	if (ierr.ne.NF_NOERR) then
440	write(,) "Problem: Could not find/open that file"
441	infile2="phisinit.nc"
442	write(,) " Trying file ",trim(infile2)
443	ierr=NF_OPEN(infile2,NF_NOWRITE,infid2)
444	if (ierr.ne.NF_NOERR) then
445	write(,) "Error: Could not open that file either"
446	stop "Might as well stop here"
447	endif
448	endif
449
450	! Get ID for phisinit
451	ierr=NF_INQ_VARID(infid2,"phisinit",tmpvarid)
452	if (ierr.ne.NF_NOERR) then
453	stop "Error: Failed to get phisinit ID"
454	endif
455	! Get physinit
456	ierr=NF_GET_VAR_REAL(infid2,tmpvarid,phisinit)
457	if (ierr.ne.NF_NOERR) then
458	stop "Error: Failed reading phisinit"
459	endif
460	! Close file
461	write(,) 'OK, got phisinit'
462	ierr=NF_CLOSE(infid2)
463	else
464	ierr=NF_GET_VAR_REAL(infid,tmpvarid,phisinit)
465	if (ierr.ne.NF_NOERR) then
466	stop "Error: Failed reading phisinit"
467	endif
468	endif
469
470	!===============================================================================
471	! 1.4 Choose and build the new vertical coordinate
472	!===============================================================================
473
474	write(,) ""
475	write(,) "Which vertical coordinate should the output be in?"
476	ztype=0
477	do while ((ztype.lt.1).or.(ztype.gt.3))
478	write(,) "(1: pressure, 2: above areoid altitude 3: above local surface)"
479	read(,)ztype
480	enddo
481
482	text="dummy" ! dummy initialization
483	auto_vert_levels=.true. ! dummy initialization to get rid of compiler warning
484	do while ((trim(text).ne."yes").and.(trim(text).ne."no"))
485	write(,) ""
486	write(,) "Automatic generation of vertical levels distribution? (yes/no)"
487	write(,) "(yes: you only provide min, max and number of levels)"
488	write(,) "(no: you provide values for each level)"
489	read(*,'(a64)') text
490	if (trim(text).eq."yes") then
491	auto_vert_levels=.true.
492	else
493	auto_vert_levels=.false.
494	endif
495	enddo
496
497	if (auto_vert_levels) then
498	! ask for # of points and end values for pressure or above areoid cases
499	write(,) ""
500	if (ztype.le.2) then
501	write(,) "Enter min and max of vertical coordinate (Pa or m)"
502	write(,) " (in that order and on the same line)"
503	if (ztype.eq.1) then ! pressure coordinate
504	read(,) pmin,pmax
505	else ! above areoid coordinate
506	read(,) zamin,zamax
507	endif
508	endif
509
510	! Build corresponding vertical coordinates
511	if (ztype.eq.1) then ! pressure coordinate
512	write(,) "Number of levels along vertical coordinate?"
513	read(,) nblev
514	allocate(plevel(nblev))
515	if (nblev.eq.1) then ! in case only one level is asked for
516	plevel(nblev)=pmin
517	else
518	do i=1,nblev
519	! build exponentially spread layers
520	plevel(i)=exp(log(pmax)+(log(pmin)-log(pmax))* &
521	((real(i)-1.0)/(real(nblev)-1.0)))
522	enddo
523	endif
524	else if (ztype.eq.2) then ! above areoid heights
525	write(,) "Number of levels along vertical coordinate?"
526	read(,) nblev
527	allocate(zareoid(nblev))
528	if (nblev.eq.1) then ! in case only one level is asked for
529	zareoid(nblev)=zamin
530	else
531	do i=1,nblev
532	zareoid(i)=zamin+(real(i)-1.0)*((zamax-zamin)/(real(nblev)-1.0))
533	enddo
534	endif
535	else ! above local surface
536	! set nblev to # of levels from input data files
537	nblev=altlength
538	allocate(zsurface(nblev))
539	! build specific above local surface altitudes
540	call build_zs(nblev,have_sigma,sigma,aps,bps,zsurface)
541	endif
542	else ! auto_vert_levels=.false. ; user provides values
543	! ask for # of points along the vertical
544	write(,) ""
545	write(,) "Number of levels along vertical coordinate?"
546	read(,) nblev
547	if (ztype.eq.1) then ! pressure coordinate
548	allocate(plevel(nblev))
549	write(,) "Enter Pressure (Pa) of levels, ordered"
550	write(,) " from max (near-surface) to min (top of atmosphere),"
551	write(,) " (one value per line)"
552	do i=1,nblev
553	read(,) plevel(i)
554	enddo
555	else if (ztype.eq.2) then ! above areoid heights
556	allocate(zareoid(nblev))
557	write(,) "Enter altitude (m) above areoid of levels, ordered"
558	write(,) " from min to max (one value per line)"
559	do i=1,nblev
560	read(,) zareoid(i)
561	enddo
562	else ! above local surface
563	allocate(zsurface(nblev))
564	write(,) "Enter altitude (m) above surface of levels, ordered"
565	write(,) " from min to max (one value per line)"
566	do i=1,nblev
567	read(,) zsurface(i)
568	enddo
569	endif
570	endif ! of if (auto_vert_levels)
571
572	!===============================================================================
573	! 1.5 Get output file name
574	!===============================================================================
575	write(,) ""
576	!write(,) "Enter output file name"
577	!read(,) outfile
578	if (ztype.eq.1) then ! pressure coordinate
579	outfile=infile(1:len_trim(infile)-3)//"_P.nc"
580	else if (ztype.eq.2) then ! above areoid coordinate
581	outfile=infile(1:len_trim(infile)-3)//"_A.nc"
582	else ! above local surface
583	outfile=infile(1:len_trim(infile)-3)//"_S.nc"
584	endif
585	write(,) "Output file name is: "//trim(outfile)
586
587	!===============================================================================
588	! 2.1 Build/store GCM fields which will be used later
589	!===============================================================================
590
591	!===============================================================================
592	! 2.1.1 Surface pressure
593	!===============================================================================
594	ierr=NF_INQ_VARID(infid,"ps",tmpvarid)
595	if (ierr.ne.NF_NOERR) then
596	stop "Error: Failed to get ps ID"
597	else
598	allocate(ps(lonlength,latlength,timelength))
599	ierr=NF_GET_VAR_REAL(infid,tmpvarid,ps)
600	if (ierr.ne.NF_NOERR) then
601	stop "Error: Failed reading surface pressure"
602	endif
603	endif
604
605	!===============================================================================
606	! 2.1.2 Atmospheric pressure
607	!===============================================================================
608	allocate(press(lonlength,latlength,altlength,timelength))
609
610	if (have_sigma) then ! sigma coordinate
611	do itim=1,timelength
612	do ilev=1,altlength
613	do ilat=1,latlength
614	do ilon=1,lonlength
615	press(ilon,ilat,ilev,itim)=sigma(ilev)*ps(ilon,ilat,itim)
616	enddo
617	enddo
618	enddo
619	enddo
620	else ! hybrid coordinates
621	do itim=1,timelength
622	do ilev=1,altlength
623	do ilat=1,latlength
624	do ilon=1,lonlength
625	press(ilon,ilat,ilev,itim)=aps(ilev)+bps(ilev)*ps(ilon,ilat,itim)
626	enddo
627	enddo
628	enddo
629	enddo
630	endif
631
632	!===============================================================================
633	! 2.1.3 Atmospheric temperature
634	!===============================================================================
635	allocate(temp(lonlength,latlength,altlength,timelength))
636
637	ierr=NF_INQ_VARID(infid,"temp",tmpvarid)
638	if (ierr.ne.NF_NOERR) then
639	! stop "Error: Failed to get temp ID"
640	! try "t" for temperature
641	ierr=NF_INQ_VARID(infid,"t",tmpvarid)
642	if (ierr.ne.NF_NOERR) then
643	stop "Error: Failed to get t ID"
644	else
645	ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
646	if (ierr.ne.NF_NOERR) then
647	stop "Error: Failed reading atmospheric temperature"
648	endif
649	endif
650	else
651	ierr=NF_GET_VAR_REAL(infid,tmpvarid,temp)
652	if (ierr.ne.NF_NOERR) then
653	stop "Error: Failed reading atmospheric temperature"
654	endif
655	endif
656
657	!===============================================================================
658	! 2.1.4 Atmospheric density
659	!===============================================================================
660
661	ierr=NF_INQ_VARID(infid,"rho",tmpvarid)
662	if (ierr.ne.NF_NOERR) then
663	write(,) "Warning: Failed to get rho ID"
664	have_rho=.false.
665	else
666	have_rho=.true.
667	allocate(rho(lonlength,latlength,altlength,timelength))
668	ierr=NF_GET_VAR_REAL(infid,tmpvarid,rho)
669	if (ierr.ne.NF_NOERR) then
670	stop "Error: Failed reading atmospheric density"
671	endif
672	endif
673
674	!===============================================================================
675	! 2.2 Build GCM Above areoid (or above surface) altitudes of GCM nodes
676	!===============================================================================
677
678	if (have_rho) then
679	if (ztype.le.2) then ! above areoid altitudes (also needed for ztype=1)
680	allocate(za_gcm(lonlength,latlength,altlength,timelength))
681	call build_gcm_za(lonlength,latlength,altlength,timelength, &
682	phisinit,ps,press,temp,rho,za_gcm)
683	else ! above local surface altitudes
684	allocate(zs_gcm(lonlength,latlength,altlength,timelength))
685	call build_gcm_zs(lonlength,latlength,altlength,timelength, &
686	phisinit,ps,press,temp,rho,zs_gcm)
687	endif
688	else
689	write(,)"Warning: Using constant R to integrate hydrostatic equation"
690	if (ztype.le.2) then ! above areoid altitudes (also needed for ztype=1)
691	allocate(za_gcm(lonlength,latlength,altlength,timelength))
692	call crude_gcm_za(lonlength,latlength,altlength,timelength, &
693	phisinit,ps,press,temp,za_gcm)
694	else ! above local surface altitudes
695	allocate(zs_gcm(lonlength,latlength,altlength,timelength))
696	call crude_gcm_zs(lonlength,latlength,altlength,timelength, &
697	phisinit,ps,press,temp,zs_gcm)
698	endif
699	endif
700
701	!===============================================================================
702	! 3. Create output file and initialize definitions of variables and dimensions
703	!===============================================================================
704
705	!===============================================================================
706	! 3.1. Output file
707	!===============================================================================
708
709	! Create output file
710	ierr=NF_CREATE(outfile,NF_CLOBBER,outfid)
711	if (ierr.ne.NF_NOERR) then
712	write(,)"Error: could not create file ",outfile
713	stop
714	endif
715
716	!===============================================================================
717	! 3.2. Define dimensions
718	!===============================================================================
719	! longitude
720	ierr=NF_DEF_DIM(outfid,"longitude",lonlength,lon_dimid)
721	if (ierr.ne.NF_NOERR) then
722	stop "Error: Could not define longitude dimension"
723	endif
724
725	! latitude
726	ierr=NF_DEF_DIM(outfid,"latitude",latlength,lat_dimid)
727	if (ierr.ne.NF_NOERR) then
728	stop "Error: Could not define latitude dimension"
729	endif
730
731	! altitude
732	ierr=NF_DEF_DIM(outfid,"altitude",nblev,alt_dimid)
733	if (ierr.ne.NF_NOERR) then
734	stop "Error: Could not define altitude dimension"
735	endif
736
737	! time
738	ierr=NF_DEF_DIM(outfid,"Time",timelength,time_dimid)
739	if (ierr.ne.NF_NOERR) then
740	stop "Error: Could not define latitude dimension"
741	endif
742
743	! GCM layers (for sigma or aps and bps)
744	ierr=NF_DEF_DIM(outfid,"GCM_layers",altlength,gcm_layers_dimid)
745	if (ierr.ne.NF_NOERR) then
746	stop "Error: Could not define GCM_layers dimension"
747	endif
748
749
750	!===============================================================================
751	! 3.3. Define variables and their attributes
752	!===============================================================================
753
754	! 3.3.1 Define 1D variables
755
756	! longitude
757	datashape(1)=lon_dimid
758	ierr=NF_DEF_VAR(outfid,"longitude",NF_REAL,1,datashape(1),lon_varid)
759	if (ierr.ne.NF_NOERR) then
760	stop "Error: Could not define longitude variable"
761	endif
762
763	! longitude attributes
764	text='east longitude'
765	ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'long_name',len_trim(text),text)
766	if (ierr.ne.NF_NOERR) then
767	stop "Error: Problem writing long_name for longitude"
768	endif
769	text='degrees_east'
770	ierr=NF_PUT_ATT_TEXT(outfid,lon_varid,'units',len_trim(text),text)
771	if (ierr.ne.NF_NOERR) then
772	stop "Error: Problem writing units for longitude"
773	endif
774
775	! latitude
776	datashape(2)=lat_dimid
777	ierr=NF_DEF_VAR(outfid,"latitude",NF_REAL,1,datashape(2),lat_varid)
778	if (ierr.ne.NF_NOERR) then
779	stop "Error: Could not define latitude variable"
780	endif
781
782	! latitude attributes
783	text='north latitude'
784	ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'long_name',len_trim(text),text)
785	if (ierr.ne.NF_NOERR) then
786	stop "Error: Problem writing long_name for latitude"
787	endif
788	text='degrees_north'
789	ierr=NF_PUT_ATT_TEXT(outfid,lat_varid,'units',len_trim(text),text)
790	if (ierr.ne.NF_NOERR) then
791	stop "Error: Problem writing units for latitude"
792	endif
793
794	! altitude
795	datashape(3)=alt_dimid
796	ierr=NF_DEF_VAR(outfid,"altitude",NF_REAL,1,datashape(3),alt_varid)
797	if (ierr.ne.NF_NOERR) then
798	stop "Error: Could not define altitude variable"
799	endif
800
801	!altitude attributes
802	if (ztype.eq.1) then ! pressure vertical coordinate
803	text='Pressure levels'
804	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text)
805	if (ierr.ne.NF_NOERR) then
806	stop "Error: Problem writing long_name for altitude"
807	endif
808	text='Pa'
809	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text)
810	if (ierr.ne.NF_NOERR) then
811	stop "Error: Problem writing units for altitude"
812	endif
813	text='down'
814	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text)
815	if (ierr.ne.NF_NOERR) then
816	stop "Error: Problem writing positive for altitude"
817	endif
818	else if (ztype.eq.2) then ! above areoid vertical coordinate
819	text='Altitude above areoid'
820	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text)
821	if (ierr.ne.NF_NOERR) then
822	stop "Error: Problem writing long_name for altitude"
823	endif
824	text='m'
825	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text)
826	if (ierr.ne.NF_NOERR) then
827	stop "Error: Problem writing units for altitude"
828	endif
829	text='up'
830	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text)
831	if (ierr.ne.NF_NOERR) then
832	stop "Error: Problem writing positive for altitude"
833	endif
834	else ! above surface vertical coordinate
835	text='Altitude above local surface'
836	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'long_name',len_trim(text),text)
837	if (ierr.ne.NF_NOERR) then
838	stop "Error: Problem writing long_name for altitude"
839	endif
840	text='m'
841	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'units',len_trim(text),text)
842	if (ierr.ne.NF_NOERR) then
843	stop "Error: Problem writing units for altitude"
844	endif
845	text='up'
846	ierr=NF_PUT_ATT_TEXT(outfid,alt_varid,'positive',len_trim(text),text)
847	if (ierr.ne.NF_NOERR) then
848	stop "Error: Problem writing positive for altitude"
849	endif
850	endif ! of if (have_sigma)
851
852	! GCM_layers
853	!ierr=NF_DEF_VAR(outfid,"gcm_layers",NF_REAL,1,,sigma_varid)
854
855	! sigma levels or hybrid coordinates
856	if (have_sigma) then
857	ierr=NF_DEF_VAR(outfid,"sigma",NF_REAL,1,gcm_layers_dimid,sigma_varid)
858	if (ierr.ne.NF_NOERR) then
859	stop "Error: Could not define sigma variable"
860	endif
861	else ! hybrid coordinates
862	ierr=NF_DEF_VAR(outfid,"aps",NF_REAL,1,gcm_layers_dimid,aps_varid)
863	if (ierr.ne.NF_NOERR) then
864	stop "Error: Could not define aps variable"
865	endif
866	ierr=NF_DEF_VAR(outfid,"bps",NF_REAL,1,gcm_layers_dimid,bps_varid)
867	if (ierr.ne.NF_NOERR) then
868	stop "Error: Could not define bps variable"
869	endif
870	endif
871
872	! sigma levels (or hybrid coordinates) attributes
873	if (have_sigma) then
874	text="sigma levels"
875	ierr=NF_PUT_ATT_TEXT(outfid,sigma_varid,'long_name',len_trim(text),text)
876	if (ierr.ne.NF_NOERR) then
877	stop "Error: Problem writing long_name for sigma"
878	endif
879	else ! hybrid coordinates
880	text="hybrid pressure at midlayers"
881	ierr=NF_PUT_ATT_TEXT(outfid,aps_varid,'long_name',len_trim(text),text)
882	if (ierr.ne.NF_NOERR) then
883	stop "Error: Problem writing long_name for aps"
884	endif
885	text="hybrid sigma at midlayers"
886	ierr=NF_PUT_ATT_TEXT(outfid,bps_varid,'long_name',len_trim(text),text)
887	if (ierr.ne.NF_NOERR) then
888	stop "Error: Problem writing long_name for bps"
889	endif
890	endif ! of if (have_sigma)
891
892	! time
893	datashape(4)=time_dimid
894	ierr=NF_DEF_VAR(outfid,"Time",NF_REAL,1,datashape(4),time_varid)
895	if (ierr.ne.NF_NOERR) then
896	stop "Error: Could not define Time variable"
897	endif
898
899	! time attributes
900	text='Time'
901	ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'long_name',len_trim(text),text)
902	if (ierr.ne.NF_NOERR) then
903	stop "Error: Problem writing long_name for Time"
904	endif
905	text='days since 0000-01-1 00:00:00'
906	ierr=NF_PUT_ATT_TEXT(outfid,time_varid,'units',len_trim(text),text)
907	if (ierr.ne.NF_NOERR) then
908	stop "Error: Problem writing units for Time"
909	endif
910
911	! 3.3.2 Define 3D variables (ie: surface+time variables)
912
913	! Surface pressure
914	surfdatashape(1)=lon_dimid
915	surfdatashape(2)=lat_dimid
916	surfdatashape(3)=time_dimid
917	ierr=NF_DEF_VAR(outfid,"ps",NF_REAL,3,surfdatashape,ps_varid)
918	if (ierr.ne.NF_NOERR) then
919	stop "Error: Could not define ps variable"
920	endif
921
922	! Surface pressure attributes
923	text='Surface pressure'
924	ierr=NF_PUT_ATT_TEXT(outfid,ps_varid,'long_name',len_trim(text),text)
925	if (ierr.ne.NF_NOERR) then
926	stop "Error: Problem writing long_name for surface pressure"
927	endif
928	text='Pa'
929	ierr=NF_PUT_ATT_TEXT(outfid,ps_varid,'units',len_trim(text),text)
930	if (ierr.ne.NF_NOERR) then
931	stop "Error: Problem writing units for surface pressure"
932	endif
933
934	! 3.3.3 Define 4D variables
935
936	! add pressure or zareoid
937	if (ztype.eq.1) then ! pressure vertical coordinate
938	! zareoid dataset
939	ierr=NF_DEF_VAR(outfid,"zareoid",NF_REAL,4,datashape,za_varid)
940	if (ierr.ne.NF_NOERR) then
941	stop "Error: Could not define zareoid variable"
942	endif
943	! zareoid attributes
944	text='altitude above areoid'
945	ierr=NF_PUT_ATT_TEXT(outfid,za_varid,'long_name',len_trim(text),text)
946	if (ierr.ne.NF_NOERR) then
947	stop "Error: Problem writing long_name for zareoid"
948	endif
949	text='m'
950	ierr=NF_PUT_ATT_TEXT(outfid,za_varid,'units',len_trim(text),text)
951	if (ierr.ne.NF_NOERR) then
952	stop "Error: Problem writing units for zareoid"
953	endif
954	! zareoid missing value
955	ierr=NF_PUT_ATT_REAL(outfid,za_varid,'missing_value',NF_REAL,1,miss_val)
956	if (ierr.ne.NF_NOERR) then
957	stop "Error: Problem writing missing_value for zareoid"
958	endif
959	else ! above areoid or above local surface vertical coordinate
960	! pressure dataset
961	ierr=NF_DEF_VAR(outfid,"pressure",NF_REAL,4,datashape,p_varid)
962	if (ierr.ne.NF_NOERR) then
963	stop "Error: Could not define pressure variable"
964	endif
965	! pressure attributes
966	text='Atmospheric pressure'
967	ierr=NF_PUT_ATT_TEXT(outfid,p_varid,'long_name',len_trim(text),text)
968	if (ierr.ne.NF_NOERR) then
969	stop "Error: Problem writing long_name for pressure"
970	endif
971	text='Pa'
972	ierr=NF_PUT_ATT_TEXT(outfid,p_varid,'units',len_trim(text),text)
973	if (ierr.ne.NF_NOERR) then
974	stop "Error: Problem writing units for pressure"
975	endif
976	! pressure missing value
977	ierr=NF_PUT_ATT_REAL(outfid,p_varid,'missing_value',NF_REAL,1,miss_val)
978	if (ierr.ne.NF_NOERR) then
979	stop "Error: Problem writing missing_value for pressure"
980	endif
981	endif
982
983	! add zs_gcm
984	if (ztype.eq.3) then
985	endif
986
987	! variables requested by user
988	allocate(var_id(nbvar))
989	do i=1,nbvar
990	write(,) ""
991	write(,) "Creating ",trim(var(i))
992	! define the variable
993	ierr=NF_DEF_VAR(outfid,var(i),NF_REAL,4,datashape,var_id(i))
994	if (ierr.ne.NF_NOERR) then
995	write(,) 'Error, could not define variable ',trim(var(i))
996	stop ""
997	endif
998
999	! Get the (input file) ID for the variable and
1000	! the # of attributes associated to that variable
1001	ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
1002	if (ierr.ne.NF_NOERR) then
1003	write(,) 'Error, failed to get ID for variable ',trim(var(i))
1004	stop ""
1005	endif
1006	ierr=NF_INQ_VARNATTS(infid,tmpvarid,nbattr)
1007	if (ierr.ne.NF_NOERR) then
1008	write(,) 'Error, could not get number of attributes for variable ',&
1009	trim(var(i))
1010	stop ""
1011	endif
1012	! inititialize j == number of attributes written to output
1013	j=0
1014
1015	! look for a "long_name" attribute
1016	text=' '
1017	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'long_name',text)
1018	if (ierr.ne.NF_NOERR) then ! no long_name attribute
1019	! try to find an equivalent 'title' attribute
1020	text=' '
1021	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'title',text)
1022	if (ierr.eq.NF_NOERR) then ! found 'title' attribute
1023	write(,) "Found title ",trim(text)
1024	j=j+1
1025	! write it as a 'long_name' attribute
1026	ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
1027	if (ierr.ne.NF_NOERR) then
1028	write(,) "Error failed to copy title attribute:",trim(text)
1029	stop ""
1030	endif
1031	else ! no 'title' attribute
1032	! try to find a "Physics_diagnostic" attribute (UK GCM outputs)
1033	text=' '
1034	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'Physics_diagnostic',text)
1035	if (ierr.eq.NF_NOERR) then ! found 'Physics_diagnostic' attribute
1036	write(,) "Found Physics_diagnostic ",trim(text)
1037	j=j+1
1038	! write it as a 'long_name' attribute
1039	ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
1040	if (ierr.ne.NF_NOERR) then
1041	write(,) "Error failed to copy Physics_diagnostic attribute:",trim(text)
1042	stop
1043	endif
1044	endif
1045	endif
1046	else ! found long_name; write it to outfile
1047	write(,) "Found long_name ",trim(text)
1048	j=j+1
1049	ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'long_name',len_trim(text),text)
1050	if (ierr.ne.NF_NOERR) then
1051	write(,) "Error failed to copy long_name attribute:",trim(text)
1052	stop""
1053	endif
1054	endif
1055
1056	! look for a "units" attribute
1057	text=' '
1058	ierr=NF_GET_ATT_TEXT(infid,tmpvarid,'units',text)
1059	if (ierr.eq.NF_NOERR) then ! found 'units' attribute
1060	write(,) "Found units ",trim(text)
1061	j=j+1
1062	! write it to output
1063	ierr=NF_PUT_ATT_TEXT(outfid,var_id(i),'units',len_trim(text),text)
1064	if (ierr.ne.NF_NOERR) then
1065	write(,) "Error failed to copy units attribute:",trim(text)
1066	stop""
1067	endif
1068	endif
1069
1070	! look for a "missing_value" attribute
1071	ierr=NF_GET_ATT_REAL(infid,tmpvarid,"missing_value",miss_val)
1072	if (ierr.eq.NF_NOERR) then ! found 'missing_value' attribute
1073	write(,) "Found missing_value ",miss_val
1074	j=j+1
1075	else ! no 'missing_value' attribute, set miss_val to default
1076	miss_val=miss_val_def
1077	endif
1078
1079	! write the missing_value attribute to output
1080	ierr=NF_PUT_ATT_REAL(outfid,var_id(i),'missing_value',NF_REAL,1,miss_val)
1081	if (ierr.ne.NF_NOERR) then
1082	stop "Error, failed to write missing_value attribute"
1083	endif
1084
1085	! warn if some attributes were missed
1086	if (j.ne.nbattr) then
1087	write(,)'Warning, it seems some attributes of variable ',trim(var(i))
1088	write(,)"were not transfered to the new file"
1089	write(,)'nbattr:',nbattr,' j:',j
1090	endif
1091
1092	enddo ! of do=1,nbvar
1093
1094
1095	!===============================================================================
1096	! 3.4. Write dimensions (and 1D varaiables)
1097	!===============================================================================
1098	! Switch out of NetCDF define mode
1099	ierr=NF_ENDDEF(outfid)
1100	if (ierr.ne.NF_NOERR) then
1101	stop "Error: Could not switch out of define mode"
1102	endif
1103
1104	! Write longitude
1105	ierr=NF_PUT_VAR_REAL(outfid,lon_varid,lon)
1106	if (ierr.ne.NF_NOERR) then
1107	stop "Error: Could not write longitude data to output file"
1108	endif
1109
1110	! Write latitude
1111	ierr=NF_PUT_VAR_REAL(outfid,lat_varid,lat)
1112	if (ierr.ne.NF_NOERR) then
1113	stop "Error: Could not write latitude data to output file"
1114	endif
1115
1116	! Write altitude
1117	if (ztype.eq.1) then ! pressure vertical coordinate
1118	ierr=NF_PUT_VAR_REAL(outfid,alt_varid,plevel)
1119	if (ierr.ne.NF_NOERR) then
1120	stop "Error: Could not write altitude data to output file"
1121	endif
1122	else if (ztype.eq.2) then ! above areoid altitude
1123	ierr=NF_PUT_VAR_REAL(outfid,alt_varid,zareoid)
1124	if (ierr.ne.NF_NOERR) then
1125	stop "Error: Could not write altitude data to output file"
1126	endif
1127	else ! above local surface
1128	ierr=NF_PUT_VAR_REAL(outfid,alt_varid,zsurface)
1129	if (ierr.ne.NF_NOERR) then
1130	stop "Error: Could not write altitude data to output file"
1131	endif
1132	endif
1133
1134	! Write sigma levels (or hybrid coordinates)
1135	if (have_sigma) then
1136	ierr=NF_PUT_VAR_REAL(outfid,sigma_varid,sigma)
1137	if (ierr.ne.NF_NOERR) then
1138	stop "Error: Could not write sigma data to output file"
1139	endif
1140	else ! hybrid coordinates
1141	ierr=NF_PUT_VAR_REAL(outfid,aps_varid,aps)
1142	if (ierr.ne.NF_NOERR) then
1143	stop "Error: Could not write aps data to output file"
1144	endif
1145	ierr=NF_PUT_VAR_REAL(outfid,bps_varid,bps)
1146	if (ierr.ne.NF_NOERR) then
1147	stop "Error: Could not write bps data to output file"
1148	endif
1149	endif
1150
1151	! Write time
1152	ierr=NF_PUT_VAR_REAL(outfid,time_varid,time)
1153	if (ierr.ne.NF_NOERR) then
1154	stop "Error: Could not write Time data to output file"
1155	endif
1156
1157	!===============================================================================
1158	! 3.5 Write 3D variables
1159	!===============================================================================
1160
1161	! Write surface pressure
1162
1163	ierr=NF_PUT_VAR_REAL(outfid,ps_varid,ps)
1164	if (ierr.ne.NF_NOERR) then
1165	stop "Error: Could not write ps data to output file"
1166	endif
1167
1168	!===============================================================================
1169	! 4. Interpolate and write 4D variables
1170	!===============================================================================
1171
1172	! 4.0 Allocations
1173	!indata() to store input (on GCM grid) data
1174	allocate(indata(lonlength,latlength,altlength,timelength))
1175	! outdata() to store output (on new vertical grid) data
1176	allocate(outdata(lonlength,latlength,nblev,timelength))
1177
1178	! 4.1 If output is in pressure coordinate
1179	if (ztype.eq.1) then
1180	do i=1,nbvar ! loop on 4D variable to process
1181	! identify and read a dataset
1182	ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
1183	if (ierr.ne.NF_NOERR) then
1184	write(,) 'Error, failed to get ID for variable ',var(i)
1185	stop
1186	endif
1187	ierr=NF_GET_VAR_REAL(infid,tmpvarid,indata)
1188	if (ierr.ne.NF_NOERR) then
1189	write(,) 'Error, failed to load variable ',var(i)
1190	stop
1191	endif
1192
1193	! interpolate dataset onto new grid
1194	call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
1195	miss_val,ps,press,indata,plevel,outdata)
1196
1197	! write new dataset to output file
1198	ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata)
1199	if (ierr.ne.NF_NOERR) then
1200	write(,)'Error, Failed to write variable ',var(i)
1201	stop
1202	endif
1203
1204	enddo ! of do i=1,nbvar
1205
1206	! interpolate zareoid onto new grid
1207	call p_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
1208	miss_val,ps,press,za_gcm,plevel,outdata)
1209	! write result to output file
1210	ierr=NF_PUT_VAR_REAL(outfid,za_varid,outdata)
1211	if (ierr.ne.NF_NOERR) then
1212	stop "Error, Failed to write zareoid to output file"
1213	endif
1214	endif ! of if (ztype.eq.1)
1215
1216	! 4.2 If output is in above areoid altitude
1217	if (ztype.eq.2) then
1218	do i=1,nbvar ! loop on 4D variable to process
1219	! identify and read a dataset
1220	ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
1221	if (ierr.ne.NF_NOERR) then
1222	write(,) 'Error, failed to get ID for variable ',var(i)
1223	stop
1224	endif
1225	ierr=NF_GET_VAR_REAL(infid,tmpvarid,indata)
1226	if (ierr.ne.NF_NOERR) then
1227	write(,) 'Error, failed to load variable ',var(i)
1228	stop
1229	endif
1230
1231	! interpolate dataset onto new grid
1232	! check if variable is "rho" (to set flag for interpolation below)
1233	if (var(i).eq.'rho') then
1234	j=1
1235	else
1236	j=0
1237	endif
1238
1239	call z_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
1240	miss_val,za_gcm,indata,j,zareoid,outdata)
1241
1242	! write new dataset to output file
1243	ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata)
1244	if (ierr.ne.NF_NOERR) then
1245	write(,)'Error, Failed to write variable ',var(i)
1246	stop
1247	endif
1248
1249	enddo ! of do i=1,nbvar
1250
1251	! interpolate pressure onto new grid
1252	write(,) ""
1253	write(,) "Processing pressure"
1254	call z_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
1255	miss_val,za_gcm,press,1,zareoid,outdata)
1256
1257	! write result to output file
1258	ierr=NF_PUT_VAR_REAL(outfid,p_varid,outdata)
1259	if (ierr.ne.NF_NOERR) then
1260	stop "Error, Failed to write pressure to output file"
1261	endif
1262	endif ! of if (ztype.eq.2)
1263
1264	! 4.3 If output is in above local surface altitude
1265	if (ztype.eq.3) then
1266	do i=1,nbvar ! loop on 4D variable to process
1267	write(,) " "
1268	write(,) "Processing "//trim(var(i))
1269	! identify and read a dataset
1270	ierr=NF_INQ_VARID(infid,var(i),tmpvarid)
1271	if (ierr.ne.NF_NOERR) then
1272	write(,) 'Error, failed to get ID for variable ',var(i)
1273	stop
1274	endif
1275	ierr=NF_GET_VAR_REAL(infid,tmpvarid,indata)
1276	if (ierr.ne.NF_NOERR) then
1277	write(,) 'Error, failed to load variable ',var(i)
1278	stop
1279	endif
1280
1281	! interpolate dataset onto new grid
1282	! check if variable is "rho" (to set flag for interpolation below)
1283	if (var(i).eq.'rho') then
1284	j=1
1285	else
1286	j=0
1287	endif
1288
1289	call zs_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
1290	miss_val,zs_gcm,indata,phisinit,press,temp,rho,&
1291	have_rho,j,zsurface,outdata)
1292
1293	! write new dataset to output file
1294	ierr=NF_PUT_VAR_REAL(outfid,var_id(i),outdata)
1295	if (ierr.ne.NF_NOERR) then
1296	write(,)'Error, Failed to write variable ',var(i)
1297	stop
1298	endif
1299	enddo ! of do i=1,nbvar
1300
1301	! interpolate pressure onto new grid
1302	write(,) ""
1303	write(,) "Processing pressure"
1304	call zs_coord_interp(lonlength,latlength,altlength,timelength,nblev, &
1305	miss_val,zs_gcm,press,phisinit,press,temp,rho,&
1306	have_rho,1,zsurface,outdata)
1307
1308	! write result to output file
1309	ierr=NF_PUT_VAR_REAL(outfid,p_varid,outdata)
1310	if (ierr.ne.NF_NOERR) then
1311	stop "Error, Failed to write pressure to output file"
1312	endif
1313
1314	endif ! of if (ztype.eq.3)
1315
1316	! 4.4 Close output file
1317	ierr=NF_CLOSE(outfid)
1318	if (ierr.ne.NF_NOERR) then
1319	write(,) 'Error, failed to close output file ',outfile
1320	endif
1321
1322	end program
1323
1324
1325	!===============================================================================
1326
1327	subroutine build_gcm_zs(lonlength,latlength,altlength,timelength, &
1328	phis,ps,press,temp,rho,zs_gcm)
1329	!==============================================================================
1330	! Purpose: Integrate hydrostatic equation in order to build the "above local
1331	! surface altitudes" corresponding to GCM atmospheric levels
1332	!==============================================================================
1333	implicit none
1334	!==============================================================================
1335	! Arguments:
1336	!==============================================================================
1337	integer,intent(in) :: lonlength ! # of points along longitude
1338	integer,intent(in) :: latlength ! # of points along latitude
1339	integer,intent(in) :: altlength ! # of points along altitude
1340	integer,intent(in) :: timelength ! # of stored time steps
1341	real,dimension(lonlength,latlength),intent(in) :: phis
1342	! phis(:,:) is the ground geopotential
1343	real,dimension(lonlength,latlength,timelength),intent(in) :: ps
1344	! ps(:,:) is the surface pressure
1345	real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press
1346	! press(:,:,:,:) is atmospheric pressure
1347	real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp
1348	! temp(:,:,:,:) is atmospheric temperature
1349	real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: rho
1350	! rho(:,:,:,:) is atmospheric density
1351
1352	real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: zs_gcm
1353	! zs_gcm(:,:,:,:) are the above local surface altitudes of GCM levels
1354
1355	!===============================================================================
1356	! Local variables:
1357	!===============================================================================
1358	real,dimension(:),allocatable :: sigma ! GCM sigma levels
1359	real,dimension(:,:,:,:),allocatable :: R ! molecular gas constant
1360	!real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude
1361	real :: Rmean ! mean value of R for a given level
1362	real :: Tmean ! "mean" value of temperature for a given level
1363	integer iloop,jloop,kloop,tloop
1364
1365	! Parameters needed to integrate hydrostatic equation:
1366	real,parameter :: g0=3.7257964
1367	!g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
1368	real,parameter :: a0=3396.E3
1369	!a0: 'mean' Mars radius=3396.km
1370	real :: gz
1371	! gz: gravitational acceleration at a given (zareoid) altitude
1372
1373	!===============================================================================
1374	! 1. Various initialisations
1375	!===============================================================================
1376	allocate(sigma(altlength))
1377	allocate(R(lonlength,latlength,altlength,timelength))
1378	!allocate(zlocal(lonlength,latlength,altlength,timelength))
1379
1380	!==============================================================================
1381	! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point
1382	! --later used to integrate the hydrostatic equation--
1383	!==============================================================================
1384	do tloop=1,timelength
1385	do kloop=1,altlength
1386	do jloop=1,latlength
1387	do iloop=1,lonlength
1388	R(iloop,jloop,kloop,tloop)=press(iloop,jloop,kloop,tloop)/ &
1389	(rho(iloop,jloop,kloop,tloop)* &
1390	temp(iloop,jloop,kloop,tloop))
1391	enddo
1392	enddo
1393	enddo
1394	enddo
1395
1396	!===============================================================================
1397	! 3. Integrate hydrostatic equation to compute zlocal and za_gcm
1398	!===============================================================================
1399	do tloop=1,timelength
1400	do jloop=1,latlength
1401	do iloop=1,lonlength
1402	! handle case of first altitude level
1403	sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop)
1404	zs_gcm(iloop,jloop,1,tloop)=-log(sigma(1))R(iloop,jloop,1,tloop) &
1405	temp(iloop,jloop,1,tloop)/g0
1406	! za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ &
1407	! phis(iloop,jloop)/g0
1408	do kloop=2,altlength
1409	! compute sigma level of layer
1410	sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
1411
1412	! compute "mean" temperature of the layer
1413	if(temp(iloop,jloop,kloop,tloop).eq. &
1414	temp(iloop,jloop,kloop-1,tloop)) then
1415	Tmean=temp(iloop,jloop,kloop,tloop)
1416	else
1417	Tmean=(temp(iloop,jloop,kloop,tloop)- &
1418	temp(iloop,jloop,kloop-1,tloop))/ &
1419	log(temp(iloop,jloop,kloop,tloop)/ &
1420	temp(iloop,jloop,kloop-1,tloop))
1421	endif
1422
1423	! compute mean value of R of the layer
1424	Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop))
1425
1426	! compute gravitational acceleration (at altitude zaeroid(kloop-1))
1427	! NB: zareoid=zsurface+phis/g0
1428	gz=g0(a0*2)/ &
1429	(a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2
1430
1431	! compute zs_gcm(iloop,jloop,kloop,tloop)
1432	zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- &
1433	log(sigma(kloop)/sigma(kloop-1))RmeanTmean/gz
1434
1435	! compute za_gcm(kloop)
1436	! za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
1437	! phis(iloop,jloop)/g0
1438	enddo ! kloop
1439	enddo ! iloop
1440	enddo ! jloop
1441	enddo ! tloop
1442
1443	! Cleanup
1444	deallocate(sigma)
1445	deallocate(R)
1446	!deallocate(zlocal)
1447
1448	end subroutine build_gcm_zs
1449
1450	!===============================================================================
1451
1452	!#include"build_gcm_alt.F90"
1453	subroutine build_gcm_za(lonlength,latlength,altlength,timelength, &
1454	phis,ps,press,temp,rho,za_gcm)
1455	!==============================================================================
1456	! Purpose: Integrate hydrostatic equation in order to build the "above areoid
1457	! altitudes" corresponding to GCM atmospheric levels
1458	!==============================================================================
1459	implicit none
1460	!==============================================================================
1461	! Arguments:
1462	!==============================================================================
1463	integer,intent(in) :: lonlength ! # of points along longitude
1464	integer,intent(in) :: latlength ! # of points along latitude
1465	integer,intent(in) :: altlength ! # of points along altitude
1466	integer,intent(in) :: timelength ! # of stored time steps
1467	real,dimension(lonlength,latlength),intent(in) :: phis
1468	! phis(:,:) is the ground geopotential
1469	real,dimension(lonlength,latlength,timelength),intent(in) :: ps
1470	! ps(:,:) is the surface pressure
1471	real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press
1472	! press(:,:,:,:) is atmospheric pressure
1473	real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp
1474	! temp(:,:,:,:) is atmospheric temperature
1475	real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: rho
1476	! rho(:,:,:,:) is atmospheric density
1477
1478	real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: za_gcm
1479	! za_gcm(:,:,:,:) are the above aroid heights of GCM levels
1480
1481	!===============================================================================
1482	! Local variables:
1483	!===============================================================================
1484	real,dimension(:),allocatable :: sigma ! GCM sigma levels
1485	real,dimension(:,:,:,:),allocatable :: R ! molecular gas constant
1486	real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude
1487	real :: Rmean ! mean value of R for a given level
1488	real :: Tmean ! "mean" value of temperature for a given level
1489	integer iloop,jloop,kloop,tloop
1490
1491	! Parameters needed to integrate hydrostatic equation:
1492	real,parameter :: g0=3.7257964
1493	!g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
1494	real,parameter :: a0=3396.E3
1495	!a0: 'mean' Mars radius=3396.km
1496	real :: gz
1497	! gz: gravitational acceleration at a given (zareoid) altitude
1498
1499	!===============================================================================
1500	! 1. Various initialisations
1501	!===============================================================================
1502	allocate(sigma(altlength))
1503	allocate(R(lonlength,latlength,altlength,timelength))
1504	allocate(zlocal(lonlength,latlength,altlength,timelength))
1505
1506	!==============================================================================
1507	! 2. Compute Molecular Gas Constant (J kg-1 K-1) at every grid point
1508	! --later used to integrate the hydrostatic equation--
1509	!==============================================================================
1510	do tloop=1,timelength
1511	do kloop=1,altlength
1512	do jloop=1,latlength
1513	do iloop=1,lonlength
1514	R(iloop,jloop,kloop,tloop)=press(iloop,jloop,kloop,tloop)/ &
1515	(rho(iloop,jloop,kloop,tloop)* &
1516	temp(iloop,jloop,kloop,tloop))
1517	enddo
1518	enddo
1519	enddo
1520	enddo
1521
1522	!===============================================================================
1523	! 3. Integrate hydrostatic equation to compute zlocal and za_gcm
1524	!===============================================================================
1525	do tloop=1,timelength
1526	do jloop=1,latlength
1527	do iloop=1,lonlength
1528	! handle case of first altitude level
1529	sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop)
1530	zlocal(iloop,jloop,1,tloop)=-log(sigma(1))R(iloop,jloop,1,tloop) &
1531	temp(iloop,jloop,1,tloop)/g0
1532	za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ &
1533	phis(iloop,jloop)/g0
1534	do kloop=2,altlength
1535	! compute sigma level of layer
1536	sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
1537
1538	! compute "mean" temperature of the layer
1539	if(temp(iloop,jloop,kloop,tloop).eq. &
1540	temp(iloop,jloop,kloop-1,tloop)) then
1541	Tmean=temp(iloop,jloop,kloop,tloop)
1542	else
1543	Tmean=(temp(iloop,jloop,kloop,tloop)- &
1544	temp(iloop,jloop,kloop-1,tloop))/ &
1545	log(temp(iloop,jloop,kloop,tloop)/ &
1546	temp(iloop,jloop,kloop-1,tloop))
1547	endif
1548
1549	! compute mean value of R of the layer
1550	Rmean=0.5*(R(iloop,jloop,kloop,tloop)+R(iloop,jloop,kloop-1,tloop))
1551
1552	! compute gravitational acceleration (at altitude zaeroid(kloop-1))
1553	gz=g0(a02)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))*2
1554
1555	! compute zlocal(kloop)
1556	zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
1557	log(sigma(kloop)/sigma(kloop-1))RmeanTmean/gz
1558
1559	! compute za_gcm(kloop)
1560	za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
1561	phis(iloop,jloop)/g0
1562	enddo ! kloop
1563	enddo ! iloop
1564	enddo ! jloop
1565	enddo ! tloop
1566
1567	! Cleanup
1568	deallocate(sigma)
1569	deallocate(R)
1570	deallocate(zlocal)
1571
1572	end subroutine build_gcm_za
1573
1574	!===============================================================================
1575
1576	subroutine build_zs(altlength,have_sigma,sigma,aps,bps,zsurface)
1577	!==============================================================================
1578	! Build generic above surface altitudes, using either sigma levels
1579	! or hybrid vertical coordinates.
1580	! In order to do so, we need to use scale heights that vary with altitude.
1581	! The scale heights distribution is then set to vary from a min. to a max.
1582	! following a tanh() law over an imposed transition range (see below).
1583	! Similarily, a reference mean surface pressure (610 Pa) is used to
1584	! compute generic above surface altitudes.
1585	!==============================================================================
1586	implicit none
1587	!==============================================================================
1588	! Arguments:
1589	!==============================================================================
1590	integer,intent(in) :: altlength ! # of points along altitude
1591	logical,intent(in) :: have_sigma ! true if sigma(:) are known
1592	! have_sigma is false if aps() and bps(:) are given instead
1593	real,dimension(altlength),intent(in) :: sigma ! sigma levels
1594	real,dimension(altlength),intent(in) :: aps ! hybrid pressure levels
1595	real,dimension(altlength),intent(in) :: bps ! hybrid sigma levels
1596	real,dimension(altlength),intent(out) :: zsurface ! altitudes (m)
1597
1598	!===============================================================================
1599	! Local variables:
1600	!===============================================================================
1601	real,dimension(:),allocatable :: H ! scale heights
1602	real,parameter :: H_low=9650 ! scale height at low altitudes
1603	real,parameter :: H_high=15000 ! scale height at high altitudes
1604	real,parameter :: trans_window=10 ! # of levels over which H(:) goes
1605	! from H_low to H_high
1606	real,parameter :: lev_trans=32+trans_window/2
1607	! lev_trans: level at which H(lev_trans)=(H_low+H_high)/2
1608	! N.B.: keep lev_trans and lev_trans as reals to avoid truncation issues
1609
1610	real,parameter :: P_ref=610 ! reference surface pressure used to build zsurface
1611	integer i
1612
1613	! 1. Build scale heights
1614	allocate(H(altlength))
1615	do i=1,altlength
1616	H(i)=H_low+(H_high-H_low)0.5(1.0+tanh(6.*(i-lev_trans)/trans_window))
1617	enddo
1618
1619	! 2. Compute zsurface(:)
1620	if (have_sigma) then ! use sigma levels
1621	do i=1,altlength
1622	zsurface(i)=-H(i)log(sigma(i)P_ref)
1623	enddo
1624	else ! use hybrid coordinates
1625	do i=1,altlength
1626	zsurface(i)=-H(i)*log((aps(i)/P_ref)+bps(i))
1627	enddo
1628	endif
1629
1630	! Cleanup
1631	deallocate(H)
1632
1633	end subroutine build_zs
1634
1635	!===============================================================================
1636
1637	subroutine crude_gcm_zs(lonlength,latlength,altlength,timelength, &
1638	phis,ps,press,temp,zs_gcm)
1639	!==============================================================================
1640	! Purpose: Integrate hydrostatic equation in order to build the "above local
1641	! surface altitudes" corresponding to GCM atmospheric levels
1642	!==============================================================================
1643	implicit none
1644	!==============================================================================
1645	! Arguments:
1646	!==============================================================================
1647	integer,intent(in) :: lonlength ! # of points along longitude
1648	integer,intent(in) :: latlength ! # of points along latitude
1649	integer,intent(in) :: altlength ! # of points along altitude
1650	integer,intent(in) :: timelength ! # of stored time steps
1651	real,dimension(lonlength,latlength),intent(in) :: phis
1652	! phis(:,:) is the ground geopotential
1653	real,dimension(lonlength,latlength,timelength),intent(in) :: ps
1654	! ps(:,:) is the surface pressure
1655	real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press
1656	! press(:,:,:,:) is atmospheric pressure
1657	real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp
1658	! temp(:,:,:,:) is atmospheric temperature
1659
1660	real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: zs_gcm
1661	! zs_gcm(:,:,:,:) are the above local surface altitudes of GCM levels
1662
1663	!===============================================================================
1664	! Local variables:
1665	!===============================================================================
1666	real,dimension(:),allocatable :: sigma ! GCM sigma levels
1667	real,parameter :: R=191 ! molecular gas constant
1668	!real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude
1669	real :: Tmean ! "mean" value of temperature for a given level
1670	integer iloop,jloop,kloop,tloop
1671
1672	! Parameters needed to integrate hydrostatic equation:
1673	real,parameter :: g0=3.7257964
1674	!g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
1675	real,parameter :: a0=3396.E3
1676	!a0: 'mean' Mars radius=3396.km
1677	real :: gz
1678	! gz: gravitational acceleration at a given (zareoid) altitude
1679
1680	!===============================================================================
1681	! 1. Various initialisations
1682	!===============================================================================
1683	allocate(sigma(altlength))
1684	!allocate(zlocal(lonlength,latlength,altlength,timelength))
1685
1686	!===============================================================================
1687	! 2. Integrate hydrostatic equation to compute zlocal and za_gcm
1688	!===============================================================================
1689	do tloop=1,timelength
1690	do jloop=1,latlength
1691	do iloop=1,lonlength
1692	! handle case of first altitude level
1693	sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop)
1694	zs_gcm(iloop,jloop,1,tloop)=-log(sigma(1))* &
1695	R*temp(iloop,jloop,1,tloop)/g0
1696	! za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ &
1697	! phis(iloop,jloop)/g0
1698	do kloop=2,altlength
1699	! compute sigma level of layer
1700	sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
1701
1702	! compute "mean" temperature of the layer
1703	if(temp(iloop,jloop,kloop,tloop).eq. &
1704	temp(iloop,jloop,kloop-1,tloop)) then
1705	Tmean=temp(iloop,jloop,kloop,tloop)
1706	else
1707	Tmean=(temp(iloop,jloop,kloop,tloop)- &
1708	temp(iloop,jloop,kloop-1,tloop))/ &
1709	log(temp(iloop,jloop,kloop,tloop)/ &
1710	temp(iloop,jloop,kloop-1,tloop))
1711	endif
1712
1713	! compute gravitational acceleration (at altitude zaeroid(kloop-1))
1714	! NB: zareoid=zsurface+phis/g0
1715	gz=g0(a0*2)/ &
1716	(a0+zs_gcm(iloop,jloop,kloop-1,tloop)+phis(iloop,jloop)/g0)**2
1717
1718	! compute zs_gcm(iloop,jloop,kloop,tloop)
1719	zs_gcm(iloop,jloop,kloop,tloop)=zs_gcm(iloop,jloop,kloop-1,tloop)- &
1720	log(sigma(kloop)/sigma(kloop-1))RTmean/gz
1721
1722	! compute za_gcm(kloop)
1723	! za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
1724	! phis(iloop,jloop)/g0
1725	enddo ! kloop
1726	enddo ! iloop
1727	enddo ! jloop
1728	enddo ! tloop
1729
1730	! Cleanup
1731	deallocate(sigma)
1732	!deallocate(zlocal)
1733
1734	end subroutine crude_gcm_zs
1735
1736	!===============================================================================
1737
1738	!#include"crude_gcm_alt.F90"
1739	subroutine crude_gcm_za(lonlength,latlength,altlength,timelength, &
1740	phis,ps,press,temp,za_gcm)
1741	!==============================================================================
1742	! Purpose: Integrate hydrostatic equation in order to build the "above areoid
1743	! altitudes" corresponding to GCM atmospheric levels
1744	!==============================================================================
1745	implicit none
1746	!==============================================================================
1747	! Arguments:
1748	!==============================================================================
1749	integer,intent(in) :: lonlength ! # of points along longitude
1750	integer,intent(in) :: latlength ! # of points along latitude
1751	integer,intent(in) :: altlength ! # of points along altitude
1752	integer,intent(in) :: timelength ! # of stored time steps
1753	real,dimension(lonlength,latlength),intent(in) :: phis
1754	! phis(:,:) is the ground geopotential
1755	real,dimension(lonlength,latlength,timelength),intent(in) :: ps
1756	! ps(:,:) is the surface pressure
1757	real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: press
1758	! press(:,:,:,:) is atmospheric pressure
1759	real,dimension(lonlength,latlength,altlength,timelength),intent(in) :: temp
1760	! temp(:,:,:,:) is atmospheric temperature
1761
1762	real,dimension(lonlength,latlength,altlength,timelength),intent(out) :: za_gcm
1763	! za_gcm(:,:,:,:) are the above aroid heights of GCM levels
1764
1765	!===============================================================================
1766	! Local variables:
1767	!===============================================================================
1768	real,dimension(:),allocatable :: sigma ! GCM sigma levels
1769	real,parameter :: R=191 ! molecular gas constant
1770	real,dimension(:,:,:,:),allocatable :: zlocal ! local above surface altitude
1771	real :: Tmean ! "mean" value of temperature for a given level
1772	integer iloop,jloop,kloop,tloop
1773
1774	! Parameters needed to integrate hydrostatic equation:
1775	real,parameter :: g0=3.7257964
1776	!g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
1777	real,parameter :: a0=3396.E3
1778	!a0: 'mean' Mars radius=3396.km
1779	real :: gz
1780	! gz: gravitational acceleration at a given (zareoid) altitude
1781
1782	!===============================================================================
1783	! 1. Various initialisations
1784	!===============================================================================
1785	allocate(sigma(altlength))
1786	allocate(zlocal(lonlength,latlength,altlength,timelength))
1787
1788	!===============================================================================
1789	! 2. Integrate hydrostatic equation to compute zlocal and za_gcm
1790	!===============================================================================
1791	do tloop=1,timelength
1792	do jloop=1,latlength
1793	do iloop=1,lonlength
1794	! handle case of first altitude level
1795	sigma(1)=press(iloop,jloop,1,tloop)/ps(iloop,jloop,tloop)
1796	zlocal(iloop,jloop,1,tloop)=-log(sigma(1))* &
1797	R*temp(iloop,jloop,1,tloop)/g0
1798	za_gcm(iloop,jloop,1,tloop)=zlocal(iloop,jloop,1,tloop)+ &
1799	phis(iloop,jloop)/g0
1800	do kloop=2,altlength
1801	! compute sigma level of layer
1802	sigma(kloop)=press(iloop,jloop,kloop,tloop)/ps(iloop,jloop,tloop)
1803
1804	! compute "mean" temperature of the layer
1805	if(temp(iloop,jloop,kloop,tloop).eq. &
1806	temp(iloop,jloop,kloop-1,tloop)) then
1807	Tmean=temp(iloop,jloop,kloop,tloop)
1808	else
1809	Tmean=(temp(iloop,jloop,kloop,tloop)- &
1810	temp(iloop,jloop,kloop-1,tloop))/ &
1811	log(temp(iloop,jloop,kloop,tloop)/ &
1812	temp(iloop,jloop,kloop-1,tloop))
1813	endif
1814
1815	! compute gravitational acceleration (at altitude zaeroid(kloop-1))
1816	gz=g0(a02)/(a0+za_gcm(iloop,jloop,kloop-1,tloop))*2
1817
1818	! compute zlocal(kloop)
1819	zlocal(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop-1,tloop)- &
1820	log(sigma(kloop)/sigma(kloop-1))RTmean/gz
1821
1822	! compute za_gcm(kloop)
1823	za_gcm(iloop,jloop,kloop,tloop)=zlocal(iloop,jloop,kloop,tloop)+ &
1824	phis(iloop,jloop)/g0
1825	enddo ! kloop
1826	enddo ! iloop
1827	enddo ! jloop
1828	enddo ! tloop
1829
1830	! Cleanup
1831	deallocate(sigma)
1832	deallocate(zlocal)
1833
1834	end subroutine crude_gcm_za
1835
1836	!===============================================================================
1837
1838	subroutine p_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, &
1839	missing,ps,press,gcmdata,plevels,newdata)
1840	!==============================================================================
1841	! Purpose:
1842	! Recast a 4D (spatio-temporal) GCM dataset, which has a vertical coordinate
1843	! in pseudo-altitude, into a dataset which has a vertical coordinate at given
1844	! pressure levels
1845	!==============================================================================
1846	implicit none
1847	!==============================================================================
1848	! Arguments:
1849	!==============================================================================
1850	integer,intent(in) :: lonlen ! # of points along longitude (GCM dataset)
1851	integer,intent(in) :: latlen ! # of points along latitude (GCM dataset)
1852	integer,intent(in) :: altlen ! # of points along altitude (GCM dataset)
1853	integer,intent(in) :: tlen ! # of stored time steps (GCM dataset)
1854	integer,intent(in) :: newaltlen ! # of points along altitude
1855	real,intent(in) :: missing ! missing value
1856	real,dimension(lonlen,latlen,tlen),intent(in) :: ps ! GCM surface pressure
1857	real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: press ! GCM pressure
1858	real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset
1859	real,dimension(newaltlen),intent(in) :: plevels
1860	! plevels(:) pressure levels at which gcmdata is to be interpolated
1861
1862	real,dimension(lonlen,latlen,newaltlen,tlen),intent(out) :: newdata
1863	! newdata(:,:,:,:) gcmdata recasted along vertical coordinate
1864	!==============================================================================
1865	! Local variables:
1866	!==============================================================================
1867	real,dimension(:),allocatable :: lnp
1868	! lnp(:): used to store log(pressure) values
1869	real,dimension(:),allocatable :: q
1870	! q(:): used to store values along vertical direction
1871	real :: x,y
1872	! x,y: temporary variables
1873	integer :: iloop,jloop,kloop,tloop
1874
1875	allocate(lnp(altlen))
1876	allocate(q(altlen))
1877
1878	do tloop=1,tlen
1879	do jloop=1,latlen
1880	do iloop=1,lonlen
1881	! build lnp() and q() along vertical coordinate
1882	do kloop=1,altlen
1883	lnp(kloop)=-log(press(iloop,jloop,kloop,tloop))
1884	q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
1885	enddo
1886
1887	! Interpolation
1888	do kloop=1,newaltlen
1889	! compute, by interpolation, value at pressure level plevels(kloop)
1890	if ((plevels(kloop).le.ps(iloop,jloop,tloop)).and. &
1891	(plevels(kloop).ge.press(iloop,jloop,altlen,tloop))) then
1892	x=-log(plevels(kloop))
1893	call interpolf(x,y,missing,lnp,q,altlen)
1894	newdata(iloop,jloop,kloop,tloop) = y
1895	else ! if plevels(kloop) is out of range,
1896	! assign a "missing_value" at this node
1897	newdata(iloop,jloop,kloop,tloop) = missing
1898	endif
1899	enddo
1900
1901	enddo !iloop
1902	enddo !jloop
1903	enddo !tloop
1904
1905	! Cleanup
1906	deallocate(lnp)
1907	deallocate(q)
1908
1909	end subroutine p_coord_interp
1910
1911	!===============================================================================
1912
1913	!#include"za_coord_interp.F90"
1914	subroutine z_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, &
1915	missing,z_gcm,gcmdata,flag,z_new,newdata)
1916	!==============================================================================
1917	! Purpose:
1918	! Recast a 4D (spatio-temporal) GCM dataset 'gcmdata', for which corresponding
1919	! grid values of altitude are known 'z_gcm', onto a new altitude grid 'z_new'.
1920	! "Altitudes" can be above areoid or above local surface altitudes, as long as
1921	! both 'z_gcm' and 'znew' refer to altitudes of a same type.
1922	! Note: If altitudes in 'znew' fall outside of the range of altitudes
1923	! in 'z_gcm' then corresponding 'newdata' value is set to 'missing'.
1924	! 'z_gcm' and 'znew' altitudes must be monotone increasing sequences.
1925	!==============================================================================
1926	implicit none
1927	!==============================================================================
1928	! Arguments:
1929	!==============================================================================
1930	integer,intent(in) :: lonlen ! # of points along longitude (GCM dataset)
1931	integer,intent(in) :: latlen ! # of points along latitude (GCM dataset)
1932	integer,intent(in) :: altlen ! # of points along altitude (GCM dataset)
1933	integer,intent(in) :: tlen ! # of stored time steps (GCM dataset)
1934	integer,intent(in) :: newaltlen ! # of points along altitude
1935	real,intent(in) :: missing ! missing value
1936	real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: z_gcm
1937	!z_gcm(:,:,:,) GCM grid points altitude (above areoid or above surface)
1938	real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset
1939	integer,intent(in) :: flag ! flag (==1 for 'log' interpolation)
1940	! flag==0 (standard linear interpolation)
1941	real,dimension(newaltlen),intent(in) :: z_new
1942	! z_new(:) altitudes (above areoid or surface) at which data must be recast
1943
1944	real,dimension(lonlen,latlen,newaltlen,tlen),intent(out) :: newdata
1945	! newdata(:,:,:,:) gcmdata recasted along vertical coordinate
1946
1947	!==============================================================================
1948	! Local variables:
1949	!==============================================================================
1950	real,dimension(:),allocatable :: za,q
1951	real,dimension(:),allocatable :: logq
1952	! za(:): used to store z_gcm at a given location
1953	! q(:): used to store field values along the vertical direction
1954	real :: x,y ! temporary variables
1955	integer :: iloop,jloop,kloop,tloop
1956
1957	allocate(za(altlen))
1958	allocate(q(altlen))
1959	allocate(logq(altlen))
1960
1961	if (flag.eq.0) then
1962	do tloop=1,tlen
1963	do jloop=1,latlen
1964	do iloop=1,lonlen
1965	do kloop=1,altlen
1966	! extract the vertical coordinates
1967	za(kloop)=z_gcm(iloop,jloop,kloop,tloop)
1968	! store values along altitude
1969	q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
1970	enddo !kloop
1971
1972	! Interpolation
1973	do kloop=1,newaltlen
1974	! Check if z_new(kloop) is within range
1975	if ((z_new(kloop).ge.z_gcm(iloop,jloop,1,tloop)).and. &
1976	(z_new(kloop).le.z_gcm(iloop,jloop,altlen,tloop))) then
1977	! z_new(kloop) is within range
1978	x=z_new(kloop)
1979	call interpolf(x,y,missing,za,q,altlen)
1980	newdata(iloop,jloop,kloop,tloop)=y
1981	else ! z_new(kloop) is out of range
1982	newdata(iloop,jloop,kloop,tloop)=missing
1983	endif
1984	enddo !kloop
1985	enddo !iloop
1986	enddo !jloop
1987	enddo !tloop
1988	else ! case when flag=1 (i.e.: rho)
1989	do tloop=1,tlen
1990	do jloop=1,latlen
1991	do iloop=1,lonlen
1992	do kloop=1,altlen
1993	! extract the vertical coordinates
1994	za(kloop)=z_gcm(iloop,jloop,kloop,tloop)
1995	! store log values along altitude
1996	logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop))
1997	enddo !kloop
1998
1999	! Interpolation
2000	do kloop=1,newaltlen
2001	! Check if z_new(kloop) is within range
2002	if ((z_new(kloop).ge.z_gcm(iloop,jloop,1,tloop)).and. &
2003	(z_new(kloop).le.z_gcm(iloop,jloop,altlen,tloop))) then
2004	! z_new(kloop) is within range
2005	x=z_new(kloop)
2006	call interpolf(x,y,missing,za,logq,altlen)
2007	newdata(iloop,jloop,kloop,tloop)=exp(y)
2008	else ! z_new(kloop) is out of range
2009	newdata(iloop,jloop,kloop,tloop)=missing
2010	endif
2011	enddo !kloop
2012	enddo !iloop
2013	enddo !jloop
2014	enddo !tloop
2015	endif
2016
2017	! Cleanup
2018	deallocate(za)
2019	deallocate(q)
2020	deallocate(logq)
2021
2022	end subroutine z_coord_interp
2023
2024	!===============================================================================
2025
2026	subroutine zs_coord_interp(lonlen,latlen,altlen,tlen,newaltlen, &
2027	missing,z_gcm,gcmdata,phis,press,temp,rho, &
2028	have_rho,flag,z_new,newdata)
2029	!==============================================================================
2030	! Purpose:
2031	! Recast a 4D (spatio-temporal) GCM dataset 'gcmdata', for which corresponding
2032	! grid values of altitude are known 'z_gcm', onto a new altitude grid 'z_new'.
2033	! "Altitudes" must be above local surface altitudes.
2034	! Notes:
2035	! 'z_gcm' and 'znew' altitudes must be monotone increasing sequences.
2036	! If altitudes in 'znew(i)' fall below those in 'z_gcm(:,:,1,:)', then
2037	! 'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,1,:)' if
2038	! flag=0, and extrapolated (exponentially) if flag=1.
2039	! If altitudes in 'znew(i)' are above those in 'z_gcm(:,:,altlen,:)', then
2040	! 'newdata(:,:,i,:)' is set to constant value 'gcmdata(:,:,altlen,:)'
2041	! if flag=0, and extrapolated (exponentially) if flag=1.
2042	!==============================================================================
2043	implicit none
2044	!==============================================================================
2045	! Arguments:
2046	!==============================================================================
2047	integer,intent(in) :: lonlen ! # of points along longitude (GCM dataset)
2048	integer,intent(in) :: latlen ! # of points along latitude (GCM dataset)
2049	integer,intent(in) :: altlen ! # of points along altitude (GCM dataset)
2050	integer,intent(in) :: tlen ! # of stored time steps (GCM dataset)
2051	integer,intent(in) :: newaltlen ! # of points along altitude
2052	real,intent(in) :: missing ! missing value
2053	real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: z_gcm
2054	!z_gcm(:,:,:,) GCM grid points altitude (above areoid or above surface)
2055	real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: gcmdata ! GCM dataset
2056	real,dimension(lonlen,latlen),intent(in) :: phis
2057	! phis(:,:) is the ground geopotential
2058	real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: press
2059	! press(:,:,:,:) is atmospheric pressure on GCM levels
2060	real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: temp
2061	! temp(:,:,:,:) is atmospheric temperature on GCM levels
2062	real,dimension(lonlen,latlen,altlen,tlen),intent(in) :: rho
2063	! rho(:,:,:,:) is atmospheric density on GCM levels
2064	logical,intent(in) :: have_rho ! trueif we have density at hand
2065	integer,intent(in) :: flag ! flag (==1 for 'log' interpolation)
2066	! flag==0 (standard linear interpolation)
2067	real,dimension(newaltlen),intent(in) :: z_new
2068	! z_new(:) altitudes (above areoid or surface) at which data must be recast
2069
2070	real,dimension(lonlen,latlen,newaltlen,tlen),intent(out) :: newdata
2071	! newdata(:,:,:,:) gcmdata recasted along vertical coordinate
2072
2073	!==============================================================================
2074	! Local variables:
2075	!==============================================================================
2076	real,dimension(:),allocatable :: z,q
2077	real,dimension(:),allocatable :: logq
2078	! z(:): used to store z_gcm at a given location
2079	! q(:): used to store field values along the vertical direction
2080	real,dimension(:,:,:),allocatable :: Rbottom,Rtop
2081	! values of R (gas constant) below and above GCM layers
2082	real :: x,y ! temporary variables
2083	integer :: iloop,jloop,kloop,tloop
2084	real,parameter :: g0=3.7257964
2085	!g0: exact mean gravity at radius=3396.km (Lemoine et al. 2001)
2086	real,parameter :: a0=3396.E3
2087	!a0: 'mean' Mars radius=3396.km
2088	real,parameter :: Rmean=191 ! molecular gas constant
2089	real :: gz
2090	! gz: gravitational acceleration at a given (above areoid) altitude
2091
2092	allocate(z(altlen))
2093	allocate(q(altlen))
2094	allocate(logq(altlen))
2095
2096	! 1. Build Rbottom and Rtop (only necessary if flag=1)
2097	if (flag.eq.1) then
2098	allocate(Rbottom(lonlen,latlen,tlen))
2099	allocate(Rtop(lonlen,latlen,tlen))
2100	if (have_rho) then
2101	do iloop=1,lonlen
2102	do jloop=1,latlen
2103	do tloop=1,tlen
2104	Rbottom(iloop,jloop,tloop)=press(iloop,jloop,1,tloop)/ &
2105	(rho(iloop,jloop,1,tloop)* &
2106	temp(iloop,jloop,1,tloop))
2107	Rtop(iloop,jloop,tloop)=press(iloop,jloop,altlen,tloop)/ &
2108	(rho(iloop,jloop,altlen,tloop)* &
2109	temp(iloop,jloop,altlen,tloop))
2110	enddo
2111	enddo
2112	enddo
2113	else ! we don't have density at hand; use mean molecular gas constant value
2114	Rbottom(:,:,:)=Rmean
2115	Rtop(:,:,:)=Rmean
2116	endif
2117	endif ! if (flag.eq.1)
2118
2119	! 2. Interpolation
2120	if (flag.eq.0) then
2121	do tloop=1,tlen
2122	do jloop=1,latlen
2123	do iloop=1,lonlen
2124	! preliminary stuff
2125	do kloop=1,altlen
2126	! extract the vertical coordinates
2127	z(kloop)=z_gcm(iloop,jloop,kloop,tloop)
2128	! store values along altitude
2129	q(kloop)=gcmdata(iloop,jloop,kloop,tloop)
2130	enddo !kloop
2131
2132	! Interpolation
2133	do kloop=1,newaltlen
2134	if (z_new(kloop).lt.z_gcm(iloop,jloop,1,tloop)) then
2135	! z_new(kloop) is below lowest GCM level
2136	newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,1,tloop)
2137	else if (z_new(kloop).gt.z_gcm(iloop,jloop,altlen,tloop)) then
2138	! z_new(kloop) is above highest GCM level
2139	newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,altlen,tloop)
2140	else ! z_new(kloop) is within range
2141	x=z_new(kloop)
2142	call interpolf(x,y,missing,z,q,altlen)
2143	newdata(iloop,jloop,kloop,tloop)=y
2144	endif
2145	enddo !kloop
2146	enddo !iloop
2147	enddo !jloop
2148	enddo !tloop
2149	else ! case when flag=1 (i.e.: rho or pressure)
2150	do tloop=1,tlen
2151	do jloop=1,latlen
2152	do iloop=1,lonlen
2153	! preliminary stuff
2154	do kloop=1,altlen
2155	! extract the vertical coordinates
2156	z(kloop)=z_gcm(iloop,jloop,kloop,tloop)
2157	! store log values along altitude
2158	logq(kloop)=log(gcmdata(iloop,jloop,kloop,tloop))
2159	enddo !kloop
2160
2161	! Interpolation
2162	do kloop=1,newaltlen
2163	if (z_new(kloop).lt.z_gcm(iloop,jloop,1,tloop)) then
2164	! z_new(kloop) is below lowest GCM level
2165	newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,1,tloop)* &
2166	exp(((z_gcm(iloop,jloop,1,tloop)-z_new(kloop))*g0)/ &
2167	(Rbottom(iloop,jloop,tloop)* &
2168	temp(iloop,jloop,1,tloop)))
2169	else if (z_new(kloop).gt.z_gcm(iloop,jloop,altlen,tloop)) then
2170	! z_new(kloop) is above highest GCM level
2171	! NB: zareoid=zsurface+phis/g0
2172	gz=g0a0a0/(a0+z_new(kloop)+phis(iloop,jloop)/g0)**2
2173	newdata(iloop,jloop,kloop,tloop)=gcmdata(iloop,jloop,altlen,tloop)* &
2174	exp(((z_gcm(iloop,jloop,altlen,tloop)-z_new(kloop))*gz)/ &
2175	(Rtop(iloop,jloop,tloop)* &
2176	temp(iloop,jloop,altlen,tloop)))
2177	else ! z_new(kloop) is within range
2178	x=z_new(kloop)
2179	call interpolf(x,y,missing,z,logq,altlen)
2180	newdata(iloop,jloop,kloop,tloop)=exp(y)
2181	endif
2182	enddo !kloop
2183	enddo !iloop
2184	enddo !jloop
2185	enddo !tloop
2186	endif
2187
2188	! Cleanup
2189	deallocate(z)
2190	deallocate(q)
2191	deallocate(logq)
2192	if (flag.eq.1) then
2193	deallocate(Rbottom)
2194	deallocate(Rtop)
2195	endif
2196
2197	end subroutine zs_coord_interp
2198
2199	!===============================================================================
2200
2201	subroutine interpolf(x,y,missing,xd,yd,nd)
2202	!==============================================================================
2203	! Purpose:
2204	! Yield y=f(x), where xd() end yd() are arrays of known values,
2205	! using linear interpolation
2206	! If x is not included in the interval spaned by xd(), then y is set
2207	! to a default value 'missing'
2208	! Note:
2209	! Array xd() should contain ordered (either increasing or decreasing) abscissas
2210	!==============================================================================
2211	implicit none
2212	!==============================================================================
2213	! Arguments:
2214	!==============================================================================
2215	real,intent(in) :: x ! abscissa at which interpolation is to be done
2216	real,intent(in) :: missing ! missing value (if no interpolation is performed)
2217	integer :: nd ! size of arrays
2218	real,dimension(nd),intent(in) :: xd ! array of known absissas
2219	real,dimension(nd),intent(in) :: yd ! array of correponding values
2220
2221	real,intent(out) :: y ! interpolated value
2222	!==============================================================================
2223	! Local variables:
2224	!==============================================================================
2225	integer :: i
2226
2227	! default: set y to 'missing'
2228	y=missing
2229
2230	do i=1,nd-1
2231	if (((x.ge.xd(i)).and.(x.le.xd(i+1))).or.&
2232	((x.le.xd(i)).and.(x.ge.xd(i+1)))) then
2233	y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
2234	exit
2235	endif
2236	enddo
2237
2238
2239	end subroutine interpolf

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

Download in other formats:

Original Format