Context Navigation

physiq.F @ 936

Last change on this file since 936 was 933, checked in by lmdzadmin, 17 years ago
Bug ecrit_hf2mth interpolation niveaux standard de pression, cas guide IM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 124.5 KB

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

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

Download in other formats:

Original Format