Context Navigation

physiq.F @ 952

Last change on this file since 952 was 952, checked in by lmdzadmin, 17 years ago
On deplace variables SAVE de physiq dans phys_state_var_mod IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 109.9 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
[864]	1076	#ifdef INCA_AER
	1077	#ifdef CPP_COUPLE
	1078	! Aerosol optical properties by INCA model
	1079	CHARACTER*4 :: rfname(9)
	1080
	1081	#endif
	1082	#endif
[524]	1083	REAL aerindex(klon) ! POLDER aerosol index
	1084
	1085	! Parameters
	1086	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
	1087	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
[559]	1088	cym
	1089	SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
[766]	1090	c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
[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.
[864]	1211	#ifdef INCA_AER
	1212	#ifdef CPP_COUPLE
	1213	tau_inca(:,:,:,:) = 0.
	1214	piz_inca(:,:,:,:) = 0.
	1215	cg_inca(:,:,:,:) = 0.
	1216	ccm(:,:,:) = 0.
	1217	topswai_inca(:) = 0.
	1218	topswad_inca(:) = 0.
	1219	topswad0_inca(:) = 0.
	1220	topsw_inca(:,:) = 0.
	1221	topsw0_inca(:,:) = 0.
	1222	solswai_inca(:) = 0.
	1223	solswad_inca(:) = 0.
	1224	solswad0_inca(:) = 0.
	1225	solsw_inca(:,:) = 0.
	1226	solsw0_inca(:,:) = 0.
	1227	#endif
	1228	#endif
[645]	1229	!rv
[864]	1230	!ACo
[645]	1231	rnebcon0(:,:) = 0.0
	1232	clwcon0(:,:) = 0.0
	1233	rnebcon(:,:) = 0.0
	1234	clwcon(:,:) = 0.0
[888]	1235	! fin anne
[645]	1236
[687]	1237	cIM
	1238	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
[524]	1239	c
	1240	c appel a la lecture du run.def physique
	1241	c
	1242	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
[883]	1243	. ok_instan, ok_hf,
[889]	1244	. solarlong0,seuil_inversion,
[879]	1245	. fact_cldcon, facttemps,ok_newmicro,iflag_radia,
[878]	1246	. iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,
[524]	1247	. ok_ade, ok_aie,
[541]	1248	. bl95_b0, bl95_b1,
[879]	1249	. iflag_thermals,nsplit_thermals,
	1250	cnv flags pour la convection et les poches froides
	1251	. iflag_coupl,iflag_clos,iflag_wake)
[524]	1252
[879]	1253	print*,'iflag_coupl,iflag_clos,iflag_wake',
	1254	. iflag_coupl,iflag_clos,iflag_wake
	1255
[524]	1256	c
	1257	c
	1258	c Initialiser les compteurs:
	1259	c
	1260	itap = 0
	1261	itaprad = 0
[782]	1262
[878]	1263	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1264	!! Un petit travail à faire ici.
	1265	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1266
	1267	if (iflag_pbl>1) then
	1268	PRINT*, "Using method MELLOR&YAMADA"
	1269	endif
	1270
[524]	1271	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
[782]	1272	. rlat,rlon,pctsrf, ftsol,
	1273	. ocean, ok_veget,
[888]	1274	. falb1, falb2, rain_fall,snow_fall,
[782]	1275	. solsw, sollw,
	1276	. radsol,clesphy0,
[524]	1277	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
[878]	1278	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon,
[937]	1279	. pbl_tke, zmax0, f0, ema_work1, ema_work2)
[524]	1280
[878]	1281
	1282
	1283
	1284	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1285
	1286
[766]	1287	DO i=1,klon
	1288	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
	1289	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
	1290	$ THEN
[782]	1291	WRITE(,)
	1292	$ 'physiq apres lecture de restart: pb sous surface au point ',
	1293	$ i, pctsrf(i, 1 : nbsrf)
[766]	1294	ENDIF
	1295	ENDDO
[782]	1296
[524]	1297	radpas = NINT( 86400./dtime/nbapp_rad)
	1298	c
	1299	C on remet le calendrier a zero
	1300	c
	1301	IF (raz_date .eq. 1) THEN
	1302	itau_phy = 0
	1303	ENDIF
	1304
[644]	1305	cIM cf. AM 081204 BEG
	1306	PRINT*,'cycle_diurne3 =',cycle_diurne
	1307	cIM cf. AM 081204 END
[524]	1308	c
[782]	1309	CALL printflag( tabcntr0,radpas,ok_journe,
[524]	1310	, ok_instan, ok_region )
	1311	c
	1312	IF (ABS(dtime-pdtphys).GT.0.001) THEN
	1313	WRITE(lunout,*) 'Pas physique n est pas correct',dtime,
	1314	. pdtphys
	1315	abort_message='Pas physique n est pas correct '
[878]	1316	! call abort_gcm(modname,abort_message,1)
	1317	dtime=pdtphys
[524]	1318	ENDIF
	1319	IF (nlon .NE. klon) THEN
	1320	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon,
	1321	. klon
	1322	abort_message='nlon et klon ne sont pas coherents'
	1323	call abort_gcm(modname,abort_message,1)
	1324	ENDIF
	1325	IF (nlev .NE. klev) THEN
	1326	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev,
	1327	. klev
	1328	abort_message='nlev et klev ne sont pas coherents'
	1329	call abort_gcm(modname,abort_message,1)
	1330	ENDIF
	1331	c
	1332	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
	1333	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	1334	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	1335	abort_message='Nbre d appels au rayonnement insuffisant'
	1336	call abort_gcm(modname,abort_message,1)
	1337	ENDIF
	1338	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
	1339	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=",
	1340	. ok_cvl
	1341	c
	1342	cKE43
	1343	c Initialisation pour la convection de K.E. (sb):
	1344	IF (iflag_con.GE.3) THEN
	1345
	1346	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
[687]	1347	WRITE(lunout,*)
	1348	. "On va utiliser le melange convectif des traceurs qui"
	1349	WRITE(lunout,*)"est calcule dans convect4.3"
	1350	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
[524]	1351
	1352	DO i = 1, klon
	1353	ema_cbmf(i) = 0.
	1354	ema_pcb(i) = 0.
	1355	ema_pct(i) = 0.
	1356	ema_workcbmf(i) = 0.
	1357	ENDDO
	1358	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
	1359	DO i = 1, klon
	1360	ibas_con(i) = 1
[619]	1361	itop_con(i) = 1
[524]	1362	ENDDO
	1363	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
[879]	1364	c===============================================================================
	1365	cCR:04.12.07: initialisations poches froides
	1366	c Controle de ALE et ALP pour la fermeture convective (jyg)
	1367	if (iflag_wake.eq.1) then
	1368	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr
	1369	s ,alp_bl_prescr, ale_bl_prescr)
	1370	c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
	1371	c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
	1372	endif
[524]	1373
[879]	1374	do i = 1,klon
	1375	wake_s(i) = 0.
	1376	wake_fip(i) = 0.
	1377	wake_cstar(i) = 0.
	1378	DO k=1,klev
	1379	wake_deltat(i,k)=0.
	1380	wake_deltaq(i,k)=0.
	1381	ENDDO
	1382	enddo
	1383	c================================================================================
	1384
[524]	1385	ENDIF
	1386
[878]	1387	rugoro=0.
[524]	1388	c34EK
	1389	IF (ok_orodr) THEN
[878]	1390
	1391	rugoro=0.
	1392
	1393	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1394	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
	1395	! justement quand ok_orodr = false.
	1396	! ce rugoro est utilise par la couche limite et fait double emploi
	1397	! avec les paramétrisations spécifiques de Francois Lott.
	1398	! DO i=1,klon
	1399	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1400	! ENDDO
	1401	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1402
[782]	1403	CALL SUGWD(klon,klev,paprs,pplay)
	1404	DO i=1,klon
	1405	zuthe(i)=0.
	1406	zvthe(i)=0.
	1407	if(zstd(i).gt.10.)then
	1408	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	1409	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	1410	endif
	1411	ENDDO
[524]	1412	ENDIF
	1413	c
	1414	c
	1415	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
	1416	WRITE(lunout,*)'La frequence de lecture surface est de ',
	1417	. lmt_pas
	1418	c
[644]	1419	cIM200505 ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
	1420	c IF (ok_mensuel) THEN
	1421	c WRITE(lunout,*)'La frequence de sortie mensuelle est de ',
	1422	c . ecrit_mth
	1423	c ENDIF
	1424	c ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
	1425	c IF (ok_journe) THEN
	1426	c WRITE(lunout,*)'La frequence de sortie journaliere est de ',
	1427	c . ecrit_day
	1428	c ENDIF
	1429	cIM 130904 BEG
	1430	cIM 080205 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
	1431	cIM 170305
	1432	c ecrit_hf = 86400./dtime/12. ! toutes les 2h
	1433	cIM 230305
	1434	cIM200505 ecrit_hf = 86400./dtime *0.25 ! toutes les 6h
	1435	c
	1436	cIM200505 ecrit_hf2mth = ecrit_day/ecrit_hf*30
	1437	c
	1438	cIM200505 IF (ok_journe) THEN
	1439	cIM200505 WRITE(lunout,*)'La frequence de sortie hf est de ',
	1440	cIM200505 . ecrit_hf
	1441	cIM200505 ENDIF
	1442	cIM 130904 END
[524]	1443	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
	1444	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
[644]	1445	c ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
	1446	c ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
	1447	cIM200505 ecrit_ins = NINT(86400./dtime/8.) ! toutes les trois heures
	1448	cIM200505 IF (ok_instan) THEN
	1449	cIM200505 WRITE(lunout,*)'La frequence de sortie instant. est de ',
	1450	cIM200505 . ecrit_ins
	1451	cIM200505 ENDIF
	1452	cIM200505 ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
	1453	cIM200505 IF (ok_region) THEN
	1454	cIM200505 WRITE(lunout,*)'La frequence de sortie region est de ',
	1455	cIM200505 . ecrit_reg
	1456	cIM200505 ENDIF
[687]	1457	cIM 030306 BEG
	1458	cIM ecrit_hf2mth = nombre de pas de temps de calcul de hf par mois apres lequel on ecrit
	1459	cIM : ne pas modifier ecrit_hf2mth
[524]	1460	c
[933]	1461	cIM 250308bad guide ecrit_hf2mth = 30*1/ecrit_hf
	1462	ecrit_hf2mth = ecrit_mth/ecrit_hf
[702]	1463	c ecrit_ins en secondes, chaque pas de temps de la physique
[687]	1464	ecrit_ins = dtime
	1465	cIM on passe les frequences de jours en secondes : ecrit_ins, ecrit_hf, ecrit_day, ecrit_mth, ecrit_tra, ecrit_reg
	1466	ecrit_hf = ecrit_hf * un_jour
