Context Navigation

gcm.F @ 2167

Last change on this file since 2167 was 516, checked in by lmdzadmin, 21 years ago
Deplacement de commentaires pour raisons d'optimisation LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 24.6 KB

Rev	Line
[2]	1
	2	PROGRAM gcm
	3
[57]	4	USE IOIPSL
	5
[2]	6	IMPLICIT NONE
	7
	8	c ...... Version du 10/01/98 ..........
	9
	10	c avec coordonnees verticales hybrides
	11	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
	12
	13	c=======================================================================
	14	c
	15	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
	16	c -------
	17	c
	18	c Objet:
	19	c ------
	20	c
	21	c GCM LMD nouvelle grille
	22	c
	23	c=======================================================================
[232]	24	c
	25	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
	26	c et possibilite d'appeler une fonction f(y) a derivee tangente
	27	c hyperbolique a la place de la fonction a derivee sinusoidale.
[2]	28
[232]	29	c ... Possibilite de choisir le shema de Van-leer pour l'advection de
	30	c q , en faisant iadv = 3 dans traceur (29/04/97) .
[2]	31	c
[232]	32	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
	33	c
[2]	34	c-----------------------------------------------------------------------
	35	c Declarations:
	36	c -------------
	37
	38	#include "dimensions.h"
	39	#include "paramet.h"
	40	#include "comconst.h"
	41	#include "comdissnew.h"
	42	#include "comvert.h"
	43	#include "comgeom.h"
	44	#include "logic.h"
	45	#include "temps.h"
	46	#include "control.h"
	47	#include "ener.h"
	48	#include "netcdf.inc"
	49	#include "description.h"
	50	#include "serre.h"
	51	#include "tracstoke.h"
	52
	53
	54	INTEGER longcles
	55	PARAMETER ( longcles = 20 )
	56	REAL clesphy0( longcles )
	57	SAVE clesphy0
	58
	59	INTEGER*4 iday ! jour julien
	60	REAL time ! Heure de la journee en fraction d'1 jour
	61	REAL zdtvr
	62	INTEGER nbetatmoy, nbetatdem,nbetat
	63
	64	c variables dynamiques
	65	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
	66	REAL teta(ip1jmp1,llm) ! temperature potentielle
	67	REAL q(ip1jmp1,llm,nqmx) ! champs advectes
	68	REAL ps(ip1jmp1) ! pression au sol
	69	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
	70	REAL pks(ip1jmp1) ! exner au sol
	71	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
	72	REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
	73	REAL masse(ip1jmp1,llm) ! masse d'air
	74	REAL phis(ip1jmp1) ! geopotentiel au sol
	75	REAL phi(ip1jmp1,llm) ! geopotentiel
	76	REAL w(ip1jmp1,llm) ! vitesse verticale
	77
	78	c variables dynamiques intermediaire pour le transport
	79	REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
	80
	81	c variables dynamiques au pas -1
	82	REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
	83	REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
	84	REAL massem1(ip1jmp1,llm)
	85
	86	c tendances dynamiques
	87	REAL dv(ip1jm,llm),du(ip1jmp1,llm)
	88	REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
	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),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
	97
	98	c variables pour le fichier histoire
	99	REAL dtav ! intervalle de temps elementaire
	100
	101	REAL tppn(iim),tpps(iim),tpn,tps
	102	c
	103	INTEGER iadv(nqmx) ! indice schema de transport pour le traceur iq
	104
	105	INTEGER itau,itaufinp1,iav
	106
	107	EXTERNAL caldyn, traceur
	108	EXTERNAL dissip,geopot,iniconst,inifilr
	109	EXTERNAL integrd,SCOPY
	110	EXTERNAL iniav,writeav,writehist
	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,numvanle,iapp_tracvl
	119
	120
[57]	121	INTEGER fluxid, fluxvid,fluxdid
[2]	122	integer histid, histvid, histaveid
	123	real time_step, t_wrt, t_ops
	124
	125	REAL rdayvrai,rdaym_ini,rday_ecri
	126	LOGICAL first
	127	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
