Context Navigation

physiq.F @ 894

Last change on this file since 894 was 892, checked in by slebonnois, 12 years ago
SL: Important commit ! Adaptation of Venus physics to parallel computation / template for arch on the LMD servers using ifort / documentation for 1D column physics and for parallel computations
File size: 39.7 KB

Rev	Line
[3]	1	!
	2	! $Header: /home/cvsroot/LMDZ4/libf/phylmd/physiq.F,v 1.8 2005/02/24 09:58:18 fairhead Exp $
	3	!
	4	c
	5	SUBROUTINE physiq (nlon,nlev,nqmax,
	6	. debut,lafin,rjourvrai,gmtime,pdtphys,
[97]	7	. paprs,pplay,ppk,pphi,pphis,presnivs,
[3]	8	. u,v,t,qx,
	9	. omega,
	10	. d_u, d_v, d_t, d_qx, d_ps)
	11
	12	c======================================================================
	13	c
	14	c Modifications pour la physique de Venus
	15	c S. Lebonnois (LMD/CNRS) Septembre 2005
	16	c
	17	c ---------------------------------------------------------------------
	18	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	19	c
	20	c Objet: Moniteur general de la physique du modele
	21	cAA Modifications quant aux traceurs :
	22	cAA - uniformisation des parametrisations ds phytrac
	23	cAA - stockage des moyennes des champs necessaires
	24	cAA en mode traceur off-line
	25	c modif ( P. Le Van , 12/10/98 )
	26	c
	27	c Arguments:
	28	c
	29	c nlon----input-I-nombre de points horizontaux
	30	c nlev----input-I-nombre de couches verticales
	31	c nqmax---input-I-nombre de traceurs
	32	c debut---input-L-variable logique indiquant le premier passage
	33	c lafin---input-L-variable logique indiquant le dernier passage
	34	c rjour---input-R-numero du jour de l'experience
	35	c gmtime--input-R-temps universel dans la journee (0 a RDAY s)
	36	c pdtphys-input-R-pas d'integration pour la physique (seconde)
	37	c paprs---input-R-pression pour chaque inter-couche (en Pa)
	38	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	39	c ppk ---input-R-fonction d'Exner au milieu de couche
	40	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
	41	c pphis---input-R-geopotentiel du sol
	42	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
	43	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
	44	c v-------input-R-vitesse Y (de S a N) en m/s
	45	c t-------input-R-temperature (K)
	46	c qx------input-R-mass mixing ratio traceurs (kg/kg)
	47	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
	48	c omega---input-R-vitesse verticale en Pa/s
	49	c
	50	c d_u-----output-R-tendance physique de "u" (m/s/s)
	51	c d_v-----output-R-tendance physique de "v" (m/s/s)
	52	c d_t-----output-R-tendance physique de "t" (K/s)
	53	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
	54	c d_ps----output-R-tendance physique de la pression au sol
	55	c======================================================================
[101]	56	USE ioipsl
[849]	57	! USE histcom ! not needed; histcom is included in ioipsl
[101]	58	USE infotrac
	59	USE control_mod
	60	use dimphy
	61	USE comgeomphy
[892]	62	USE mod_phys_lmdz_para, only : is_parallel,jj_nb
	63	USE phys_state_var_mod ! Variables sauvegardees de la physique
	64	USE iophy
[101]	65	IMPLICIT none
	66	c======================================================================
	67	c CLEFS CPP POUR LES IO
	68	c =====================
	69	c#define histhf
	70	#define histday
	71	#define histmth
	72	#define histins
	73	c======================================================================
[3]	74	#include "dimensions.h"
	75	integer jjmp1
	76	parameter (jjmp1=jjm+1-1/jjm)
	77	#include "dimsoil.h"
	78	#include "clesphys.h"
	79	#include "temps.h"
	80	#include "iniprint.h"
	81	#include "timerad.h"
	82	#include "logic.h"
[815]	83	#include "tabcontrol.h"
[3]	84	c======================================================================
	85	LOGICAL ok_journe ! sortir le fichier journalier
	86	save ok_journe
	87	c PARAMETER (ok_journe=.true.)
	88	c
	89	LOGICAL ok_mensuel ! sortir le fichier mensuel
	90	save ok_mensuel
	91	c PARAMETER (ok_mensuel=.true.)
	92	c
	93	LOGICAL ok_instan ! sortir le fichier instantane
	94	save ok_instan
	95	c PARAMETER (ok_instan=.true.)
	96	c
	97	c======================================================================
	98	c
	99	c Variables argument:
	100	c
	101	INTEGER nlon
	102	INTEGER nlev
	103	INTEGER nqmax
	104	REAL rjourvrai
	105	REAL gmtime
	106	REAL pdtphys
	107	LOGICAL debut, lafin
	108	REAL paprs(klon,klev+1)
	109	REAL pplay(klon,klev)
	110	REAL pphi(klon,klev)
	111	REAL pphis(klon)
	112	REAL presnivs(klev)
	113
	114	! ADAPTATION GCM POUR CP(T)
	115	REAL ppk(klon,klev)
	116
	117	REAL u(klon,klev)
	118	REAL v(klon,klev)
	119	REAL t(klon,klev)
	120	REAL qx(klon,klev,nqmax)
	121
	122	REAL d_u_dyn(klon,klev)
	123	REAL d_t_dyn(klon,klev)
	124
	125	REAL omega(klon,klev)
	126
	127	REAL d_u(klon,klev)
	128	REAL d_v(klon,klev)
	129	REAL d_t(klon,klev)
	130	REAL d_qx(klon,klev,nqmax)
	131	REAL d_ps(klon)
	132
	133	logical ok_hf
	134	real ecrit_hf
	135	integer nid_hf
	136	save ok_hf, ecrit_hf, nid_hf
	137
	138	#ifdef histhf
	139	data ok_hf,ecrit_hf/.true.,0.25/
	140	#else
	141	data ok_hf/.false./
	142	#endif
	143
	144	c Variables propres a la physique
	145	c
[101]	146	REAL,save,allocatable :: rlev(:,:) ! altitude a chaque niveau (interface inferieure de la couche)
	147	INTEGER,save :: itap ! compteur pour la physique
