Context Navigation

offlinenc.F @ 2523

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

Line
1	c
2	c $Header$
3	c
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
71	logical debut,lectstart,rnpb,lafin
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'
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.
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
213	c premiere cnitialisation, eventuellement bidon
214	cALL iniconst
215	CALL inigeom
216	C
217	print*,'ENTREE DANS redecoupenc ou lectfluxnc
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
229	call lectfluxnc(0,masse,pbaru,pbarv,w,teta,phi,
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
383	call lectfluxnc(irec,masse,pbaru,pbarv,w,teta,phi,
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
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))
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
406	c CALL initdynav(dynhistave_file,day_ini,anne_ini,dtav,
407	c . t_ops, t_wrt, nq,mode, histaveid)
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
847	lafin=.false. ! en attendant mieux.
848	print*,'dtphys avant phytrac ',dtphys
849	print*,'TESTPHYS: APPEL A PHYTRAC IT=',itau
850
851	call phytrac(rnpb,
852	I debutphy, lafin,
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
873	c CALL writedynav(histaveid, nq,mode, itauav,vcov ,
874	c , ucov,teta,pk,phi,q,masse,ps,phis)
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