Context Navigation

gcm.F @ 477

Last change on this file since 477 was 477, checked in by lmdzadmin, 21 years ago
Quand on remet a zero le calendrier, il faut aussi reinitialiser day_ini qui sert a la determination du jour a lire dans le fichier de limite et au calcul de l'orbite LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 24.3 KB

Line
1
2	PROGRAM gcm
3
4	USE IOIPSL
5
6	IMPLICIT NONE
7
8	c ...... Version du 10/01/98 ..........
9
10	c avec coordonnees verticales hybrides
11	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
12
13	c=======================================================================
14	c
15	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
16	c -------
17	c
18	c Objet:
19	c ------
20	c
21	c GCM LMD nouvelle grille
22	c
23	c=======================================================================
24	c
25	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
26	c et possibilite d'appeler une fonction f(y) a derivee tangente
27	c hyperbolique a la place de la fonction a derivee sinusoidale.
28
29	c ... Possibilite de choisir le shema de Van-leer pour l'advection de
30	c q , en faisant iadv = 3 dans traceur (29/04/97) .
31	c
32	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
33	c
34	c-----------------------------------------------------------------------
35	c Declarations:
36	c -------------
37
38	#include "dimensions.h"
39	#include "paramet.h"
40	#include "comconst.h"
41	#include "comdissnew.h"
42	#include "comvert.h"
43	#include "comgeom.h"
44	#include "logic.h"
45	#include "temps.h"
46	#include "control.h"
47	#include "ener.h"
48	#include "netcdf.inc"
49	#include "description.h"
50	#include "serre.h"
51	#include "tracstoke.h"
52
53
54	INTEGER longcles
55	PARAMETER ( longcles = 20 )
56	REAL clesphy0( longcles )
57	SAVE clesphy0
58
59	INTEGER*4 iday ! jour julien
60	REAL time ! Heure de la journee en fraction d'1 jour
61	REAL zdtvr
62	INTEGER nbetatmoy, nbetatdem,nbetat
63
64	c variables dynamiques
65	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
66	REAL teta(ip1jmp1,llm) ! temperature potentielle
67	REAL q(ip1jmp1,llm,nqmx) ! champs advectes
68	REAL ps(ip1jmp1) ! pression au sol
69	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
70	REAL pks(ip1jmp1) ! exner au sol
71	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
72	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
73	REAL masse(ip1jmp1,llm) ! masse d'air
74	REAL phis(ip1jmp1) ! geopotentiel au sol
75	REAL phi(ip1jmp1,llm) ! geopotentiel
76	REAL w(ip1jmp1,llm) ! vitesse verticale
77
78	c variables dynamiques intermediaire pour le transport
79	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
80
81	c variables dynamiques au pas -1
82	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
83	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
84	REAL massem1(ip1jmp1,llm)
85
86	c tendances dynamiques
87	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
88	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
89
90	c tendances de la dissipation
91	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
92	REAL dhdis(ip1jmp1,llm)
93
94	c tendances physiques
95	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
96	REAL dhfi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
97
98	c variables pour le fichier histoire
99	REAL dtav ! intervalle de temps elementaire
100
101	REAL tppn(iim),tpps(iim),tpn,tps
102	c
103	INTEGER iadv(nqmx) ! indice schema de transport pour le traceur iq
104
105	INTEGER itau,itaufinp1,iav
106
107	EXTERNAL caldyn, traceur
108	EXTERNAL dissip,geopot,iniconst,inifilr
109	EXTERNAL integrd,SCOPY
110	EXTERNAL iniav,writeav,writehist
111	EXTERNAL inigeom
112	EXTERNAL exner_hyb,addit
113	EXTERNAL defrun_new, test_period
114	REAL SSUM
115	REAL time_0 , finvmaold(ip1jmp1,llm)
116
117	LOGICAL lafin
118	INTEGER ij,iq,l,numvanle,iapp_tracvl
119
120
121	INTEGER fluxid, fluxvid,fluxdid
122	integer histid, histvid, histaveid
123	real time_step, t_wrt, t_ops
124
125	REAL rdayvrai,rdaym_ini,rday_ecri
126	LOGICAL first
127	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
128	c+jld variables test conservation energie
129	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
130	C Tendance de la temp. potentiel d (theta)/ d t due a la
131	C tansformation d'energie cinetique en energie thermique
132	C cree par la dissipation
133	REAL dhecdt(ip1jmp1,llm)
134	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
135	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
136	CHARACTER*15 ztit
137	INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
138	SAVE ip_ebil_dyn
139	DATA ip_ebil_dyn/0/
140	c-jld
141
142	LOGICAL offline ! Controle du stockage ds "fluxmass"
143	PARAMETER (offline=.false.)
144
145	character*80 dynhist_file, dynhistave_file
146	character*20 modname
147	character*80 abort_message
148
149	C Calendrier
150	LOGICAL true_calendar
151	PARAMETER (true_calendar = .false.)
152	C Run nudge
153	LOGICAL ok_nudge
154	PARAMETER (ok_nudge = .false.)
155
156	c-----------------------------------------------------------------------
157	c Initialisations:
158	c ----------------
159
160	abort_message = 'last timestep reached'
161	modname = 'gcm'
162	descript = 'Run GCM LMDZ'
163	lafin = .FALSE.
164	dynhist_file = 'dyn_hist.nc'
165	dynhistave_file = 'dyn_hist_ave.nc'
166
167	if (true_calendar) then
168	call ioconf_calendar('gregorian')
169	else
170	call ioconf_calendar('360d')
171	endif
172	c-----------------------------------------------------------------------
173	c
174	c .... Choix des shemas d'advection pour l'eau et les traceurs ...
175	c ...................................................................
176	c
177	c iadv = 1 shema transport type "humidite specifique LMD"
178	c iadv = 2 shema amont
179	c iadv = 3 shema Van-leer
180	c iadv = 4 schema Van-leer + humidite specifique
181	c Modif F.Codron
182	c
183	c dans le tableau q(ij,l,iq) , iq = 1 pour l'eau vapeur
184	c , iq = 2 pour l'eau liquide
185	c et eventuellement , iq = 3, nqmx pour les autres traceurs
186	c
187	DO iq = 1, nqmx
188	iadv( iq ) = 3
189	ENDDO
190	c
191	c iadv(1): choix pour l'eau vap. et iadv(2) : choix pour l'eau liq.
192	iadv( 1 ) = 4
193	c
194	IF( iadv(1).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
195	* ,' pour la vapeur d''eau'
196	IF( iadv(1).EQ.2 ) PRINT *,' Choix du shema amont',' pour la'
197	* ,' vapeur d''eau '
198	IF( iadv(1).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour'
199	* ,'la vapeur d''eau'
200	IF( iadv(1).EQ.4 ) PRINT *,' Choix du shema Van-Leer + humidite'
201	* ,' specifique pour la vapeur d''eau'
202	c
203	IF( iadv(1).LE.0.OR.iadv(1).GT.4 ) THEN
204	PRINT *,' Erreur dans le choix de iadv (1).Corriger et repasser
205	* , car iadv(1) = ', iadv(1)
206	CALL ABORT
207	ENDIF
208
209	DO iq = 2, nqmx
210	IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
211	* ,' pour le traceur no ', iq
212	IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema amont',' pour le'
213	* ,' traceur no ', iq
214	IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour'
215	* ,'le traceur no ', iq
216
217	IF( iadv(iq).EQ.4 ) THEN
218	PRINT *,' Le shema Van-Leer + humidite specifique ',
219	* ' est uniquement pour la vapeur d eau .'
220	PRINT *,' Corriger iadv( ',iq, ') et repasser ! '
221	CALL ABORT
222	ENDIF
223
224	IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 ) THEN
225	PRINT *,' Erreur dans le choix de iadv . Corriger et repasser
226	* . '
227	CALL ABORT
228	ENDIF
229	ENDDO
230	c
231	first = .TRUE.
232	numvanle = nqmx + 1
233	DO iq = 1, nqmx
234	IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN
235	numvanle = iq
236	first = .FALSE.
237	ENDIF
238	ENDDO
239	c
240	DO iq = 1, nqmx
241	IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle ) THEN
242	PRINT *,' Il y a discontinuite dans le choix du shema de ',
243	* 'Van-leer pour les traceurs . Corriger et repasser . '
244	CALL ABORT
245	ENDIF
246	IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 ) THEN
247	PRINT *,' Le choix de iadv est errone pour le traceur ',
248	* iq
249	STOP
250	ENDIF
251	ENDDO
252	c
253
254	CALL dynetat0("start.nc",nqmx,vcov,ucov,
255	. teta,q,masse,ps,phis, time_0)
256
257
258	c OPEN( 99,file ='run.def',status='old',form='formatted')
259	c CALL defrun_new( 99, .TRUE. , clesphy0 )
260	CALL conf_gcm( 99, .TRUE. , clesphy0 )
261	c CLOSE(99)
262
263	c on recalcule eventuellement le pas de temps
264
265	IF(MOD(day_step,iperiod).NE.0)
266	* STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
267
268	IF(MOD(day_step,iphysiq).NE.0)
269	* STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
270
271	zdtvr = daysec/FLOAT(day_step)
272	IF(dtvr.NE.zdtvr) THEN
273	PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
274	ENDIF
275	C
276	C on remet le calendrier à zero
277	c
278	if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
279	write(,)' Attention les dates initiales lues dans le fichier'
280	write(,)' restart ne correspondent pas a celles lues dans '
281	write(,)' gcm.def'
282	if (raz_date .ne. 1) then
283	write(,)' On garde les dates du fichier restart'
284	else
285	annee_ref = anneeref
286	day_ref = dayref
287	day_ini = dayref
288	itau_dyn = 0
289	itau_phy = 0
290	time_0 = 0.
291	write(,)' On reinitialise a la date lue dans gcm.def'
292	endif
293	ELSE
294	raz_date = 0
295	endif
296
297	c nombre d'etats dans les fichiers demarrage et histoire
298	nbetatdem = nday / iecri
299	nbetatmoy = nday / periodav + 1
300
301	dtvr = zdtvr
302	CALL iniconst
303	CALL inigeom
304
305	CALL inifilr
306	c
307	c ...... P.Le Van ( modif le 29/04/97 ) .........
308	c
309	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
310	* tetagdiv, tetagrot , tetatemp )
311	c+jld
312	C initialisation constantes thermo utilisees dans diagedyn
313	IF (ip_ebil_dyn.ge.1 ) THEN
314	CALL suphec
315	ENDIF
316	c-jld
317	c
318
319	CALL pression ( ip1jmp1, ap, bp, ps, p )
320	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
321	c
322
323	c numero de stockage pour les fichiers de redemarrage:
324
325	c-----------------------------------------------------------------------
326	c temps de depart et de fin:
327	c --------------------------
328
329	itau = 0
330	iday = day_ini+itau/day_step
331	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
332	IF(time.GT.1.) THEN
333	time = time-1.
334	iday = iday+1
335	ENDIF
336	itaufin = nday*day_step
337	itaufinp1 = itaufin +1
338
339	day_end = day_ini + nday
340	PRINT 300, itau,itaufin,day_ini,day_end
341
342	CALL dynredem0("restart.nc", day_end, phis, nqmx)
343
344	ecripar = .TRUE.
345
346	time_step = zdtvr
347	t_ops = iecri * daysec
348	t_wrt = iecri * daysec
349	C CALL inithist(dynhist_file,day_ref,annee_ref,time_step,
350	c . t_ops, t_wrt, nqmx, histid, histvid)
351
352	t_ops = iperiod * time_step
353	t_wrt = periodav * daysec
354	CALL initdynav(dynhistave_file,day_ref,annee_ref,time_step,
355	. t_ops, t_wrt, nqmx, histaveid)
356
357	dtav = iperiod*dtvr/daysec
358
359
360	c Quelques initialisations pour les traceurs
361	call initial0(ijp1llm*nqmx,dq)
362	istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
363	istphy=istdyn/iphysiq
364
365	c-----------------------------------------------------------------------
366	c Debut de l'integration temporelle:
367	c ----------------------------------
368
369	1 CONTINUE
370
371
372
373	if (ok_nudge) then
374	call nudge(itau,ucov,vcov,teta,masse,ps)
375	endif
376	c
377	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
378	CALL test_period ( ucov,vcov,teta,q,p,phis )
379	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
380	c ENDIF
381	c
382	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
383	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
384	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
385	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
386	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
387
388	forward = .TRUE.
389	leapf = .FALSE.
390	dt = dtvr
391
392	c ... P.Le Van .26/04/94 ....
393
394	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
395	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
396
397
398	2 CONTINUE
399
400	c-----------------------------------------------------------------------
401
402	c date:
403	c -----
404
405
406	c gestion des appels de la physique et des dissipations:
407	c ------------------------------------------------------
408	c
409	c ... P.Le Van ( 6/02/95 ) ....
410
411	apphys = .FALSE.
412	statcl = .FALSE.
413	conser = .FALSE.
414	apdiss = .FALSE.
415
416	IF( purmats ) THEN
417	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
418	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
419	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
420	$ .AND. physic ) apphys = .TRUE.
421	ELSE
422	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
423	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
424	IF( MOD(itau+1,iphysiq).EQ.0. AND. physic ) apphys = .TRUE.
425	END IF
426
427	c-----------------------------------------------------------------------
428	c calcul des tendances dynamiques:
429	c --------------------------------
430
431	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
432	c
433	c
434	CALL caldyn
435	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
436	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
437
438	c-----------------------------------------------------------------------
439	c calcul des tendances advection des traceurs (dont l'humidite)
440	c -------------------------------------------------------------
441
442	IF( forward. OR . leapf ) THEN
443	c
444	DO 15 iq = 1, nqmx
445	c
446	IF( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 ) THEN
447	CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
448
449	ELSE IF( iq.EQ. nqmx ) THEN
450	c
451	iapp_tracvl = 5
452	c
453	cccc iapp_tracvl est la frequence en pas du groupement des flux
454	cccc de masse pour Van-Leer dans la routine tracvl .
455	c
456	CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
457	* p, masse, dq, iadv(1), teta, pk )
458	c
459	c ... Modif F.Codron ....
460	c
461	ENDIF
462	c
463	15 CONTINUE
464	C
465	IF (offline) THEN
466	Cmaf stokage du flux de masse pour traceurs OFF-LINE
467
468	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
469	. time_step, itau)
470
471
472	ENDIF
473	c
474	ENDIF
475
476
477	c-----------------------------------------------------------------------
478	c integrations dynamique et traceurs:
479	c ----------------------------------
480
481
482	CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
483	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
484	$ finvmaold )
485
486	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
487	c
488	c-----------------------------------------------------------------------
489	c calcul des tendances physiques:
490	c -------------------------------
491	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
492	c
493	IF( purmats ) THEN
494	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
495	ELSE
496	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
497	ENDIF
498	c
499	c
500	IF( apphys ) THEN
501	c
502	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
503	c
504
505	CALL pression ( ip1jmp1, ap, bp, ps, p )
506	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
507
508	rdaym_ini = itau * dtvr / daysec
509	rdayvrai = rdaym_ini + day_ini
510
511	IF ( ecritphy.LT.1. ) THEN
512	rday_ecri = rdaym_ini
513	ELSE
514	rday_ecri = NINT( rdayvrai )
515	ENDIF
516	c
517
518	c rajout debug
519	c lafin = .true.
520	c+jld
521	IF (ip_ebil_dyn.ge.1 ) THEN
522	ztit='bil dyn'
523	CALL diagedyn(ztit,2,1,1,dtphys
524	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
525	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
526	ENDIF
527	c-jld
528	CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time ,
529	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
530	$ du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi )
531
532	c ajout des tendances physiques:
533	c ------------------------------
534	CALL addfi( nqmx, dtphys, leapf, forward ,
535	$ ucov, vcov, teta , q ,ps ,
536	$ dufi, dvfi, dhfi , dqfi ,dpfi )
537	c
538	c+jld
539	IF (ip_ebil_dyn.ge.1 ) THEN
540	ztit='bil phys'
541	CALL diagedyn(ztit,2,1,1,dtphys
542	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
543	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
544	ENDIF
545	c-jld
546	ENDIF
547
548	CALL pression ( ip1jmp1, ap, bp, ps, p )
549	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
550
551	c-----------------------------------------------------------------------
552	c
553	c
554	c dissipation horizontale et verticale des petites echelles:
555	c ----------------------------------------------------------
556
557	IF(apdiss) THEN
558	c+jld
559	IF (ip_ebil_dyn.ge.2 ) THEN
560	ztit='avant dissip'
561	CALL diagedyn(ztit,2,2,2,dtvr
562	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
563	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
564	ENDIF
565	CALL covcont ( llm , ucov , vcov , ucont, vcont )
566	CALL enercin ( vcov , ucov , vcont , ucont , ecin0 )
567	c-jld
568
569	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
570
571	CALL addit( ijp1llm,ucov ,dudis,ucov )
572	CALL addit( ijmllm ,vcov ,dvdis,vcov )
573	c+jld
574	CALL covcont ( llm , ucov , vcov , ucont, vcont )
575	CALL enercin ( vcov , ucov , vcont , ucont , ecin )
576	C
577	C On rajoute la tendance due a la transform. Ec -> E therm. cree
578	C lors de la dissipation
579	DO l = 1, llm
580	DO ij = 1, ip1jmp1
581	dhecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
582	dhdis(ij,l) = dhdis(ij,l) + dhecdt(ij,l)
583	ENDDO
584	ENDDO
585	c-jld
586	CALL addit( ijp1llm,teta ,dhdis,teta )
587	c+jld
588	c
589	IF (ip_ebil_dyn.ge.2 ) THEN
590	ztit='apres dissip'
591	CALL diagedyn(ztit,2,2,2,dtdiss
592	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
593	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
594	ENDIF
595	c-jld
596
597	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
598	c ... Calcul de la valeur moyenne, unique de h aux poles .....
599	c
600
601	DO l = 1, llm
602	DO ij = 1,iim
603	tppn(ij) = aire( ij ) * teta( ij ,l)
604	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
605	ENDDO
606	tpn = SSUM(iim,tppn,1)/apoln
607	tps = SSUM(iim,tpps,1)/apols
608
609	DO ij = 1, iip1
610	teta( ij ,l) = tpn
611	teta(ij+ip1jm,l) = tps
612	ENDDO
613	ENDDO
614
615	DO ij = 1,iim
616	tppn(ij) = aire( ij ) * ps ( ij )
617	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
618	ENDDO
619	tpn = SSUM(iim,tppn,1)/apoln
620	tps = SSUM(iim,tpps,1)/apols
621
622	DO ij = 1, iip1
623	ps( ij ) = tpn
624	ps(ij+ip1jm) = tps
625	ENDDO
626
627
628	END IF
629
630	c ajout debug
631	c IF( lafin ) then
632	c abort_message = 'Simulation finished'
633	c call abort_gcm(modname,abort_message,0)
634	c ENDIF
635
636	c ********************************************************************
637	c ********************************************************************
638	c .... fin de l'integration dynamique et physique pour le pas itau ..
639	c ********************************************************************
640	c ********************************************************************
641
642	c preparation du pas d'integration suivant ......
643
644	IF ( .NOT.purmats ) THEN
645	c ........................................................
646	c .............. schema matsuno + leapfrog ..............
647	c ........................................................
648
649	IF(forward. OR. leapf) THEN
650	itau= itau + 1
651	iday= day_ini+itau/day_step
652	time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
653	IF(time.GT.1.) THEN
654	time = time-1.
655	iday = iday+1
656	ENDIF
657	ENDIF
658
659
660	IF( itau. EQ. itaufinp1 ) then
661
662	abort_message = 'Simulation finished'
663	call abort_gcm(modname,abort_message,0)
664	ENDIF
665	c-----------------------------------------------------------------------
666	c ecriture du fichier histoire moyenne:
667	c -------------------------------------
668
669	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
670	IF(itau.EQ.itaufin) THEN
671	iav=1
672	ELSE
673	iav=0
674	ENDIF
675	CALL writedynav(histaveid, nqmx, itau,vcov ,
676	, ucov,teta,pk,phi,q,masse,ps,phis)
677	cccIM cf. FH
678	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
679	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
680
681	ENDIF
682
683	c-----------------------------------------------------------------------
684	c ecriture de la bande histoire:
685	c ------------------------------
686
687	IF( MOD(itau,iecri*day_step).EQ.0) THEN
688
689	nbetat = nbetatdem
690	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
691	c CALL writehist( histid, histvid, nqmx, itau,vcov ,
692	c , ucov,teta,phi,q,masse,ps,phis)
693
694
695	ENDIF
696
697	IF(itau.EQ.itaufin) THEN
698
699
700	PRINT *,' Appel test_period avant redem ', itau
701	CALL test_period ( ucov,vcov,teta,q,p,phis )
702	CALL dynredem1("restart.nc",0.0,
703	. vcov,ucov,teta,q,nqmx,masse,ps)
704
705	CLOSE(99)
706	ENDIF
707
708	c-----------------------------------------------------------------------
709	c gestion de l'integration temporelle:
710	c ------------------------------------
711
712	IF( MOD(itau,iperiod).EQ.0 ) THEN
713	GO TO 1
714	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
715
716	IF( forward ) THEN
717	c fin du pas forward et debut du pas backward
718
719	forward = .FALSE.
720	leapf = .FALSE.
721	GO TO 2
722
723	ELSE
724	c fin du pas backward et debut du premier pas leapfrog
725
726	leapf = .TRUE.
727	dt = 2.*dtvr
728	GO TO 2
729	END IF
730	ELSE
731
732	c ...... pas leapfrog .....
733
734	leapf = .TRUE.
735	dt = 2.*dtvr
736	GO TO 2
737	END IF
738
739	ELSE
740
741	c ........................................................
742	c .............. schema matsuno ...............
743	c ........................................................
744	IF( forward ) THEN
745
746	itau = itau + 1
747	iday = day_ini+itau/day_step
748	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
749
750	IF(time.GT.1.) THEN
751	time = time-1.
752	iday = iday+1
753	ENDIF
754
755	forward = .FALSE.
756	IF( itau. EQ. itaufinp1 ) then
757	abort_message = 'Simulation finished'
758	call abort_gcm(modname,abort_message,0)
759	ENDIF
760	GO TO 2
761
762	ELSE
763
764	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
765	IF(itau.EQ.itaufin) THEN
766	iav=1
767	ELSE
768	iav=0
769	ENDIF
770	CALL writedynav(histaveid, nqmx, itau,vcov ,
771	, ucov,teta,pk,phi,q,masse,ps,phis)
772	cccIM cf. FH
773	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
774	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
775
776	ENDIF
777
778	IF(MOD(itau,iecri*day_step).EQ.0) THEN
779	nbetat = nbetatdem
780	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
781	c CALL writehist( histid, histvid, nqmx, itau,vcov ,
782	c , ucov,teta,phi,q,masse,ps,phis)
783	ENDIF
784
785	IF(itau.EQ.itaufin)
786	. CALL dynredem1("restart.nc",0.0,
787	. vcov,ucov,teta,q,nqmx,masse,ps)
788
789	forward = .TRUE.
790	GO TO 1
791
792	ENDIF
793
794	END IF
795
796	300 FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x,
797	. 'c''est a dire du jour',i7,3x,'au jour',i7//)
798
799	STOP
800	END

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

Download in other formats:

Original Format