Context Navigation

leapfrog.F @ 5425

Last change on this file since 5425 was 1717, checked in by Ehouarn Millour, 12 years ago
Added "arch" files for Ada (using dynamic libraries for NetCDF, you must have export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/smplocal/pub/NetCDF/4.1.3/lib:/smplocal/pub/HDF5/1.8.9/seq/lib in your .bashrc or .bash_login or in your job to run). Also updated some sources so that gcm bench runs in "debug" mode (note that all these changes are minor and have already been implemented in LMDZ5 trunk). EM
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 1717 2013-01-25 08:26:03Z jyg $
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 - day_ref + 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 - day_ref
359	$ + int (itau * dtvr / daysec)
360	jH_cur = jH_ref + &
361	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
362	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
363	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
364	! write(lunout,*)'current date = ',an, mois, jour, secondes
365
366	c rajout debug
367	c lafin = .true.
368
369
370	c Inbterface avec les routines de phylmd (phymars ... )
371	c -----------------------------------------------------
372
373	c+jld
374
375	c Diagnostique de conservation de l'énergie : initialisation
376	IF (ip_ebil_dyn.ge.1 ) THEN
377	ztit='bil dyn'
378	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
379	IF (planet_type.eq."earth") THEN
380	CALL diagedyn(ztit,2,1,1,dtphys
381	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
382	ENDIF
383	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
384	c-jld
385	#ifdef CPP_IOIPSL
386	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
387	IF (first) THEN
388	first=.false.
389	#include "ini_paramLMDZ_dyn.h"
390	ENDIF
391	c
392	#include "write_paramLMDZ_dyn.h"
393	c
394	#endif
395	! #endif of #ifdef CPP_IOIPSL
396	CALL calfis( lafin , jD_cur, jH_cur,
397	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
398	$ du,dv,dteta,dq,
399	$ flxw,
400	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
401
402	IF (ok_strato) THEN
403	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
404	ENDIF
405
406	c ajout des tendances physiques:
407	c ------------------------------
408	CALL addfi( dtphys, leapf, forward ,
409	$ ucov, vcov, teta , q ,ps ,
410	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
411	c
412	c Diagnostique de conservation de l'énergie : difference
413	IF (ip_ebil_dyn.ge.1 ) THEN
414	ztit='bil phys'
415	IF (planet_type.eq."earth") THEN
416	CALL diagedyn(ztit,2,1,1,dtphys
417	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
418	ENDIF
419	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
420
421	ENDIF ! of IF( apphys )
422
423	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
424	c Calcul academique de la physique = Rappel Newtonien + friction
425	c --------------------------------------------------------------
426	teta(:,:)=teta(:,:)
427	s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
428	call friction(ucov,vcov,iphysiq*dtvr)
429	ENDIF
430
431
432	c-jld
433
434	CALL pression ( ip1jmp1, ap, bp, ps, p )
435	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
436
437
438	c-----------------------------------------------------------------------
439	c dissipation horizontale et verticale des petites echelles:
440	c ----------------------------------------------------------
441
442	IF(apdiss) THEN
443
444
445	c calcul de l'energie cinetique avant dissipation
446	call covcont(llm,ucov,vcov,ucont,vcont)
447	call enercin(vcov,ucov,vcont,ucont,ecin0)
448
449	c dissipation
450	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
451	ucov=ucov+dudis
452	vcov=vcov+dvdis
453	c teta=teta+dtetadis
454
455
456	c------------------------------------------------------------------------
457	if (dissip_conservative) then
458	C On rajoute la tendance due a la transform. Ec -> E therm. cree
459	C lors de la dissipation
460	call covcont(llm,ucov,vcov,ucont,vcont)
461	call enercin(vcov,ucov,vcont,ucont,ecin)
462	dtetaecdt= (ecin0-ecin)/ pk
463	c teta=teta+dtetaecdt
464	dtetadis=dtetadis+dtetaecdt
465	endif
466	teta=teta+dtetadis
467	c------------------------------------------------------------------------
468
469
470	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
471	c ... Calcul de la valeur moyenne, unique de h aux poles .....
472	c
473
474	DO l = 1, llm
475	DO ij = 1,iim
476	tppn(ij) = aire( ij ) * teta( ij ,l)
477	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
478	ENDDO
479	tpn = SSUM(iim,tppn,1)/apoln
480	tps = SSUM(iim,tpps,1)/apols
481
482	DO ij = 1, iip1
483	teta( ij ,l) = tpn
484	teta(ij+ip1jm,l) = tps
485	ENDDO
486	ENDDO
487
488	DO ij = 1,iim
489	tppn(ij) = aire( ij ) * ps ( ij )
490	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
491	ENDDO
492	tpn = SSUM(iim,tppn,1)/apoln
493	tps = SSUM(iim,tpps,1)/apols
494
495	DO ij = 1, iip1
496	ps( ij ) = tpn
497	ps(ij+ip1jm) = tps
498	ENDDO
499
500
501	END IF ! of IF(apdiss)
502
503	c ajout debug
504	c IF( lafin ) then
505	c abort_message = 'Simulation finished'
506	c call abort_gcm(modname,abort_message,0)
507	c ENDIF
508
509	c ********************************************************************
510	c ********************************************************************
511	c .... fin de l'integration dynamique et physique pour le pas itau ..
512	c ********************************************************************
513	c ********************************************************************
514
515	c preparation du pas d'integration suivant ......
516
517	IF ( .NOT.purmats ) THEN
518	c ........................................................
519	c .............. schema matsuno + leapfrog ..............
520	c ........................................................
521
522	IF(forward. OR. leapf) THEN
523	itau= itau + 1
524	c$$$ iday= day_ini+itau/day_step
525	c$$$ time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
526	c$$$ IF(time.GT.1.) THEN
527	c$$$ time = time-1.
528	c$$$ iday = iday+1
529	c$$$ ENDIF
530	ENDIF
531
532
533	IF( itau. EQ. itaufinp1 ) then
534	if (flag_verif) then
535	write(79,*) 'ucov',ucov
536	write(80,*) 'vcov',vcov
537	write(81,*) 'teta',teta
538	write(82,*) 'ps',ps
539	write(83,*) 'q',q
540	WRITE(85,*) 'q1 = ',q(:,:,1)
541	WRITE(86,*) 'q3 = ',q(:,:,3)
542	endif
543
544	abort_message = 'Simulation finished'
545
546	call abort_gcm(modname,abort_message,0)
547	ENDIF
548	c-----------------------------------------------------------------------
549	c ecriture du fichier histoire moyenne:
550	c -------------------------------------
551
552	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
553	IF(itau.EQ.itaufin) THEN
554	iav=1
555	ELSE
556	iav=0
557	ENDIF
558
559	IF (ok_dynzon) THEN
560	#ifdef CPP_IOIPSL
561	! CALL writedynav(histaveid, itau,vcov ,
562	! , ucov,teta,pk,phi,q,masse,ps,phis)
563	CALL bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
564	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
565	#endif
566	END IF
567
568	ENDIF
569
570	c-----------------------------------------------------------------------
571	c ecriture de la bande histoire:
572	c ------------------------------
573
574	IF( MOD(itau,iecri ).EQ.0) THEN
575	c IF( MOD(itau,iecri*day_step).EQ.0) THEN
576
577	! nbetat = nbetatdem
578	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
579	unat=0.
580	do l=1,llm
581	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
582	vnat(:,l)=vcov(:,l)/cv(:)
583	enddo
584	#ifdef CPP_IOIPSL
585	c CALL writehist(histid,histvid,itau,vcov,
586	c & ucov,teta,phi,q,masse,ps,phis)
587	#endif
588	! For some Grads outputs of fields
589	if (output_grads_dyn) then
590	#include "write_grads_dyn.h"
591	endif
592
593	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
594
595	IF(itau.EQ.itaufin) THEN
596
597
598	if (planet_type.eq."earth") then
599	! Write an Earth-format restart file
600	CALL dynredem1("restart.nc",0.0,
601	& vcov,ucov,teta,q,masse,ps)
602	endif ! of if (planet_type.eq."earth")
603
604	CLOSE(99)
605	ENDIF ! of IF (itau.EQ.itaufin)
606
607	c-----------------------------------------------------------------------
608	c gestion de l'integration temporelle:
609	c ------------------------------------
610
611	IF( MOD(itau,iperiod).EQ.0 ) THEN
612	GO TO 1
613	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
614
615	IF( forward ) THEN
616	c fin du pas forward et debut du pas backward
617
618	forward = .FALSE.
619	leapf = .FALSE.
620	GO TO 2
621
622	ELSE
623	c fin du pas backward et debut du premier pas leapfrog
624
625	leapf = .TRUE.
626	dt = 2.*dtvr
627	GO TO 2
628	END IF ! of IF (forward)
629	ELSE
630
631	c ...... pas leapfrog .....
632
633	leapf = .TRUE.
634	dt = 2.*dtvr
635	GO TO 2
636	END IF ! of IF (MOD(itau,iperiod).EQ.0)
637	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
638
639	ELSE ! of IF (.not.purmats)
640
641	c ........................................................
642	c .............. schema matsuno ...............
643	c ........................................................
644	IF( forward ) THEN
645
646	itau = itau + 1
647	c$$$ iday = day_ini+itau/day_step
648	c$$$ time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
649	c$$$
650	c$$$ IF(time.GT.1.) THEN
651	c$$$ time = time-1.
652	c$$$ iday = iday+1
653	c$$$ ENDIF
654
655	forward = .FALSE.
656	IF( itau. EQ. itaufinp1 ) then
657	abort_message = 'Simulation finished'
658	call abort_gcm(modname,abort_message,0)
659	ENDIF
660	GO TO 2
661
662	ELSE ! of IF(forward)
663
664	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
665	IF(itau.EQ.itaufin) THEN
666	iav=1
667	ELSE
668	iav=0
669	ENDIF
670
671	IF (ok_dynzon) THEN
672	#ifdef CPP_IOIPSL
673	! CALL writedynav(histaveid, itau,vcov ,
674	! , ucov,teta,pk,phi,q,masse,ps,phis)
675	CALL bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
676	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
677	#endif
678	END IF
679
680	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
681
682	IF(MOD(itau,iecri ).EQ.0) THEN
683	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
684	! nbetat = nbetatdem
685	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
686	unat=0.
687	do l=1,llm
688	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
689	vnat(:,l)=vcov(:,l)/cv(:)
690	enddo
691	#ifdef CPP_IOIPSL
692	c CALL writehist( histid, histvid, itau,vcov ,
693	c & ucov,teta,phi,q,masse,ps,phis)
694	#endif
695	! For some Grads outputs
696	if (output_grads_dyn) then
697	#include "write_grads_dyn.h"
698	endif
699
700	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
701
702	IF(itau.EQ.itaufin) THEN
703	if (planet_type.eq."earth") then
704	CALL dynredem1("restart.nc",0.0,
705	& vcov,ucov,teta,q,masse,ps)
706	endif ! of if (planet_type.eq."earth")
707	ENDIF ! of IF(itau.EQ.itaufin)
708
709	forward = .TRUE.
710	GO TO 1
711
712	ENDIF ! of IF (forward)
713
714	END IF ! of IF(.not.purmats)
715
716	STOP
717	END

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

Download in other formats:

Original Format