Context Navigation

physiq.F @ 956

Last change on this file since 956 was 956, checked in by lmdzadmin, 17 years ago

Nettoyage du controle des parametres physiques. FH

Les parametres cycle_diurne, soil_model, new_oliq, ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad et iflag_con
sont maintenant geres par la physique uniquement.
ecritphy est elimine.
dimphy.F90 et clesphys.h ne sont plus utilises par le code dynamique.
Le test academique obtenu en compilant avec
makegcm -p nophys gcm
fonctionne. FH
IM

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 110.1 KB

Rev	Line
[524]	1	!
	2	! $Header$
	3	!
	4	c
[766]	5	c#define IO_DEBUG
	6
[524]	7	SUBROUTINE physiq (nlon,nlev,nqmax,
	8	. debut,lafin,rjourvrai,gmtime,pdtphys,
	9	. paprs,pplay,pphi,pphis,presnivs,clesphy0,
	10	. u,v,t,qx,
	11	. omega,
[616]	12	#ifdef INCA
[524]	13	. flxmass_w,
	14	#endif
[644]	15	. d_u, d_v, d_t, d_qx, d_ps
	16	. , dudyn
	17	. , PVteta)
[524]	18
	19	USE ioipsl
[766]	20	USE comgeomphy
	21	USE write_field_phy
	22	USE dimphy
[776]	23	USE mod_grid_phy_lmdz
	24	USE mod_phys_lmdz_para
[766]	25	USE iophy
	26	USE misc_mod, mydebug=>debug
	27	USE vampir
[782]	28	USE pbl_surface_mod, ONLY : pbl_surface
[904]	29	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
[913]	30	USE phys_state_var_mod ! Variables sauvegardees de la physique
[782]	31
	32
	33	USE ocean_slab_mod, ONLY : ocean_slab_get_vars
	34	USE ocean_cpl_mod, ONLY : ocean_cpl_get_vars
	35	USE ocean_forced_mod, ONLY : ocean_forced_get_vars
	36	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
[909]	37	USE phys_output_mod
[782]	38
[524]	39	IMPLICIT none
	40	c======================================================================
	41	c
	42	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	43	c
	44	c Objet: Moniteur general de la physique du modele
	45	cAA Modifications quant aux traceurs :
	46	cAA - uniformisation des parametrisations ds phytrac
	47	cAA - stockage des moyennes des champs necessaires
	48	cAA en mode traceur off-line
	49	c======================================================================
	50	c CLEFS CPP POUR LES IO
	51	c =====================
[766]	52	c#define histhf
[800]	53	#define histday
[524]	54	#define histmth
[933]	55	c#define histmthNMC
[766]	56	c#define histins
	57	c#define histISCCP
[524]	58	c======================================================================
	59	c modif ( P. Le Van , 12/10/98 )
	60	c
	61	c Arguments:
	62	c
	63	c nlon----input-I-nombre de points horizontaux
	64	c nlev----input-I-nombre de couches verticales
	65	c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
	66	c debut---input-L-variable logique indiquant le premier passage
	67	c lafin---input-L-variable logique indiquant le dernier passage
	68	c rjour---input-R-numero du jour de l'experience
	69	c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
	70	c pdtphys-input-R-pas d'integration pour la physique (seconde)
	71	c paprs---input-R-pression pour chaque inter-couche (en Pa)
	72	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	73	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
	74	c pphis---input-R-geopotentiel du sol
	75	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
	76	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
	77	c v-------input-R-vitesse Y (de S a N) en m/s
	78	c t-------input-R-temperature (K)
	79	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
	80	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
	81	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
	82	c omega---input-R-vitesse verticale en Pa/s
	83	c d_u-----output-R-tendance physique de "u" (m/s/s)
	84	c d_v-----output-R-tendance physique de "v" (m/s/s)
	85	c d_t-----output-R-tendance physique de "t" (K/s)
	86	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
	87	c d_ps----output-R-tendance physique de la pression au sol
[776]	88	cIM
[644]	89	c PVteta--output-R-vorticite potentielle a des thetas constantes
[524]	90	c======================================================================
	91	#include "dimensions.h"
	92	integer jjmp1
	93	parameter (jjmp1=jjm+1-1/jjm)
[766]	94	integer iip1
	95	parameter (iip1=iim+1)
[782]	96
[524]	97	#include "regdim.h"
	98	#include "indicesol.h"
	99	#include "dimsoil.h"
	100	#include "clesphys.h"
	101	#include "control.h"
[616]	102	#include "logic.h"
[524]	103	#include "temps.h"
[766]	104	cym#include "comgeomphy.h"
[524]	105	#include "advtrac.h"
	106	#include "iniprint.h"
[541]	107	#include "thermcell.h"
[524]	108	c======================================================================
	109	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
	110	PARAMETER (ok_cvl=.TRUE.)
	111	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
	112	PARAMETER (ok_gust=.FALSE.)
[879]	113	integer iflag_radia ! active ou non le rayonnement (MPL)
	114	save iflag_radia
[524]	115	c======================================================================
	116	LOGICAL check ! Verifier la conservation du modele en eau
	117	PARAMETER (check=.FALSE.)
	118	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
	119	PARAMETER (ok_stratus=.FALSE.)
	120	c======================================================================
	121	logical rnpb
	122	#ifdef INCA
	123	parameter(rnpb=.false.)
	124	#else
	125	parameter(rnpb=.true.)
	126	#endif
[782]	127
[644]	128	cIM "slab" ocean
	129	REAL tslab(klon) !Temperature du slab-ocean
	130	REAL seaice(klon) !glace de mer (kg/m2)
	131	REAL fluxo(klon) !flux turbulents ocean-glace de mer
	132	REAL fluxg(klon) !flux turbulents ocean-atmosphere
[687]	133	REAL amn, amx
[879]	134	INTEGER igout
[524]	135	c======================================================================
	136	c Clef controlant l'activation du cycle diurne:
	137	ccc LOGICAL cycle_diurne
	138	ccc PARAMETER (cycle_diurne=.FALSE.)
	139	c======================================================================
	140	c Modele thermique du sol, a activer pour le cycle diurne:
	141	ccc LOGICAL soil_model
	142	ccc PARAMETER (soil_model=.FALSE.)
	143	c======================================================================
	144	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
	145	c le calcul du rayonnement est celle apres la precipitation des nuages.
	146	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
	147	c la condensation et la precipitation. Cette cle augmente les impacts
	148	c radiatifs des nuages.
	149	ccc LOGICAL new_oliq
	150	ccc PARAMETER (new_oliq=.FALSE.)
	151	c======================================================================
	152	c Clefs controlant deux parametrisations de l'orographie:
	153	cc LOGICAL ok_orodr
	154	ccc PARAMETER (ok_orodr=.FALSE.)
	155	ccc LOGICAL ok_orolf
	156	ccc PARAMETER (ok_orolf=.FALSE.)
	157	c======================================================================
	158	LOGICAL ok_journe ! sortir le fichier journalier
	159	save ok_journe
[766]	160	c$OMP THREADPRIVATE(ok_journe)
[524]	161	c
	162	LOGICAL ok_mensuel ! sortir le fichier mensuel
	163	save ok_mensuel
[766]	164	c$OMP THREADPRIVATE(ok_mensuel)
[524]	165	c
	166	LOGICAL ok_instan ! sortir le fichier instantane
	167	save ok_instan
[766]	168	c$OMP THREADPRIVATE(ok_instan)
[524]	169	c
	170	LOGICAL ok_region ! sortir le fichier regional
	171	PARAMETER (ok_region=.FALSE.)
	172	c======================================================================
[878]	173	real weak_inversion(klon),dthmin(klon)
	174	real seuil_inversion
	175	save seuil_inversion
	176	c$OMP THREADPRIVATE(seuil_inversion)
	177	integer iflag_ratqs
	178	save iflag_ratqs
	179	c$OMP THREADPRIVATE(iflag_ratqs)
	180
	181	integer lmax_th(klon)
	182	integer limbas(klon)
	183	real ratqscth(klon,klev)
	184	real ratqsdiff(klon,klev)
	185	real zqsatth(klon,klev)
	186
[541]	187	c======================================================================
[524]	188	c
	189	INTEGER ivap ! indice de traceurs pour vapeur d'eau
	190	PARAMETER (ivap=1)
	191	INTEGER iliq ! indice de traceurs pour eau liquide
	192	PARAMETER (iliq=2)
	193
	194	c
	195	c
	196	c Variables argument:
	197	c
	198	INTEGER nlon
	199	INTEGER nlev
	200	INTEGER nqmax
	201	REAL rjourvrai
	202	REAL gmtime
	203	REAL pdtphys
	204	LOGICAL debut, lafin
	205	REAL paprs(klon,klev+1)
	206	REAL pplay(klon,klev)
	207	REAL pphi(klon,klev)
	208	REAL pphis(klon)
	209	REAL presnivs(klev)
	210	REAL znivsig(klev)
[644]	211	real pir
[719]	212
[524]	213	REAL u(klon,klev)
	214	REAL v(klon,klev)
[879]	215	REAL t(klon,klev),theta(klon,klev)
[524]	216	REAL qx(klon,klev,nqmax)
	217
	218	REAL omega(klon,klev)
	219
[616]	220	#ifdef INCA
[524]	221	REAL flxmass_w(klon,klev)
	222	#endif
	223	REAL d_u(klon,klev)
	224	REAL d_v(klon,klev)
	225	REAL d_t(klon,klev)
	226	REAL d_qx(klon,klev,nqmax)
	227	REAL d_ps(klon)
[619]	228	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
[644]	229	c
	230	cIM Amip2 PV a theta constante
	231	c
	232	INTEGER nbteta
	233	PARAMETER(nbteta=3)
	234	CHARACTER*3 ctetaSTD(nbteta)
	235	DATA ctetaSTD/'350','380','405'/
[766]	236	SAVE ctetaSTD
	237	c$OMP THREADPRIVATE(ctetaSTD)
[644]	238	REAL rtetaSTD(nbteta)
	239	DATA rtetaSTD/350., 380., 405./
[766]	240	SAVE rtetaSTD
	241	c$OMP THREADPRIVATE(rtetaSTD)
[644]	242	c
	243	REAL PVteta(klon,nbteta)
	244	REAL zx_tmp_3dte(iim,jjmp1,nbteta)
	245	c
	246	cMI Amip2 PV a theta constante
[524]	247
[766]	248	cym INTEGER klevp1, klevm1
	249	cym PARAMETER(klevp1=klev+1,klevm1=klev-1)
	250	cym#include "raddim.h"
[524]	251	c
	252	c
[644]	253	cIM Amip2
	254	c variables a une pression donnee
[524]	255	c
	256	real rlevSTD(nlevSTD)
	257	DATA rlevSTD/100000., 92500., 85000., 70000.,
	258	.60000., 50000., 40000., 30000., 25000., 20000.,
	259	.15000., 10000., 7000., 5000., 3000., 2000., 1000./
[766]	260	SAVE rlevstd
[644]	261	CHARACTER*4 clevSTD(nlevSTD)
[524]	262	DATA clevSTD/'1000','925 ','850 ','700 ','600 ',
	263	.'500 ','400 ','300 ','250 ','200 ','150 ','100 ',
	264	.'70 ','50 ','30 ','20 ','10 '/
[766]	265	SAVE clevSTD
[524]	266	c
[901]	267	CHARACTER*4 bb2
[644]	268	CHARACTER*2 bb3
	269	c
[524]	270	real tlevSTD(klon,nlevSTD), qlevSTD(klon,nlevSTD)
	271	real rhlevSTD(klon,nlevSTD), philevSTD(klon,nlevSTD)
	272	real ulevSTD(klon,nlevSTD), vlevSTD(klon,nlevSTD)
[644]	273	real wlevSTD(klon,nlevSTD)
[524]	274	c
[644]	275	c nout : niveau de output des variables a une pression donnee
	276	logical oknondef(klon,nlevSTD,nout)
	277	c
	278	c les produits uvSTD, vqSTD, .., T2STD sont calcules
	279	c a partir des valeurs instantannees toutes les 6 h
	280	c qui sont moyennees sur le mois
	281	c
	282	real uvSTD(klon,nlevSTD)
	283	real vqSTD(klon,nlevSTD)
	284	real vTSTD(klon,nlevSTD)
	285	real wqSTD(klon,nlevSTD)
	286	c
	287	real vphiSTD(klon,nlevSTD)
	288	real wTSTD(klon,nlevSTD)
	289	real u2STD(klon,nlevSTD)
	290	real v2STD(klon,nlevSTD)
	291	real T2STD(klon,nlevSTD)
	292	c
	293	#include "radepsi.h"
	294	#include "radopt.h"
	295	c
	296	c
[524]	297	c prw: precipitable water
	298	real prw(klon)
	299
	300	REAL convliq(klon,klev) ! eau liquide nuageuse convective
	301	REAL convfra(klon,klev) ! fraction nuageuse convective
	302
	303	REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
	304	REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
	305	REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
	306	REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
	307
[766]	308	INTEGER linv, kp1
[524]	309	c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
	310	c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
	311	REAL flwp(klon), fiwp(klon)
	312	REAL flwc(klon,klev), fiwc(klon,klev)
	313	REAL flwp_c(klon), fiwp_c(klon)
	314	REAL flwc_c(klon,klev), fiwc_c(klon,klev)
	315	REAL flwp_s(klon), fiwp_s(klon)
	316	REAL flwc_s(klon,klev), fiwc_s(klon,klev)
	317
[644]	318	cIM ISCCP simulator v3.4
[524]	319	c dans clesphys.h top_height, overlap
	320	cv3.4
	321	INTEGER debug, debugcol
[766]	322	cym INTEGER npoints
	323	cym PARAMETER(npoints=klon)
[524]	324	c
	325	INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
	326	INTEGER nregISCtot
	327	PARAMETER(nregISCtot=1)
	328	c
	329	c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
	330	c y compris pour 1 point
	331	c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
	332	c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
	333	INTEGER imin_debut, nbpti
	334	INTEGER jmin_debut, nbptj
[687]	335	cIM parametres ISCCP BEG
[828]	336	INTEGER nbapp_isccp
	337	! INTEGER nbapp_isccp,isccppas
	338	! PARAMETER(isccppas=6) !appel du simulateurs tous les 6pas de temps de la physique
	339	! !i.e. toutes les 3 heures
[952]	340	INTEGER n
[687]	341	INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
	342	DATA ifreq_isccp/3/
	343	SAVE ifreq_isccp
[766]	344	c$OMP THREADPRIVATE(ifreq_isccp)
[687]	345	CHARACTER*5 typinout(napisccp)
	346	DATA typinout/'i3od'/
[766]	347	SAVE typinout
	348	c$OMP THREADPRIVATE(typinout)
[687]	349	cIM verif boxptop BEG
	350	CHARACTER*1 verticaxe(napisccp)
	351	DATA verticaxe/'1'/
[766]	352	SAVE verticaxe
	353	c$OMP THREADPRIVATE(verticaxe)
[687]	354	cIM verif boxptop END
	355	INTEGER nvlev(napisccp)
	356	c INTEGER nvlev
	357	REAL t1, aa
	358	REAL seed_re(klon,napisccp)
[766]	359	cym !!!! A voir plus tard
	360	cym INTEGER iphy(iim,jjmp1)
[687]	361	cIM parametres ISCCP END
[524]	362	c
	363	c ncol = nb. de sous-colonnes pour chaque maille du GCM
[687]	364	c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM
[766]	365	c INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
	366	INTEGER,SAVE :: ncol(napisccp)
	367	INTEGER ncolmx, seed(klon,napisccp)
[687]	368	REAL nbsunlit(nregISCtot,klon,napisccp) !nbsunlit : moyenne de sunlit
[828]	369	c PARAMETER(ncolmx=1500)
	370	PARAMETER(ncolmx=300)
[687]	371	c
	372	cIM verif boxptop BEG
	373	REAL vertlev(ncolmx,napisccp)
	374	cIM verif boxptop END
	375	c
[766]	376	REAL,SAVE :: tautab_omp(0:255),tautab(0:255)
	377	INTEGER,SAVE :: invtau_omp(-20:45000),invtau(-20:45000)
	378	c$OMP THREADPRIVATE(tautab,invtau)
[524]	379	REAL emsfc_lw
	380	PARAMETER(emsfc_lw=0.99)
