Context Navigation

leapfrog.F @ 101

Last change on this file since 101 was 101, checked in by slebonnois, 14 years ago
SL: modifications pour arriver a compiler le gcm VENUS ! Ca marche ! A noter: modifs de makelmdz
File size: 27.1 KB

Line
1	!
2	! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
3	!
4	c
5	c
6	SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,
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	real zqmin,zqmax
70	INTEGER nbetatmoy, nbetatdem,nbetat
71
72	c variables dynamiques
73	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
74	REAL teta(ip1jmp1,llm) ! temperature potentielle
75	REAL q(ip1jmp1,llm,nqtot) ! champs advectes
76	REAL ps(ip1jmp1) ! pression au sol
77	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
78	REAL pks(ip1jmp1) ! exner au sol
79	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
80	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
81	REAL masse(ip1jmp1,llm) ! masse d'air
82	REAL phis(ip1jmp1) ! geopotentiel au sol
83	REAL phi(ip1jmp1,llm) ! geopotentiel
84	REAL w(ip1jmp1,llm) ! vitesse verticale
85	! ADAPTATION GCM POUR CP(T)
86	REAL temp(ip1jmp1,llm) ! temperature
87	REAL tsurpk(ip1jmp1,llm) ! cpp*T/pk
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 en */s
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 en */s
102	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
103	REAL dtetadis(ip1jmp1,llm)
104
105	c tendances de la couche superieure */s
106	REAL dvtop(ip1jm,llm),dutop(ip1jmp1,llm)
107	REAL dtetatop(ip1jmp1,llm)
108	REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
109
110	c tendances physiques */s
111	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
112	REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
113
114	real :: duspg(ip1jmp1,llm) ! for bilan_dyn
115
116	c variables pour le fichier histoire
117	REAL dtav ! intervalle de temps elementaire
118
119	REAL tppn(iim),tpps(iim),tpn,tps
120	c
121	INTEGER itau,itaufinp1,iav
122	! INTEGER iday ! jour julien
123	REAL time
124
125	REAL SSUM
126	REAL time_0 , finvmaold(ip1jmp1,llm)
127
128	cym LOGICAL lafin
129	LOGICAL :: lafin=.false.
130	INTEGER ij,iq,l
131	INTEGER ik
132
133	real time_step, t_wrt, t_ops
134
135	! REAL rdayvrai,rdaym_ini
136	! jD_cur: jour julien courant
137	! jH_cur: heure julienne courante
138	REAL :: jD_cur, jH_cur
139	INTEGER :: an, mois, jour
140	REAL :: secondes
141
142	LOGICAL first,callinigrads
143	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
144	save first
145	data first/.true./
146	real dt_cum
147	character*10 infile
148	integer zan, tau0, thoriid
149	integer nid_ctesGCM
150	save nid_ctesGCM
151	real degres
152	real rlong(iip1), rlatg(jjp1)
153	real zx_tmp_2d(iip1,jjp1)
154	integer ndex2d(iip1*jjp1)
155	logical ok_sync
156	parameter (ok_sync = .true.)
157
158	data callinigrads/.true./
159	character*10 string10
160
161	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
162	REAL :: flxw(ip1jmp1,llm) ! flux de masse verticale
163
164	c+jld variables test conservation energie
165	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
166	C Tendance de la temp. potentiel d (theta)/ d t due a la
167	C tansformation d'energie cinetique en energie thermique
168	C cree par la dissipation
169	REAL dtetaecdt(ip1jmp1,llm)
170	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
171	REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
172	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
173	CHARACTER*15 ztit
174	!IM INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
175	!IM SAVE ip_ebil_dyn
176	!IM DATA ip_ebil_dyn/0/
177	c-jld
178
179	integer :: itau_w ! for write_paramLMDZ_dyn.h
180
181	character*80 dynhist_file, dynhistave_file
182	character(len=20) :: modname
183	character*80 abort_message
184
185	logical dissip_conservative
186	save dissip_conservative
187	data dissip_conservative/.true./
188
189	INTEGER testita
190	PARAMETER (testita = 9)
191
192	logical , parameter :: flag_verif = .false.
193
194	! for CP(T)
195	real :: dtec
196	real,external :: cpdet
197	real :: ztetaec(ip1jmp1,llm)
198
199	c dummy: sinon cette routine n'est jamais compilee...
200	if(1.eq.0) then
201	CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
202	endif
203
204	itaufin = nday*day_step
205	if (less1day) then
206	c MODIF VENUS: to run less than one day:
207	itaufin = int(fractday*day_step)
208	endif
209	itaufinp1 = itaufin +1
210	modname="leapfrog"
211
212	c INITIALISATIONS
213	dudis(:,:) =0.
214	dvdis(:,:) =0.
215	dtetadis(:,:)=0.
216	dutop(:,:) =0.
217	dvtop(:,:) =0.
218	dtetatop(:,:)=0.
219	dqtop(:,:,:) =0.
220	dptop(:) =0.
221	dufi(:,:) =0.
222	dvfi(:,:) =0.
223	dtetafi(:,:)=0.
224	dqfi(:,:,:) =0.
225	dpfi(:) =0.
226
227	itau = 0
228	c iday = day_ini+itau/day_step
229	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
230	c IF(time.GT.1.) THEN
231	c time = time-1.
232	c iday = iday+1
233	c ENDIF
234
235
236	c-----------------------------------------------------------------------
237	c On initialise la pression et la fonction d'Exner :
238	c --------------------------------------------------
239
240	dq(:,:,:)=0.
241	CALL pression ( ip1jmp1, ap, bp, ps, p )
242	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
243	c------------------
244	c TEST PK MONOTONE
245	c------------------
246	write(,) "Test PK"
247	do ij=1,ip1jmp1
248	do l=2,llm
249	if(pk(ij,l).gt.pk(ij,l-1)) then
250	c write(,) ij,l,pk(ij,l)
251	abort_message = 'PK non strictement decroissante'
252	call abort_gcm(modname,abort_message,1)
253	c write(,) "ATTENTION, Test PK deconnecté..."
254	endif
255	enddo
256	enddo
257	write(,) "Fin Test PK"
258	c stop
259	c------------------
260
261	c-----------------------------------------------------------------------
262	c Debut de l'integration temporelle:
263	c ----------------------------------
264
265	1 CONTINUE
266
267	jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
268	jH_cur = jH_ref + &
269	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
270
271
272	#ifdef CPP_IOIPSL
273	if (ok_guide) then
274	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
275	endif
276	#endif
277
278
279	c
280	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
281	c CALL test_period ( ucov,vcov,teta,q,p,phis )
282	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
283	c ENDIF
284	c
285
286	! Save fields obtained at previous time step as '...m1'
287	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
288	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
289	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
290	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
291	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
292
293	forward = .TRUE.
294	leapf = .FALSE.
295	dt = dtvr
296
297	c ... P.Le Van .26/04/94 ....
298
299	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
300	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
301
302	2 CONTINUE
303
304	c-----------------------------------------------------------------------
305
306	c date:
307	c -----
308
309
310	c gestion des appels de la physique et des dissipations:
311	c ------------------------------------------------------
312	c
313	c ... P.Le Van ( 6/02/95 ) ....
314
315	apphys = .FALSE.
316	statcl = .FALSE.
317	conser = .FALSE.
318	apdiss = .FALSE.
319
320	IF( purmats ) THEN
321	! Purely Matsuno time stepping
322	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
323	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
324	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
325	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
326	ELSE
327	! Leapfrog/Matsuno time stepping
328	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
329	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
330	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
331	END IF
332
333	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
334	! supress dissipation step
335	if (llm.eq.1) then
336	apdiss=.false.
337	endif
338
339	c-----------------------------------------------------------------------
340	c calcul des tendances dynamiques:
341	c --------------------------------
342
343	! ADAPTATION GCM POUR CP(T)
344	call tpot2t(ijp1llm,teta,temp,pk)
345	tsurpk = cpp*temp/pk
346	CALL geopot ( ip1jmp1, tsurpk , pk , pks, phis , phi )
347
348	time = jD_cur + jH_cur
349	CALL caldyn
350	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,tsurpk,phis ,
351	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
352
353
354	c-----------------------------------------------------------------------
355	c calcul des tendances advection des traceurs (dont l'humidite)
356	c -------------------------------------------------------------
357
358	IF( forward. OR . leapf ) THEN
359
360	CALL caladvtrac(q,pbaru,pbarv,
361	* p, masse, dq, teta,
362	. flxw, pk)
363
364	IF (offline) THEN
365	Cmaf stokage du flux de masse pour traceurs OFF-LINE
366
367	#ifdef CPP_IOIPSL
368	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
369	. dtvr, itau)
370	#endif
371
372
373	ENDIF ! of IF (offline)
374	c
375	ENDIF ! of IF( forward. OR . leapf )
376
377
378	c-----------------------------------------------------------------------
379	c integrations dynamique et traceurs:
380	c ----------------------------------
381
382
383	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
384	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
385	$ finvmaold )
386
387
388	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
389	c
390	c-----------------------------------------------------------------------
391	c calcul des tendances physiques:
392	c -------------------------------
393	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
394	c
395	IF( purmats ) THEN
396	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
397	ELSE
398	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
399	ENDIF
400	c
401	c
402	IF( apphys ) THEN
403	c
404	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
405	c
406
407	CALL pression ( ip1jmp1, ap, bp, ps, p )
408	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
409
410	! rdaym_ini = itau * dtvr / daysec
411	! rdayvrai = rdaym_ini + day_ini
412	jD_cur = jD_ref + day_ini - day_ref
413	$ + int (itau * dtvr / daysec)
414	jH_cur = jH_ref + &
415	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
416	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
417	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
418	! write(lunout,*)'current date = ',an, mois, jour, secondes
419
420	c rajout debug
421	c lafin = .true.
422
423
424	c Interface avec les routines de phylmd (phymars ... )
425	c -----------------------------------------------------
426
427	c+jld
428
429	c Diagnostique de conservation de l'énergie : initialisation
430	IF (ip_ebil_dyn.ge.1 ) THEN
431	ztit='bil dyn'
432	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
433	IF (planet_type.eq."earth") THEN
434	CALL diagedyn(ztit,2,1,1,dtphys
435	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
436	ENDIF
437	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
438	c-jld
439	#ifdef CPP_IOIPSL
440	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
441	IF (first) THEN
442	#include "ini_paramLMDZ_dyn.h"
443	ENDIF
444	c
445	#include "write_paramLMDZ_dyn.h"
446	c
447	#endif
448	! #endif of #ifdef CPP_IOIPSL
449
450	CALL calfis( lafin , jD_cur, jH_cur,
451	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
452	$ du,dv,dteta,dq,
453	$ flxw,
454	$ dufi,dvfi,dtetafi,dqfi,dpfi )
455
456	c ajout des tendances physiques:
457	c ------------------------------
458	CALL addfi( dtphys, leapf, forward ,
459	$ ucov, vcov, teta , q ,ps ,
460	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
461
462	c Couche superieure :
463	c -------------------
464	IF (ok_strato) THEN
465	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
466	ENDIF
467	c dqtop=0, dptop=0
468
469	c ajout des tendances physiques:
470	c ------------------------------
471	CALL addfi( dtphys, leapf, forward ,
472	$ ucov, vcov, teta , q ,ps ,
473	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
474	c
475	c Diagnostique de conservation de l'énergie : difference
476	IF (ip_ebil_dyn.ge.1 ) THEN
477	ztit='bil phys'
478	IF (planet_type.eq."earth") THEN
479	CALL diagedyn(ztit,2,1,1,dtphys
480	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
481	ENDIF
482	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
483
484	ENDIF ! of IF( apphys )
485
486	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
487	! Academic case : Simple friction and Newtonan relaxation
488	! -------------------------------------------------------
489	DO l=1,llm
490	DO ij=1,ip1jmp1
491	teta(ij,l)=teta(ij,l)-dtvr*
492	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
493	ENDDO
494	ENDDO ! of DO l=1,llm
495	call friction(ucov,vcov,dtvr)
496
497	if (planet_type.eq."giant") then
498	! Intrinsic heat flux
499	! Aymeric -- for giant planets
500	if (ihf .gt. 1.e-6) then
501	!print , '* INTRINSIC HEAT FLUX **', ihf
502	DO ij=1,ip1jmp1
503	teta(ij,1) = teta(ij,1)
504	& + dtvr * aire(ij) * ihf / cpp / masse(ij,1)
505	ENDDO
506	!print , '*** d teta '
507	!print , dtvr aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
508	endif
509	endif
510
511
512	! Sponge layer (if any)
513	IF (ok_strato) THEN
514	dutop(:,:)=0.
515	dvtop(:,:)=0.
516	dtetatop(:,:)=0.
517	dqtop(:,:,:)=0.
518	dptop(:)=0.
519	CALL top_bound(vcov,ucov,teta,masse,dutop,dvtop,dtetatop)
520	CALL addfi( dtvr, leapf, forward ,
521	$ ucov, vcov, teta , q ,ps ,
522	$ dutop, dvtop, dtetatop , dqtop ,dptop )
523	ENDIF ! of IF (ok_strato)
524	ENDIF ! of IF (iflag_phys.EQ.2)
525
526
527	c-jld
528
529	CALL pression ( ip1jmp1, ap, bp, ps, p )
530	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
531
532
533	c-----------------------------------------------------------------------
534	c dissipation horizontale et verticale des petites echelles:
535	c ----------------------------------------------------------
536
537	IF(apdiss) THEN
538
539
540	c calcul de l'energie cinetique avant dissipation
541	call covcont(llm,ucov,vcov,ucont,vcont)
542	call enercin(vcov,ucov,vcont,ucont,ecin0)
543
544	c dissipation
545	! ADAPTATION GCM POUR CP(T)
546	call tpot2t(ijp1llm,teta,temp,pk)
547
548	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
549	ucov=ucov+dudis
550	vcov=vcov+dvdis
551	dudis=dudis/dtdiss ! passage en (m/s)/s
552	dvdis=dvdis/dtdiss ! passage en (m/s)/s
553
554	c------------------------------------------------------------------------
555	if (dissip_conservative) then
556	C On rajoute la tendance due a la transform. Ec -> E therm. cree
557	C lors de la dissipation
558	call covcont(llm,ucov,vcov,ucont,vcont)
559	call enercin(vcov,ucov,vcont,ucont,ecin)
560	! ADAPTATION GCM POUR CP(T)
561	do ij=1,ip1jmp1
562	do l=1,llm
563	dtec = (ecin0(ij,l)-ecin(ij,l))/cpdet(temp(ij,l))
564	temp(ij,l) = temp(ij,l) + dtec
565	enddo
566	enddo
567	call t2tpot(ijp1llm,temp,ztetaec,pk)
568	dtetaecdt=ztetaec-teta
569	dtetadis=dtetadis+dtetaecdt
570	endif
571	teta=teta+dtetadis
572	dtetadis=dtetadis/dtdiss ! passage en K/s
573	c------------------------------------------------------------------------
574
575
576	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
577	c ... Calcul de la valeur moyenne, unique de h aux poles .....
578	c
579
580	DO l = 1, llm
581	DO ij = 1,iim
582	tppn(ij) = aire( ij ) * teta( ij ,l)
583	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
584	ENDDO
585	tpn = SSUM(iim,tppn,1)/apoln
586	tps = SSUM(iim,tpps,1)/apols
587
588	DO ij = 1, iip1
589	teta( ij ,l) = tpn
590	teta(ij+ip1jm,l) = tps
591	ENDDO
592	ENDDO
593
594	DO ij = 1,iim
595	tppn(ij) = aire( ij ) * ps ( ij )
596	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
597	ENDDO
598	tpn = SSUM(iim,tppn,1)/apoln
599	tps = SSUM(iim,tpps,1)/apols
600
601	DO ij = 1, iip1
602	ps( ij ) = tpn
603	ps(ij+ip1jm) = tps
604	ENDDO
605
606
607	END IF ! of IF(apdiss)
608
609	c ajout debug
610	c IF( lafin ) then
611	c abort_message = 'Simulation finished'
612	c call abort_gcm(modname,abort_message,0)
613	c ENDIF
614
615	c ********************************************************************
616	c ********************************************************************
617	c .... fin de l'integration dynamique et physique pour le pas itau ..
618	c ********************************************************************
619	c ********************************************************************
620
621	c preparation du pas d'integration suivant ......
622
623	IF ( .NOT.purmats ) THEN
624	c ........................................................
625	c .............. schema matsuno + leapfrog ..............
626	c ........................................................
627
628	IF(forward. OR. leapf) THEN
629	itau= itau + 1
630	c iday= day_ini+itau/day_step
631	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
632	c IF(time.GT.1.) THEN
633	c time = time-1.
634	c iday = iday+1
635	c ENDIF
636	ENDIF
637
638
639	IF( itau. EQ. itaufinp1 ) then
640	if (flag_verif) then
641	write(79,*) 'ucov',ucov
642	write(80,*) 'vcov',vcov
643	write(81,*) 'teta',teta
644	write(82,*) 'ps',ps
645	write(83,*) 'q',q
646	WRITE(85,*) 'q1 = ',q(:,:,1)
647	WRITE(86,*) 'q3 = ',q(:,:,3)
648	endif
649
650	abort_message = 'Simulation finished'
651
652	call abort_gcm(modname,abort_message,0)
653	ENDIF
654	c-----------------------------------------------------------------------
655	c ecriture du fichier histoire moyenne:
656	c -------------------------------------
657
658	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
659	IF(itau.EQ.itaufin) THEN
660	iav=1
661	ELSE
662	iav=0
663	ENDIF
664
665	IF (ok_dynzon) THEN
666	#ifdef CPP_IOIPSL
667	c les traceurs ne sont pas sortis, trop lourd.
668	c Peut changer eventuellement si besoin.
669	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
670	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
671	& du,dudis,duspg,dufi)
672	#endif
673	END IF
674	IF (ok_dyn_ave) THEN
675	#ifdef CPP_IOIPSL
676	CALL writedynav(itau,vcov,
677	& ucov,teta,pk,phi,q,masse,ps,phis)
678	#endif
679	ENDIF
680
681	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
682
683	c-----------------------------------------------------------------------
684	c ecriture de la bande histoire:
685	c ------------------------------
686
687	IF( MOD(itau,iecri).EQ.0) THEN
688	! Ehouarn: output only during LF or Backward Matsuno
689	if (leapf.or.(.not.leapf.and.(.not.forward))) then
690	nbetat = nbetatdem
691	! ADAPTATION GCM POUR CP(T)
692	call tpot2t(ijp1llm,teta,temp,pk)
693	tsurpk = cpp*temp/pk
694	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
695	unat=0.
696	do l=1,llm
697	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
698	vnat(:,l)=vcov(:,l)/cv(:)
699	enddo
700	#ifdef CPP_IOIPSL
701	if (ok_dyn_ins) then
702	! write(lunout,*) "leapfrog: call writehist, itau=",itau
703	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
704	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
705	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
706	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
707	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
708	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
709	endif ! of if (ok_dyn_ins)
710	#endif
711	! For some Grads outputs of fields
712	if (output_grads_dyn) then
713	#include "write_grads_dyn.h"
714	endif
715	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
716	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
717
718	IF(itau.EQ.itaufin) THEN
719
720
721	if (planet_type.eq."mars") then
722	! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
723	abort_message = 'dynredem1_mars A FAIRE'
724	call abort_gcm(modname,abort_message,0)
725	else
726	CALL dynredem1("restart.nc",0.0,
727	& vcov,ucov,teta,q,masse,ps)
728	endif ! of if (planet_type.eq."mars")
729
730	CLOSE(99)
731	ENDIF ! of IF (itau.EQ.itaufin)
732
733	c-----------------------------------------------------------------------
734	c gestion de l'integration temporelle:
735	c ------------------------------------
736
737	IF( MOD(itau,iperiod).EQ.0 ) THEN
738	GO TO 1
739	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
740
741	IF( forward ) THEN
742	c fin du pas forward et debut du pas backward
743
744	forward = .FALSE.
745	leapf = .FALSE.
746	GO TO 2
747
748	ELSE
749	c fin du pas backward et debut du premier pas leapfrog
750
751	leapf = .TRUE.
752	dt = 2.*dtvr
753	GO TO 2
754	END IF ! of IF (forward)
755	ELSE
756
757	c ...... pas leapfrog .....
758
759	leapf = .TRUE.
760	dt = 2.*dtvr
761	GO TO 2
762	END IF ! of IF (MOD(itau,iperiod).EQ.0)
763	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
764
765	ELSE ! of IF (.not.purmats)
766
767	c ........................................................
768	c .............. schema matsuno ...............
769	c ........................................................
770	IF( forward ) THEN
771
772	itau = itau + 1
773	c iday = day_ini+itau/day_step
774	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
775	c
776	c IF(time.GT.1.) THEN
777	c time = time-1.
778	c iday = iday+1
779	c ENDIF
780
781	forward = .FALSE.
782	IF( itau. EQ. itaufinp1 ) then
783	abort_message = 'Simulation finished'
784	call abort_gcm(modname,abort_message,0)
785	ENDIF
786	GO TO 2
787
788	ELSE ! of IF(forward) i.e. backward step
789
790	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
791	IF(itau.EQ.itaufin) THEN
792	iav=1
793	ELSE
794	iav=0
795	ENDIF
796
797	IF (ok_dynzon) THEN
798	#ifdef CPP_IOIPSL
799	c les traceurs ne sont pas sortis, trop lourd.
800	c Peut changer eventuellement si besoin.
801	CALL bilan_dyn(dtvriperiod,dtvrday_step*periodav,
802	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
803	& du,dudis,duspg,dufi)
804	#endif
805	ENDIF
806	IF (ok_dyn_ave) THEN
807	#ifdef CPP_IOIPSL
808	CALL writedynav(itau,vcov,
809	& ucov,teta,pk,phi,q,masse,ps,phis)
810	#endif
811	ENDIF
812
813	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
814
815	IF(MOD(itau,iecri ).EQ.0) THEN
816	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
817	nbetat = nbetatdem
818	! ADAPTATION GCM POUR CP(T)
819	call tpot2t(ijp1llm,teta,temp,pk)
820	tsurpk = cpp*temp/pk
821	CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
822	unat=0.
823	do l=1,llm
824	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
825	vnat(:,l)=vcov(:,l)/cv(:)
826	enddo
827	#ifdef CPP_IOIPSL
828	if (ok_dyn_ins) then
829	! write(lunout,*) "leapfrog: call writehist (b)",
830	! & itau,iecri
831	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
832	endif ! of if (ok_dyn_ins)
833	#endif
834	! For some Grads outputs
835	if (output_grads_dyn) then
836	#include "write_grads_dyn.h"
837	endif
838
839	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
840
841	IF(itau.EQ.itaufin) THEN
842	if (planet_type.eq."mars") then
843	! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
844	abort_message = 'dynredem1_mars A FAIRE'
845	call abort_gcm(modname,abort_message,0)
846	else
847	CALL dynredem1("restart.nc",0.0,
848	& vcov,ucov,teta,q,masse,ps)
849	endif ! of if (planet_type.eq."mars")
850	ENDIF ! of IF(itau.EQ.itaufin)
851
852	forward = .TRUE.
853	GO TO 1
854
855	ENDIF ! of IF (forward)
856
857	END IF ! of IF(.not.purmats)
858
859	first=.false.
860
861	STOP
862	END

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

Download in other formats:

Original Format