Context Navigation

anl_stats-diag.F @ 937

Last change on this file since 937 was 410, checked in by acolaitis, 13 years ago
Modified NCDF norm of our files from classic to 64bit offset to support variables indices of more than integer*4 maximum length. This is a priori retrocompatible with classic format. I have tested it, it works for all file outputs, newstart (on classic or 64-bit offset files). One can check the format of his .nc with ncdump -k file.nc
File size: 44.2 KB

Line
1	PROGRAM anldiag_nc
2	IMPLICIT NONE
3	c======================================================================
4	c
5	c
6	c Program designed to read output files from the MARS gcm
7	c or the Mars Climate database.
8	c
9	c Francois Forget 1999
10	c Version used by monica to interpolate in zaeroid and p coordinate
11
12	c
13	c=======================================================================
14	c declarations:
15	c -------------
16
17	#include "dimensions.h"
18	#include "paramet.h"
19	#include "comconst.h"
20	#include "comdissip.h"
21	#include "comvert.h"
22	#include "comgeom2.h"
23	#include "logic.h"
24	#include "temps.h"
25	#include "control.h"
26	#include "ener.h"
27	#include "description.h"
28	#include "netcdf.inc"
29
30	INTEGER itau,nbpas,nbpasmx, coor,sor,varoutsize
31	integer varspecsize
32	integer varout(6),varspec(nqmx)
33	integer itau0
34	PARAMETER(nbpasmx=1000000)
35	REAL temps(nbpasmx),y
36	INTEGER unitlec
37	INTEGER i,j,l,irec,itautot,ierr,iq,n,ff,k
38	real, dimension(:),allocatable :: lat,lon,alt
39	integer:: nout,latdimout,londimout,altdimout,timedimout,timevarout
40	integer:: latdim,londim,altdim,dimid
41	integer:: latvar,lonvar,altvar,timevar
42	integer:: latlen,lonlen,altlen,timelen
43	integer, dimension(4) :: edges,corner,id
44	integer, dimension(3) :: edges2,corner2
45
46	INTEGER unit,nvarid,nid,varid
47	REAL mugaz
48
49	REAL missing
50	PARAMETER(missing=1E+20)
51	REAL valid_range(2)
52	DATA valid_range /0., 300.0/
53	c variables meteo dans la grille verticale GCM
54	c --------------------------------------------
55	REAL v(iip1,jjp1,llm),u(iip1,jjp1,llm),w(iip1,jjp1,llm)
56	REAL t(iip1,jjp1,llm),rho(iip1,jjp1,llm),phi(iip1,jjp1,llm)
57	REAL ps(iip1,jjp1) , tsurf(iip1,jjp1)
58	REAL Rnew(iip1,jjp1,llm)
59	real euv(iip1,jjp1,llm),con(iip1,jjp1,llm),nir(iip1,jjp1,llm)
60	real nspec(iip1,jjp1,llm,nqmx),spec(iip1,jjp1,llm)
61
62	REAL v_sd(iip1,jjp1,llm),u_sd(iip1,jjp1,llm)
63	real w_sd(iip1,jjp1,llm),t_sd(iip1,jjp1,llm)
64	real rho_sd(iip1,jjp1,llm)
65	real euv_sd(iip1,jjp1,llm),con_sd(iip1,jjp1,llm)
66	real nir_sd(iip1,jjp1,llm)
67	real nspec_sd(iip1,jjp1,llm,nqmx)
68
69	c ALtitude of the GCM levels (+1 dummy "thermosphere level")
70	c -------------------------
71	REAL zlocal(iip1,jjp1,llm+1) ! altitude above the local surface (m)
72	REAL zareoid(iip1,jjp1,llm+1) ! altitude above the mola areoid (m)
73	real tmean ! use to integrate hydrostatic equation
74	real sig_therm ! dummy "thermosphere level" as in atmemcd.F
75	real T_therm
76
77	c variables meteo en coordonnee de pression
78	c --------------------------------------------
79	REAL vp(iip1,jjp1,llm),up(iip1,jjp1,llm),wp(iip1,jjp1,llm)
80	REAL tp(iip1,jjp1,llm),rhop(iip1,jjp1,llm)
81	REAL nspecp(iip1,jjp1,llm,nqmx)
82	REAL phip(iip1,jjp1,llm)
83
84	REAL vp_sd(iip1,jjp1,llm),up_sd(iip1,jjp1,llm)
85	real wp_sd(iip1,jjp1,llm),tp_sd(iip1,jjp1,llm)
86	real rhop_sd(iip1,jjp1,llm),nspecp_sd(iip1,jjp1,llm,nqmx)
87
88	c Niveaux de pression (Pa)
89	real pref(llm)
90
91	c Version 32 layers : ************
92	data pref/
93	&1000.,884., 783. , 610., 475, 370, 288, 224, 174, 140.,
94	&115 , 80.6481, 48.0445, 26.6881, 14.1462, 7.29102, 3.70004,
95	&1.86241, 0.933516, 0.466924, 0.233296, 0.116503, 5.81631E-02,
96	& 2.90335E-02,1.44919E-02, 7.23328E-03, 3.61027E-03, 1.80193E-03,
97	& 8.99364E-04,4.48881E-04, 2.15769E-04, 1.3E-04/
98
99	c Version 50 layers : ************
100	c data pref/
101	c &1000.,884., 783. , 610., 475, 370, 288, 224, 174, 140.,
102	c &115 , 80.6481, 48.0445, 26.6881, 14.1462, 7.29102, 3.70004,
103	c &1.86241, 0.933516, 0.466924, 0.233296, 0.116503, 5.81631E-02,
104	c & 2.90335E-02,1.44919E-02, 7.23328E-03, 3.61027E-03, 1.80193E-03,
105	c & 8.99364E-04,4.48881E-04, 2.15769E-04, 1.3E-04,
106	c & 7.23328e-05, 3.61027e-05, 1.80193e-05, 8.99364e-06,
107	c & 4.48881e-06, 2.15769e-06, 1.3e-06,
108	c & 7.23328e-07, 3.61027e-07, 1.80193e-07, 8.99364e-08,
109	c & 4.48881e-08, 2.15769e-08, 1.3e-08,
110	c & 7.23328e-09, 3.61027e-09, 1.80193e-09, 8.99364e-10/
111
112	c data pref/
113	c &1000.,884., 783. , 610., 475, 370, 288, 224, 174, 140.,
114	c &115 , 80.6481, 48.0445, 26.6881, 14.1462, 7.29102, 3.70004,
115	c &1.86241, 0.933516, 0.466924, 0.233296, 0.116503, 5.81631E-02,
116	c & 2.90335E-02,1.44919E-02 /
117	c data pref/
118	c & 101.4369,94.4369,87.4369,80.4369,73.4369,66.4369,59.4369,52.4369,
119	c & 45.4369,38.4369,31.5527,25.0626,18.9666,13.5138,8.97780,
120	c & 5.54010,3.19040,1.73630,0.907000,0.460400,0.228200,0.109700,
121	c & 0.0500,0.0200,0.0050/
122	c variables meteo en coordonnee de zareoid
123	c --------------------------------------------
124	REAL va(iip1,jjp1,llm),ua(iip1,jjp1,llm),wa(iip1,jjp1,llm)
125	REAL ta(iip1,jjp1,llm),rhoa(iip1,jjp1,llm),ra(iip1,jjp1,llm)
126	real euva(iip1,jjp1,llm),cona(iip1,jjp1,llm),nira(iip1,jjp1,llm)
127	real nspeca(iip1,jjp1,llm,nqmx)
128
129	REAL va_sd(iip1,jjp1,llm),ua_sd(iip1,jjp1,llm)
130	real wa_sd(iip1,jjp1,llm),ta_sd(iip1,jjp1,llm)
131	real rhoa_sd(iip1,jjp1,llm)
132	real euva_sd(iip1,jjp1,llm),cona_sd(iip1,jjp1,llm)
133	real nira_sd(iip1,jjp1,llm)
134	real nspeca_sd(iip1,jjp1,llm,nqmx)
135
136	REAL vij(llm),uij(llm),zaij(llm),wij(llm),rij(llm)
137	REAL tij(llm),rhoij(llm+1),sigij(llm)
138	real euvij(llm),conij(llm),nirij(llm)
139	real nspecij(llm+1,nqmx)
140
141	REAL vij_sd(llm),uij_sd(llm)
142	real wij_sd(llm),tij_sd(llm),rhoij_sd(llm+1)
143	real euvij_sd(llm),conij_sd(llm),nirij_sd(llm)
144	real nspecij_sd(llm+1,nqmx)
145
146	c Niveaux de zareoide (m)
147	real zaref(llm)
148	real p_zaref,p_zaij(llm+1) ! used to interpolate rho
149
150	c Surface data (sometime needed for some analysis):
151	c ------------------------------------------------
152	character relief*3
153	real phis(iip1,jjp1)
154	real alb(iip1,jjp1),ith(iip1,jjp1)
155	real zmeaS(iip1,jjp1),zstdS(iip1,jjp1)
156	real zsigS(iip1,jjp1),zgamS(iip1,jjp1),ztheS(iip1,jjp1)
157
158
159	real utime
160	real localtime(iip1)
161	common/temporaire/localtime
162
163	INTEGER*4 day0
164	integer nmoy(jjp1,llm),tp1,tp2,pass
165
166	CHARACTER*1 yes,yescomp
167
168	integer iformat
169
170	c declarations de l'interface avec mywrite:
171	c -----------------------------------------
172
173	CHARACTER file*80
174	CHARACTER nomfich*60
175	CHARACTER nomfile(2)60,nomfileout(2)60
176	integer fichsize,fstats
177	CHARACTER fdiag*60
178
179	c externe:
180	c --------
181
182	EXTERNAL iniconst,inigeom,covcont,mywrite
183	EXTERNAL inifilr,exner
184	EXTERNAL solarlong,coordij,moy2
185	EXTERNAL SSUM
186	REAL SSUM
187	EXTERNAL lnblnk
188	INTEGER lnblnk
189
190	c Dust
191	c ----
192
193	#include "fxyprim.h"
194
195	character*10 noms(nqmx),varnoms(6)
196	data noms/"co2","co","o","o1d","o2","o3","h","h2",
197	$ "oh","ho2","h2o2", "n2", "h2o"/
198	data varnoms/"T","U","V","W","RHO","Species"/
199
200
201	c-----------------------------------------------------------------------
202	c initialisations:
203	c ----------------
204
205	unitlec=11
206	itautot=0
207
208
209	c Lecture du fichier a lire
210	c -------------------------
211
212	print*,' '
213	PRINT*,'File to process: 1)DIAGFI 2)STATS 3)BOTH'
214	READ(5,'(i)') fichsize
215	if(fichsize.ne.2)then
216	i=1
217	PRINT*,' enter name of diagfi file (w/o .nc)'
218	READ(5,'(a)') nomfile(i)
219	fdiag=nomfile(i)
220	PRINT*,' do you want to rename the output file?
221	& 1) Yes 2) No'
222	READ(5,'(i)') sor
223	if (sor.eq.1) then
224	PRINT*,' enter label for output file'
225	READ(5,'(a)') nomfileout(i)
226	elseif(sor.eq.2) then
227	nomfileout(i)=nomfile(i)
228	endif
229	endif
230	if(fichsize.ne.1)then
231	i=i+1
232	PRINT*,' enter name of stats file (w/o .nc)'
233	READ(5,'(a)') nomfile(i)
234	PRINT*,' do you want to rename the output file?
235	& 1) Yes 2) No'
236	READ(5,'(i)') sor
237	if (sor.eq.1) then
238	PRINT*,' enter label for output file'
239	READ(5,'(a)') nomfileout(i)
240	elseif(sor.eq.2) then
241	nomfileout(i)=nomfile(i)
242	endif
243	if(fichsize.eq.2) then
244	PRINT*,' enter name of diagfi file for controle (w/o .nc)'
245	read(5,'(a)') fdiag
246	endif
247	fstats=1
248	endif
249	if (fichsize.ge.2) fichsize=fichsize-1
250	print*,' '
251	print*,'Output in: 1) netcdf 2) IDL 3) both ?'
252	READ(5,'(i)') sor
253	print*,' '
254	print*,'Coordinates in: 1) pressure 2) Zareoid 3) both ?'
255	READ(5,'(i)') coor
256	print*,' '
257	print*,'Fields: T U V W Rho Species All ?'
258	print*,'Type total number of desired outputs or 6 for All '
259	READ(5,'(i)') varoutsize
260	if(varoutsize.lt.6) then
261	print*,' '
262	print*,'Fields: T U V W Rho Species '
263	print*,' # : 1 2 3 4 5 6'
264	print*,' '
265	do i=1,varoutsize
266	print*,'Type number of variable ',i
267	READ(5,'(i)') varout(i)
268	enddo
269	endif
270	if(varoutsize.eq.6) then
271	do i=1,6
272	varout(i)=i
273	enddo
274	endif
275	do i=1,varoutsize
276	if (varout(i).eq.6) then
277	print*,'Type total number of desired species or 13 for All'
278	READ(5,'(i2)') varspecsize
279	if(varspecsize.ne.13) then
280	print*,' '
281	print*,'Select Species: CO2 CO O O(1D) O2 O3'
282	print*,' 1 2 3 4 5 6'
283	print*,' H H2 OH HO2 H2O2 N2 H2O'
284	print*,' 7 8 9 10 11 12 13'
285	do n=1,varspecsize
286	print*,'Type number of species ',i
287	read(5,'(i)') varspec(n)
288	enddo
289	else
290	do n=1,varspecsize
291	varspec(n)=n
292	enddo
293	endif
294	endif
295	enddo
296	if(varoutsize.gt.6) print*,'Value must be <= 6'
297
298
299	! LOOP ON FILES
300	! =================================================================
301	do ff=1,fichsize
302
303	file=nomfile(ff)
304	nomfich=nomfile(ff)
305	file=file(1:lnblnk(file))
306	PRINT*,'file',file
307
308	c Ouverture fichiers :
309	c ------------------
310	write(,) "opening "//trim(nomfich)//"..."
311	ierr = NF_OPEN(trim(adjustl(nomfich))//".nc",NF_NOWRITE,nid)
312	if (ierr.NE.NF_NOERR) then
313	write(,) 'ERROR: Pb opening file '//nomfich
314	stop ""
315	endif
316
317	! Control & lecture on dimensions
318	! =================================================================
319	ierr=NF_INQ_DIMID(nid,"latitude",latdim)
320	ierr=NF_INQ_VARID(nid,"latitude",latvar)
321	if (ierr.NE.NF_NOERR) then
322	write(,) 'ERROR: Field <latitude> is missing'
323	stop ""
324	endif
325	ierr=NF_INQ_DIMLEN(nid,latdim,latlen)
326	! write(,) "latlen: ",latlen
327
328	ierr=NF_INQ_DIMID(nid,"longitude",londim)
329	ierr=NF_INQ_VARID(nid,"longitude",lonvar)
330	if (ierr.NE.NF_NOERR) then
331	write(,) 'ERROR: Field <longitude> is lacking'
332	stop ""
333	endif
334	ierr=NF_INQ_DIMLEN(nid,londim,lonlen)
335	! write(,) "lonlen: ",lonlen
336
337	ierr=NF_INQ_DIMID(nid,"altitude",altdim)
338	ierr=NF_INQ_VARID(nid,"altitude",altvar)
339	if (ierr.NE.NF_NOERR) then
340	write(,) 'ERROR: Field <altitude> is lacking'
341	stop ""
342	endif
343
344	ierr=NF_INQ_DIMLEN(nid,altdim,altlen)
345	! write(,) "altlen: ",altlen
346
347	if (ff.eq.1) then
348	allocate(lat(latlen))
349	allocate(lon(lonlen))
350	allocate(alt(altlen))
351	endif
352
353	#ifdef NC_DOUBLE
354	ierr = NF_GET_VAR_DOUBLE(nid,latvar,lat)
355	ierr = NF_GET_VAR_DOUBLE(nid,lonvar,lon)
356	ierr = NF_GET_VAR_DOUBLE(nid,altvar,alt)
357	#else
358	ierr = NF_GET_VAR_REAL(nid,latvar,lat)
359	ierr = NF_GET_VAR_REAL(nid,lonvar,lon)
360	ierr = NF_GET_VAR_REAL(nid,altvar,alt)
361	#endif
362
363	c Lecture Time :
364
365	ierr= NF_INQ_DIMID (nid,"Time",dimid)
366	IF (ierr.NE.NF_NOERR) THEN
367	PRINT*, 'anl_NC: Le champ <Time> est absent'
368	CALL abort
369	ENDIF
370
371	ierr= NF_INQ_DIMLEN (nid,dimid,nbpas)
372	ierr = NF_INQ_VARID (nid, "Time", nvarid)
373	#ifdef NC_DOUBLE
374	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, temps)
375	#else
376	ierr = NF_GET_VAR_REAL(nid, nvarid, temps)
377	#endif
378	IF (ierr.NE.NF_NOERR) THEN
379	PRINT*, 'anl_NC: Lecture echouee pour <Time>'
380	CALL abort
381	ENDIF
382	PRINT*,'temps',(temps(itau),itau=1,10)
383	! =================================================================
384
385	! SETTING OUTPUT FILES
386
387	file=nomfileout(ff)
388	nomfich=nomfileout(ff)
389	file=file(1:lnblnk(file))
390
391	c layers in zareoid : ************
392	do l=0,llm-1
393	zaref(l+1)=4500.*l
394	alt(l+1)=zaref(l+1)/1000.
395	enddo
396
397	! =================================================================
398	! Output in netcdf
399	! =================================================================
400	if (sor.ne.2) then
401
402	if(coor.ne.2) then ! PRESSURE COORDINATES
403
404	ierr = NF_CREATE(trim(adjustl(nomfich))//"_P.nc",
405	& IOR(NF_CLOBBER,NF_64BIT_OFFSET), nout)
406	if (ierr.NE.NF_NOERR) THEN
407	write(,)' Pb d''ouverture du fichier '
408	write(,)' ierr = ', ierr
409	STOP
410	ENDIF
411	ierr = NF_PUT_ATT_TEXT(nout, NF_GLOBAL, "title", 24,
412	& "Pressure coordinates ")
413	ierr = NF_DEF_DIM(nout,"latitude", size(lat), latdimout)
414	ierr = NF_DEF_DIM(nout,"longitude", size(lon), londimout)
415	ierr = NF_DEF_DIM(nout,"altitude", size(pref), altdimout)
416	ierr = NF_DEF_DIM(nout,"Time", NF_UNLIMITED, timedimout)
417	ierr = NF_ENDDEF(nout)
418
419	endif
420
421	if(coor.ne.1) then ! ZAEROID COORDINATES
422
423	ierr = NF_CREATE(trim(adjustl(nomfich))//"_Z.nc",
424	& IOR(NF_CLOBBER,NF_64BIT_OFFSET), nout)
425	IF(ierr.NE.NF_NOERR) THEN
426	write(,)' Pb d''ouverture du fichier '
427	write(,)' ierr = ', ierr
428	STOP
429	ENDIF
430	ierr = NF_PUT_ATT_TEXT(nout, NF_GLOBAL, "title", 24,
431	& "zaeroid coordinates")
432	ierr = NF_DEF_DIM(nout,"latitude", size(lat), latdimout)
433	ierr = NF_DEF_DIM(nout,"longitude", size(lon), londimout)
434	ierr = NF_DEF_DIM(nout,"altitude ", size(pref), altdimout)
435	ierr = NF_DEF_DIM(nout,"Time", NF_UNLIMITED, timedimout)
436	ierr = NF_ENDDEF(nout)
437
438	endif
439
440	call def_var(nout,"Time","Time","days since 0000-00-0
441	& 00:00:00",1,(/timedimout/),timevarout,ierr)
442
443	ierr = NF_REDEF (nout)
444	call def_var(nout,"latitude","latitude","degrees_north",1,
445	& (/latdimout/),nvarid,ierr)
446
447	#ifdef NC_DOUBLE
448	ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lat)
449	#else
450	ierr = NF_PUT_VAR_REAL (nout,nvarid,lat)
451	#endif
452	ierr = NF_REDEF (nout)
453	call def_var(nout,"longitude","East longitude","degrees_east",1,
454	& (/londimout/),nvarid,ierr)
455	#ifdef NC_DOUBLE
456	ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,lon)
457	#else
458	ierr = NF_PUT_VAR_REAL (nout,nvarid,lon)
459	#endif
460	ierr = NF_REDEF (nout)
461
462	if(coor.ne.2) then ! PRESSURE CCOORDINATES
463	#ifdef NC_DOUBLE
464	ierr = NF_DEF_VAR (nout,"altitude",NF_DOUBLE,1,
465	& (/altdimout/), nvarid)
466	#else
467	ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,
468	& (/altdimout/),nvarid)
469	#endif
470
471	ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",22,"altitude
472	& pression")
473	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,"Pa")
474	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',4,"down")
475	ierr = NF_ENDDEF(nout)
476	#ifdef NC_DOUBLE
477	ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,pref)
478	#else
479	ierr = NF_PUT_VAR_REAL (nout,nvarid,pref)
480	#endif
481	endif
482
483	if(coor.ne.1) then ! ZAEROID COORDINATES
484
485	#ifdef NC_DOUBLE
486	ierr = NF_DEF_VAR (nout,"altitude",NF_DOUBLE,1,
487	& (/altdimout/), nvarid)
488	#else
489	ierr = NF_DEF_VAR (nout,"altitude",NF_FLOAT,1,
490	& (/altdimout/),nvarid)
491	#endif
492	ierr = NF_PUT_ATT_TEXT (nout,nvarid,"long_name",20,"altitude
493	& zaeroid")
494	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'units',2,"km")
495	ierr = NF_PUT_ATT_TEXT (nout,nvarid,'positive',2,"up")
496	ierr = NF_ENDDEF(nout)
497	#ifdef NC_DOUBLE
498	ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,alt)
499	#else
500	ierr = NF_PUT_VAR_REAL (nout,nvarid,alt)
501	#endif
502	endif
503
504	endif !end if netcdf
505
506	! Output in idl
507	! =================================================================
508	if (sor.ne.1) then
509
510	if(coor.ne.2) then ! PRESSURE COORDINATES
511	do i=1,varoutsize
512	j=121+i
513	if(varout(i).eq.6) then
514	do iq=1,varspecsize
515	file=trim(nomfileout(ff))//'_P_'
516	& //trim(noms(iq))//'.dat'
517	j=121+i+(iq-1)
518	open(j,file=trim(file), status='unknown')
519	write(j,117) noms(iq)
520	write(j,120) iim,jjp1,llm
521	write(j,119) lon
522	write(j,119) lat
523	write(j,119) pref
524	enddo
525	else
526	file=trim(nomfileout(ff))//'_P_'
527	& //trim(varnoms(varout(i)))//'.dat'
528	open(j,file=trim(file), status='unknown')
529	write(j,117) varnoms(varout(i))
530	write(j,120) iim,jjp1,llm
531	write(j,119) lon
532	write(j,119) lat
533	write(j,119) pref
534	endif
535	enddo
536	endif
537
538	if(coor.ne.1) then ! ZAEROID COORDINATES
539	do i=1,varoutsize
540	j=1221+i
541	if(varout(i).eq.6) then
542	do iq=1,varspecsize
543	file=trim(nomfileout(ff))//'_Z_'
544	& //trim(noms(iq))//'.dat'
545	j=1221+i+(iq-1)
546	open(j,file=trim(file), status='unknown')
547	write(j,117) noms(iq)
548	write(j,120) iim,jjp1,llm
549	write(j,119) lon
550	write(j,119) lat
551	write(j,119) alt
552	enddo
553	else
554	file=trim(nomfileout(ff))//'_Z_'
555	& //trim(varnoms(varout(i)))//'.dat'
556	open(j,file=trim(file), status='unknown')
557	write(j,117) varnoms(varout(i))
558	write(j,120) iim,jjp1,llm
559	write(j,119) lon
560	write(j,119) lat
561	write(j,119) alt
562	endif
563	enddo
564
565	c file=trim(nomfileout(ff))
566	c j=j+1
567	c open(j,file=file(1:lnblnk(file))//'_Z_EUV.dat',
568	c & status='unknown')
569	c write(j,117) 'EUV'
570	c write(j,120) iim,jjp1,llm
571	c write(j,119) lon
572	c write(j,119) lat
573	c write(j,119) alt
574	c j=j+1
575	c open(j,file=file(1:lnblnk(file))//'_Z_CON.dat',
576	c & status='unknown')
577	c write(j,117) 'CONDUCTION'
578	c write(j,120) iim,jjp1,llm
579	c write(j,119) lon
580	c write(j,119) lat
581	c write(j,119) alt
582	c j=j+1
583	c open(j,file=file(1:lnblnk(file))//'_Z_NIR.dat',
584	c & status='unknown')
585	c write(j,117) 'NIR'
586	c write(j,120) iim,jjp1,llm
587	c write(j,119) lon
588	c write(j,119) lat
589	c write(j,119) alt
590	endif
591
592	endif
593	120 format(3i3)
594	119 format(12(e10.4,1x))
595	118 format(2i3)
596	117 format(6(a7))
597	116 format(13(a6))
598	! =================================================================
599
600
601	800 continue ! LOOP SUR LES FICHIERS
602
603
604	! =================================================================
605	! INITIALIZATION OF PARAMETERS
606	! =================================================================
607
608	rad=3397200. ! rayon de mars (m) ~3397200 m
609	daysec=88775. ! duree du sol (s) ~88775 s
610	omeg=4.*asin(1.)/(daysec) ! vitesse de rotation (rad.s-1)
611	g=3.72 ! gravite (m.s-2) ~3.72
612	mugaz=43.49 ! Masse molaire de l'atm (g.mol-1) ~43.49
613	kappa=.256793 ! = r/cp ~0.256793
614	r = 191.18213
615	pi=2.*asin(1.)
616	ecritphy =1
617	iphysiq=1
618	day_ini=0.
619	day_step=1
620
621	CALL readhead_NC(trim(adjustl(fdiag))//".nc",day0,phis,R)
622	WRITE (,) 'day0 = ' , day0
623
624	CALL iniconst
625	CALL inigeom
626	CALL inifilr
627
628	c Dummy "thermosphere level" as in atmemcd (extrapol above top GCM level)
629	sig_therm=bps(llm)*exp(-20./7.)
630	T_therm = 200. ! temperature (K) of the "thermosphere level"
631
632
633	c If needed : getting topography, albedo, thermal inertia...
634	c -------------------------------------------------------
635	c "phis" is the surface geopotential = topography (m) * g
636	relief='mola' ! Topography MOLA par defaut
637
638	c CALL datareadnc(relief,phis,alb,ith,zmeaS,zstdS,zsigS,zgamS,
639	c . ztheS)
640
641
642	! =================================================================
643	! LOOP IN TIME
644	! =================================================================
645
646	DO itau=1,nbpas
647
648	! FILE READING
649	! =================================================================
650	print*,'Timestep = ',itau,'/',nbpas
651
652	c call lect_var(nid,'tsurf',3,ierr,itau, tsurf)
653	call lect_var(nid,'ps',3,ierr,itau, ps)
654
655	call lect_var(nid,'temp',4,ierr,itau, t)
656	call lect_var(nid,'u',4,ierr,itau, u)
657	call lect_var(nid,'v',4,ierr,itau, v)
658	c call lect_var(nid,'w',4,ierr,itau, w)
659	c call lect_var(nid,'rho',4,ierr,itau, rho)
660	do l=1,llm
661	do j=1,jjp1
662	do i=1,iip1
663	Rnew(i,j,l)=8.314/mugaz*1.e3
664	! A enregistrer lors du run pour pouvoir le lire avec la ligne dessous
665	enddo
666	enddo
667	enddo
668	c call lect_var(nid,'r',4,ierr,itau, Rnew)
669	c call lect_var(nid,'euv',4,ierr,itau, euv)
670	c call lect_var(nid,'cond',4,ierr,itau, con)
671	c call lect_var(nid,'nir',4,ierr,itau, nir)
672	do n=1,varspecsize
673	i=varspec(n)
674	call lect_var(nid,'n_'//trim(noms(i)),4,ierr,itau, spec)
675	do l=1,llm
676	do j=1,jjp1
677	do i=1,iip1
678	nspec(i,j,l,n)=spec(i,j,l)
679	enddo
680	enddo
681	enddo
682	enddo
683
684	if(fstats*ff-fichsize .eq. 0) then
685	call lect_var(nid,'temp_sd',4,ierr,itau, t_sd)
686	call lect_var(nid,'u_sd',4,ierr,itau, u_sd)
687	call lect_var(nid,'v_sd',4,ierr,itau, v_sd)
688	c call lect_var(nid,'w_sd',4,ierr,itau, w_sd)
689	c call lect_var(nid,'rho_sd',4,ierr,itau, rho_sd)
690	do n=1,varspecsize
691	i=varspec(n)
692	call lect_var(nid,'n_'//trim(noms(i))//'_sd',4,ierr,
693	& itau, spec)
694	do l=1,llm
695	do j=1,jjp1
696	do i=1,iip1
697	nspec_sd(i,j,l,n)=spec(i,j,l)
698	enddo
699	enddo
700	enddo
701	enddo
702	endif
703	! =================================================================
704
705	! COMPUTING ALTITUDE OF GCM LEVELS IN GCM GRID
706	! =================================================================
707	c zlocal ! altitude above the local surface (m)
708	c zareoid ! altitude above the mola areoid (m)
709
710	do j=1,jjp1
711	do i=1,iip1
712	sigij(1)=aps(1)/ps(i,j)+bps(1)
713	zlocal(i,j,1)=-log(sigij(1))* Rnew(i,j,1)*t(i,j,1)/g
714	c phis(i,j)=0.0
715	c phi(i,j,1)=phis(i,j)
716	c rho(i,j,1)=sigij(1)ps(i,j)/(Rnew(i,j,1)t(i,j,1))
717	zareoid(i,j,1) = zlocal(i,j,1) + phis(i,j)/g
718	do l=2,llm
719	sigij(l)=aps(l)/ps(i,j)+bps(l)
720	tmean=t(i,j,l)
721	if(t(i,j,l).ne.t(i,j,l-1))
722	& tmean=(t(i,j,l)-t(i,j,l-1))/log(t(i,j,l)/t(i,j,l-1))
723	zlocal(i,j,l)=zlocal(i,j,l-1)
724	& -log(sigij(l)/sigij(l-1))Rnew(i,j,l-1)tmean/g
725	zareoid(i,j,l) = zlocal(i,j,l) + phis(i,j)/g
726	c rho(i,j,l)=sigij(l)ps(i,j)/(Rnew(i,j,l)t(i,j,l))
727	phi(i,j,l)=zareoid(i,j,l)*g
728	enddo
729
730	c Altitude of the dummy thermosphere level (as in atmemcd)
731	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
732
733	tmean=(T_therm-t(i,j,llm))/log(T_therm/t(i,j,llm))
734	zlocal(i,j,llm+1)=zlocal(i,j,llm)
735	& -log(sig_therm/sigij(llm))Rnew(i,j,llm)tmean/g
736	zareoid(i,j,llm+1) = zlocal(i,j,llm+1) + phis(i,j)/g
737
738	enddo
739	enddo
740
741	c do l=1,llm
742	c print*,zareoid(32,24,l),zlocal(32,24,l),zaref(l)
743	c print*,zareoid(32,24,l),zaref(l)
744	c enddo
745
746	c Local time
747	c ----------
748	c universal time (in stats oe MCD files)
749	utime = (itau-1)*2.
750
751	c local time (in stats file)
752	DO i = 1,iim
753	localtime(i)=utime + 12.*rlonv(i)/pi
754	if(localtime(i).gt.24) localtime(i)=localtime(i)-24.
755	if(localtime(i).lt.0) localtime(i)=localtime(i)+24.
756	ENDDO
757
758	c c Passage en niveaux de pression
759	c ------------------------------
760	if(coor.ne.2) then ! PRESSURE COORDINATES
761	call xsig2p (ps,T,pref,llm,Tp)
762	call xsig2p (ps,u,pref,llm,up)
763	call xsig2p (ps,v,pref,llm,vp)
764	call xsig2p (ps,w,pref,llm,wp)
765	call xsig2p (ps,rho,pref,llm,rhop)
766	do iq=1,varspecsize
767	call xsig2p (ps,nspec(1,1,1,iq),pref,llm,nspecp(1,1,1,iq))
768	enddo
769	if(fstats*ff-fichsize .eq. 0) then
770	call xsig2p (ps,u_sd,pref,llm,up_sd)
771	call xsig2p (ps,v_sd,pref,llm,vp_sd)
772	call xsig2p(ps,t_sd,pref,llm,Tp_sd)
773	call xsig2p (ps,w_sd,pref,llm,wp_sd)
774	call xsig2p (ps,rho_sd,pref,llm,rhop_sd)
775	do iq=1,varspecsize
776	call xsig2p (ps,nspec_sd(1,1,1,iq),pref,llm,
777	$ nspecp_sd(1,1,1,iq))
778	enddo
779	endif
780	endif
781
782	c c Passage en niveaux de zareoid
783	c ------------------------------
784	if(coor.ne.1) then ! ZAEROID COORDINATES
785	do j=1,jjp1
786	do i=1,iip1
787	do l=1,llm
788	Tij(l)=T(i,j,l)
789	uij(l)=u(i,j,l)
790	vij(l)=v(i,j,l)
791	wij(l)=w(i,j,l)
792	rhoij(l)=rho(i,j,l)
793	do iq=1,varspecsize
794	nspecij(l,iq)=nspec(i,j,l,iq)
795	if(l.eq.llm) nspecij(l+1,iq)=nspecij(l,iq)*1.e-3
796	enddo
797	c euvij(l)=euv(i,j,l)
798	c conij(l)=con(i,j,l)
799	c nirij(l)=nir(i,j,l)
800	zaij(l)=zareoid(i,j,l)
801	p_zaij(l)=exp(-zareoid(i,j,l)/1.e4) ! used to interpolate rho
802	if(fstats*ff-fichsize .eq. 0) then
803	Tij_sd(l)=T_sd(i,j,l)
804	uij_sd(l)=u_sd(i,j,l)
805	vij_sd(l)=v_sd(i,j,l)
806	wij_sd(l)=w_sd(i,j,l)
807	rhoij_sd(l)=rho_sd(i,j,l)
808	do iq=1,varspecsize
809	nspecij_sd(l,iq)=nspec_sd(i,j,l,iq)
810	if(l.eq.llm) nspecij_sd(l+1,iq)=nspecij_sd(l,iq)*1.e-3
811	enddo
812	endif
813	enddo
814	p_zaij(l+1)=exp(-zareoid(i,j,l+1)/1.e4)
815	c rhoij(l+1)=sig_thermps(i,j)/(RT_therm)!in dummy thermosphere
816	rhoij(l+1)=rhoij(llm)*1.e-3
817	rhoij_sd(l+1)=rhoij_sd(llm)*1.e-3
818
819	do l=1,llm
820	call interpolf(zaref(l),y,zaij,Tij,llm)
821	Ta(i,j,l)=y
822	call interpolf(zaref(l),y,zaij,uij,llm)
823	ua(i,j,l)=y
824	call interpolf(zaref(l),y,zaij,vij,llm)
825	va(i,j,l)=y
826	call interpolf(zaref(l),y,zaij,wij,llm)
827	wa(i,j,l)=y
828	c call interpolf(zaref(l),y,zaij,euvij,llm)
829	c euva(i,j,l)=y
830	c call interpolf(zaref(l),y,zaij,conij,llm)
831	c cona(i,j,l)=y
832	c call interpolf(zaref(l),y,zaij,nirij,llm)
833	c nira(i,j,l)=y
834
835	p_zaref=exp(-zaref(l)/1.e4)
836	call interpolf(p_zaref,y,p_zaij,rhoij,llm+1)
837	rhoa(i,j,l)=y
838
839	do iq=1,varspecsize
840	call interpolf(p_zaref,y,p_zaij,nspecij(1,iq),llm+1)
841	c call interpolf(zaref(l),y,zaij,nspecij(1,iq),llm)
842	nspeca(i,j,l,iq)=y
843	enddo
844
845	if(fstats*ff-fichsize .eq. 0) then
846	call interpolf(zaref(l),y,zaij,Tij_sd,llm)
847	Ta_sd(i,j,l)=y
848	call interpolf(zaref(l),y,zaij,uij_sd,llm)
849	ua_sd(i,j,l)=y
850	call interpolf(zaref(l),y,zaij,vij_sd,llm)
851	va_sd(i,j,l)=y
852	call interpolf(zaref(l),y,zaij,wij_sd,llm)
853	wa_sd(i,j,l)=y
854	call interpolf(p_zaref,y,p_zaij,rhoij_sd,llm+1)
855	rhoa_sd(i,j,l)=y
856	do iq=1,varspecsize
857	call interpolf(p_zaref,y,p_zaij,nspecij_sd(1,iq),
858	& llm+1)
859	c call interpolf(zaref(l),y,zaij,nspecij_sd(1,iq),llm)
860	nspeca_sd(i,j,l,iq)=y
861	enddo
862	endif
863
864	enddo
865	enddo
866	enddo
867
868	c do l=1,llm
869	c print*,T(32,24,l),Ta(32,24,l)
870	c print*,nspec(32,24,l,1),nspeca(32,24,l,1)
871	c enddo
872
873	endif
874
875	c Output at each time-step, in netcdf format
876	! =================================================================
877	if (sor .ne. 2) then
878
879	ierr= NF_INQ_VARID(nout,"Time",nvarid)
880	#ifdef NC_DOUBLE
881	ierr= NF_PUT_VARA_DOUBLE(nout,nvarid,itau,1,temps(itau))
882	#else
883	ierr= NF_PUT_VARA_REAL(nout,nvarid,itau,1,temps(itau))
884	#endif
885	if (ierr.ne.NF_NOERR) then
886	write(,) "time problem ",NF_STRERROR(ierr)
887	stop
888	endif
889
890	if (varoutsize.eq.6) n=1
891	if (varoutsize.ne.6) n=varoutsize
892	do i=1,n
893	if(varout(i).eq.1 .or. varoutsize.eq.6) then
894	if(coor.ne.1)
895	& call put_var(nout,'temp','temp','K',4,ierr,itau,Ta)
896	if(coor.ne.2)
897	& call put_var(nout,'temp','temp','K',4,ierr,itau,Tp)
898	if(fstats*ff-fichsize.eq.0) then
899	if(coor.ne.1) call put_var(nout,'temp_sd','temp_sd','K',
900	& 4,ierr,itau,Ta_sd)
901	if(coor.ne.2) call put_var(nout,'temp_sd','temp_sd','K',
902	& 4,ierr,itau,Tp_sd)
903	endif
904	endif
905	if(varout(i).eq.2 .or. varoutsize.eq.6) then
906	if(coor.ne.1)
907	& call put_var(nout,'u','u','m s-1',4,ierr,itau,ua)
908	if(coor.ne.2)
909	& call put_var(nout,'u','u','m s-1',4,ierr,itau,up)
910	if(fstats*ff-fichsize.eq.0) then
911	if(coor.ne.1) call put_var(nout,'u_sd','u_sd','m s-1',
912	& 4,ierr,itau,ua_sd)
913	if(coor.ne.2) call put_var(nout,'u_sd','u_sd','m s-1',
914	& 4,ierr,itau,up_sd)
915	endif
916	endif
917	if(varout(i).eq.3 .or. varoutsize.eq.6) then
918	if(coor.ne.1)
919	& call put_var(nout,'v','v','m s-1',4,ierr,itau,va)
920	if(coor.ne.2)
921	& call put_var(nout,'v','v','m s-1',4,ierr,itau,vp)
922	if(fstats*ff-fichsize.eq.0) then
923	if(coor.ne.1) call put_var(nout,'v_sd','v_sd','m s-1',
924	& 4,ierr,itau,va_sd)
925	if(coor.ne.2) call put_var(nout,'v_sd','v_sd','m s-1',
926	& 4,ierr,itau,vp_sd)
927	endif
928	endif
929	if(varout(i).eq.4 .or. varoutsize.eq.6) then
930	if(coor.ne.1)
931	& call put_var(nout,'w','w','m s-1',4,ierr,itau,wa)
932	if(coor.ne.2)
933	& call put_var(nout,'w','w','m s-1',4,ierr,itau,wp)
934	if(fstats*ff-fichsize.eq.0) then
935	if(coor.ne.1) call put_var(nout,'w_sd','w_sd','m s-1',
936	& 4,ierr,itau,wa_sd)
937	if(coor.ne.2) call put_var(nout,'w_sd','w_sd','m s-1',
938	& 4,ierr,itau,wp_sd)
939	endif
940	endif
941	if(varout(i).eq.5 .or. varoutsize.eq.6) then
942	if(coor.ne.1) call put_var(nout,'rho','rho','',4,ierr,
943	& itau,rhoa)
944	if(coor.ne.2) call put_var(nout,'rho','rho','',4,ierr,
945	& itau,rhop)
946	if(fstats*ff-fichsize.eq.0) then
947	if(coor.ne.1) call put_var(nout,'rho_sd','rho_sd','',
948	& 4,ierr,itau,rhoa_sd)
949	if(coor.ne.2) call put_var(nout,'rho_sd','rho_sd','',
950	& 4,ierr,itau,rhop_sd)
951	endif
952	endif
953	if(varout(i).eq.6 .or. varoutsize.eq.6) then
954	do iq=1,varspecsize
955	j=varspec(iq)
956	if(coor.ne.1) call put_var(nout,'n_'//trim(noms(j)),
957	& 'n_'//trim(noms(j)),'cm-3',4,ierr,itau,nspeca(1,1,1,iq))
958	if(coor.ne.2) call put_var(nout,'n_'//trim(noms(j)),
959	& 'n_'//trim(noms(j)),'cm-3',4,ierr,itau,nspecp(1,1,1,iq))
960	if(fstats*ff-fichsize.eq.0) then
961	if(coor.ne.1) call put_var(nout,
962	& 'n_'//trim(noms(j))//'_sd',
963	& 'n_'//trim(noms(j))//'_sd','cm-3',4,ierr,
964	& itau,nspeca_sd(1,1,1,iq))
965	if(coor.ne.2) call put_var(nout,
966	& 'n_'//trim(noms(j))//'_sd',
967	& 'n_'//trim(noms(j))//'_sd','cm-3',4,ierr,
968	& itau,nspecp_sd(1,1,1,iq))
969	endif
970	enddo
971	endif
972	enddo
973	endif !endif writing data on Netcdf file
974
975	c Output at each time-step, ASCII file for IDL
976	! =================================================================
977	c In pressure and zareoid coordinates:
978
979	if (sor .ne. 1) then
980	do n=1,varoutsize
981	if(varout(n).eq.1) then
982	do l=1,llm
983	if(coor.ne.2) then
984	k=121+n
985	write(k,124) ((Tp(i,j,l),i=1,iim),j=1,jjp1)
986	elseif(coor.ne.1) then
987	k=1221+n
988	write(k,124) ((Ta(i,j,l),i=1,iim),j=1,jjp1)
989	endif
990	enddo
991	endif
992	if(varout(n).eq.2) then
993	do l=1,llm
994	if(coor.ne.2) then
995	k=121+n
996	write(k,124) ((up(i,j,l),i=1,iim),j=1,jjp1)
997	elseif(coor.ne.1)then
998	k=1221+n
999	write(k,124) ((ua(i,j,l),i=1,iim),j=1,jjp1)
1000	endif
1001	enddo
1002	endif
1003	if(varout(n).eq.3) then
1004	do l=1,llm
1005	if(coor.ne.2) then
1006	k=121+n
1007	write(k,124) ((vp(i,j,l),i=1,iim),j=1,jjp1)
1008	elseif(coor.ne.1) then
1009	k=1221+n
1010	write(k,124) ((va(i,j,l),i=1,iim),j=1,jjp1)
1011	endif
1012	enddo
1013	endif
1014	if(varout(n).eq.4) then
1015	do l=1,llm
1016	if(coor.ne.2) then
1017	k=121+n
1018	write(122,124) ((wp(i,j,l),i=1,iim),j=1,jjp1)
1019	elseif(coor.ne.1) then
1020	k=1221+n
1021	write(k,124) ((wa(i,j,l),i=1,iim),j=1,jjp1)
1022	endif
1023	enddo
1024	endif
1025	if(varout(n).eq.5) then
1026	do l=1,llm
1027	if(coor.ne.2) then
1028	k=121+n
1029	write(k,124) ((rhop(i,j,l),i=1,iim),j=1,jjp1)
1030	elseif(coor.ne.1) then
1031	k=1221+n
1032	write(k,124) ((rhoa(i,j,l),i=1,iim),j=1,jjp1)
1033	endif
1034	enddo
1035	endif
1036	if(varout(n).eq.6) then
1037	do iq=1,varspecsize
1038	do l=1,llm
1039	if(coor.ne.1) then
1040	k=1221+n+(iq-1)
1041	write(k,124) ((nspeca(i,j,l,iq),i=1,iim),j=1,jjp1)
1042	endif
1043	enddo
1044	enddo
1045	endif
1046	enddo
1047	c k=k+1
1048	c do l=1,llm
1049	c if(coor.ne.1)write(k,124) ((euva(i,j,l),i=1,iim),j=1,jjp1)
1050	c enddo
1051	c k=k+1
1052	c do l=1,llm
1053	c if(coor.ne.1)write(k,124) ((cona(i,j,l),i=1,iim),j=1,jjp1)
1054	c enddo
1055	c k=k+1
1056	c do l=1,llm
1057	c if(coor.ne.1)write(k,124) ((nira(i,j,l),i=1,iim),j=1,jjp1)
1058	c enddo
1059	endif
1060	124 format(12(e10.4,1x))
1061
1062
1063	ENDDO !LOOP in TIME
1064	! =================================================================
1065
1066	900 continue
1067
1068
1069	c ECRITURE FINALE si besoin
1070	c *************************
1071	ierr= NF_CLOSE(nid)
1072	ierr=NF_CLOSE(nout)
1073	close (33)
1074	close (122)
1075	close (1222)
1076
1077	ENDDO !LOOP in FILES
1078	! =================================================================
1079
1080	9999 PRINT*,'Fin '
1081	1000 format(a5,3x,i4,i3,x,a39)
1082	7777 FORMAT ('latitude/longitude',4f7.1)
1083
1084	END
1085
1086	c*******************************************************
1087
1088	subroutine xsig2p (ps,qsig,pref,npref,qp)
1089	IMPLICIT NONE
1090	c=======================================================================
1091	c
1092	c Francois Forget Avril 1996
1093	c
1094	c Cette subroutine calcule les champs interpole en niveaux de pression
1095	c different niveaux de pression pref
1096	c
1097	c=======================================================================
1098	c-----------------------------------------------------------------------
1099	c declarations:
1100	c -------------
1101	#include "dimensions.h"
1102	#include "paramet.h"
1103	#include "comgeom.h"
1104	#include "comvert.h"
1105	#include "comconst.h"
1106
1107	c
1108	c ARGUMENTS
1109	c ---------
1110
1111	c Inputs:
1112
1113	REAL qsig(iip1,jjp1,llm)
1114	REAL ps(iip1,jjp1)
1115	integer npref
1116	real pref(npref)
1117
1118	c outputs
1119	REAL qp(iip1,jjp1,npref)
1120
1121
1122	c Variables du modele
1123	c -------------------
1124
1125	REAL lnp(llm) , q(llm)
1126	REAL x,y
1127	INTEGER i,j,l ,n
1128
1129
1130	logical firstcall
1131	save firstcall
1132	data firstcall/.true./
1133
1134	#include "fxyprim.h"
1135	c**********************************************************************
1136	if (firstcall) then
1137	write(,) 'kappa = ', kappa
1138	write(,) 'cpp = ', cpp
1139	firstcall = .false.
1140	end if
1141
1142	DO j=1,jjp1
1143	DO i=1,iip1
1144	do l=1,llm
1145	lnp(l) = -log(aps(l)+bps(l)*ps(i,j))
1146	q(l) = qsig(i,j,l)
1147	end do
1148	do n =1,npref
1149	if ( (pref(n).le.ps(i,j)) .and.
1150	& (pref(n).ge.(ps(i,j)*bps(llm)+aps(llm))) ) then
1151	x = -log(pref(n))
1152	call interpolf(x,y,lnp,q,llm)
1153	qp(i,j,n) = y
1154	else
1155	qp(i,j,n) = 1E+20
1156	end if
1157	end do
1158	ENDDO
1159	ENDDO
1160
1161	return
1162	end
1163
1164	subroutine missing_value(nout,nvarid,valid_range,missing)
1165	IMPLICIT NONE
1166	c=======================================================================
1167	c
1168	c=======================================================================
1169	c-----------------------------------------------------------------------
1170	c declarations:
1171	c -------------
1172	include "netcdf.inc"
1173
1174	INTEGER nout,nvarid,ierr
1175	REAL missing
1176	REAL valid_range(2)
1177	c
1178	c ARGUMENTS
1179	c ---------
1180
1181	c Inputs:
1182
1183	ierr= NF_PUT_ATT_REAL(nout,nvarid,'valid_range',
1184	$ NF_FLOAT,2,valid_range)
1185	IF (ierr.NE.NF_NOERR) THEN
1186	PRINT*, 'anl_NC: valid_range attribution failed'
1187	WRITE(,) 'NF_NOERR', NF_NOERR
1188	CALL abort
1189	ENDIF
1190
1191	#ifdef NC_DOUBLE
1192	ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',
1193	$ NF_DOUBLE,1,missing)
1194	#else
1195	ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',
1196	$ NF_FLOAT,1,missing)
1197
1198	#endif
1199
1200	IF (ierr.NE.NF_NOERR) THEN
1201	PRINT*, 'anl_NC: missing value attribution failed'
1202	WRITE(,) 'NF_NOERR', NF_NOERR
1203	CALL abort
1204	ENDIF
1205	return
1206	end
1207	**********************************************************
1208	subroutine lect_var(nid,name,nbdim,ierr,itau,
1209	$ var)
1210	IMPLICIT NONE
1211	c=======================================================================
1212	c
1213	c=======================================================================
1214	c-----------------------------------------------------------------------
1215	c declarations:
1216	c -------------
1217	#include "netcdf.inc"
1218	#include "dimensions.h"
1219	#include "paramet.h"
1220	#include "comgeom.h"
1221	#include "comvert.h"
1222	#include "comconst.h"
1223
1224	c inputs
1225	character (len=*) :: name
1226	INTEGER nid,nbdim,nvarid,ierr,itau
1227	integer, dimension(nbdim) :: edges,corner
1228	c output var
1229	integer, dimension(nbdim) :: var
1230	c ---------
1231	if (nbdim.eq.4) then
1232	corner=(/1,1,1,itau/)
1233	edges=(/iip1,jjp1,llm,1/)
1234	endif
1235
1236	if (nbdim.eq.3) then
1237	corner=(/1,1,itau/)
1238	edges=(/iip1,jjp1,1/)
1239	endif
1240
1241
1242	ierr = NF_INQ_VARID (nid,adjustl(name), nvarid)
1243	#ifdef NC_DOUBLE
1244	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,corner,edges, var)
1245	#else
1246	ierr = NF_GET_VARA_REAL(nid, nvarid,corner,edges, var)
1247	#endif
1248	IF (ierr.NE.NF_NOERR) THEN
1249	write(,) 'anl_NC: reading failed for ', name
1250	CALL abort
1251	ENDIF
1252
1253	return
1254	end
1255
1256
1257	subroutine put_var(nout,name,title,units,nbdim,ierr,itau,
1258	$ var)
1259	IMPLICIT NONE
1260	c=======================================================================
1261	c
1262	c=======================================================================
1263	c-----------------------------------------------------------------------
1264	c declarations:
1265	c -------------
1266	#include "netcdf.inc"
1267	#include "dimensions.h"
1268	#include "paramet.h"
1269	#include "comgeom.h"
1270	#include "comvert.h"
1271	#include "comconst.h"
1272
1273	character (len=*) :: title,units,name
1274	INTEGER nout,nbdim,nvarid,ierr,itau
1275	integer, dimension(nbdim) :: edges,corner,id,var
1276	REAL valid_range(2)
1277	DATA valid_range /0., 300.0/
1278	c ---------
1279	if (nbdim.eq.4) then
1280	ierr= NF_INQ_DIMID(nout,"longitude",id(1))
1281	ierr= NF_INQ_DIMID(nout,"latitude",id(2))
1282	ierr= NF_INQ_DIMID(nout,"altitude",id(3))
1283	ierr= NF_INQ_DIMID(nout,"Time",id(4))
1284	corner=(/1,1,1,itau/)
1285	edges=(/iip1,jjp1,llm,1/)
1286	endif
1287
1288	if (nbdim.eq.3) then
1289	ierr= NF_INQ_DIMID(nout,"longitude",id(1))
1290	ierr= NF_INQ_DIMID(nout,"latitude",id(2))
1291	ierr= NF_INQ_DIMID(nout,"Time",id(3))
1292	corner=(/1,1,itau/)
1293	edges=(/iip1,jjp1,1/)
1294	endif
1295
1296
1297	ierr = NF_INQ_VARID (nout,adjustl(name), nvarid)
1298	if (ierr /= NF_NOERR) then
1299	! choix du nom des coordonnees
1300	! Creation de la variable si elle n'existait pas
1301	write (,) "====================="
1302	write (,) "creation de ",name
1303	call def_var(nout,adjustl(title),adjustl(name),adjustl(units),
1304	$ nbdim,id,nvarid,ierr)
1305	if (name.eq.'temp') then
1306	ierr = NF_REDEF (nout)
1307	call missing_value(nout,nvarid,valid_range,1E+20)
1308	ierr = NF_ENDDEF(nout)
1309	endif
1310	endif
1311	#ifdef NC_DOUBLE
1312	ierr = NF_PUT_VARA_DOUBLE(nout,nvarid,corner,edges,var)
1313	#else
1314	ierr = NF_PUT_VARA_REAL(nout,nvarid,corner,edges,var)
1315	#endif
1316	IF (ierr.NE.NF_NOERR) THEN
1317	write(,) 'anl_NC: writing failed for ',name
1318	CALL abort
1319	ENDIF
1320
1321
1322	return
1323	end
1324
1325
1326	Subroutine interpolf(x,y,xd,yd,nd)
1327
1328	c******************************************************
1329	c SUBROUTINE (interpol)
1330	c interpolation, give y = f(x) with array xd,yd known, size nd
1331
1332	c Version with CONSTANT values oustide limits
1333	c**********************************************************
1334
1335
1336	c Variable declaration
1337	c --------------------
1338	c Arguments :
1339	real x,y
1340	real xd(),yd()
1341	integer nd
1342	c internal
1343	integer i,j
1344	real y_undefined
1345
1346	c run
1347	c ---
1348	c y_undefined=-999.999
1349	y_undefined=1.e20
1350
1351	y=0.
1352	if ((x.le.xd(1)).and.(x.le.xd(nd))) then
1353	if (xd(1).lt.xd(nd)) y = y_undefined ! yd(1)
1354	if (xd(1).ge.xd(nd)) y = y_undefined ! yd(nd)
1355	else if ((x.ge.xd(1)).and.(x.ge.xd(nd))) then
1356	if (xd(1).lt.xd(nd)) y = y_undefined ! yd(nd)
1357	if (xd(1).ge.xd(nd)) y = y_undefined ! yd(1)
1358	c y = yd(nd)
1359	else
1360	do i=1,nd-1
1361	if ( ( (x.ge.xd(i)).and.(x.lt.xd(i+1)) )
1362	& .or. ( (x.le.xd(i)).and.(x.gt.xd(i+1)) ) ) then
1363	y=yd(i)+(x-xd(i))*(yd(i+1)-yd(i))/(xd(i+1)-xd(i))
1364	goto 99
1365	end if
1366	end do
1367	end if
1368
1369	c write (,)' x , y' , x,y
1370	c do i=1,nd
1371	c write (,)' i, xd , yd' , xd(i),yd(i)
1372	c end do
1373	c stop
1374
1375	99 continue
1376
1377	end

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

Download in other formats:

Original Format