Context Navigation

offlinenc.F @ 331

Last change on this file since 331 was 220, checked in by lmdz, 23 years ago
Erreur dans les parametres d'appel a phytrac (P.Bousquet) LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.5 KB

Rev	Line
[204]	1	c
[207]	2	c $Header$
[204]	3	c
[187]	4	PROGRAM offlinenc
	5	USE ioipsl
	6	IMPLICIT NONE
	7
	8	#include "dimensions.h"
	9	#include "paramet.h"
	10	#include "comconst.h"
	11	#include "comdissnew.h"
	12	#include "comvert.h"
	13	#include "comgeom2.h"
	14	#include "logic.h"
	15	#include "temps.h"
	16	#include "control.h"
	17	#include "ener.h"
	18	#include "description.h"
	19	#include "serre.h"
	20	#include "tracstoke.h"
	21	integer ngridmx,nbsrf
	22	PARAMETER (ngridmx=iim*(jjm-1)+2,nbsrf=4)
	23	c-----------------------------------------------------------------------
	24
	25	INTEGER*4 iday ! jour julien
	26	REAL time ! Heure de la journee en fraction d'1 jour
	27	REAL dtdyn,dtav
	28	INTEGER iday_step
	29
	30	INTEGER longcles
	31	PARAMETER ( longcles = 20 )
	32	REAL clesphy0( longcles )
	33	logical clefinit
	34
	35	c variables dynamiques
	36	REAL vcov(iip1,jjm,llm),ucov(iip1,jjp1,llm) ! vents covariants
	37	REAL q(iip1,jjp1,llm,nqmx) ! champs advectes
	38	real qmoy(iip1,jjp1,nqmx)
	39	REAL pks(iip1,jjp1) ! exner (f pour filtre)
	40	REAL phis(iip1,jjp1) ! geopotentiel au sol
	41	REAL phi(iip1,jjp1,llm) ! geopotentiel
	42	REAL w(iip1,jjp1,llm) ! vitesse verticale
	43	C
	44	C variables physiques
	45	real zmfd(ngridmx,llm),zde_d(ngridmx,llm),zen_d(ngridmx,llm)
	46	REAL qfi(ngridmx,llm,nqmx) ! champs advectes
	47	real zmfu(ngridmx,llm),zde_u(ngridmx,llm),zen_u(ngridmx,llm)
	48	real t(ngridmx,llm)
	49	real coefkz(ngridmx,llm)
	50	real frac_impa(ngridmx,llm),frac_nucl(ngridmx,llm)
	51	real pphis (ngridmx)
	52	REAL airefi(ngridmx)
	53	REAL latfi(ngridmx),lonfi(ngridmx) ! latitude et long pour chaque point
	54	real yu1(ngridmx), yv1(ngridmx)
	55	real ftsol(ngridmx,nbsrf),pctsrf(ngridmx,nbsrf)
	56	C
	57	character*10 file
	58	c variables dynamiques intermediaire pour le transport
	59	REAL pbaru(iip1,jjp1,llm),pbarv(iip1,jjm,llm) !flux de masse
	60	REAL masse(iip1,jjp1,llm)
	61	REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
	62	REAL teta(iip1,jjp1,llm)
	63	INTEGER itau,irec,nrec,ldec,itauav
	64	integer recmin,recmax,recstkmin,recstkmax
	65	character*1 injecte,stoke
	66	character*10 nom
	67
	68	real tinj2,tinj1
	69	integer iiinj
	70	integer jour0,isplit,nsplit_dyn,nsplit_phy,nsplit
