Context Navigation

physiq.F @ 715

Last change on this file since 715 was 702, checked in by lmdzadmin, 19 years ago
Frequence sortie traceurs (ecrit_tra) en jours. MAF CVS
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 102.1 KB

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

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

Download in other formats:

Original Format