Context Navigation

leapfrog.F @ 4400

Last change on this file since 4400 was 1793, checked in by Ehouarn Millour, 11 years ago

Improved sponge layer:

Sponge tendencies are now computed analytically, instead than using a Forward Euler approximation.
Sponge tendencies are added within top_bound, and the sponge is applied after physics tendencies have been taken into account.

These changes imply that GCM results (when using sponge layer) will be differentwrt bench test cases using previous revisions.
EM

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

File size: 25.4 KB

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

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

Download in other formats:

Original Format