[644]	381	c REAL ran0 ! type for random number fuction
[524]	382	c
	383	REAL cldtot(klon,klev)
	384	c variables de haut en bas pour le simulateur ISCCP
	385	REAL dtau_s(klon,klev) !tau nuages startiformes
	386	REAL dtau_c(klon,klev) !tau nuages convectifs
	387	REAL dem_s(klon,klev) !emissivite nuages startiformes
	388	REAL dem_c(klon,klev) !emissivite nuages convectifs
	389	c
	390	c variables de haut en bas pour le simulateur ISCCP
	391	REAL pfull(klon,klev)
	392	REAL phalf(klon,klev+1)
	393	REAL qv(klon,klev)
	394	REAL cc(klon,klev)
	395	REAL conv(klon,klev)
	396	REAL dtau_sH2B(klon,klev)
	397	REAL dtau_cH2B(klon,klev)
	398	REAL at(klon,klev)
	399	REAL dem_sH2B(klon,klev)
	400	REAL dem_cH2B(klon,klev)
	401	c
[687]	402	INTEGER kmax, lmax, lmax3
	403	PARAMETER(kmax=8, lmax=8, lmax3=3)
[524]	404	INTEGER kmaxm1, lmaxm1
	405	PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
	406	INTEGER iimx7, jjmx7, jjmp1x7
	407	PARAMETER(iimx7=iimkmaxm1, jjmx7=jjmlmaxm1,
	408	.jjmp1x7=jjmp1*lmaxm1)
	409	c
[687]	410	c output from ISCCP simulator
	411	REAL fq_isccp(klon,kmaxm1,lmaxm1,napisccp)
	412	REAL fq_is_true(klon,kmaxm1,lmaxm1,napisccp)
	413	REAL totalcldarea(klon,napisccp)
	414	REAL meanptop(klon,napisccp)
	415	REAL meantaucld(klon,napisccp)
	416	REAL boxtau(klon,ncolmx,napisccp)
	417	REAL boxptop(klon,ncolmx,napisccp)
	418	REAL zx_tmp_fi3d_bx(klon,ncolmx)
	419	REAL zx_tmp_3d_bx(iim,jjmp1,ncolmx)
	420	c
	421	REAL cld_fi3d(klon,lmax3)
	422	REAL cld_3d(iim,jjmp1,lmax3)
	423	c
[524]	424	INTEGER iw, iwmax
	425	REAL wmin, pas_w
[766]	426	c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
	427	cIM 051005 PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
[687]	428	PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
[524]	429	REAL o500(klon)
	430	c
	431
	432	c sorties ISCCP
	433
	434	integer nid_isccp
[644]	435	save nid_isccp
[766]	436	c$OMP THREADPRIVATE(nid_isccp)
[524]	437
	438	REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
	439	DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
[766]	440	SAVE zx_tau
[687]	441	DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
[766]	442	SAVE zx_pc
	443	c$OMP THREADPRIVATE(zx_tau,zx_pc)
[524]	444	c cldtopres pression au sommet des nuages
[687]	445	REAL cldtopres(lmaxm1), cldtopres3(lmax3)
	446	DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
	447	DATA cldtopres3/440., 680., 1000./
[766]	448	SAVE cldtopres,cldtopres3
	449	c$OMP THREADPRIVATE(cldtopres,cldtopres3)
	450	cIM 051005 BEG
[524]	451	INTEGER komega, nhoriRD
	452
	453	c taulev: numero du niveau de tau dans les sorties ISCCP
	454	CHARACTER *4 taulev(kmaxm1)
[644]	455	c DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
	456	DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
	457	CHARACTER *3 pclev(lmaxm1)
	458	DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
[766]	459	SAVE taulev,pclev
	460	c$OMP THREADPRIVATE(taulev,pclev)
[644]	461	c
	462	c cnameisccp
	463	CHARACTER *27 cnameisccp(lmaxm1,kmaxm1)
[687]	464	cIM bad 151205 DATA cnameisccp/'pc< 50hPa, tau< 0.3',
	465	DATA cnameisccp/'pc= 50-180hPa, tau< 0.3',
[644]	466	. 'pc= 180-310hPa, tau< 0.3',
	467	. 'pc= 310-440hPa, tau< 0.3',
	468	. 'pc= 440-560hPa, tau< 0.3',
	469	. 'pc= 560-680hPa, tau< 0.3',
	470	. 'pc= 680-800hPa, tau< 0.3',
[687]	471	. 'pc= 800-1000hPa, tau< 0.3',
[644]	472	. 'pc= 50-180hPa, tau= 0.3-1.3',
	473	. 'pc= 180-310hPa, tau= 0.3-1.3',
	474	. 'pc= 310-440hPa, tau= 0.3-1.3',
	475	. 'pc= 440-560hPa, tau= 0.3-1.3',
	476	. 'pc= 560-680hPa, tau= 0.3-1.3',
	477	. 'pc= 680-800hPa, tau= 0.3-1.3',
[687]	478	. 'pc= 800-1000hPa, tau= 0.3-1.3',
[644]	479	. 'pc= 50-180hPa, tau= 1.3-3.6',
	480	. 'pc= 180-310hPa, tau= 1.3-3.6',
	481	. 'pc= 310-440hPa, tau= 1.3-3.6',
	482	. 'pc= 440-560hPa, tau= 1.3-3.6',
	483	. 'pc= 560-680hPa, tau= 1.3-3.6',
	484	. 'pc= 680-800hPa, tau= 1.3-3.6',
[687]	485	. 'pc= 800-1000hPa, tau= 1.3-3.6',
[644]	486	. 'pc= 50-180hPa, tau= 3.6-9.4',
	487	. 'pc= 180-310hPa, tau= 3.6-9.4',
	488	. 'pc= 310-440hPa, tau= 3.6-9.4',
	489	. 'pc= 440-560hPa, tau= 3.6-9.4',
	490	. 'pc= 560-680hPa, tau= 3.6-9.4',
	491	. 'pc= 680-800hPa, tau= 3.6-9.4',
[687]	492	. 'pc= 800-1000hPa, tau= 3.6-9.4',
[644]	493	. 'pc= 50-180hPa, tau= 9.4-23',
	494	. 'pc= 180-310hPa, tau= 9.4-23',
	495	. 'pc= 310-440hPa, tau= 9.4-23',
	496	. 'pc= 440-560hPa, tau= 9.4-23',
	497	. 'pc= 560-680hPa, tau= 9.4-23',
	498	. 'pc= 680-800hPa, tau= 9.4-23',
[687]	499	. 'pc= 800-1000hPa, tau= 9.4-23',
[644]	500	. 'pc= 50-180hPa, tau= 23-60',
	501	. 'pc= 180-310hPa, tau= 23-60',
	502	. 'pc= 310-440hPa, tau= 23-60',
	503	. 'pc= 440-560hPa, tau= 23-60',
	504	. 'pc= 560-680hPa, tau= 23-60',
	505	. 'pc= 680-800hPa, tau= 23-60',
[687]	506	. 'pc= 800-1000hPa, tau= 23-60',
[644]	507	. 'pc= 50-180hPa, tau> 60.',
	508	. 'pc= 180-310hPa, tau> 60.',
	509	. 'pc= 310-440hPa, tau> 60.',
	510	. 'pc= 440-560hPa, tau> 60.',
	511	. 'pc= 560-680hPa, tau> 60.',
[687]	512	. 'pc= 680-800hPa, tau> 60.',
	513	. 'pc= 800-1000hPa, tau> 60.'/
[766]	514	SAVE cnameisccp
	515	c$OMP THREADPRIVATE(cnameisccp)
[644]	516	c
	517	c REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
	518	c INTEGER nhorix7
[524]	519	cIM: region='3d' <==> sorties en global
	520	CHARACTER*3 region
	521	PARAMETER(region='3d')
	522	c
[644]	523	cIM ISCCP simulator v3.4
	524	c
[524]	525	logical ok_hf
[644]	526	c
[524]	527	integer nid_hf, nid_hf3d
[644]	528	save ok_hf, nid_hf, nid_hf3d
[766]	529	c$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
[524]	530	c QUESTION : noms de variables ?
	531
[828]	532	c#ifdef histhf
	533	c data ok_hf/.true./
	534	c#else
	535	c data ok_hf/.false./
	536	c#endif
[524]	537	INTEGER longcles
	538	PARAMETER ( longcles = 20 )
	539	REAL clesphy0( longcles )
	540	c
	541	c Variables quasi-arguments
	542	c
	543	REAL xjour
	544	SAVE xjour
[766]	545	c$OMP THREADPRIVATE(xjour)
[524]	546	c
	547	c
	548	c Variables propres a la physique
	549	c
	550	INTEGER radpas
	551	SAVE radpas ! frequence d'appel rayonnement
[766]	552	c$OMP THREADPRIVATE(radpas)
[524]	553	c
	554	cc INTEGER iflag_con
	555	cc SAVE iflag_con ! indicateur de la convection
	556	c
	557	INTEGER itap
	558	SAVE itap ! compteur pour la physique
[766]	559	c$OMP THREADPRIVATE(itap)
[524]	560	c
	561	real slp(klon) ! sea level pressure
	562	c
[782]	563	REAL fevap(klon,nbsrf)
	564	REAL fluxlat(klon,nbsrf)
[524]	565	c
[782]	566	REAL qsol(klon)
[883]	567	REAL,save :: solarlong0
[524]	568	c
	569	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
	570	c
[644]	571	cIM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
	572	REAL zulow(klon),zvlow(klon)
[524]	573	c
	574	INTEGER igwd,idx(klon),itest(klon)
	575	c
[782]	576	REAL agesno(klon,nbsrf)
[524]	577	c
[782]	578	c REAL,allocatable,save :: run_off_lic_0(:)
	579	cc$OMP THREADPRIVATE(run_off_lic_0)
[766]	580	cym SAVE run_off_lic_0
[524]	581	cKE43
	582	c Variables liees a la convection de K. Emanuel (sb):
	583	c
	584	REAL bas, top ! cloud base and top levels
	585	SAVE bas
	586	SAVE top
[766]	587	c$OMP THREADPRIVATE(bas, top)
[524]	588
	589	REAL wdn(klon), tdn(klon), qdn(klon)
[879]	590	c
	591	c=================================================================================================
	592	cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
	593	c Variables liées à la poche froide (jyg)
	594
	595	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
	596	REAL Vprecip(klon,klev) ! precipitation vertical profile
	597	c
	598	REAL wape_prescr, fip_prescr
	599	INTEGER it_wape_prescr
	600	SAVE wape_prescr, fip_prescr, it_wape_prescr
	601	c
	602	c variables supplementaires de concvl
	603	REAL Tconv(klon,klev)
	604	REAL ment(klon,klev,klev),sij(klon,klev,klev)
	605	REAL dd_t(klon,klev),dd_q(klon,klev)
	606	real alp_bl_prescr
	607	save alp_bl_prescr
	608	real ale_bl_prescr
	609	save ale_bl_prescr
	610	real ale_wake(klon)
	611	real alp_wake(klon)
	612	cRC
	613	c Variables liées à la poche froide (jyg et rr)
	614	c Version diagnostique pour l'instant : pas de rétroaction sur la convection
	615
	616	REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
	617
	618	REAL wake_dth(klon,klev) ! wake : temp pot difference
	619
	620	REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
	621	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta transported by LS omega
	622	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of large scale omega
	623	REAL wake_dtKE(klon,klev) ! Wake : differential heating (wake - unpertubed) CONV
	624	REAL wake_dqKE(klon,klev) ! Wake : differential moistening (wake - unpertubed) CONV
	625	REAL wake_dtPBL(klon,klev) ! Wake : differential heating (wake - unpertubed) PBL
	626	REAL wake_dqPBL(klon,klev) ! Wake : differential moistening (wake - unpertubed) PBL
	627	REAL wake_omg(klon,klev) ! Wake : velocity difference (wake - unpertubed)
	628	REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
	629	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
	630	REAL wake_spread(klon,klev) ! spreading term in wake_delt
[952]	631	c
[879]	632	cpourquoi y'a pas de save??
	633	REAL wake_h(klon) ! Wake : hauteur de la poche froide
[952]	634	c
[879]	635	INTEGER wake_k(klon) ! Wake sommet
	636	c
	637	REAL t_undi(klon,klev) ! temperature moyenne dans la zone non perturbee
	638	REAL q_undi(klon,klev) ! humidite moyenne dans la zone non perturbee
	639	c
	640	REAL wake_pe(klon) ! Wake potential energy - WAPE
	641
	642	REAL wake_gfl(klon) ! Gust Front Length
	643	REAL wake_dens(klon)
	644	c
	645	c
	646	REAL dt_dwn(klon,klev)
	647	REAL dq_dwn(klon,klev)
	648	REAL wdt_PBL(klon,klev)
	649	REAL udt_PBL(klon,klev)
	650	REAL wdq_PBL(klon,klev)
	651	REAL udq_PBL(klon,klev)
	652	REAL M_dwn(klon,klev)
	653	REAL M_up(klon,klev)
	654	REAL dt_a(klon,klev)
	655	REAL dq_a(klon,klev)
	656	c
	657	cRR:fin declarations poches froides
	658	c=======================================================================================================
	659
[524]	660	c Variables locales pour la couche limite (al1):
	661	c
	662	cAl1 REAL pblh(klon) ! Hauteur de couche limite
	663	cAl1 SAVE pblh
	664	c34EK
	665	c
	666	c Variables locales:
	667	c
	668	REAL cdragh(klon) ! drag coefficient pour T and Q
	669	REAL cdragm(klon) ! drag coefficient pour vent
	670	cAA
	671	cAA Pour phytrac
	672	cAA
	673	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
	674	REAL yu1(klon) ! vents dans la premiere couche U
	675	REAL yv1(klon) ! vents dans la premiere couche V
[782]	676
[766]	677	REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
[524]	678
[766]	679	c@$$ LOGICAL offline ! Controle du stockage ds "physique"
	680	c@$$ PARAMETER (offline=.false.)
	681	c@$$ INTEGER physid
[524]	682	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
	683	REAL frac_nucl(klon,klev) ! idem (nucleation)
	684	#ifdef INCA
[567]	685	INTEGER :: iii
[524]	686	REAL :: calday
	687	#endif
	688
[644]	689	cIM cf FH pour Tiedtke 080604
	690	REAL rain_tiedtke(klon),snow_tiedtke(klon)
	691	c
[766]	692	cIM 050204 END
[524]	693	REAL evap(klon), devap(klon) ! evaporation et sa derivee
	694	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
[782]	695
[524]	696	REAL bils(klon) ! bilan de chaleur au sol
[687]	697	REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
	698	C ! type de sous-surface et pondere par la fraction
[524]	699	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
	700	C ! type de sous-surface et pondere par la fraction
[782]	701	REAL fder(klon)
[524]	702	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
	703	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
	704	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
	705	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
	706	c
[782]	707	REAL frugs(klon,nbsrf)
[524]	708	REAL zxrugs(klon) ! longueur de rugosite
	709	c
	710	c Conditions aux limites
	711	c
	712	INTEGER julien
	713	c
	714	INTEGER lmt_pas
	715	SAVE lmt_pas ! frequence de mise a jour
[766]	716	c$OMP THREADPRIVATE(lmt_pas)
[524]	717	cIM
	718	REAL pctsrf_new(klon,nbsrf) !pourcentage surfaces issus d'ORCHIDEE
	719
[766]	720	cym SAVE pctsrf ! sous-fraction du sol
	721
[687]	722	cIM sorties
	723	REAL un_jour
	724	PARAMETER(un_jour=86400.)
[524]	725	c======================================================================
	726	c
	727	c Declaration des procedures appelees
	728	c
	729	EXTERNAL angle ! calculer angle zenithal du soleil
	730	EXTERNAL alboc ! calculer l'albedo sur ocean
	731	EXTERNAL ajsec ! ajustement sec
	732	EXTERNAL conlmd ! convection (schema LMD)
	733	cKE43
	734	EXTERNAL conema3 ! convect4.3
	735	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
	736	cAA
	737	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
	738	c ! stockage des coefficients necessaires au
	739	c ! lessivage OFF-LINE et ON-LINE
	740	EXTERNAL hgardfou ! verifier les temperatures
	741	EXTERNAL nuage ! calculer les proprietes radiatives
	742	EXTERNAL o3cm ! initialiser l'ozone
	743	EXTERNAL orbite ! calculer l'orbite terrestre
	744	EXTERNAL ozonecm ! prescrire l'ozone
	745	EXTERNAL phyetat0 ! lire l'etat initial de la physique
	746	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
	747	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