[204]	71	logical debut,lectstart,rnpb,lafin
[187]	72
	73	EXTERNAL inidissip,iniconst,inifilr
	74	EXTERNAL inigeom
	75	EXTERNAL exner_hyb,pression
	76	EXTERNAL defrun_new
	77
	78	REAL time_0 ! temps de debut de simulation (% journee)
	79
	80	character*2 str2
	81	INTEGER iq,j,i,l,nq,mode
	82	real paprs(ngridmx,llm+1),pplay(ngridmx,llm),tempfi(ngridmx,llm)
	83	real p(iip1,jjp1,llmp1),pk(iip1,jjp1,llm)
	84	real lon_stat(nqmx),lat_stat(nqmx),zlon,zlat
	85	integer i_stat(nqmx),j_stat(nqmx)
	86	real intensite(nqmx)
	87	integer iqmin,iqmax
	88	integer npasinject
	89	data npasinject/1/
	90	integer histid, histvid, histaveid
	91	real time_step, t_wrt, t_ops
	92	integer itime_step
	93	character*80 dynhist_file, dynhistave_file
	94
	95	logical avant
	96	logical debutphy,redecoupage
	97	DATA debutphy/.true./
	98
	99	real unskap,pksurcp,smass(iip1,jjp1)
	100	integer ig0,ii,jj,iii,ij
	101
	102	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
	103	REAL ps(iip1,jjp1),pkf(ip1jmp1,llm)
	104
	105	integer ntau, iinjec2,iij
	106	real zttt
	107	c-----------------------------------------------------------------------
	108	c Initialisations:
	109	c ----------------
	110	descript = 'Run GCM LMDZ'
	111	tinj1=0.
	112	tinj2=0.
	113	rad = 6400000
	114	omeg = 7.272205e-05
	115	g = 9.8
	116	kappa = 0.285716
	117	daysec = 86400
	118	cpp = 1004.70885
	119
	120	Cmaf a verifier....
	121	preff = 101325.
	122	pa= 50000.
	123	C
	124	print*,'PARAMETRES A ENTRER DANS offline.def'
	125	print*,'Nombre effectif de traceurs'
	126	read(,) nq
	127	write(,) nq
	128	print*,'Splitting du pas de temps dynamique'
	129	read(,) nsplit_dyn
	130	write(,) nsplit_dyn
	131	print*,'Splitting du pas de temps physique'
	132	read(,) nsplit_phy
	133	write(,) nsplit_phy
	134	print*,'jour0, par rapport au debut du fichier fluxmass'
	135	read(,) jour0
	136	write(,) jour0
	137	print*,'redecoupage (T) ou standard (F)'
	138	read(*,'(l1)') redecoupage
	139	write(,) redecoupage
	140	print*,'Lecture fichier start (T)?'
	141	read(*,'(l1)') lectstart
	142	write(*,'(l1)') lectstart
	143	print*,'Traceurs 1 et 2 = Rn222 et Pb'
	144	read(*,'(l1)') rnpb
	145	write(*,'(l1)') rnpb
	146	print*,'Avant T ou arriere F?'
	147	read(*,'(l1)') avant
	148	write(*,'(l1)') avant
	149	print*,'Index min et max des traceurs a initaliser'
	150	read(,) iqmin,iqmax
	151	write(,) iqmin,iqmax
	152	print*,'localisation (lon, lat, intensite)'
	153
	154	if(iqmax.le.0) then
	155	iqmin=iqmax+1
	156	endif
	157
	158	do iq=iqmin,iqmax
	159	read(,) lon_stat(iq),lat_stat(iq),intensite(iq)
	160	write(,) lon_stat(iq),lat_stat(iq),intensite(iq)
	161	enddo
	162
	163	print*,'Records min et max pour le stokage et couche '
	164	s ,'verticale (specifique reseaux), def 0 0 1'
	165	read(,) recstkmin,recstkmax,ldec
	166	write(,) recstkmin,recstkmax,ldec
	167
	168	c redefinition du splitting
	169
	170	if(mod(nsplit_dyn,nsplit_phy).eq.0) then
	171	nsplit=nsplit_phy
	172	elseif (mod(nsplit_phy,nsplit_dyn).eq.0) then
	173	nsplit=nsplit_dyn
	174	else
	175	stop'prendre un des deux nsplit doit diviser l autre'
	176	endif
	177
	178	nsplit_dyn=nsplit_dyn/nsplit
	179	nsplit_phy=nsplit_phy/nsplit
	180
	181	print*,'nsplit,nsplit_dyn,nsplit_phy'
	182	s ,nsplit,nsplit_dyn,nsplit_phy
	183
	184	c Lecture eventuelle
	185	time_0=0.
	186	if(nq.gt.nqmx) stop'surdimensioner nqmx'
	187	print*,'av initial0'
	188	call initial0(nq*ijp1llm,q)
	189	print*,'av lectstart'
	190
	191	if(lectstart) then
	192	CALL dynetat0("start.nc",nq,vcov,ucov,
	193	. teta,q,masse,ps,phis, time_0)
	194	print*,'Lecture du start'
[204]	195	c On zappe le radon et le plomb.
	196	q(:,:,:,1)=q(:,:,:,3)
	197	q(:,:,:,2)=q(:,:,:,4)
	198	c Initialisation d'un traceur ` 1 pour tester l'impact du lessivage.
	199	q(:,:,:,3)=1.
[187]	200	else
	201	day_ini=0
	202	anne_ini=0
	203	endif
	204
	205	day_end=day_ini+nday
	206	C print*,'av defrun'
	207	c open (99,file='run.def',status='old',form='formatted')
	208	clefinit=.true.
	209	CALL defrun_new( 99, clefinit, clesphy0 )
	210	c close(99)
	211	C print*,'av iniconst'
	212	c lecture du jour de demarrage