[828]	1467	!IM
	1468	IF(ecrit_day.LE.1.) THEN
	1469	ecrit_day = ecrit_day * un_jour !en secondes
	1470	ENDIF
	1471	!IM
[687]	1472	ecrit_mth = ecrit_mth * un_jour
	1473	ecrit_reg = ecrit_reg * un_jour
[702]	1474	ecrit_tra = ecrit_tra * un_jour
[828]	1475	ecrit_ISCCP = ecrit_ISCCP * un_jour
	1476	c
[933]	1477	PRINT*,'physiq ecrit_ hf day mth reg tra ISCCP hf2mth',
	1478	. ecrit_hf,ecrit_day,ecrit_mth,ecrit_reg,ecrit_tra,ecrit_ISCCP,
	1479	. ecrit_hf2mth
[687]	1480	cIM 030306 END
[782]	1481
[524]	1482	capemaxcels = 't_max(X)'
	1483	t2mincels = 't_min(X)'
	1484	t2maxcels = 't_max(X)'
[644]	1485	tinst = 'inst(X)'
	1486	tave = 'ave(X)'
	1487	cIM cf. AM 081204 BEG
	1488	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
	1489	cIM cf. AM 081204 END
[524]	1490	c
	1491	c=============================================================
	1492	c Initialisation des sorties
	1493	c=============================================================
	1494
	1495	#ifdef CPP_IOIPSL
	1496
[909]	1497	c Commente par abderrahmane 11 2 08
	1498	c#ifdef histhf
	1499	c#include "ini_histhf.h"
	1500	c#endif
[524]	1501
[909]	1502	c#ifdef histday
	1503	c#include "ini_histday.h"
[644]	1504	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
	1505	c#include "ini_bilKP_ins.h"
	1506	c#include "ini_bilKP_ave.h"
[909]	1507	c#endif
[524]	1508
[909]	1509	c#ifdef histmth
	1510	c#include "ini_histmth.h"
	1511	c#endif
[524]	1512
[909]	1513	c#ifdef histins
	1514	c#include "ini_histins.h"
	1515	c#endif
[524]	1516
[929]	1517	call phys_output_open(jjmp1,nqmax,nlevSTD,clevSTD,nbteta,
	1518	& ctetaSTD,dtime,presnivs,ok_veget,
	1519	& ocean,iflag_pbl,ok_mensuel,ok_journe,
	1520	& ok_hf,ok_instan,nid_files)
[909]	1521
[524]	1522	#ifdef histISCCP
	1523	#include "ini_histISCCP.h"
	1524	#endif
	1525
	1526	#ifdef histmthNMC
	1527	#include "ini_histmthNMC.h"
	1528	#endif
	1529
[687]	1530	#include "ini_histday_seri.h"
[524]	1531
[687]	1532	#include "ini_paramLMDZ_phy.h"
[524]	1533
[644]	1534	#endif
	1535
[524]	1536	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
	1537	date0 = zjulian
	1538	C date0 = day_ini
	1539	WRITE(,) 'physiq date0 : ',date0
	1540	c
	1541	c
	1542	c
	1543	c Prescrire l'ozone dans l'atmosphere
	1544	c
	1545	c
	1546	cc DO i = 1, klon
	1547	cc DO k = 1, klev
	1548	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
	1549	cc ENDDO
	1550	cc ENDDO
	1551	c
	1552	#ifdef INCA
[766]	1553	call VTe(VTphysiq)
	1554	call VTb(VTinca)
[524]	1555	iii = MOD(NINT(xjour),360)
	1556	calday = FLOAT(iii) + gmtime
	1557	WRITE(lunout,*) 'initial time ', xjour, calday
	1558	#ifdef INCAINFO
	1559	WRITE(lunout,*) 'Appel CHEMINI ...'
	1560	#endif
[766]	1561	CALL chemini(
[524]	1562	$ rg,
	1563	$ ra,
	1564	$ airephy,
	1565	$ rlat,
	1566	$ rlon,
	1567	$ presnivs,
	1568	$ calday,
	1569	$ klon,
	1570	$ nqmax,
	1571	$ pdtphys,
[567]	1572	$ annee_ref,
[524]	1573	$ day_ini)
	1574	#ifdef INCAINFO
	1575	WRITE(lunout,*) 'OK.'
	1576	#endif
[766]	1577	call VTe(VTinca)
	1578	call VTb(VTphysiq)
[524]	1579	#endif
	1580	c
	1581	ENDIF
	1582	c
	1583	c ************** Fin de IF ( debut ) *************
	1584	c
[904]	1585
	1586	! Tendances bidons pour les processus qui n'affectent pas certaines
	1587	! variables.
	1588	du0(:,:)=0.
	1589	dv0(:,:)=0.
	1590	dq0(:,:)=0.
	1591	dql0(:,:)=0.
[524]	1592	c
	1593	c Mettre a zero des variables de sortie (pour securite)
	1594	c
	1595	DO i = 1, klon
	1596	d_ps(i) = 0.0
	1597	ENDDO
	1598	DO k = 1, klev
	1599	DO i = 1, klon
	1600	d_t(i,k) = 0.0
	1601	d_u(i,k) = 0.0
	1602	d_v(i,k) = 0.0
	1603	ENDDO
	1604	ENDDO
	1605	DO iq = 1, nqmax
	1606	DO k = 1, klev
	1607	DO i = 1, klon
	1608	d_qx(i,k,iq) = 0.0
	1609	ENDDO
	1610	ENDDO
	1611	ENDDO
[660]	1612	da(:,:)=0.
	1613	mp(:,:)=0.
	1614	phi(:,:,:)=0.
[524]	1615	c
	1616	c Ne pas affecter les valeurs entrees de u, v, h, et q
	1617	c
	1618	DO k = 1, klev
	1619	DO i = 1, klon
	1620	t_seri(i,k) = t(i,k)
	1621	u_seri(i,k) = u(i,k)
	1622	v_seri(i,k) = v(i,k)
	1623	q_seri(i,k) = qx(i,k,ivap)
	1624	ql_seri(i,k) = qx(i,k,iliq)
	1625	qs_seri(i,k) = 0.
	1626	ENDDO
	1627	ENDDO
	1628	IF (nqmax.GE.3) THEN
	1629	DO iq = 3, nqmax
	1630	DO k = 1, klev
	1631	DO i = 1, klon
	1632	tr_seri(i,k,iq-2) = qx(i,k,iq)
	1633	ENDDO
	1634	ENDDO
	1635	ENDDO
	1636	ELSE
	1637	DO k = 1, klev
	1638	DO i = 1, klon
	1639	tr_seri(i,k,1) = 0.0
	1640	ENDDO
	1641	ENDDO
	1642	ENDIF
	1643	C
	1644	DO i = 1, klon
	1645	ztsol(i) = 0.
	1646	ENDDO
	1647	DO nsrf = 1, nbsrf
	1648	DO i = 1, klon
	1649	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
	1650	ENDDO
	1651	ENDDO
[687]	1652	cIM
	1653	IF (ip_ebil_phy.ge.1) THEN
[524]	1654	ztit='after dynamic'
[687]	1655	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
[524]	1656	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1657	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1658	C Comme les tendances de la physique sont ajoute dans la dynamique,
	1659	C on devrait avoir que la variation d'entalpie par la dynamique
	1660	C est egale a la variation de la physique au pas de temps precedent.
	1661	C Donc la somme de ces 2 variations devrait etre nulle.
[687]	1662	call diagphy(airephy,ztit,ip_ebil_phy
[524]	1663	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1664	e , zero_v, zero_v, zero_v, ztsol
	1665	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
	1666	s , fs_bound, fq_bound )
	1667	END IF
	1668
	1669	c Diagnostiquer la tendance dynamique
	1670	c
	1671	IF (ancien_ok) THEN
	1672	DO k = 1, klev
	1673	DO i = 1, klon
	1674	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	1675	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
	1676	ENDDO
	1677	ENDDO
	1678	ELSE
	1679	DO k = 1, klev
	1680	DO i = 1, klon
	1681	d_t_dyn(i,k) = 0.0
	1682	d_q_dyn(i,k) = 0.0
	1683	ENDDO
	1684	ENDDO
	1685	ancien_ok = .TRUE.
	1686	ENDIF
	1687	c
	1688	c Ajouter le geopotentiel du sol:
	1689	c
	1690	DO k = 1, klev
	1691	DO i = 1, klon
	1692	zphi(i,k) = pphi(i,k) + pphis(i)
	1693	ENDDO
	1694	ENDDO
	1695	c
	1696	c Verifier les temperatures
	1697	c
[687]	1698	cIM BEG
	1699	IF (check) THEN
	1700	amn=MIN(ftsol(1,is_ter),1000.)
	1701	amx=MAX(ftsol(1,is_ter),-1000.)
	1702	DO i=2, klon
	1703	amn=MIN(ftsol(i,is_ter),amn)
	1704	amx=MAX(ftsol(i,is_ter),amx)
	1705	ENDDO
	1706	c
	1707	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
	1708	ENDIF !(check) THEN
	1709	cIM END
	1710	c
[524]	1711	CALL hgardfou(t_seri,ftsol,'debutphy')
	1712	c
[687]	1713	cIM BEG
	1714	IF (check) THEN
	1715	amn=MIN(ftsol(1,is_ter),1000.)
	1716	amx=MAX(ftsol(1,is_ter),-1000.)
	1717	DO i=2, klon
	1718	amn=MIN(ftsol(i,is_ter),amn)
	1719	amx=MAX(ftsol(i,is_ter),amx)
	1720	ENDDO
	1721	c
	1722	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
	1723	ENDIF !(check) THEN
	1724	cIM END
	1725	c
[524]	1726	c Incrementer le compteur de la physique
	1727	c
	1728	itap = itap + 1
	1729	julien = MOD(NINT(xjour),360)
	1730	if (julien .eq. 0) julien = 360
	1731	c
	1732	c Mettre en action les conditions aux limites (albedo, sst, etc.).
	1733	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
	1734	c
	1735	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
	1736	WRITE(lunout,*)' PHYS cond julien ',julien
	1737	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
	1738	ENDIF
	1739	c
	1740	c Re-evaporer l'eau liquide nuageuse
	1741	c
	1742	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
	1743	DO i = 1, klon
	1744	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1745	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1746	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1747	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
	1748	zb = MAX(0.0,ql_seri(i,k))
	1749	za = - MAX(0.0,ql_seri(i,k))
	1750	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
	1751	t_seri(i,k) = t_seri(i,k) + za
	1752	q_seri(i,k) = q_seri(i,k) + zb
	1753	ql_seri(i,k) = 0.0
	1754	d_t_eva(i,k) = za
	1755	d_q_eva(i,k) = zb
	1756	ENDDO
	1757	ENDDO
[687]	1758	cIM
	1759	IF (ip_ebil_phy.ge.2) THEN
