Context Navigation

physiq.F @ 735

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

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

Download in other formats:

Original Format