Context Navigation

anl_stats-diag.F @ 937

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

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

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

Download in other formats:

Original Format