Context Navigation

leapfrog.F @ 1549

Last change on this file since 1549 was 1529, checked in by Laurent Fairhead, 13 years ago

Necessary modifications to use the model in an aquaplanet or terraplanet configuration

Modifications nécessaires pour l'utilisation du moèle en configuration aquaplanète
ou terraplanète

Dans la dynamique:
Changement pour le paramètre iflag_phys:

auparavant : 0 pas de physique, 1 phyique, 2 rappel
ici on ajoute iflag_phys>=100 pour les aqua et terra

Du coup, dans leapfrog.F et gcm.F on appelle la physique si
iflag_phys=1.or.iflag_phys>=100
Dans iniacademic, on initialise si iflag_phys>=2 au lieu de =2
Dans gcm.F, on appelle en plus iniaqua (sous une clef CPP_EARTH)
Dans iniacademic, on met ps=108080 pour les aqua et terra pour répondre
à une specification CFMIP.

Dans la physique:
On ajoute phyaqua.F qui contient :
iniaqua : initialise les startphy.nc et limit.nc pour la phyique
zenang_an : calcule un ensolleillement moyen sur l'année en fonction de

la latiude (A. Campoy).

profil_sst : qui calcule différents profils latitudinaux de SSTs.
writelim : écrit le fichier limit.nc

Dans physiq.F
On ajoute la possibilité d'appeller un calcul d'ensoleillement moyen
sur l'année quand solarlong0=1000.

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

File size: 24.4 KB

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

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

Download in other formats:

Original Format