[524]	1760	ztit='after reevap'
[687]	1761	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime
[524]	1762	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1763	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	1764	call diagphy(airephy,ztit,ip_ebil_phy
[524]	1765	e , zero_v, zero_v, zero_v, zero_v, zero_v
	1766	e , zero_v, zero_v, zero_v, ztsol
	1767	e , d_h_vcol, d_qt, d_ec
	1768	s , fs_bound, fq_bound )
	1769	C
	1770	END IF
[782]	1771
[524]	1772	c
[883]	1773	c=========================================================================
	1774	! Calculs de l'orbite.
	1775	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
	1776	! doit donc etre placé avant radlwsw et pbl_surface
	1777
	1778	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
	1779	! solarlong0
	1780
	1781	if (solarlong0<-999.) then
	1782	CALL orbite(FLOAT(julien),zlongi,dist)
	1783	else
	1784	zlongi=solarlong0 ! longitude solaire vraie
	1785	dist=1. ! distance au soleil / moyenne
	1786	endif
	1787
[942]	1788	if(prt_level.ge.1) print*,'Longitude solaire ',zlongi,solarlong0
[883]	1789
	1790	! Avec ou sans cycle diurne
[524]	1791	IF (cycle_diurne) THEN
	1792	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
	1793	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
	1794	ELSE
	1795	rmu0 = -999.999
	1796	ENDIF
	1797
[766]	1798	if (mydebug) then
	1799	call writefield_phy('u_seri',u_seri,llm)
	1800	call writefield_phy('v_seri',v_seri,llm)
	1801	call writefield_phy('t_seri',t_seri,llm)
	1802	call writefield_phy('q_seri',q_seri,llm)
	1803	endif
[782]	1804
	1805	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	1806	c Appel au pbl_surface : Planetary Boudary Layer et Surface
	1807	c Cela implique tous les interactions des sous-surfaces et la partie diffusion
	1808	c turbulent du couche limit.
	1809	c
	1810	c Certains varibales de sorties de pbl_surface sont utiliser que pour
	1811	c ecriture des fihiers hist_XXXX.nc, ces sont :
	1812	c qsol, zq2m, s_pblh, s_lcl,
	1813	c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
	1814	c s_therm, s_trmb1, s_trmb2, s_trmb3,
	1815	c zxrugs, zu10m, zv10m, fder,
	1816	c zxqsurf, rh2m, zxfluxu, zxfluxv,
	1817	c frugs, agesno, fsollw, fsolsw,
	1818	c d_ts, fevap, fluxlat, t2m,
	1819	c wfbils, wfbilo, fluxt, fluxu, fluxv,
[687]	1820	c
[782]	1821	c Certains ne sont pas utiliser du tout :
	1822	c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
[687]	1823	c
[883]	1824
[782]	1825	CALL pbl_surface(
	1826	e dtime, date0, itap, julien,
	1827	e debut, lafin,
	1828	e rlon, rlat, rugoro, rmu0,
	1829	e rain_fall, snow_fall, solsw, sollw,
	1830	e t_seri, q_seri, u_seri, v_seri,
	1831	e pplay, paprs, pctsrf,
[888]	1832	+ ftsol, falb1, falb2, u10m, v10m,
[782]	1833	s sollwdown, cdragh, cdragm, yu1, yv1,
[888]	1834	s albsol1, albsol2, sens, evap,
[782]	1835	s zxtsol, zxfluxlat, zt2m, qsat2m,
	1836	s d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf,
	1837	s ycoefh, pctsrf_new,
	1838	d qsol, zq2m, s_pblh, s_lcl,
	1839	d s_capCL, s_oliqCL, s_cteiCL,s_pblT,
	1840	d s_therm, s_trmb1, s_trmb2, s_trmb3,
	1841	d zxrugs, zu10m, zv10m, fder,
	1842	d zxqsurf, rh2m, zxfluxu, zxfluxv,
	1843	d frugs, agesno, fsollw, fsolsw,
	1844	d d_ts, fevap, fluxlat, t2m,
	1845	d wfbils, wfbilo, fluxt, fluxu, fluxv,
	1846	- dsens, devap, zxsnow,
[878]	1847	- zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
[524]	1848	c
[782]	1849	c
	1850	pctsrf(:,:) = pctsrf_new(:,:)
	1851
[904]	1852	!-----------------------------------------------------------------------------------------
	1853	! ajout des tendances de la diffusion turbulente
	1854	CALL add_phys_tend(d_u_vdf,d_v_vdf,d_t_vdf,d_q_vdf,dql0,'vdf')
	1855	!-----------------------------------------------------------------------------------------
[766]	1856
	1857	if (mydebug) then
	1858	call writefield_phy('u_seri',u_seri,llm)
	1859	call writefield_phy('v_seri',v_seri,llm)
	1860	call writefield_phy('t_seri',t_seri,llm)
	1861	call writefield_phy('q_seri',q_seri,llm)
	1862	endif
	1863
	1864
[687]	1865	IF (ip_ebil_phy.ge.2) THEN
[782]	1866	ztit='after surface_main'
[687]	1867	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	1868	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	1869	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	1870	call diagphy(airephy,ztit,ip_ebil_phy
[524]	1871	e , zero_v, zero_v, zero_v, zero_v, sens
	1872	e , evap , zero_v, zero_v, ztsol
	1873	e , d_h_vcol, d_qt, d_ec
	1874	s , fs_bound, fq_bound )
	1875	END IF
	1876
[881]	1877	c =================================================================== c
	1878	c Calcul de Qsat
	1879
	1880	DO k = 1, klev
	1881	DO i = 1, klon
	1882	zx_t = t_seri(i,k)
	1883	IF (thermcep) THEN
	1884	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	1885	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	1886	zx_qs = MIN(0.5,zx_qs)
	1887	zcor = 1./(1.-retv*zx_qs)
	1888	zx_qs = zx_qs*zcor
	1889	ELSE
	1890	IF (zx_t.LT.t_coup) THEN
	1891	zx_qs = qsats(zx_t)/pplay(i,k)
	1892	ELSE
	1893	zx_qs = qsatl(zx_t)/pplay(i,k)
	1894	ENDIF
	1895	ENDIF
	1896	zqsat(i,k)=zx_qs
	1897	ENDDO
	1898	ENDDO
	1899
[942]	1900	if (prt_level.ge.1) then
[881]	1901	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
	1902	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
	1903	endif
[524]	1904	c
	1905	c Appeler la convection (au choix)
	1906	c
	1907	DO k = 1, klev
	1908	DO i = 1, klon
	1909	conv_q(i,k) = d_q_dyn(i,k)
	1910	. + d_q_vdf(i,k)/dtime
	1911	conv_t(i,k) = d_t_dyn(i,k)
	1912	. + d_t_vdf(i,k)/dtime
	1913	ENDDO
	1914	ENDDO
	1915	IF (check) THEN
	1916	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	1917	WRITE(lunout,*) "avantcon=", za
	1918	ENDIF
	1919	zx_ajustq = .FALSE.
	1920	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
	1921	IF (zx_ajustq) THEN
	1922	DO i = 1, klon
	1923	z_avant(i) = 0.0
	1924	ENDDO
	1925	DO k = 1, klev
	1926	DO i = 1, klon
	1927	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
	1928	. *(paprs(i,k)-paprs(i,k+1))/RG
	1929	ENDDO
	1930	ENDDO
	1931	ENDIF
	1932	IF (iflag_con.EQ.1) THEN
	1933	stop'reactiver le call conlmd dans physiq.F'
	1934	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
	1935	c . d_t_con, d_q_con,
	1936	c . rain_con, snow_con, ibas_con, itop_con)
	1937	ELSE IF (iflag_con.EQ.2) THEN
	1938	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
[782]	1939	e conv_t, conv_q, -evap, omega,
[524]	1940	s d_t_con, d_q_con, rain_con, snow_con,
	1941	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	1942	s kcbot, kctop, kdtop, pmflxr, pmflxs)
	1943	WHERE (rain_con < 0.) rain_con = 0.
	1944	WHERE (snow_con < 0.) snow_con = 0.
	1945	DO i = 1, klon
	1946	ibas_con(i) = klev+1 - kcbot(i)
	1947	itop_con(i) = klev+1 - kctop(i)
	1948	ENDDO
	1949	ELSE IF (iflag_con.GE.3) THEN
	1950	c nb of tracers for the KE convection:
[619]	1951	c MAF la partie traceurs est faite dans phytrac
	1952	c on met ntra=1 pour limiter les appels mais on peut
	1953	c supprimer les calculs / ftra.
	1954	ntra = 1
[879]	1955
	1956	c=====================================================================================
	1957	cajout pour la parametrisation des poches froides:
	1958	ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
	1959	do k=1,klev
	1960	do i=1,klon
	1961	if (iflag_wake.eq.1) then
	1962	t_wake(i,k) = t_seri(i,k)
	1963	. +(1-wake_s(i))*wake_deltat(i,k)
	1964	q_wake(i,k) = q_seri(i,k)
	1965	. +(1-wake_s(i))*wake_deltaq(i,k)
	1966	t_undi(i,k) = t_seri(i,k)
	1967	. -wake_s(i)*wake_deltat(i,k)
	1968	q_undi(i,k) = q_seri(i,k)
	1969	. -wake_s(i)*wake_deltaq(i,k)
	1970	else
	1971	t_wake(i,k) = t_seri(i,k)
	1972	q_wake(i,k) = q_seri(i,k)
	1973	t_undi(i,k) = t_seri(i,k)
	1974	q_undi(i,k) = q_seri(i,k)
	1975	endif
	1976	enddo
	1977	enddo
	1978
	1979	cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
	1980	cc-- pour le soulevement des particules dans le modele convectif
	1981	c
	1982	do i = 1,klon
	1983	ALE(i) = 0.
	1984	ALP(i) = 0.
	1985	enddo
	1986	c
	1987	ccalcul de ale_wake et alp_wake
	1988	do i = 1,klon
	1989	if (iflag_wake.eq.1) then
	1990	ale_wake(i) = 0.5wake_cstar(i)*2
	1991	alp_wake(i) = wake_fip(i)
	1992	else
	1993	ale_wake(i) = 0.
	1994	alp_wake(i) = 0.
	1995	endif
	1996	enddo
	1997	ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
	1998	cdans le thermique sinon
	1999	if (iflag_coupl.eq.0) then
[942]	2000	if (debut) print*,'ALE et ALP imposes'
[879]	2001	do i = 1,klon
	2002	con ne couple que ale
	2003	c ALE(i) = max(ale_wake(i),Ale_bl(i))
	2004	ALE(i) = max(ale_wake(i),ale_bl_prescr)
	2005	con ne couple que alp
	2006	c ALP(i) = alp_wake(i) + Alp_bl(i)
	2007	ALP(i) = alp_wake(i) + alp_bl_prescr
	2008	enddo
	2009	else
	2010	print*,'ALE et ALP couples au thermique'
	2011	do i = 1,klon
	2012	ALE(i) = max(ale_wake(i),Ale_bl(i))
	2013	ALP(i) = alp_wake(i) + Alp_bl(i)
	2014	c write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
	2015	c write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
	2016	enddo
	2017	endif
	2018
	2019	cfin calcul ale et alp
	2020	c=================================================================================================
	2021
	2022
