Context Navigation

physiq.F @ 5319

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

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

Download in other formats:

Original Format