[879]	748	EXTERNAL suphel ! initialiser certaines constantes
[524]	749	EXTERNAL transp ! transport total de l'eau et de l'energie
	750	EXTERNAL ecribina ! ecrire le fichier binaire global
	751	EXTERNAL ecribins ! ecrire le fichier binaire global
	752	EXTERNAL ecrirega ! ecrire le fichier binaire regional
	753	EXTERNAL ecriregs ! ecrire le fichier binaire regional
	754	cIM
	755	EXTERNAL haut2bas !variables de haut en bas
	756	INTEGER lnblnk1
	757	EXTERNAL lnblnk1 !enleve les blancs a la fin d'une variable de type
	758	!caracter
[644]	759	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
	760	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
	761	c EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
	762	c EXTERNAL moyglo_aire !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
	763	c !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
[524]	764	c
	765	c Variables locales
	766	c
	767	REAL rhcl(klon,klev) ! humiditi relative ciel clair
	768	REAL dialiq(klon,klev) ! eau liquide nuageuse
	769	REAL diafra(klon,klev) ! fraction nuageuse
	770	REAL cldliq(klon,klev) ! eau liquide nuageuse
	771	REAL cldfra(klon,klev) ! fraction nuageuse
	772	REAL cldtau(klon,klev) ! epaisseur optique
	773	REAL cldemi(klon,klev) ! emissivite infrarouge
	774	c
	775	CXXX PB
	776	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
	777	REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
	778	REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
	779	REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
	780	c
	781	REAL zxfluxt(klon, klev)
	782	REAL zxfluxq(klon, klev)
	783	REAL zxfluxu(klon, klev)
	784	REAL zxfluxv(klon, klev)
	785	CXXX
[952]	786	c
[524]	787	REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
	788	REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
	789	c Le rayonnement n'est pas calcule tous les pas, il faut donc
	790	c sauvegarder les sorties du rayonnement
[766]	791	cym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
	792	cym SAVE sollwdownclr, toplwdown, toplwdownclr
	793	cym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
[524]	794	c
	795	INTEGER itaprad
	796	SAVE itaprad
[766]	797	c$OMP THREADPRIVATE(itaprad)
[524]	798	c
	799	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
	800	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
	801	c
	802	REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
	803	REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
	804	c
	805	REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
[782]	806	REAL zxsnow_dummy(klon)
[524]	807	c
	808	REAL dist, rmu0(klon), fract(klon)
	809	REAL zdtime, zlongi
	810	c
	811	CHARACTER*2 str2
	812	CHARACTER*2 iqn
	813	c
	814	REAL qcheck
	815	REAL z_avant(klon), z_apres(klon), z_factor(klon)
	816	LOGICAL zx_ajustq
	817	c
	818	REAL za, zb
	819	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
	820	real zqsat(klon,klev)
[909]	821	INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq, iff
[524]	822	REAL t_coup
	823	PARAMETER (t_coup=234.0)
	824	c
	825	REAL zphi(klon,klev)
[766]	826	cym A voir plus tard !!
	827	cym REAL zx_relief(iim,jjmp1)
	828	cym REAL zx_aire(iim,jjmp1)
[644]	829	c
[782]	830	c Grandeurs de sorties
[644]	831	REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
	832	REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
	833	REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
	834	REAL s_trmb3(klon)
[524]	835	cKE43
	836	c Variables locales pour la convection de K. Emanuel (sb):
	837	c
	838	REAL upwd(klon,klev) ! saturated updraft mass flux
	839	REAL dnwd(klon,klev) ! saturated downdraft mass flux
	840	REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
	841	REAL tvp(klon,klev) ! virtual temp of lifted parcel
	842	CHARACTER*40 capemaxcels !max(CAPE)
	843
	844	REAL rflag(klon) ! flag fonctionnement de convect
	845	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
	846	c -- convect43:
[644]	847	INTEGER ntra ! nb traceurs pour convect4.3
[524]	848	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
	849	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
	850	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
	851	REAL dplcldt(klon), dplcldr(klon)
	852	c? . condm_con(klon,klev),conda_con(klon,klev),
	853	c? . mr_con(klon,klev),ep_con(klon,klev)
	854	c? . ,sadiab(klon,klev),wadiab(klon,klev)
	855	c --
	856	c34EK
	857	c
	858	c Variables du changement
	859	c
	860	c con: convection
	861	c lsc: condensation a grande echelle (Large-Scale-Condensation)
	862	c ajs: ajustement sec
	863	c eva: evaporation de l'eau liquide nuageuse
	864	c vdf: couche limite (Vertical DiFfusion)
	865	REAL rneb(klon,klev)
[904]	866
	867	! tendance nulles
	868	REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
	869
[524]	870	c
[644]	871	*********************************************************
	872	* declarations
	873
	874	*********************************************************
	875	cIM 081204 END
	876	c
[524]	877	REAL pmfu(klon,klev), pmfd(klon,klev)
	878	REAL pen_u(klon,klev), pen_d(klon,klev)
	879	REAL pde_u(klon,klev), pde_d(klon,klev)
	880	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
	881	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
	882	REAL prfl(klon,klev+1), psfl(klon,klev+1)
	883	c
[766]	884	REAL rain_lsc(klon)
	885	REAL snow_lsc(klon)
[524]	886	c
[766]	887	REAL ratqss(klon,klev),ratqsc(klon,klev)
[524]	888	real ratqsbas,ratqshaut
[766]	889	save ratqsbas,ratqshaut
[524]	890	real zpt_conv(klon,klev)
	891
	892	c Parametres lies au nouveau schema de nuages (SB, PDF)
	893	real fact_cldcon
	894	real facttemps
	895	logical ok_newmicro
	896	save ok_newmicro
[766]	897	c$OMP THREADPRIVATE(ok_newmicro)
[524]	898	save fact_cldcon,facttemps
[766]	899	c$OMP THREADPRIVATE(fact_cldcon,facttemps)
[524]	900	real facteur
	901
	902	integer iflag_cldcon
	903	save iflag_cldcon
[766]	904	c$OMP THREADPRIVATE(iflag_cldcon)
[524]	905	logical ptconv(klon,klev)
[644]	906	cIM cf. AM 081204 BEG
	907	logical ptconvth(klon,klev)
	908	cIM cf. AM 081204 END
[524]	909	c
	910	c Variables liees a l'ecriture de la bande histoire physique
	911	c
[644]	912	c======================================================================
[524]	913	c
[644]	914	cIM cf. AM 081204 BEG
	915	c declarations pour sortir sur une sous-region
	916	integer imin_ins,imax_ins,jmin_ins,jmax_ins
	917	save imin_ins,imax_ins,jmin_ins,jmax_ins
[766]	918	c$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
[644]	919	c real lonmin_ins,lonmax_ins,latmin_ins
	920	c s ,latmax_ins
	921	c data lonmin_ins,lonmax_ins,latmin_ins
	922	c s ,latmax_ins/
	923	c valeurs initiales s -5.,20.,41.,55./
	924	c s 100.,130.,-20.,20./
	925	c s -180.,180.,-90.,90./
	926	c======================================================================
	927	cIM cf. AM 081204 END
	928
[524]	929	c
	930	integer itau_w ! pas de temps ecriture = itap + itau_phy
	931	c
	932	c
	933	c Variables locales pour effectuer les appels en serie
	934	c
	935	REAL zx_rh(klon,klev)
[687]	936	cIM RH a 2m (la surface)
	937	REAL rh2m(klon), qsat2m(klon)
	938	REAL tpot(klon), tpote(klon)
	939	REAL Lheat
[524]	940
	941	INTEGER length
	942	PARAMETER ( length = 100 )
	943	REAL tabcntr0( length )
	944	c
	945	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
[687]	946	cIM
	947	INTEGER ndex2d1(iwmax)
[644]	948	c
	949	cIM AMIP2 BEG
	950	REAL moyglo, mountor
	951	cIM 141004 BEG
	952	REAL zustrdr(klon), zvstrdr(klon)
	953	REAL zustrli(klon), zvstrli(klon)
	954	REAL zustrph(klon), zvstrph(klon)
	955	REAL aam, torsfc
	956	cIM 141004 END
	957	cIM 190504 BEG
	958	INTEGER ij, imp1jmp1
	959	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
[766]	960	cym A voir plus tard
[644]	961	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
	962	REAL padyn(iim+1,jjmp1,klev+1)
	963	REAL dudyn(iim+1,jjmp1,klev)
	964	REAL rlatdyn(iim+1,jjmp1)
	965	cIM 190504 END
	966	LOGICAL ok_msk
	967	REAL msk(klon)
	968	cIM
	969	REAL airetot, pi
[766]	970	cym A voir plus tard
	971	cym REAL zm_wo(jjmp1, klev)
[644]	972	cIM AMIP2 END
	973	c
[524]	974	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
	975	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
[766]	976	c#ifdef histmthNMC
	977	cym A voir plus tard !!!!
	978	cym REAL zx_tmp_NC(iim,jjmp1,nlevSTD)
[694]	979	REAL zx_tmp_fiNC(klon,nlevSTD)
[766]	980	c#endif
[644]	981	REAL*8 zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
[524]	982	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
	983	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
	984	c
[644]	985	INTEGER nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
[687]	986	INTEGER nid_ctesGCM
[644]	987	SAVE nid_day, nid_mth, nid_ins, nid_nmc, nid_day_seri
[687]	988	SAVE nid_ctesGCM
[766]	989	c$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins, nid_nmc)
	990	c$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
[524]	991	c
[644]	992	cIM 280405 BEG
	993	INTEGER nid_bilKPins, nid_bilKPave
	994	SAVE nid_bilKPins, nid_bilKPave
[766]	995	c$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
[644]	996	c
	997	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
	998	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
	999	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
	1000	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
	1001	c
	1002	cIM 280405 END
	1003	c
[687]	1004	INTEGER nhori, nvert, nvert1, nvert3
	1005	REAL zsto, zsto1, zsto2
	1006	REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
	1007	REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
	1008	REAL zout_isccp(napisccp)
	1009	SAVE zcals, zcalh, zoutj, zout_isccp
[766]	1010	c$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
[687]	1011
[524]	1012	real zjulian
	1013	save zjulian
[766]	1014	c$OMP THREADPRIVATE(zjulian)
[524]	1015
	1016	character*20 modname
	1017	character*80 abort_message
	1018	logical ok_sync
	1019	real date0
	1020	integer idayref
	1021
	1022	C essai writephys
	1023	integer fid_day, fid_mth, fid_ins
	1024	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
	1025	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
	1026	parameter (prof2d_on = 1, prof3d_on = 2,
	1027	. prof2d_av = 3, prof3d_av = 4)
	1028	character*30 nom_fichier
	1029	character*10 varname
	1030	character*40 vartitle
	1031	character*20 varunits
	1032	C Variables liees au bilan d'energie et d'enthalpi
	1033	REAL ztsol(klon)
	1034	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	1035	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	1036	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
	1037	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
[766]	1038	c$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
	1039	c$OMP+ h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
[524]	1040	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
	1041	REAL d_h_vcol_phy
	1042	REAL fs_bound, fq_bound
	1043	SAVE d_h_vcol_phy
[766]	1044	c$OMP THREADPRIVATE(d_h_vcol_phy)
[524]	1045	REAL zero_v(klon)
	1046	CHARACTER*15 ztit
[766]	1047	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
	1048	SAVE ip_ebil
	1049	DATA ip_ebil/0/
	1050	c$OMP THREADPRIVATE(ip_ebil)
	1051	INTEGER if_ebil ! level for energy conserv. dignostics
	1052	SAVE if_ebil
	1053	c$OMP THREADPRIVATE(if_ebil)
[524]	1054	c+jld ec_conser
	1055	REAL ZRCPD
	1056	c-jld ec_conser
[782]	1057	REAL t2m(klon,nbsrf) ! temperature a 2m
	1058	REAL q2m(klon,nbsrf) ! humidite a 2m
	1059
[524]	1060	cIM: t2m, q2m, u10m, v10m et t2mincels, t2maxcels
	1061	REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
	1062	REAL zu10m(klon), zv10m(klon) !vents a 10m moyennes s/1 maille
	1063	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
[644]	1064	CHARACTER*40 tinst, tave, typeval
[524]	1065	cjq Aerosol effects (Johannes Quaas, 27/11/2003)
	1066	REAL sulfate(klon, klev) ! SO4 aerosol concentration [ug/m3]
	1067
	1068	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
	1069
	1070	REAL re(klon, klev) ! Cloud droplet effective radius
	1071	REAL fl(klon, klev) ! denominator of re
	1072
	1073	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
	1074
	1075	! Aerosol optical properties
[955]	1076	#ifdef INCA
[864]	1077	! Aerosol optical properties by INCA model
	1078	CHARACTER*4 :: rfname(9)
	1079	#endif
[524]	1080	REAL aerindex(klon) ! POLDER aerosol index
	1081
	1082	! Parameters
	1083	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
	1084	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
[559]	1085	cym
	1086	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
[766]	1087	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
[955]	1088	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
	1089	! false : lecture des aerosol dans un fichier
	1090	c$OMP THREADPRIVATE(aerosol_couple)
[674]	1091
[524]	1092	c
	1093	c Declaration des constantes et des fonctions thermodynamiques
	1094	c
[766]	1095	LOGICAL,SAVE :: first=.true.
	1096	c$OMP THREADPRIVATE(first)
[524]	1097	#include "YOMCST.h"
	1098	#include "YOETHF.h"
	1099	#include "FCTTRE.h"
[687]	1100	cIM 100106 BEG : pouvoir sortir les ctes de la physique
	1101	#include "conema3.h"
	1102	#include "fisrtilp.h"
	1103	#include "nuage.h"
	1104	#include "compbl.h"
	1105	cIM 100106 END : pouvoir sortir les ctes de la physique
	1106	c
[524]	1107	c======================================================================
[879]	1108	! Ecriture eventuelle d'un profil verticale en entree de la physique.
	1109	! Utilise notamment en 1D mais peut etre active egalement en 3D
	1110	! en imposant la valeur de igout.
	1111	c======================================================================
[766]	1112
[942]	1113	if (prt_level.ge.1) then
[950]	1114	igout=klon/2+1/klon
[879]	1115	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	1116	write(lunout,*)
	1117	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys'
	1118	write(lunout,*)
	1119	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys
	1120
	1121	write(lunout,*) 'papers, play, phi, u, v, t, omega'
	1122	do k=1,nlev
	1123	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k),
	1124	s u(igout,k),v(igout,k),t(igout,k),omega(igout,k)
	1125	enddo
	1126	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
	1127	do k=1,nlev
	1128	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
	1129	enddo
	1130	endif
	1131
	1132	c======================================================================
	1133
[766]	1134	cym => necessaire pour iflag_con != 2
	1135	pmfd(:,:) = 0.
	1136	pen_u(:,:) = 0.
	1137	pen_d(:,:) = 0.
	1138	pde_d(:,:) = 0.
	1139	pde_u(:,:) = 0.
	1140	aam=0.
	1141	torsfc=0.
	1142
	1143	if (first) then
	1144
[879]	1145	cCR:nvelles variables convection/poches froides
[766]	1146
[909]	1147	print*, '================================================='
	1148	print*, 'Allocation des variables locales et sauvegardees'
	1149	call phys_local_var_init
	1150	call phys_state_var_init
	1151	print*, '================================================='
	1152
[766]	1153	paire_ter(:)=0.
	1154	clwcon(:,:)=0.
	1155	rnebcon(:,:)=0.
	1156	ratqs(:,:)=0.
	1157	sollw(:)=0.
	1158	cym Attention pbase pas initialise dans concvl !!!!
	1159	pbase(:)=0
	1160
	1161	first=.false.
	1162
[904]	1163	endif ! fisrt
	1164
[766]	1165	modname = 'physiq'
[687]	1166	cIM
	1167	IF (ip_ebil_phy.ge.1) THEN
[524]	1168	DO i=1,klon
	1169	zero_v(i)=0.
	1170	END DO
	1171	END IF
	1172	ok_sync=.TRUE.
	1173	IF (nqmax .LT. 2) THEN
	1174	abort_message = 'eaux vapeur et liquide sont indispensables'
	1175	CALL abort_gcm (modname,abort_message,1)
	1176	ENDIF
	1177	IF (debut) THEN
[879]	1178	CALL suphel ! initialiser constantes et parametres phys.
[644]	1179	ENDIF
	1180