[403]	128	c+jld variables test conservation energie
	129	REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
	130	C Tendance de la temp. potentiel d (theta)/ d t due a la
	131	C tansformation d'energie cinetique en energie thermique
	132	C cree par la dissipation
	133	REAL dhecdt(ip1jmp1,llm)
	134	REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
	135	REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
	136	CHARACTER*15 ztit
[459]	137	INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
	138	SAVE ip_ebil_dyn
	139	DATA ip_ebil_dyn/0/
[403]	140	c-jld
[2]	141
	142	LOGICAL offline ! Controle du stockage ds "fluxmass"
[99]	143	PARAMETER (offline=.false.)
[2]	144
	145	character*80 dynhist_file, dynhistave_file
	146	character*20 modname
	147	character*80 abort_message
	148
[57]	149	C Calendrier
	150	LOGICAL true_calendar
	151	PARAMETER (true_calendar = .false.)
[232]	152	C Run nudge
	153	LOGICAL ok_nudge
	154	PARAMETER (ok_nudge = .false.)
[495]	155	c
[57]	156
[2]	157	c-----------------------------------------------------------------------
	158	c Initialisations:
	159	c ----------------
	160
[99]	161	abort_message = 'last timestep reached'
[2]	162	modname = 'gcm'
	163	descript = 'Run GCM LMDZ'
	164	lafin = .FALSE.
	165	dynhist_file = 'dyn_hist.nc'
	166	dynhistave_file = 'dyn_hist_ave.nc'
[57]	167
	168	if (true_calendar) then
	169	call ioconf_calendar('gregorian')
	170	else
	171	call ioconf_calendar('360d')
	172	endif
[2]	173	c-----------------------------------------------------------------------
	174	c
	175	c .... Choix des shemas d'advection pour l'eau et les traceurs ...
	176	c ...................................................................
	177	c
	178	c iadv = 1 shema transport type "humidite specifique LMD"
	179	c iadv = 2 shema amont
	180	c iadv = 3 shema Van-leer
[232]	181	c iadv = 4 schema Van-leer + humidite specifique
	182	c Modif F.Codron
[2]	183	c
	184	c dans le tableau q(ij,l,iq) , iq = 1 pour l'eau vapeur
	185	c , iq = 2 pour l'eau liquide
	186	c et eventuellement , iq = 3, nqmx pour les autres traceurs
	187	c
	188	DO iq = 1, nqmx
	189	iadv( iq ) = 3
	190	ENDDO
	191	c
[468]	192	c iadv(1): choix pour l'eau vap. et iadv(2) : choix pour l'eau liq.
	193	iadv( 1 ) = 4
	194	c
[232]	195	IF( iadv(1).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
	196	* ,' pour la vapeur d''eau'
	197	IF( iadv(1).EQ.2 ) PRINT *,' Choix du shema amont',' pour la'
	198	* ,' vapeur d''eau '
	199	IF( iadv(1).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour'
	200	* ,'la vapeur d''eau'
	201	IF( iadv(1).EQ.4 ) PRINT *,' Choix du shema Van-Leer + humidite'
	202	* ,' specifique pour la vapeur d''eau'
	203	c
	204	IF( iadv(1).LE.0.OR.iadv(1).GT.4 ) THEN
	205	PRINT *,' Erreur dans le choix de iadv (1).Corriger et repasser
	206	* , car iadv(1) = ', iadv(1)
	207	CALL ABORT
	208	ENDIF
	209
	210	DO iq = 2, nqmx
[2]	211	IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
	212	* ,' pour le traceur no ', iq
	213	IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema amont',' pour le'
	214	* ,' traceur no ', iq
	215	IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema Van-Leer ',' pour'
	216	* ,'le traceur no ', iq
[232]	217
	218	IF( iadv(iq).EQ.4 ) THEN
	219	PRINT *,' Le shema Van-Leer + humidite specifique ',
	220	* ' est uniquement pour la vapeur d eau .'
	221	PRINT *,' Corriger iadv( ',iq, ') et repasser ! '
	222	CALL ABORT
	223	ENDIF
	224
	225	IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 ) THEN
[2]	226	PRINT *,' Erreur dans le choix de iadv . Corriger et repasser
	227	* . '
[232]	228	CALL ABORT
[2]	229	ENDIF
	230	ENDDO
	231	c
	232	first = .TRUE.
	233	numvanle = nqmx + 1
	234	DO iq = 1, nqmx
[232]	235	IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN
[2]	236	numvanle = iq
	237	first = .FALSE.
	238	ENDIF
	239	ENDDO
	240	c
	241	DO iq = 1, nqmx
[232]	242	IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle ) THEN
[2]	243	PRINT *,' Il y a discontinuite dans le choix du shema de ',
	244	* 'Van-leer pour les traceurs . Corriger et repasser . '
