Context Navigation

offlinenc.F @ 195

Last change on this file since 195 was 187, checked in by lmdzadmin, 24 years ago
Ajout Idelkadi pour le offline LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.0 KB

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

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

Download in other formats:

Original Format