Context Navigation

gcm.F @ 403

Last change on this file since 403 was 403, checked in by lmdzadmin, 22 years ago
Conservation de l'energie dans la dynamique (JLD) + diagnostics transports (FH) IM/LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 23.6 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 ! PRINT level for energy conserv. diag.
138	SAVE ip_ebil
139	DATA ip_ebil/1/
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	c iadv(1): choix pour l'eau vap. et iadv(2) : choix pour l'eau liq.
188	c
189	DO iq = 1, nqmx
190	iadv( iq ) = 3
191	ENDDO
192	c
193	IF( iadv(1).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
194	* ,' pour la vapeur d''eau'
195	IF( iadv(1).EQ.2 ) PRINT *,' Choix du shema amont',' pour la'
196	* ,' vapeur d''eau '
197	IF( iadv(1).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour'
198	* ,'la vapeur d''eau'
199	IF( iadv(1).EQ.4 ) PRINT *,' Choix du shema Van-Leer + humidite'
200	* ,' specifique pour la vapeur d''eau'
201	c
202	IF( iadv(1).LE.0.OR.iadv(1).GT.4 ) THEN
203	PRINT *,' Erreur dans le choix de iadv (1).Corriger et repasser
204	* , car iadv(1) = ', iadv(1)
205	CALL ABORT
206	ENDIF
207
208	DO iq = 2, nqmx
209	IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
210	* ,' pour le traceur no ', iq
211	IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema amont',' pour le'
212	* ,' traceur no ', iq
213	IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour'
214	* ,'le traceur no ', iq
215
216	IF( iadv(iq).EQ.4 ) THEN
217	PRINT *,' Le shema Van-Leer + humidite specifique ',
218	* ' est uniquement pour la vapeur d eau .'
219	PRINT *,' Corriger iadv( ',iq, ') et repasser ! '
220	CALL ABORT
221	ENDIF
222
223	IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 ) THEN
224	PRINT *,' Erreur dans le choix de iadv . Corriger et repasser
225	* . '
226	CALL ABORT
227	ENDIF
228	ENDDO
229	c
230	first = .TRUE.
231	numvanle = nqmx + 1
232	DO iq = 1, nqmx
233	IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN
234	numvanle = iq
235	first = .FALSE.
236	ENDIF
237	ENDDO
238	c
239	DO iq = 1, nqmx
240	IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle ) THEN
241	PRINT *,' Il y a discontinuite dans le choix du shema de ',
242	* 'Van-leer pour les traceurs . Corriger et repasser . '
243	CALL ABORT
244	ENDIF
245	IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 ) THEN
246	PRINT *,' Le choix de iadv est errone pour le traceur ',
247	* iq
248	STOP
249	ENDIF
250	ENDDO
251	c
252
253	CALL dynetat0("start.nc",nqmx,vcov,ucov,
254	. teta,q,masse,ps,phis, time_0)
255
256
257	c OPEN( 99,file ='run.def',status='old',form='formatted')
258	c CALL defrun_new( 99, .TRUE. , clesphy0 )
259	CALL conf_gcm( 99, .TRUE. , clesphy0 )
260	c CLOSE(99)
261
262	c on recalcule eventuellement le pas de temps
263
264	IF(MOD(day_step,iperiod).NE.0)
265	* STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
266
267	IF(MOD(day_step,iphysiq).NE.0)
268	* STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
269
270	zdtvr = daysec/FLOAT(day_step)
271	IF(dtvr.NE.zdtvr) THEN
272	PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
273	ENDIF
274
275	c nombre d'etats dans les fichiers demarrage et histoire
276	nbetatdem = nday / iecri
277	nbetatmoy = nday / periodav + 1
278
279	dtvr = zdtvr
280	CALL iniconst
281	CALL inigeom
282
283	CALL inifilr
284	c
285	c ...... P.Le Van ( modif le 29/04/97 ) .........
286	c
287	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
288	* tetagdiv, tetagrot , tetatemp )
289	c+jld
290	C initialisation constantes thermo utilisees dans diagedyn
291	IF (ip_ebil.ge.1 ) THEN
292	CALL suphec
293	ENDIF
294	c-jld
295	c
296
297	CALL pression ( ip1jmp1, ap, bp, ps, p )
298	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
299	c
300
301	c numero de stockage pour les fichiers de redemarrage:
302
303	c-----------------------------------------------------------------------
304	c temps de depart et de fin:
305	c --------------------------
306
307	itau = 0
308	iday = day_ini+itau/day_step
309	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
310	IF(time.GT.1.) THEN
311	time = time-1.
312	iday = iday+1
313	ENDIF
314	itaufin = nday*day_step
315	itaufinp1 = itaufin +1
316
317	day_end = day_ini + nday
318	PRINT 300, itau,itaufin,day_ini,day_end
319
320	CALL dynredem0("restart.nc", day_end, phis, nqmx)
321
322	ecripar = .TRUE.
323
324	time_step = zdtvr
325	t_ops = iecri * daysec
326	t_wrt = iecri * daysec
327	C CALL inithist(dynhist_file,day_ref,annee_ref,time_step,
328	c . t_ops, t_wrt, nqmx, histid, histvid)
329
330	t_ops = iperiod * time_step
331	t_wrt = periodav * daysec
332	CALL initdynav(dynhistave_file,day_ref,annee_ref,time_step,
333	. t_ops, t_wrt, nqmx, histaveid)
334
335	dtav = iperiod*dtvr/daysec
336
337
338	c Quelques initialisations pour les traceurs
339	call initial0(ijp1llm*nqmx,dq)
340	istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
341	istphy=istdyn/iphysiq
342
343	c-----------------------------------------------------------------------
344	c Debut de l'integration temporelle:
345	c ----------------------------------
346
347	1 CONTINUE
348
349
350
351	if (ok_nudge) then
352	call nudge(itau,ucov,vcov,teta,masse,ps)
353	endif
354	c
355	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
356	CALL test_period ( ucov,vcov,teta,q,p,phis )
357	PRINT *,' ---- Test_period apres continue OK ! -----', itau
358	c ENDIF
359	c
360	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
361	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
362	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
363	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
364	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
365
366	forward = .TRUE.
367	leapf = .FALSE.
368	dt = dtvr
369
370	c ... P.Le Van .26/04/94 ....
371
372	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
373	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
374
375
376	2 CONTINUE
377
378	c-----------------------------------------------------------------------
379
380	c date:
381	c -----
382
383
384	c gestion des appels de la physique et des dissipations:
385	c ------------------------------------------------------
386	c
387	c ... P.Le Van ( 6/02/95 ) ....
388
389	apphys = .FALSE.
390	statcl = .FALSE.
391	conser = .FALSE.
392	apdiss = .FALSE.
393
394	IF( purmats ) THEN
395	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
396	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
397	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
398	$ .AND. physic ) apphys = .TRUE.
399	ELSE
400	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
401	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
402	IF( MOD(itau+1,iphysiq).EQ.0. AND. physic ) apphys = .TRUE.
403	END IF
404
405	c-----------------------------------------------------------------------
406	c calcul des tendances dynamiques:
407	c --------------------------------
408
409	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
410	c
411	c
412	CALL caldyn
413	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
414	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
415
416	c-----------------------------------------------------------------------
417	c calcul des tendances advection des traceurs (dont l'humidite)
418	c -------------------------------------------------------------
419
420	IF( forward. OR . leapf ) THEN
421	c
422	DO 15 iq = 1, nqmx
423	c
424	IF( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 ) THEN
425	CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
426
427	ELSE IF( iq.EQ. nqmx ) THEN
428	c
429	iapp_tracvl = 5
430	c
431	cccc iapp_tracvl est la frequence en pas du groupement des flux
432	cccc de masse pour Van-Leer dans la routine tracvl .
433	c
434	CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
435	* p, masse, dq, iadv(1), teta, pk )
436	c
437	c ... Modif F.Codron ....
438	c
439	ENDIF
440	c
441	15 CONTINUE
442	C
443	IF (offline) THEN
444	Cmaf stokage du flux de masse pour traceurs OFF-LINE
445
446	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
447	. time_step, itau)
448
449
450	ENDIF
451	c
452	ENDIF
453
454
455	c-----------------------------------------------------------------------
456	c integrations dynamique et traceurs:
457	c ----------------------------------
458
459
460	CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
461	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
462	$ finvmaold )
463
464	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
465	c
466	c-----------------------------------------------------------------------
467	c calcul des tendances physiques:
468	c -------------------------------
469	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
470	c
471	IF( purmats ) THEN
472	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
473	ELSE
474	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
475	ENDIF
476	c
477	c
478	IF( apphys ) THEN
479	c
480	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
481	c
482
483	CALL pression ( ip1jmp1, ap, bp, ps, p )
484	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
485
486	rdaym_ini = itau * dtvr / daysec
487	rdayvrai = rdaym_ini + day_ini
488
489	IF ( ecritphy.LT.1. ) THEN
490	rday_ecri = rdaym_ini
491	ELSE
492	rday_ecri = NINT( rdayvrai )
493	ENDIF
494	c
495
496	c rajout debug
497	c lafin = .true.
498	c+jld
499	IF (ip_ebil.ge.1 ) THEN
500	ztit='bil dyn'
501	CALL diagedyn(ztit,2,1,1,dtphys
502	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
503	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
504	ENDIF
505	c-jld
506	CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time ,
507	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
508	$ du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi )
509
510	c ajout des tendances physiques:
511	c ------------------------------
512	CALL addfi( nqmx, dtphys, leapf, forward ,
513	$ ucov, vcov, teta , q ,ps ,
514	$ dufi, dvfi, dhfi , dqfi ,dpfi )
515	c
516	c+jld
517	IF (ip_ebil.ge.1 ) THEN
518	ztit='bil phys'
519	CALL diagedyn(ztit,2,1,1,dtphys
520	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
521	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
522	ENDIF
523	c-jld
524	ENDIF
525
526	CALL pression ( ip1jmp1, ap, bp, ps, p )
527	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
528
529	c-----------------------------------------------------------------------
530	c
531	c
532	c dissipation horizontale et verticale des petites echelles:
533	c ----------------------------------------------------------
534
535	IF(apdiss) THEN
536	c+jld
537	IF (ip_ebil.ge.2 ) THEN
538	ztit='avant dissip'
539	CALL diagedyn(ztit,2,2,2,dtvr
540	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
541	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
542	ENDIF
543	CALL covcont ( llm , ucov , vcov , ucont, vcont )
544	CALL enercin ( vcov , ucov , vcont , ucont , ecin0 )
545	c-jld
546
547	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
548
549	CALL addit( ijp1llm,ucov ,dudis,ucov )
550	CALL addit( ijmllm ,vcov ,dvdis,vcov )
551	c+jld
552	CALL covcont ( llm , ucov , vcov , ucont, vcont )
553	CALL enercin ( vcov , ucov , vcont , ucont , ecin )
554	C
555	C On rajoute la tendance due a la transform. Ec -> E therm. cree
556	C lors de la dissipation
557	DO l = 1, llm
558	DO ij = 1, ip1jmp1
559	dhecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
560	dhdis(ij,l) = dhdis(ij,l) + dhecdt(ij,l)
561	ENDDO
562	ENDDO
563	c-jld
564	CALL addit( ijp1llm,teta ,dhdis,teta )
565	c+jld
566	c
567	IF (ip_ebil.ge.2 ) THEN
568	ztit='apres dissip'
569	CALL diagedyn(ztit,2,2,2,dtdiss
570	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
571	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
572	ENDIF
573	c-jld
574
575	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
576	c ... Calcul de la valeur moyenne, unique de h aux poles .....
577	c
578
579	DO l = 1, llm
580	DO ij = 1,iim
581	tppn(ij) = aire( ij ) * teta( ij ,l)
582	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
583	ENDDO
584	tpn = SSUM(iim,tppn,1)/apoln
585	tps = SSUM(iim,tpps,1)/apols
586
587	DO ij = 1, iip1
588	teta( ij ,l) = tpn
589	teta(ij+ip1jm,l) = tps
590	ENDDO
591	ENDDO
592
593	DO ij = 1,iim
594	tppn(ij) = aire( ij ) * ps ( ij )
595	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
596	ENDDO
597	tpn = SSUM(iim,tppn,1)/apoln
598	tps = SSUM(iim,tpps,1)/apols
599
600	DO ij = 1, iip1
601	ps( ij ) = tpn
602	ps(ij+ip1jm) = tps
603	ENDDO
604
605
606	END IF
607
608	c ajout debug
609	c IF( lafin ) then
610	c abort_message = 'Simulation finished'
611	c call abort_gcm(modname,abort_message,0)
612	c ENDIF
613
614	c ********************************************************************
615	c ********************************************************************
616	c .... fin de l'integration dynamique et physique pour le pas itau ..
617	c ********************************************************************
618	c ********************************************************************
619
620	c preparation du pas d'integration suivant ......
621
622	IF ( .NOT.purmats ) THEN
623	c ........................................................
624	c .............. schema matsuno + leapfrog ..............
625	c ........................................................
626
627	IF(forward. OR. leapf) THEN
628	itau= itau + 1
629	iday= day_ini+itau/day_step
630	time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
631	IF(time.GT.1.) THEN
632	time = time-1.
633	iday = iday+1
634	ENDIF
635	ENDIF
636
637
638	IF( itau. EQ. itaufinp1 ) then
639
640	abort_message = 'Simulation finished'
641	call abort_gcm(modname,abort_message,0)
642	ENDIF
643	c-----------------------------------------------------------------------
644	c ecriture du fichier histoire moyenne:
645	c -------------------------------------
646
647	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
648	IF(itau.EQ.itaufin) THEN
649	iav=1
650	ELSE
651	iav=0
652	ENDIF
653	CALL writedynav(histaveid, nqmx, itau,vcov ,
654	, ucov,teta,pk,phi,q,masse,ps,phis)
655	cccIM cf. FH
656	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
657	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
658
659	ENDIF
660
661	c-----------------------------------------------------------------------
662	c ecriture de la bande histoire:
663	c ------------------------------
664
665	IF( MOD(itau,iecri*day_step).EQ.0) THEN
666
667	nbetat = nbetatdem
668	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
669	c CALL writehist( histid, histvid, nqmx, itau,vcov ,
670	c , ucov,teta,phi,q,masse,ps,phis)
671
672
673	ENDIF
674
675	IF(itau.EQ.itaufin) THEN
676
677
678	PRINT *,' Appel test_period avant redem ', itau
679	CALL test_period ( ucov,vcov,teta,q,p,phis )
680	CALL dynredem1("restart.nc",0.0,
681	. vcov,ucov,teta,q,nqmx,masse,ps)
682
683	CLOSE(99)
684	ENDIF
685
686	c-----------------------------------------------------------------------
687	c gestion de l'integration temporelle:
688	c ------------------------------------
689
690	IF( MOD(itau,iperiod).EQ.0 ) THEN
691	GO TO 1
692	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
693
694	IF( forward ) THEN
695	c fin du pas forward et debut du pas backward
696
697	forward = .FALSE.
698	leapf = .FALSE.
699	GO TO 2
700
701	ELSE
702	c fin du pas backward et debut du premier pas leapfrog
703
704	leapf = .TRUE.
705	dt = 2.*dtvr
706	GO TO 2
707	END IF
708	ELSE
709
710	c ...... pas leapfrog .....
711
712	leapf = .TRUE.
713	dt = 2.*dtvr
714	GO TO 2
715	END IF
716
717	ELSE
718
719	c ........................................................
720	c .............. schema matsuno ...............
721	c ........................................................
722	IF( forward ) THEN
723
724	itau = itau + 1
725	iday = day_ini+itau/day_step
726	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
727
728	IF(time.GT.1.) THEN
729	time = time-1.
730	iday = iday+1
731	ENDIF
732
733	forward = .FALSE.
734	IF( itau. EQ. itaufinp1 ) then
735	abort_message = 'Simulation finished'
736	call abort_gcm(modname,abort_message,0)
737	ENDIF
738	GO TO 2
739
740	ELSE
741
742	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
743	IF(itau.EQ.itaufin) THEN
744	iav=1
745	ELSE
746	iav=0
747	ENDIF
748	CALL writedynav(histaveid, nqmx, itau,vcov ,
749	, ucov,teta,pk,phi,q,masse,ps,phis)
750	cccIM cf. FH
751	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
752	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
753
754	ENDIF
755
756	IF(MOD(itau,iecri*day_step).EQ.0) THEN
757	nbetat = nbetatdem
758	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
759	c CALL writehist( histid, histvid, nqmx, itau,vcov ,
760	c , ucov,teta,phi,q,masse,ps,phis)
761	ENDIF
762
763	IF(itau.EQ.itaufin)
764	. CALL dynredem1("restart.nc",0.0,
765	. vcov,ucov,teta,q,nqmx,masse,ps)
766
767	forward = .TRUE.
768	GO TO 1
769
770	ENDIF
771
772	END IF
773
774	300 FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x,
775	. 'c''est a dire du jour',i7,3x,'au jour',i7//)
776
777	STOP
778	END

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

Download in other formats:

Original Format