[524]	2023	c sb, oct02:
	2024	c Schema de convection modularise et vectorise:
	2025	c (driver commun aux versions 3 et 4)
	2026	c
	2027	IF (ok_cvl) THEN ! new driver for convectL
	2028
[879]	2029	CALL concvl (iflag_con,iflag_clos,
	2030	. dtime,paprs,pplay,t_undi,q_undi,
	2031	. t_wake,q_wake,
	2032	. u_seri,v_seri,tr_seri,nbtr,
	2033	. ALE,ALP,
[524]	2034	. ema_work1,ema_work2,
	2035	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
[879]	2036	. rain_con, snow_con, ibas_con, itop_con, sigd,
[524]	2037	. upwd,dnwd,dnwd0,
[879]	2038	. Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl,
[619]	2039	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd,
[879]	2040	. pmflxr,pmflxs,da,phi,mp,
	2041	. ftd,fqd,lalim_conv,wght_th)
[619]	2042
[524]	2043	cIM cf. FH
	2044	clwcon0=qcondc
[619]	2045	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
[524]	2046
	2047	ELSE ! ok_cvl
[619]	2048	c MAF conema3 ne contient pas les traceurs
[524]	2049	CALL conema3 (dtime,
	2050	. paprs,pplay,t_seri,q_seri,
[619]	2051	. u_seri,v_seri,tr_seri,ntra,
[524]	2052	. ema_work1,ema_work2,
	2053	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
	2054	. rain_con, snow_con, ibas_con, itop_con,
	2055	. upwd,dnwd,dnwd0,bas,top,
	2056	. Ma,cape,tvp,rflag,
	2057	. pbase
	2058	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
	2059	. ,clwcon0)
	2060
	2061	ENDIF ! ok_cvl
	2062
[766]	2063	c
	2064	c Correction precip
	2065	rain_con = rain_con * cvl_corr
	2066	snow_con = snow_con * cvl_corr
	2067	c
	2068
[524]	2069	IF (.NOT. ok_gust) THEN
	2070	do i = 1, klon
	2071	wd(i)=0.0
	2072	enddo
	2073	ENDIF
	2074
	2075	c =================================================================== c
	2076	c Calcul des proprietes des nuages convectifs
	2077	c
	2078
	2079	c calcul des proprietes des nuages convectifs
	2080	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
	2081	call clouds_gno
	2082	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
	2083
	2084	c =================================================================== c
	2085
	2086	DO i = 1, klon
	2087	ema_pcb(i) = pbase(i)
	2088	ENDDO
	2089	DO i = 1, klon
[879]	2090
	2091	! L'idicage de itop_con peut cacher un pb potentiel
	2092	! FH sous la dictee de JYG, CR
	2093	ema_pct(i) = paprs(i,itop_con(i)+1)
	2094
[878]	2095	if (itop_con(i).gt.klev-3) then
	2096	print*,'La convection monte trop haut '
	2097	print*,'itop_con(,',i,',)=',itop_con(i)
	2098	endif
[524]	2099	ENDDO
	2100	DO i = 1, klon
	2101	ema_cbmf(i) = ema_workcbmf(i)
	2102	ENDDO
[881]	2103	ELSE IF (iflag_con.eq.0) THEN
	2104	write(lunout,*) 'On n appelle pas la convection'
	2105	clwcon0=0.
	2106	rnebcon0=0.
	2107	d_t_con=0.
	2108	d_q_con=0.
	2109	d_u_con=0.
	2110	d_v_con=0.
	2111	rain_con=0.
	2112	snow_con=0.
	2113	bas=1
	2114	top=1
[524]	2115	ELSE
	2116	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
	2117	CALL abort
	2118	ENDIF
	2119
	2120	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
	2121	c . d_u_con, d_v_con)
	2122
[904]	2123	!-----------------------------------------------------------------------------------------
	2124	! ajout des tendances de la diffusion turbulente
	2125	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
	2126	!-----------------------------------------------------------------------------------------
[766]	2127
	2128	if (mydebug) then
	2129	call writefield_phy('u_seri',u_seri,llm)
	2130	call writefield_phy('v_seri',v_seri,llm)
	2131	call writefield_phy('t_seri',t_seri,llm)
	2132	call writefield_phy('q_seri',q_seri,llm)
	2133	endif
	2134
[687]	2135	cIM
	2136	IF (ip_ebil_phy.ge.2) THEN
[524]	2137	ztit='after convect'
[687]	2138	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2139	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2140	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	2141	call diagphy(airephy,ztit,ip_ebil_phy
[524]	2142	e , zero_v, zero_v, zero_v, zero_v, zero_v
	2143	e , zero_v, rain_con, snow_con, ztsol
	2144	e , d_h_vcol, d_qt, d_ec
	2145	s , fs_bound, fq_bound )
	2146	END IF
	2147	C
	2148	IF (check) THEN
	2149	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2150	WRITE(lunout,*)"aprescon=", za
	2151	zx_t = 0.0
	2152	za = 0.0
	2153	DO i = 1, klon
	2154	za = za + airephy(i)/FLOAT(klon)
	2155	zx_t = zx_t + (rain_con(i)+
	2156	. snow_con(i))*airephy(i)/FLOAT(klon)
	2157	ENDDO
	2158	zx_t = zx_t/za*dtime
	2159	WRITE(lunout,*)"Precip=", zx_t
	2160	ENDIF
	2161	IF (zx_ajustq) THEN
	2162	DO i = 1, klon
	2163	z_apres(i) = 0.0
	2164	ENDDO
	2165	DO k = 1, klev
	2166	DO i = 1, klon
	2167	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
	2168	. *(paprs(i,k)-paprs(i,k+1))/RG
	2169	ENDDO
	2170	ENDDO
	2171	DO i = 1, klon
	2172	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
	2173	. /z_apres(i)
	2174	ENDDO
	2175	DO k = 1, klev
	2176	DO i = 1, klon
	2177	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
	2178	. z_factor(i).LT.(1.0-1.0E-08)) THEN
	2179	q_seri(i,k) = q_seri(i,k) * z_factor(i)
	2180	ENDIF
	2181	ENDDO
	2182	ENDDO
	2183	ENDIF
	2184	zx_ajustq=.FALSE.
[879]	2185
[524]	2186	c
[879]	2187	c=============================================================================
	2188	cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
	2189	cpour la couche limite diffuse pour l instant
	2190	c
	2191	if (iflag_wake.eq.1) then
	2192	DO k=1,klev
	2193	DO i=1,klon
	2194	dt_dwn(i,k) = ftd(i,k)
	2195	wdt_PBL(i,k) = 0.
	2196	dq_dwn(i,k) = fqd(i,k)
	2197	wdq_PBL(i,k) = 0.
	2198	M_dwn(i,k) = dnwd0(i,k)
	2199	M_up(i,k) = upwd(i,k)
	2200	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
	2201	udt_PBL(i,k) = 0.
	2202	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
	2203	udq_PBL(i,k) = 0.
	2204	ENDDO
	2205	ENDDO
	2206	c
	2207	ccalcul caracteristiques de la poche froide
	2208	call calWAKE (paprs,pplay,dtime
	2209	: ,t_seri,q_seri,omega,ibas_con
	2210	: ,dt_dwn,dq_dwn,M_dwn,M_up
	2211	: ,dt_a,dq_a,sigd
	2212	: ,wdt_PBL,wdq_PBL
	2213	: ,udt_PBL,udq_PBL
	2214	o ,wake_deltat,wake_deltaq,wake_dth
	2215	o ,wake_h,wake_s,wake_dens
	2216	o ,wake_pe,wake_fip,wake_gfl
	2217	o ,dt_wake,dq_wake
	2218	o ,wake_k, t_undi,q_undi
	2219	o ,wake_omgbdth,wake_dp_omgb
	2220	o ,wake_dtKE,wake_dqKE
	2221	o ,wake_dtPBL,wake_dqPBL
	2222	o ,wake_omg,wake_dp_deltomg
	2223	o ,wake_spread,wake_Cstar,wake_d_deltat_gw
	2224	o ,wake_ddeltat,wake_ddeltaq)
	2225	c
[904]	2226	!-----------------------------------------------------------------------------------------
	2227	! ajout des tendances des poches froides
	2228	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
	2229	! coherent avec les autres d_t_...
	2230	d_t_wake(:,:)=dt_wake(:,:)*dtime
	2231	d_q_wake(:,:)=dq_wake(:,:)*dtime
	2232	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
	2233	!-----------------------------------------------------------------------------------------
[879]	2234
	2235	endif
	2236	c print*,'apres callwake iflag_cldcon=', iflag_cldcon
	2237	c
[541]	2238	c===================================================================
	2239	c Convection seche (thermiques ou ajustement)
	2240	c===================================================================
[524]	2241	c
[878]	2242	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri
	2243	s ,seuil_inversion,weak_inversion,dthmin)
	2244
	2245
	2246
	2247	d_t_ajsb(:,:)=0.
	2248	d_q_ajsb(:,:)=0.
[541]	2249	d_t_ajs(:,:)=0.
	2250	d_u_ajs(:,:)=0.
	2251	d_v_ajs(:,:)=0.
	2252	d_q_ajs(:,:)=0.
[878]	2253	clwcon0th(:,:)=0.
[541]	2254	c
[557]	2255	IF(prt_level>9)WRITE(lunout,*)
	2256	. 'AVANT LA CONVECTION SECHE , iflag_thermals='
[541]	2257	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
	2258	if(iflag_thermals.lt.0) then
	2259	c Rien
	2260	c ====
[557]	2261	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
[541]	2262
[878]	2263
[541]	2264	else
[878]	2265
[541]	2266	c Thermiques
	2267	c ==========
[557]	2268	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals='
[541]	2269	s ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
[878]	2270
	2271
	2272	if (iflag_thermals.gt.1) then
