Context Navigation

leapfrog.F @ 1616

Last change on this file since 1616 was 1616, checked in by Ehouarn Millour, 12 years ago
Some cleanup around what is done during the integration step of dynamical tendencies and namely removed computation of (unused) finvmaold, thereby saving us the expense of a call to the (costly) filter at every dynamical time step. Checked (on Vargas, in seq, omp, mpi and mixed mode) that this doesn't change the GCM results, as expected. EM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.3 KB

Line
1	!
2	! $Id: leapfrog.F 1616 2012-02-17 11:59:00Z emillour $
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 (disvert_type==1) then
215	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
216	else ! we assume that we are in the disvert_type==2 case
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	& int (itau * dtvr / daysec)
228	jH_cur = jH_ref + start_time + &
229	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
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 (disvert_type==1) then
372	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
373	else ! we assume that we are in the disvert_type==2 case
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	& int (itau * dtvr / daysec)
385	jH_cur = jH_ref + start_time + &
386	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
387	jD_cur = jD_cur + int(jH_cur)
388	jH_cur = jH_cur - int(jH_cur)
389	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
390	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
391	! write(lunout,*)'current date = ',an, mois, jour, secondes
392
393	c rajout debug
394	c lafin = .true.
395
396
397	c Inbterface avec les routines de phylmd (phymars ... )
398	c -----------------------------------------------------
399
400	c+jld
401
402	c Diagnostique de conservation de l'énergie : initialisation
403	IF (ip_ebil_dyn.ge.1 ) THEN
404	ztit='bil dyn'
405	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
406	IF (planet_type.eq."earth") THEN
407	#ifdef CPP_EARTH
408	CALL diagedyn(ztit,2,1,1,dtphys
409	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
410	#endif
411	ENDIF
412	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
413	c-jld
414	#ifdef CPP_IOIPSL
415	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
416	IF (first) THEN
417	first=.false.
418	#include "ini_paramLMDZ_dyn.h"
419	ENDIF
420	c
421	#include "write_paramLMDZ_dyn.h"
422	c
423	#endif
424	! #endif of #ifdef CPP_IOIPSL
425	CALL calfis( lafin , jD_cur, jH_cur,
426	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
427	$ du,dv,dteta,dq,
428	$ flxw,
429	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
430
431	IF (ok_strato) THEN
432	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
433	ENDIF
434
435	c ajout des tendances physiques:
436	c ------------------------------
437	CALL addfi( dtphys, leapf, forward ,
438	$ ucov, vcov, teta , q ,ps ,
439	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
440	c
441	c Diagnostique de conservation de l'énergie : difference
442	IF (ip_ebil_dyn.ge.1 ) THEN
443	ztit='bil phys'
444	IF (planet_type.eq."earth") THEN
445	CALL diagedyn(ztit,2,1,1,dtphys
446	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
447	ENDIF
448	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
449
450	ENDIF ! of IF( apphys )
451
452	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
453	! Academic case : Simple friction and Newtonan relaxation
454	! -------------------------------------------------------
455	DO l=1,llm
456	DO ij=1,ip1jmp1
457	teta(ij,l)=teta(ij,l)-dtvr*
458	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
459	ENDDO
460	ENDDO ! of DO l=1,llm
461
462	if (planet_type.eq."giant") then
463	! add an intrinsic heat flux at the base of the atmosphere
464	teta(:,1)=teta(:,1)+dtvraire(:)ihf/cpp/masse(:,1)
465	endif
466
467	call friction(ucov,vcov,dtvr)
468
469	! Sponge layer (if any)
470	IF (ok_strato) THEN
471	dufi(:,:)=0.
472	dvfi(:,:)=0.
473	dtetafi(:,:)=0.
474	dqfi(:,:,:)=0.
475	dpfi(:)=0.
476	CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
477	CALL addfi( dtvr, leapf, forward ,
478	$ ucov, vcov, teta , q ,ps ,
479	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
480	ENDIF ! of IF (ok_strato)
481	ENDIF ! of IF (iflag_phys.EQ.2)
482
483
484	c-jld
485
486	CALL pression ( ip1jmp1, ap, bp, ps, p )
487	if (disvert_type==1) then
488	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
489	else ! we assume that we are in the disvert_type==2 case
490	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
491	endif
492
493
494	c-----------------------------------------------------------------------
495	c dissipation horizontale et verticale des petites echelles:
496	c ----------------------------------------------------------
497
498	IF(apdiss) THEN
499
500
501	c calcul de l'energie cinetique avant dissipation
502	call covcont(llm,ucov,vcov,ucont,vcont)
503	call enercin(vcov,ucov,vcont,ucont,ecin0)
504
505	c dissipation
506	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
507	ucov=ucov+dudis
508	vcov=vcov+dvdis
509	c teta=teta+dtetadis
510
511
512	c------------------------------------------------------------------------
513	if (dissip_conservative) then
514	C On rajoute la tendance due a la transform. Ec -> E therm. cree
515	C lors de la dissipation
516	call covcont(llm,ucov,vcov,ucont,vcont)
517	call enercin(vcov,ucov,vcont,ucont,ecin)
518	dtetaecdt= (ecin0-ecin)/ pk
519	c teta=teta+dtetaecdt
520	dtetadis=dtetadis+dtetaecdt
521	endif
522	teta=teta+dtetadis
523	c------------------------------------------------------------------------
524
525
526	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
527	c ... Calcul de la valeur moyenne, unique de h aux poles .....
528	c
529
530	DO l = 1, llm
531	DO ij = 1,iim
532	tppn(ij) = aire( ij ) * teta( ij ,l)
533	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
534	ENDDO
535	tpn = SSUM(iim,tppn,1)/apoln
536	tps = SSUM(iim,tpps,1)/apols
537
538	DO ij = 1, iip1
539	teta( ij ,l) = tpn
540	teta(ij+ip1jm,l) = tps
541	ENDDO
542	ENDDO
543
544	if (1 == 0) then
545	!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
546	!!! 2) should probably not be here anyway
547	!!! but are kept for those who would want to revert to previous behaviour
548	DO ij = 1,iim
549	tppn(ij) = aire( ij ) * ps ( ij )
550	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
551	ENDDO
552	tpn = SSUM(iim,tppn,1)/apoln
553	tps = SSUM(iim,tpps,1)/apols
554
555	DO ij = 1, iip1
556	ps( ij ) = tpn
557	ps(ij+ip1jm) = tps
558	ENDDO
559	endif ! of if (1 == 0)
560
561	END IF ! of IF(apdiss)
562
563	c ajout debug
564	c IF( lafin ) then
565	c abort_message = 'Simulation finished'
566	c call abort_gcm(modname,abort_message,0)
567	c ENDIF
568
569	c ********************************************************************
570	c ********************************************************************
571	c .... fin de l'integration dynamique et physique pour le pas itau ..
572	c ********************************************************************
573	c ********************************************************************
574
575	c preparation du pas d'integration suivant ......
576
577	IF ( .NOT.purmats ) THEN
578	c ........................................................
579	c .............. schema matsuno + leapfrog ..............
580	c ........................................................
581
582	IF(forward. OR. leapf) THEN
583	itau= itau + 1
584	c iday= day_ini+itau/day_step
585	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
586	c IF(time.GT.1.) THEN
587	c time = time-1.
588	c iday = iday+1
589	c ENDIF
590	ENDIF
591
592
593	IF( itau. EQ. itaufinp1 ) then
594	if (flag_verif) then
595	write(79,*) 'ucov',ucov
596	write(80,*) 'vcov',vcov
597	write(81,*) 'teta',teta
598	write(82,*) 'ps',ps
599	write(83,*) 'q',q
600	WRITE(85,*) 'q1 = ',q(:,:,1)
601	WRITE(86,*) 'q3 = ',q(:,:,3)
602	endif
603
604	abort_message = 'Simulation finished'
605
606	call abort_gcm(modname,abort_message,0)
607	ENDIF
608	c-----------------------------------------------------------------------
609	c ecriture du fichier histoire moyenne:
610	c -------------------------------------
611
612	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
613	IF(itau.EQ.itaufin) THEN
614	iav=1
615	ELSE
616	iav=0
617	ENDIF
618
619	IF (ok_dynzon) THEN
620	#ifdef CPP_IOIPSL
621	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
622	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
623	#endif
624	END IF
625	IF (ok_dyn_ave) THEN
626	#ifdef CPP_IOIPSL
627	CALL writedynav(itau,vcov,
628	& ucov,teta,pk,phi,q,masse,ps,phis)
629	#endif
630	ENDIF
631
632	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
633
634	c-----------------------------------------------------------------------
635	c ecriture de la bande histoire:
636	c ------------------------------
637
638	IF( MOD(itau,iecri).EQ.0) THEN
639	! Ehouarn: output only during LF or Backward Matsuno
640	if (leapf.or.(.not.leapf.and.(.not.forward))) then
641	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
642	unat=0.
643	do l=1,llm
644	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
645	vnat(:,l)=vcov(:,l)/cv(:)
646	enddo
647	#ifdef CPP_IOIPSL
648	if (ok_dyn_ins) then
649	! write(lunout,*) "leapfrog: call writehist, itau=",itau
650	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
651	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
652	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
653	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
654	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
655	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
656	endif ! of if (ok_dyn_ins)
657	#endif
658	! For some Grads outputs of fields
659	if (output_grads_dyn) then
660	#include "write_grads_dyn.h"
661	endif
662	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
663	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
664
665	IF(itau.EQ.itaufin) THEN
666
667
668	! if (planet_type.eq."earth") then
669	! Write an Earth-format restart file
670	CALL dynredem1("restart.nc",start_time,
671	& vcov,ucov,teta,q,masse,ps)
672	! endif ! of if (planet_type.eq."earth")
673
674	CLOSE(99)
675	!!! Ehouarn: Why not stop here and now?
676	ENDIF ! of IF (itau.EQ.itaufin)
677
678	c-----------------------------------------------------------------------
679	c gestion de l'integration temporelle:
680	c ------------------------------------
681
682	IF( MOD(itau,iperiod).EQ.0 ) THEN
683	GO TO 1
684	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
685
686	IF( forward ) THEN
687	c fin du pas forward et debut du pas backward
688
689	forward = .FALSE.
690	leapf = .FALSE.
691	GO TO 2
692
693	ELSE
694	c fin du pas backward et debut du premier pas leapfrog
695
696	leapf = .TRUE.
697	dt = 2.*dtvr
698	GO TO 2
699	END IF ! of IF (forward)
700	ELSE
701
702	c ...... pas leapfrog .....
703
704	leapf = .TRUE.
705	dt = 2.*dtvr
706	GO TO 2
707	END IF ! of IF (MOD(itau,iperiod).EQ.0)
708	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
709
710	ELSE ! of IF (.not.purmats)
711
712	c ........................................................
713	c .............. schema matsuno ...............
714	c ........................................................
715	IF( forward ) THEN
716
717	itau = itau + 1
718	c iday = day_ini+itau/day_step
719	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
720	c
721	c IF(time.GT.1.) THEN
722	c time = time-1.
723	c iday = iday+1
724	c ENDIF
725
726	forward = .FALSE.
727	IF( itau. EQ. itaufinp1 ) then
728	abort_message = 'Simulation finished'
729	call abort_gcm(modname,abort_message,0)
730	ENDIF
731	GO TO 2
732
733	ELSE ! of IF(forward) i.e. backward step
734
735	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
736	IF(itau.EQ.itaufin) THEN
737	iav=1
738	ELSE
739	iav=0
740	ENDIF
741
742	IF (ok_dynzon) THEN
743	#ifdef CPP_IOIPSL
744	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
745	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
746	#endif
747	ENDIF
748	IF (ok_dyn_ave) THEN
749	#ifdef CPP_IOIPSL
750	CALL writedynav(itau,vcov,
751	& ucov,teta,pk,phi,q,masse,ps,phis)
752	#endif
753	ENDIF
754
755	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
756
757	IF(MOD(itau,iecri ).EQ.0) THEN
758	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
759	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
760	unat=0.
761	do l=1,llm
762	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
763	vnat(:,l)=vcov(:,l)/cv(:)
764	enddo
765	#ifdef CPP_IOIPSL
766	if (ok_dyn_ins) then
767	! write(lunout,*) "leapfrog: call writehist (b)",
768	! & itau,iecri
769	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
770	endif ! of if (ok_dyn_ins)
771	#endif
772	! For some Grads outputs
773	if (output_grads_dyn) then
774	#include "write_grads_dyn.h"
775	endif
776
777	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
778
779	IF(itau.EQ.itaufin) THEN
780	! if (planet_type.eq."earth") then
781	CALL dynredem1("restart.nc",start_time,
782	& vcov,ucov,teta,q,masse,ps)
783	! endif ! of if (planet_type.eq."earth")
784	ENDIF ! of IF(itau.EQ.itaufin)
785
786	forward = .TRUE.
787	GO TO 1
788
789	ENDIF ! of IF (forward)
790
791	END IF ! of IF(.not.purmats)
792
793	STOP
794	END

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

Download in other formats:

Original Format