[232]	245	CALL ABORT
	246	ENDIF
	247	IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 ) THEN
[2]	248	PRINT *,' Le choix de iadv est errone pour le traceur ',
	249	* iq
	250	STOP
	251	ENDIF
	252	ENDDO
	253	c
	254
	255	CALL dynetat0("start.nc",nqmx,vcov,ucov,
	256	. teta,q,masse,ps,phis, time_0)
	257
	258
	259	c OPEN( 99,file ='run.def',status='old',form='formatted')
[228]	260	c CALL defrun_new( 99, .TRUE. , clesphy0 )
	261	CALL conf_gcm( 99, .TRUE. , clesphy0 )
[2]	262	c CLOSE(99)
	263
	264	c on recalcule eventuellement le pas de temps
	265
	266	IF(MOD(day_step,iperiod).NE.0)
	267	* STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
	268
	269	IF(MOD(day_step,iphysiq).NE.0)
	270	* STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
	271
	272	zdtvr = daysec/FLOAT(day_step)
	273	IF(dtvr.NE.zdtvr) THEN
	274	PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
	275	ENDIF
[431]	276	C
	277	C on remet le calendrier à zero
	278	c
	279	if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
	280	write(,)' Attention les dates initiales lues dans le fichier'
	281	write(,)' restart ne correspondent pas a celles lues dans '
	282	write(,)' gcm.def'
	283	if (raz_date .ne. 1) then
	284	write(,)' On garde les dates du fichier restart'
	285	else
	286	annee_ref = anneeref
	287	day_ref = dayref
[477]	288	day_ini = dayref
[431]	289	itau_dyn = 0
	290	itau_phy = 0
[474]	291	time_0 = 0.
[431]	292	write(,)' On reinitialise a la date lue dans gcm.def'
	293	endif
[474]	294	ELSE
	295	raz_date = 0
[431]	296	endif
[2]	297
	298	c nombre d'etats dans les fichiers demarrage et histoire
	299	nbetatdem = nday / iecri
	300	nbetatmoy = nday / periodav + 1
	301
	302	dtvr = zdtvr
	303	CALL iniconst
	304	CALL inigeom
	305
	306	CALL inifilr
	307	c
	308	c ...... P.Le Van ( modif le 29/04/97 ) .........
	309	c
	310	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
	311	* tetagdiv, tetagrot , tetatemp )
[403]	312	c+jld
	313	C initialisation constantes thermo utilisees dans diagedyn
[459]	314	IF (ip_ebil_dyn.ge.1 ) THEN
[403]	315	CALL suphec
	316	ENDIF
	317	c-jld
[2]	318	c
	319
	320	CALL pression ( ip1jmp1, ap, bp, ps, p )
	321	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	322	c
	323
	324	c numero de stockage pour les fichiers de redemarrage:
	325
	326	c-----------------------------------------------------------------------
	327	c temps de depart et de fin:
	328	c --------------------------
	329
	330	itau = 0
	331	iday = day_ini+itau/day_step
	332	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	333	IF(time.GT.1.) THEN
	334	time = time-1.
	335	iday = iday+1
	336	ENDIF
	337	itaufin = nday*day_step
	338	itaufinp1 = itaufin +1
	339
	340	day_end = day_ini + nday
	341	PRINT 300, itau,itaufin,day_ini,day_end
	342
