Context Navigation

anl_mcdstats.F @ 620

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

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

Download in other formats:

Original Format