Context Navigation

leapfrog.F @ 1201

Last change on this file since 1201 was 1201, checked in by Laurent Fairhead, 15 years ago
Modifications nécessaires a l'inclusion d'un calendrier réaliste. La date courante est calculée dans leapfrog.F et exprimée en Jour Julien (modifié). On en a profité pour faire un peu de ménage dans la gestion des dates du modèle. Dans la physique, on utilise les routines de passages entre calendrier Julien et Gregorien incluses dans IOIPSL pour calculer le nombre de jours écoulés depuis le 1er janvier (pour les conditions aux limites) ou l'equinoxe (pour le calcul de la longitude solaire). Le calcul de l'orbite reprend celui du gcm planétaire (codé par FH) On décide du calendrier à utiliser à l'aide du paramètre calend du run.def. Par défaut celui-ci est à earth_360d LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.0 KB

Line
1	!
2	! $Id: leapfrog.F 1201 2009-07-07 14:01:00Z fairhead $
3	!
4	c
5	c
6	SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
7	& time_0)
8
9
10	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
11	#ifdef CPP_IOIPSL
12	use IOIPSL
13	#endif
14	USE infotrac
15	USE guide_mod, ONLY : guide_main
16	USE write_field
17	IMPLICIT NONE
18
19	c ...... Version du 10/01/98 ..........
20
21	c avec coordonnees verticales hybrides
22	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
23
24	c=======================================================================
25	c
26	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
27	c -------
28	c
29	c Objet:
30	c ------
31	c
32	c GCM LMD nouvelle grille
33	c
34	c=======================================================================
35	c
36	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
37	c et possibilite d'appeler une fonction f(y) a derivee tangente
38	c hyperbolique a la place de la fonction a derivee sinusoidale.
39
40	c ... Possibilite de choisir le shema pour l'advection de
41	c q , en modifiant iadv dans traceur.def (10/02) .
42	c
43	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
44	c Pour Van-Leer iadv=10
45	c
46	c-----------------------------------------------------------------------
47	c Declarations:
48	c -------------
49
50	#include "dimensions.h"
51	#include "paramet.h"
52	#include "comconst.h"
53	#include "comdissnew.h"
54	#include "comvert.h"
55	#include "comgeom.h"
56	#include "logic.h"
57	#include "temps.h"
58	#include "control.h"
59	#include "ener.h"
60	#include "description.h"
61	#include "serre.h"
62	#include "com_io_dyn.h"
63	#include "iniprint.h"
64	#include "academic.h"
65
66	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
67	! #include "clesphys.h"
68
69	INTEGER longcles
70	PARAMETER ( longcles = 20 )
71	REAL clesphy0( longcles )
72
73	real zqmin,zqmax
74	INTEGER nbetatmoy, nbetatdem,nbetat
75
76	c variables dynamiques
77	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
78	REAL teta(ip1jmp1,llm) ! temperature potentielle
79	REAL q(ip1jmp1,llm,nqtot) ! champs advectes
80	REAL ps(ip1jmp1) ! pression au sol
81	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
82	REAL pks(ip1jmp1) ! exner au sol
83	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
84	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
85	REAL masse(ip1jmp1,llm) ! masse d'air
86	REAL phis(ip1jmp1) ! geopotentiel au sol
87	REAL phi(ip1jmp1,llm) ! geopotentiel
88	REAL w(ip1jmp1,llm) ! vitesse verticale
89
90	c variables dynamiques intermediaire pour le transport
91	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
92
93	c variables dynamiques au pas -1
94	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
95	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
96	REAL massem1(ip1jmp1,llm)
97
98	c tendances dynamiques
99	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
100	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
101
102	c tendances de la dissipation
103	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
104	REAL dtetadis(ip1jmp1,llm)
105
106	c tendances physiques
107	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
108	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
109
110	c variables pour le fichier histoire
111	REAL dtav ! intervalle de temps elementaire
112
113	REAL tppn(iim),tpps(iim),tpn,tps
114	c
115	INTEGER itau,itaufinp1,iav
116	! INTEGER iday ! jour julien
117	REAL time
118
119	REAL SSUM
120	REAL time_0 , finvmaold(ip1jmp1,llm)
121
122	cym LOGICAL lafin
123	LOGICAL :: lafin=.false.
124	INTEGER ij,iq,l
125	INTEGER ik
126
127	real time_step, t_wrt, t_ops
128
129	! REAL rdayvrai,rdaym_ini
130	! jD_cur: jour julien courant
131	! jH_cur: heure julienne courante
132	REAL :: jD_cur, jH_cur
133	INTEGER :: an, mois, jour
134	REAL :: secondes
135
136	LOGICAL first,callinigrads
137	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
138	save first
139	data first/.true./
140	real dt_cum
141	character*10 infile
142	integer zan, tau0, thoriid
143	integer nid_ctesGCM
144	save nid_ctesGCM
145	real degres
146	real rlong(iip1), rlatg(jjp1)
147	real zx_tmp_2d(iip1,jjp1)
148	integer ndex2d(iip1*jjp1)
149	logical ok_sync
150	parameter (ok_sync = .true.)
151
152	data callinigrads/.true./
153	character*10 string10
154
155	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
156	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
157
158	c+jld variables test conservation energie
159	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
160	C Tendance de la temp. potentiel d (theta)/ d t due a la
161	C tansformation d'energie cinetique en energie thermique
162	C cree par la dissipation
163	REAL dtetaecdt(ip1jmp1,llm)
164	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
165	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
166	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
167	CHARACTER*15 ztit
168	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
169	!IM SAVE ip_ebil_dyn
170	!IM DATA ip_ebil_dyn/0/
171	c-jld
172
173	character*80 dynhist_file, dynhistave_file
174	character(len=20) :: modname
175	character*80 abort_message
176
177	logical dissip_conservative
178	save dissip_conservative
179	data dissip_conservative/.true./
180
181	LOGICAL prem
182	save prem
183	DATA prem/.true./
184	INTEGER testita
185	PARAMETER (testita = 9)
186
187	logical , parameter :: flag_verif = .false.
188
189
190	integer itau_w ! pas de temps ecriture = itap + itau_phy
191
192
193	itaufin = nday*day_step
194	itaufinp1 = itaufin +1
195	modname="leapfrog"
196
197
198	itau = 0
199	c$$$ iday = day_ini+itau/day_step
200	c$$$ time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
201	c$$$ IF(time.GT.1.) THEN
202	c$$$ time = time-1.
203	c$$$ iday = iday+1
204	c$$$ ENDIF
205
206
207	c-----------------------------------------------------------------------
208	c On initialise la pression et la fonction d'Exner :
209	c --------------------------------------------------
210
211	dq=0.
212	CALL pression ( ip1jmp1, ap, bp, ps, p )
213	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
214
215	c-----------------------------------------------------------------------
216	c Debut de l'integration temporelle:
217	c ----------------------------------
218
219	1 CONTINUE
220
221	jD_cur = jD_ref + (day_ini - 1) + int (itau * dtvr / daysec)
222	jH_cur = jH_ref + &
223	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
224
225
226	#ifdef CPP_IOIPSL
227	if (ok_guide) then
228	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
229	endif
230	#endif
231
232
233	c
234	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
235	c CALL test_period ( ucov,vcov,teta,q,p,phis )
236	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
237	c ENDIF
238	c
239
240	! Save fields obtained at previous time step as '...m1'
241	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
242	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
243	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
244	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
245	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
246
247	forward = .TRUE.
248	leapf = .FALSE.
249	dt = dtvr
250
251	c ... P.Le Van .26/04/94 ....
252
253	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
254	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
255
256	! Ehouarn: what is this for? zqmin & zqmax are not used anyway ...
257	! call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
258
259	2 CONTINUE
260
261	c-----------------------------------------------------------------------
262
263	c date:
264	c -----
265
266
267	c gestion des appels de la physique et des dissipations:
268	c ------------------------------------------------------
269	c
270	c ... P.Le Van ( 6/02/95 ) ....
271
272	apphys = .FALSE.
273	statcl = .FALSE.
274	conser = .FALSE.
275	apdiss = .FALSE.
276
277	IF( purmats ) THEN
278	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
279	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
280	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
281	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
282	ELSE
283	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
284	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
285	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
286	END IF
287
288	c-----------------------------------------------------------------------
289	c calcul des tendances dynamiques:
290	c --------------------------------
291
292	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
293
294	time = jD_cur + jH_cur
295	CALL caldyn
296	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
297	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
298
299
300	c-----------------------------------------------------------------------
301	c calcul des tendances advection des traceurs (dont l'humidite)
302	c -------------------------------------------------------------
303
304	IF( forward. OR . leapf ) THEN
305
306	CALL caladvtrac(q,pbaru,pbarv,
307	* p, masse, dq, teta,
308	. flxw, pk)
309
310	IF (offline) THEN
311	Cmaf stokage du flux de masse pour traceurs OFF-LINE
312
313	#ifdef CPP_IOIPSL
314	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
315	. dtvr, itau)
316	#endif
317
318
319	ENDIF ! of IF (offline)
320	c
321	ENDIF ! of IF( forward. OR . leapf )
322
323
324	c-----------------------------------------------------------------------
325	c integrations dynamique et traceurs:
326	c ----------------------------------
327
328
329	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
330	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
331	$ finvmaold )
332
333
334	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
335	c
336	c-----------------------------------------------------------------------
337	c calcul des tendances physiques:
338	c -------------------------------
339	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
340	c
341	IF( purmats ) THEN
342	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
343	ELSE
344	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
345	ENDIF
346	c
347	c
348	IF( apphys ) THEN
349	c
350	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
351	c
352
353	CALL pression ( ip1jmp1, ap, bp, ps, p )
354	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
355
356	! rdaym_ini = itau * dtvr / daysec
357	! rdayvrai = rdaym_ini + day_ini
358	jD_cur = jD_ref + (day_ini - 1) + int (itau * dtvr / daysec)
359	jH_cur = jH_ref + &
360	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
361	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
362	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
363	! write(lunout,*)'current date = ',an, mois, jour, secondes
364
365	c rajout debug
366	c lafin = .true.
367
368
369	c Inbterface avec les routines de phylmd (phymars ... )
370	c -----------------------------------------------------
371
372	c+jld
373
374	c Diagnostique de conservation de l'énergie : initialisation
375	IF (ip_ebil_dyn.ge.1 ) THEN
376	ztit='bil dyn'
377	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
378	IF (planet_type.eq."earth") THEN
379	CALL diagedyn(ztit,2,1,1,dtphys
380	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
381	ENDIF
382	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
383	c-jld
384	#ifdef CPP_IOIPSL
385	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
386	IF (first) THEN
387	first=.false.
388	#include "ini_paramLMDZ_dyn.h"
389	ENDIF
390	c
391	#include "write_paramLMDZ_dyn.h"
392	c
393	#endif
394	! #endif of #ifdef CPP_IOIPSL
395	CALL calfis( lafin , jD_cur, jH_cur,
396	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
397	$ du,dv,dteta,dq,
398	$ flxw,
399	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
400
401	IF (ok_strato) THEN
402	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
403	ENDIF
404
405	c ajout des tendances physiques:
406	c ------------------------------
407	CALL addfi( dtphys, leapf, forward ,
408	$ ucov, vcov, teta , q ,ps ,
409	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
410	c
411	c Diagnostique de conservation de l'énergie : difference
412	IF (ip_ebil_dyn.ge.1 ) THEN
413	ztit='bil phys'
414	IF (planet_type.eq."earth") THEN
415	CALL diagedyn(ztit,2,1,1,dtphys
416	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
417	ENDIF
418	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
419
420	ENDIF ! of IF( apphys )
421
422	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
423	c Calcul academique de la physique = Rappel Newtonien + friction
424	c --------------------------------------------------------------
425	teta(:,:)=teta(:,:)
426	s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
427	call friction(ucov,vcov,iphysiq*dtvr)
428	ENDIF
429
430
431	c-jld
432
433	CALL pression ( ip1jmp1, ap, bp, ps, p )
434	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
435
436
437	c-----------------------------------------------------------------------
438	c dissipation horizontale et verticale des petites echelles:
439	c ----------------------------------------------------------
440
441	IF(apdiss) THEN
442
443
444	c calcul de l'energie cinetique avant dissipation
445	call covcont(llm,ucov,vcov,ucont,vcont)
446	call enercin(vcov,ucov,vcont,ucont,ecin0)
447
448	c dissipation
449	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
450	ucov=ucov+dudis
451	vcov=vcov+dvdis
452	c teta=teta+dtetadis
453
454
455	c------------------------------------------------------------------------
456	if (dissip_conservative) then
457	C On rajoute la tendance due a la transform. Ec -> E therm. cree
458	C lors de la dissipation
459	call covcont(llm,ucov,vcov,ucont,vcont)
460	call enercin(vcov,ucov,vcont,ucont,ecin)
461	dtetaecdt= (ecin0-ecin)/ pk
462	c teta=teta+dtetaecdt
463	dtetadis=dtetadis+dtetaecdt
464	endif
465	teta=teta+dtetadis
466	c------------------------------------------------------------------------
467
468
469	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
470	c ... Calcul de la valeur moyenne, unique de h aux poles .....
471	c
472
473	DO l = 1, llm
474	DO ij = 1,iim
475	tppn(ij) = aire( ij ) * teta( ij ,l)
476	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
477	ENDDO
478	tpn = SSUM(iim,tppn,1)/apoln
479	tps = SSUM(iim,tpps,1)/apols
480
481	DO ij = 1, iip1
482	teta( ij ,l) = tpn
483	teta(ij+ip1jm,l) = tps
484	ENDDO
485	ENDDO
486
487	DO ij = 1,iim
488	tppn(ij) = aire( ij ) * ps ( ij )
489	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
490	ENDDO
491	tpn = SSUM(iim,tppn,1)/apoln
492	tps = SSUM(iim,tpps,1)/apols
493
494	DO ij = 1, iip1
495	ps( ij ) = tpn
496	ps(ij+ip1jm) = tps
497	ENDDO
498
499
500	END IF ! of IF(apdiss)
501
502	c ajout debug
503	c IF( lafin ) then
504	c abort_message = 'Simulation finished'
505	c call abort_gcm(modname,abort_message,0)
506	c ENDIF
507
508	c ********************************************************************
509	c ********************************************************************
510	c .... fin de l'integration dynamique et physique pour le pas itau ..
511	c ********************************************************************
512	c ********************************************************************
513
514	c preparation du pas d'integration suivant ......
515
516	IF ( .NOT.purmats ) THEN
517	c ........................................................
518	c .............. schema matsuno + leapfrog ..............
519	c ........................................................
520
521	IF(forward. OR. leapf) THEN
522	itau= itau + 1
523	c$$$ iday= day_ini+itau/day_step
524	c$$$ time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
525	c$$$ IF(time.GT.1.) THEN
526	c$$$ time = time-1.
527	c$$$ iday = iday+1
528	c$$$ ENDIF
529	ENDIF
530
531
532	IF( itau. EQ. itaufinp1 ) then
533	if (flag_verif) then
534	write(79,*) 'ucov',ucov
535	write(80,*) 'vcov',vcov
536	write(81,*) 'teta',teta
537	write(82,*) 'ps',ps
538	write(83,*) 'q',q
539	WRITE(85,*) 'q1 = ',q(:,:,1)
540	WRITE(86,*) 'q3 = ',q(:,:,3)
541	endif
542
543	abort_message = 'Simulation finished'
544
545	call abort_gcm(modname,abort_message,0)
546	ENDIF
547	c-----------------------------------------------------------------------
548	c ecriture du fichier histoire moyenne:
549	c -------------------------------------
550
551	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
552	IF(itau.EQ.itaufin) THEN
553	iav=1
554	ELSE
555	iav=0
556	ENDIF
557
558	IF (ok_dynzon) THEN
559	#ifdef CPP_IOIPSL
560	CALL writedynav(histaveid, itau,vcov ,
561	, ucov,teta,pk,phi,q,masse,ps,phis)
562	CALL bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
563	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
564	#endif
565	END IF
566
567	ENDIF
568
569	c-----------------------------------------------------------------------
570	c ecriture de la bande histoire:
571	c ------------------------------
572
573	IF( MOD(itau,iecri ).EQ.0) THEN
574	c IF( MOD(itau,iecri*day_step).EQ.0) THEN
575
576	nbetat = nbetatdem
577	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
578	unat=0.
579	do l=1,llm
580	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
581	vnat(:,l)=vcov(:,l)/cv(:)
582	enddo
583	#ifdef CPP_IOIPSL
584	c CALL writehist(histid,histvid,itau,vcov,
585	c & ucov,teta,phi,q,masse,ps,phis)
586	#endif
587	! For some Grads outputs of fields
588	if (output_grads_dyn) then
589	#include "write_grads_dyn.h"
590	endif
591
592	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
593
594	IF(itau.EQ.itaufin) THEN
595
596
597	if (planet_type.eq."earth") then
598	#ifdef CPP_EARTH
599	! Write an Earth-format restart file
600	CALL dynredem1("restart.nc",0.0,
601	& vcov,ucov,teta,q,masse,ps)
602	#endif
603	endif ! of if (planet_type.eq."earth")
604
605	CLOSE(99)
606	ENDIF ! of IF (itau.EQ.itaufin)
607
608	c-----------------------------------------------------------------------
609	c gestion de l'integration temporelle:
610	c ------------------------------------
611
612	IF( MOD(itau,iperiod).EQ.0 ) THEN
613	GO TO 1
614	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
615
616	IF( forward ) THEN
617	c fin du pas forward et debut du pas backward
618
619	forward = .FALSE.
620	leapf = .FALSE.
621	GO TO 2
622
623	ELSE
624	c fin du pas backward et debut du premier pas leapfrog
625
626	leapf = .TRUE.
627	dt = 2.*dtvr
628	GO TO 2
629	END IF ! of IF (forward)
630	ELSE
631
632	c ...... pas leapfrog .....
633
634	leapf = .TRUE.
635	dt = 2.*dtvr
636	GO TO 2
637	END IF ! of IF (MOD(itau,iperiod).EQ.0)
638	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
639
640	ELSE ! of IF (.not.purmats)
641
642	c ........................................................
643	c .............. schema matsuno ...............
644	c ........................................................
645	IF( forward ) THEN
646
647	itau = itau + 1
648	c$$$ iday = day_ini+itau/day_step
649	c$$$ time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
650	c$$$
651	c$$$ IF(time.GT.1.) THEN
652	c$$$ time = time-1.
653	c$$$ iday = iday+1
654	c$$$ ENDIF
655
656	forward = .FALSE.
657	IF( itau. EQ. itaufinp1 ) then
658	abort_message = 'Simulation finished'
659	call abort_gcm(modname,abort_message,0)
660	ENDIF
661	GO TO 2
662
663	ELSE ! of IF(forward)
664
665	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
666	IF(itau.EQ.itaufin) THEN
667	iav=1
668	ELSE
669	iav=0
670	ENDIF
671
672	IF (ok_dynzon) THEN
673	#ifdef CPP_IOIPSL
674	CALL writedynav(histaveid, itau,vcov ,
675	, ucov,teta,pk,phi,q,masse,ps,phis)
676	CALL bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
677	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
678	#endif
679	END IF
680
681	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
682
683	IF(MOD(itau,iecri ).EQ.0) THEN
684	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
685	nbetat = nbetatdem
686	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
687	unat=0.
688	do l=1,llm
689	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
690	vnat(:,l)=vcov(:,l)/cv(:)
691	enddo
692	#ifdef CPP_IOIPSL
693	c CALL writehist( histid, histvid, itau,vcov ,
694	c & ucov,teta,phi,q,masse,ps,phis)
695	#endif
696	! For some Grads outputs
697	if (output_grads_dyn) then
698	#include "write_grads_dyn.h"
699	endif
700
701	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
702
703	IF(itau.EQ.itaufin) THEN
704	if (planet_type.eq."earth") then
705	#ifdef CPP_EARTH
706	CALL dynredem1("restart.nc",0.0,
707	& vcov,ucov,teta,q,masse,ps)
708	#endif
709	endif ! of if (planet_type.eq."earth")
710	ENDIF ! of IF(itau.EQ.itaufin)
711
712	forward = .TRUE.
713	GO TO 1
714
715	ENDIF ! of IF (forward)
716
717	END IF ! of IF(.not.purmats)
718
719	STOP
720	END

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

Download in other formats:

Original Format