Context Navigation

physiq.F @ 969

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

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

Download in other formats:

Original Format