[3]	148	REAL delp(klon,klev) ! epaisseur d'une couche
	149
	150	INTEGER igwd,idx(klon),itest(klon)
	151	c
	152	c Diagnostiques 2D de drag_noro, lift_noro et gw_nonoro
	153
	154	REAL zulow(klon),zvlow(klon)
	155	REAL zustrdr(klon), zvstrdr(klon)
	156	REAL zustrli(klon), zvstrli(klon)
	157	REAL zustrhi(klon), zvstrhi(klon)
	158
	159	c Pour calcul GW drag oro et nonoro: CALCUL de N2:
	160	real zdzlev(klon,klev)
	161	real ztlev(klon,klev),zpklev(klon,klev)
	162	real ztetalay(klon,klev),ztetalev(klon,klev)
	163	real zdtetalev(klon,klev)
	164	real zn2(klon,klev) ! BV^2 at plev
	165
	166	c Pour les bilans de moment angulaire,
	167	integer bilansmc
	168	c Pour le transport de ballons
	169	integer ballons
	170	c j'ai aussi besoin
	171	c du stress de couche limite a la surface:
	172
	173	REAL zustrcl(klon),zvstrcl(klon)
	174
	175	c et du stress total c de la physique:
	176
	177	REAL zustrph(klon),zvstrph(klon)
	178
	179	c Variables locales:
	180	c
	181	REAL cdragh(klon) ! drag coefficient pour T and Q
	182	REAL cdragm(klon) ! drag coefficient pour vent
	183	c
	184	cAA Pour TRACEURS
	185	cAA
[105]	186	REAL,save,allocatable :: source(:,:)
[3]	187	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
	188	REAL yu1(klon) ! vents dans la premiere couche U
	189	REAL yv1(klon) ! vents dans la premiere couche V
	190
	191	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
	192	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
	193	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
	194	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
	195	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
	196	c
	197
	198	c======================================================================
	199	c
	200	c Declaration des procedures appelees
	201	c
	202	EXTERNAL ajsec ! ajustement sec
	203	EXTERNAL clmain ! couche limite
	204	EXTERNAL hgardfou ! verifier les temperatures
	205	c EXTERNAL orbite ! calculer l'orbite
	206	EXTERNAL phyetat0 ! lire l'etat initial de la physique
	207	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
	208	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
	209	EXTERNAL suphec ! initialiser certaines constantes
	210	EXTERNAL transp ! transport total de l'eau et de l'energie
	211	EXTERNAL abort_gcm
	212	EXTERNAL printflag
	213	EXTERNAL zenang
	214	EXTERNAL diagetpq
	215	EXTERNAL conf_phys
	216	EXTERNAL diagphy
	217	EXTERNAL mucorr
	218	EXTERNAL phytrac
	219	c
	220	c Variables locales
	221	c
	222	CXXX PB
	223	REAL fluxt(klon,klev) ! flux turbulent de chaleur
	224	REAL fluxu(klon,klev) ! flux turbulent de vitesse u
	225	REAL fluxv(klon,klev) ! flux turbulent de vitesse v
	226	c
	227	REAL flux_dyn(klon,klev) ! flux de chaleur produit par la dynamique
	228	REAL flux_ajs(klon,klev) ! flux de chaleur ajustement sec
	229	REAL flux_ec(klon,klev) ! flux de chaleur Ec
	230	c
	231	REAL tmpout(klon,klev) ! K s-1
	232
	233	INTEGER itaprad
	234	SAVE itaprad
	235	c
	236	REAL dist, rmu0(klon), fract(klon)
	237	REAL zdtime, zlongi
	238	c
	239	INTEGER i, k, iq, ig, j, ll
	240	c
	241	REAL zphi(klon,klev)
	242	c
	243	c Variables du changement
	244	c
	245	c ajs: ajustement sec
	246	c vdf: couche limite (Vertical DiFfusion)
	247	REAL d_t_ajs(klon,klev), d_tr_ajs(klon,klev,nqmax)
	248	REAL d_u_ajs(klon,klev), d_v_ajs(klon,klev)
	249	c
	250	REAL d_ts(klon)
	251	c
	252	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
	253	REAL d_t_vdf(klon,klev), d_tr_vdf(klon,klev,nqmax)
	254	c
	255	CMOD LOTT: Tendances Orography Sous-maille
	256	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
	257	REAL d_t_oro(klon,klev)
	258	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
	259	REAL d_t_lif(klon,klev)
	260	C Tendances Ondes de G non oro (runs strato).
	261	REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
	262	REAL d_t_hin(klon,klev)
	263
	264	c
	265	c Variables liees a l'ecriture de la bande histoire physique
	266	c
	267	INTEGER ecrit_mth
	268	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
	269	c
	270	INTEGER ecrit_day
	271	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
	272	c
	273	INTEGER ecrit_ins
	274	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
	275	c
	276	integer itau_w ! pas de temps ecriture = itap + itau_phy
	277
	278	c Variables locales pour effectuer les appels en serie
	279	c
	280	REAL t_seri(klon,klev)
	281	REAL u_seri(klon,klev), v_seri(klon,klev)
	282	c
	283	REAL tr_seri(klon,klev,nqmax)
	284	REAL d_tr(klon,klev,nqmax)
	285	c
	286	c pour ioipsl
	287	INTEGER nid_day, nid_mth, nid_ins
	288	SAVE nid_day, nid_mth, nid_ins
	289	INTEGER nhori, nvert, idayref
[97]	290	REAL zsto, zout, zsto1, zsto2, zero
[3]	291	parameter (zero=0.0e0)
	292	real zjulian
	293	save zjulian
	294
	295	CHARACTER*2 str2
	296	character*20 modname
	297	character*80 abort_message
	298	logical ok_sync
	299
	300	character*30 nom_fichier
	301	character*10 varname
	302	character*40 vartitle
	303	character*20 varunits
	304	C Variables liees au bilan d'energie et d'enthalpi
	305	REAL ztsol(klon)
	306	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	307	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	308	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	309	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	310	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
	311	REAL d_h_vcol_phy
	312	REAL fs_bound, fq_bound
	313	SAVE d_h_vcol_phy
	314	REAL zero_v(klon),zero_v2(klon,klev)
	315	CHARACTER*15 ztit
	316	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
	317	SAVE ip_ebil
	318	DATA ip_ebil/2/
	319	INTEGER if_ebil ! level for energy conserv. dignostics
	320	SAVE if_ebil
	321	c+jld ec_conser
	322	REAL d_t_ec(klon,klev) ! tendance du a la conversion Ec -> E thermique
	323	c-jld ec_conser
	324
	325	c TEST VENUS...
	326	REAL mang(klon,klev) ! moment cinetique
	327	REAL mangtot ! moment cinetique total
	328
	329	c Declaration des constantes et des fonctions thermodynamiques
	330	c
	331	#include "YOMCST.h"
	332
	333	c======================================================================
	334	c INITIALISATIONS
	335	c================
	336
	337	modname = 'physiq'
	338	ok_sync=.TRUE.
	339
