Context Navigation

leapfrog.F @ 98

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

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

Download in other formats:

Original Format