Context Navigation

gcm.F @ 1523

Last change on this file since 1523 was 1523, checked in by emillour, 10 years ago

All models: More updates to make planetary codes (+Earth) setups converge.

in dyn3d_common:
convmas.F => convmas.F90
enercin.F => enercin.F90
flumass.F => flumass.F90
massbar.F => massbar.F90
tourpot.F => tourpot.F90
vitvert.F => vitvert.F90

in misc:
move "q_sat" from "dyn3d_common" to "misc" (in Earth model, it is also called by the physics)
move "write_field" from "dyn3d_common" to "misc"(may be called from physics or dynamics and depends on neither).

in phy_common:
move "write_field_phy" here since it may be called from any physics package)
add module "regular_lonlat_mod" to store global information on lon-lat grid

in dynlonlat_phylonlat/phy*:
turn "iniphysiq.F90" into module "iniphysiq_mod.F90" (and of course adapt gcm.F[90] and 1D models accordingly)

File size: 22.2 KB

Line
1	PROGRAM gcm
2
3	use infotrac, only: iniadvtrac, nqtot, iadv
4	use control_mod, only: day_step, iperiod, iphysiq, ndynstep,
5	& nday_r, idissip, iconser, ecritstart,
6	& ecritphy
7	use filtreg_mod, only: inifilr
8	! use comgeomphy, only: initcomgeomphy
9	USE comvert_mod, ONLY: ap,bp
10	USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss,rad,g,r,cpp
11	USE logic_mod, ONLY: ecripar,forward,leapf,apphys,statcl,conser,
12	. apdiss,purmats,physic,apphys
13	USE temps_mod, ONLY: day_ini,day_end,dt,itaufin
14	USE iniphysiq_mod, ONLY: iniphysiq
15	IMPLICIT NONE
16
17	c ...... Version du 10/01/98 ..........
18
19	c avec coordonnees verticales hybrides
20	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
21
22	c=======================================================================
23	c
24	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
25	c -------
26	c
27	c Objet:
28	c ------
29	c
30	c GCM LMD nouvelle grille
31	c
32	c=======================================================================
33	c
34	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
35	c et possibilite d'appeler une fonction f(y) a derivee tangente
36	c hyperbolique a la place de la fonction a derivee sinusoidale.
37
38	c ... Possibilite de choisir le shema de Van-leer pour l'advection de
39	c q , en faisant iadv = 3 dans traceur (29/04/97) .
40	c
41	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
42	c
43	c-----------------------------------------------------------------------
44	c Declarations:
45	c -------------
46
47	#include "dimensions.h"
48	#include "paramet.h"
49	#include "comdissnew.h"
50	#include "comgeom.h"
51	!#include "control.h"
52	#include "netcdf.inc"
53	#include "tracstoke.h"
54	#include "sponge.h"
55	!#include"advtrac.h"
56
57	INTEGER*4 iday ! jour julien
58	REAL time ! Heure de la journee en fraction d''1 jour
59	REAL zdtvr
60
61	c variables dynamiques
62	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
63	real, dimension(ip1jmp1,llm) :: teta ! temperature potentielle
64	REAL,allocatable :: q(:,:,:) ! champs advectes
65	REAL ps(ip1jmp1) ! pression au sol
66	REAL pext(ip1jmp1) ! pression extensive
67	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
68	REAL pks(ip1jmp1) ! exner au sol
69	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
70	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
71	REAL masse(ip1jmp1,llm) ! masse d''air
72	REAL phis(ip1jmp1) ! geopotentiel au sol
73	REAL phi(ip1jmp1,llm) ! geopotentiel
74	REAL w(ip1jmp1,llm) ! vitesse verticale
75
76	c variables dynamiques intermediaire pour le transport
77	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
78
79	c variables dynamiques au pas -1
80	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
81
82	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
83	REAL massem1(ip1jmp1,llm)
84
85	c tendances dynamiques
86	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
87	REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
88	REAL,ALLOCATABLE :: dq(:,:,:)
89
90	c tendances de la dissipation
91	REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
92	REAL dhdis(ip1jmp1,llm)
93
94	c tendances physiques
95	REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
96	REAL dhfi(ip1jmp1,llm),dpfi(ip1jmp1)
97	REAL,ALLOCATABLE :: dqfi(:,:,:)
98
99	c variables pour le fichier histoire
100	REAL dtav ! intervalle de temps elementaire
101
102	REAL tppn(iim),tpps(iim),tpn,tps
103	c
104
105	INTEGER itau,itaufinp1,iav
106
107
108	EXTERNAL caldyn, traceur
109	EXTERNAL dissip,geopot,iniconst
110	EXTERNAL integrd,SCOPY
111	EXTERNAL inigeom
112	EXTERNAL exner_hyb,addit
113	EXTERNAL defrun_new, test_period
114	REAL SSUM
115	REAL time_0 , finvmaold(ip1jmp1,llm)
116
117	LOGICAL lafin
118	INTEGER ij,iq,l,ierr,numvanle,iapp_tracvl
119
120	REAL rdayvrai,rdaym_ini,rday_ecri
121	! LOGICAL first
122	REAL beta(ip1jmp1,llm)
123
124	LOGICAL offline ! Controle du stockage ds "fluxmass"
125	PARAMETER (offline=.false.)
126
127	character*20 modname
128	character*80 abort_message
129
130	LOGICAL tracer
131	data tracer/.true./
132	! INTEGER nq
133
134	C Calendrier
135	LOGICAL true_calendar
136	PARAMETER (true_calendar = .false.)
137
138	! flag to set/remove calls to groupeun
139	logical callgroupeun
140	parameter (callgroupeun = .false.)
141
142	c-----------------------------------------------------------------------
143	c variables pour l'initialisation de la physique :
144	c ------------------------------------------------
145	! INTEGER ngridmx
146	! PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
147	! REAL zcufi(ngridmx),zcvfi(ngridmx)
148	! REAL latfi(ngridmx),lonfi(ngridmx)
149	! REAL airefi(ngridmx)
150	! SAVE latfi, lonfi, airefi
151	! INTEGER i,j
152
153	c-----------------------------------------------------------------------
154	c Initialisations:
155	c ----------------
156
157	modname = 'gcm'
158	lafin = .FALSE.
159
160	c-----------------------------------------------------------------------
161	CALL defrun_new( 99, .TRUE. )
162
163	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
164	! FH 2008/05/02
165	! A nettoyer. On ne veut qu'une ou deux routines d'interface
166	! dynamique -> physique pour l'initialisation
167	!#ifdef CPP_PHYS
168	CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
169	! call initcomgeomphy ! now done in iniphysiq
170	!#endif
171	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
172
173	! Initialize tracers
174	call iniadvtrac(nqtot,numvanle)
175	! Allocation de la tableau q : champs advectes
176	allocate(q(ip1jmp1,llm,nqtot))
177	allocate(dq(ip1jmp1,llm,nqtot))
178	allocate(dqfi(ip1jmp1,llm,nqtot))
179
180	CALL dynetat0("start.nc",vcov,ucov,
181	. teta,q,masse,ps,phis,time_0)
182
183	! in case time_0 (because of roundoffs) is close to zero,
184	! set it to zero to avoid roundoff propagation issues
185	if ((time_0.gt.0.).and.(time_0.lt.(1./day_step))) then
186	write(,)"GCM: In start.nc, time=",time_0
187	write(,)" but day_step=",day_step
188	write(,)" and 1./day_step=",1./day_step
189	write(,)" fix this drift by setting time=0"
190	time_0=0.
191	endif
192
193
194	c on recalcule eventuellement le pas de temps
195
196	IF(MOD(day_step,iperiod).NE.0)
197	* STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
198
199	IF(MOD(day_step,iphysiq).NE.0)
200	* STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
201
202	zdtvr = daysec/REAL(day_step)
203	IF(dtvr.NE.zdtvr) THEN
204	PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
205	ENDIF
206
207	c nombre d'etats dans les fichiers demarrage et histoire
208
209	dtvr = zdtvr
210	CALL iniconst
211	CALL inigeom
212
213	CALL inifilr
214
215	c
216	c ...... P.Le Van ( modif le 29/04/97 ) .........
217	c
218
219	CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh ,
220	* tetagdiv, tetagrot , tetatemp )
221	c
222
223	c-----------------------------------------------------------------------
224	c Initialisation de la physique :
225	c -------------------------------
226
227	! IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
228	! latfi(1)=rlatu(1)
229	! lonfi(1)=0.
230	! zcufi(1) = cu(1)
231	! zcvfi(1) = cv(1)
232	! DO j=2,jjm
233	! DO i=1,iim
234	! latfi((j-2)*iim+1+i)= rlatu(j)
235	! lonfi((j-2)*iim+1+i)= rlonv(i)
236	! zcufi((j-2)iim+1+i) = cu((j-1)iip1+i)
237	! zcvfi((j-2)iim+1+i) = cv((j-1)iip1+i)
238	! ENDDO
239	! ENDDO
240	! latfi(ngridmx)= rlatu(jjp1)
241	! lonfi(ngridmx)= 0.
242	! zcufi(ngridmx) = cu(ip1jm+1)
243	! zcvfi(ngridmx) = cv(ip1jm-iim)
244	!
245	! ! build airefi(), mesh area on physics grid
246	! CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
247	! ! Poles are single points on physics grid
248	! airefi(1)=airefi(1)*iim
249	! airefi(ngridmx)=airefi(ngridmx)*iim
250
251	! Initialisation de la physique: pose probleme quand on tourne
252	! SANS physique, car iniphysiq.F est dans le repertoire phy[]...
253	! Il faut une cle CPP_PHYS
254	!#ifdef CPP_PHYS
255	! CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys,
256	CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys,
257	& rlatu,rlatv,rlonu,rlonv,
258	& aire,cu,cv,rad,g,r,cpp,
259	& 1)
260	! & iflag_phys)
261	!#endif
262	! call_iniphys=.false.
263	! ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
264	! call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
265	! . latfi,lonfi,airefi,rad,g,r,cpp)
266
267	call dump2d(iip1,jjp1,ps,'PRESSION SURFACE')
268
269	c
270	c numero de stockage pour les fichiers de redemarrage:
271
272	c-----------------------------------------------------------------------
273	c temps de depart et de fin:
274	c --------------------------
275
276	itau = 0
277	iday = day_ini+itau/day_step
278	time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
279	IF(time.GT.1.) THEN
280	time = time-1.
281	iday = iday+1
282	ENDIF
283	if (ndynstep .gt. 0) then
284	itaufin = ndynstep
285	else
286	itaufin=nint(nday_r*day_step) ! nint() to avoid problematic roundoffs
287	endif
288	! check that this is compatible with call sequence dyn/phys/dissip
289	! i.e. that itaufin is a multiple of iphysiq and idissip
290	if ((modulo(itaufin,iphysiq).ne.0).or.
291	& (modulo(itaufin,idissip).ne.0)) then
292	if (ndynstep .gt. 0) then
293	write(*,'(A,I5)')
294	& "gcm: Problem: incompatibility between ndynstep=",ndynstep
295	else
296	write(*,'((A,F9.2),2(A,I5))')
297	& "gcm: Problem: incompatibility between nday=",nday_r,
298	& " day_step=",day_step," which imply itaufin=",itaufin
299	endif
300	write(*,'(2(A,I5))')
301	& " whereas iphysiq=",iphysiq," and idissip=",
302	& idissip
303	stop
304	endif
305	! write(,)"gcm: itaufin=",itaufin
306	c ********************************
307	c itaufin = 120 ! temporaire !!
308	c ********************************
309	itaufinp1 = itaufin +1
310
311	if (ndynstep .gt. 0) then
312	day_end = day_ini
313	& + floor(float(ndynstep)/float(day_step)+time_0)
314	else
315	day_end = day_ini + floor(nday_r+time_0)
316	endif
317	PRINT 300, itau,itaufin,day_ini,day_end
318	300 FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x,
319	. 'c''est a dire du jour',i7,3x,'au jour',i7//)
320
321	CALL dynredem0("restart.nc",day_ini,phis)
322
323	ecripar = .TRUE.
324
325	dtav = iperiod*dtvr/daysec
326
327
328	c Quelques initialisations pour les traceurs
329	dq(:,:,:)=0
330	c istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
331	c istphy=istdyn/iphysiq
332
333	write(,) "gcm: callgroupeun set to:",callgroupeun
334	c-----------------------------------------------------------------------
335	c Debut de l'integration temporelle:
336	c ----------------------------------
337
338	1 CONTINUE
339	c
340	c TN 09/2012. To ensure "1+1=2" in dynamical core :
341	c update atmospheric pressure IN the main loop
342	CALL pression ( ip1jmp1, ap, bp, ps, p )
343	CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf )
344
345	IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
346	CALL test_period ( ucov,vcov,teta,q,p,phis )
347	write(,)' GCM ---- Test_period apres continue OK ! -----',
348	& ' itau: ',itau
349	ENDIF
350
351	if (callgroupeun) then
352	call groupeun(jjp1,llm,ucov,.true.)
353	call groupeun(jjm,llm,vcov,.true.)
354	call groupeun(jjp1,llm,teta,.true.)
355	call groupeun(jjp1,llm,masse,.true.)
356	call groupeun(jjp1,1,ps,.false.)
357	endif
358
359	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
360	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
361	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
362	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
363	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
364
365	forward = .TRUE.
366	leapf = .FALSE.
367	dt = dtvr
368
369	c ... P.Le Van .26/04/94 ....
370
371	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
372	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
373
374
375	2 CONTINUE
376
377	c-----------------------------------------------------------------------
378
379	c date:
380	c -----
381
382	! write(,) 'GCM: itau=',itau
383
384	c gestion des appels de la physique et des dissipations:
385	c ------------------------------------------------------
386	c
387	c ... P.Le Van ( 6/02/95 ) ....
388
389	apphys = .FALSE.
390	statcl = .FALSE.
391	conser = .FALSE.
392	apdiss = .FALSE.
393
394	IF( purmats ) THEN
395	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
396	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
397	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
398	$ .AND. physic ) apphys = .TRUE.
399	ELSE
400	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
401	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
402	IF( MOD(itau+1,iphysiq).EQ.0. AND. physic ) apphys = .TRUE.
403	END IF
404
405	c-----------------------------------------------------------------------
406	c calcul des tendances dynamiques:
407	c --------------------------------
408
409	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
410	c
411	c
412	CALL caldyn
413	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
414	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
415
416
417	c-----------------------------------------------------------------------
418	c calcul des tendances advection des traceurs (dont l'humidite)
419	c -------------------------------------------------------------
420
421	if (tracer) then
422	IF( forward. OR . leapf ) THEN
423
424	DO iq = 1, nqtot
425	c
426	IF ( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 ) THEN
427	CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
428
429	ELSE IF( iq.EQ. nqtot ) THEN
430	c
431	iapp_tracvl = 5
432	c
433	cccc iapp_tracvl est la frequence en pas du groupement des flux
434	cccc de masse pour Van-Leer dans la routine tracvl .
435	c
436
437	CALL vanleer(numvanle,iapp_tracvl,nqtot,q,pbaru,pbarv,
438	* p, masse, dq, iadv(1), teta, pk )
439
440	c
441	c ... Modif F.Codron ....
442	c
443	ENDIF
444	c
445	ENDDO
446	C
447	c IF (offline) THEN
448	C maf stokage du flux de masse pour traceurs OFF-LINE
449
450	c CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
451	c . time_step, itau)
452
453	c ENDIF
454	c
455	ENDIF
456	END IF ! tracer
457
458
459	c-----------------------------------------------------------------------
460	c integrations dynamique et traceurs:
461	c ----------------------------------
462
463	CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
464	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
465	$ finvmaold )
466
467	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
468	c
469	c-----------------------------------------------------------------------
470	c calcul des tendances physiques:
471	c -------------------------------
472	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
473	c
474	IF( purmats ) THEN
475	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
476	ELSE
477	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
478	ENDIF
479	c
480	c
481	IF( apphys ) THEN
482	c
483	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
484	c
485
486	CALL pression ( ip1jmp1, ap, bp, ps, p )
487	CALL exner_hyb( ip1jmp1, ps, p,beta,pks, pk, pkf )
488
489	rdaym_ini = itau * dtvr / daysec + time_0
490	rdayvrai = rdaym_ini + day_ini
491
492	IF ( ecritphy.LT.1. ) THEN
493	rday_ecri = rdaym_ini
494	ELSE
495	rday_ecri = INT( rdayvrai )
496	ENDIF
497	c
498	CALL calfis( nqtot, lafin ,rdayvrai,rday_ecri,time ,
499	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
500	$ du,dv,dteta,dq,w, dufi,dvfi,dhfi,dqfi,dpfi,tracer)
501
502
503	c ajout des tendances physiques:
504	c ------------------------------
505	CALL addfi( nqtot, dtphys, leapf, forward ,
506	$ ucov, vcov, teta , q ,ps , masse,
507	$ dufi, dvfi, dhfi , dqfi ,dpfi )
508	c
509	ENDIF
510
511	CALL pression ( ip1jmp1, ap, bp, ps, p )
512
513	CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf )
514	c ----------------------------------------------------------
515
516	c
517	c
518	c dissipation horizontale et verticale des petites echelles:
519	c ----------------------------------------------------------
520
521
522	IF(apdiss) THEN
523
524	c Sponge layer
525	c ~~~~~~~~~~~~
526	DO ij=1, ip1jmp1
527	pext(ij)=ps(ij)*aire(ij)
528	ENDDO
529	IF (callsponge) THEN
530	CALL sponge(ucov,vcov,teta,pext,dtdiss,mode_sponge)
531	ENDIF
532
533	c Dissipation horizontale
534	c ~~~~~~~~~~~~~~~~~~~~~~~
535	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
536
537	CALL addit( ijp1llm,ucov ,dudis,ucov )
538	CALL addit( ijmllm ,vcov ,dvdis,vcov )
539	CALL addit( ijp1llm,teta ,dhdis,teta )
540
541
542	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
543	c ... Calcul de la valeur moyenne, unique de h aux poles .....
544	c
545
546	DO l = 1, llm
547	DO ij = 1,iim
548	tppn(ij) = aire( ij ) * teta( ij ,l)
549	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
550	ENDDO
551
552	tpn = SSUM(iim,tppn,1)/apoln
553	tps = SSUM(iim,tpps,1)/apols
554
555	DO ij = 1, iip1
556	teta( ij ,l) = tpn
557	teta(ij+ip1jm,l) = tps
558	ENDDO
559	ENDDO
560
561	DO ij = 1,iim
562	tppn(ij) = aire( ij ) * ps ( ij )
563	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
564	ENDDO
565	tpn = SSUM(iim,tppn,1)/apoln
566
567	tps = SSUM(iim,tpps,1)/apols
568
569	DO ij = 1, iip1
570	ps( ij ) = tpn
571	ps(ij+ip1jm) = tps
572	ENDDO
573
574
575	END IF
576
577	c ********************************************************************
578	c ********************************************************************
579	c .... fin de l'integration dynamique et physique pour le pas itau ..
580	c ********************************************************************
581	c ********************************************************************
582
583	c preparation du pas d'integration suivant ......
584
585	IF ( .NOT.purmats ) THEN
586	c ........................................................
587	c .............. schema matsuno + leapfrog ..............
588	c ........................................................
589
590	IF(forward. OR. leapf) THEN
591	itau= itau + 1
592	iday= day_ini+itau/day_step
593	time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
594	!IF(time.GT.1.) THEN
595	IF(time.GE.1.) THEN
596	time = time-1.
597	iday = iday+1
598	ENDIF
599	ENDIF
600
601
602	IF( itau. EQ. itaufinp1 ) then
603	abort_message = 'Simulation finished'
604	call abort_gcm(modname,abort_message,0)
605	ENDIF
606	c-----------------------------------------------------------------------
607	c ecriture du fichier histoire moyenne:
608	c -------------------------------------
609
610	c IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
611	c IF(itau.EQ.itaufin) THEN
612	c iav=1
613	c ELSE
614	c iav=0
615	c ENDIF
616	c CALL writedynav(histaveid, nqtot, itau,vcov ,
617	c , ucov,teta,pk,phi,q,masse,ps,phis)
618	c ENDIF
619
620	c-----------------------------------------------------------------------
621
622	IF( (ecritstart.GT.0) .and. (MOD(itau,ecritstart).EQ.0)
623	. .or. (itau.EQ.itaufin) ) THEN
624
625	!write(,)' GCM: Appel test_period avant redem ; itau=',itau
626	CALL test_period ( ucov,vcov,teta,q,p,phis )
627
628	write(*,'(A,I7,A,F12.5)')
629	. 'GCM: Ecriture du fichier restart ; itau= ',itau,
630	. ' date=',REAL(itau)/REAL(day_step)
631	CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
632	. vcov,ucov,teta,q,masse,ps)
633
634	CLOSE(99)
635
636	ENDIF
637
638
639	c-----------------------------------------------------------------------
640	c gestion de l'integration temporelle:
641	c ------------------------------------
642
643	IF( MOD(itau,iperiod).EQ.0 ) THEN
644	GO TO 1
645	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
646
647	IF( forward ) THEN
648	c fin du pas forward et debut du pas backward
649
650	forward = .FALSE.
651	leapf = .FALSE.
652	GO TO 2
653
654	ELSE
655	c fin du pas backward et debut du premier pas leapfrog
656
657	leapf = .TRUE.
658	dt = 2.*dtvr
659	GO TO 2
660	END IF
661	ELSE
662
663	c ...... pas leapfrog .....
664
665	leapf = .TRUE.
666	dt = 2.*dtvr
667	GO TO 2
668	END IF
669
670	ELSE
671
672	c ........................................................
673	c .............. schema matsuno ...............
674	c ........................................................
675	IF( forward ) THEN
676
677	itau = itau + 1
678	iday = day_ini+itau/day_step
679	time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
680
681	IF(time.GE.1.) THEN
682	time = time-1.
683	iday = iday+1
684	ENDIF
685
686	forward = .FALSE.
687	IF( itau. EQ. itaufinp1 ) then
688	abort_message = 'Simulation finished'
689	call abort_gcm(modname,abort_message,0)
690	ENDIF
691	GO TO 2
692
693	ELSE
694
695	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
696	IF(itau.EQ.itaufin) THEN
697	iav=1
698	ELSE
699	iav=0
700	ENDIF
701	c CALL writedynav(histaveid, nqtot, itau,vcov ,
702	c , ucov,teta,pk,phi,q,masse,ps,phis)
703
704	ENDIF
705
706
707	IF( (ecritstart.GT.0) .and. (MOD(itau,ecritstart).EQ.0)
708	. .or. (itau.EQ.itaufin) ) THEN
709
710	CALL dynredem1("restart.nc",
711	. REAL(itau)/REAL(day_step),
712	. vcov,ucov,teta,q,masse,ps)
713
714	forward = .TRUE.
715	GO TO 1
716
717	ENDIF
718
719
720	ENDIF
721
722	END IF
723
724	STOP
725	END

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

Download in other formats:

Original Format