[119]	340	bilansmc = 0
[3]	341	ballons = 0
[892]	342	! NE FONCTIONNENT PAS ENCORE EN PARALLELE !!!
	343	if (is_parallel) then
	344	bilansmc = 0
	345	ballons = 0
	346	endif
[3]	347
	348	IF (if_ebil.ge.1) THEN
	349	DO i=1,klon
	350	zero_v(i)=0.
	351	END DO
	352	DO i=1,klon
	353	DO j=1,klev
	354	zero_v2(i,j)=0.
	355	END DO
	356	END DO
	357	END IF
	358
	359	c PREMIER APPEL SEULEMENT
	360	c========================
	361	IF (debut) THEN
[892]	362	allocate(rlev(klon,klevp1))
[105]	363	allocate(source(klon,nqmax))
[3]	364
	365	CALL suphec ! initialiser constantes et parametres phys.
	366
	367	IF (if_ebil.ge.1) d_h_vcol_phy=0.
	368	c
	369	c appel a la lecture du physiq.def
	370	c
	371	call conf_phys(ok_journe, ok_mensuel,
	372	. ok_instan,
	373	. if_ebil)
	374
[892]	375	call phys_state_var_init
[3]	376	c
	377	c Initialiser les compteurs:
	378	c
	379	itap = 0
	380	itaprad = 0
	381	c
	382	c Lecture startphy.nc :
	383	c
[892]	384	CALL phyetat0 ("startphy.nc")
[3]	385
[815]	386	c dtime est defini dans tabcontrol.h et lu dans startphy
[150]	387	c pdtphys est calcule a partir des nouvelles conditions:
	388	c Reinitialisation du pas de temps physique quand changement
	389	IF (ABS(dtime-pdtphys).GT.0.001) THEN
	390	WRITE(lunout,*) 'Pas physique a change',dtime,
	391	. pdtphys
	392	c abort_message='Pas physique n est pas correct '
	393	c call abort_gcm(modname,abort_message,1)
[152]	394	c----------------
	395	c pour initialiser convenablement le time_counter, il faut tenir compte
	396	c du changement de dtime en changeant itau_phy (point de depart)
	397	itau_phy = NINT(itau_phy*dtime/pdtphys)
	398	c----------------
[150]	399	dtime=pdtphys
	400	ENDIF
	401
[3]	402	radpas = NINT( RDAY/pdtphys/nbapp_rad)
	403
[815]	404	CALL printflag( ok_journe,ok_instan )
[3]	405	c
	406	c---------
	407	c FLOTT
	408	IF (ok_orodr) THEN
	409	DO i=1,klon
	410	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	411	ENDDO
	412	CALL SUGWD(klon,klev,paprs,pplay)
	413	DO i=1,klon
	414	zuthe(i)=0.
	415	zvthe(i)=0.
	416	if(zstd(i).gt.10.)then
	417	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	418	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	419	endif
	420	ENDDO
	421	ENDIF
	422
	423	if (bilansmc.eq.1) then
	424	C OUVERTURE D'UN FICHIER FORMATTE POUR STOCKER LES COMPOSANTES
	425	C DU BILAN DE MOMENT ANGULAIRE.
	426	open(27,file='aaam_bud.out',form='formatted')
	427	open(28,file='fields_2d.out',form='formatted')
	428	write(,)'Ouverture de aaam_bud.out (FL Vous parle)'
	429	write(,)'Ouverture de fields_2d.out (FL Vous parle)'
	430	endif !bilansmc
	431
	432	c--------------SLEBONNOIS
	433	C OUVERTURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
	434	C DES BALLONS
	435	if (ballons.eq.1) then
	436	open(30,file='ballons-lat.out',form='formatted')
	437	open(31,file='ballons-lon.out',form='formatted')
	438	open(32,file='ballons-u.out',form='formatted')
	439	open(33,file='ballons-v.out',form='formatted')
	440	open(34,file='ballons-alt.out',form='formatted')
	441	write(,)'Ouverture des ballons*.out'
	442	endif !ballons
	443	c-------------
	444
	445	c---------
	446	C TRACEURS
	447	C source dans couche limite
	448	source = 0.0 ! pas de source, pour l'instant
	449	C
	450	c---------
	451	c
	452	c Verifications:
	453	c
	454	IF (nlon .NE. klon) THEN
	455	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
	456	. klon
	457	abort_message='nlon et klon ne sont pas coherents'
	458	call abort_gcm(modname,abort_message,1)
	459	ENDIF
	460	IF (nlev .NE. klev) THEN
	461	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
	462	. klev
	463	abort_message='nlev et klev ne sont pas coherents'
	464	call abort_gcm(modname,abort_message,1)
	465	ENDIF
	466	c
[808]	467	IF (dtimeREAL(radpas).GT.(RDAY0.25).AND.cycle_diurne)
[3]	468	$ THEN
	469	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	470	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	471	abort_message='Nbre d appels au rayonnement insuffisant'
	472	call abort_gcm(modname,abort_message,1)
	473	ENDIF
	474	c
	475	WRITE(lunout,*)"Clef pour la convection seche, iflag_ajs=",
	476	. iflag_ajs
	477	c
[97]	478	ecrit_mth = NINT(RDAY/dtime*ecritphy) ! tous les ecritphy jours
[3]	479	IF (ok_mensuel) THEN
	480	WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
	481	. ecrit_mth
	482	ENDIF
[97]	483
[3]	484	ecrit_day = NINT(RDAY/dtime *1.0) ! tous les jours
	485	IF (ok_journe) THEN
	486	WRITE(lunout,*)'La frequence de sortie journaliere est de ',
	487	. ecrit_day
	488	ENDIF