[541]	2273	call calltherm(pdtphys
[878]	2274	s ,pplay,paprs,pphi,weak_inversion
	2275	s ,u_seri,v_seri,t_seri,q_seri,zqsat,debut
[541]	2276	s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
[879]	2277	s ,fm_therm,entr_therm,zqasc,clwcon0th,lmax_th,ratqscth
	2278	s ,ratqsdiff,zqsatth
	2279	con rajoute ale et alp, et les caracteristiques de la couche alim
[927]	2280	s ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0)
[878]	2281	endif
	2282
	2283	! call calltherm(pdtphys
	2284	! s ,pplay,paprs,pphi
	2285	! s ,u_seri,v_seri,t_seri,q_seri
	2286	! s ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs
	2287	! s ,fm_therm,entr_therm)
	2288
	2289	c Ajustement sec
	2290	c ==============
	2291
	2292	! Dans le cas où on active les thermiques, on fait partir l'ajustement
	2293	! a partir du sommet des thermiques.
	2294	! Dans le cas contraire, on demarre au niveau 1.
	2295
	2296	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
	2297
	2298	if(iflag_thermals.eq.0) then
	2299	IF(prt_level>9)WRITE(lunout,*)'ajsec'
	2300	limbas(:)=1
	2301	else
	2302	limbas(:)=lmax_th(:)
	2303	endif
	2304
	2305	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
	2306	! pour des test de convergence numerique.
	2307	! Le nouveau ajsec est a priori mieux, meme pour le cas
	2308	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
	2309	! non nulles numeriquement pour des mailles non concernees.
	2310
	2311	if (iflag_thermals.eq.0) then
	2312	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri
	2313	s , d_t_ajsb, d_q_ajsb)
	2314	else
	2315	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas
	2316	s , d_t_ajsb, d_q_ajsb)
	2317	endif
	2318
[904]	2319	!-----------------------------------------------------------------------------------------
	2320	! ajout des tendances de l'ajustement sec ou des thermiques
	2321	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
[878]	2322	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
	2323	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
	2324
[904]	2325	!-----------------------------------------------------------------------------------------
	2326
[878]	2327	endif
	2328
[541]	2329	endif
	2330	c
	2331	c===================================================================
[687]	2332	cIM
	2333	IF (ip_ebil_phy.ge.2) THEN
[524]	2334	ztit='after dry_adjust'
[687]	2335	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2336	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2337	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2338	END IF
	2339
	2340
	2341	c-------------------------------------------------------------------------
	2342	c Caclul des ratqs
	2343	c-------------------------------------------------------------------------
	2344
	2345	c print*,'calcul des ratqs'
	2346	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
	2347	c ----------------
	2348	c on ecrase le tableau ratqsc calcule par clouds_gno
	2349	if (iflag_cldcon.eq.1) then
	2350	do k=1,klev
	2351	do i=1,klon
	2352	if(ptconv(i,k)) then
	2353	ratqsc(i,k)=ratqsbas
	2354	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
	2355	else
	2356	ratqsc(i,k)=0.
	2357	endif
	2358	enddo
	2359	enddo
[878]	2360	else if (iflag_cldcon.eq.4) then
	2361	ptconvth(:,:)=.false.
	2362	ratqsc(:,:)=0.
[942]	2363	if(prt_level.ge.9) print*,'avant clouds_gno thermique'
[878]	2364	call clouds_gno
	2365	s (klon,klev,q_seri,zqsat,clwcon0th,ptconvth,ratqsc,rnebcon0th)
[942]	2366	if(prt_level.ge.9) print*,' CLOUDS_GNO OK'
[524]	2367	endif
	2368
	2369	c ratqs stables
	2370	c -------------
	2371
[878]	2372	if (iflag_ratqs.eq.0) then
[524]	2373
[878]	2374	! Le cas iflag_ratqs=0 correspond a la version IPCC 2005 du modele.
	2375	do k=1,klev
	2376	do i=1, klon
	2377	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
	2378	s min((paprs(i,1)-pplay(i,k))/(paprs(i,1)-30000.),1.)
	2379	enddo
	2380	enddo
	2381
	2382	! Pour iflag_ratqs=1 ou 2, le ratqs est constant au dessus de
	2383	! 300 hPa (ratqshaut), varie lineariement en fonction de la pression
	2384	! entre 600 et 300 hPa et est soit constant (ratqsbas) pour iflag_ratqs=1
	2385	! soit lineaire (entre 0 a la surface et ratqsbas) pour iflag_ratqs=2
	2386	! Il s'agit de differents tests dans la phase de reglage du modele
	2387	! avec thermiques.
	2388
	2389	else if (iflag_ratqs.eq.1) then
	2390
	2391	do k=1,klev
	2392	do i=1, klon
	2393	if (pplay(i,k).ge.60000.) then
	2394	ratqss(i,k)=ratqsbas
	2395	else if ((pplay(i,k).ge.30000.).and.
	2396	s (pplay(i,k).lt.60000.)) then
	2397	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
	2398	s (60000.-pplay(i,k))/(60000.-30000.)
	2399	else
	2400	ratqss(i,k)=ratqshaut
	2401	endif
	2402	enddo
	2403	enddo
	2404
	2405	else
	2406
	2407	do k=1,klev
	2408	do i=1, klon
	2409	if (pplay(i,k).ge.60000.) then
	2410	ratqss(i,k)=ratqsbas
	2411	s *(paprs(i,1)-pplay(i,k))/(paprs(i,1)-60000.)
	2412	else if ((pplay(i,k).ge.30000.).and.
	2413	s (pplay(i,k).lt.60000.)) then
	2414	ratqss(i,k)=ratqsbas+(ratqshaut-ratqsbas)*
	2415	s (60000.-pplay(i,k))/(60000.-30000.)
	2416	else
	2417	ratqss(i,k)=ratqshaut
	2418	endif
	2419	enddo
	2420	enddo
	2421	endif
	2422
	2423
	2424
	2425
[524]	2426	c ratqs final
	2427	c -----------
[878]	2428
	2429	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2
	2430	s .or.iflag_cldcon.eq.4) then
	2431
	2432	! On ajoute une constante au ratqsc*2 pour tenir compte de
	2433	! fluctuations turbulentes de petite echelle
	2434
	2435	do k=1,klev
	2436	do i=1,klon
	2437	if ((fm_therm(i,k).gt.1.e-10)) then
	2438	ratqsc(i,k)=sqrt(ratqsc(i,k)2+0.052)
	2439	endif
	2440	enddo
	2441	enddo
	2442
	2443	! les ratqs sont une conbinaison de ratqss et ratqsc
	2444	! 1800s, en dur pour le moment, est le temps de
	2445	! relaxation des ratqs
	2446
	2447	facteur=exp(-pdtphys/1800.)
	2448
	2449	print*,'WARNING ratqs a revoir '
	2450	ratqs(:,:)=ratqsc(:,:)(1.-facteur)+ratqs(:,:)facteur
	2451	ratqs(:,:)=sqrt(ratqs(:,:)2+ratqss(:,:)2)
	2452
[524]	2453	else
[878]	2454	! on ne prend que le ratqs stable pour fisrtilp
[524]	2455	ratqs(:,:)=ratqss(:,:)
	2456	endif
	2457
	2458
	2459	c
	2460	c Appeler le processus de condensation a grande echelle
	2461	c et le processus de precipitation
	2462	c-------------------------------------------------------------------------
	2463	CALL fisrtilp(dtime,paprs,pplay,
	2464	. t_seri, q_seri,ptconv,ratqs,
	2465	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
	2466	. rain_lsc, snow_lsc,
	2467	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
	2468	. frac_impa, frac_nucl,
	2469	. prfl, psfl, rhcl)
	2470
	2471	WHERE (rain_lsc < 0) rain_lsc = 0.
	2472	WHERE (snow_lsc < 0) snow_lsc = 0.
[904]	2473	!-----------------------------------------------------------------------------------------
	2474	! ajout des tendances de la diffusion turbulente
	2475	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
	2476	!-----------------------------------------------------------------------------------------
[524]	2477	DO k = 1, klev
	2478	DO i = 1, klon
	2479	cldfra(i,k) = rneb(i,k)
	2480	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
	2481	ENDDO
	2482	ENDDO
	2483	IF (check) THEN
	2484	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2485	WRITE(lunout,*)"apresilp=", za
	2486	zx_t = 0.0
	2487	za = 0.0
	2488	DO i = 1, klon
	2489	za = za + airephy(i)/FLOAT(klon)
	2490	zx_t = zx_t + (rain_lsc(i)
	2491	. + snow_lsc(i))*airephy(i)/FLOAT(klon)
	2492	ENDDO
	2493	zx_t = zx_t/za*dtime
	2494	WRITE(lunout,*)"Precip=", zx_t
	2495	ENDIF
[687]	2496	cIM
	2497	IF (ip_ebil_phy.ge.2) THEN
[524]	2498	ztit='after fisrt'
[687]	2499	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2500	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2501	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	2502	call diagphy(airephy,ztit,ip_ebil_phy
[524]	2503	e , zero_v, zero_v, zero_v, zero_v, zero_v
	2504	e , zero_v, rain_lsc, snow_lsc, ztsol
	2505	e , d_h_vcol, d_qt, d_ec
	2506	s , fs_bound, fq_bound )
	2507	END IF
[766]	2508
	2509	if (mydebug) then
	2510	call writefield_phy('u_seri',u_seri,llm)
	2511	call writefield_phy('v_seri',v_seri,llm)
	2512	call writefield_phy('t_seri',t_seri,llm)
	2513	call writefield_phy('q_seri',q_seri,llm)
	2514	endif
	2515
[524]	2516	c
	2517	c-------------------------------------------------------------------
	2518	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
	2519	c-------------------------------------------------------------------
	2520
	2521	c 1. NUAGES CONVECTIFS
	2522	c
[644]	2523	cIM cf FH
	2524	c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
[878]	2525	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
[644]	2526	snow_tiedtke=0.
	2527	c print*,'avant calcul de la pseudo precip '
	2528	c print*,'iflag_cldcon',iflag_cldcon
	2529	if (iflag_cldcon.eq.-1) then
	2530	rain_tiedtke=rain_con
	2531	else
	2532	c print*,'calcul de la pseudo precip '
	2533	rain_tiedtke=0.
	2534	c print*,'calcul de la pseudo precip 0'
	2535	do k=1,klev
	2536	do i=1,klon
	2537	if (d_q_con(i,k).lt.0.) then
	2538	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys
	2539	s *(paprs(i,k)-paprs(i,k+1))/rg
	2540	endif
	2541	enddo
	2542	enddo
	2543	endif
	2544	c
	2545	c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
	2546	c
[524]	2547
	2548	c Nuages diagnostiques pour Tiedtke
	2549	CALL diagcld1(paprs,pplay,
[644]	2550	cIM cf FH . rain_con,snow_con,ibas_con,itop_con,
	2551	. rain_tiedtke,snow_tiedtke,ibas_con,itop_con,
[524]	2552	. diafra,dialiq)
	2553	DO k = 1, klev
	2554	DO i = 1, klon
	2555	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	2556	cldliq(i,k) = dialiq(i,k)
	2557	cldfra(i,k) = diafra(i,k)
	2558	ENDIF
	2559	ENDDO
	2560	ENDDO
	2561