[353]	343	CALL dynredem0("restart.nc", day_end, phis, nqmx)
[2]	344
	345	ecripar = .TRUE.
	346
	347	time_step = zdtvr
	348	t_ops = iecri * daysec
	349	t_wrt = iecri * daysec
[353]	350	C CALL inithist(dynhist_file,day_ref,annee_ref,time_step,
[295]	351	c . t_ops, t_wrt, nqmx, histid, histvid)
[2]	352
	353	t_ops = iperiod * time_step
	354	t_wrt = periodav * daysec
[353]	355	CALL initdynav(dynhistave_file,day_ref,annee_ref,time_step,
[2]	356	. t_ops, t_wrt, nqmx, histaveid)
	357
	358	dtav = iperiod*dtvr/daysec
	359
	360
	361	c Quelques initialisations pour les traceurs
	362	call initial0(ijp1llm*nqmx,dq)
	363	istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
	364	istphy=istdyn/iphysiq
[516]	365	print,'WARNING*: Pour des rasions d optimisations'
	366	print,'WARNING*: iapp_tracvl = iperiod'
	367	print,'WARNING*: Appel vanleer avec 2 traceurs'
	368	print*,' et non nqmx'
[2]	369
	370	c-----------------------------------------------------------------------
	371	c Debut de l'integration temporelle:
	372	c ----------------------------------
	373
[495]	374
[2]	375	1 CONTINUE
[232]	376
[403]	377
	378
[232]	379	if (ok_nudge) then
	380	call nudge(itau,ucov,vcov,teta,masse,ps)
	381	endif
[2]	382	c
	383	c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
	384	CALL test_period ( ucov,vcov,teta,q,p,phis )
[462]	385	c PRINT *,' ---- Test_period apres continue OK ! -----', itau
[2]	386	c ENDIF
	387	c
	388	CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
	389	CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
	390	CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
	391	CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
	392	CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
	393
	394	forward = .TRUE.
	395	leapf = .FALSE.
	396	dt = dtvr
	397
	398	c ... P.Le Van .26/04/94 ....
	399
	400	CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
	401	CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
	402
	403
	404	2 CONTINUE
	405
	406	c-----------------------------------------------------------------------
	407
	408	c date:
	409	c -----
	410
	411
	412	c gestion des appels de la physique et des dissipations:
	413	c ------------------------------------------------------
	414	c
	415	c ... P.Le Van ( 6/02/95 ) ....
	416
	417	apphys = .FALSE.
	418	statcl = .FALSE.
	419	conser = .FALSE.
	420	apdiss = .FALSE.
	421
	422	IF( purmats ) THEN
	423	IF( MOD(itau,iconser) .EQ.0.AND. forward ) conser = .TRUE.
	424	IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
	425	IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
	426	$ .AND. physic ) apphys = .TRUE.
	427	ELSE
	428	IF( MOD(itau ,iconser) .EQ. 0 ) conser = .TRUE.
	429	IF( MOD(itau+1,idissip) .EQ. 0 ) apdiss = .TRUE.
	430	IF( MOD(itau+1,iphysiq).EQ.0. AND. physic ) apphys = .TRUE.
	431	END IF
	432
	433	c-----------------------------------------------------------------------
	434	c calcul des tendances dynamiques:
	435	c --------------------------------
	436
	437	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
	438	c
	439	c
	440	CALL caldyn
	441	$ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
	442	$ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
	443
	444	c-----------------------------------------------------------------------
	445	c calcul des tendances advection des traceurs (dont l'humidite)
	446	c -------------------------------------------------------------
	447
	448	IF( forward. OR . leapf ) THEN
	449	c
	450	DO 15 iq = 1, nqmx
	451	c
	452	IF( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 ) THEN
	453	CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
	454
	455	ELSE IF( iq.EQ. nqmx ) THEN
	456	c
[495]	457	c iapp_tracvl = 5
	458	iapp_tracvl = iperiod
[2]	459	c
	460	cccc iapp_tracvl est la frequence en pas du groupement des flux
	461	cccc de masse pour Van-Leer dans la routine tracvl .
	462	c
[495]	463	c CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
	464	CALL vanleer(numvanle,iapp_tracvl,2,q,pbaru,pbarv,
[232]	465	* p, masse, dq, iadv(1), teta, pk )
[2]	466	c
[232]	467	c ... Modif F.Codron ....
	468	c
[2]	469	ENDIF
	470	c
	471	15 CONTINUE
	472	C
	473	IF (offline) THEN
	474	Cmaf stokage du flux de masse pour traceurs OFF-LINE
[57]	475
	476	CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
[232]	477	. time_step, itau)
[57]	478
[232]	479
[2]	480	ENDIF
	481	c
	482	ENDIF
	483
	484
	485	c-----------------------------------------------------------------------
	486	c integrations dynamique et traceurs:
	487	c ----------------------------------
	488
	489
	490	CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
	491	$ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
	492	$ finvmaold )
	493