[204]	213	c premiere cnitialisation, eventuellement bidon
	214	cALL iniconst
[187]	215	CALL inigeom
	216	C
[204]	217	print*,'ENTREE DANS redecoupenc ou lectfluxnc
[187]	218	s pour irec=0'
	219
	220	c Lecture des entetes des fichiers de flux
	221
	222	if(redecoupage) then
	223	call redecoupenc(0,masse,pbaru,pbarv,w,teta,phi,
	224	s nrec,avant,airefi,pphis,
	225	s t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
	226	s yu1,yv1,ftsol,pctsrf,
	227	s frac_impa,frac_nucl,phis)
	228	else
[204]	229	call lectfluxnc(0,masse,pbaru,pbarv,w,teta,phi,
[187]	230	s nrec,avant,airefi,pphis,
	231	s t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
	232	s yu1,yv1,ftsol,pctsrf,
	233	s frac_impa,frac_nucl,phis)
	234	call inigeom
	235	endif
	236
	237	c controle sur les records
	238
	239	if(avant) then
	240	recmin=2+(jour0-1)*iday_step
	241	recmax=1+(jour0+nday-1)*iday_step
	242	else
	243	recmin=(jour0-nday)*iday_step+2
	244	recmax=jour0*iday_step+1
	245	endif
	246
	247	if(recmin.le.1.or.recmax.ge.nrec) then
	248	print*,'Probleme de lecture '
	249	print*,'recmin,recmax,nrec',recmin,recmax,nrec
	250	stop
	251	endif
	252
	253	C on recalcule le nombre de pas de temps dans une journee
	254	c istdyn est la frequence de stokage en "pas de temps dynamique" dans
	255	c le modele on-line. dtvr est le pas de temps du meme modele en s.
	256	c iday-step est le nombre de lectures par jour (grand pas de temps).
	257
	258	iday_step=daysec/(istdyn*dtvr)
	259	C
	260	CALL iniconst
	261	c on calcule le pas de temps
	262	dtvr = daysec/FLOAT(iday_step)
	263	dtphys=dtvr/(nsplit*nsplit_phy) ! a mettre apres iniconst!
	264	print*,'Pas de temps physique ',dtphys
	265	C
	266	call gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
	267	print *,'Pas de temps dynamique', dtvr
	268	print*,'Pas de temps physique ',dtphys
	269
	270	c-----------------------------------------------------------------------
	271
	272	C maf est ce utile ?
	273	call suphec
	274	CALL inifilr
	275	print*,'Pas de temps physique ',dtphys
	276	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
	277	* tetagdiv, tetagrot , tetatemp )
	278	print*,'Pas de temps physique ',dtphys
	279
	280	c-----------------------------------------------------------------------
	281	c temps de depart et de fin:
	282	c --------------------------
	283	itau = 0
	284	iday = dayref
	285	time = time_0
	286	IF(time.GT.1.) THEN
	287	time = time-1.
	288	iday = iday+1
	289	ENDIF
	290	print*,'Pas de temps physique ',dtphys
	291	c=======================================================================
	292	c Initialisation des fichiers de sortie
	293	c=======================================================================
	294	print*,'Initialisation de wrgras'
	295	do i=1,iip1
	296	print*,'rlonu(',i,')=',rlonu(i)
	297	enddo
	298	do i=1,jjp1
	299	print*,'rlatu(',i,')=',rlatu(i)
	300	enddo
	301	do i=1,iip1
	302	print*,'rlonv(',i,')=',rlonv(i)
	303	enddo
	304	do i=1,jjm
	305	print*,'rlatv(',i,')=',rlatv(i)
	306	enddo
	307
	308	file='sol'
	309	call inigrads(1,iip1
	310	s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
	311	s ,llm,presnivs,1.
	312	s ,dtphys,file,'gcmq2 ')
	313
	314	file='atm'
	315	print*,'Pas de temps physique ',dtphys
	316	call inigrads(2,iip1
	317	s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
	318	s ,llm,presnivs,1.
	319	s ,dtphys,file,'gcmq2 ')
	320
	321	file='pbu'
	322	call inigrads(3,iip1
	323	s ,rlonu,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
	324	s ,llm,presnivs,1.
	325	s ,dtphys,file,'gcmq2 ')
	326
	327	file='pbv'
	328	call inigrads(4,iip1
	329	s ,rlonv,180./pi,-180.,180.,jjm,rlatv,-90.,90.,180./pi
	330	s ,llm,presnivs,1.
	331	s ,dtphys,file,'gcmq2 ')
	332
	333	file='qadv'
	334	call inigrads(5,iip1
	335	s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
	336	s ,llm,presnivs,1.
	337	s ,dtphys,file,'gcmq2 ')
	338
	339	print*,'dtphys ap inigrads ',dtphys
	340
	341	print*,'Ecriture du fichier de conditions aux limites'
	342	open (98,file='limit',form='unformatted')
	343	write(98) float(im),float(jm),float(nq),float(nday)
	344	write(98) (rlonv(i)*180./pi,i=1,iip1)
	345	write(98) (rlatu(j)*180./pi,j=1,jjp1)
	346	write(98) ((phis(i,j)/g,i=1,iip1),j=1,jjp1)
	347	close(98)
	348
	349	CALL dynredem0("restart.nc",day_end,anne_ini,phis,nq)
	350
	351	c-----------------------------------------------------------------------
	352	c Debut de la boucle en temps:
	353	c ----------------------------
	354
	355	do itau=1,nday*iday_step
	356	print*,'ITAU ITAU ITAU =',itau
	357	injecte='.'
	358	stoke='.'
	359	c Gestion du temps
	360	if (avant) then
	361	dtdyn=dtvr/float(nsplit*nsplit_dyn)
	362	c irec=1+(jour0-1)*iday_step+itau
	363	irec=(jour0-1)*iday_step+itau
	364	else
	365	dtdyn=-dtvr/float(nsplit*nsplit_dyn)
	366	irec=1+jour0*iday_step-itau+1
	367	endif
	368
	369	debut=itau.eq.1
	370
	371	c-----------------------------------------------------------------------
	372	c Lecture des fichiers fluxmass et physique:
	373	c -----------------------------------------------------
	374	print*,'CCCCC Appel a REDECOUPENC ou LECTFLUXNC
	375	s au pas itau=',itau
	376	if(redecoupage) then
	377	call redecoupenc(irec,masse,pbaru,pbarv,w,teta,phi,
	378	s nrec,avant,airefi,pphis,
	379	s t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
	380	s yu1,yv1,ftsol,pctsrf,
	381	s frac_impa,frac_nucl,phis)
	382	else