[97]	489
[101]	490	ecrit_ins = NINT(RDAY/dtime*ecritphy) ! Fraction de jour reglable
[3]	491	IF (ok_instan) THEN
	492	WRITE(lunout,*)'La frequence de sortie instant. est de ',
	493	. ecrit_ins
	494	ENDIF
	495
	496	c Initialisation des sorties
	497	c===========================
	498
	499	#ifdef CPP_IOIPSL
	500
	501	#ifdef histhf
	502	#include "ini_histhf.h"
	503	#endif
	504
	505	#ifdef histday
	506	#include "ini_histday.h"
	507	#endif
	508
	509	#ifdef histmth
	510	#include "ini_histmth.h"
	511	#endif
	512
	513	#ifdef histins
	514	#include "ini_histins.h"
	515	#endif
	516
	517	#endif
	518
	519	c
	520	c Initialiser les valeurs de u pour calculs tendances
	521	c (pour T, c'est fait dans phyetat0)
	522	c
	523	DO k = 1, klev
	524	DO i = 1, klon
	525	u_ancien(i,k) = u(i,k)
	526	ENDDO
	527	ENDDO
	528
	529
	530	ENDIF ! debut
	531	c====================================================================
	532	c======================================================================
	533
	534	c Mettre a zero des variables de sortie (pour securite)
	535	c
	536	DO i = 1, klon
	537	d_ps(i) = 0.0
	538	ENDDO
	539	DO k = 1, klev
	540	DO i = 1, klon
	541	d_t(i,k) = 0.0
	542	d_u(i,k) = 0.0
	543	d_v(i,k) = 0.0
	544	ENDDO
	545	ENDDO
	546	DO iq = 1, nqmax
	547	DO k = 1, klev
	548	DO i = 1, klon
	549	d_qx(i,k,iq) = 0.0
	550	ENDDO
	551	ENDDO
	552	ENDDO
	553	c
	554	c Ne pas affecter les valeurs entrees de u, v, h, et q
	555	c
	556	DO k = 1, klev
	557	DO i = 1, klon
	558	t_seri(i,k) = t(i,k)
	559	u_seri(i,k) = u(i,k)
	560	v_seri(i,k) = v(i,k)
	561	ENDDO
	562	ENDDO
	563	DO iq = 1, nqmax
	564	DO k = 1, klev
	565	DO i = 1, klon
	566	tr_seri(i,k,iq) = qx(i,k,iq)
	567	ENDDO
	568	ENDDO
	569	ENDDO
	570	C
	571	DO i = 1, klon
	572	ztsol(i) = ftsol(i)
	573	ENDDO
	574	C
	575	IF (if_ebil.ge.1) THEN
	576	ztit='after dynamic'
	577	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
	578	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	579	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	580	C Comme les tendances de la physique sont ajoute dans la dynamique,
	581	C on devrait avoir que la variation d'entalpie par la dynamique
	582	C est egale a la variation de la physique au pas de temps precedent.
	583	C Donc la somme de ces 2 variations devrait etre nulle.
	584	call diagphy(airephy,ztit,ip_ebil
	585	e , zero_v, zero_v, zero_v, zero_v, zero_v
	586	e , zero_v, zero_v, zero_v, ztsol
	587	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
	588	s , fs_bound, fq_bound )
	589	END IF
	590
	591	c====================================================================
	592	c Diagnostiquer la tendance dynamique
	593	c
	594	IF (ancien_ok) THEN
	595	DO k = 1, klev
	596	DO i = 1, klon
	597	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
	598	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	599	ENDDO
	600	ENDDO
	601
	602	! ADAPTATION GCM POUR CP(T)
	603	do i=1,klon
	604	flux_dyn(i,1) = 0.0
	605	do j=2,klev
	606	flux_dyn(i,j) = flux_dyn(i,j-1)
	607	. +cpdet(t_seri(i,j-1))/RGd_t_dyn(i,j-1)(paprs(i,j-1)-paprs(i,j))
	608	enddo
	609	enddo
	610
	611	ELSE
	612	DO k = 1, klev
	613	DO i = 1, klon
	614	d_u_dyn(i,k) = 0.0
	615	d_t_dyn(i,k) = 0.0
	616	ENDDO
	617	ENDDO
	618	ancien_ok = .TRUE.
	619	ENDIF
	620	c====================================================================
	621	c
	622	c Ajouter le geopotentiel du sol:
	623	c
	624	DO k = 1, klev
	625	DO i = 1, klon
	626	zphi(i,k) = pphi(i,k) + pphis(i)
	627	ENDDO
	628	ENDDO
	629	c====================================================================
	630	c
	631	c Verifier les temperatures
	632	c
	633	CALL hgardfou(t_seri,ftsol,'debutphy')
	634	c====================================================================
	635	c
	636	c Incrementer le compteur de la physique
	637	c
	638	itap = itap + 1
	639
	640	c====================================================================
	641	c Orbite et eclairement
	642	c====================================================================
	643
	644	c Pour VENUS, on fixe l'obliquite a 0 et l'eccentricite a 0.
	645	c donc pas de variations de Ls, ni de dist.
	646	c La seule chose qui compte, c'est la rotation de la planete devant
	647	c le Soleil...
	648
	649	zlongi = 0.0
	650	dist = 0.72 ! en UA
	651
	652	c Si on veut remettre l'obliquite a 3 degres et/ou l'eccentricite
	653	c a sa valeur, et prendre en compte leur evolution,
	654	c il faudra refaire un orbite.F...
	655	c CALL orbite(zlongi,dist)
	656
	657	IF (cycle_diurne) THEN
[808]	658	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[3]	659	CALL zenang(zlongi,gmtime,zdtime,rlatd,rlond,rmu0,fract)
	660	ELSE
	661	call mucorr(klon,zlongi,rlatd,rmu0,fract)
	662	ENDIF
	663
	664	c====================================================================
	665	c Calcul des tendances traceurs
	666	c====================================================================
	667
	668	if (iflag_trac.eq.1) then
	669	call phytrac (
	670	I itap, gmtime,
	671	I debut,lafin,
	672	I nqmax,
	673	I nlon,nlev,dtime,
	674	I u,v,t,paprs,pplay,
	675	I rlatd,
	676	I rlond,presnivs,pphis,pphi,
	677	I falbe,
	678	O tr_seri)
	679	endif
	680
	681	c
	682	c====================================================================
	683	c Appeler la diffusion verticale (programme de couche limite)
	684	c====================================================================
	685
	686	c-------------------------------
	687	c VENUS TEST: on ne tient pas compte des calculs de clmain mais on force
	688	c l'equilibre radiatif du sol
	689	if (1.eq.0) then
	690	if (debut) then
	691	print*,"ATTENTION, CLMAIN SHUNTEE..."
	692	endif
	693
	694	DO i = 1, klon
	695	sens(i) = 0.0e0 ! flux de chaleur sensible au sol
	696	fder(i) = 0.0e0
	697	dlw(i) = 0.0e0
	698	ENDDO
	699
	700	c Incrementer la temperature du sol
	701	c
	702	DO i = 1, klon
	703	d_ts(i) = dtime * radsol(i)/22000. !valeur calculee par GCM pour I=200
	704	ftsol(i) = ftsol(i) + d_ts(i)
	705	do j=1,nsoilmx
	706	ftsoil(i,j)=ftsol(i)
	707	enddo
	708	ENDDO
	709
	710	c-------------------------------
	711	else
	712	c-------------------------------
	713
	714	fder = dlw
	715
	716	c print*,"radsol avant clmain=",radsol(klon/2)
	717	c print*,"solsw avant clmain=",solsw(klon/2)
	718	c print*,"sollw avant clmain=",sollw(klon/2)
	719
	720	! ADAPTATION GCM POUR CP(T)
[808]	721
[3]	722	CALL clmain(dtime,itap,
	723	e t_seri,u_seri,v_seri,
	724	e rmu0,
	725	e ftsol,
	726	$ ftsoil,
	727	$ paprs,pplay,ppk,radsol,falbe,
	728	e solsw, sollw, sollwdown, fder,
	729	e rlond, rlatd, cuphy, cvphy,
	730	e debut, lafin,
	731	s d_t_vdf,d_u_vdf,d_v_vdf,d_ts,
	732	s fluxt,fluxu,fluxv,cdragh,cdragm,
	733	s dsens,
	734	s ycoefh,yu1,yv1)
	735
	736	c print*,"radsol apres clmain=",radsol(klon/2)
	737	c print*,"solsw apres clmain=",solsw(klon/2)
	738	c print*,"sollw apres clmain=",sollw(klon/2)
	739
	740	CXXX Incrementation des flux
	741	DO i = 1, klon
	742	sens(i) = - fluxt(i,1) ! flux de chaleur sensible au sol
	743	fder(i) = dlw(i) + dsens(i)
	744	ENDDO
	745	CXXX
	746
	747	DO k = 1, klev
	748	DO i = 1, klon
	749	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
	750	d_t_vdf(i,k)= d_t_vdf(i,k)/dtime ! K/s
	751	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
	752	d_u_vdf(i,k)= d_u_vdf(i,k)/dtime ! (m/s)/s
	753	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
	754	d_v_vdf(i,k)= d_v_vdf(i,k)/dtime ! (m/s)/s
	755	ENDDO
	756	ENDDO
	757	c
	758	c print,"d_t_vdf1=",d_t_vdf(1,:)dtime
	759	c print,"d_t_vdf2=",d_t_vdf(klon/2,:)dtime
	760	c print,"d_t_vdf3=",d_t_vdf(klon,:)dtime
	761	c print,"d_u_vdf=",d_u_vdf(klon/2,:)dtime
	762	c print,"d_v_vdf=",d_v_vdf(klon/2,:)dtime
	763
	764	C TRACEURS
	765
	766	if (iflag_trac.eq.1) then
[101]	767	DO k = 1, klev
	768	DO i = 1, klon
	769	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	770	ENDDO
	771	ENDDO
[3]	772	DO iq=1, nqmax
	773	CALL cltrac(dtime,ycoefh,t_seri,
	774	s tr_seri(1,1,iq),source,
	775	e paprs, pplay,delp,
	776	s d_tr_vdf(1,1,iq))
	777	tr_seri(:,:,iq) = tr_seri(:,:,iq) + d_tr_vdf(:,:,iq)
	778	d_tr_vdf(:,:,iq)= d_tr_vdf(:,:,iq)/dtime ! /s
	779	ENDDO
	780	endif
	781
	782	IF (if_ebil.ge.2) THEN
	783	ztit='after clmain'
	784	CALL diagetpq(airephy,ztit,ip_ebil,2,1,dtime
	785	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	786	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	787	call diagphy(airephy,ztit,ip_ebil
	788	e , zero_v, zero_v, zero_v, zero_v, sens
	789	e , zero_v, zero_v, zero_v, ztsol
	790	e , d_h_vcol, d_qt, d_ec
	791	s , fs_bound, fq_bound )
	792	END IF
	793	C
	794	c
	795	c Incrementer la temperature du sol
	796	c
	797	c print*,'Tsol avant clmain:',ftsol(1)
	798	DO i = 1, klon
	799	ftsol(i) = ftsol(i) + d_ts(i)
	800	ENDDO
	801	c print*,'Tsol apres clmain:',ftsol(1)
	802
	803	c Calculer la derive du flux infrarouge
	804	c
	805	DO i = 1, klon
	806	dlw(i) = - 4.0RSIGMAftsol(i)**3
	807	ENDDO
	808
	809	c-------------------------------
	810	endif ! fin du VENUS TEST
	811
	812	c
	813	c Appeler l'ajustement sec
	814	c
	815	c===================================================================
	816	c Convection seche
	817	c===================================================================
	818	c
	819	d_t_ajs(:,:)=0.
	820	d_u_ajs(:,:)=0.
	821	d_v_ajs(:,:)=0.
	822	d_tr_ajs(:,:,:)=0.
	823	c
	824	IF(prt_level>9)WRITE(lunout,*)
	825	. 'AVANT LA CONVECTION SECHE , iflag_ajs='
	826	s ,iflag_ajs
	827
	828	if(iflag_ajs.eq.0) then
	829	c Rien
	830	c ====
	831	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
	832
	833	else if(iflag_ajs.eq.1) then
	834
	835	c Ajustement sec
	836	c ==============
	837	IF(prt_level>9)WRITE(lunout,*)'ajsec'
	838
	839	! ADAPTATION GCM POUR CP(T)
	840	CALL ajsec(paprs, pplay, ppk, t_seri, u_seri, v_seri, nqmax,
	841	. tr_seri, d_t_ajs, d_u_ajs, d_v_ajs, d_tr_ajs)
	842
	843	! ADAPTATION GCM POUR CP(T)
	844	do i=1,klon
	845	flux_ajs(i,1) = 0.0
	846	do j=2,klev
	847	flux_ajs(i,j) = flux_ajs(i,j-1)
	848	. + cpdet(t_seri(i,j-1))/RG*d_t_ajs(i,j-1)/dtime
	849	. *(paprs(i,j-1)-paprs(i,j))
	850	enddo
	851	enddo
	852
	853	t_seri(:,:) = t_seri(:,:) + d_t_ajs(:,:)
	854	d_t_ajs(:,:)= d_t_ajs(:,:)/dtime ! K/s
	855	u_seri(:,:) = u_seri(:,:) + d_u_ajs(:,:)
	856	d_u_ajs(:,:)= d_u_ajs(:,:)/dtime ! (m/s)/s
	857	v_seri(:,:) = v_seri(:,:) + d_v_ajs(:,:)
	858	d_v_ajs(:,:)= d_v_ajs(:,:)/dtime ! (m/s)/s
	859	if (iflag_trac.eq.1) then
	860	tr_seri(:,:,:) = tr_seri(:,:,:) + d_tr_ajs(:,:,:)
	861	d_tr_ajs(:,:,:)= d_tr_ajs(:,:,:)/dtime ! /s
	862	endif
	863
	864	c print,"d_t_ajs1=",d_t_ajs(1,:)dtime
	865	c print,"d_t_ajs2=",d_t_ajs(klon/2,:)dtime
	866	c print,"d_t_ajs3=",d_t_ajs(klon,:)dtime
	867	c print,"d_u_ajs=",d_u_ajs(klon/2,:)dtime
	868	c print,"d_v_ajs=",d_v_ajs(klon/2,:)dtime
	869
	870	endif
	871	c
	872	IF (if_ebil.ge.2) THEN
	873	ztit='after dry_adjust'
	874	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	875	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	876	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	877	call diagphy(airephy,ztit,ip_ebil
	878	e , zero_v, zero_v, zero_v, zero_v, sens
	879	e , zero_v, zero_v, zero_v, ztsol
	880	e , d_h_vcol, d_qt, d_ec
	881	s , fs_bound, fq_bound )
	882	END IF
	883
	884	c====================================================================
	885	c RAYONNEMENT
	886	c====================================================================
	887
	888	IF (MOD(itaprad,radpas).EQ.0) THEN
	889	c print*,'RAYONNEMENT ',
	890	c $ ' (itaprad=',itaprad,'/radpas=',radpas,')'
	891
[808]	892	dtimerad = dtime*REAL(radpas) ! pas de temps du rayonnement (s)
[3]	893	c PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
	894
	895	c print*,"radsol avant radlwsw=",radsol(klon/2)
	896	c print*,"solsw avant radlwsw=",solsw(klon/2)
	897	c print*,"sollw avant radlwsw=",sollw(klon/2)
	898	c print*,"avant radlwsw"
[808]	899
[3]	900	CALL radlwsw
	901	e (dist, rmu0, fract,
	902	e paprs, pplay,ftsol, t_seri,
	903	s heat,cool,radsol,
	904	s topsw,toplw,solsw,sollw,
	905	s sollwdown,
	906	s lwnet, swnet)
[808]	907
[3]	908	c print*,"apres radlwsw"
	909	c print*,"radsol apres radlwsw=",radsol(klon/2)
	910	c print*,"solsw apres radlwsw=",solsw(klon/2)
	911	c print*,"sollw apres radlwsw=",sollw(klon/2)
	912	itaprad = 0
	913	DO k = 1, klev
	914	DO i = 1, klon
	915	dtrad(i,k) = heat(i,k)-cool(i,k)
	916	dtrad(i,k) = dtrad(i,k)/RDAY !K/s
	917	ENDDO
	918	ENDDO
	919	c print,"dtrad1=",dtrad(1,:)dtime
	920	c print,"dtrad2=",dtrad(klon/2,:)dtime
	921	c print,"dtrad3=",dtrad(klon,:)dtime
	922
	923	ENDIF
	924	itaprad = itaprad + 1
	925	c====================================================================
	926	c
	927	c Ajouter la tendance des rayonnements (tous les pas)
	928	c
	929	DO k = 1, klev
	930	DO i = 1, klon
	931	t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime
	932	ENDDO
	933	ENDDO
	934
	935	IF (if_ebil.ge.2) THEN
	936	ztit='after rad'
	937	CALL diagetpq(airephy,ztit,ip_ebil,2,2,dtime
	938	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	939	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	940	call diagphy(airephy,ztit,ip_ebil
	941	e , topsw, toplw, solsw, sollw, zero_v
	942	e , zero_v, zero_v, zero_v, ztsol
	943	e , d_h_vcol, d_qt, d_ec
	944	s , fs_bound, fq_bound )
	945	END IF
	946	c
	947
	948	c====================================================================
	949	c Calcul des gravity waves FLOTT
	950	c====================================================================
	951	c
	952	if (ok_orodr.or.ok_gw_nonoro) then
	953	c CALCUL DE N2
	954	do i=1,klon
	955	do k=2,klev
	956	ztlev(i,k) = (t_seri(i,k)+t_seri(i,k-1))/2.
	957	zpklev(i,k) = sqrt(ppk(i,k)*ppk(i,k-1))
	958	enddo
	959	enddo
	960	call t2tpot(klon*klev,ztlev, ztetalev,zpklev)
	961	call t2tpot(klon*klev,t_seri,ztetalay,ppk)
	962	do i=1,klon
	963	do k=2,klev
	964	zdtetalev(i,k) = ztetalay(i,k)-ztetalay(i,k-1)
	965	zdzlev(i,k) = (zphi(i,k)-zphi(i,k-1))/RG
	966	zn2(i,k) = RGzdtetalev(i,k)/(ztetalev(i,k)zdzlev(i,k))
	967	zn2(i,k) = max(zn2(i,k),1.e-12) ! securite
	968	enddo
[808]	969	zn2(i,1) = 1.e-12 ! securite
[3]	970	enddo
	971
	972	endif
	973
	974	c ----------------------------ORODRAG
	975	IF (ok_orodr) THEN
	976	c
	977	c selection des points pour lesquels le shema est actif:
	978	igwd=0
	979	DO i=1,klon
	980	itest(i)=0
	981	c IF ((zstd(i).gt.10.0)) THEN
	982	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	983	itest(i)=1
	984	igwd=igwd+1
	985	idx(igwd)=i
	986	ENDIF
	987	ENDDO
	988	c igwdim=MAX(1,igwd)
	989	c
	990	c A ADAPTER POUR VENUS!!!
	991	CALL drag_noro(klon,klev,dtime,paprs,pplay,zphi,zn2,
	992	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
	993	e igwd,idx,itest,
	994	e t_seri, u_seri, v_seri,
	995	s zulow, zvlow, zustrdr, zvstrdr,
	996	s d_t_oro, d_u_oro, d_v_oro)
	997
[808]	998	c print*,"d_u_oro=",d_u_oro(klon/2,:)
[3]	999	c ajout des tendances
	1000	t_seri(:,:) = t_seri(:,:) + d_t_oro(:,:)
	1001	d_t_oro(:,:)= d_t_oro(:,:)/dtime ! K/s
	1002	u_seri(:,:) = u_seri(:,:) + d_u_oro(:,:)
	1003	d_u_oro(:,:)= d_u_oro(:,:)/dtime ! (m/s)/s
	1004	v_seri(:,:) = v_seri(:,:) + d_v_oro(:,:)
	1005	d_v_oro(:,:)= d_v_oro(:,:)/dtime ! (m/s)/s
	1006	c
	1007	ELSE
	1008	d_t_oro = 0.
	1009	d_u_oro = 0.
	1010	d_v_oro = 0.
	1011	zustrdr = 0.
	1012	zvstrdr = 0.
	1013	c
	1014	ENDIF ! fin de test sur ok_orodr
	1015	c
	1016	c ----------------------------OROLIFT
	1017	IF (ok_orolf) THEN
[808]	1018	print*,"ok_orolf NOT IMPLEMENTED !"
	1019	stop
[3]	1020	c
	1021	c selection des points pour lesquels le shema est actif:
	1022	igwd=0
	1023	DO i=1,klon
	1024	itest(i)=0
	1025	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	1026	itest(i)=1
	1027	igwd=igwd+1
	1028	idx(igwd)=i
	1029	ENDIF
	1030	ENDDO
	1031	c igwdim=MAX(1,igwd)
	1032	c
	1033	c A ADAPTER POUR VENUS!!!
	1034	c CALL lift_noro(klon,klev,dtime,paprs,pplay,
	1035	c e rlatd,zmea,zstd,zpic,zgam,zthe,zpic,zval,
	1036	c e igwd,idx,itest,
	1037	c e t_seri, u_seri, v_seri,
	1038	c s zulow, zvlow, zustrli, zvstrli,
	1039	c s d_t_lif, d_u_lif, d_v_lif )
	1040
	1041	c
	1042	c ajout des tendances
	1043	t_seri(:,:) = t_seri(:,:) + d_t_lif(:,:)
	1044	d_t_lif(:,:)= d_t_lif(:,:)/dtime ! K/s
	1045	u_seri(:,:) = u_seri(:,:) + d_u_lif(:,:)
	1046	d_u_lif(:,:)= d_u_lif(:,:)/dtime ! (m/s)/s
	1047	v_seri(:,:) = v_seri(:,:) + d_v_lif(:,:)
	1048	d_v_lif(:,:)= d_v_lif(:,:)/dtime ! (m/s)/s
	1049	c
	1050	ELSE
	1051	d_t_lif = 0.
	1052	d_u_lif = 0.
	1053	d_v_lif = 0.
	1054	zustrli = 0.
	1055	zvstrli = 0.
	1056	c
	1057	ENDIF ! fin de test sur ok_orolf
	1058
	1059	c ---------------------------- NON-ORO GRAVITY WAVES
	1060	IF(ok_gw_nonoro) then
	1061
	1062	call flott_gwd_ran(klon,klev,dtime,pplay,zn2,
	1063	e t_seri, u_seri, v_seri,
	1064	o zustrhi,zvstrhi,
	1065	o d_t_hin, d_u_hin, d_v_hin)
	1066
	1067	c ajout des tendances
	1068
	1069	t_seri(:,:) = t_seri(:,:) + d_t_hin(:,:)
	1070	d_t_hin(:,:)= d_t_hin(:,:)/dtime ! K/s
	1071	u_seri(:,:) = u_seri(:,:) + d_u_hin(:,:)
	1072	d_u_hin(:,:)= d_u_hin(:,:)/dtime ! (m/s)/s
	1073	v_seri(:,:) = v_seri(:,:) + d_v_hin(:,:)
	1074	d_v_hin(:,:)= d_v_hin(:,:)/dtime ! (m/s)/s
	1075
	1076	ELSE
	1077	d_t_hin = 0.
	1078	d_u_hin = 0.
	1079	d_v_hin = 0.
	1080	zustrhi = 0.
	1081	zvstrhi = 0.
	1082
	1083	ENDIF ! fin de test sur ok_gw_nonoro
	1084
	1085	c====================================================================
	1086	c Transport de ballons
	1087	c====================================================================
	1088	if (ballons.eq.1) then
	1089	CALL ballon(30,pdtphys,rjourvrai,gmtime,rlatd,rlond,
	1090	c C t,pplay,u,v,pphi) ! alt above surface (smoothed for GCM)
	1091	C t,pplay,u,v,zphi) ! alt above planet average radius
	1092	endif !ballons
	1093
	1094	c====================================================================
	1095	c Bilan de mmt angulaire
	1096	c====================================================================
	1097	if (bilansmc.eq.1) then
	1098	CMODDEB FLOTT
	1099	C CALCULER LE BILAN DE MOMENT ANGULAIRE (DIAGNOSTIQUE)
	1100	C STRESS NECESSAIRES: COUCHE LIMITE ET TOUTE LA PHYSIQUE
	1101
	1102	DO i = 1, klon
	1103	zustrph(i)=0.
	1104	zvstrph(i)=0.
	1105	zustrcl(i)=0.
	1106	zvstrcl(i)=0.
	1107	ENDDO
	1108	DO k = 1, klev
	1109	DO i = 1, klon
	1110	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
	1111	c (paprs(i,k)-paprs(i,k+1))/rg
	1112	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
	1113	c (paprs(i,k)-paprs(i,k+1))/rg
	1114	zustrcl(i)=zustrcl(i)+d_u_vdf(i,k)*
	1115	c (paprs(i,k)-paprs(i,k+1))/rg
	1116	zvstrcl(i)=zvstrcl(i)+d_v_vdf(i,k)*
	1117	c (paprs(i,k)-paprs(i,k+1))/rg
	1118	ENDDO
	1119	ENDDO
	1120
	1121	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
	1122	C ra,rg,romega,
	1123	C rlatd,rlond,pphis,
	1124	C zustrdr,zustrli,zustrcl,
	1125	C zvstrdr,zvstrli,zvstrcl,
	1126	C paprs,u,v)
	1127
	1128	CCMODFIN FLOTT
	1129	endif !bilansmc
	1130
	1131	c====================================================================
	1132	c====================================================================
	1133	c Calculer le transport de l'eau et de l'energie (diagnostique)
	1134	c
	1135	c A REVOIR POUR VENUS...
	1136	c
	1137	c CALL transp (paprs,ftsol,
	1138	c e t_seri, q_seri, u_seri, v_seri, zphi,
	1139	c s ve, vq, ue, uq)
	1140	c
	1141	c====================================================================
	1142	c+jld ec_conser
	1143	DO k = 1, klev
	1144	DO i = 1, klon
	1145	d_t_ec(i,k)=0.5/cpdet(t_seri(i,k))
	1146	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
	1147	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
	1148	d_t_ec(i,k) = d_t_ec(i,k)/dtime
	1149	END DO
	1150	END DO
	1151	do i=1,klon
	1152	flux_ec(i,1) = 0.0
	1153	do j=2,klev
	1154	flux_ec(i,j) = flux_ec(i,j-1)
	1155	. +cpdet(t_seri(i,j-1))/RGd_t_ec(i,j-1)(paprs(i,j-1)-paprs(i,j))
	1156	enddo
	1157	enddo
	1158
	1159	c-jld ec_conser
	1160	c====================================================================
	1161	IF (if_ebil.ge.1) THEN
	1162	ztit='after physic'
	1163	CALL diagetpq(airephy,ztit,ip_ebil,1,1,dtime
	1164	e , t_seri,zero_v2,zero_v2,zero_v2,u_seri,v_seri,paprs,pplay
	1165	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1166	C Comme les tendances de la physique sont ajoute dans la dynamique,
	1167	C on devrait avoir que la variation d'entalpie par la dynamique
	1168	C est egale a la variation de la physique au pas de temps precedent.
	1169	C Donc la somme de ces 2 variations devrait etre nulle.
	1170	call diagphy(airephy,ztit,ip_ebil
	1171	e , topsw, toplw, solsw, sollw, sens
	1172	e , zero_v, zero_v, zero_v, ztsol
	1173	e , d_h_vcol, d_qt, d_ec
	1174	s , fs_bound, fq_bound )
	1175	C
	1176	d_h_vcol_phy=d_h_vcol
	1177	C
	1178	END IF
	1179	C
	1180	c=======================================================================
	1181	c SORTIES
	1182	c=======================================================================
	1183
	1184	c Convertir les incrementations en tendances
	1185	c
	1186	DO k = 1, klev
	1187	DO i = 1, klon
	1188	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	1189	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	1190	d_t(i,k) = ( t_seri(i,k) - t(i,k) ) / dtime
	1191	ENDDO
	1192	ENDDO
	1193	c
	1194	DO iq = 1, nqmax
	1195	DO k = 1, klev
	1196	DO i = 1, klon
	1197	d_qx(i,k,iq) = ( tr_seri(i,k,iq) - qx(i,k,iq) ) / dtime
	1198	ENDDO
	1199	ENDDO
	1200	ENDDO
	1201
	1202	c------------------------
	1203	c Calcul moment cinetique
	1204	c------------------------
	1205	c TEST VENUS...
	1206	c mangtot = 0.0
	1207	c DO k = 1, klev
	1208	c DO i = 1, klon
	1209	c mang(i,k) = RAcos(rlatd(i)RPI/180.)
	1210	c . (u_seri(i,k)+RAcos(rlatd(i)RPI/180.)ROMEGA)
	1211	c . airephy(i)(paprs(i,k)-paprs(i,k+1))/RG
	1212	c mangtot=mangtot+mang(i,k)
	1213	c ENDDO
	1214	c ENDDO
	1215	c print*,"Moment cinetique total = ",mangtot
	1216	c
	1217	c------------------------
	1218	c
	1219	c Sauvegarder les valeurs de t et u a la fin de la physique:
	1220	c
	1221	DO k = 1, klev
	1222	DO i = 1, klon
	1223	u_ancien(i,k) = u_seri(i,k)
	1224	t_ancien(i,k) = t_seri(i,k)
	1225	ENDDO
	1226	ENDDO
	1227	c
	1228	c=============================================================
	1229	c Ecriture des sorties
	1230	c=============================================================
	1231
	1232	#ifdef CPP_IOIPSL
	1233
	1234	#ifdef histhf
	1235	#include "write_histhf.h"
	1236	#endif
	1237
	1238	#ifdef histday
	1239	#include "write_histday.h"
	1240	#endif
	1241
	1242	#ifdef histmth
	1243	#include "write_histmth.h"
	1244	#endif
	1245
	1246	#ifdef histins
	1247	#include "write_histins.h"
	1248	#endif
	1249
	1250	#endif
	1251
	1252	c====================================================================
	1253	c Si c'est la fin, il faut conserver l'etat de redemarrage
	1254	c====================================================================
	1255	c
	1256	IF (lafin) THEN
	1257	itau_phy = itau_phy + itap
[892]	1258	CALL phyredem ("restartphy.nc")
[3]	1259
	1260	c--------------FLOTT
	1261	CMODEB LOTT
	1262	C FERMETURE DU FICHIER FORMATTE CONTENANT LES COMPOSANTES
	1263	C DU BILAN DE MOMENT ANGULAIRE.
	1264	if (bilansmc.eq.1) then
	1265	write(,)'Fermeture de aaam_bud.out (FL Vous parle)'
	1266	close(27)
	1267	close(28)
	1268	endif !bilansmc
	1269	CMODFIN
	1270	c-------------
	1271	c--------------SLEBONNOIS
	1272	C FERMETURE DES FICHIERS FORMATTES CONTENANT LES POSITIONS ET VITESSES
	1273	C DES BALLONS
	1274	if (ballons.eq.1) then
	1275	write(,)'Fermeture des ballons*.out'
	1276	close(30)
	1277	close(31)
	1278	close(32)
	1279	close(33)
	1280	close(34)
	1281	endif !ballons
	1282	c-------------
	1283	ENDIF
	1284
	1285	RETURN
	1286	END
	1287
	1288
	1289
	1290	***********************************************************************
	1291	***********************************************************************
	1292	***********************************************************************
	1293	***********************************************************************
	1294	***********************************************************************
	1295	***********************************************************************
	1296	***********************************************************************
	1297	***********************************************************************
	1298
	1299	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	1300	IMPLICIT none
	1301	c
	1302	c Tranformer une variable de la grille physique a
	1303	c la grille d'ecriture
	1304	c
	1305	INTEGER nfield,nlon,iim,jjmp1, jjm
	1306	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
	1307	c
	1308	INTEGER i, n, ig
	1309	c
	1310	jjm = jjmp1 - 1
	1311	DO n = 1, nfield
	1312	DO i=1,iim
	1313	ecrit(i,n) = fi(1,n)
	1314	ecrit(i+jjm*iim,n) = fi(nlon,n)
	1315	ENDDO
	1316	DO ig = 1, nlon - 2
	1317	ecrit(iim+ig,n) = fi(1+ig,n)
	1318	ENDDO
	1319	ENDDO
	1320	RETURN
	1321	END

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

Download in other formats:

Original Format