[495]	494
[2]	495	c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
	496	c
	497	c-----------------------------------------------------------------------
	498	c calcul des tendances physiques:
	499	c -------------------------------
	500	c ######## P.Le Van ( Modif le 6/02/95 ) ###########
	501	c
	502	IF( purmats ) THEN
	503	IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
	504	ELSE
	505	IF( itau+1. EQ. itaufin ) lafin = .TRUE.
	506	ENDIF
	507	c
	508	c
	509	IF( apphys ) THEN
	510	c
	511	c ....... Ajout P.Le Van ( 17/04/96 ) ...........
	512	c
	513
	514	CALL pression ( ip1jmp1, ap, bp, ps, p )
	515	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
	516
	517	rdaym_ini = itau * dtvr / daysec
	518	rdayvrai = rdaym_ini + day_ini
	519
	520	IF ( ecritphy.LT.1. ) THEN
	521	rday_ecri = rdaym_ini
	522	ELSE
	523	rday_ecri = NINT( rdayvrai )
	524	ENDIF
	525	c
[99]	526
	527	c rajout debug
[113]	528	c lafin = .true.
[403]	529	c+jld
[459]	530	IF (ip_ebil_dyn.ge.1 ) THEN
[403]	531	ztit='bil dyn'
	532	CALL diagedyn(ztit,2,1,1,dtphys
	533	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
	534	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
	535	ENDIF
	536	c-jld
[2]	537	CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time ,
	538	$ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
	539	$ du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi )
	540
	541	c ajout des tendances physiques:
	542	c ------------------------------
	543	CALL addfi( nqmx, dtphys, leapf, forward ,
	544	$ ucov, vcov, teta , q ,ps ,
	545	$ dufi, dvfi, dhfi , dqfi ,dpfi )
	546	c