[878]	2562	ELSE IF (iflag_cldcon.ge.3) THEN
[524]	2563	c On prend pour les nuages convectifs le max du calcul de la
[766]	2564	c convection et du calcul du pas de temps precedent diminue d'un facteur
[524]	2565	c facttemps
	2566	facteur = pdtphys *facttemps
	2567	do k=1,klev
	2568	do i=1,klon
	2569	rnebcon(i,k)=rnebcon(i,k)*facteur
	2570	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
	2571	s then
	2572	rnebcon(i,k)=rnebcon0(i,k)
	2573	clwcon(i,k)=clwcon0(i,k)
	2574	endif
	2575	enddo
	2576	enddo
	2577
[644]	2578	c
[766]	2579	cjq - introduce the aerosol direct and first indirect radiative forcings
	2580	cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
	2581	IF (ok_ade.OR.ok_aie) THEN
[864]	2582	#if defined(CPP_COUPLE) && !defined(INCA_AER)
[766]	2583	! Get sulfate aerosol distribution
	2584	CALL readsulfate(rjourvrai, debut, sulfate)
	2585	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
	2586
	2587	! Calculate aerosol optical properties (Olivier Boucher)
	2588	CALL aeropt(pplay, paprs, t_seri, sulfate, rhcl,
	2589	. tau_ae, piz_ae, cg_ae, aerindex)
[864]	2590	#endif
	2591	#if !defined(CPP_COUPLE)
	2592	! Get sulfate aerosol distribution
	2593	CALL readsulfate(rjourvrai, debut, sulfate)
	2594	CALL readsulfate_preind(rjourvrai, debut, sulfate_pi)
	2595
	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	cym
	2601	ELSE
	2602	tau_ae(:,:,:)=0.0
	2603	piz_ae(:,:,:)=0.0
	2604	cg_ae(:,:,:)=0.0
	2605	cym
	2606	ENDIF
	2607
	2608	c
[524]	2609	cIM calcul nuages par le simulateur ISCCP
[644]	2610	c
[839]	2611	#ifdef histISCCP
[524]	2612	IF (ok_isccp) THEN
[828]	2613	cIM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
	2614	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
[644]	2615	#include "calcul_simulISCCP.h"
[828]	2616	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
[524]	2617	ENDIF !ok_isccp
[839]	2618	#endif
[524]	2619
	2620	c On prend la somme des fractions nuageuses et des contenus en eau
	2621	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	2622	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
	2623
	2624	ENDIF
	2625
	2626	c
	2627	c 2. NUAGES STARTIFORMES
	2628	c
	2629	IF (ok_stratus) THEN
	2630	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
	2631	DO k = 1, klev
	2632	DO i = 1, klon
	2633	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	2634	cldliq(i,k) = dialiq(i,k)
	2635	cldfra(i,k) = diafra(i,k)
	2636	ENDIF
	2637	ENDDO
	2638	ENDDO
	2639	ENDIF
	2640	c
	2641	c Precipitation totale
	2642	c
	2643	DO i = 1, klon
	2644	rain_fall(i) = rain_con(i) + rain_lsc(i)
	2645	snow_fall(i) = snow_con(i) + snow_lsc(i)
	2646	ENDDO
[687]	2647	cIM
	2648	IF (ip_ebil_phy.ge.2) THEN
[524]	2649	ztit="after diagcld"
[687]	2650	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2651	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2652	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2653	END IF
	2654	c
	2655	c Calculer l'humidite relative pour diagnostique
	2656	c
	2657	DO k = 1, klev
	2658	DO i = 1, klon
	2659	zx_t = t_seri(i,k)
	2660	IF (thermcep) THEN
	2661	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	2662	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	2663	zx_qs = MIN(0.5,zx_qs)
	2664	zcor = 1./(1.-retv*zx_qs)
	2665	zx_qs = zx_qs*zcor
	2666	ELSE
	2667	IF (zx_t.LT.t_coup) THEN
	2668	zx_qs = qsats(zx_t)/pplay(i,k)
	2669	ELSE
	2670	zx_qs = qsatl(zx_t)/pplay(i,k)
	2671	ENDIF
	2672	ENDIF
	2673	zx_rh(i,k) = q_seri(i,k)/zx_qs
	2674	zqsat(i,k)=zx_qs
	2675	ENDDO
	2676	ENDDO
[782]	2677
[687]	2678	cIM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
	2679	c equivalente a 2m (tpote) pour diagnostique
	2680	c
	2681	DO i = 1, klon
	2682	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
	2683	IF (thermcep) THEN
	2684	IF(zt2m(i).LT.RTT) then
	2685	Lheat=RLSTT
	2686	ELSE
	2687	Lheat=RLVTT
	2688	ENDIF
	2689	ELSE
	2690	IF (zt2m(i).LT.RTT) THEN
	2691	Lheat=RLSTT
	2692	ELSE
	2693	Lheat=RLVTT
	2694	ENDIF
	2695	ENDIF
	2696	tpote(i) = tpot(i)*
	2697	. EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
	2698	ENDDO
[524]	2699
	2700
	2701	#ifdef INCA
[766]	2702	call VTe(VTphysiq)
	2703	call VTb(VTinca)
[524]	2704	calday = FLOAT(julien) + gmtime
	2705
	2706	#ifdef INCA_AER
[766]	2707	call AEROSOL_METEO_CALC(calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs
	2708	& ,prfl,psfl,pctsrf,airephy,xjour,rlat,rlon,u10m,v10m)
[524]	2709	#endif
	2710
	2711	#ifdef INCAINFO
	2712	WRITE(lunout,*)'Appel CHEMHOOK_BEGIN ...'
	2713	#endif
[625]	2714
[782]	2715	zxsnow_dummy(:) = 0.0
	2716
[524]	2717	CALL chemhook_begin (calday,
[678]	2718	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
	2719	$ julien,
	2720	$ gmtime,
	2721	#endif
[593]	2722	$ pctsrf(1,1),
[524]	2723	$ rlat,
	2724	$ rlon,
	2725	$ airephy,
	2726	$ paprs,
	2727	$ pplay,
	2728	$ ycoefh,
	2729	$ pphi,
	2730	$ t_seri,
	2731	$ u,
	2732	$ v,
	2733	$ wo,
	2734	$ q_seri,
	2735	$ zxtsol,
[782]	2736	$ zxsnow_dummy,
[524]	2737	$ solsw,
[888]	2738	$ albsol1,
[524]	2739	$ rain_fall,
	2740	$ snow_fall,
	2741	$ itop_con,
	2742	$ ibas_con,
	2743	$ cldfra,
	2744	$ iim,
	2745	$ jjm,
[616]	2746	#ifdef INCA_AER
	2747	$ tr_seri,
	2748	$ ftsol,
	2749	$ paprs,
	2750	$ cdragh,
	2751	$ cdragm,
	2752	$ pctsrf,
	2753	$ pdtphys,
	2754	$ itap)
	2755	#else
	2756	$ tr_seri)
	2757	#endif
	2758
[625]	2759
[524]	2760	#ifdef INCAINFO
	2761	WRITE(lunout,*)'OK.'
	2762	#endif
[766]	2763	call VTe(VTinca)
	2764	call VTb(VTphysiq)
[524]	2765	#endif
	2766	c
	2767	c Calculer les parametres optiques des nuages et quelques
	2768	c parametres pour diagnostiques:
	2769	c
	2770	if (ok_newmicro) then
	2771	CALL newmicro (paprs, pplay,ok_newmicro,
	2772	. t_seri, cldliq, cldfra, cldtau, cldemi,
	2773	. cldh, cldl, cldm, cldt, cldq,
	2774	. flwp, fiwp, flwc, fiwc,
	2775	e ok_aie,
[864]	2776	#if defined(CPP_COUPLE) && defined(INCA_AER)
	2777	e ccm(:,:,1), ccm(:,:,2),
	2778	#else
[524]	2779	e sulfate, sulfate_pi,
[864]	2780	#endif
[524]	2781	e bl95_b0, bl95_b1,
	2782	s cldtaupi, re, fl)
	2783	else
	2784	CALL nuage (paprs, pplay,
	2785	. t_seri, cldliq, cldfra, cldtau, cldemi,
	2786	. cldh, cldl, cldm, cldt, cldq,
	2787	e ok_aie,
	2788	e sulfate, sulfate_pi,
	2789	e bl95_b0, bl95_b1,
	2790	s cldtaupi, re, fl)
	2791
	2792	endif
	2793	c
	2794	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
	2795	c
	2796	IF (MOD(itaprad,radpas).EQ.0) THEN
[782]	2797
[524]	2798	DO i = 1, klon
[888]	2799	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce)
	2800	. + falb1(i,is_lic) * pctsrf(i,is_lic)
	2801	. + falb1(i,is_ter) * pctsrf(i,is_ter)
	2802	. + falb1(i,is_sic) * pctsrf(i,is_sic)
	2803	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce)
	2804	. + falb2(i,is_lic) * pctsrf(i,is_lic)
	2805	. + falb2(i,is_ter) * pctsrf(i,is_ter)
	2806	. + falb2(i,is_sic) * pctsrf(i,is_sic)
[524]	2807	ENDDO
[766]	2808
	2809	if (mydebug) then
	2810	call writefield_phy('u_seri',u_seri,llm)
	2811	call writefield_phy('v_seri',v_seri,llm)
	2812	call writefield_phy('t_seri',t_seri,llm)
	2813	call writefield_phy('q_seri',q_seri,llm)
	2814	endif
	2815
[864]	2816
	2817	#if defined(CPP_COUPLE) && defined(INCA_AER)
	2818	CALL radlwsw_inca ! nouveau rayonnement (compatible Arpege-IFS)
	2819	e (kdlon,kflev,dist, rmu0, fract,
	2820	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
	2821	e wo,
	2822	e cldfra, cldemi, cldtau,
	2823	s heat,heat0,cool,cool0,radsol,albpla,
	2824	s topsw,toplw,solsw,sollw,
	2825	s sollwdown,
	2826	s topsw0,toplw0,solsw0,sollw0,
	2827	s lwdn0, lwdn, lwup0, lwup,
	2828	s swdn0, swdn, swup0, swup,
	2829	e ok_ade, ok_aie, ! new for aerosol radiative effects
	2830	e tau_inca, piz_inca, cg_inca, ! ="=
	2831	s topswad_inca, solswad_inca, ! ="=
	2832	s topswad0_inca, solswad0_inca, ! ="=
	2833	s topsw_inca, topsw0_inca,
	2834	s solsw_inca, solsw0_inca,
	2835	e cldtaupi, ! ="=
	2836	s topswai_inca, solswai_inca) ! ="=
	2837	#else
[524]	2838	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
	2839	e (dist, rmu0, fract,
[888]	2840	e paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri,
[524]	2841	e wo,
	2842	e cldfra, cldemi, cldtau,
	2843	s heat,heat0,cool,cool0,radsol,albpla,
	2844	s topsw,toplw,solsw,sollw,
	2845	s sollwdown,
	2846	s topsw0,toplw0,solsw0,sollw0,
	2847	s lwdn0, lwdn, lwup0, lwup,
	2848	s swdn0, swdn, swup0, swup,
	2849	e ok_ade, ok_aie, ! new for aerosol radiative effects
	2850	e tau_ae, piz_ae, cg_ae, ! ="=
	2851	s topswad, solswad, ! ="=
	2852	e cldtaupi, ! ="=
	2853	s topswai, solswai) ! ="=