[204]	383	call lectfluxnc(irec,masse,pbaru,pbarv,w,teta,phi,
[187]	384	s nrec,avant,airefi,pphis,
	385	s t,zmfu, zmfd, zen_u, zde_u,zen_d, zde_d, coefkz,
	386	s yu1,yv1,ftsol,pctsrf,
	387	s frac_impa,frac_nucl,phis)
	388	endif
	389
[204]	390
	391	print,'TESTPHYS: ON PREND LA PUIS ',1./float(nsplitnsplit_phy)
	392	s ,' DES FRAC A IT=',itau
	393	frac_impa(:,:)=frac_impa(:,:)*(1./float(nsplitnsplit_phy))
	394	frac_nucl(:,:)=frac_nucl(:,:)*(1./float(nsplitnsplit_phy))
[187]	395	c ... ouverture du fichier de stockage netcdf ...
	396	C
	397	IF (debut) then
	398	dynhistave_file = 'histmoy.nc'
	399	day_ini=0
	400	anne_ini=0
	401	t_ops =(1./48.)*daysec
	402	t_wrt =(1./48.)*daysec
	403	dtav=dtvr/float(nsplit)
	404	mode=1
	405
[204]	406	c CALL initdynav(dynhistave_file,day_ini,anne_ini,dtav,
	407	c . t_ops, t_wrt, nq,mode, histaveid)
