Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

anl_mcdstats.F @ 259

Last change on this file since 259 was 38, checked in by emillour, 14 years ago
Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs. EM
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	& NF_CLOBBER, 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	& NF_CLOBBER, 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