[942]	1181	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
[644]	1182
[878]	1183
[524]	1184	c======================================================================
	1185	xjour = rjourvrai
	1186	c
	1187	c Si c'est le debut, il faut initialiser plusieurs choses
	1188	c ********
	1189	c
	1190	IF (debut) THEN
	1191	C
[645]	1192	!rv
[879]	1193	cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
	1194	cde la convection a partir des caracteristiques du thermique
	1195	wght_th(:,:)=1.
	1196	lalim_conv(:)=1
	1197	cRC
[645]	1198	u10m(:,:)=0.
	1199	v10m(:,:)=0.
	1200	piz_ae(:,:,:)=0.
	1201	tau_ae(:,:,:)=0.
	1202	cg_ae(:,:,:)=0.
	1203	rain_con(:)=0.
	1204	snow_con(:)=0.
	1205	bl95_b0=0.
	1206	bl95_b1=0.
	1207	topswai(:)=0.
	1208	topswad(:)=0.
	1209	solswai(:)=0.
	1210	solswad(:)=0.
[955]	1211	#ifdef INCA
[864]	1212	tau_inca(:,:,:,:) = 0.
	1213	piz_inca(:,:,:,:) = 0.
	1214	cg_inca(:,:,:,:) = 0.
	1215	ccm(:,:,:) = 0.
	1216	topswai_inca(:) = 0.
	1217	topswad_inca(:) = 0.
	1218	topswad0_inca(:) = 0.
	1219	topsw_inca(:,:) = 0.
	1220	topsw0_inca(:,:) = 0.
	1221	solswai_inca(:) = 0.
	1222	solswad_inca(:) = 0.
	1223	solswad0_inca(:) = 0.
	1224	solsw_inca(:,:) = 0.
	1225	solsw0_inca(:,:) = 0.
	1226	#endif
[645]	1227	!rv
[864]	1228	!ACo
[645]	1229	rnebcon0(:,:) = 0.0
	1230	clwcon0(:,:) = 0.0
	1231	rnebcon(:,:) = 0.0
	1232	clwcon(:,:) = 0.0
[888]	1233	! fin anne
[645]	1234
[687]	1235	cIM
	1236	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
[524]	1237	c
	1238	c appel a la lecture du run.def physique
	1239	c
	1240	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
[883]	1241	. ok_instan, ok_hf,
[889]	1242	. solarlong0,seuil_inversion,
[879]	1243	. fact_cldcon, facttemps,ok_newmicro,iflag_radia,
[878]	1244	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
[955]	1245	. ok_ade, ok_aie, aerosol_couple,
[541]	1246	. bl95_b0, bl95_b1,
[879]	1247	. iflag_thermals,nsplit_thermals,
	1248	cnv flags pour la convection et les poches froides
	1249	. iflag_coupl,iflag_clos,iflag_wake)
[524]	1250
[879]	1251	print*,'iflag_coupl,iflag_clos,iflag_wake',
	1252	. iflag_coupl,iflag_clos,iflag_wake
[956]	1253	print*,'CYCLE_DIURNE', cycle_diurne
[879]	1254
[524]	1255	c
	1256	c
	1257	c Initialiser les compteurs:
	1258	c
	1259	itap = 0
	1260	itaprad = 0
[782]	1261
[878]	1262	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1263	!! Un petit travail à faire ici.
	1264	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1265
	1266	if (iflag_pbl>1) then
	1267	PRINT*, "Using method MELLOR&YAMADA"
	1268	endif
	1269
[956]	1270	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1271	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
	1272	! dyn3d
	1273	! Attention : la version precedente n'etait pas tres propre.
	1274	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
	1275	! pour obtenir le meme resultat.
	1276	dtime=pdtphys
	1277	radpas = NINT( 86400./dtime/nbapp_rad)
	1278	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1279
	1280
[524]	1281	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
[782]	1282	. rlat,rlon,pctsrf, ftsol,
	1283	. ocean, ok_veget,
[888]	1284	. falb1, falb2, rain_fall,snow_fall,
[782]	1285	. solsw, sollw,
	1286	. radsol,clesphy0,
[524]	1287	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
[878]	1288	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
[937]	1289	. pbl_tke, zmax0, f0, ema_work1, ema_work2)
[524]	1290
[878]	1291
	1292
	1293
	1294	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1295
	1296
[766]	1297	DO i=1,klon
	1298	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
	1299	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
	1300	$ THEN
[782]	1301	WRITE(,)
	1302	$ 'physiq apres lecture de restart: pb sous surface au point ',
	1303	$ i, pctsrf(i, 1 : nbsrf)
[766]	1304	ENDIF
	1305	ENDDO
[782]	1306
[524]	1307	c
	1308	C on remet le calendrier a zero
	1309	c
	1310	IF (raz_date .eq. 1) THEN
	1311	itau_phy = 0
	1312	ENDIF
	1313
[644]	1314	cIM cf. AM 081204 BEG
	1315	PRINT*,'cycle_diurne3 =',cycle_diurne
	1316	cIM cf. AM 081204 END
[524]	1317	c
[782]	1318	CALL printflag( tabcntr0,radpas,ok_journe,
[524]	1319	, ok_instan, ok_region )
	1320	c
	1321	IF (ABS(dtime-pdtphys).GT.0.001) THEN
	1322	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
	1323	. pdtphys
	1324	abort_message='Pas physique n est pas correct '
[878]	1325	! call abort_gcm(modname,abort_message,1)
	1326	dtime=pdtphys
[524]	1327	ENDIF
	1328	IF (nlon .NE. klon) THEN
	1329	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
	1330	. klon
	1331	abort_message='nlon et klon ne sont pas coherents'
	1332	call abort_gcm(modname,abort_message,1)
	1333	ENDIF
	1334	IF (nlev .NE. klev) THEN
	1335	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
	1336	. klev
	1337	abort_message='nlev et klev ne sont pas coherents'
	1338	call abort_gcm(modname,abort_message,1)
	1339	ENDIF
	1340	c
	1341	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
	1342	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	1343	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	1344	abort_message='Nbre d appels au rayonnement insuffisant'
	1345	call abort_gcm(modname,abort_message,1)
	1346	ENDIF
	1347	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
	1348	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
	1349	. ok_cvl
	1350	c
	1351	cKE43
	1352	c Initialisation pour la convection de K.E. (sb):
	1353	IF (iflag_con.GE.3) THEN
	1354
	1355	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
[687]	1356	WRITE(lunout,*)
	1357	. "On va utiliser le melange convectif des traceurs qui"
	1358	WRITE(lunout,*)"est calcule dans convect4.3"
	1359	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
[524]	1360
	1361	DO i = 1, klon
	1362	ema_cbmf(i) = 0.
	1363	ema_pcb(i) = 0.
	1364	ema_pct(i) = 0.
	1365	ema_workcbmf(i) = 0.
	1366	ENDDO
	1367	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
	1368	DO i = 1, klon
	1369	ibas_con(i) = 1
[619]	1370	itop_con(i) = 1
[524]	1371	ENDDO
	1372	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
[879]	1373	c===============================================================================
	1374	cCR:04.12.07: initialisations poches froides
	1375	c Controle de ALE et ALP pour la fermeture convective (jyg)
	1376	if (iflag_wake.eq.1) then
	1377	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
	1378	s ,alp_bl_prescr, ale_bl_prescr)
	1379	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
	1380	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
	1381	endif
[524]	1382
[879]	1383	do i = 1,klon
	1384	wake_s(i) = 0.
	1385	wake_fip(i) = 0.
	1386	wake_cstar(i) = 0.
	1387	DO k=1,klev
	1388	wake_deltat(i,k)=0.
	1389	wake_deltaq(i,k)=0.
	1390	ENDDO
	1391	enddo
	1392	c================================================================================
	1393
[524]	1394	ENDIF
	1395
[878]	1396	rugoro=0.
[524]	1397	c34EK
	1398	IF (ok_orodr) THEN
[878]	1399
	1400	rugoro=0.
	1401
	1402	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1403	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
	1404	! justement quand ok_orodr = false.
	1405	! ce rugoro est utilise par la couche limite et fait double emploi
	1406	! avec les paramétrisations spécifiques de Francois Lott.
	1407	! DO i=1,klon
	1408	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1409	! ENDDO
	1410	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1411
[782]	1412	CALL SUGWD(klon,klev,paprs,pplay)
	1413	DO i=1,klon
	1414	zuthe(i)=0.
	1415	zvthe(i)=0.
	1416	if(zstd(i).gt.10.)then
	1417	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	1418	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	1419	endif
	1420	ENDDO
[524]	1421	ENDIF
	1422	c
	1423	c
	1424	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
	1425	WRITE(lunout,*)'La frequence de lecture surface est de ',
	1426	. lmt_pas
	1427	c
[644]	1428	cIM200505 ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
	1429	c IF (ok_mensuel) THEN
	1430	c WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
	1431	c . ecrit_mth
	1432	c ENDIF
	1433	c ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
	1434	c IF (ok_journe) THEN
	1435	c WRITE(lunout,*)'La frequence de sortie journaliere est de ',
	1436	c . ecrit_day
	1437	c ENDIF
	1438	cIM 130904 BEG
	1439	cIM 080205 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
	1440	cIM 170305
	1441	c ecrit_hf = 86400./dtime/12. ! toutes les 2h
	1442	cIM 230305
	1443	cIM200505 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
	1444	c
	1445	cIM200505 ecrit_hf2mth = ecrit_day/ecrit_hf*30
	1446	c
	1447	cIM200505 IF (ok_journe) THEN
	1448	cIM200505 WRITE(lunout,*)'La frequence de sortie hf est de ',
	1449	cIM200505 . ecrit_hf
	1450	cIM200505 ENDIF
	1451	cIM 130904 END
[524]	1452	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
	1453	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
[644]	1454	c ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
	1455	c ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
	1456	cIM200505 ecrit_ins = NINT(86400./dtime/8.) ! toutes les trois heures
	1457	cIM200505 IF (ok_instan) THEN
	1458	cIM200505 WRITE(lunout,*)'La frequence de sortie instant. est de ',
	1459	cIM200505 . ecrit_ins
	1460	cIM200505 ENDIF
	1461	cIM200505 ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
	1462	cIM200505 IF (ok_region) THEN
	1463	cIM200505 WRITE(lunout,*)'La frequence de sortie region est de ',
	1464	cIM200505 . ecrit_reg
	1465	cIM200505 ENDIF
[687]	1466	cIM 030306 BEG
	1467	cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
	1468	cIM : ne pas modifier ecrit_hf2mth
[524]	1469	c
[933]	1470	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
	1471	ecrit_hf2mth = ecrit_mth/ecrit_hf
[702]	1472	c ecrit_ins en secondes, chaque pas de temps de la physique
[687]	1473	ecrit_ins = dtime
	1474	cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
	1475	ecrit_hf = ecrit_hf * un_jour
[828]	1476	!IM
	1477	IF(ecrit_day.LE.1.) THEN
	1478	ecrit_day = ecrit_day * un_jour !en secondes
	1479	ENDIF
	1480	!IM
[687]	1481	ecrit_mth = ecrit_mth * un_jour
	1482	ecrit_reg = ecrit_reg * un_jour
[702]	1483	ecrit_tra = ecrit_tra * un_jour
[828]	1484	ecrit_ISCCP = ecrit_ISCCP * un_jour
	1485	c
[933]	1486	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
	1487	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
	1488	. ecrit_hf2mth
[687]	1489	cIM 030306 END
[782]	1490
[524]	1491	capemaxcels = 't_max(X)'
	1492	t2mincels = 't_min(X)'
	1493	t2maxcels = 't_max(X)'
[644]	1494	tinst = 'inst(X)'
	1495	tave = 'ave(X)'
	1496	cIM cf. AM 081204 BEG
	1497	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
	1498	cIM cf. AM 081204 END
[524]	1499	c
	1500	c=============================================================
	1501	c Initialisation des sorties
	1502	c=============================================================
	1503
	1504	#ifdef CPP_IOIPSL
	1505
[909]	1506	c Commente par abderrahmane 11 2 08
	1507	c#ifdef histhf
	1508	c#include "ini_histhf.h"
	1509	c#endif
[524]	1510
[909]	1511	c#ifdef histday
	1512	c#include "ini_histday.h"
[644]	1513	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
	1514	c#include "ini_bilKP_ins.h"
	1515	c#include "ini_bilKP_ave.h"
[909]	1516	c#endif
[524]	1517
[909]	1518	c#ifdef histmth
	1519	c#include "ini_histmth.h"
	1520	c#endif
[524]	1521
[909]	1522	c#ifdef histins
	1523	c#include "ini_histins.h"
	1524	c#endif
[524]	1525
[929]	1526	call phys_output_open(jjmp1,nqmax,nlevSTD,clevSTD,nbteta,
	1527	& ctetaSTD,dtime,presnivs,ok_veget,
	1528	& ocean,iflag_pbl,ok_mensuel,ok_journe,
	1529	& ok_hf,ok_instan,nid_files)
[909]	1530
[524]	1531	#ifdef histISCCP
	1532	#include "ini_histISCCP.h"
	1533	#endif
	1534
	1535	#ifdef histmthNMC
	1536	#include "ini_histmthNMC.h"
	1537	#endif
	1538
[687]	1539	#include "ini_histday_seri.h"
[524]	1540
[687]	1541	#include "ini_paramLMDZ_phy.h"
[524]	1542
[644]	1543	#endif
	1544
[524]	1545	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
	1546	date0 = zjulian
	1547	C date0 = day_ini
	1548	WRITE(,) 'physiq date0 : ',date0
	1549	c
	1550	c
	1551	c
	1552	c Prescrire l'ozone dans l'atmosphere
	1553	c
	1554	c
	1555	cc DO i = 1, klon
	1556	cc DO k = 1, klev
	1557	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
	1558	cc ENDDO
	1559	cc ENDDO
	1560	c
	1561	#ifdef INCA
[766]	1562	call VTe(VTphysiq)
	1563	call VTb(VTinca)
[524]	1564	iii = MOD(NINT(xjour),360)
	1565	calday = FLOAT(iii) + gmtime
	1566	WRITE(lunout,*) 'initial time ', xjour, calday
	1567	#ifdef INCAINFO
	1568	WRITE(lunout,*) 'Appel CHEMINI ...'
	1569	#endif
[766]	1570	CALL chemini(
[524]	1571	$ rg,
	1572	$ ra,
	1573	$ airephy,
	1574	$ rlat,
	1575	$ rlon,
	1576	$ presnivs,
	1577	$ calday,
	1578	$ klon,
	1579	$ nqmax,
	1580	$ pdtphys,
[567]	1581	$ annee_ref,
[524]	1582	$ day_ini)
	1583	#ifdef INCAINFO
	1584	WRITE(lunout,*) 'OK.'
	1585	#endif
[766]	1586	call VTe(VTinca)
	1587	call VTb(VTphysiq)
[524]	1588	#endif
	1589	c
	1590	ENDIF
	1591	c
	1592	c ************** Fin de IF ( debut ) *************
	1593	c
[904]	1594
	1595	! Tendances bidons pour les processus qui n'affectent pas certaines
	1596	! variables.
	1597	du0(:,:)=0.
	1598	dv0(:,:)=0.
	1599	dq0(:,:)=0.
	1600	dql0(:,:)=0.
[524]	1601	c
	1602	c Mettre a zero des variables de sortie (pour securite)
	1603	c
	1604	DO i = 1, klon
	1605	d_ps(i) = 0.0
	1606	ENDDO
	1607	DO k = 1, klev
	1608	DO i = 1, klon
	1609	d_t(i,k) = 0.0
	1610	d_u(i,k) = 0.0
	1611	d_v(i,k) = 0.0
	1612	ENDDO
	1613	ENDDO
	1614	DO iq = 1, nqmax
	1615	DO k = 1, klev
	1616	DO i = 1, klon
	1617	d_qx(i,k,iq) = 0.0
	1618	ENDDO
	1619	ENDDO
	1620	ENDDO
[660]	1621	da(:,:)=0.
	1622	mp(:,:)=0.
	1623	phi(:,:,:)=0.
