Context Navigation

leapfrog.F @ 131

Last change on this file since 131 was 130, checked in by slebonnois, 14 years ago

Sebastien:

correction de bugs dans phytitan suite a compil avec -debug sur gnome
elimination de nbetat* dans gcm.F et leapfrog.F (et dans le parallele)
correction rday -> rjour dans sortvarc.F
ajustement des arch-GNOME*

File size: 27.3 KB

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

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

Download in other formats:

Original Format