Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

xvik.F @ 57

Last change on this file since 57 was 57, checked in by aslmd, 14 years ago
mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)
File size: 15.3 KB

Rev	Line
[57]	1	PROGRAM xvik
	2	IMPLICIT NONE
	3	c=======================================================================
	4	c
	5	c Pression au site Viking
	6	c
	7	c=======================================================================
	8	c-----------------------------------------------------------------------
	9	c declarations:
	10	c -------------
	11
	12
	13	#include "dimensions.h"
	14	#include "paramet.h"
	15	#include "comconst.h"
	16	#include "comdissip.h"
	17	#include "comvert.h"
	18	#include "comgeom2.h"
	19	#include "logic.h"
	20	#include "temps.h"
	21	#include "control.h"
	22	#include "ener.h"
	23	#include "description.h"
	24	#include "netcdf.inc"
	25
	26
	27	INTEGER itau,nbpas,nbpasmx
	28	PARAMETER(nbpasmx=1000000)
	29	REAL temps(nbpasmx)
	30	INTEGER unitlec
	31	INTEGER i,j,l,jj
	32	REAL constR
	33
	34	c Declarations NCDF:
	35	c -----------------
	36	CHARACTER*100 varname
	37	INTEGER ierr,nid,nvarid,dimid
	38	LOGICAL nc
	39	INTEGER start_ps(3),start_temp(4),start_co2ice(3)
	40	INTEGER count_ps(3),count_temp(4),count_co2ice(3)
	41
	42	c declarations pour les points viking:
	43	c ------------------------------------
	44	INTEGER ivik(2),jvik(2),ifile(2),iv
	45	REAL lonvik(2),latvik(2),phivik(2),phisim(2)
	46	REAL unanj
	47
	48	c variables meteo:
	49	c ----------------
	50	REAL vnat(iip1,jjm,llm),unat(iip1,jjp1,llm)
	51	REAL t(iip1,jjp1,llm),ps(iip1,jjp1),pstot, phis(iip1,jjp1)
	52	REAL co2ice(iip1,jjp1), captotN,captotS
	53
	54	REAL zp1,zp2,zp2_sm,zu,zv,zw(0:1,0:1,2),zalpha,zbeta
	55
	56	LOGICAL firstcal,lcal,latcal,lvent,day_ls
	57	INTEGER*4 day0
	58
	59	REAL ziceco2(iip1,jjp1)
	60	REAL day,zt,sollong,sol,dayw
	61	REAL airtot1,gh
	62
	63	INTEGER ii,iyear,kyear
	64
	65	CHARACTER*1 chr2
	66
	67
	68	c declarations de l'interface avec mywrite:
	69	c -----------------------------------------
	70
	71	CHARACTER file*80
	72	CHARACTER pathchmp80,pathsor80,nomfich*80
	73
	74	c externe:
	75	c --------
	76
	77	EXTERNAL iniconst,inigeom,covcont,mywrite
	78	EXTERNAL inifilr,exner,pbar
	79	EXTERNAL solarlong,coordij,moy2
	80	EXTERNAL SSUM
	81	REAL SSUM
	82	EXTERNAL lnblnk
	83	INTEGER lnblnk
	84
	85	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	86
	87	c-----------------------------------------------------------------------
	88	c initialisations:
	89	c ----------------
	90
	91	unanj=667.9
	92	print*,'WARNING!!!',unanj,'Jours/an'
	93	nc=.true.
	94	lcal=.true.
	95	latcal=.true.
	96	lvent=.false.
	97	day_ls=.true.
	98
	99	c lecture du fichier xvik.def
	100
	101	phivik(1)=-3627
	102	phivik(2)=-4505
	103
	104
	105
	106	OPEN(99,file='xvik.def',form='formatted')
	107
	108	READ(99,*)
	109	READ(99,*,iostat=ierr) phivik
	110	IF(ierr.NE.0) GOTO 105
	111
	112	READ(99,*,END=105)
	113	READ(99,'(a)',END=105) pathchmp
	114	READ(99,*,END=105)
	115	READ(99,'(a)',END=105) pathsor
	116	READ(99,*,END=105)
	117	c READ(99,'(l1)',END=105) day_ls
	118	READ(99,'(l1)',END=105)
	119	READ(99,'(l1)',END=105) lcal
	120	READ(99,'(l1)',END=105)
	121	READ(99,'(l1)',END=105) lvent
	122	READ(99,'(l1)',END=105)
	123	READ(99,'(l1)',END=105) latcal
	124
	125	105 CONTINUE
	126	CLOSE(99)
	127	write (,)'>>>>>>>>>>>>>>>>', phivik,g
	128	DO iv=1,2
	129	phivik(iv)=phivik(iv)*3.73
	130	END DO
	131
	132
	133	c-----------------------------------------------------------------------
	134	c-----------------------------------------------------------------------
	135	c ouverture des fichiers xgraph:
	136	c ------------------------------
	137
	138	ifile(1)=12
	139	ifile(2)=13
	140	kyear=-1
	141	c OPEN(77,file='xlongday',form='formatted')
	142
	143	unitlec=11
	144	c
	145	PRINT*,'entrer le nom du fichier NC'
	146	READ(5,'(a)') nomfich
	147
	148	PRINT *,'nomfich : ',nomfich
	149
	150
	151	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	152	c grande boucle sur les fichiers histoire:
	153	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	154
	155	firstcal=.true.
	156	DO WHILE(lnblnk(nomfich).GT.0.AND.lnblnk(nomfich).LT.50)
	157	PRINT *,'nomfich : ',nomfich
	158
	159	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	160
	161	file=pathchmp(1:lnblnk(pathchmp))//'/'//
	162	s nomfich(1:lnblnk(nomfich))
	163	PRINT*,'file.nc: ', file(1:lnblnk(file))//'.nc'
	164	PRINT*,'timestep ',dtvr
	165
	166	IF(nc) THEN
	167	ierr= NF_OPEN(file(1:lnblnk(file))//'.nc',NF_NOWRITE,nid)
	168	ELSE
	169	PRINT*,'Ouverture binaire ',file
	170	OPEN(unitlec,file=file,status='old',form='unformatted',
	171	. iostat=ierr)
	172	ENDIF
	173
	174	c----------------------------------------------------------------------
	175	c initialisation de la physique:
	176	c ------------------------------
	177
	178	CALL readhead_NC(file(1:lnblnk(file))//'.nc',day0,phis,constR)
	179
	180	WRITE (,) 'day0 = ' , day0
	181
	182	CALL iniconst
	183	CALL inigeom
	184	CALL inifilr
	185
	186
	187	c Lecture temps :
	188
	189	ierr= NF_INQ_DIMID (nid,"Time",dimid)
	190	IF (ierr.NE.NF_NOERR) THEN
	191	PRINT*, 'xvik: Le champ <Time> est absent'
	192	CALL abort
	193	ENDIF
	194
	195	ierr= NF_INQ_DIMLEN (nid,dimid,nbpas)
	196
	197	ierr = NF_INQ_VARID (nid, "Time", nvarid)
	198	#ifdef NC_DOUBLE
	199	ierr = NF_GET_VAR_DOUBLE(nid, nvarid, temps)
	200	#else
	201	ierr = NF_GET_VAR_REAL(nid, nvarid, temps)
	202	#endif
	203	IF (ierr.NE.NF_NOERR) THEN
	204	PRINT*, 'xvik: Lecture echouee pour <Time>'
	205	CALL abort
	206	ENDIF
	207
	208	PRINT*,'temps',(temps(itau),itau=1,10)
	209
	210	c-----------------------------------------------------------------------
	211	c coordonnees des point Viking:
	212	c -----------------------------
	213
	214	latvik(1)=22.27*pi/180.
	215	lonvik(1)=-47.9*pi/180.
	216	latvik(2)=47.67*pi/180.
	217	lonvik(2)=(360.-225.71)*pi/180.
	218
	219	c ponderations pour les 4 points autour de Viking
	220	DO iv=1,2
	221	CALL coordij(lonvik(iv),latvik(iv),ivik(iv),jvik(iv))
	222	IF(lonvik(iv).lt.rlonv(ivik(iv))) THEN
	223	ivik(iv)=ivik(iv)-1
	224	ENDIF
	225	IF(latvik(iv).gt.rlatu(jvik(iv))) THEN
	226	jvik(iv)=jvik(iv)-1
	227	ENDIF
	228	zalpha=(lonvik(iv)-rlonv(ivik(iv)))/
	229	s (rlonv(ivik(iv)+1)-rlonv(ivik(iv)))
	230	zbeta=(latvik(iv)-rlatu(jvik(iv)))/
	231	s (rlatu(jvik(iv)+1)-rlatu(jvik(iv)))
	232	zw(0,0,iv)=(1.-zalpha)*(1.-zbeta)
	233	zw(1,0,iv)=zalpha*(1.-zbeta)
	234	zw(0,1,iv)=(1.-zalpha)*zbeta
	235	zw(1,1,iv)=zalpha*zbeta
	236	ENDDO
	237
	238	c altitude reelle et modele aux points Viking
	239	DO iv=1,2
	240	phisim(iv)=0.
	241	DO jj=0,1
	242	j=jvik(iv)+jj
	243	DO ii=0,1
	244	i=ivik(iv)+ii
	245	phisim(iv)=phisim(iv)+zw(ii,jj,iv)*phis(i,j)
	246	ENDDO
	247	ENDDO
	248	ENDDO
	249	PRINT*,'relief aux points Viking pour les sorties:',phivik
	250
	251	c----------------------------------------------------------------------
	252	c lectures des etats:
	253	c -------------------
	254
	255	airtot1=1./(SSUM(ip1jmp1,aire,1)-SSUM(jjp1,aire,iip1))
	256
	257	c======================================================================
	258	c debut de la boucle sur les etats dans un fichier histoire:
	259	c======================================================================
	260	count_ps=(/iip1,jjp1,1/)
	261	count_co2ice=(/iip1,jjp1,1/)
	262	count_temp=(/iip1,jjp1,llm,1/)
	263
	264	DO itau=1,nbpas
	265
	266	start_ps=(/1,1,itau/)
	267	start_co2ice=(/1,1,itau/)
	268	start_temp=(/1,1,1,itau/)
	269	c lecture drs des champs:
	270	c -----------------------
	271	c varname='u'
	272	c ierr=drsread (unitlec,varname,unat,itau)
	273	c PRINT*,'unat',unat(iip1/2,jjp1/2,llm/2)
	274	c varname='v'
	275	c ierr=drsread (unitlec,varname,vnat,itau)
	276	c PRINT*,'vnat',vnat(iip1/2,jjp1/2,llm/2)
	277
	278	ccccccccc LECTURE Ps ccccccccccccccccccccccccccc
	279	ierr = NF_INQ_VARID (nid, "ps", nvarid)
	280	#ifdef NC_DOUBLE
	281	ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start_ps,count_ps, ps)
	282	#else
	283	ierr = NF_GET_VARA_REAL(nid, nvarid,start_ps,count_ps, ps)
	284	#endif
	285	IF (ierr.NE.NF_NOERR) THEN
	286	PRINT*, 'xvik: Lecture echouee pour <ps>'
	287	CALL abort
	288	ENDIF
	289
	290	PRINT*,'ps',ps(iip1/2,jjp1/2)
	291
	292	ccccccccc LECTURE Temperature ccccccccccccccccccccccccccc
	293	ierr = NF_INQ_VARID (nid, "temp", nvarid)
	294	#ifdef NC_DOUBLE
	295	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_temp,count_temp, t)
	296	#else
	297	ierr = NF_GET_VARA_REAL(nid,nvarid,start_temp,count_temp, t)
	298	#endif
	299	IF (ierr.NE.NF_NOERR) THEN
	300	PRINT*, 'xvik: Lecture echouee pour <temp>'
	301	CALL abort
	302	ENDIF
	303
	304	c PRINT*,'t',t(iip1/2,jjp1/2,llm/2)
	305
	306	ccccccccc LECTURE co2ice ccccccccccccccccccccccccccc
	307	ierr = NF_INQ_VARID (nid, "co2ice", nvarid)
	308	#ifdef NC_DOUBLE
	309	ierr = NF_GET_VARA_DOUBLE(nid,nvarid,start_co2ice,
	310	& count_co2ice, co2ice)
	311	#else
	312	ierr = NF_GET_VARA_REAL(nid, nvarid,start_co2ice,
	313	& count_co2ice, co2ice)
	314	#endif
	315	IF (ierr.NE.NF_NOERR) THEN
	316	PRINT*, 'xvik: Lecture echouee pour <co2ice>'
	317	CALL abort
	318	ENDIF
	319
	320
	321	c Gestion du temps
	322	c ----------------
	323	day=temps(itau)
	324	PRINT*,'day ',day
	325	CALL solarlong(day+day0,sollong)
	326	sol=day+day0+461.
	327	iyear=sol/unanj
	328	WRITE (,) 'iyear',iyear
	329	sol=sol-iyear*unanj
	330	c
	331	c Ouverture / fermeture des fichiers
	332	c ----------------------------------
	333	IF (iyear.NE.kyear) THEN
	334	WRITE(chr2(1:1),'(i1)') iyear+1
	335	WRITE (,) 'iyear bis',iyear
	336	WRITE (,) 'chr2'
	337	WRITE (,) chr2
	338	IF(iyear.GE.9) WRITE(chr2,'(i2)') iyear+1
	339	kyear=iyear
	340	DO ii=1,2
	341	CLOSE(10+ifile(ii))
	342	CLOSE(2+ifile(ii))
	343	CLOSE(4+ifile(ii))
	344	CLOSE(6+ifile(ii))
	345	CLOSE(8+ifile(ii))
	346	CLOSE(16+ifile(ii))
	347	CLOSE(12+ifile(ii))
	348	CLOSE(14+ifile(ii))
	349	CLOSE(97)
	350	CLOSE(98)
	351	ENDDO
	352	CLOSE(5+ifile(1))
	353	OPEN(ifile(1)+10,file='xpsol1'//chr2,form='formatted')
	354	OPEN(ifile(2)+10,file='xpsol2'//chr2,form='formatted')
	355	c OPEN(ifile(1)+8,file='xbpsol1'//chr2,form='formatted')
	356	c OPEN(ifile(2)+8,file='xbpsol2'//chr2,form='formatted')
	357	c OPEN(ifile(1)+2,file='xlps1'//chr2,form='formatted')
	358	c OPEN(ifile(2)+2,file='xlps2'//chr2,form='formatted')
	359	IF(lcal) THEN
	360	c OPEN(ifile(2)+4,file='xpressud'//chr2,form='formatted')
	361	c OPEN(ifile(1)+4,file='xpresnord'//chr2,form='formatted')
	362	c OPEN(ifile(1)+6,file='xpm2'//chr2,form='formatted')
	363	ENDIF
	364	IF(latcal) THEN
	365	c OPEN(ifile(2)+14,file='xlats'//chr2,form='formatted')
	366	c OPEN(ifile(1)+14,file='xlatn'//chr2,form='formatted')
	367	ENDIF
	368	IF(lvent) THEN
	369	c OPEN(ifile(1)+16,file='xu1'//chr2,form='formatted')
	370	c OPEN(ifile(2)+16,file='xu2'//chr2,form='formatted')
	371	c OPEN(ifile(1)+12,file='xv1'//chr2,form='formatted')
	372	c OPEN(ifile(2)+12,file='xv2'//chr2,form='formatted')
	373	ENDIF
	374	OPEN(97,file='xprestot'//chr2,form='formatted')
	375	c OPEN(98,file='xlat37_'//chr2,form='formatted')
	376	WRITE(98,'(f5.1,16f7.1)') 0.,(rlonv(i)*180./pi,i=1,iim,4)
	377	ENDIF
	378
	379
	380	sollong=sollong*180./pi
	381	IF(day_ls) THEN
	382	dayw=sol
	383	write(,) 'dayw', dayw
	384	ELSE
	385	dayw=sollong
	386	ENDIF
	387
	388	c Calcul de la moyenne planetaire
	389	c -------------------------------
	390	pstot=0.
	391	captotS=0.
	392	captotN=0.
	393	DO j=1,jjp1
	394	DO i=1,iim
	395	pstot=pstot+aire(i,j)*ps(i,j)
	396	ENDDO
	397	ENDDO
	398
	399	DO j=1,jjp1/2
	400	DO i=1,iim
	401	captotN = captotN +aire(i,j)*co2ice(i,j)
	402	ENDDO
	403	ENDDO
	404	DO j=jjp1/2+1, jjp1
	405	DO i=1,iim
	406	captotS = captotS +aire(i,j)*co2ice(i,j)
	407	ENDDO
	408	ENDDO
	409	WRITE(97,'(4e16.6)') dayw,pstot*airtot1
	410	& , captotNgairtot1, captotSgairtot1
	411
	412	IF(.NOT.firstcal) THEN
	413	WRITE(98,'(f5.1,16f7.3)')
	414	s dayw,(ps(i,37),i=1,iim,4)
	415
	416	c boucle sur les sites vikings:
	417	c ----------------------------
	418
	419	DO iv=1,2
	420
	421	c interpolation de la temperature dans la 7eme couche, de la pression
	422	c de surface et des vents aux points viking.
	423
	424	zp1=0.
	425	zp2=0.
	426	zp2_sm=0.
	427	zt=0.
	428	zu=0.
	429	zv=0.
	430	DO jj=0,1
	431	j=jvik(iv)+jj
	432	DO ii=0,1
	433	i=ivik(iv)+ii
	434	zt=zt+zw(ii,jj,iv)*t(i,j,7)
	435	zp1=zp1+zw(ii,jj,iv)*ps(i,j)
	436	WRITE (,) 'ps autour iv',ps(i,j),iv
	437	zu=zu+zw(ii,jj,iv)*unat(i,j,1)/cu(i,j)
	438	zv=zv+zw(ii,jj,iv)*vnat(i,j,1)/cv(i,j)
	439	ENDDO
	440	ENDDO
	441
	442
	443	c pression au sol extrapolee a partir de la temp. 7eme couche
	444	WRITE (,) 'constR ',constR
	445	WRITE (,) 'zt ',zt
	446	gh=constR*zt
	447	zp2=zp1*exp(-(phivik(iv)-phisim(iv))/gh)
	448
	449	WRITE (,) 'iv,pstot,zp2, zp1, phivik(iv),phisim(iv),gh'
	450	WRITE (,) iv,pstot*airtot1,zp2,zp1,phivik(iv),phisim(iv),gh
	451	WRITE(ifile(iv)+10,'(2e15.5)') dayw,zp1
	452
	453
	454	c sorties eventuelles de vent
	455	IF(lvent) THEN
	456	WRITE(ifile(iv)+16,'(2e15.5)')
	457	s dayw,zu
	458	WRITE(ifile(iv)+12,'(2e15.5)')
	459	s dayw,zv
	460	ENDIF
	461	ENDDO
	462	c IF (lcal) THEN
	463	c WRITE(ifile(1)+4,'(2e15.6)') dayw,airtot1g.01*
	464	c s (SSUM(ip1jmp1/2,ziceco2,1)-SSUM(jjp1/2,ziceco2,iip1))
	465	c WRITE(ifile(2)+4,'(2e15.6)') dayw,airtot1g.01*
	466	c s (SSUM(iip1*jjm/2,ziceco2(1,jjm/2+2),1)-
	467	c s SSUM(jjm/2,ziceco2(1,jjm/2+2),iip1))
	468	c ENDIF
	469	c IF(latcal) THEN
	470	c CALL icelat(iim,jjm,ziceco2,rlatv,zicelat)
	471	c WRITE(ifile(1)+14,'(2e15.6)') dayw,zicelat(1)*180./pi
	472	c WRITE(ifile(2)+14,'(2e15.6)') dayw,zicelat(2)*180./pi
	473	c ENDIF
	474	ENDIF
	475	firstcal=.false.
	476
	477	c======================================================================
	478	c Fin de la boucle sur les etats du fichier histoire:
	479	c======================================================================
	480	ENDDO
	481
	482	ierr= NF_CLOSE(nid)
	483
	484	PRINT*,'Fin du fichier',nomfich
	485	print*,'Entrer un nouveau fichier ou return pour finir'
	486	READ(5,'(a)',err=9999) nomfich
	487
	488	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	489	c Fin de la boucle sur les fichiers histoire:
	490	c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	491
	492	ENDDO
	493
	494	PRINT,'relief du point V1',.001phis(ivik(1),jvik(1))/g
	495	PRINT,'relief du point V2',.001phis(ivik(2),jvik(2))/g
	496	DO iv=1,2
	497	PRINT*,'Viking',iv,' i=',ivik(iv),'j =',jvik(iv)
	498	WRITE(6,7777)
	499	s (rlonv(i)*180./pi,i=ivik(iv)-1,ivik(iv)+2)
	500	print*
	501	DO j=jvik(iv)-1,jvik(iv)+2
	502	WRITE(6,'(f8.1,10x,5f7.1)')
	503	s rlatu(j)180./pi,(phis(i,j)/(g1000.),i=ivik(iv)-1,ivik(iv)+2)
	504	ENDDO
	505	print*
	506	print*,'zw'
	507	write(6,'(2(2f10.4/))') ((zw(ii,jj,iv),ii=0,1),jj=0,1)
	508	print*,'altitude interpolee (km) ',phisim(iv)/1000./g
	509	ENDDO
	510	PRINT*,'R=',r
	511	9999 PRINT*,'Fin '
	512
	513	7777 FORMAT ('latitude/longitude',4f7.1)
	514	END

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

Download in other formats:

Original Format