Context Navigation

zrecast.F90 @ 937

Last change on this file since 937 was 816, checked in by slebonnois, 12 years ago
SL: tools for postprocessing (Veznus and Titan); see DOC/documentation/vt-tools.pdf
File size: 85.7 KB

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

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

Download in other formats:

Original Format