Context Navigation

physiq.F @ 1301

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

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

Download in other formats:

Original Format