Context Navigation

physiq.F @ 974

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

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

Download in other formats:

Original Format