Context Navigation

leapfrog.F @ 3285

Last change on this file since 3285 was 1676, checked in by aslmd, 12 years ago

Repare le probleme "un jour sans fin" des versions precedentes pour le modele generique. Pas d'impact sur modele terrestre.

NB: changements egalement faits sur dyn3dmem en prevision de son utilisation

Fix the bug "Groundhog Day" which was a problem for previous versions using generic physics. No impact whatsoever on terrestrial model.

NB: changes also done on dyn3dmem which is planned to be used soon

Author : AS equipe planeto

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 25.4 KB

Line
1	!
2	! $Id: leapfrog.F 1676 2012-11-08 20:03:30Z 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	USE control_mod
18	IMPLICIT NONE
19
20	c ...... Version du 10/01/98 ..........
21
22	c avec coordonnees verticales hybrides
23	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
24
25	c=======================================================================
26	c
27	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
28	c -------
29	c
30	c Objet:
31	c ------
32	c
33	c GCM LMD nouvelle grille
34	c
35	c=======================================================================
36	c
37	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
38	c et possibilite d'appeler une fonction f(y) a derivee tangente
39	c hyperbolique a la place de la fonction a derivee sinusoidale.
40
41	c ... Possibilite de choisir le shema pour l'advection de
42	c q , en modifiant iadv dans traceur.def (10/02) .
43	c
44	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
45	c Pour Van-Leer iadv=10
46	c
47	c-----------------------------------------------------------------------
48	c Declarations:
49	c -------------
50
51	#include "dimensions.h"
52	#include "paramet.h"
53	#include "comconst.h"
54	#include "comdissnew.h"
55	#include "comvert.h"
56	#include "comgeom.h"
57	#include "logic.h"
58	#include "temps.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
75	c variables dynamiques
76	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
77	REAL teta(ip1jmp1,llm) ! temperature potentielle
78	REAL q(ip1jmp1,llm,nqtot) ! champs advectes
79	REAL ps(ip1jmp1) ! pression au sol
80	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
81	REAL pks(ip1jmp1) ! exner au sol
82	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
83	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
84	REAL masse(ip1jmp1,llm) ! masse d'air
85	REAL phis(ip1jmp1) ! geopotentiel au sol
86	REAL phi(ip1jmp1,llm) ! geopotentiel
87	REAL w(ip1jmp1,llm) ! vitesse verticale
88
89	c variables dynamiques intermediaire pour le transport
90	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
91
92	c variables dynamiques au pas -1
93	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
94	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
95	REAL massem1(ip1jmp1,llm)
96
97	c tendances dynamiques
98	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
99	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
100
101	c tendances de la dissipation
102	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
103	REAL dtetadis(ip1jmp1,llm)
104
105	c tendances physiques
106	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
107	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
108
109	c variables pour le fichier histoire
110	REAL dtav ! intervalle de temps elementaire
111
112	REAL tppn(iim),tpps(iim),tpn,tps
113	c
114	INTEGER itau,itaufinp1,iav
115	! INTEGER iday ! jour julien
116	REAL time
117
118	REAL SSUM
119	REAL time_0
120	! REAL 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	logical physic
152
153	data callinigrads/.true./
154	character*10 string10
155
156	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
157	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
158
159	c+jld variables test conservation energie
160	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
161	C Tendance de la temp. potentiel d (theta)/ d t due a la
162	C tansformation d'energie cinetique en energie thermique
163	C cree par la dissipation
164	REAL dtetaecdt(ip1jmp1,llm)
165	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
166	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
167	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
168	CHARACTER*15 ztit
169	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
170	!IM SAVE ip_ebil_dyn
171	!IM DATA ip_ebil_dyn/0/
172	c-jld
173
174	character*80 dynhist_file, dynhistave_file
175	character(len=*),parameter :: modname="leapfrog"
176	character*80 abort_message
177
178	logical dissip_conservative
179	save dissip_conservative
180	data dissip_conservative/.true./
181
182	LOGICAL prem
183	save prem
184	DATA prem/.true./
185	INTEGER testita
186	PARAMETER (testita = 9)
187
188	logical , parameter :: flag_verif = .false.
189
190
191	integer itau_w ! pas de temps ecriture = itap + itau_phy
192
193
194	itaufin = nday*day_step
195	itaufinp1 = itaufin +1
196	itau = 0
197	physic=.true.
198	if (iflag_phys==0.or.iflag_phys==2) physic=.false.
199
200	c iday = day_ini+itau/day_step
201	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
202	c IF(time.GT.1.) THEN
203	c time = time-1.
204	c iday = iday+1
205	c ENDIF
206
207
208	c-----------------------------------------------------------------------
209	c On initialise la pression et la fonction d'Exner :
210	c --------------------------------------------------
211
212	dq(:,:,:)=0.
213	CALL pression ( ip1jmp1, ap, bp, ps, p )
214	if (pressure_exner) then
215	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
216	else
217	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
218	endif
219
220	c-----------------------------------------------------------------------
221	c Debut de l'integration temporelle:
222	c ----------------------------------
223
224	1 CONTINUE ! Matsuno Forward step begins here
225
226	jD_cur = jD_ref + day_ini - day_ref + &
227	& itau/day_step
228	jH_cur = jH_ref + start_time + &
229	& mod(itau,day_step)/float(day_step)
230	jD_cur = jD_cur + int(jH_cur)
231	jH_cur = jH_cur - int(jH_cur)
232
233
234	#ifdef CPP_IOIPSL
235	if (ok_guide) then
236	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
237	endif
238	#endif
239
240
241	c
242	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
243	c CALL test_period ( ucov,vcov,teta,q,p,phis )
244	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
245	c ENDIF
246	c
247
248	! Save fields obtained at previous time step as '...m1'
249	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
250	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
251	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
252	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
253	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
254
255	forward = .TRUE.
256	leapf = .FALSE.
257	dt = dtvr
258
259	c ... P.Le Van .26/04/94 ....
260	! Ehouarn: finvmaold is actually not used
261	! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
262	! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
263
264	2 CONTINUE ! Matsuno backward or leapfrog step begins here
265
266	c-----------------------------------------------------------------------
267
268	c date:
269	c -----
270
271
272	c gestion des appels de la physique et des dissipations:
273	c ------------------------------------------------------
274	c
275	c ... P.Le Van ( 6/02/95 ) ....
276
277	apphys = .FALSE.
278	statcl = .FALSE.
279	conser = .FALSE.
280	apdiss = .FALSE.
281
282	IF( purmats ) THEN
283	! Purely Matsuno time stepping
284	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
285	IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
286	s apdiss = .TRUE.
287	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
288	s .and. physic ) apphys = .TRUE.
289	ELSE
290	! Leapfrog/Matsuno time stepping
291	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
292	IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
293	s apdiss = .TRUE.
294	IF( MOD(itau+1,iphysiq).EQ.0.AND.physic ) apphys=.TRUE.
295	END IF
296
297	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
298	! supress dissipation step
299	if (llm.eq.1) then
300	apdiss=.false.
301	endif
302
303	c-----------------------------------------------------------------------
304	c calcul des tendances dynamiques:
305	c --------------------------------
306
307	! compute geopotential phi()
308	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
309
310	time = jD_cur + jH_cur
311	CALL caldyn
312	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
313	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
314
315
316	c-----------------------------------------------------------------------
317	c calcul des tendances advection des traceurs (dont l'humidite)
318	c -------------------------------------------------------------
319
320	IF( forward. OR . leapf ) THEN
321	! Ehouarn: NB: at this point p with ps are not synchronized
322	! (whereas mass and ps are...)
323	CALL caladvtrac(q,pbaru,pbarv,
324	* p, masse, dq, teta,
325	. flxw, pk)
326
327	IF (offline) THEN
328	Cmaf stokage du flux de masse pour traceurs OFF-LINE
329
330	#ifdef CPP_IOIPSL
331	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
332	. dtvr, itau)
333	#endif
334
335
336	ENDIF ! of IF (offline)
337	c
338	ENDIF ! of IF( forward. OR . leapf )
339
340
341	c-----------------------------------------------------------------------
342	c integrations dynamique et traceurs:
343	c ----------------------------------
344
345
346	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
347	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
348	! $ finvmaold )
349
350
351	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
352	c
353	c-----------------------------------------------------------------------
354	c calcul des tendances physiques:
355	c -------------------------------
356	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
357	c
358	IF( purmats ) THEN
359	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
360	ELSE
361	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
362	ENDIF
363	c
364	c
365	IF( apphys ) THEN
366	c
367	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
368	c
369
370	CALL pression ( ip1jmp1, ap, bp, ps, p )
371	if (pressure_exner) then
372	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
373	else
374	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
375	endif
376
377	! rdaym_ini = itau * dtvr / daysec
378	! rdayvrai = rdaym_ini + day_ini
379	! jD_cur = jD_ref + day_ini - day_ref
380	! $ + int (itau * dtvr / daysec)
381	! jH_cur = jH_ref + &
382	! & (itau * dtvr / daysec - int(itau * dtvr / daysec))
383	jD_cur = jD_ref + day_ini - day_ref + &
384	& itau/day_step
385
386	IF (planet_type .eq."generic") THEN
387	! AS: we make jD_cur to be pday
388	jD_cur = int(day_ini + itau/day_step)
389	ENDIF
390
391	jH_cur = jH_ref + start_time + &
392	& mod(itau,day_step)/float(day_step)
393	jD_cur = jD_cur + int(jH_cur)
394	jH_cur = jH_cur - int(jH_cur)
395	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
396	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
397	! write(lunout,*)'current date = ',an, mois, jour, secondes
398
399	c rajout debug
400	c lafin = .true.
401
402
403	c Inbterface avec les routines de phylmd (phymars ... )
404	c -----------------------------------------------------
405
406	c+jld
407
408	c Diagnostique de conservation de l'énergie : initialisation
409	IF (ip_ebil_dyn.ge.1 ) THEN
410	ztit='bil dyn'
411	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
412	IF (planet_type.eq."earth") THEN
413	#ifdef CPP_EARTH
414	CALL diagedyn(ztit,2,1,1,dtphys
415	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
416	#endif
417	ENDIF
418	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
419	c-jld
420	#ifdef CPP_IOIPSL
421	cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
422	cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
423	c IF (first) THEN
424	c first=.false.
425	c#include "ini_paramLMDZ_dyn.h"
426	c ENDIF
427	c
428	c#include "write_paramLMDZ_dyn.h"
429	c
430	#endif
431	! #endif of #ifdef CPP_IOIPSL
432	CALL calfis( lafin , jD_cur, jH_cur,
433	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
434	$ du,dv,dteta,dq,
435	$ flxw,
436	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
437
438	IF (ok_strato) THEN
439	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
440	ENDIF
441
442	c ajout des tendances physiques:
443	c ------------------------------
444	CALL addfi( dtphys, leapf, forward ,
445	$ ucov, vcov, teta , q ,ps ,
446	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
447	c
448	c Diagnostique de conservation de l'énergie : difference
449	IF (ip_ebil_dyn.ge.1 ) THEN
450	ztit='bil phys'
451	IF (planet_type.eq."earth") THEN
452	CALL diagedyn(ztit,2,1,1,dtphys
453	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
454	ENDIF
455	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
456
457	ENDIF ! of IF( apphys )
458
459	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
460	! Academic case : Simple friction and Newtonan relaxation
461	! -------------------------------------------------------
462	DO l=1,llm
463	DO ij=1,ip1jmp1
464	teta(ij,l)=teta(ij,l)-dtvr*
465	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
466	ENDDO
467	ENDDO ! of DO l=1,llm
468
469	if (planet_type.eq."giant") then
470	! add an intrinsic heat flux at the base of the atmosphere
471	teta(:,1)=teta(:,1)+dtvraire(:)ihf/cpp/masse(:,1)
472	endif
473
474	call friction(ucov,vcov,dtvr)
475
476	! Sponge layer (if any)
477	IF (ok_strato) THEN
478	dufi(:,:)=0.
479	dvfi(:,:)=0.
480	dtetafi(:,:)=0.
481	dqfi(:,:,:)=0.
482	dpfi(:)=0.
483	CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
484	CALL addfi( dtvr, leapf, forward ,
485	$ ucov, vcov, teta , q ,ps ,
486	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
487	ENDIF ! of IF (ok_strato)
488	ENDIF ! of IF (iflag_phys.EQ.2)
489
490
491	c-jld
492
493	CALL pression ( ip1jmp1, ap, bp, ps, p )
494	if (pressure_exner) then
495	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
496	else
497	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
498	endif
499
500
501	c-----------------------------------------------------------------------
502	c dissipation horizontale et verticale des petites echelles:
503	c ----------------------------------------------------------
504
505	IF(apdiss) THEN
506
507
508	c calcul de l'energie cinetique avant dissipation
509	call covcont(llm,ucov,vcov,ucont,vcont)
510	call enercin(vcov,ucov,vcont,ucont,ecin0)
511
512	c dissipation
513	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
514	ucov=ucov+dudis
515	vcov=vcov+dvdis
516	c teta=teta+dtetadis
517
518
519	c------------------------------------------------------------------------
520	if (dissip_conservative) then
521	C On rajoute la tendance due a la transform. Ec -> E therm. cree
522	C lors de la dissipation
523	call covcont(llm,ucov,vcov,ucont,vcont)
524	call enercin(vcov,ucov,vcont,ucont,ecin)
525	dtetaecdt= (ecin0-ecin)/ pk
526	c teta=teta+dtetaecdt
527	dtetadis=dtetadis+dtetaecdt
528	endif
529	teta=teta+dtetadis
530	c------------------------------------------------------------------------
531
532
533	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
534	c ... Calcul de la valeur moyenne, unique de h aux poles .....
535	c
536
537	DO l = 1, llm
538	DO ij = 1,iim
539	tppn(ij) = aire( ij ) * teta( ij ,l)
540	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
541	ENDDO
542	tpn = SSUM(iim,tppn,1)/apoln
543	tps = SSUM(iim,tpps,1)/apols
544
545	DO ij = 1, iip1
546	teta( ij ,l) = tpn
547	teta(ij+ip1jm,l) = tps
548	ENDDO
549	ENDDO
550
551	if (1 == 0) then
552	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
553	!!! 2) should probably not be here anyway
554	!!! but are kept for those who would want to revert to previous behaviour
555	DO ij = 1,iim
556	tppn(ij) = aire( ij ) * ps ( ij )
557	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
558	ENDDO
559	tpn = SSUM(iim,tppn,1)/apoln
560	tps = SSUM(iim,tpps,1)/apols
561
562	DO ij = 1, iip1
563	ps( ij ) = tpn
564	ps(ij+ip1jm) = tps
565	ENDDO
566	endif ! of if (1 == 0)
567
568	END IF ! of IF(apdiss)
569
570	c ajout debug
571	c IF( lafin ) then
572	c abort_message = 'Simulation finished'
573	c call abort_gcm(modname,abort_message,0)
574	c ENDIF
575
576	c ********************************************************************
577	c ********************************************************************
578	c .... fin de l'integration dynamique et physique pour le pas itau ..
579	c ********************************************************************
580	c ********************************************************************
581
582	c preparation du pas d'integration suivant ......
583
584	IF ( .NOT.purmats ) THEN
585	c ........................................................
586	c .............. schema matsuno + leapfrog ..............
587	c ........................................................
588
589	IF(forward. OR. leapf) THEN
590	itau= itau + 1
591	c iday= day_ini+itau/day_step
592	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
593	c IF(time.GT.1.) THEN
594	c time = time-1.
595	c iday = iday+1
596	c ENDIF
597	ENDIF
598
599
600	IF( itau. EQ. itaufinp1 ) then
601	if (flag_verif) then
602	write(79,*) 'ucov',ucov
603	write(80,*) 'vcov',vcov
604	write(81,*) 'teta',teta
605	write(82,*) 'ps',ps
606	write(83,*) 'q',q
607	WRITE(85,*) 'q1 = ',q(:,:,1)
608	WRITE(86,*) 'q3 = ',q(:,:,3)
609	endif
610
611	abort_message = 'Simulation finished'
612
613	call abort_gcm(modname,abort_message,0)
614	ENDIF
615	c-----------------------------------------------------------------------
616	c ecriture du fichier histoire moyenne:
617	c -------------------------------------
618
619	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
620	IF(itau.EQ.itaufin) THEN
621	iav=1
622	ELSE
623	iav=0
624	ENDIF
625
626	IF (ok_dynzon) THEN
627	#ifdef CPP_IOIPSL
628	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
629	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
630	#endif
631	END IF
632	IF (ok_dyn_ave) THEN
633	#ifdef CPP_IOIPSL
634	CALL writedynav(itau,vcov,
635	& ucov,teta,pk,phi,q,masse,ps,phis)
636	#endif
637	ENDIF
638
639	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
640
641	c-----------------------------------------------------------------------
642	c ecriture de la bande histoire:
643	c ------------------------------
644
645	IF( MOD(itau,iecri).EQ.0) THEN
646	! Ehouarn: output only during LF or Backward Matsuno
647	if (leapf.or.(.not.leapf.and.(.not.forward))) then
648	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
649	unat=0.
650	do l=1,llm
651	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
652	vnat(:,l)=vcov(:,l)/cv(:)
653	enddo
654	#ifdef CPP_IOIPSL
655	if (ok_dyn_ins) then
656	! write(lunout,*) "leapfrog: call writehist, itau=",itau
657	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
658	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
659	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
660	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
661	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
662	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
663	endif ! of if (ok_dyn_ins)
664	#endif
665	! For some Grads outputs of fields
666	if (output_grads_dyn) then
667	#include "write_grads_dyn.h"
668	endif
669	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
670	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
671
672	IF(itau.EQ.itaufin) THEN
673
674
675	! if (planet_type.eq."earth") then
676	! Write an Earth-format restart file
677	CALL dynredem1("restart.nc",start_time,
678	& vcov,ucov,teta,q,masse,ps)
679	! endif ! of if (planet_type.eq."earth")
680
681	CLOSE(99)
682	!!! Ehouarn: Why not stop here and now?
683	ENDIF ! of IF (itau.EQ.itaufin)
684
685	c-----------------------------------------------------------------------
686	c gestion de l'integration temporelle:
687	c ------------------------------------
688
689	IF( MOD(itau,iperiod).EQ.0 ) THEN
690	GO TO 1
691	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
692
693	IF( forward ) THEN
694	c fin du pas forward et debut du pas backward
695
696	forward = .FALSE.
697	leapf = .FALSE.
698	GO TO 2
699
700	ELSE
701	c fin du pas backward et debut du premier pas leapfrog
702
703	leapf = .TRUE.
704	dt = 2.*dtvr
705	GO TO 2
706	END IF ! of IF (forward)
707	ELSE
708
709	c ...... pas leapfrog .....
710
711	leapf = .TRUE.
712	dt = 2.*dtvr
713	GO TO 2
714	END IF ! of IF (MOD(itau,iperiod).EQ.0)
715	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
716
717	ELSE ! of IF (.not.purmats)
718
719	c ........................................................
720	c .............. schema matsuno ...............
721	c ........................................................
722	IF( forward ) THEN
723
724	itau = itau + 1
725	c iday = day_ini+itau/day_step
726	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
727	c
728	c IF(time.GT.1.) THEN
729	c time = time-1.
730	c iday = iday+1
731	c 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 ! of IF(forward) i.e. backward step
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
749	IF (ok_dynzon) THEN
750	#ifdef CPP_IOIPSL
751	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
752	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
753	#endif
754	ENDIF
755	IF (ok_dyn_ave) THEN
756	#ifdef CPP_IOIPSL
757	CALL writedynav(itau,vcov,
758	& ucov,teta,pk,phi,q,masse,ps,phis)
759	#endif
760	ENDIF
761
762	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
763
764	IF(MOD(itau,iecri ).EQ.0) THEN
765	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
766	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
767	unat=0.
768	do l=1,llm
769	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
770	vnat(:,l)=vcov(:,l)/cv(:)
771	enddo
772	#ifdef CPP_IOIPSL
773	if (ok_dyn_ins) then
774	! write(lunout,*) "leapfrog: call writehist (b)",
775	! & itau,iecri
776	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
777	endif ! of if (ok_dyn_ins)
778	#endif
779	! For some Grads outputs
780	if (output_grads_dyn) then
781	#include "write_grads_dyn.h"
782	endif
783
784	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
785
786	IF(itau.EQ.itaufin) THEN
787	! if (planet_type.eq."earth") then
788	CALL dynredem1("restart.nc",start_time,
789	& vcov,ucov,teta,q,masse,ps)
790	! endif ! of if (planet_type.eq."earth")
791	ENDIF ! of IF(itau.EQ.itaufin)
792
793	forward = .TRUE.
794	GO TO 1
795
796	ENDIF ! of IF (forward)
797
798	END IF ! of IF(.not.purmats)
799
800	STOP
801	END

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

Download in other formats:

Original Format