Context Navigation

physiq.F @ 1098

Last change on this file since 1098 was 864, checked in by lsce, 17 years ago
ACo : ajout clef CPP (COUPLE + AER) pour le modele couple avec les aerosol (IPSL_ESM_V1)
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 112.7 KB

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

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

Download in other formats:

Original Format