Context Navigation

physiq.F @ 717

Last change on this file since 717 was 717, checked in by Laurent Fairhead, 18 years ago
Modifs pour compilation sur Brodie LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 117.1 KB

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

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

Download in other formats:

Original Format