Context Navigation

leapfrog.F @ 493

Last change on this file since 493 was 492, checked in by emillour, 14 years ago
Common dynamics: updates to keep up with LMDZ5 Earth (rev 1605) See file "DOC/chantiers/commit_importants.log" for details. EM
File size: 27.5 KB

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

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

Download in other formats:

Original Format