Context Navigation

physiq.F @ 5434

Last change on this file since 5434 was 954, checked in by lsce, 17 years ago

ACo+JG

Ajout du flag aerosol_couple :
false=lecture des sulfates dans un fichier(par defaut) - configuration existante
true=calcul des aerosol par INCA

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 108.1 KB

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

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

Download in other formats:

Original Format