Context Navigation

physiq.F @ 674

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

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

Download in other formats:

Original Format