Context Navigation

leapfrog.F @ 4682

Last change on this file since 4682 was 1357, checked in by Ehouarn Millour, 15 years ago

Some cleanup and fixing the possibility to output fields in the dynamics, on the dynamical grids.

CLEANUPS:

arch-PW6_VARGAS.fcm : add potentially benefic compiling options
removed obsolete "control.h" in dyn3d/dyn3dpar (module control_mod.F90 is used instead)

OUTPUTS in the dynamics (3 sets of files, one for each grid: scalar, u, v):

removed "com_io_dyn.h" common; use module "com_io_dyn_mod.F90" instead
updated "initdynav.F","inithist.F","writehist.F" and "writedynav.F" in bibio: which field will be written is hard coded there.
flags "ok_dyn_ins" and "ok_dyn_ave" (loaded via conf_gcm.F) trigger output of fields in the dynamics: if ok_dyn_ins is true, then files "dyn_hist.nc", "dyn_histu.nc" and "dyn_histv.nc" are written (the frequency of the outputs is given by 'iecri' in run.def; values are written every 'iecri' dynamical step). if ok_dyn_ave is true then files "dyn_hist_ave.nc", "dyn_histu_ave.nc" and "dyn_histv_ave.nc" are written (the rate at which averages and made/written, in days, is given by 'periodav' in run.def).

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

File size: 23.1 KB