[864]	2854	#endif
[524]	2855	itaprad = 0
	2856	ENDIF
	2857	itaprad = itaprad + 1
[879]	2858
	2859	if (iflag_radia.eq.0) then
	2860	print *,'--------------------------------------------------'
	2861	print *,'>>>> ATTENTION rayonnement desactive pour ce cas'
	2862	print *,'>>>> heat et cool mis a zero '
	2863	print *,'--------------------------------------------------'
	2864	heat=0.
	2865	cool=0.
	2866	endif
	2867
	2868
[524]	2869	c
	2870	c Ajouter la tendance des rayonnements (tous les pas)
	2871	c
	2872	DO k = 1, klev
	2873	DO i = 1, klon
	2874	t_seri(i,k) = t_seri(i,k)
	2875	. + (heat(i,k)-cool(i,k)) * dtime/86400.
	2876	ENDDO
	2877	ENDDO
[766]	2878	c
	2879	if (mydebug) then
	2880	call writefield_phy('u_seri',u_seri,llm)
	2881	call writefield_phy('v_seri',v_seri,llm)
	2882	call writefield_phy('t_seri',t_seri,llm)
	2883	call writefield_phy('q_seri',q_seri,llm)
	2884	endif
	2885
[687]	2886	cIM
	2887	IF (ip_ebil_phy.ge.2) THEN
[524]	2888	ztit='after rad'
[687]	2889	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	2890	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	2891	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
[687]	2892	call diagphy(airephy,ztit,ip_ebil_phy
[524]	2893	e , topsw, toplw, solsw, sollw, zero_v
	2894	e , zero_v, zero_v, zero_v, ztsol
	2895	e , d_h_vcol, d_qt, d_ec
	2896	s , fs_bound, fq_bound )
	2897	END IF
	2898	c
	2899	c
	2900	c Calculer l'hydrologie de la surface
	2901	c
	2902	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
	2903	c . agesno, ftsol,fqsurf,fsnow, ruis)
	2904	c
[782]	2905
[524]	2906	c
	2907	c Calculer le bilan du sol et la derive de temperature (couplage)
	2908	c
	2909	DO i = 1, klon
	2910	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
	2911	c a la demande de JLD
	2912	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
	2913	ENDDO
	2914	c
	2915	cmoddeblott(jan95)
	2916	c Appeler le programme de parametrisation de l'orographie
	2917	c a l'echelle sous-maille:
	2918	c
	2919	IF (ok_orodr) THEN
	2920	c
	2921	c selection des points pour lesquels le shema est actif:
	2922	igwd=0
	2923	DO i=1,klon
	2924	itest(i)=0
	2925	c IF ((zstd(i).gt.10.0)) THEN
	2926	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	2927	itest(i)=1
	2928	igwd=igwd+1
	2929	idx(igwd)=i
	2930	ENDIF
	2931	ENDDO
	2932	c igwdim=MAX(1,igwd)
	2933	c
	2934	CALL drag_noro(klon,klev,dtime,paprs,pplay,
	2935	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
	2936	e igwd,idx,itest,
	2937	e t_seri, u_seri, v_seri,
[644]	2938	cIM 141004 s zulow, zvlow, zustr, zvstr,
	2939	s zulow, zvlow, zustrdr, zvstrdr,
[524]	2940	s d_t_oro, d_u_oro, d_v_oro)
	2941	c
	2942	c ajout des tendances
[904]	2943	!-----------------------------------------------------------------------------------------
	2944	! ajout des tendances de la trainee de l'orographie
	2945	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
	2946	!-----------------------------------------------------------------------------------------
[524]	2947	c
	2948	ENDIF ! fin de test sur ok_orodr
	2949	c
[766]	2950	if (mydebug) then
	2951	call writefield_phy('u_seri',u_seri,llm)
	2952	call writefield_phy('v_seri',v_seri,llm)
	2953	call writefield_phy('t_seri',t_seri,llm)
	2954	call writefield_phy('q_seri',q_seri,llm)
	2955	endif
	2956
[524]	2957	IF (ok_orolf) THEN
	2958	c
	2959	c selection des points pour lesquels le shema est actif:
	2960	igwd=0
	2961	DO i=1,klon
	2962	itest(i)=0
	2963	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	2964	itest(i)=1
	2965	igwd=igwd+1
	2966	idx(igwd)=i
	2967	ENDIF
	2968	ENDDO
	2969	c igwdim=MAX(1,igwd)
	2970	c
	2971	CALL lift_noro(klon,klev,dtime,paprs,pplay,
	2972	e rlat,zmea,zstd,zpic,
	2973	e itest,
	2974	e t_seri, u_seri, v_seri,
[644]	2975	s zulow, zvlow, zustrli, zvstrli,
[524]	2976	s d_t_lif, d_u_lif, d_v_lif)
	2977	c
[904]	2978	!-----------------------------------------------------------------------------------------
	2979	! ajout des tendances de la portance de l'orographie
	2980	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
	2981	!-----------------------------------------------------------------------------------------
[524]	2982	c
	2983	ENDIF ! fin de test sur ok_orolf
	2984	c
[644]	2985	cIM cf. FLott BEG
	2986	C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
	2987
[766]	2988	if (mydebug) then
	2989	call writefield_phy('u_seri',u_seri,llm)
	2990	call writefield_phy('v_seri',v_seri,llm)
	2991	call writefield_phy('t_seri',t_seri,llm)
	2992	call writefield_phy('q_seri',q_seri,llm)
	2993	endif
	2994
[644]	2995	DO i = 1, klon
	2996	zustrph(i)=0.
	2997	zvstrph(i)=0.
	2998	ENDDO
	2999	DO k = 1, klev
	3000	DO i = 1, klon
	3001	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime*
	3002	c (paprs(i,k)-paprs(i,k+1))/rg
	3003	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime*
	3004	c (paprs(i,k)-paprs(i,k+1))/rg
	3005	ENDDO
	3006	ENDDO
	3007	c
	3008	cIM calcul composantes axiales du moment angulaire et couple des montagnes
	3009	c
[776]	3010	IF (is_sequential) THEN
[766]	3011
	3012	CALL aaam_bud (27,klon,klev,rjourvrai,gmtime,
	3013	C ra,rg,romega,
	3014	C rlat,rlon,pphis,
	3015	C zustrdr,zustrli,zustrph,
	3016	C zvstrdr,zvstrli,zvstrph,
	3017	C paprs,u,v,
	3018	C aam, torsfc)
	3019	ENDIF
[644]	3020	cIM cf. FLott END
[687]	3021	cIM
	3022	IF (ip_ebil_phy.ge.2) THEN
[524]	3023	ztit='after orography'
[687]	3024	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime
[524]	3025	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	3026	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3027	END IF
	3028	c
	3029	c
	3030	cAA
	3031	cAA Installation de l'interface online-offline pour traceurs
	3032	cAA
	3033	c====================================================================
	3034	c Calcul des tendances traceurs
	3035	c====================================================================
	3036	C
[658]	3037	call phytrac ( rnpb,
[625]	3038	I itap,
	3039	I julien,
	3040	I gmtime,
	3041	I debut,
	3042	I lafin,
[524]	3043	I nqmax-2,
[625]	3044	I nlon,
	3045	I nlev,
	3046	I dtime,
	3047	I u,
	3048	I v,
	3049	I t,
	3050	I paprs,
	3051	I pplay,
	3052	I pmfu,
	3053	I pmfd,
	3054	I pen_u,
	3055	I pde_u,
	3056	I pen_d,
	3057	I pde_d,
	3058	I ycoefh,
	3059	I fm_therm,
	3060	I entr_therm,
	3061	I yu1,
	3062	I yv1,
	3063	I ftsol,
	3064	I pctsrf,
	3065	I rlat,
	3066	I frac_impa,
	3067	I frac_nucl,
	3068	I rlon,
	3069	I presnivs,
	3070	I pphis,
	3071	I pphi,
[888]	3072	I albsol1,
[625]	3073	I qx(1,1,1),
	3074	I rhcl,
	3075	I cldfra,
	3076	I rneb,
	3077	I diafra,
	3078	I cldliq,
	3079	I itop_con,
[524]	3080	I ibas_con,
[625]	3081	I pmflxr,
	3082	I pmflxs,
	3083	I prfl,
	3084	I psfl,
	3085	I da,
	3086	I phi,
	3087	I mp,
	3088	I upwd,
	3089	I dnwd,
	3090	#ifdef INCA
[524]	3091	I flxmass_w,
[864]	3092	#if defined(INCA_AER) && defined(CPP_COUPLE)
	3093	I tau_inca,
	3094	I piz_inca,
	3095	I cg_inca,
	3096	I ccm,
	3097	I rfname,
[524]	3098	#endif
[864]	3099	#endif
[524]	3100	O tr_seri)
	3101
	3102	IF (offline) THEN
	3103
[541]	3104	print*,'Attention on met a 0 les thermiques pour phystoke'
[524]	3105	call phystokenc (
	3106	I nlon,nlev,pdtphys,rlon,rlat,
	3107	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
[541]	3108	I fm_therm,entr_therm,
[524]	3109	I ycoefh,yu1,yv1,ftsol,pctsrf,
	3110	I frac_impa, frac_nucl,
	3111	I pphis,airephy,dtime,itap)
	3112
	3113
	3114	ENDIF
	3115
	3116	c
	3117	c Calculer le transport de l'eau et de l'energie (diagnostique)
	3118	c
	3119	CALL transp (paprs,zxtsol,
	3120	e t_seri, q_seri, u_seri, v_seri, zphi,
	3121	s ve, vq, ue, uq)
	3122	c
[687]	3123	cIM global posePB BEG
	3124	IF(1.EQ.0) THEN
[524]	3125	c
[644]	3126	CALL transp_lay (paprs,zxtsol,
	3127	e t_seri, q_seri, u_seri, v_seri, zphi,
	3128	s ve_lay, vq_lay, ue_lay, uq_lay)
[524]	3129	c
[687]	3130	ENDIF !(1.EQ.0) THEN
	3131	cIM global posePB END
[644]	3132	c Accumuler les variables a stocker dans les fichiers histoire:
[524]	3133	c
	3134	c+jld ec_conser
	3135	DO k = 1, klev
	3136	DO i = 1, klon
	3137	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
	3138	d_t_ec(i,k)=0.5/ZRCPD
	3139	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
	3140	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
	3141	d_t_ec(i,k) = d_t_ec(i,k)/dtime
	3142	END DO
	3143	END DO
	3144	c-jld ec_conser
