Context Navigation

leapfrog.F @ 1520

Last change on this file since 1520 was 1520, checked in by Ehouarn Millour, 14 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

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

File size: 24.3 KB

Line
1	!
2	! $Id: leapfrog.F 1520 2011-05-23 11:37:09Z emillour $
3	!
4	c
5	c
6	SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
7	& time_0)
8
9
10	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
11	#ifdef CPP_IOIPSL
12	use IOIPSL
13	#endif
14	USE infotrac
15	USE guide_mod, ONLY : guide_main
16	USE write_field
17	USE control_mod
18	IMPLICIT NONE
19
20	c ...... Version du 10/01/98 ..........
21
22	c avec coordonnees verticales hybrides
23	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
24
25	c=======================================================================
26	c
27	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
28	c -------
29	c
30	c Objet:
31	c ------
32	c
33	c GCM LMD nouvelle grille
34	c
35	c=======================================================================
36	c
37	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
38	c et possibilite d'appeler une fonction f(y) a derivee tangente
39	c hyperbolique a la place de la fonction a derivee sinusoidale.
40
41	c ... Possibilite de choisir le shema pour l'advection de
42	c q , en modifiant iadv dans traceur.def (10/02) .
43	c
44	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
45	c Pour Van-Leer iadv=10
46	c
47	c-----------------------------------------------------------------------
48	c Declarations:
49	c -------------
50
51	#include "dimensions.h"
52	#include "paramet.h"
53	#include "comconst.h"
54	#include "comdissnew.h"
55	#include "comvert.h"
56	#include "comgeom.h"
57	#include "logic.h"
58	#include "temps.h"
59	#include "ener.h"
60	#include "description.h"
61	#include "serre.h"
62	!#include "com_io_dyn.h"
63	#include "iniprint.h"
64	#include "academic.h"
65
66	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
67	! #include "clesphys.h"
68
69	INTEGER longcles
70	PARAMETER ( longcles = 20 )
71	REAL clesphy0( longcles )
72
73	real zqmin,zqmax
74	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=*),parameter :: modname="leapfrog"
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
196
197	itau = 0
198	c iday = day_ini+itau/day_step
199	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
200	c IF(time.GT.1.) THEN
201	c time = time-1.
202	c iday = iday+1
203	c ENDIF
204
205
206	c-----------------------------------------------------------------------
207	c On initialise la pression et la fonction d'Exner :
208	c --------------------------------------------------
209
210	dq(:,:,:)=0.
211	CALL pression ( ip1jmp1, ap, bp, ps, p )
212	if (disvert_type==1) then
213	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
214	else ! we assume that we are in the disvert_type==2 case
215	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
216	endif
217
218	c-----------------------------------------------------------------------
219	c Debut de l'integration temporelle:
220	c ----------------------------------
221
222	1 CONTINUE
223
224	jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
225	jH_cur = jH_ref + &
226	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
227
228
229	#ifdef CPP_IOIPSL
230	if (ok_guide) then
231	call guide_main(itau,ucov,vcov,teta,q,masse,ps)
232	endif
233	#endif
234
235
236	c
237	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
238	c CALL test_period ( ucov,vcov,teta,q,p,phis )
239	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
240	c ENDIF
241	c
242
243	! Save fields obtained at previous time step as '...m1'
244	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
245	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
246	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
247	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
248	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
249
250	forward = .TRUE.
251	leapf = .FALSE.
252	dt = dtvr
253
254	c ... P.Le Van .26/04/94 ....
255
256	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
257	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
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,dissip_period ).EQ.0.AND..NOT.forward )
281	s apdiss = .TRUE.
282	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
283	s .and. iflag_phys.EQ.1 ) apphys = .TRUE.
284	ELSE
285	! Leapfrog/Matsuno time stepping
286	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
287	IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
288	s apdiss = .TRUE.
289	IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
290	END IF
291
292	! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
293	! supress dissipation step
294	if (llm.eq.1) then
295	apdiss=.false.
296	endif
297
298	c-----------------------------------------------------------------------
299	c calcul des tendances dynamiques:
300	c --------------------------------
301
302	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
303
304	time = jD_cur + jH_cur
305	CALL caldyn
306	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
307	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
308
309
310	c-----------------------------------------------------------------------
311	c calcul des tendances advection des traceurs (dont l'humidite)
312	c -------------------------------------------------------------
313
314	IF( forward. OR . leapf ) THEN
315
316	CALL caladvtrac(q,pbaru,pbarv,
317	* p, masse, dq, teta,
318	. flxw, pk)
319
320	IF (offline) THEN
321	Cmaf stokage du flux de masse pour traceurs OFF-LINE
322
323	#ifdef CPP_IOIPSL
324	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
325	. dtvr, itau)
326	#endif
327
328
329	ENDIF ! of IF (offline)
330	c
331	ENDIF ! of IF( forward. OR . leapf )
332
333
334	c-----------------------------------------------------------------------
335	c integrations dynamique et traceurs:
336	c ----------------------------------
337
338
339	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
340	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
341	$ finvmaold )
342
343
344	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
345	c
346	c-----------------------------------------------------------------------
347	c calcul des tendances physiques:
348	c -------------------------------
349	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
350	c
351	IF( purmats ) THEN
352	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
353	ELSE
354	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
355	ENDIF
356	c
357	c
358	IF( apphys ) THEN
359	c
360	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
361	c
362
363	CALL pression ( ip1jmp1, ap, bp, ps, p )
364	if (disvert_type==1) then
365	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
366	else ! we assume that we are in the disvert_type==2 case
367	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
368	endif
369
370	! rdaym_ini = itau * dtvr / daysec
371	! rdayvrai = rdaym_ini + day_ini
372	jD_cur = jD_ref + day_ini - day_ref
373	$ + int (itau * dtvr / daysec)
374	jH_cur = jH_ref + &
375	& (itau * dtvr / daysec - int(itau * dtvr / daysec))
376	! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
377	! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
378	! write(lunout,*)'current date = ',an, mois, jour, secondes
379
380	c rajout debug
381	c lafin = .true.
382
383
384	c Inbterface avec les routines de phylmd (phymars ... )
385	c -----------------------------------------------------
386
387	c+jld
388
389	c Diagnostique de conservation de l'énergie : initialisation
390	IF (ip_ebil_dyn.ge.1 ) THEN
391	ztit='bil dyn'
392	! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
393	IF (planet_type.eq."earth") THEN
394	CALL diagedyn(ztit,2,1,1,dtphys
395	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
396	ENDIF
397	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
398	c-jld
399	#ifdef CPP_IOIPSL
400	cIM : pour sortir les param. du modele dans un fis. netcdf 110106
401	IF (first) THEN
402	first=.false.
403	#include "ini_paramLMDZ_dyn.h"
404	ENDIF
405	c
406	#include "write_paramLMDZ_dyn.h"
407	c
408	#endif
409	! #endif of #ifdef CPP_IOIPSL
410	CALL calfis( lafin , jD_cur, jH_cur,
411	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
412	$ du,dv,dteta,dq,
413	$ flxw,
414	$ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
415
416	IF (ok_strato) THEN
417	CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
418	ENDIF
419
420	c ajout des tendances physiques:
421	c ------------------------------
422	CALL addfi( dtphys, leapf, forward ,
423	$ ucov, vcov, teta , q ,ps ,
424	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
425	c
426	c Diagnostique de conservation de l'énergie : difference
427	IF (ip_ebil_dyn.ge.1 ) THEN
428	ztit='bil phys'
429	IF (planet_type.eq."earth") THEN
430	CALL diagedyn(ztit,2,1,1,dtphys
431	& , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
432	ENDIF
433	ENDIF ! of IF (ip_ebil_dyn.ge.1 )
434
435	ENDIF ! of IF( apphys )
436
437	IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
438	! Academic case : Simple friction and Newtonan relaxation
439	! -------------------------------------------------------
440	DO l=1,llm
441	DO ij=1,ip1jmp1
442	teta(ij,l)=teta(ij,l)-dtvr*
443	& (teta(ij,l)-tetarappel(ij,l))(knewt_g+knewt_t(l)clat4(ij))
444	ENDDO
445	ENDDO ! of DO l=1,llm
446
447	if (planet_type.eq."giant") then
448	! add an intrinsic heat flux at the base of the atmosphere
449	teta(:,1)=teta(:,1)+dtvraire(:)ihf/cpp/masse(:,1)
450	endif
451
452	call friction(ucov,vcov,dtvr)
453
454	! Sponge layer (if any)
455	IF (ok_strato) THEN
456	dufi(:,:)=0.
457	dvfi(:,:)=0.
458	dtetafi(:,:)=0.
459	dqfi(:,:,:)=0.
460	dpfi(:)=0.
461	CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
462	CALL addfi( dtvr, leapf, forward ,
463	$ ucov, vcov, teta , q ,ps ,
464	$ dufi, dvfi, dtetafi , dqfi ,dpfi )
465	ENDIF ! of IF (ok_strato)
466	ENDIF ! of IF (iflag_phys.EQ.2)
467
468
469	c-jld
470
471	CALL pression ( ip1jmp1, ap, bp, ps, p )
472	if (disvert_type==1) then
473	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
474	else ! we assume that we are in the disvert_type==2 case
475	CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
476	endif
477
478
479	c-----------------------------------------------------------------------
480	c dissipation horizontale et verticale des petites echelles:
481	c ----------------------------------------------------------
482
483	IF(apdiss) THEN
484
485
486	c calcul de l'energie cinetique avant dissipation
487	call covcont(llm,ucov,vcov,ucont,vcont)
488	call enercin(vcov,ucov,vcont,ucont,ecin0)
489
490	c dissipation
491	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
492	ucov=ucov+dudis
493	vcov=vcov+dvdis
494	c teta=teta+dtetadis
495
496
497	c------------------------------------------------------------------------
498	if (dissip_conservative) then
499	C On rajoute la tendance due a la transform. Ec -> E therm. cree
500	C lors de la dissipation
501	call covcont(llm,ucov,vcov,ucont,vcont)
502	call enercin(vcov,ucov,vcont,ucont,ecin)
503	dtetaecdt= (ecin0-ecin)/ pk
504	c teta=teta+dtetaecdt
505	dtetadis=dtetadis+dtetaecdt
506	endif
507	teta=teta+dtetadis
508	c------------------------------------------------------------------------
509
510
511	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
512	c ... Calcul de la valeur moyenne, unique de h aux poles .....
513	c
514
515	DO l = 1, llm
516	DO ij = 1,iim
517	tppn(ij) = aire( ij ) * teta( ij ,l)
518	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
519	ENDDO
520	tpn = SSUM(iim,tppn,1)/apoln
521	tps = SSUM(iim,tpps,1)/apols
522
523	DO ij = 1, iip1
524	teta( ij ,l) = tpn
525	teta(ij+ip1jm,l) = tps
526	ENDDO
527	ENDDO
528
529	DO ij = 1,iim
530	tppn(ij) = aire( ij ) * ps ( ij )
531	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
532	ENDDO
533	tpn = SSUM(iim,tppn,1)/apoln
534	tps = SSUM(iim,tpps,1)/apols
535
536	DO ij = 1, iip1
537	ps( ij ) = tpn
538	ps(ij+ip1jm) = tps
539	ENDDO
540
541
542	END IF ! of IF(apdiss)
543
544	c ajout debug
545	c IF( lafin ) then
546	c abort_message = 'Simulation finished'
547	c call abort_gcm(modname,abort_message,0)
548	c ENDIF
549
550	c ********************************************************************
551	c ********************************************************************
552	c .... fin de l'integration dynamique et physique pour le pas itau ..
553	c ********************************************************************
554	c ********************************************************************
555
556	c preparation du pas d'integration suivant ......
557
558	IF ( .NOT.purmats ) THEN
559	c ........................................................
560	c .............. schema matsuno + leapfrog ..............
561	c ........................................................
562
563	IF(forward. OR. leapf) THEN
564	itau= itau + 1
565	c iday= day_ini+itau/day_step
566	c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
567	c IF(time.GT.1.) THEN
568	c time = time-1.
569	c iday = iday+1
570	c ENDIF
571	ENDIF
572
573
574	IF( itau. EQ. itaufinp1 ) then
575	if (flag_verif) then
576	write(79,*) 'ucov',ucov
577	write(80,*) 'vcov',vcov
578	write(81,*) 'teta',teta
579	write(82,*) 'ps',ps
580	write(83,*) 'q',q
581	WRITE(85,*) 'q1 = ',q(:,:,1)
582	WRITE(86,*) 'q3 = ',q(:,:,3)
583	endif
584
585	abort_message = 'Simulation finished'
586
587	call abort_gcm(modname,abort_message,0)
588	ENDIF
589	c-----------------------------------------------------------------------
590	c ecriture du fichier histoire moyenne:
591	c -------------------------------------
592
593	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
594	IF(itau.EQ.itaufin) THEN
595	iav=1
596	ELSE
597	iav=0
598	ENDIF
599
600	IF (ok_dynzon) THEN
601	#ifdef CPP_IOIPSL
602	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
603	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
604	#endif
605	END IF
606	IF (ok_dyn_ave) THEN
607	#ifdef CPP_IOIPSL
608	CALL writedynav(itau,vcov,
609	& ucov,teta,pk,phi,q,masse,ps,phis)
610	#endif
611	ENDIF
612
613	ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
614
615	c-----------------------------------------------------------------------
616	c ecriture de la bande histoire:
617	c ------------------------------
618
619	IF( MOD(itau,iecri).EQ.0) THEN
620	! Ehouarn: output only during LF or Backward Matsuno
621	if (leapf.or.(.not.leapf.and.(.not.forward))) then
622	nbetat = nbetatdem
623	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
624	unat=0.
625	do l=1,llm
626	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
627	vnat(:,l)=vcov(:,l)/cv(:)
628	enddo
629	#ifdef CPP_IOIPSL
630	if (ok_dyn_ins) then
631	! write(lunout,*) "leapfrog: call writehist, itau=",itau
632	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
633	! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
634	! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
635	! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
636	! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
637	! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
638	endif ! of if (ok_dyn_ins)
639	#endif
640	! For some Grads outputs of fields
641	if (output_grads_dyn) then
642	#include "write_grads_dyn.h"
643	endif
644	endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
645	ENDIF ! of IF(MOD(itau,iecri).EQ.0)
646
647	IF(itau.EQ.itaufin) THEN
648
649
650	! if (planet_type.eq."earth") then
651	! Write an Earth-format restart file
652	CALL dynredem1("restart.nc",0.0,
653	& vcov,ucov,teta,q,masse,ps)
654	! endif ! of if (planet_type.eq."earth")
655
656	CLOSE(99)
657	ENDIF ! of IF (itau.EQ.itaufin)
658
659	c-----------------------------------------------------------------------
660	c gestion de l'integration temporelle:
661	c ------------------------------------
662
663	IF( MOD(itau,iperiod).EQ.0 ) THEN
664	GO TO 1
665	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
666
667	IF( forward ) THEN
668	c fin du pas forward et debut du pas backward
669
670	forward = .FALSE.
671	leapf = .FALSE.
672	GO TO 2
673
674	ELSE
675	c fin du pas backward et debut du premier pas leapfrog
676
677	leapf = .TRUE.
678	dt = 2.*dtvr
679	GO TO 2
680	END IF ! of IF (forward)
681	ELSE
682
683	c ...... pas leapfrog .....
684
685	leapf = .TRUE.
686	dt = 2.*dtvr
687	GO TO 2
688	END IF ! of IF (MOD(itau,iperiod).EQ.0)
689	! ELSEIF (MOD(itau-1,iperiod).EQ.0)
690
691	ELSE ! of IF (.not.purmats)
692
693	c ........................................................
694	c .............. schema matsuno ...............
695	c ........................................................
696	IF( forward ) THEN
697
698	itau = itau + 1
699	c iday = day_ini+itau/day_step
700	c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
701	c
702	c IF(time.GT.1.) THEN
703	c time = time-1.
704	c iday = iday+1
705	c ENDIF
706
707	forward = .FALSE.
708	IF( itau. EQ. itaufinp1 ) then
709	abort_message = 'Simulation finished'
710	call abort_gcm(modname,abort_message,0)
711	ENDIF
712	GO TO 2
713
714	ELSE ! of IF(forward) i.e. backward step
715
716	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
717	IF(itau.EQ.itaufin) THEN
718	iav=1
719	ELSE
720	iav=0
721	ENDIF
722
723	IF (ok_dynzon) THEN
724	#ifdef CPP_IOIPSL
725	CALL bilan_dyn(2,dtvriperiod,dtvrday_step*periodav,
726	& ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
727	#endif
728	ENDIF
729	IF (ok_dyn_ave) THEN
730	#ifdef CPP_IOIPSL
731	CALL writedynav(itau,vcov,
732	& ucov,teta,pk,phi,q,masse,ps,phis)
733	#endif
734	ENDIF
735
736	ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
737
738	IF(MOD(itau,iecri ).EQ.0) THEN
739	c IF(MOD(itau,iecri*day_step).EQ.0) THEN
740	nbetat = nbetatdem
741	CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
742	unat=0.
743	do l=1,llm
744	unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
745	vnat(:,l)=vcov(:,l)/cv(:)
746	enddo
747	#ifdef CPP_IOIPSL
748	if (ok_dyn_ins) then
749	! write(lunout,*) "leapfrog: call writehist (b)",
750	! & itau,iecri
751	CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
752	endif ! of if (ok_dyn_ins)
753	#endif
754	! For some Grads outputs
755	if (output_grads_dyn) then
756	#include "write_grads_dyn.h"
757	endif
758
759	ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
760
761	IF(itau.EQ.itaufin) THEN
762	! if (planet_type.eq."earth") then
763	CALL dynredem1("restart.nc",0.0,
764	& vcov,ucov,teta,q,masse,ps)
765	! endif ! of if (planet_type.eq."earth")
766	ENDIF ! of IF(itau.EQ.itaufin)
767
768	forward = .TRUE.
769	GO TO 1
770
771	ENDIF ! of IF (forward)
772
773	END IF ! of IF(.not.purmats)
774
775	STOP
776	END

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

Download in other formats:

Original Format