[524]	1624	c
	1625	c Ne pas affecter les valeurs entrees de u, v, h, et q
	1626	c
	1627	DO k = 1, klev
	1628	DO i = 1, klon
	1629	t_seri(i,k) = t(i,k)
	1630	u_seri(i,k) = u(i,k)
	1631	v_seri(i,k) = v(i,k)
	1632	q_seri(i,k) = qx(i,k,ivap)
	1633	ql_seri(i,k) = qx(i,k,iliq)
	1634	qs_seri(i,k) = 0.
	1635	ENDDO
	1636	ENDDO
	1637	IF (nqmax.GE.3) THEN
	1638	DO iq = 3, nqmax
	1639	DO k = 1, klev
	1640	DO i = 1, klon
	1641	tr_seri(i,k,iq-2) = qx(i,k,iq)
	1642	ENDDO
	1643	ENDDO
	1644	ENDDO
	1645	ELSE
	1646	DO k = 1, klev
	1647	DO i = 1, klon
	1648	tr_seri(i,k,1) = 0.0
	1649	ENDDO
	1650	ENDDO
	1651	ENDIF
	1652	C
	1653	DO i = 1, klon
	1654	ztsol(i) = 0.
	1655	ENDDO
	1656	DO nsrf = 1, nbsrf
	1657	DO i = 1, klon
	1658	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
	1659	ENDDO
	1660	ENDDO
[687]	1661	cIM
	1662	IF (ip_ebil_phy.ge.1) THEN
[524]	1663	ztit='after dynamic'
[687]	1664	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
[524]	1665	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1666	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1667	C Comme les tendances de la physique sont ajoute dans la dynamique,
	1668	C on devrait avoir que la variation d'entalpie par la dynamique
	1669	C est egale a la variation de la physique au pas de temps precedent.
	1670	C Donc la somme de ces 2 variations devrait etre nulle.
[687]	1671	call diagphy(airephy,ztit,ip_ebil_phy
[524]	1672	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1673	e , zero_v, zero_v, zero_v, ztsol
	1674	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
	1675	s , fs_bound, fq_bound )
	1676	END IF
	1677
	1678	c Diagnostiquer la tendance dynamique
	1679	c
	1680	IF (ancien_ok) THEN
	1681	DO k = 1, klev
	1682	DO i = 1, klon
	1683	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	1684	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
	1685	ENDDO
	1686	ENDDO
	1687	ELSE
	1688	DO k = 1, klev
	1689	DO i = 1, klon
	1690	d_t_dyn(i,k) = 0.0
	1691	d_q_dyn(i,k) = 0.0
	1692	ENDDO
	1693	ENDDO
	1694	ancien_ok = .TRUE.
	1695	ENDIF
	1696	c
	1697	c Ajouter le geopotentiel du sol:
	1698	c
	1699	DO k = 1, klev
	1700	DO i = 1, klon
	1701	zphi(i,k) = pphi(i,k) + pphis(i)
	1702	ENDDO
	1703	ENDDO
	1704	c
	1705	c Verifier les temperatures
	1706	c
[687]	1707	cIM BEG
	1708	IF (check) THEN
	1709	amn=MIN(ftsol(1,is_ter),1000.)
	1710	amx=MAX(ftsol(1,is_ter),-1000.)
	1711	DO i=2, klon
	1712	amn=MIN(ftsol(i,is_ter),amn)
	1713	amx=MAX(ftsol(i,is_ter),amx)
	1714	ENDDO
	1715	c
	1716	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
	1717	ENDIF !(check) THEN
	1718	cIM END
	1719	c
[524]	1720	CALL hgardfou(t_seri,ftsol,'debutphy')
	1721	c
[687]	1722	cIM BEG
	1723	IF (check) THEN
	1724	amn=MIN(ftsol(1,is_ter),1000.)
	1725	amx=MAX(ftsol(1,is_ter),-1000.)
	1726	DO i=2, klon
	1727	amn=MIN(ftsol(i,is_ter),amn)
	1728	amx=MAX(ftsol(i,is_ter),amx)
	1729	ENDDO
	1730	c
	1731	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
	1732	ENDIF !(check) THEN
	1733	cIM END
	1734	c
[524]	1735	c Incrementer le compteur de la physique
	1736	c
	1737	itap = itap + 1
	1738	julien = MOD(NINT(xjour),360)
	1739	if (julien .eq. 0) julien = 360
	1740	c
	1741	c Mettre en action les conditions aux limites (albedo, sst, etc.).
	1742	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
	1743	c
	1744	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
	1745	WRITE(lunout,*)' PHYS cond julien ',julien
	1746	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
	1747	ENDIF
	1748	c
	1749	c Re-evaporer l'eau liquide nuageuse
	1750	c
	1751	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
	1752	DO i = 1, klon
	1753	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1754	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1755	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1756	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
	1757	zb = MAX(0.0,ql_seri(i,k))
	1758	za = - MAX(0.0,ql_seri(i,k))
	1759	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
	1760	t_seri(i,k) = t_seri(i,k) + za
	1761	q_seri(i,k) = q_seri(i,k) + zb
	1762	ql_seri(i,k) = 0.0
	1763	d_t_eva(i,k) = za
	1764	d_q_eva(i,k) = zb
	1765	ENDDO
	1766	ENDDO
[687]	1767	cIM
	1768	IF (ip_ebil_phy.ge.2) THEN
[524]	1769	ztit='after reevap'
[687]	1770	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
[524]	1771	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1772	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	1773	call diagphy(airephy,ztit,ip_ebil_phy
[524]	1774	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1775	e , zero_v, zero_v, zero_v, ztsol
	1776	e , d_h_vcol, d_qt, d_ec
	1777	s , fs_bound, fq_bound )
	1778	C
	1779	END IF
[782]	1780
[524]	1781	c
[883]	1782	c=========================================================================
	1783	! Calculs de l'orbite.
	1784	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
	1785	! doit donc etre placé avant radlwsw et pbl_surface
	1786
	1787	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
	1788	! solarlong0
	1789
	1790	if (solarlong0<-999.) then
	1791	CALL orbite(FLOAT(julien),zlongi,dist)
	1792	else
	1793	zlongi=solarlong0 ! longitude solaire vraie
	1794	dist=1. ! distance au soleil / moyenne
	1795	endif
	1796
[942]	1797	if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
[883]	1798
	1799	! Avec ou sans cycle diurne
[524]	1800	IF (cycle_diurne) THEN
	1801	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
	1802	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
	1803	ELSE
	1804	rmu0 = -999.999
	1805	ENDIF
	1806
[766]	1807	if (mydebug) then
	1808	call writefield_phy('u_seri',u_seri,llm)
	1809	call writefield_phy('v_seri',v_seri,llm)
	1810	call writefield_phy('t_seri',t_seri,llm)
	1811	call writefield_phy('q_seri',q_seri,llm)
	1812	endif
[782]	1813
	1814	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1815	c Appel au pbl_surface : Planetary Boudary Layer et Surface
	1816	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
	1817	c turbulent du couche limit.
	1818	c
	1819	c Certains varibales de sorties de pbl_surface sont utiliser que pour
	1820	c ecriture des fihiers hist_XXXX.nc, ces sont :
	1821	c qsol, zq2m, s_pblh, s_lcl,
	1822	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
	1823	c s_therm, s_trmb1, s_trmb2, s_trmb3,
	1824	c zxrugs, zu10m, zv10m, fder,
	1825	c zxqsurf, rh2m, zxfluxu, zxfluxv,
	1826	c frugs, agesno, fsollw, fsolsw,
	1827	c d_ts, fevap, fluxlat, t2m,
	1828	c wfbils, wfbilo, fluxt, fluxu, fluxv,
[687]	1829	c
[782]	1830	c Certains ne sont pas utiliser du tout :
	1831	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
[687]	1832	c
[883]	1833
[782]	1834	CALL pbl_surface(
	1835	e dtime, date0, itap, julien,
	1836	e debut, lafin,
	1837	e rlon, rlat, rugoro, rmu0,
	1838	e rain_fall, snow_fall, solsw, sollw,
	1839	e t_seri, q_seri, u_seri, v_seri,
	1840	e pplay, paprs, pctsrf,
[888]	1841	+ ftsol, falb1, falb2, u10m, v10m,
[782]	1842	s sollwdown, cdragh, cdragm, yu1, yv1,
[888]	1843	s albsol1, albsol2, sens, evap,
[782]	1844	s zxtsol, zxfluxlat, zt2m, qsat2m,
	1845	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
	1846	s ycoefh, pctsrf_new,
	1847	d qsol, zq2m, s_pblh, s_lcl,
	1848	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
	1849	d s_therm, s_trmb1, s_trmb2, s_trmb3,
	1850	d zxrugs, zu10m, zv10m, fder,
	1851	d zxqsurf, rh2m, zxfluxu, zxfluxv,
	1852	d frugs, agesno, fsollw, fsolsw,
	1853	d d_ts, fevap, fluxlat, t2m,
	1854	d wfbils, wfbilo, fluxt, fluxu, fluxv,
	1855	- dsens, devap, zxsnow,
[878]	1856	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
[524]	1857	c
[782]	1858	c
	1859	pctsrf(:,:) = pctsrf_new(:,:)
	1860
[904]	1861	!-----------------------------------------------------------------------------------------
	1862	! ajout des tendances de la diffusion turbulente
	1863	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
	1864	!-----------------------------------------------------------------------------------------
[766]	1865
	1866	if (mydebug) then
	1867	call writefield_phy('u_seri',u_seri,llm)
	1868	call writefield_phy('v_seri',v_seri,llm)
	1869	call writefield_phy('t_seri',t_seri,llm)
	1870	call writefield_phy('q_seri',q_seri,llm)
	1871	endif
	1872
	1873
[687]	1874	IF (ip_ebil_phy.ge.2) THEN
[782]	1875	ztit='after surface_main'
[687]	1876	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	1877	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1878	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	1879	call diagphy(airephy,ztit,ip_ebil_phy
[524]	1880	e , zero_v, zero_v, zero_v, zero_v, sens
	1881	e , evap , zero_v, zero_v, ztsol
	1882	e , d_h_vcol, d_qt, d_ec
	1883	s , fs_bound, fq_bound )
	1884	END IF
	1885
[881]	1886	c =================================================================== c
	1887	c Calcul de Qsat
	1888
	1889	DO k = 1, klev
	1890	DO i = 1, klon
	1891	zx_t = t_seri(i,k)
	1892	IF (thermcep) THEN
	1893	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	1894	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	1895	zx_qs = MIN(0.5,zx_qs)
	1896	zcor = 1./(1.-retv*zx_qs)
	1897	zx_qs = zx_qs*zcor
	1898	ELSE
	1899	IF (zx_t.LT.t_coup) THEN
	1900	zx_qs = qsats(zx_t)/pplay(i,k)
	1901	ELSE
	1902	zx_qs = qsatl(zx_t)/pplay(i,k)
	1903	ENDIF
	1904	ENDIF
	1905	zqsat(i,k)=zx_qs
	1906	ENDDO
	1907	ENDDO
	1908
[942]	1909	if (prt_level.ge.1) then
[881]	1910	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
	1911	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
	1912	endif
[524]	1913	c
	1914	c Appeler la convection (au choix)
	1915	c
	1916	DO k = 1, klev
	1917	DO i = 1, klon
	1918	conv_q(i,k) = d_q_dyn(i,k)
	1919	. + d_q_vdf(i,k)/dtime
	1920	conv_t(i,k) = d_t_dyn(i,k)
	1921	. + d_t_vdf(i,k)/dtime
	1922	ENDDO
	1923	ENDDO
	1924	IF (check) THEN
	1925	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	1926	WRITE(lunout,*) "avantcon=", za
	1927	ENDIF
	1928	zx_ajustq = .FALSE.
	1929	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
	1930	IF (zx_ajustq) THEN
	1931	DO i = 1, klon
	1932	z_avant(i) = 0.0
	1933	ENDDO
	1934	DO k = 1, klev
	1935	DO i = 1, klon
	1936	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
	1937	. *(paprs(i,k)-paprs(i,k+1))/RG
	1938	ENDDO
	1939	ENDDO
	1940	ENDIF
	1941	IF (iflag_con.EQ.1) THEN
	1942	stop'reactiver le call conlmd dans physiq.F'
	1943	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
	1944	c . d_t_con, d_q_con,
	1945	c . rain_con, snow_con, ibas_con, itop_con)
	1946	ELSE IF (iflag_con.EQ.2) THEN
	1947	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
[782]	1948	e conv_t, conv_q, -evap, omega,
[524]	1949	s d_t_con, d_q_con, rain_con, snow_con,
	1950	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	1951	s kcbot, kctop, kdtop, pmflxr, pmflxs)
	1952	WHERE (rain_con < 0.) rain_con = 0.
	1953	WHERE (snow_con < 0.) snow_con = 0.
	1954	DO i = 1, klon
	1955	ibas_con(i) = klev+1 - kcbot(i)
	1956	itop_con(i) = klev+1 - kctop(i)
	1957	ENDDO
	1958	ELSE IF (iflag_con.GE.3) THEN
	1959	c nb of tracers for the KE convection:
[619]	1960	c MAF la partie traceurs est faite dans phytrac
	1961	c on met ntra=1 pour limiter les appels mais on peut
	1962	c supprimer les calculs / ftra.
	1963	ntra = 1
[879]	1964
	1965	c=====================================================================================
	1966	cajout pour la parametrisation des poches froides:
	1967	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
	1968	do k=1,klev
	1969	do i=1,klon
	1970	if (iflag_wake.eq.1) then
	1971	t_wake(i,k) = t_seri(i,k)
	1972	. +(1-wake_s(i))*wake_deltat(i,k)
	1973	q_wake(i,k) = q_seri(i,k)
	1974	. +(1-wake_s(i))*wake_deltaq(i,k)
	1975	t_undi(i,k) = t_seri(i,k)
	1976	. -wake_s(i)*wake_deltat(i,k)
	1977	q_undi(i,k) = q_seri(i,k)
	1978	. -wake_s(i)*wake_deltaq(i,k)
	1979	else
	1980	t_wake(i,k) = t_seri(i,k)
	1981	q_wake(i,k) = q_seri(i,k)
	1982	t_undi(i,k) = t_seri(i,k)
	1983	q_undi(i,k) = q_seri(i,k)
	1984	endif
	1985	enddo
	1986	enddo
	1987
	1988	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
	1989	cc-- pour le soulevement des particules dans le modele convectif
	1990	c
	1991	do i = 1,klon
	1992	ALE(i) = 0.
	1993	ALP(i) = 0.
	1994	enddo
	1995	c
	1996	ccalcul de ale_wake et alp_wake
	1997	do i = 1,klon
	1998	if (iflag_wake.eq.1) then
	1999	ale_wake(i) = 0.5wake_cstar(i)*2
	2000	alp_wake(i) = wake_fip(i)
	2001	else
	2002	ale_wake(i) = 0.
	2003	alp_wake(i) = 0.
	2004	endif
	2005	enddo
	2006	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
	2007	cdans le thermique sinon
	2008	if (iflag_coupl.eq.0) then
[942]	2009	if (debut) print*,'ALE et ALP imposes'
[879]	2010	do i = 1,klon
	2011	con ne couple que ale
	2012	c ALE(i) = max(ale_wake(i),Ale_bl(i))
	2013	ALE(i) = max(ale_wake(i),ale_bl_prescr)
	2014	con ne couple que alp
	2015	c ALP(i) = alp_wake(i) + Alp_bl(i)
	2016	ALP(i) = alp_wake(i) + alp_bl_prescr
	2017	enddo
	2018	else
	2019	print*,'ALE et ALP couples au thermique'
	2020	do i = 1,klon
	2021	ALE(i) = max(ale_wake(i),Ale_bl(i))
	2022	ALP(i) = alp_wake(i) + Alp_bl(i)
	2023	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
	2024	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
	2025	enddo
	2026	endif
	2027
	2028	cfin calcul ale et alp
	2029	c=================================================================================================
	2030
	2031
[524]	2032	c sb, oct02:
	2033	c Schema de convection modularise et vectorise:
	2034	c (driver commun aux versions 3 et 4)
	2035	c
	2036	IF (ok_cvl) THEN ! new driver for convectL
	2037
[879]	2038	CALL concvl (iflag_con,iflag_clos,
	2039	. dtime,paprs,pplay,t_undi,q_undi,
	2040	. t_wake,q_wake,
	2041	. u_seri,v_seri,tr_seri,nbtr,
	2042	. ALE,ALP,
[524]	2043	. ema_work1,ema_work2,
	2044	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
[879]	2045	. rain_con, snow_con, ibas_con, itop_con, sigd,
[524]	2046	. upwd,dnwd,dnwd0,
[879]	2047	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
[619]	2048	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
[879]	2049	. pmflxr,pmflxs,da,phi,mp,
	2050	. ftd,fqd,lalim_conv,wght_th)
[619]	2051
[524]	2052	cIM cf. FH
	2053	clwcon0=qcondc
[619]	2054	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
[524]	2055
	2056	ELSE ! ok_cvl