[403]	547	c+jld
[459]	548	IF (ip_ebil_dyn.ge.1 ) THEN
[403]	549	ztit='bil phys'
	550	CALL diagedyn(ztit,2,1,1,dtphys
	551	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
	552	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
	553	ENDIF
	554	c-jld
[2]	555	ENDIF
	556
	557	CALL pression ( ip1jmp1, ap, bp, ps, p )
	558	CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
	559
	560	c-----------------------------------------------------------------------
	561	c
	562	c
	563	c dissipation horizontale et verticale des petites echelles:
	564	c ----------------------------------------------------------
	565
	566	IF(apdiss) THEN
[403]	567	c+jld
[459]	568	IF (ip_ebil_dyn.ge.2 ) THEN
[403]	569	ztit='avant dissip'
	570	CALL diagedyn(ztit,2,2,2,dtvr
	571	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
	572	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
	573	ENDIF
	574	CALL covcont ( llm , ucov , vcov , ucont, vcont )
	575	CALL enercin ( vcov , ucov , vcont , ucont , ecin0 )
	576	c-jld
[2]	577
	578	CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
	579
	580	CALL addit( ijp1llm,ucov ,dudis,ucov )
	581	CALL addit( ijmllm ,vcov ,dvdis,vcov )
[403]	582	c+jld
	583	CALL covcont ( llm , ucov , vcov , ucont, vcont )
	584	CALL enercin ( vcov , ucov , vcont , ucont , ecin )
	585	C
	586	C On rajoute la tendance due a la transform. Ec -> E therm. cree
	587	C lors de la dissipation
	588	DO l = 1, llm
	589	DO ij = 1, ip1jmp1
	590	dhecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
	591	dhdis(ij,l) = dhdis(ij,l) + dhecdt(ij,l)
	592	ENDDO
	593	ENDDO
	594	c-jld
[2]	595	CALL addit( ijp1llm,teta ,dhdis,teta )
[403]	596	c+jld
	597	c
[459]	598	IF (ip_ebil_dyn.ge.2 ) THEN
[403]	599	ztit='apres dissip'
	600	CALL diagedyn(ztit,2,2,2,dtdiss
	601	e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
	602	s , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
	603	ENDIF
	604	c-jld
[2]	605
	606	c ....... P. Le Van ( ajout le 17/04/96 ) ...........
	607	c ... Calcul de la valeur moyenne, unique de h aux poles .....
	608	c
	609
	610	DO l = 1, llm
	611	DO ij = 1,iim
	612	tppn(ij) = aire( ij ) * teta( ij ,l)
	613	tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
	614	ENDDO
	615	tpn = SSUM(iim,tppn,1)/apoln
	616	tps = SSUM(iim,tpps,1)/apols
	617
	618	DO ij = 1, iip1
	619	teta( ij ,l) = tpn
	620	teta(ij+ip1jm,l) = tps
	621	ENDDO
	622	ENDDO
	623
[49]	624	DO ij = 1,iim
	625	tppn(ij) = aire( ij ) * ps ( ij )
	626	tpps(ij) = aire(ij+ip1jm) * ps (ij+ip1jm)
	627	ENDDO
	628	tpn = SSUM(iim,tppn,1)/apoln
	629	tps = SSUM(iim,tpps,1)/apols
[2]	630
[49]	631	DO ij = 1, iip1
	632	ps( ij ) = tpn
	633	ps(ij+ip1jm) = tps
	634	ENDDO
	635
	636
[2]	637	END IF
[99]	638
	639	c ajout debug
[113]	640	c IF( lafin ) then
	641	c abort_message = 'Simulation finished'
	642	c call abort_gcm(modname,abort_message,0)
	643	c ENDIF
[2]	644
	645	c ********************************************************************
	646	c ********************************************************************
	647	c .... fin de l'integration dynamique et physique pour le pas itau ..
	648	c ********************************************************************
	649	c ********************************************************************
	650
	651	c preparation du pas d'integration suivant ......
	652
	653	IF ( .NOT.purmats ) THEN
	654	c ........................................................
	655	c .............. schema matsuno + leapfrog ..............
	656	c ........................................................
	657
	658	IF(forward. OR. leapf) THEN
	659	itau= itau + 1
	660	iday= day_ini+itau/day_step
	661	time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	662	IF(time.GT.1.) THEN
	663	time = time-1.
	664	iday = iday+1
	665	ENDIF
	666	ENDIF
	667
	668
	669	IF( itau. EQ. itaufinp1 ) then
[403]	670
[495]	671
[2]	672	abort_message = 'Simulation finished'
	673	call abort_gcm(modname,abort_message,0)
	674	ENDIF
	675	c-----------------------------------------------------------------------
	676	c ecriture du fichier histoire moyenne:
	677	c -------------------------------------
	678
	679	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	680	IF(itau.EQ.itaufin) THEN
	681	iav=1
	682	ELSE
	683	iav=0
	684	ENDIF
	685	CALL writedynav(histaveid, nqmx, itau,vcov ,
	686	, ucov,teta,pk,phi,q,masse,ps,phis)
[403]	687	cccIM cf. FH
	688	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
	689	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
	690
[2]	691	ENDIF
	692
	693	c-----------------------------------------------------------------------
	694	c ecriture de la bande histoire:
	695	c ------------------------------
	696
	697	IF( MOD(itau,iecri*day_step).EQ.0) THEN
	698
	699	nbetat = nbetatdem
	700	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
[295]	701	c CALL writehist( histid, histvid, nqmx, itau,vcov ,
	702	c , ucov,teta,phi,q,masse,ps,phis)
[2]	703
	704
	705	ENDIF
	706
	707	IF(itau.EQ.itaufin) THEN
	708
	709
	710	PRINT *,' Appel test_period avant redem ', itau
	711	CALL test_period ( ucov,vcov,teta,q,p,phis )
	712	CALL dynredem1("restart.nc",0.0,
	713	. vcov,ucov,teta,q,nqmx,masse,ps)
	714
	715	CLOSE(99)
	716	ENDIF
	717
	718	c-----------------------------------------------------------------------
	719	c gestion de l'integration temporelle:
	720	c ------------------------------------
	721
	722	IF( MOD(itau,iperiod).EQ.0 ) THEN
	723	GO TO 1
	724	ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
	725
	726	IF( forward ) THEN
	727	c fin du pas forward et debut du pas backward
	728
	729	forward = .FALSE.
	730	leapf = .FALSE.
	731	GO TO 2
	732
	733	ELSE
	734	c fin du pas backward et debut du premier pas leapfrog
	735
	736	leapf = .TRUE.
	737	dt = 2.*dtvr
	738	GO TO 2
	739	END IF
	740	ELSE
	741
	742	c ...... pas leapfrog .....
	743
	744	leapf = .TRUE.
	745	dt = 2.*dtvr
	746	GO TO 2
	747	END IF
	748
	749	ELSE
	750
	751	c ........................................................
	752	c .............. schema matsuno ...............
	753	c ........................................................
	754	IF( forward ) THEN
	755
	756	itau = itau + 1
	757	iday = day_ini+itau/day_step
	758	time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
	759
	760	IF(time.GT.1.) THEN
	761	time = time-1.
	762	iday = iday+1
	763	ENDIF
	764
	765	forward = .FALSE.
	766	IF( itau. EQ. itaufinp1 ) then
	767	abort_message = 'Simulation finished'
	768	call abort_gcm(modname,abort_message,0)
	769	ENDIF
	770	GO TO 2
	771
	772	ELSE
	773
	774	IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
	775	IF(itau.EQ.itaufin) THEN
	776	iav=1
	777	ELSE
	778	iav=0
	779	ENDIF
	780	CALL writedynav(histaveid, nqmx, itau,vcov ,
	781	, ucov,teta,pk,phi,q,masse,ps,phis)
[403]	782	cccIM cf. FH
	783	call bilan_dyn (2,dtvriperiod,dtvrday_step*periodav,
	784	, ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
[2]	785
	786	ENDIF
	787
	788	IF(MOD(itau,iecri*day_step).EQ.0) THEN
	789	nbetat = nbetatdem
	790	CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
[295]	791	c CALL writehist( histid, histvid, nqmx, itau,vcov ,
	792	c , ucov,teta,phi,q,masse,ps,phis)
[2]	793	ENDIF
	794
	795	IF(itau.EQ.itaufin)
	796	. CALL dynredem1("restart.nc",0.0,
	797	. vcov,ucov,teta,q,nqmx,masse,ps)
	798
	799	forward = .TRUE.
	800	GO TO 1
	801
	802	ENDIF
	803
	804	END IF
	805
	806	300 FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x,
	807	. 'c''est a dire du jour',i7,3x,'au jour',i7//)
	808
	809	STOP
	810	END

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

Download in other formats:

Original Format