[687]	3145	cIM
	3146	IF (ip_ebil_phy.ge.1) THEN
[524]	3147	ztit='after physic'
[687]	3148	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime
[524]	3149	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
	3150	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3151	C Comme les tendances de la physique sont ajoute dans la dynamique,
	3152	C on devrait avoir que la variation d'entalpie par la dynamique
	3153	C est egale a la variation de la physique au pas de temps precedent.
	3154	C Donc la somme de ces 2 variations devrait etre nulle.
[687]	3155	call diagphy(airephy,ztit,ip_ebil_phy
[524]	3156	e , topsw, toplw, solsw, sollw, sens
	3157	e , evap, rain_fall, snow_fall, ztsol
	3158	e , d_h_vcol, d_qt, d_ec
	3159	s , fs_bound, fq_bound )
	3160	C
	3161	d_h_vcol_phy=d_h_vcol
	3162	C
	3163	END IF
	3164	C
	3165	c=======================================================================
	3166	c SORTIES
	3167	c=======================================================================
	3168
[644]	3169	cIM Interpolation sur les niveaux de pression du NMC
	3170	c -------------------------------------------------
[524]	3171	c
[644]	3172	#include "calcul_STDlev.h"
[524]	3173	c
	3174	c slp sea level pressure
	3175	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
	3176	c
	3177	ccc prw = eau precipitable
	3178	DO i = 1, klon
	3179	prw(i) = 0.
	3180	DO k = 1, klev
	3181	prw(i) = prw(i) +
	3182	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
	3183	ENDDO
	3184	ENDDO
	3185	c
[644]	3186	cIM initialisation + calculs divers diag AMIP2
[524]	3187	c
[644]	3188	#include "calcul_divers.h"
	3189	c
[655]	3190	#ifdef INCA
[766]	3191	call VTe(VTphysiq)
	3192	call VTb(VTinca)
[655]	3193	#ifdef INCAINFO
	3194	WRITE(lunout,*)'Appel CHEMHOOK_END ...'
	3195	#endif
	3196	CALL chemhook_end (calday,
	3197	$ dtime,
	3198	$ pplay,
	3199	$ t_seri,
	3200	$ tr_seri,
	3201	$ nbtr,
	3202	$ paprs,
	3203	$ q_seri,
	3204	$ annee_ref,
	3205	$ day_ini,
[791]	3206	$ airephy,
[655]	3207	#ifdef INCA_AER
	3208	$ xjour,
	3209	$ pphi,
	3210	$ pphis,
[766]	3211	$ zx_rh)
[655]	3212	#else
	3213	$ xjour)
	3214	#endif
	3215	#ifdef INCAINFO
	3216	WRITE(lunout,*)'OK.'
	3217	#endif
[766]	3218	call VTe(VTinca)
	3219	call VTb(VTphysiq)
[655]	3220	#endif
	3221
[524]	3222	c=============================================================
	3223	c
	3224	c Convertir les incrementations en tendances
	3225	c
[766]	3226	if (mydebug) then
	3227	call writefield_phy('u_seri',u_seri,llm)
	3228	call writefield_phy('v_seri',v_seri,llm)
	3229	call writefield_phy('t_seri',t_seri,llm)
	3230	call writefield_phy('q_seri',q_seri,llm)
	3231	endif
	3232
[524]	3233	DO k = 1, klev
	3234	DO i = 1, klon
	3235	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	3236	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	3237	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
	3238	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
	3239	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
	3240	ENDDO
	3241	ENDDO
	3242	c
	3243	IF (nqmax.GE.3) THEN
	3244	DO iq = 3, nqmax
	3245	DO k = 1, klev
	3246	DO i = 1, klon
	3247	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
	3248	ENDDO
	3249	ENDDO
	3250	ENDDO
	3251	ENDIF
	3252	c
[644]	3253	cIM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
[687]	3254	cIM global posePB#include "write_bilKP_ins.h"
	3255	cIM global posePB#include "write_bilKP_ave.h"
[644]	3256	c
[524]	3257	c Sauvegarder les valeurs de t et q a la fin de la physique:
	3258	c
	3259	DO k = 1, klev
	3260	DO i = 1, klon
	3261	t_ancien(i,k) = t_seri(i,k)
	3262	q_ancien(i,k) = q_seri(i,k)
	3263	ENDDO
	3264	ENDDO
	3265	c
[879]	3266	!==========================================================================
	3267	! Sorties des tendances pour un point particulier
	3268	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
	3269	! pour le debug
	3270	! La valeur de igout est attribuee plus haut dans le programme
	3271	!==========================================================================
	3272
[942]	3273	if (prt_level.ge.1) then
[879]	3274	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	3275	write(lunout,*)
[929]	3276	s 'nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys pct tlos'
[879]	3277	write(lunout,*)
[929]	3278	s nlon,nlev,nqmax,debut,lafin,rjourvrai,gmtime,pdtphys,
[930]	3279	s pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce),
[929]	3280	s pctsrf(igout,is_sic)
[879]	3281	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
	3282	do k=1,nlev
	3283	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k),
	3284	s d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k),
	3285	s d_t_eva(igout,k)
	3286	enddo
	3287	write(lunout,*) 'cool,heat'
	3288	do k=1,nlev
	3289	write(lunout,*) cool(igout,k),heat(igout,k)
	3290	enddo
	3291
	3292	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
	3293	do k=1,nlev
	3294	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k),
	3295	s d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
	3296	enddo
	3297
	3298	write(lunout,*) 'd_ps ',d_ps(igout)
	3299	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
	3300	do k=1,nlev
	3301	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k),
	3302	s d_qx(igout,k,1),d_qx(igout,k,2)
	3303	enddo
	3304	endif
	3305
	3306	!==========================================================================
	3307
	3308	c============================================================
	3309	c Calcul de la temperature potentielle
	3310	c============================================================
	3311	DO k = 1, klev
	3312	DO i = 1, klon
	3313	theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
	3314	ENDDO
	3315	ENDDO
	3316	c
	3317
[524]	3318	c 22.03.04 BEG
	3319	c=============================================================
	3320	c Ecriture des sorties
	3321	c=============================================================
	3322	#ifdef CPP_IOIPSL
[782]	3323
	3324	c Recupere des varibles calcule dans differents modules
	3325	c pour ecriture dans histxxx.nc
[524]	3326
[888]	3327	! Get some variables from module fonte_neige_mod
[782]	3328	CALL fonte_neige_get_vars(pctsrf,
	3329	. zxfqcalving, zxfqfonte, zxffonte)
	3330
	3331	IF (ocean == 'slab') THEN
[888]	3332	! Get some variables from module ocean_slab_mod
[782]	3333	CALL ocean_slab_get_vars(tslab, seaice, fluxo, fluxg)
	3334	ELSEIF (ocean == 'couple') THEN
[888]	3335	! Get some variables from module ocean_cpl_mod
[782]	3336	CALL ocean_cpl_get_vars(fluxo, fluxg)
	3337	ELSE
[888]	3338	! Get some variables from module ocean_forced_mod
[782]	3339	CALL ocean_forced_get_vars(fluxo, fluxg)
	3340	ENDIF
	3341
[524]	3342
[909]	3343	c Commente par abderrahmane le 11 2 08
	3344	c#ifdef histhf
	3345	c#include "write_histhf.h"
	3346	c#endif
[524]	3347
[909]	3348	c#ifdef histday
	3349	c#include "write_histday.h"
	3350	c#endif
[524]	3351
[909]	3352	c#ifdef histmth
	3353	c#include "write_histmth.h"
	3354	c#endif
[524]	3355
[909]	3356	c#ifdef histins
	3357	c#include "write_histins.h"
	3358	c#endif
	3359
	3360	#include "phys_output_write.h"
	3361
[524]	3362	#ifdef histISCCP
	3363	#include "write_histISCCP.h"
	3364	#endif
	3365
	3366	#ifdef histmthNMC
	3367	#include "write_histmthNMC.h"
	3368	#endif
	3369
[687]	3370	#include "write_histday_seri.h"
	3371
	3372	#include "write_paramLMDZ_phy.h"
	3373
[524]	3374	#endif
	3375
	3376	c 22.03.04 END
	3377	c
	3378	c====================================================================
	3379	c Si c'est la fin, il faut conserver l'etat de redemarrage
	3380	c====================================================================
	3381	c
[782]	3382
	3383
[524]	3384	IF (lafin) THEN
	3385	itau_phy = itau_phy + itap
[782]	3386	CALL phyredem ("restartphy.nc",dtime,radpas,ocean,
	3387	. rlat, rlon, pctsrf, ftsol,
[888]	3388	. falb1, falb2, rain_fall,
[782]	3389	. snow_fall,
	3390	. solsw, sollw,
	3391	. radsol,
[524]	3392	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
[878]	3393	. t_ancien, q_ancien, rnebcon, ratqs, clwcon,
[937]	3394	. pbl_tke, zmax0, f0, ema_work1, ema_work2)
[904]	3395	open(97,form="unformatted",file="finbin")
	3396	write(97) u_seri,v_seri,t_seri,q_seri
	3397	close(97)
[524]	3398	ENDIF
	3399
	3400
	3401	RETURN
	3402	END
	3403	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
	3404	IMPLICIT none
	3405	c
	3406	c Calculer et imprimer l'eau totale. A utiliser pour verifier
	3407	c la conservation de l'eau
	3408	c
	3409	#include "YOMCST.h"
	3410	INTEGER klon,klev
	3411	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
	3412	REAL aire(klon)
	3413	REAL qtotal, zx, qcheck
	3414	INTEGER i, k
	3415	c
	3416	zx = 0.0
	3417	DO i = 1, klon
	3418	zx = zx + aire(i)
	3419	ENDDO
	3420	qtotal = 0.0
	3421	DO k = 1, klev
	3422	DO i = 1, klon
	3423	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
	3424	. *(paprs(i,k)-paprs(i,k+1))/RG
	3425	ENDDO
	3426	ENDDO
	3427	c
	3428	qcheck = qtotal/zx
	3429	c
	3430	RETURN
	3431	END
	3432	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	3433	IMPLICIT none
	3434	c
	3435	c Tranformer une variable de la grille physique a
	3436	c la grille d'ecriture
	3437	c
	3438	INTEGER nfield,nlon,iim,jjmp1, jjm
	3439	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
	3440	c
	3441	INTEGER i, n, ig
	3442	c
	3443	jjm = jjmp1 - 1
	3444	DO n = 1, nfield
	3445	DO i=1,iim
	3446	ecrit(i,n) = fi(1,n)
	3447	ecrit(i+jjm*iim,n) = fi(nlon,n)
	3448	ENDDO
	3449	DO ig = 1, nlon - 2
	3450	ecrit(iim+ig,n) = fi(1+ig,n)
	3451	ENDDO
	3452	ENDDO
	3453	RETURN
	3454	END

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

Download in other formats:

Original Format