[619]	2057	c MAF conema3 ne contient pas les traceurs
[524]	2058	CALL conema3 (dtime,
	2059	. paprs,pplay,t_seri,q_seri,
[619]	2060	. u_seri,v_seri,tr_seri,ntra,
[524]	2061	. ema_work1,ema_work2,
	2062	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
	2063	. rain_con, snow_con, ibas_con, itop_con,
	2064	. upwd,dnwd,dnwd0,bas,top,
	2065	. Ma,cape,tvp,rflag,
	2066	. pbase
	2067	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
	2068	. ,clwcon0)
	2069
	2070	ENDIF ! ok_cvl
	2071
[766]	2072	c
	2073	c Correction precip
	2074	rain_con = rain_con * cvl_corr
	2075	snow_con = snow_con * cvl_corr
	2076	c
	2077
[524]	2078	IF (.NOT. ok_gust) THEN
	2079	do i = 1, klon
	2080	wd(i)=0.0
	2081	enddo
	2082	ENDIF
	2083
	2084	c =================================================================== c
	2085	c Calcul des proprietes des nuages convectifs
	2086	c
	2087
	2088	c calcul des proprietes des nuages convectifs
	2089	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
	2090	call clouds_gno
	2091	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
	2092
	2093	c =================================================================== c
	2094
	2095	DO i = 1, klon
	2096	ema_pcb(i) = pbase(i)
	2097	ENDDO
	2098	DO i = 1, klon
[879]	2099
	2100	! L'idicage de itop_con peut cacher un pb potentiel
	2101	! FH sous la dictee de JYG, CR
	2102	ema_pct(i) = paprs(i,itop_con(i)+1)
	2103
[878]	2104	if (itop_con(i).gt.klev-3) then
	2105	print*,'La convection monte trop haut '
	2106	print*,'itop_con(,',i,',)=',itop_con(i)
	2107	endif
[524]	2108	ENDDO
	2109	DO i = 1, klon
	2110	ema_cbmf(i) = ema_workcbmf(i)
	2111	ENDDO
[881]	2112	ELSE IF (iflag_con.eq.0) THEN
	2113	write(lunout,*) 'On n appelle pas la convection'
	2114	clwcon0=0.
	2115	rnebcon0=0.
	2116	d_t_con=0.
	2117	d_q_con=0.
	2118	d_u_con=0.
	2119	d_v_con=0.
	2120	rain_con=0.
	2121	snow_con=0.
	2122	bas=1
	2123	top=1
[524]	2124	ELSE
	2125	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
	2126	CALL abort
	2127	ENDIF
	2128
	2129	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
	2130	c . d_u_con, d_v_con)
	2131
[904]	2132	!-----------------------------------------------------------------------------------------
	2133	! ajout des tendances de la diffusion turbulente
	2134	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
	2135	!-----------------------------------------------------------------------------------------
[766]	2136
	2137	if (mydebug) then
	2138	call writefield_phy('u_seri',u_seri,llm)
	2139	call writefield_phy('v_seri',v_seri,llm)
	2140	call writefield_phy('t_seri',t_seri,llm)
	2141	call writefield_phy('q_seri',q_seri,llm)
	2142	endif
	2143
[687]	2144	cIM
	2145	IF (ip_ebil_phy.ge.2) THEN
[524]	2146	ztit='after convect'
[687]	2147	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2148	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2149	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	2150	call diagphy(airephy,ztit,ip_ebil_phy
[524]	2151	e , zero_v, zero_v, zero_v, zero_v, zero_v
	2152	e , zero_v, rain_con, snow_con, ztsol
	2153	e , d_h_vcol, d_qt, d_ec
	2154	s , fs_bound, fq_bound )
	2155	END IF
	2156	C
	2157	IF (check) THEN
	2158	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2159	WRITE(lunout,*)"aprescon=", za
	2160	zx_t = 0.0
	2161	za = 0.0
	2162	DO i = 1, klon
	2163	za = za + airephy(i)/FLOAT(klon)
	2164	zx_t = zx_t + (rain_con(i)+
	2165	. snow_con(i))*airephy(i)/FLOAT(klon)
	2166	ENDDO
	2167	zx_t = zx_t/za*dtime
	2168	WRITE(lunout,*)"Precip=", zx_t
	2169	ENDIF
	2170	IF (zx_ajustq) THEN
	2171	DO i = 1, klon
	2172	z_apres(i) = 0.0
	2173	ENDDO
	2174	DO k = 1, klev
	2175	DO i = 1, klon
	2176	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
	2177	. *(paprs(i,k)-paprs(i,k+1))/RG
	2178	ENDDO
	2179	ENDDO
	2180	DO i = 1, klon
	2181	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
	2182	. /z_apres(i)
	2183	ENDDO
	2184	DO k = 1, klev
	2185	DO i = 1, klon
	2186	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
	2187	. z_factor(i).LT.(1.0-1.0E-08)) THEN
	2188	q_seri(i,k) = q_seri(i,k) * z_factor(i)
	2189	ENDIF
	2190	ENDDO
	2191	ENDDO
	2192	ENDIF
	2193	zx_ajustq=.FALSE.
[879]	2194
[524]	2195	c
[879]	2196	c=============================================================================
	2197	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
	2198	cpour la couche limite diffuse pour l instant
	2199	c
	2200	if (iflag_wake.eq.1) then
	2201	DO k=1,klev
	2202	DO i=1,klon
	2203	dt_dwn(i,k) = ftd(i,k)
	2204	wdt_PBL(i,k) = 0.
	2205	dq_dwn(i,k) = fqd(i,k)
	2206	wdq_PBL(i,k) = 0.
	2207	M_dwn(i,k) = dnwd0(i,k)
	2208	M_up(i,k) = upwd(i,k)
	2209	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
	2210	udt_PBL(i,k) = 0.
	2211	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
	2212	udq_PBL(i,k) = 0.
	2213	ENDDO
	2214	ENDDO
	2215	c
	2216	ccalcul caracteristiques de la poche froide
	2217	call calWAKE (paprs,pplay,dtime
[953]	2218	: ,t_seri,q_seri,omega
[879]	2219	: ,dt_dwn,dq_dwn,M_dwn,M_up
	2220	: ,dt_a,dq_a,sigd
	2221	: ,wdt_PBL,wdq_PBL
	2222	: ,udt_PBL,udq_PBL
	2223	o ,wake_deltat,wake_deltaq,wake_dth
	2224	o ,wake_h,wake_s,wake_dens
	2225	o ,wake_pe,wake_fip,wake_gfl
	2226	o ,dt_wake,dq_wake
	2227	o ,wake_k, t_undi,q_undi
	2228	o ,wake_omgbdth,wake_dp_omgb
	2229	o ,wake_dtKE,wake_dqKE
	2230	o ,wake_dtPBL,wake_dqPBL
	2231	o ,wake_omg,wake_dp_deltomg
	2232	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
	2233	o ,wake_ddeltat,wake_ddeltaq)
	2234	c
[904]	2235	!-----------------------------------------------------------------------------------------
	2236	! ajout des tendances des poches froides
	2237	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
	2238	! coherent avec les autres d_t_...
	2239	d_t_wake(:,:)=dt_wake(:,:)*dtime
	2240	d_q_wake(:,:)=dq_wake(:,:)*dtime
	2241	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
	2242	!-----------------------------------------------------------------------------------------
[879]	2243
	2244	endif
	2245	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
	2246	c
[541]	2247	c===================================================================
	2248	c Convection seche (thermiques ou ajustement)
	2249	c===================================================================
[524]	2250	c
[878]	2251	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
	2252	s ,seuil_inversion,weak_inversion,dthmin)
	2253
	2254
	2255
	2256	d_t_ajsb(:,:)=0.
	2257	d_q_ajsb(:,:)=0.
[541]	2258	d_t_ajs(:,:)=0.
	2259	d_u_ajs(:,:)=0.
	2260	d_v_ajs(:,:)=0.
	2261	d_q_ajs(:,:)=0.
[878]	2262	clwcon0th(:,:)=0.
[541]	2263	c
[557]	2264	IF(prt_level>9)WRITE(lunout,*)
	2265	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
[541]	2266	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
	2267	if(iflag_thermals.lt.0) then
	2268	c Rien
	2269	c ====
[557]	2270	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
[541]	2271
[878]	2272
[541]	2273	else
[878]	2274
[541]	2275	c Thermiques
	2276	c ==========
[557]	2277	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
[541]	2278	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
[878]	2279
	2280
	2281	if (iflag_thermals.gt.1) then
[541]	2282	call calltherm(pdtphys
[878]	2283	s ,pplay,paprs,pphi,weak_inversion
	2284	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
[541]	2285	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
[879]	2286	s ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth
	2287	s ,ratqsdiff,zqsatth
	2288	con rajoute ale et alp, et les caracteristiques de la couche alim
[927]	2289	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0)
[878]	2290	endif
	2291
	2292	! call calltherm(pdtphys
	2293	! s ,pplay,paprs,pphi
	2294	! s ,u_seri,v_seri,t_seri,q_seri
	2295	! s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
	2296	! s ,fm_therm,entr_therm)
	2297
	2298	c Ajustement sec
	2299	c ==============
	2300
	2301	! Dans le cas où on active les thermiques, on fait partir l'ajustement
	2302	! a partir du sommet des thermiques.
	2303	! Dans le cas contraire, on demarre au niveau 1.
	2304
	2305	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
	2306
	2307	if(iflag_thermals.eq.0) then
	2308	IF(prt_level>9)WRITE(lunout,*)'ajsec'
	2309	limbas(:)=1
	2310	else
	2311	limbas(:)=lmax_th(:)
	2312	endif
	2313
	2314	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
	2315	! pour des test de convergence numerique.
	2316	! Le nouveau ajsec est a priori mieux, meme pour le cas
	2317	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
	2318	! non nulles numeriquement pour des mailles non concernees.
	2319
	2320	if (iflag_thermals.eq.0) then
	2321	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
	2322	s , d_t_ajsb, d_q_ajsb)
	2323	else
	2324	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
	2325	s , d_t_ajsb, d_q_ajsb)
	2326	endif
	2327
[904]	2328	!-----------------------------------------------------------------------------------------
	2329	! ajout des tendances de l'ajustement sec ou des thermiques
	2330	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
[878]	2331	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
	2332	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
	2333
[904]	2334	!-----------------------------------------------------------------------------------------
	2335
[878]	2336	endif
	2337
[541]	2338	endif
	2339	c
	2340	c===================================================================
[687]	2341	cIM
	2342	IF (ip_ebil_phy.ge.2) THEN
[524]	2343	ztit='after dry_adjust'
[687]	2344	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2345	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2346	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2347	END IF
	2348
	2349
	2350	c-------------------------------------------------------------------------
	2351	c Caclul des ratqs
	2352	c-------------------------------------------------------------------------
	2353
	2354	c print*,'calcul des ratqs'
	2355	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
	2356	c ----------------
	2357	c on ecrase le tableau ratqsc calcule par clouds_gno
	2358	if (iflag_cldcon.eq.1) then
	2359	do k=1,klev
	2360	do i=1,klon
	2361	if(ptconv(i,k)) then
	2362	ratqsc(i,k)=ratqsbas
	2363	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
	2364	else
	2365	ratqsc(i,k)=0.
	2366	endif
	2367	enddo
	2368	enddo
[878]	2369	else if (iflag_cldcon.eq.4) then
	2370	ptconvth(:,:)=.false.
	2371	ratqsc(:,:)=0.
[942]	2372	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
[878]	2373	call clouds_gno
	2374	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
[942]	2375	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
[524]	2376	endif
	2377
	2378	c ratqs stables
	2379	c -------------
	2380
[878]	2381	if (iflag_ratqs.eq.0) then
[524]	2382
[878]	2383	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
	2384	do k=1,klev
	2385	do i=1, klon
	2386	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
	2387	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
	2388	enddo
	2389	enddo
	2390
	2391	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
	2392	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
	2393	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
	2394	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
	2395	! Il s'agit de differents tests dans la phase de reglage du modele
	2396	! avec thermiques.
	2397
	2398	else if (iflag_ratqs.eq.1) then
	2399
	2400	do k=1,klev
	2401	do i=1, klon
	2402	if (pplay(i,k).ge.60000.) then
	2403	ratqss(i,k)=ratqsbas
	2404	else if ((pplay(i,k).ge.30000.).and.
	2405	s (pplay(i,k).lt.60000.)) then
	2406	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
	2407	s (60000.-pplay(i,k))/(60000.-30000.)
	2408	else
	2409	ratqss(i,k)=ratqshaut
	2410	endif
	2411	enddo
	2412	enddo
	2413
	2414	else
	2415
	2416	do k=1,klev
	2417	do i=1, klon
	2418	if (pplay(i,k).ge.60000.) then
	2419	ratqss(i,k)=ratqsbas
	2420	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
	2421	else if ((pplay(i,k).ge.30000.).and.
	2422	s (pplay(i,k).lt.60000.)) then
	2423	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
	2424	s (60000.-pplay(i,k))/(60000.-30000.)
	2425	else
	2426	ratqss(i,k)=ratqshaut
	2427	endif
	2428	enddo
	2429	enddo
	2430	endif
	2431
	2432
	2433
	2434
[524]	2435	c ratqs final
	2436	c -----------
[878]	2437
	2438	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
	2439	s .or.iflag_cldcon.eq.4) then
	2440
	2441	! On ajoute une constante au ratqsc*2 pour tenir compte de
	2442	! fluctuations turbulentes de petite echelle
	2443
	2444	do k=1,klev
	2445	do i=1,klon
	2446	if ((fm_therm(i,k).gt.1.e-10)) then
	2447	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
	2448	endif
	2449	enddo
	2450	enddo
	2451
	2452	! les ratqs sont une conbinaison de ratqss et ratqsc
	2453	! 1800s, en dur pour le moment, est le temps de
	2454	! relaxation des ratqs
	2455
	2456	facteur=exp(-pdtphys/1800.)
	2457
	2458	print*,'WARNING ratqs a revoir '
	2459	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
	2460	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
	2461
[524]	2462	else
[878]	2463	! on ne prend que le ratqs stable pour fisrtilp
[524]	2464	ratqs(:,:)=ratqss(:,:)
	2465	endif
	2466
	2467
	2468	c
	2469	c Appeler le processus de condensation a grande echelle
	2470	c et le processus de precipitation
	2471	c-------------------------------------------------------------------------
	2472	CALL fisrtilp(dtime,paprs,pplay,
	2473	. t_seri, q_seri,ptconv,ratqs,
	2474	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
	2475	. rain_lsc, snow_lsc,
	2476	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
	2477	. frac_impa, frac_nucl,
	2478	. prfl, psfl, rhcl)
	2479
	2480	WHERE (rain_lsc < 0) rain_lsc = 0.
	2481	WHERE (snow_lsc < 0) snow_lsc = 0.
[904]	2482	!-----------------------------------------------------------------------------------------
	2483	! ajout des tendances de la diffusion turbulente
	2484	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
	2485	!-----------------------------------------------------------------------------------------
