Context Navigation

physiq.F @ 3947

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

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

Download in other formats:

Original Format