Context Navigation

anl_mcdstats.F @ 664

Last change on this file since 664 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

Rev	Line
[38]	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",
[410]	414	& IOR(NF_CLOBBER,NF_64BIT_OFFSET), nout)
[38]	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",
[410]	433	& IOR(NF_CLOBBER,NF_64BIT_OFFSET), nout)
[38]	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