[524]	2486	DO k = 1, klev
	2487	DO i = 1, klon
	2488	cldfra(i,k) = rneb(i,k)
	2489	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
	2490	ENDDO
	2491	ENDDO
	2492	IF (check) THEN
	2493	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2494	WRITE(lunout,*)"apresilp=", za
	2495	zx_t = 0.0
	2496	za = 0.0
	2497	DO i = 1, klon
	2498	za = za + airephy(i)/FLOAT(klon)
	2499	zx_t = zx_t + (rain_lsc(i)
	2500	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
	2501	ENDDO
	2502	zx_t = zx_t/za*dtime
	2503	WRITE(lunout,*)"Precip=", zx_t
	2504	ENDIF
[687]	2505	cIM
	2506	IF (ip_ebil_phy.ge.2) THEN
[524]	2507	ztit='after fisrt'
[687]	2508	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2509	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2510	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	2511	call diagphy(airephy,ztit,ip_ebil_phy
[524]	2512	e , zero_v, zero_v, zero_v, zero_v, zero_v
	2513	e , zero_v, rain_lsc, snow_lsc, ztsol
	2514	e , d_h_vcol, d_qt, d_ec
	2515	s , fs_bound, fq_bound )
	2516	END IF
[766]	2517
	2518	if (mydebug) then
	2519	call writefield_phy('u_seri',u_seri,llm)
	2520	call writefield_phy('v_seri',v_seri,llm)
	2521	call writefield_phy('t_seri',t_seri,llm)
	2522	call writefield_phy('q_seri',q_seri,llm)
	2523	endif
	2524
[524]	2525	c
	2526	c-------------------------------------------------------------------
	2527	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
	2528	c-------------------------------------------------------------------
	2529
	2530	c 1. NUAGES CONVECTIFS
	2531	c
[644]	2532	cIM cf FH
	2533	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
[878]	2534	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
[644]	2535	snow_tiedtke=0.
	2536	c print*,'avant calcul de la pseudo precip '
	2537	c print*,'iflag_cldcon',iflag_cldcon
	2538	if (iflag_cldcon.eq.-1) then
	2539	rain_tiedtke=rain_con
	2540	else
	2541	c print*,'calcul de la pseudo precip '
	2542	rain_tiedtke=0.
	2543	c print*,'calcul de la pseudo precip 0'
	2544	do k=1,klev
	2545	do i=1,klon
	2546	if (d_q_con(i,k).lt.0.) then
	2547	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
	2548	s *(paprs(i,k)-paprs(i,k+1))/rg
	2549	endif
	2550	enddo
	2551	enddo
	2552	endif
	2553	c
	2554	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
	2555	c
[524]	2556
	2557	c Nuages diagnostiques pour Tiedtke
	2558	CALL diagcld1(paprs,pplay,
[644]	2559	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
	2560	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
[524]	2561	. diafra,dialiq)
	2562	DO k = 1, klev
	2563	DO i = 1, klon
	2564	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	2565	cldliq(i,k) = dialiq(i,k)
	2566	cldfra(i,k) = diafra(i,k)
	2567	ENDIF
	2568	ENDDO
	2569	ENDDO
	2570
[878]	2571	ELSE IF (iflag_cldcon.ge.3) THEN
[524]	2572	c On prend pour les nuages convectifs le max du calcul de la
[766]	2573	c convection et du calcul du pas de temps precedent diminue d'un facteur
[524]	2574	c facttemps
	2575	facteur = pdtphys *facttemps
	2576	do k=1,klev
	2577	do i=1,klon
	2578	rnebcon(i,k)=rnebcon(i,k)*facteur
	2579	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
	2580	s then
	2581	rnebcon(i,k)=rnebcon0(i,k)
	2582	clwcon(i,k)=clwcon0(i,k)
	2583	endif
	2584	enddo
	2585	enddo
	2586
[644]	2587	c
[766]	2588	cjq - introduce the aerosol direct and first indirect radiative forcings
	2589	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
	2590	IF (ok_ade.OR.ok_aie) THEN
[955]	2591	IF ( .NOT. aerosol_couple ) THEN
	2592	! Get sulfate aerosol distribution
	2593	CALL readsulfate(rjourvrai, debut, sulfate)
	2594	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
[766]	2595
[955]	2596	! Calculate aerosol optical properties (Olivier Boucher)
	2597	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
	2598	. tau_ae, piz_ae, cg_ae, aerindex)
	2599	ENDIF
[766]	2600	ELSE
	2601	tau_ae(:,:,:)=0.0
	2602	piz_ae(:,:,:)=0.0
	2603	cg_ae(:,:,:)=0.0
	2604	ENDIF
	2605
	2606	c
[524]	2607	cIM calcul nuages par le simulateur ISCCP
[644]	2608	c
[839]	2609	#ifdef histISCCP
[524]	2610	IF (ok_isccp) THEN
[828]	2611	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
	2612	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
[644]	2613	#include "calcul_simulISCCP.h"
[828]	2614	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
[524]	2615	ENDIF !ok_isccp
[839]	2616	#endif
[524]	2617
	2618	c On prend la somme des fractions nuageuses et des contenus en eau
	2619	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	2620	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
	2621
	2622	ENDIF
	2623
	2624	c
	2625	c 2. NUAGES STARTIFORMES
	2626	c
	2627	IF (ok_stratus) THEN
	2628	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
	2629	DO k = 1, klev
	2630	DO i = 1, klon
	2631	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	2632	cldliq(i,k) = dialiq(i,k)
	2633	cldfra(i,k) = diafra(i,k)
	2634	ENDIF
	2635	ENDDO
	2636	ENDDO
	2637	ENDIF
	2638	c
	2639	c Precipitation totale
	2640	c
	2641	DO i = 1, klon
	2642	rain_fall(i) = rain_con(i) + rain_lsc(i)
	2643	snow_fall(i) = snow_con(i) + snow_lsc(i)
	2644	ENDDO
[687]	2645	cIM
	2646	IF (ip_ebil_phy.ge.2) THEN
[524]	2647	ztit="after diagcld"
[687]	2648	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2649	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2650	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2651	END IF
	2652	c
	2653	c Calculer l'humidite relative pour diagnostique
	2654	c
	2655	DO k = 1, klev
	2656	DO i = 1, klon
	2657	zx_t = t_seri(i,k)
	2658	IF (thermcep) THEN
	2659	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	2660	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	2661	zx_qs = MIN(0.5,zx_qs)
	2662	zcor = 1./(1.-retv*zx_qs)
	2663	zx_qs = zx_qs*zcor
	2664	ELSE
	2665	IF (zx_t.LT.t_coup) THEN
	2666	zx_qs = qsats(zx_t)/pplay(i,k)
	2667	ELSE
	2668	zx_qs = qsatl(zx_t)/pplay(i,k)
	2669	ENDIF
	2670	ENDIF
	2671	zx_rh(i,k) = q_seri(i,k)/zx_qs
	2672	zqsat(i,k)=zx_qs
	2673	ENDDO
	2674	ENDDO
[782]	2675
[687]	2676	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
	2677	c equivalente a 2m (tpote) pour diagnostique
	2678	c
	2679	DO i = 1, klon
	2680	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
	2681	IF (thermcep) THEN
	2682	IF(zt2m(i).LT.RTT) then
	2683	Lheat=RLSTT
	2684	ELSE
	2685	Lheat=RLVTT
	2686	ENDIF
	2687	ELSE
	2688	IF (zt2m(i).LT.RTT) THEN
	2689	Lheat=RLSTT
	2690	ELSE
	2691	Lheat=RLVTT
	2692	ENDIF
	2693	ENDIF
	2694	tpote(i) = tpot(i)*
	2695	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
	2696	ENDDO
[524]	2697
	2698
	2699	#ifdef INCA
[766]	2700	call VTe(VTphysiq)
	2701	call VTb(VTinca)
[524]	2702	calday = FLOAT(julien) + gmtime
	2703
	2704	#ifdef INCA_AER
[766]	2705	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs
	2706	& ,prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
[524]	2707	#endif
	2708
	2709	#ifdef INCAINFO
	2710	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
	2711	#endif
[625]	2712
[782]	2713	zxsnow_dummy(:) = 0.0
	2714
[524]	2715	CALL chemhook_begin (calday,
[678]	2716	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
	2717	$ julien,
	2718	$ gmtime,
	2719	#endif
[593]	2720	$ pctsrf(1,1),
[524]	2721	$ rlat,
	2722	$ rlon,
	2723	$ airephy,
	2724	$ paprs,
	2725	$ pplay,
	2726	$ ycoefh,
	2727	$ pphi,
	2728	$ t_seri,
	2729	$ u,
	2730	$ v,
	2731	$ wo,
	2732	$ q_seri,
	2733	$ zxtsol,
[782]	2734	$ zxsnow_dummy,
[524]	2735	$ solsw,
[888]	2736	$ albsol1,
[524]	2737	$ rain_fall,
	2738	$ snow_fall,
	2739	$ itop_con,
	2740	$ ibas_con,
	2741	$ cldfra,
	2742	$ iim,
	2743	$ jjm,
[616]	2744	#ifdef INCA_AER
	2745	$ tr_seri,
	2746	$ ftsol,
	2747	$ paprs,
	2748	$ cdragh,
	2749	$ cdragm,
	2750	$ pctsrf,
	2751	$ pdtphys,
	2752	$ itap)
	2753	#else
	2754	$ tr_seri)
	2755	#endif
	2756
[625]	2757
[524]	2758	#ifdef INCAINFO
	2759	WRITE(lunout,*)'OK.'
	2760	#endif
[766]	2761	call VTe(VTinca)
	2762	call VTb(VTphysiq)
[524]	2763	#endif
	2764	c
	2765	c Calculer les parametres optiques des nuages et quelques
	2766	c parametres pour diagnostiques:
	2767	c
[955]	2768	IF (aerosol_couple ) THEN
	2769	#ifdef INCA
	2770	sulfate(:,:) = ccm(:,:,1)
	2771	sulfate_pi(:,:) = ccm(:,:,2)
	2772	#endif
	2773	ENDIF
	2774
[524]	2775	if (ok_newmicro) then
	2776	CALL newmicro (paprs, pplay,ok_newmicro,
	2777	. t_seri, cldliq, cldfra, cldtau, cldemi,
	2778	. cldh, cldl, cldm, cldt, cldq,
	2779	. flwp, fiwp, flwc, fiwc,
	2780	e ok_aie,
	2781	e sulfate, sulfate_pi,
	2782	e bl95_b0, bl95_b1,
	2783	s cldtaupi, re, fl)
	2784	else
	2785	CALL nuage (paprs, pplay,
	2786	. t_seri, cldliq, cldfra, cldtau, cldemi,
	2787	. cldh, cldl, cldm, cldt, cldq,
	2788	e ok_aie,
	2789	e sulfate, sulfate_pi,
	2790	e bl95_b0, bl95_b1,
	2791	s cldtaupi, re, fl)
	2792
	2793	endif
	2794	c
	2795	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
	2796	c
	2797	IF (MOD(itaprad,radpas).EQ.0) THEN
[782]	2798
[524]	2799	DO i = 1, klon
[888]	2800	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
	2801	. + falb1(i,is_lic) * pctsrf(i,is_lic)
	2802	. + falb1(i,is_ter) * pctsrf(i,is_ter)
	2803	. + falb1(i,is_sic) * pctsrf(i,is_sic)
	2804	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
	2805	. + falb2(i,is_lic) * pctsrf(i,is_lic)
	2806	. + falb2(i,is_ter) * pctsrf(i,is_ter)
	2807	. + falb2(i,is_sic) * pctsrf(i,is_sic)
[524]	2808	ENDDO
[766]	2809
	2810	if (mydebug) then
	2811	call writefield_phy('u_seri',u_seri,llm)
	2812	call writefield_phy('v_seri',v_seri,llm)
	2813	call writefield_phy('t_seri',t_seri,llm)
	2814	call writefield_phy('q_seri',q_seri,llm)
	2815	endif
	2816
[955]	2817	IF (aerosol_couple) THEN
	2818	#ifdef INCA_AER
	2819	CALL radlwsw_inca
	2820	e (kdlon,kflev,dist, rmu0, fract, solaire,
[864]	2821	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
	2822	e wo,
	2823	e cldfra, cldemi, cldtau,
	2824	s heat,heat0,cool,cool0,radsol,albpla,
	2825	s topsw,toplw,solsw,sollw,
	2826	s sollwdown,
	2827	s topsw0,toplw0,solsw0,sollw0,
	2828	s lwdn0, lwdn, lwup0, lwup,
	2829	s swdn0, swdn, swup0, swup,
[955]	2830	e ok_ade, ok_aie,
	2831	e tau_inca, piz_inca, cg_inca,
	2832	s topswad_inca, solswad_inca,
	2833	s topswad0_inca, solswad0_inca,
[864]	2834	s topsw_inca, topsw0_inca,
	2835	s solsw_inca, solsw0_inca,
[955]	2836	e cldtaupi,
	2837	s topswai_inca, solswai_inca)
	2838	#endif
	2839	ELSE
	2840	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
[524]	2841	e (dist, rmu0, fract,
[888]	2842	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
[524]	2843	e wo,
	2844	e cldfra, cldemi, cldtau,
	2845	s heat,heat0,cool,cool0,radsol,albpla,
	2846	s topsw,toplw,solsw,sollw,
	2847	s sollwdown,
	2848	s topsw0,toplw0,solsw0,sollw0,
	2849	s lwdn0, lwdn, lwup0, lwup,
	2850	s swdn0, swdn, swup0, swup,
	2851	e ok_ade, ok_aie, ! new for aerosol radiative effects
	2852	e tau_ae, piz_ae, cg_ae, ! ="=
	2853	s topswad, solswad, ! ="=
	2854	e cldtaupi, ! ="=
	2855	s topswai, solswai) ! ="=
[955]	2856	ENDIF
	2857
[524]	2858	itaprad = 0
	2859	ENDIF
	2860	itaprad = itaprad + 1
[879]	2861
	2862	if (iflag_radia.eq.0) then
	2863	print *,'--------------------------------------------------'
	2864	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
	2865	print *,'>>>> heat et cool mis a zero '
	2866	print *,'--------------------------------------------------'
	2867	heat=0.
	2868	cool=0.
	2869	endif
	2870
	2871
[524]	2872	c
	2873	c Ajouter la tendance des rayonnements (tous les pas)
	2874	c
	2875	DO k = 1, klev
	2876	DO i = 1, klon
	2877	t_seri(i,k) = t_seri(i,k)
	2878	. + (heat(i,k)-cool(i,k)) * dtime/86400.
	2879	ENDDO
	2880	ENDDO
[766]	2881	c
	2882	if (mydebug) then
	2883	call writefield_phy('u_seri',u_seri,llm)
	2884	call writefield_phy('v_seri',v_seri,llm)
	2885	call writefield_phy('t_seri',t_seri,llm)
	2886	call writefield_phy('q_seri',q_seri,llm)
	2887	endif
	2888
[687]	2889	cIM
	2890	IF (ip_ebil_phy.ge.2) THEN
[524]	2891	ztit='after rad'
[687]	2892	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2893	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2894	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	2895	call diagphy(airephy,ztit,ip_ebil_phy
[524]	2896	e , topsw, toplw, solsw, sollw, zero_v
	2897	e , zero_v, zero_v, zero_v, ztsol
	2898	e , d_h_vcol, d_qt, d_ec
	2899	s , fs_bound, fq_bound )
	2900	END IF
	2901	c
	2902	c
	2903	c Calculer l'hydrologie de la surface
	2904	c
	2905	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
	2906	c . agesno, ftsol,fqsurf,fsnow, ruis)
	2907	c
[782]	2908
[524]	2909	c
	2910	c Calculer le bilan du sol et la derive de temperature (couplage)
	2911	c
	2912	DO i = 1, klon
	2913	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
	2914	c a la demande de JLD
	2915	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
	2916	ENDDO
	2917	c
	2918	cmoddeblott(jan95)
	2919	c Appeler le programme de parametrisation de l'orographie
	2920	c a l'echelle sous-maille:
	2921	c
	2922	IF (ok_orodr) THEN
	2923	c
	2924	c selection des points pour lesquels le shema est actif:
	2925	igwd=0
	2926	DO i=1,klon
	2927	itest(i)=0
	2928	c IF ((zstd(i).gt.10.0)) THEN
	2929	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	2930	itest(i)=1
	2931	igwd=igwd+1
	2932	idx(igwd)=i
	2933	ENDIF
	2934	ENDDO
	2935	c igwdim=MAX(1,igwd)
	2936	c
	2937	CALL drag_noro(klon,klev,dtime,paprs,pplay,
	2938	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
	2939	e igwd,idx,itest,
	2940	e t_seri, u_seri, v_seri,
[644]	2941	cIM 141004 s zulow, zvlow, zustr, zvstr,
	2942	s zulow, zvlow, zustrdr, zvstrdr,
[524]	2943	s d_t_oro, d_u_oro, d_v_oro)
	2944	c
	2945	c ajout des tendances
[904]	2946	!-----------------------------------------------------------------------------------------
	2947	! ajout des tendances de la trainee de l'orographie
	2948	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
	2949	!-----------------------------------------------------------------------------------------
[524]	2950	c
	2951	ENDIF ! fin de test sur ok_orodr
	2952	c
[766]	2953	if (mydebug) then
	2954	call writefield_phy('u_seri',u_seri,llm)
	2955	call writefield_phy('v_seri',v_seri,llm)
	2956	call writefield_phy('t_seri',t_seri,llm)
	2957	call writefield_phy('q_seri',q_seri,llm)
	2958	endif
	2959