[187]	408
	409	pi=2.*asin(1.)
	410	do iq=iqmin,iqmax
	411	do l=1,llm
	412	do j=1,jjp1
	413	do i=1,iip1
	414	q(i,j,l,iq)=0.
	415	enddo
	416	enddo
	417	enddo
	418
	419	zlon=lon_stat(iq)*pi/180.
	420	zlat=lat_stat(iq)*pi/180.
	421
	422	CALL coordij(zlon,zlat,i_stat(iq),j_stat(iq))
	423
	424	enddo
	425	ENDIF
	426
	427	C calcul de la pression de surface.
	428
	429	do j=1,jjp1
	430	do i=1,iip1
	431	smass(i,j)=0.
	432	enddo
	433	enddo
	434
	435	do l=1,llm
	436	do j=1,jjp1
	437	do i=1,iip1
	438	smass(i,j)=smass(i,j)+masse(i,j,l)
	439	enddo
	440	enddo
	441	enddo
	442
	443	do j=1,jjp1
	444	do i=1,iip1
	445	ps (i,j)=smass(i,j)/aire(i,j)*g
	446	end do
	447	end do
	448	C
	449	C calcul de la pression et de pk (fonction d'Exner)
	450	C
	451	CALL pression ( ip1jmp1, ap, bp, ps, p )
	452	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	453	C
	454	do jj=1,jjp1
	455	do ii=1,iip1
	456	phis(ii,jj)= phi(ii,jj,1)-teta(ii,jj,1)*
	457	s (pks(ii,jj)-pk(ii,jj,1))
	458	end do
	459	end do
	460	c=======================================================================
	461	c TERMES SOURCES OU PUITS
	462	do isplit=1,nsplit
	463	c=======================================================================
	464	c injection par Fred
	465	c Print*,Verifier en fct de iapp_tracvl la boucle inj
	466	goto 333
	467	c Inj pour le stokage ttes les 0.5h
	468	if (itau.gt.32.and.itau.le.56)then
	469	print*,'Masse de la maille inj :',
	470	s masse(i_stat(1),j_stat(1),1)
	471	do iiinj=1,20
	472	q(i_stat(1),j_stat(1),1,1)=
	473	s q(i_stat(1),j_stat(1),1,1)
	474	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	475	s (float(nsplit)24.20.)
	476
	477	tinj1=tinj1+intensite(1)/(float(nsplit)24.20.)
	478	enddo
	479
	480	print*,'Apres injection a itau ',itau
	481
	482	do iq=1,nq
	483	write(str2,'(i2.2)') iq
	484	call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
	485	enddo
	486
	487	print*,'QT TOTALE INJECTEE TR1',tinj1
	488	print*,'fin des injections'
	489
	490	endif
	491
	492	c Inj pour le stokage ttes les 1h
	493
	494	if (itau.gt.16.and.itau.le.28)then
	495	print*,'Masse de la maille inj :',
	496	s masse(i_stat(1),j_stat(1),1)
	497	do iiinj=1,20
	498	q(i_stat(1),j_stat(1),1,1)=
	499	s q(i_stat(1),j_stat(1),1,1)
	500	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	501	s (float(nsplit)12.20.)
	502
	503	tinj1=tinj1+intensite(1)/(float(nsplit)12.20.)
	504
	505	enddo
	506
	507	print*,'Apres injection a itau ',itau
	508
	509	do iq=1,nq
	510	write(str2,'(i2.2)') iq
	511	call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
	512	enddo
	513
	514	print*,'QT TOTALE INJECTEE TR1',tinj1
	515	print*,'fin des injections'
	516
	517	endif
	518
	519	c Inj pour le stokage ttes les 1.5h
	520
	521	if (itau.eq.11.and.isplit.gt.2)then
	522	do iiinj=1,20
	523	q(i_stat(1),j_stat(1),1,1)=
	524	s q(i_stat(1),j_stat(1),1,1)
	525	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	526	s (24.*20.)
	527	enddo
	528	endif
	529
	530	if (itau.ge.12.and.itau.le.18)then
	531	do iiinj=1,20
	532	q(i_stat(1),j_stat(1),1,1)=
	533	s q(i_stat(1),j_stat(1),1,1)
	534	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	535	s (24.*20.)
	536	enddo
	537	endif
	538
	539	if (itau.eq.19.and.isplit.le.2)then
	540	do iiinj=1,20
	541	q(i_stat(1),j_stat(1),1,1)=
	542	s q(i_stat(1),j_stat(1),1,1)
	543	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	544	s (24.*20.)
	545	enddo
	546
	547	endif
	548
	549	c Inj pour le stokage ttes les 4h
	550
	551	if (itau.ge.5.and.itau.le.7)then
	552	do iiinj=1,20
	553	q(i_stat(1),j_stat(1),1,1)=
	554	s q(i_stat(1),j_stat(1),1,1)
	555	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	556	s (24.*20.)
	557	enddo
	558	endif
	559
	560	c Inj pour le stokage ttes les 3h
	561
	562	if (itau.eq.6.and.isplit.gt.2)then
	563	do iiinj=1,20
	564	q(i_stat(1),j_stat(1),1,1)=
	565	s q(i_stat(1),j_stat(1),1,1)
	566	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	567	s (24.*20.)
	568	enddo
	569	endif
	570
	571	if (itau.ge.7.and.itau.le.9)then
	572	do iiinj=1,20
	573	q(i_stat(1),j_stat(1),1,1)=
	574	s q(i_stat(1),j_stat(1),1,1)
	575	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	576	s (24.*20.)
	577	enddo
	578	endif
	579
	580	if (itau.eq.10.and.isplit.le.2)then
	581	do iiinj=1,20
	582	q(i_stat(1),j_stat(1),1,1)=
	583	s q(i_stat(1),j_stat(1),1,1)
	584	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	585	s (24.*20.)
	586	enddo
	587	endif
	588
	589	c Inj pour le stokage ttes les 6h
	590
	591	if (itau.eq.3.and.isplit.gt.8)then
	592	do iiinj=1,20
	593	q(i_stat(1),j_stat(1),1,1)=
	594	s q(i_stat(1),j_stat(1),1,1)
	595	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	596	s (24.*20.)
	597	enddo
	598	endif
	599
	600	if (itau.eq.4)then
	601	do iiinj=1,20
	602	q(i_stat(1),j_stat(1),1,1)=
	603	s q(i_stat(1),j_stat(1),1,1)
	604	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	605	s (24.*20.)
	606	enddo
	607	endif
	608
	609	if (itau.eq.5.and.isplit.le.8)then
	610	do iiinj=1,20
	611	q(i_stat(1),j_stat(1),1,1)=
	612	s q(i_stat(1),j_stat(1),1,1)
	613	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	614	s (24.*20.)
	615	enddo
	616	endif
	617
	618	c Inj sur les 4 1ere couches pour le stokage ttes les 3h
	619
	620	if (itau.eq.6.and.isplit.gt.2)then
	621	do l=1,4
	622	q(i_stat(1),j_stat(1),l,1)=
	623	s q(i_stat(1),j_stat(1),l,1)
	624	s +intensite(1)/masse(i_stat(1),j_stat(1),l)/
	625	s (24.*4.)
	626	enddo
	627	endif
	628
	629	if (itau.ge.7.and.itau.le.9)then
	630	do l=1,4
	631	q(i_stat(1),j_stat(1),l,1)=
	632	s q(i_stat(1),j_stat(1),l,1)
	633	s +intensite(1)/masse(i_stat(1),j_stat(1),l)/
	634	s (24.*4.)
	635	enddo
	636	endif
	637
	638	if (itau.eq.10.and.isplit.le.2)then
	639	do l=1,4
	640	q(i_stat(1),j_stat(1),l,1)=
	641	s q(i_stat(1),j_stat(1),l,1)
	642	s +intensite(1)/masse(i_stat(1),j_stat(1),l)/
	643	s (24.*4.)
	644	enddo
	645	endif
	646
	647	c Inj pour le stokage ttes les 12h
	648
	649	if (itau.eq.2.and.isplit.gt.8)then
	650	do iiinj=1,20
	651	q(i_stat(1),j_stat(1),1,1)=
	652	s q(i_stat(1),j_stat(1),1,1)
	653	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	654	s (24.*20.)
	655	enddo
	656	endif
	657
	658	if (itau.eq.3.and.isplit.le.8)then
	659	do iiinj=1,20
	660	q(i_stat(1),j_stat(1),1,1)=
	661	s q(i_stat(1),j_stat(1),1,1)
	662	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	663	s (24.*20.)
	664	enddo
	665	endif
	666
	667	c Inj pour le stokage ttes les 24h
	668
	669	if (itau.eq.1.and.isplit.gt.32)then
	670	q(i_stat(1),j_stat(1),1,1)=
	671	s q(i_stat(1),j_stat(1),1,1)
	672	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/24.
	673	endif
	674
	675	if (itau.eq.2.and.isplit.le.8)then
	676	q(i_stat(1),j_stat(1),1,1)=
	677	s q(i_stat(1),j_stat(1),1,1)
	678	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/24.
	679	endif
	680
	681	c Inj pour le stokage ttes les 3h
	682
	683	if (itau.eq.6.and.isplit.gt.2)then
	684	do iiinj=1,20
	685	q(i_stat(1),j_stat(1),1,1)=
	686	s q(i_stat(1),j_stat(1),1,1)
	687	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	688	s (24.*20.)
	689	enddo
	690	endif
	691
	692	if (itau.ge.7.and.itau.le.9)then
	693	do iiinj=1,20
	694	q(i_stat(1),j_stat(1),1,1)=
	695	s q(i_stat(1),j_stat(1),1,1)
	696	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	697	s (24.*20.)
	698	enddo
	699	endif
	700
	701	if (itau.eq.10.and.isplit.le.2)then
	702	do iiinj=1,20
	703	q(i_stat(1),j_stat(1),1,1)=
	704	s q(i_stat(1),j_stat(1),1,1)
	705	s +intensite(1)/masse(i_stat(1),j_stat(1),1)/
	706	s (24.*20.)
	707	enddo
	708	endif
	709
	710	333 continue
	711
	712	if (itau.eq.6.and.isplit.gt.2)then
	713
	714	do l=1,2
	715	q(i_stat(1),j_stat(1),l,1)=
	716	s q(i_stat(1),j_stat(1),l,1)
	717	s +intensite(1)/masse(i_stat(1),j_stat(1),l)/
	718	s (24.*2.)
	719	enddo
	720	endif
	721
	722	if (itau.ge.7.and.itau.le.9)then
	723	do l=1,2
	724	q(i_stat(1),j_stat(1),l,1)=
	725	s q(i_stat(1),j_stat(1),l,1)
	726	s +intensite(1)/masse(i_stat(1),j_stat(1),l)/
	727	s (24.*2.)
	728	enddo
	729	endif
	730
	731	if (itau.eq.10.and.isplit.le.2)then
	732	do l=1,2
	733	q(i_stat(1),j_stat(1),l,1)=
	734	s q(i_stat(1),j_stat(1),l,1)
	735	s +intensite(1)/masse(i_stat(1),j_stat(1),l)/
	736	s (24.*2.)
	737	enddo
	738	endif
	739
	740	c=======================================================================
	741	c FIN DE LA PARTIE TERMES SOURCES
	742	c=======================================================================
	743	c-----------------------------------------------------------------------
	744	c Advection
	745	c ----------
	746	print*,'ENTREE DANS VLSPLT'
	747	c=======================================================================
	748	c do isplit=1,nsplit
	749	c=======================================================================
	750	CC MAF dans le cas du offline on ne fait pas la difference masse et massem
	751	do iii=1,nsplit_dyn
	752	nom='pbaru'
	753	call wrgrads(3,llm,pbaru(:,:,1),nom,nom)
	754	nom='pbarv'
	755	call wrgrads(4,llm,pbarv(:,:,1),nom,nom)
	756	nom='masse'
	757	call wrgrads(3,llm,masse(:,:,1),nom,nom)
	758	nom='w'
	759	call wrgrads(3,llm,w(:,:,1),nom,nom)
	760	print*,'Entree dans vlsplt a itau et au temps',itau,time
	761	print*,'Pas advect :',dtdyn
	762	print*,'Avant vlsplt'
	763	do iq=1,nq
	764	write(str2,'(i2.2)') iq
	765	call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
	766	enddo
	767	do iq=1,nq
	768	ctest abder CALL vlsplt(q(1,1,1,iq),2.,masse,w,pbaru,pbarv,dtdyn)
	769	CALL vlsplt(q(:,:,:,iq),2.,masse,w,pbaru,pbarv,dtdyn)
	770	enddo
	771
	772	do iq=iqmin,iqmax
	773	write(str2,'(i2.2)') iq
	774	nom='q'//str2
	775	call wrgrads(5,1,q(:,:,1,iq),nom,nom)
	776	enddo
	777	print*,'Apres vlsplt'
	778	do iq=1,nq
	779	write(str2,'(i2.2)') iq
	780	call diagadv(q(:,:,:,iq),masse,'Traceur '//str2)
	781	enddo
	782
	783	c mise a jour du champ de masse tenant compte de l'advection sur
	784	c le pas de temps considere.
	785
	786	enddo
	787
	788	c-----------------------------------------------------------------------
	789	c Physique (lessivage dans phytrac maintenant) MAF
	790	c ---------------------------------------------------------
	791	if(physic) then
	792	if(debut) then
	793	latfi(1)=rlatu(1)*180./pi
	794	lonfi(1)=0.
	795	DO j=2,jjm
	796	DO i=1,iim
	797	latfi((j-2)iim+1+i)= rlatu(j)180./pi
	798	lonfi((j-2)iim+1+i)= rlonv(i)180./pi
	799	ENDDO
	800	ENDDO
	801	latfi(ngridmx)= rlatu(jjp1)*180./pi
	802	lonfi(ngridmx)= 0.
	803	endif
	804	c
	805	C Calcul de paprs et pplay et de temp
	806	c -----------------------------------------------------------------
	807	c .... paprs definis aux (llm +1) interfaces des couches ....
	808	c .... pplay definis aux ( llm ) milieux des couches ....
	809	c -----------------------------------------------------------------
	810	C
	811	DO l = 1, llmp1
	812	paprs( 1,l ) = p(1,1,l)
	813	ig0 = 2
	814	DO j = 2, jjm
	815	DO i =1, iim
	816	paprs( ig0,l ) = p(i,j,l)
	817	ig0 = ig0 +1
	818	ENDDO
	819	ENDDO
	820	paprs( ngridmx,l ) = p(1,jjp1,l)
	821	ENDDO
	822	c
	823	unskap = 1./ kappa
	824	DO l=1,llm
	825	pksurcp = pk(1,1,l) / cpp
	826	pplay(1,l) = preff * pksurcp ** unskap
	827	tempfi(1,l) = teta(1,1,l) * pksurcp
	828	ig0 = 2
	829	DO j = 2, jjm
	830	DO i = 1, iim
	831	pksurcp = pk(i,j,l) / cpp
	832	pplay(ig0,l) = preff * pksurcp ** unskap
	833	tempfi(ig0,l) = teta(i,j,l) * pksurcp
	834	ig0 = ig0 + 1
	835	ENDDO
	836	ENDDO
	837	pksurcp = pk(1,jjp1,l) / cpp
	838	pplay(ig0,l) = preff * pksurcp ** unskap
	839	tempfi (ig0,l) = teta(1,jjp1,l) * pksurcp
	840	ENDDO
	841
	842	c On passe le traceur et le geopot de surf sur la grille physique
	843	call gr_dyn_fi(llm*nq,iip1,jjp1,ngridmx,q,qfi)
	844	c call gr_dyn_fi(1,iip1,jjp1,ngridmx,phis,pphis)
	845	do iii=1,nsplit_phy
	846	C
[204]	847	lafin=.false. ! en attendant mieux.
[187]	848	print*,'dtphys avant phytrac ',dtphys
[204]	849	print*,'TESTPHYS: APPEL A PHYTRAC IT=',itau
	850
[220]	851	call phytrac(rnpb,
[204]	852	I debutphy, lafin,
[187]	853	I nq,
	854	I ngridmx,llm,dtphys,
	855	I t,paprs,pplay,
	856	I zmfu, zmfd, zen_u, zde_u, zen_d, zde_d,
	857	I coefkz,yu1,yv1,ftsol,pctsrf,latfi,
	858	I frac_impa, frac_nucl,
	859	I lonfi,presnivs,airefi,pphis,
	860	O qfi)
	861	C
	862	debutphy= .FALSE.
	863	C
	864	enddo ! nsplit_phy
	865	C
	866	c On passe le traceur sur la grille dynamique.
	867	call gr_fi_dyn(llm*nq,ngridmx,iip1,jjp1,qfi,q)
	868
	869	endif
	870
	871	itauav=(itau-1)*nsplit+isplit
	872	c itauav=itau*nsplit+isplit
[204]	873	c CALL writedynav(histaveid, nq,mode, itauav,vcov ,
	874	c , ucov,teta,pk,phi,q,masse,ps,phis)
[187]	875	c qmoy(:,:,:)=qmoy(:,:,:)+q(:,:,1,:)
	876	do iq=1,nq
	877	do j=1,jjp1
	878	do i=1,iip1
	879	qmoy(i,j,iq)=qmoy(i,j,iq)+q(i,j,1,iq)
	880	enddo
	881	enddo
	882	enddo
	883	enddo ! isplit=1,nsplit
	884
	885	call histsync
	886
	887	c-----------------------------------------------------------------------
	888	c Calcul de ucov et vcov pour les sorties
	889	c -------------------------------------
	890
	891	call massbar(masse, massebx, masseby )
	892	do l=1,llm
	893	do j=1,jjp1
	894	do i=1,iip1
	895	ucov(i,j,l)=pbaru(i,j,l)/massebx(i,j,l)cu(i,j)cu(i,j)
	896	s /istdyn
	897	enddo
	898	enddo
	899
	900	do j=1,jjm
	901	do i=1,iip1
	902	vcov(i,j,l)=pbarv(i,j,l)/masseby(i,j,l)cv(i,j)cv(i,j)
	903	s /istdyn
	904	enddo
	905	enddo
	906	enddo
	907	iday= dayref+itau/iday_step
	908	time=
	909	s float(itau-(iday-dayref)*iday_step)/iday_step+time_0
	910	IF(time.GT.1.) THEN
	911	time = time-1.
	912	iday = iday+1
	913	ENDIF
	914	C
	915	c=======================================================================
	916	do iq=iqmin,iqmax
	917	write(str2,'(i2.2)') iq
	918	nom='q'//str2
	919	call wrgrads(1,1,qmoy(:,:,iq),nom,nom)
	920	call wrgrads(2,1,q(:,:,ldec,iq),nom,nom)
	921	enddo
	922
	923	do iq=1,nq
	924	do j=1,jjp1
	925	do i=1,iip1
	926	qmoy(i,j,iq)=0.
	927	enddo
	928	enddo
	929	enddo
	930
	931	print*,'Record=',irec,' stoke=',stoke,
	932	s ' injecte=',injecte
	933
	934	enddo ! itau
	935
	936	CALL dynredem1("restart.nc",0.0,
	937	. vcov,ucov,teta,q,nq,masse,ps)
	938
	939	c=======================================================================
	940	C fermeture des fichiers netcdf
	941	CALL histclo
	942	C
	943	300 FORMAT('1',/,15x,'run du pas',i7,2x,'au pas',i7,2x,
	944	s 'c"est a dire du jour',i7,3x,'au jour',i7,/,/)
	945	C----------------------------------------------------------------------
	946
	947	stop
	948	end

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

Download in other formats:

Original Format