Line
1	!
2	! $Id: leapfrog.F 1357 2010-04-14 14:03:19Z fhourdin $
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	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 = REAL(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	! Purely Matsuno time stepping
279	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
280	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
281	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
282	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
283	ELSE
284	! Leapfrog/Matsuno time stepping
285	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
286	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
287	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
288	END IF
289
290	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
291	! supress dissipation step
292	if (llm.eq.1) then
293	apdiss=.false.
294	endif
295
296	c-----------------------------------------------------------------------
297	c calcul des tendances dynamiques:
298	c --------------------------------
299
300	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
301
302	time = jD_cur + jH_cur
303	CALL caldyn
304	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
305	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
306
307
308	c-----------------------------------------------------------------------
309	c calcul des tendances advection des traceurs (dont l'humidite)
310	c -------------------------------------------------------------
311
312	IF( forward. OR . leapf ) THEN
313
314	CALL caladvtrac(q,pbaru,pbarv,
315	* p, masse, dq, teta,
316	. flxw, pk)
317
318	IF (offline) THEN
319	Cmaf stokage du flux de masse pour traceurs OFF-LINE
320
321	#ifdef CPP_IOIPSL
322	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
323	. dtvr, itau)
324	#endif
325
326
327	ENDIF ! of IF (offline)
328	c
329	ENDIF ! of IF( forward. OR . leapf )
330
331
332	c-----------------------------------------------------------------------
333	c integrations dynamique et traceurs:
334	c ----------------------------------
335
336
337	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
338	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
339	$ finvmaold )
340
341
342	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
343	c
344	c-----------------------------------------------------------------------
345	c calcul des tendances physiques:
346	c -------------------------------
347	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
348	c
349	IF( purmats ) THEN
350	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
351	ELSE
352	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
353	ENDIF
354	c
355	c
356	IF( apphys ) THEN
357	c
358	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
359	c
360
361	CALL pression ( ip1jmp1, ap, bp, ps, p )
362	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
363
364	! rdaym_ini = itau * dtvr / daysec
365	! rdayvrai = rdaym_ini + day_ini
366	jD_cur = jD_ref + day_ini - day_ref
367	$ + int (itau * dtvr / daysec)
368	jH_cur = jH_ref + &
369	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
370	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
371	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
372	! write(lunout,*)'current date = ',an, mois, jour, secondes
373
374	c rajout debug
375	c lafin = .true.
376
377
378	c Inbterface avec les routines de phylmd (phymars ... )
379	c -----------------------------------------------------
380
381	c+jld
382
383	c Diagnostique de conservation de l'énergie : initialisation
384	IF (ip_ebil_dyn.ge.1 ) THEN
385	ztit='bil dyn'
386	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
387	IF (planet_type.eq."earth") THEN
388	CALL diagedyn(ztit,2,1,1,dtphys
389	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
390	ENDIF
391	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
392	c-jld
393	#ifdef CPP_IOIPSL
394	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
395	IF (first) THEN
396	first=.false.
397	#include "ini_paramLMDZ_dyn.h"
398	ENDIF
399	c
400	#include "write_paramLMDZ_dyn.h"
401	c
402	#endif
403	! #endif of #ifdef CPP_IOIPSL
404	CALL calfis( lafin , jD_cur, jH_cur,
405	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
406	$ du,dv,dteta,dq,
407	$ flxw,
408	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
409
410	IF (ok_strato) THEN
411	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
412	ENDIF
413
414	c ajout des tendances physiques:
415	c ------------------------------
416	CALL addfi( dtphys, leapf, forward ,
417	$ ucov, vcov, teta , q ,ps ,
418	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
419	c
420	c Diagnostique de conservation de l'énergie : difference
421	IF (ip_ebil_dyn.ge.1 ) THEN
422	ztit='bil phys'
423	IF (planet_type.eq."earth") THEN
424	CALL diagedyn(ztit,2,1,1,dtphys
425	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
426	ENDIF
427	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
428
429	ENDIF ! of IF( apphys )
430
431	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
432	c Calcul academique de la physique = Rappel Newtonien + friction
433	c --------------------------------------------------------------
434	teta(:,:)=teta(:,:)
435	s -iphysiqdtvr(teta(:,:)-tetarappel(:,:))/taurappel
436	call friction(ucov,vcov,iphysiq*dtvr)
437	ENDIF
438
439
440	c-jld
441
442	CALL pression ( ip1jmp1, ap, bp, ps, p )
443	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
444
445
446	c-----------------------------------------------------------------------
447	c dissipation horizontale et verticale des petites echelles:
448	c ----------------------------------------------------------
449
450	IF(apdiss) THEN
451
452
453	c calcul de l'energie cinetique avant dissipation
454	call covcont(llm,ucov,vcov,ucont,vcont)
455	call enercin(vcov,ucov,vcont,ucont,ecin0)
456
457	c dissipation
458	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
459	ucov=ucov+dudis
460	vcov=vcov+dvdis
461	c teta=teta+dtetadis
462
463
464	c------------------------------------------------------------------------
465	if (dissip_conservative) then
466	C On rajoute la tendance due a la transform. Ec -> E therm. cree
467	C lors de la dissipation
468	call covcont(llm,ucov,vcov,ucont,vcont)
469	call enercin(vcov,ucov,vcont,ucont,ecin)
470	dtetaecdt= (ecin0-ecin)/ pk
471	c teta=teta+dtetaecdt
472	dtetadis=dtetadis+dtetaecdt
473	endif
474	teta=teta+dtetadis
475	c------------------------------------------------------------------------
476
477
478	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
479	c ... Calcul de la valeur moyenne, unique de h aux poles .....
480	c
481
482	DO l = 1, llm
483	DO ij = 1,iim
484	tppn(ij) = aire( ij ) * teta( ij ,l)
485	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
486	ENDDO
487	tpn = SSUM(iim,tppn,1)/apoln
488	tps = SSUM(iim,tpps,1)/apols
489
490	DO ij = 1, iip1
491	teta( ij ,l) = tpn
492	teta(ij+ip1jm,l) = tps
493	ENDDO
494	ENDDO
495
496	DO ij = 1,iim
497	tppn(ij) = aire( ij ) * ps ( ij )
498	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
499	ENDDO
500	tpn = SSUM(iim,tppn,1)/apoln
501	tps = SSUM(iim,tpps,1)/apols
502
503	DO ij = 1, iip1
504	ps( ij ) = tpn
505	ps(ij+ip1jm) = tps
506	ENDDO
507
508
509	END IF ! of IF(apdiss)
510
511	c ajout debug
512	c IF( lafin ) then
513	c abort_message = 'Simulation finished'
514	c call abort_gcm(modname,abort_message,0)
515	c ENDIF
516
517	c ********************************************************************
518	c ********************************************************************
519	c .... fin de l'integration dynamique et physique pour le pas itau ..
520	c ********************************************************************
521	c ********************************************************************
522
523	c preparation du pas d'integration suivant ......
524
525	IF ( .NOT.purmats ) THEN
526	c ........................................................
527	c .............. schema matsuno + leapfrog ..............
528	c ........................................................
529
530	IF(forward. OR. leapf) THEN
531	itau= itau + 1
532	c iday= day_ini+itau/day_step
533	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
534	c IF(time.GT.1.) THEN
535	c time = time-1.
536	c iday = iday+1
537	c ENDIF
538	ENDIF
539
540
541	IF( itau. EQ. itaufinp1 ) then
542	if (flag_verif) then
543	write(79,*) 'ucov',ucov
544	write(80,*) 'vcov',vcov
545	write(81,*) 'teta',teta
546	write(82,*) 'ps',ps
547	write(83,*) 'q',q
548	WRITE(85,*) 'q1 = ',q(:,:,1)
549	WRITE(86,*) 'q3 = ',q(:,:,3)
550	endif
551
552	abort_message = 'Simulation finished'
553
554	call abort_gcm(modname,abort_message,0)
555	ENDIF
556	c-----------------------------------------------------------------------
557	c ecriture du fichier histoire moyenne:
558	c -------------------------------------
559
560	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
561	IF(itau.EQ.itaufin) THEN
562	iav=1
563	ELSE
564	iav=0
565	ENDIF
566
567	IF (ok_dynzon) THEN
568	#ifdef CPP_IOIPSL
569	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
570	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
571	#endif
572	END IF
573	IF (ok_dyn_ave) THEN
574	#ifdef CPP_IOIPSL
575	CALL writedynav(itau,vcov,
576	& ucov,teta,pk,phi,q,masse,ps,phis)
577	#endif
578	ENDIF
579
580	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
581
582	c-----------------------------------------------------------------------
583	c ecriture de la bande histoire:
584	c ------------------------------
585
586	IF( MOD(itau,iecri).EQ.0) THEN
587	! Ehouarn: output only during LF or Backward Matsuno
588	if (leapf.or.(.not.leapf.and.(.not.forward))) then
589	nbetat = nbetatdem
590	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
591	unat=0.
592	do l=1,llm
593	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
594	vnat(:,l)=vcov(:,l)/cv(:)
595	enddo
596	#ifdef CPP_IOIPSL
597	if (ok_dyn_ins) then
598	! write(lunout,*) "leapfrog: call writehist, itau=",itau
599	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
600	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
601	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
602	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
603	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
604	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
605	endif ! of if (ok_dyn_ins)
606	#endif
607	! For some Grads outputs of fields
608	if (output_grads_dyn) then
609	#include "write_grads_dyn.h"
610	endif
611	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
612	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
613
614	IF(itau.EQ.itaufin) THEN
615
616
617	if (planet_type.eq."earth") then
618	! Write an Earth-format restart file
619	CALL dynredem1("restart.nc",0.0,
620	& vcov,ucov,teta,q,masse,ps)
621	endif ! of if (planet_type.eq."earth")
622
623	CLOSE(99)
624	ENDIF ! of IF (itau.EQ.itaufin)
625
626	c-----------------------------------------------------------------------
627	c gestion de l'integration temporelle:
628	c ------------------------------------
629
630	IF( MOD(itau,iperiod).EQ.0 ) THEN
631	GO TO 1
632	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
633
634	IF( forward ) THEN
635	c fin du pas forward et debut du pas backward
636
637	forward = .FALSE.
638	leapf = .FALSE.
639	GO TO 2
640
641	ELSE
642	c fin du pas backward et debut du premier pas leapfrog
643
644	leapf = .TRUE.
645	dt = 2.*dtvr
646	GO TO 2
647	END IF ! of IF (forward)
648	ELSE
649
650	c ...... pas leapfrog .....
651
652	leapf = .TRUE.
653	dt = 2.*dtvr
654	GO TO 2
655	END IF ! of IF (MOD(itau,iperiod).EQ.0)
656	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
657
658	ELSE ! of IF (.not.purmats)
659
660	c ........................................................
661	c .............. schema matsuno ...............
662	c ........................................................
663	IF( forward ) THEN
664
665	itau = itau + 1
666	c iday = day_ini+itau/day_step
667	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
668	c
669	c IF(time.GT.1.) THEN
670	c time = time-1.
671	c iday = iday+1
672	c ENDIF
673
674	forward = .FALSE.
675	IF( itau. EQ. itaufinp1 ) then
676	abort_message = 'Simulation finished'
677	call abort_gcm(modname,abort_message,0)
678	ENDIF
679	GO TO 2
680
681	ELSE ! of IF(forward) i.e. backward step
682
683	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
684	IF(itau.EQ.itaufin) THEN
685	iav=1
686	ELSE
687	iav=0
688	ENDIF
689
690	IF (ok_dynzon) THEN
691	#ifdef CPP_IOIPSL
692	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
693	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
694	#endif
695	ENDIF
696	IF (ok_dyn_ave) THEN
697	#ifdef CPP_IOIPSL
698	CALL writedynav(itau,vcov,
699	& ucov,teta,pk,phi,q,masse,ps,phis)
700	#endif
701	ENDIF
702
703	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
704
705	IF(MOD(itau,iecri ).EQ.0) THEN
706	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
707	nbetat = nbetatdem
708	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
709	unat=0.
710	do l=1,llm
711	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
712	vnat(:,l)=vcov(:,l)/cv(:)
713	enddo
714	#ifdef CPP_IOIPSL
715	if (ok_dyn_ins) then
716	! write(lunout,*) "leapfrog: call writehist (b)",
717	! & itau,iecri
718	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
719	endif ! of if (ok_dyn_ins)
720	#endif
721	! For some Grads outputs
722	if (output_grads_dyn) then
723	#include "write_grads_dyn.h"
724	endif
725
726	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
727
728	IF(itau.EQ.itaufin) THEN
729	if (planet_type.eq."earth") then
730	CALL dynredem1("restart.nc",0.0,
731	& vcov,ucov,teta,q,masse,ps)
732	endif ! of if (planet_type.eq."earth")
733	ENDIF ! of IF(itau.EQ.itaufin)
734
735	forward = .TRUE.
736	GO TO 1
737
738	ENDIF ! of IF (forward)
739
740	END IF ! of IF(.not.purmats)
741
742	STOP
743	END

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

Download in other formats:

Original Format