[524]	2960	IF (ok_orolf) THEN
	2961	c
	2962	c selection des points pour lesquels le shema est actif:
	2963	igwd=0
	2964	DO i=1,klon
	2965	itest(i)=0
	2966	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	2967	itest(i)=1
	2968	igwd=igwd+1
	2969	idx(igwd)=i
	2970	ENDIF
	2971	ENDDO
	2972	c igwdim=MAX(1,igwd)
	2973	c
	2974	CALL lift_noro(klon,klev,dtime,paprs,pplay,
	2975	e rlat,zmea,zstd,zpic,
	2976	e itest,
	2977	e t_seri, u_seri, v_seri,
[644]	2978	s zulow, zvlow, zustrli, zvstrli,
[524]	2979	s d_t_lif, d_u_lif, d_v_lif)
	2980	c
[904]	2981	!-----------------------------------------------------------------------------------------
	2982	! ajout des tendances de la portance de l'orographie
	2983	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
	2984	!-----------------------------------------------------------------------------------------
[524]	2985	c
	2986	ENDIF ! fin de test sur ok_orolf
	2987	c
[644]	2988	cIM cf. FLott BEG
	2989	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
	2990
[766]	2991	if (mydebug) then
	2992	call writefield_phy('u_seri',u_seri,llm)
	2993	call writefield_phy('v_seri',v_seri,llm)
	2994	call writefield_phy('t_seri',t_seri,llm)
	2995	call writefield_phy('q_seri',q_seri,llm)
	2996	endif
	2997
[644]	2998	DO i = 1, klon
	2999	zustrph(i)=0.
	3000	zvstrph(i)=0.
	3001	ENDDO
	3002	DO k = 1, klev
	3003	DO i = 1, klon
	3004	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
	3005	c (paprs(i,k)-paprs(i,k+1))/rg
	3006	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
	3007	c (paprs(i,k)-paprs(i,k+1))/rg
	3008	ENDDO
	3009	ENDDO
	3010	c
	3011	cIM calcul composantes axiales du moment angulaire et couple des montagnes
	3012	c
[776]	3013	IF (is_sequential) THEN
[766]	3014
	3015	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
	3016	C ra,rg,romega,
	3017	C rlat,rlon,pphis,
	3018	C zustrdr,zustrli,zustrph,
	3019	C zvstrdr,zvstrli,zvstrph,
	3020	C paprs,u,v,
	3021	C aam, torsfc)
	3022	ENDIF
[644]	3023	cIM cf. FLott END
[687]	3024	cIM
	3025	IF (ip_ebil_phy.ge.2) THEN
[524]	3026	ztit='after orography'
[687]	3027	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	3028	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	3029	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3030	END IF
	3031	c
	3032	c
	3033	cAA
	3034	cAA Installation de l'interface online-offline pour traceurs
	3035	cAA
	3036	c====================================================================
	3037	c Calcul des tendances traceurs
	3038	c====================================================================
	3039	C
[658]	3040	call phytrac ( rnpb,
[625]	3041	I itap,
	3042	I julien,
	3043	I gmtime,
	3044	I debut,
	3045	I lafin,
[524]	3046	I nqmax-2,
[625]	3047	I nlon,
	3048	I nlev,
	3049	I dtime,
	3050	I u,
	3051	I v,
	3052	I t,
	3053	I paprs,
	3054	I pplay,
	3055	I pmfu,
	3056	I pmfd,
	3057	I pen_u,
	3058	I pde_u,
	3059	I pen_d,
	3060	I pde_d,
	3061	I ycoefh,
	3062	I fm_therm,
	3063	I entr_therm,
	3064	I yu1,
	3065	I yv1,
	3066	I ftsol,
	3067	I pctsrf,
	3068	I rlat,
	3069	I frac_impa,
	3070	I frac_nucl,
	3071	I rlon,
	3072	I presnivs,
	3073	I pphis,
	3074	I pphi,
[888]	3075	I albsol1,
[625]	3076	I qx(1,1,1),
	3077	I rhcl,
	3078	I cldfra,
	3079	I rneb,
	3080	I diafra,
	3081	I cldliq,
	3082	I itop_con,
[524]	3083	I ibas_con,
[625]	3084	I pmflxr,
	3085	I pmflxs,
	3086	I prfl,
	3087	I psfl,
	3088	I da,
	3089	I phi,
	3090	I mp,
	3091	I upwd,
	3092	I dnwd,
[955]	3093	I aerosol_couple,
[625]	3094	#ifdef INCA
[524]	3095	I flxmass_w,
[864]	3096	I tau_inca,
	3097	I piz_inca,
	3098	I cg_inca,
	3099	I ccm,
	3100	I rfname,
[524]	3101	#endif
	3102	O tr_seri)
	3103
	3104	IF (offline) THEN
	3105
[541]	3106	print*,'Attention on met a 0 les thermiques pour phystoke'
[524]	3107	call phystokenc (
	3108	I nlon,nlev,pdtphys,rlon,rlat,
	3109	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
[541]	3110	I fm_therm,entr_therm,
[524]	3111	I ycoefh,yu1,yv1,ftsol,pctsrf,
	3112	I frac_impa, frac_nucl,
	3113	I pphis,airephy,dtime,itap)
	3114
	3115
	3116	ENDIF
	3117
	3118	c
	3119	c Calculer le transport de l'eau et de l'energie (diagnostique)
	3120	c
	3121	CALL transp (paprs,zxtsol,
	3122	e t_seri, q_seri, u_seri, v_seri, zphi,
	3123	s ve, vq, ue, uq)
	3124	c
[687]	3125	cIM global posePB BEG
	3126	IF(1.EQ.0) THEN
[524]	3127	c
[644]	3128	CALL transp_lay (paprs,zxtsol,
	3129	e t_seri, q_seri, u_seri, v_seri, zphi,
	3130	s ve_lay, vq_lay, ue_lay, uq_lay)
[524]	3131	c
[687]	3132	ENDIF !(1.EQ.0) THEN
	3133	cIM global posePB END
[644]	3134	c Accumuler les variables a stocker dans les fichiers histoire:
[524]	3135	c
	3136	c+jld ec_conser
	3137	DO k = 1, klev
	3138	DO i = 1, klon
	3139	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
	3140	d_t_ec(i,k)=0.5/ZRCPD
	3141	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
	3142	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
	3143	d_t_ec(i,k) = d_t_ec(i,k)/dtime
	3144	END DO
	3145	END DO
	3146	c-jld ec_conser
[687]	3147	cIM
	3148	IF (ip_ebil_phy.ge.1) THEN
[524]	3149	ztit='after physic'
[687]	3150	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
[524]	3151	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	3152	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3153	C Comme les tendances de la physique sont ajoute dans la dynamique,
	3154	C on devrait avoir que la variation d'entalpie par la dynamique
	3155	C est egale a la variation de la physique au pas de temps precedent.
	3156	C Donc la somme de ces 2 variations devrait etre nulle.
[687]	3157	call diagphy(airephy,ztit,ip_ebil_phy
[524]	3158	e , topsw, toplw, solsw, sollw, sens
	3159	e , evap, rain_fall, snow_fall, ztsol
	3160	e , d_h_vcol, d_qt, d_ec
	3161	s , fs_bound, fq_bound )
	3162	C
	3163	d_h_vcol_phy=d_h_vcol
	3164	C
	3165	END IF
	3166	C
	3167	c=======================================================================
	3168	c SORTIES
	3169	c=======================================================================
	3170
[644]	3171	cIM Interpolation sur les niveaux de pression du NMC
	3172	c -------------------------------------------------
[524]	3173	c
[644]	3174	#include "calcul_STDlev.h"
[524]	3175	c
	3176	c slp sea level pressure
	3177	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
	3178	c
	3179	ccc prw = eau precipitable
	3180	DO i = 1, klon
	3181	prw(i) = 0.
	3182	DO k = 1, klev
	3183	prw(i) = prw(i) +
	3184	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
	3185	ENDDO
	3186	ENDDO
	3187	c
[644]	3188	cIM initialisation + calculs divers diag AMIP2
[524]	3189	c
[644]	3190	#include "calcul_divers.h"
	3191	c
[655]	3192	#ifdef INCA
[766]	3193	call VTe(VTphysiq)
	3194	call VTb(VTinca)
[655]	3195	#ifdef INCAINFO
	3196	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
	3197	#endif
	3198	CALL chemhook_end (calday,
	3199	$ dtime,
	3200	$ pplay,
	3201	$ t_seri,
	3202	$ tr_seri,
	3203	$ nbtr,
	3204	$ paprs,
	3205	$ q_seri,
	3206	$ annee_ref,
	3207	$ day_ini,
[791]	3208	$ airephy,
[655]	3209	#ifdef INCA_AER
	3210	$ xjour,
	3211	$ pphi,
	3212	$ pphis,
[766]	3213	$ zx_rh)
[655]	3214	#else
	3215	$ xjour)
	3216	#endif
	3217	#ifdef INCAINFO
	3218	WRITE(lunout,*)'OK.'
	3219	#endif
[766]	3220	call VTe(VTinca)
	3221	call VTb(VTphysiq)
[655]	3222	#endif
	3223
[524]	3224	c=============================================================
	3225	c
	3226	c Convertir les incrementations en tendances
	3227	c
[766]	3228	if (mydebug) then
	3229	call writefield_phy('u_seri',u_seri,llm)
	3230	call writefield_phy('v_seri',v_seri,llm)
	3231	call writefield_phy('t_seri',t_seri,llm)
	3232	call writefield_phy('q_seri',q_seri,llm)
	3233	endif
	3234
[524]	3235	DO k = 1, klev
	3236	DO i = 1, klon
	3237	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	3238	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	3239	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
	3240	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
	3241	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
	3242	ENDDO
	3243	ENDDO
	3244	c
	3245	IF (nqmax.GE.3) THEN
	3246	DO iq = 3, nqmax
	3247	DO k = 1, klev
	3248	DO i = 1, klon
	3249	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
	3250	ENDDO
	3251	ENDDO
	3252	ENDDO
	3253	ENDIF
	3254	c
[644]	3255	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
[687]	3256	cIM global posePB#include "write_bilKP_ins.h"
	3257	cIM global posePB#include "write_bilKP_ave.h"
[644]	3258	c
[524]	3259	c Sauvegarder les valeurs de t et q a la fin de la physique:
	3260	c
	3261	DO k = 1, klev
	3262	DO i = 1, klon
	3263	t_ancien(i,k) = t_seri(i,k)
	3264	q_ancien(i,k) = q_seri(i,k)
	3265	ENDDO
	3266	ENDDO
	3267	c
[879]	3268	!==========================================================================
	3269	! Sorties des tendances pour un point particulier
	3270	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
	3271	! pour le debug
	3272	! La valeur de igout est attribuee plus haut dans le programme
	3273	!==========================================================================
	3274
[942]	3275	if (prt_level.ge.1) then
[879]	3276	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	3277	write(lunout,*)
[929]	3278	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
[879]	3279	write(lunout,*)
[929]	3280	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys,
[930]	3281	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
[929]	3282	s pctsrf(igout,is_sic)
[879]	3283	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
	3284	do k=1,nlev
	3285	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
	3286	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
	3287	s d_t_eva(igout,k)
	3288	enddo
	3289	write(lunout,*) 'cool,heat'
	3290	do k=1,nlev
	3291	write(lunout,*) cool(igout,k),heat(igout,k)
	3292	enddo
	3293
	3294	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
	3295	do k=1,nlev
	3296	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
	3297	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
	3298	enddo
	3299
	3300	write(lunout,*) 'd_ps ',d_ps(igout)
	3301	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
	3302	do k=1,nlev
	3303	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
	3304	s d_qx(igout,k,1),d_qx(igout,k,2)
	3305	enddo
	3306	endif
	3307
	3308	!==========================================================================
	3309
	3310	c============================================================
	3311	c Calcul de la temperature potentielle
	3312	c============================================================
	3313	DO k = 1, klev
	3314	DO i = 1, klon
	3315	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
	3316	ENDDO
	3317	ENDDO
	3318	c
	3319
[524]	3320	c 22.03.04 BEG
	3321	c=============================================================
	3322	c Ecriture des sorties
	3323	c=============================================================
	3324	#ifdef CPP_IOIPSL
[782]	3325
	3326	c Recupere des varibles calcule dans differents modules
	3327	c pour ecriture dans histxxx.nc
[524]	3328
[888]	3329	! Get some variables from module fonte_neige_mod
[782]	3330	CALL fonte_neige_get_vars(pctsrf,
	3331	. zxfqcalving, zxfqfonte, zxffonte)
	3332
	3333	IF (ocean == 'slab') THEN
[888]	3334	! Get some variables from module ocean_slab_mod
[782]	3335	CALL ocean_slab_get_vars(tslab, seaice, fluxo, fluxg)
	3336	ELSEIF (ocean == 'couple') THEN
[888]	3337	! Get some variables from module ocean_cpl_mod
[782]	3338	CALL ocean_cpl_get_vars(fluxo, fluxg)
	3339	ELSE
[888]	3340	! Get some variables from module ocean_forced_mod
[782]	3341	CALL ocean_forced_get_vars(fluxo, fluxg)
	3342	ENDIF
	3343
[524]	3344
[909]	3345	c Commente par abderrahmane le 11 2 08
	3346	c#ifdef histhf
	3347	c#include "write_histhf.h"
	3348	c#endif
[524]	3349
[909]	3350	c#ifdef histday
	3351	c#include "write_histday.h"
	3352	c#endif
[524]	3353
[909]	3354	c#ifdef histmth
	3355	c#include "write_histmth.h"
	3356	c#endif
[524]	3357
[909]	3358	c#ifdef histins
	3359	c#include "write_histins.h"
	3360	c#endif
	3361
	3362	#include "phys_output_write.h"
	3363
[524]	3364	#ifdef histISCCP
	3365	#include "write_histISCCP.h"
	3366	#endif
	3367
	3368	#ifdef histmthNMC
	3369	#include "write_histmthNMC.h"
	3370	#endif
	3371
[687]	3372	#include "write_histday_seri.h"
	3373
	3374	#include "write_paramLMDZ_phy.h"
	3375
[524]	3376	#endif
	3377
	3378	c 22.03.04 END
	3379	c
	3380	c====================================================================
	3381	c Si c'est la fin, il faut conserver l'etat de redemarrage
	3382	c====================================================================
	3383	c
[782]	3384
	3385
[524]	3386	IF (lafin) THEN
	3387	itau_phy = itau_phy + itap
[782]	3388	CALL phyredem ("restartphy.nc",dtime,radpas,ocean,
	3389	. rlat, rlon, pctsrf, ftsol,
[888]	3390	. falb1, falb2, rain_fall,
[782]	3391	. snow_fall,
	3392	. solsw, sollw,
	3393	. radsol,
[524]	3394	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
[878]	3395	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,
[937]	3396	. pbl_tke, zmax0, f0, ema_work1, ema_work2)
[904]	3397	open(97,form="unformatted",file="finbin")
	3398	write(97) u_seri,v_seri,t_seri,q_seri
	3399	close(97)
[524]	3400	ENDIF
	3401
	3402
	3403	RETURN
	3404	END
	3405	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
	3406	IMPLICIT none
	3407	c
	3408	c Calculer et imprimer l'eau totale. A utiliser pour verifier
	3409	c la conservation de l'eau
	3410	c
	3411	#include "YOMCST.h"
	3412	INTEGER klon,klev
	3413	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
	3414	REAL aire(klon)
	3415	REAL qtotal, zx, qcheck
	3416	INTEGER i, k
	3417	c
	3418	zx = 0.0
	3419	DO i = 1, klon
	3420	zx = zx + aire(i)
	3421	ENDDO
	3422	qtotal = 0.0
	3423	DO k = 1, klev
	3424	DO i = 1, klon
	3425	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
	3426	. *(paprs(i,k)-paprs(i,k+1))/RG
	3427	ENDDO
	3428	ENDDO
	3429	c
	3430	qcheck = qtotal/zx
	3431	c
	3432	RETURN
	3433	END
	3434	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	3435	IMPLICIT none
	3436	c
	3437	c Tranformer une variable de la grille physique a
	3438	c la grille d'ecriture
	3439	c
	3440	INTEGER nfield,nlon,iim,jjmp1, jjm
	3441	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
	3442	c
	3443	INTEGER i, n, ig
	3444	c
	3445	jjm = jjmp1 - 1
	3446	DO n = 1, nfield
	3447	DO i=1,iim
	3448	ecrit(i,n) = fi(1,n)
	3449	ecrit(i+jjm*iim,n) = fi(nlon,n)
	3450	ENDDO
	3451	DO ig = 1, nlon - 2
	3452	ecrit(iim+ig,n) = fi(1+ig,n)
	3453	ENDDO
	3454	ENDDO
	3455	RETURN
	3456	END

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

Download in other formats:

Original Format