Context Navigation

physiq.F90 @ 354

Last change on this file since 354 was 186, checked in by lfita, 10 years ago
Removing checking printings
File size: 165.2 KB

Rev	Line
[1]	1	! $Id: physiq.F 1795 2013-07-18 08:20:28Z emillour $
	2	!c#define IO_DEBUG
	3
	4	SUBROUTINE physiq (nlon,nlev, &
	5	& debut,lafin,jD_cur, jH_cur,pdtphys, &
	6	& paprs,pplay,pphi,pphis,presnivs,clesphy0, &
	7	& u,v,t,qx, &
	8	& flxmass_w, &
	9	& d_u, d_v, d_t, d_qx, d_ps &
	10	& , dudyn &
	11	& , PVteta)
	12
	13	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
	14	& histwrite, ju2ymds, ymds2ju, ioget_year_len
	15	USE comgeomphy
	16	USE phys_cal_mod
	17	USE write_field_phy
	18	USE dimphy
	19	USE infotrac
	20	USE mod_grid_phy_lmdz
	21	USE mod_phys_lmdz_para
	22	USE iophy
	23	USE misc_mod, mydebug=>debug
	24	USE vampir
	25	!L. Fita, LMD. November 2013
	26	! USE pbl_surface_mod, ONLY : pbl_surface
	27	! USE pbl_surface_mod, ONLY : pbl_surface, pbl_surface_init
	28	USE pbl_surface_mod
	29	USE change_srf_frac_mod
	30	USE surface_data, ONLY : type_ocean, ok_veget
	31	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
	32	USE phys_state_var_mod ! Variables sauvegardees de la physique
	33	USE phys_output_var_mod ! Variables pour les ecritures des sorties
	34	USE fonte_neige_mod, ONLY : fonte_neige_get_vars, fonte_neige_init
	35	USE phys_output_mod
	36	USE phys_output_ctrlout_mod
	37	USE iophy
	38	use open_climoz_m, only: open_climoz ! ozone climatology from a file
	39	use regr_pr_av_m, only: regr_pr_av
	40	use netcdf95, only: nf95_close
	41	!IM for NMC files
	42	!c use netcdf, only: nf90_fill_real
	43	use netcdf
	44	use mod_phys_lmdz_mpi_data, only: is_mpi_root
	45	USE aero_mod
	46	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
	47	use conf_phys_m, only: conf_phys
	48	use radlwsw_m, only: radlwsw
	49	USE control_mod
	50	#ifdef REPROBUS
	51	USE CHEM_REP, ONLY : Init_chem_rep_xjour
	52	#endif
	53	USE indice_sol_mod
	54
	55	!IM stations CFMIP
	56	USE CFMIP_point_locations
	57
	58	!L. Fita. July 2013
	59	USE orografi
	60	USE orografi_strato_mod
	61	USE diagphy_mod
	62	! USE orbite_mod, ONLY : angle
	63	USE phyaqua_mod, ONLY: zenang_an
	64	USE nuage_mod, ONLY : diagcld1
	65
	66	! L. Fita, LMD. November 2013.
	67	!! USE wrf_lmdz_mod
	68	!! USE module_lmdz_variables, ONLY: neige_initialize, limit_initialize
	69	!! USE physiq_limit_variables_mod
	70	!! USE physiq_limit_variables_mod
	71	USE lmdz_wrf_variables_mod
	72	! L. Fita, LMD. January 2014
	73	! Defining with SAVE statement all that variables that should appear in the output
	74	USE output_lmdz_NOmodule
	75
	76	IMPLICIT none
	77	!>======================================================================
	78	!!
	79	!! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	80	!!
	81	!! Objet: Moniteur general de la physique du modele
	82	!!AA Modifications quant aux traceurs :
	83	!!AA - uniformisation des parametrisations ds phytrac
	84	!!AA - stockage des moyennes des champs necessaires
	85	!!AA en mode traceur off-line
	86	!!======================================================================
	87	!! CLEFS CPP POUR LES IO
	88	!! =====================
	89	#define histNMC
	90	!c#define histISCCP
	91	!!======================================================================
	92	!! modif ( P. Le Van , 12/10/98 )
	93	!!
	94	!! Arguments:
	95	!!
	96	!! nlon----input-I-nombre de points horizontaux
	97	!! nlev----input-I-nombre de couches verticales, doit etre egale a klev
	98	!! debut---input-L-variable logique indiquant le premier passage
	99	!! lafin---input-L-variable logique indiquant le dernier passage
	100	!! jD_cur -R-jour courant a l'appel de la physique (jour julien)
	101	!! jH_cur -R-heure courante a l'appel de la physique (jour julien)
	102	!! pdtphys-input-R-pas d'integration pour la physique (seconde)
	103	!! paprs---input-R-pression pour chaque inter-couche (en Pa)
	104	!! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	105	!! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
	106	!! pphis---input-R-geopotentiel du sol
	107	!! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
	108	!! u-------input-R-vitesse dans la direction X (de O a E) en m/s
	109	!! v-------input-R-vitesse Y (de S a N) en m/s
	110	!! t-------input-R-temperature (K)
	111	!! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
	112	!! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
	113	!! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
	114	!! flxmass_w -input-R- flux de masse verticale
	115	!! d_u-----output-R-tendance physique de "u" (m/s/s)
	116	!! d_v-----output-R-tendance physique de "v" (m/s/s)
	117	!! d_t-----output-R-tendance physique de "t" (K/s)
	118	!! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
	119	!! d_ps----output-R-tendance physique de la pression au sol
	120	!!IM
	121	!! PVteta--output-R-vorticite potentielle a des thetas constantes
	122	!!======================================================================
	123	#include "dimensions.h"
	124	integer jjmp1
	125	parameter (jjmp1=jjm+1-1/jjm)
	126	integer iip1
	127	parameter (iip1=iim+1)
	128
	129	#include "regdim.h"
	130	#include "dimsoil.h"
	131	#include "clesphys.h"
	132	#include "temps.h"
	133	#include "iniprint.h"
	134	#include "thermcell.h"
	135	!c======================================================================
	136	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
	137	PARAMETER (ok_cvl=.TRUE.)
	138	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
	139	PARAMETER (ok_gust=.FALSE.)
	140	integer iflag_radia ! active ou non le rayonnement (MPL)
	141	save iflag_radia
	142	!$OMP THREADPRIVATE(iflag_radia)
	143	!c======================================================================
	144	LOGICAL check ! Verifier la conservation du modele en eau
	145	PARAMETER (check=.FALSE.)
	146	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
	147	PARAMETER (ok_stratus=.FALSE.)
	148	!c======================================================================
	149	REAL amn, amx
	150	INTEGER igout
	151	!c======================================================================
	152	!c Clef controlant l'activation du cycle diurne:
	153	!ccc LOGICAL cycle_diurne
	154	!ccc PARAMETER (cycle_diurne=.FALSE.)
	155	!c======================================================================
	156	!c Modele thermique du sol, a activer pour le cycle diurne:
	157	!ccc LOGICAL soil_model
	158	!ccc PARAMETER (soil_model=.FALSE.)
	159	!c======================================================================
	160	!c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
	161	!c le calcul du rayonnement est celle apres la precipitation des nuages.
	162	!c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
	163	!c la condensation et la precipitation. Cette cle augmente les impacts
	164	!c radiatifs des nuages.
	165	!ccc LOGICAL new_oliq
	166	!ccc PARAMETER (new_oliq=.FALSE.)
	167	!c======================================================================
	168	!c Clefs controlant deux parametrisations de l'orographie:
	169	!cc LOGICAL ok_orodr
	170	!ccc PARAMETER (ok_orodr=.FALSE.)
	171	!ccc LOGICAL ok_orolf
	172	!ccc PARAMETER (ok_orolf=.FALSE.)
	173	!c======================================================================
	174	LOGICAL ok_journe ! sortir le fichier journalier
	175	save ok_journe
	176	!$OMP THREADPRIVATE(ok_journe)
	177	!c
	178	LOGICAL ok_mensuel ! sortir le fichier mensuel
	179	save ok_mensuel
	180	!$OMP THREADPRIVATE(ok_mensuel)
	181	!c
	182	LOGICAL ok_instan ! sortir le fichier instantane
	183	save ok_instan
	184	!$OMP THREADPRIVATE(ok_instan)
	185	!c
	186	LOGICAL ok_LES ! sortir le fichier LES
	187	save ok_LES
	188	!$OMP THREADPRIVATE(ok_LES)
	189	!c
	190	LOGICAL callstats ! sortir le fichier stats
	191	save callstats
	192	!$OMP THREADPRIVATE(callstats)
	193	!c
	194	LOGICAL ok_region ! sortir le fichier regional
	195	PARAMETER (ok_region=.FALSE.)
	196	!c======================================================================
	197	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	198	! real weak_inversion(klon),dthmin(klon)
	199	real seuil_inversion
	200	save seuil_inversion
	201	!$OMP THREADPRIVATE(seuil_inversion)
	202	integer iflag_ratqs
	203	save iflag_ratqs
	204	!$OMP THREADPRIVATE(iflag_ratqs)
	205	real facteur
	206
	207	REAL zz,znum,zden
	208	REAL wmax_th(klon)
	209	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	210	! REAL zmax_th(klon)
	211	REAL tau_overturning_th(klon)
	212
	213	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	214	! integer lmax_th(klon)
	215	integer limbas(klon)
	216	real ratqscth(klon,klev)
	217	real ratqsdiff(klon,klev)
	218	real zqsatth(klon,klev)
	219
	220	!c======================================================================
	221	!c
	222	INTEGER ivap ! indice de traceurs pour vapeur d'eau
	223	PARAMETER (ivap=1)
	224	INTEGER iliq ! indice de traceurs pour eau liquide
	225	PARAMETER (iliq=2)
	226
	227	!c
	228	!c
	229	!c Variables argument:
	230	!c
	231	INTEGER nlon
	232	INTEGER nlev
	233	REAL, intent(in):: jD_cur, jH_cur
	234
	235	REAL pdtphys
	236	LOGICAL debut, lafin
	237	REAL paprs(klon,klev+1)
	238	REAL pplay(klon,klev)
	239	REAL pphi(klon,klev)
	240	REAL pphis(klon)
	241	REAL presnivs(klev)
	242	REAL znivsig(klev)
	243	real pir
	244
	245	REAL u(klon,klev)
	246	REAL v(klon,klev)
	247	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	248	! REAL t(klon,klev),theta(klon,klev),thetal(klon,klev)
	249	REAL t(klon,klev),thetal(klon,klev)
	250
	251	!c thetal: ligne suivante a decommenter si vous avez les fichiers MPL 20130625
	252	!c fth_fonctions.F90 et parkind1.F90
	253	!c sinon thetal=theta
	254	!c REAL fth_thetae,fth_thetav,fth_thetal
	255	REAL qx(klon,klev,nqtot)
	256	REAL flxmass_w(klon,klev)
	257	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	258	! REAL omega(klon,klev) ! vitesse verticale en Pa/s
	259	REAL d_u(klon,klev)
	260	REAL d_v(klon,klev)
	261	REAL d_t(klon,klev)
	262	REAL d_qx(klon,klev,nqtot)
	263	REAL d_ps(klon)
	264	! Variables pour le transport convectif
	265	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
	266	! Variables pour le lessivage convectif
	267	! RomP >>>
	268	real phi2(klon,klev,klev)
	269	real d1a(klon,klev),dam(klon,klev)
	270	real ev(klon,klev),ep(klon,klev)
	271	real clw(klon,klev),elij(klon,klev,klev)
	272	real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
	273	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	274	! real wdtrainA(klon,klev),wdtrainM(klon,klev)
	275
	276	! RomP <<<
	277	!IM definition dynamique o_trac dans phys_output_open
	278	! type(ctrl_out) :: o_trac(nqtot)
	279	!c
	280	!IM Amip2 PV a theta constante
	281	!c
	282	INTEGER nbteta
	283	PARAMETER(nbteta=3)
	284	CHARACTER*3 ctetaSTD(nbteta)
	285	DATA ctetaSTD/'350','380','405'/
	286	SAVE ctetaSTD
	287	!$OMP THREADPRIVATE(ctetaSTD)
	288	REAL rtetaSTD(nbteta)
	289	DATA rtetaSTD/350., 380., 405./
	290	SAVE rtetaSTD
	291	!$OMP THREADPRIVATE(rtetaSTD)
	292	!c
	293	REAL PVteta(klon,nbteta)
	294	REAL zx_tmp_3dte(iim,jjmp1,nbteta)
	295	!c
	296	!MI Amip2 PV a theta constante
	297
	298	!cym INTEGER klevp1, klevm1
	299	!cym PARAMETER(klevp1=klev+1,klevm1=klev-1)
	300	!cym#include "raddim.h"
	301	!c
	302	!c
	303	!IM Amip2
	304	!c variables a une pression donnee
	305	!c
	306	#include "declare_STDlev.h"
	307	!c
	308	CHARACTER*4 bb2
	309	CHARACTER*2 bb3
	310	!c
	311	#include "radopt.h"
	312	!c
	313	!c
	314	!c prw: precipitable water
	315	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	316	! real prw(klon)
	317
	318	REAL convliq(klon,klev) ! eau liquide nuageuse convective
	319	REAL convfra(klon,klev) ! fraction nuageuse convective
	320
	321	REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
	322	REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
	323	REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
	324	REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
	325
	326	INTEGER linv, kp1
	327	!c flwp, fiwp = Liquid Water Path & Ice Water Path (kg/m2)
	328	!c flwc, fiwc = Liquid Water Content & Ice Water Content (kg/kg)
	329	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	330	! REAL flwp(klon), fiwp(klon)
	331	! REAL flwc(klon,klev), fiwc(klon,klev)
	332	REAL flwp_c(klon), fiwp_c(klon)
	333	REAL flwc_c(klon,klev), fiwc_c(klon,klev)
	334	REAL flwp_s(klon), fiwp_s(klon)
	335	REAL flwc_s(klon,klev), fiwc_s(klon,klev)
	336
	337	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	338	! REAL evap_pot(klon,nbsrf)
	339
	340	!IM ISCCP simulator v3.4
	341	!c dans clesphys.h top_height, overlap
	342	!cv3.4
	343	INTEGER debug, debugcol
	344	!cym INTEGER npoints
	345	!cym PARAMETER(npoints=klon)
	346	!c
	347	INTEGER sunlit(klon) !sunlit=1 if day; sunlit=0 if night
	348	INTEGER nregISCtot
	349	PARAMETER(nregISCtot=1)
	350	!c
	351	!c imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
	352	!c y compris pour 1 point
	353	!c imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
	354	!c jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
	355	INTEGER imin_debut, nbpti
	356	INTEGER jmin_debut, nbptj
	357	!IM parametres ISCCP BEG
	358	INTEGER nbapp_isccp
	359	! INTEGER nbapp_isccp,isccppas
	360	! PARAMETER(isccppas=6) !appel du simulateurs tous les 6pas de temps de la physique
	361	! !i.e. toutes les 3 heures
	362	INTEGER n
	363	INTEGER ifreq_isccp(napisccp), freqin_pdt(napisccp)
	364	DATA ifreq_isccp/3/
	365	SAVE ifreq_isccp
	366	!$OMP THREADPRIVATE(ifreq_isccp)
	367	CHARACTER*5 typinout(napisccp)
	368	DATA typinout/'i3od'/
	369	SAVE typinout
	370	!$OMP THREADPRIVATE(typinout)
	371	!IM verif boxptop BEG
	372	CHARACTER*1 verticaxe(napisccp)
	373	DATA verticaxe/'1'/
	374	SAVE verticaxe
	375	!$OMP THREADPRIVATE(verticaxe)
	376	!IM verif boxptop END
	377	INTEGER nvlev(napisccp)
	378	!c INTEGER nvlev
	379	REAL t1, aa
	380	REAL seed_re(klon,napisccp)
	381	!cym !!!! A voir plus tard
	382	!cym INTEGER iphy(iim,jjmp1)
	383	!IM parametres ISCCP END
	384	!c
	385	!c ncol = nb. de sous-colonnes pour chaque maille du GCM
	386	!c ncolmx = No. max. de sous-colonnes pour chaque maille du GCM
	387	!c INTEGER ncol(napisccp), ncolmx, seed(klon,napisccp)
	388	INTEGER,SAVE :: ncol(napisccp)
	389	!$OMP THREADPRIVATE(ncol)
	390	INTEGER ncolmx, seed(klon,napisccp)
	391	REAL nbsunlit(nregISCtot,klon,napisccp) !nbsunlit : moyenne de sunlit
	392	!c PARAMETER(ncolmx=1500)
	393	PARAMETER(ncolmx=300)
	394	!c
	395	!IM verif boxptop BEG
	396	REAL vertlev(ncolmx,napisccp)
	397	!IM verif boxptop END
	398	!c
	399	REAL,SAVE :: tautab_omp(0:255),tautab(0:255)
	400	INTEGER,SAVE :: invtau_omp(-20:45000),invtau(-20:45000)
	401	!$OMP THREADPRIVATE(tautab,invtau)
	402	REAL emsfc_lw
	403	PARAMETER(emsfc_lw=0.99)
	404	!c REAL ran0 ! type for random number fuction
	405	!c
	406	REAL cldtot(klon,klev)
	407	!c variables de haut en bas pour le simulateur ISCCP
	408	REAL dtau_s(klon,klev) !tau nuages startiformes
	409	REAL dtau_c(klon,klev) !tau nuages convectifs
	410	REAL dem_s(klon,klev) !emissivite nuages startiformes
	411	REAL dem_c(klon,klev) !emissivite nuages convectifs
	412	!c
	413	!c variables de haut en bas pour le simulateur ISCCP
	414	REAL pfull(klon,klev)
	415	REAL phalf(klon,klev+1)
	416	REAL qv(klon,klev)
	417	REAL cc(klon,klev)
	418	REAL conv(klon,klev)
	419	REAL dtau_sH2B(klon,klev)
	420	REAL dtau_cH2B(klon,klev)
	421	REAL at(klon,klev)
	422	REAL dem_sH2B(klon,klev)
	423	REAL dem_cH2B(klon,klev)
	424	!c
	425	INTEGER kmax, lmax, lmax3
	426	PARAMETER(kmax=8, lmax=8, lmax3=3)
	427	INTEGER kmaxm1, lmaxm1
	428	PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
	429	INTEGER iimx7, jjmx7, jjmp1x7
	430	PARAMETER(iimx7=iimkmaxm1, jjmx7=jjmlmaxm1, &
	431	&jjmp1x7=jjmp1*lmaxm1)
	432	!c
	433	!c output from ISCCP simulator
	434	REAL fq_isccp(klon,kmaxm1,lmaxm1,napisccp)
	435	REAL fq_is_true(klon,kmaxm1,lmaxm1,napisccp)
	436	REAL totalcldarea(klon,napisccp)
	437	REAL meanptop(klon,napisccp)
	438	REAL meantaucld(klon,napisccp)
	439	REAL boxtau(klon,ncolmx,napisccp)
	440	REAL boxptop(klon,ncolmx,napisccp)
	441	REAL zx_tmp_fi3d_bx(klon,ncolmx)
	442	REAL zx_tmp_3d_bx(iim,jjmp1,ncolmx)
	443	!c
	444	REAL cld_fi3d(klon,lmax3)
	445	REAL cld_3d(iim,jjmp1,lmax3)
	446	!c
	447	INTEGER iw, iwmax
	448	REAL wmin, pas_w
	449	!c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
	450	!IM 051005 PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
	451	PARAMETER(wmin=-100.,pas_w=10.,iwmax=20)
	452	REAL o500(klon)
	453	!c
	454
	455	!c sorties ISCCP
	456
	457	integer nid_isccp
	458	save nid_isccp
	459	!$OMP THREADPRIVATE(nid_isccp)
	460
	461	REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
	462	DATA zx_tau/0.0, 0.3, 1.3, 3.6, 9.4, 23., 60./
	463	SAVE zx_tau
	464	DATA zx_pc/180., 310., 440., 560., 680., 800., 1000./
	465	SAVE zx_pc
	466	!$OMP THREADPRIVATE(zx_tau,zx_pc)
	467	!c cldtopres pression au sommet des nuages
	468	REAL cldtopres(lmaxm1), cldtopres3(lmax3)
	469	DATA cldtopres/180., 310., 440., 560., 680., 800., 1000./
	470	DATA cldtopres3/440., 680., 1000./
	471	SAVE cldtopres,cldtopres3
	472	!$OMP THREADPRIVATE(cldtopres,cldtopres3)
	473	!IM 051005 BEG
	474	INTEGER komega, nhoriRD
	475
	476	!c taulev: numero du niveau de tau dans les sorties ISCCP
	477	CHARACTER *4 taulev(kmaxm1)
	478	!c DATA taulev/'tau1','tau2','tau3','tau4','tau5','tau6','tau7'/
	479	DATA taulev/'tau0','tau1','tau2','tau3','tau4','tau5','tau6'/
	480	CHARACTER *3 pclev(lmaxm1)
	481	DATA pclev/'pc1','pc2','pc3','pc4','pc5','pc6','pc7'/
	482	SAVE taulev,pclev
	483	!$OMP THREADPRIVATE(taulev,pclev)
	484	!c
	485	!c cnameisccp
	486	CHARACTER *29 cnameisccp(lmaxm1,kmaxm1)
	487	!IM bad 151205 DATA cnameisccp/'pc< 50hPa, tau< 0.3',
	488	DATA cnameisccp/'pc= 50-180hPa, tau< 0.3', &
	489	& 'pc= 180-310hPa, tau< 0.3', &
	490	& 'pc= 310-440hPa, tau< 0.3', &
	491	& 'pc= 440-560hPa, tau< 0.3', &
	492	& 'pc= 560-680hPa, tau< 0.3', &
	493	& 'pc= 680-800hPa, tau< 0.3', &
	494	& 'pc= 800-1000hPa, tau< 0.3', &
	495	& 'pc= 50-180hPa, tau= 0.3-1.3', &
	496	& 'pc= 180-310hPa, tau= 0.3-1.3', &
	497	& 'pc= 310-440hPa, tau= 0.3-1.3', &
	498	& 'pc= 440-560hPa, tau= 0.3-1.3', &
	499	& 'pc= 560-680hPa, tau= 0.3-1.3', &
	500	& 'pc= 680-800hPa, tau= 0.3-1.3', &
	501	& 'pc= 800-1000hPa, tau= 0.3-1.3', &
	502	& 'pc= 50-180hPa, tau= 1.3-3.6', &
	503	& 'pc= 180-310hPa, tau= 1.3-3.6', &
	504	& 'pc= 310-440hPa, tau= 1.3-3.6', &
	505	& 'pc= 440-560hPa, tau= 1.3-3.6', &
	506	& 'pc= 560-680hPa, tau= 1.3-3.6', &
	507	& 'pc= 680-800hPa, tau= 1.3-3.6', &
	508	& 'pc= 800-1000hPa, tau= 1.3-3.6', &
	509	& 'pc= 50-180hPa, tau= 3.6-9.4', &
	510	& 'pc= 180-310hPa, tau= 3.6-9.4', &
	511	& 'pc= 310-440hPa, tau= 3.6-9.4', &
	512	& 'pc= 440-560hPa, tau= 3.6-9.4', &
	513	& 'pc= 560-680hPa, tau= 3.6-9.4', &
	514	& 'pc= 680-800hPa, tau= 3.6-9.4', &
	515	& 'pc= 800-1000hPa, tau= 3.6-9.4', &
	516	& 'pc= 50-180hPa, tau= 9.4-23', &
	517	& 'pc= 180-310hPa, tau= 9.4-23', &
	518	& 'pc= 310-440hPa, tau= 9.4-23', &
	519	& 'pc= 440-560hPa, tau= 9.4-23', &
	520	& 'pc= 560-680hPa, tau= 9.4-23', &
	521	& 'pc= 680-800hPa, tau= 9.4-23', &
	522	& 'pc= 800-1000hPa, tau= 9.4-23', &
	523	& 'pc= 50-180hPa, tau= 23-60', &
	524	& 'pc= 180-310hPa, tau= 23-60', &
	525	& 'pc= 310-440hPa, tau= 23-60', &
	526	& 'pc= 440-560hPa, tau= 23-60', &
	527	& 'pc= 560-680hPa, tau= 23-60', &
	528	& 'pc= 680-800hPa, tau= 23-60', &
	529	& 'pc= 800-1000hPa, tau= 23-60', &
	530	& 'pc= 50-180hPa, tau> 60.', &
	531	& 'pc= 180-310hPa, tau> 60.', &
	532	& 'pc= 310-440hPa, tau> 60.', &
	533	& 'pc= 440-560hPa, tau> 60.', &
	534	& 'pc= 560-680hPa, tau> 60.', &
	535	& 'pc= 680-800hPa, tau> 60.', &
	536	& 'pc= 800-1000hPa, tau> 60.'/
	537	SAVE cnameisccp
	538	!$OMP THREADPRIVATE(cnameisccp)
	539	!c
	540	!c REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
	541	!c INTEGER nhorix7
	542	!IM: region='3d' <==> sorties en global
	543	CHARACTER*3 region
	544	PARAMETER(region='3d')
	545	!c
	546	!IM ISCCP simulator v3.4
	547	!c
	548	logical ok_hf
	549	!c
	550	integer nid_hf, nid_hf3d
	551	save ok_hf, nid_hf, nid_hf3d
	552	!$OMP THREADPRIVATE(ok_hf, nid_hf, nid_hf3d)
	553	!c QUESTION : noms de variables ?
	554
	555	INTEGER longcles
	556	PARAMETER ( longcles = 20 )
	557	REAL clesphy0( longcles )
	558	!c
	559	!c Variables propres a la physique
	560	INTEGER itap
	561	SAVE itap ! compteur pour la physique
	562	!$OMP THREADPRIVATE(itap)
	563	!c
	564	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	565	! real slp(klon) ! sea level pressure
	566	!c
	567	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	568	! REAL fevap(klon,nbsrf)
	569	! REAL fluxlat(klon,nbsrf)
	570	!c
	571	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	572	! REAL qsol(klon)
	573	REAL,save :: solarlong0
	574	!$OMP THREADPRIVATE(solarlong0)
	575
	576	!c
	577	!c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
	578	!c
	579	!IM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
	580	REAL zulow(klon),zvlow(klon)
	581	!c
	582	INTEGER igwd,idx(klon),itest(klon)
	583	!c
	584	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	585	! REAL agesno(klon,nbsrf)
	586
	587	!c
	588	!c REAL,allocatable,save :: run_off_lic_0(:)
	589	!cc$OMP THREADPRIVATE(run_off_lic_0)
	590	!cym SAVE run_off_lic_0
	591	!cKE43
	592	!c Variables liees a la convection de K. Emanuel (sb):
	593	!c
	594	REAL bas, top ! cloud base and top levels
	595	SAVE bas
	596	SAVE top
	597	!$OMP THREADPRIVATE(bas, top)
	598
	599	REAL wdn(klon), tdn(klon), qdn(klon)
	600	!c
	601	!c=================================================================================================
	602	!cCR04.12.07: on ajoute les nouvelles variables du nouveau schema de convection avec poches froides
	603	!c Variables li\'ees \`a la poche froide (jyg)
	604
	605	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
	606	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	607	! REAL Vprecip(klon,klev+1) ! precipitation vertical profile
	608	!c
	609	REAL wape_prescr, fip_prescr
	610	INTEGER it_wape_prescr
	611	SAVE wape_prescr, fip_prescr, it_wape_prescr
	612	!$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
	613	!c
	614	!c variables supplementaires de concvl
	615	REAL Tconv(klon,klev)
	616	REAL ment(klon,klev,klev),sij(klon,klev,klev)
	617	REAL dd_t(klon,klev),dd_q(klon,klev)
	618
	619	real, save :: alp_bl_prescr=0.
	620	real, save :: ale_bl_prescr=0.
	621
	622	real, save :: ale_max=1000.
	623	real, save :: alp_max=2.
	624
	625	real, save :: wake_s_min_lsp=0.1
	626
	627	!$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
	628	!$OMP THREADPRIVATE(ale_max,alp_max)
	629	!$OMP THREADPRIVATE(wake_s_min_lsp)
	630
	631	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	632	! real ale_wake(klon)
	633	! real alp_wake(klon)
	634
	635	real ok_wk_lsp(klon)
	636
	637	!cRC
	638	!c Variables li\'ees \`a la poche froide (jyg et rr)
	639	!c Version diagnostique pour l'instant : pas de r\'etroaction sur la convection
	640
	641	REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
	642
	643	REAL wake_dth(klon,klev) ! wake : temp pot difference
	644
	645	REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to Gravity Wave (/s)
	646	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta transported by LS omega
	647	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of large scale omega
	648	REAL wake_dtKE(klon,klev) ! Wake : differential heating (wake - unpertubed) CONV
	649	REAL wake_dqKE(klon,klev) ! Wake : differential moistening (wake - unpertubed) CONV
	650	REAL wake_dtPBL(klon,klev) ! Wake : differential heating (wake - unpertubed) PBL
	651	REAL wake_dqPBL(klon,klev) ! Wake : differential moistening (wake - unpertubed) PBL
	652	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	653	! REAL wake_omg(klon,klev) ! Wake : velocity difference (wake - unpertubed)
	654	REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
	655	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
	656	REAL wake_spread(klon,klev) ! spreading term in wake_delt
	657	!c
	658	!cpourquoi y'a pas de save??
	659	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	660	! REAL wake_h(klon) ! Wake : hauteur de la poche froide
	661	!c
	662	INTEGER wake_k(klon) ! Wake sommet
	663	!c
	664	REAL t_undi(klon,klev) ! temperature moyenne dans la zone non perturbee
	665	REAL q_undi(klon,klev) ! humidite moyenne dans la zone non perturbee
	666	!c
	667	!cjyg
	668	!ccc REAL wake_pe(klon) ! Wake potential energy - WAPE
	669
	670	REAL wake_gfl(klon) ! Gust Front Length
	671	REAL wake_dens(klon)
	672	!c
	673	!c
	674	REAL dt_dwn(klon,klev)
	675	REAL dq_dwn(klon,klev)
	676	REAL wdt_PBL(klon,klev)
	677	REAL udt_PBL(klon,klev)
	678	REAL wdq_PBL(klon,klev)
	679	REAL udq_PBL(klon,klev)
	680	REAL M_dwn(klon,klev)
	681	REAL M_up(klon,klev)
	682	REAL dt_a(klon,klev)
	683	REAL dq_a(klon,klev)
	684	REAL, SAVE :: alp_offset
	685	!$OMP THREADPRIVATE(alp_offset)
	686
	687	!c
	688	!cRR:fin declarations poches froides
	689	!c=======================================================================================================
	690
	691	REAL zw2(klon,klev+1)
	692
	693	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	694	! REAL fraca(klon,klev+1)
	695	REAL ztv(klon,klev),ztva(klon,klev)
	696	REAL zpspsk(klon,klev)
	697	REAL ztla(klon,klev),zqla(klon,klev)
	698	REAL zthl(klon,klev)
	699
	700	!ccc nrlmd le 10/04/2012
	701
	702	!c--------Stochastic Boundary Layer Triggering: ALE_BL--------
	703	!c---Propri\'et\'es du thermiques au LCL
	704	real zlcl_th(klon) ! Altitude du LCL calcul\'e continument (pcon dans thermcell_main.F90)
	705	real fraca0(klon) ! Fraction des thermiques au LCL
	706	real w0(klon) ! Vitesse des thermiques au LCL
	707	real w_conv(klon) ! Vitesse verticale de grande \'echelle au LCL
	708	real tke0(klon,klev+1) ! TKE au dÃ©but du pas de temps
	709	real therm_tke_max0(klon) ! TKE dans les thermiques au LCL
	710	real env_tke_max0(klon) ! TKE dans l'environnement au LCL
	711
	712	!c---Spectre de thermiques de type 2 au LCL
	713	real n2(klon),s2(klon)
	714	real ale_bl_stat(klon)
	715
	716	!c---D\'eclenchement stochastique
	717	integer :: tau_trig(klon)
	718	real proba_notrig(klon)
	719	real random_notrig(klon)
	720
	721	!c--------Statistical Boundary Layer Closure: ALP_BL--------
	722	!c---Profils de TKE dans et hors du thermique
	723	real pbl_tke_input(klon,klev+1,nbsrf)
	724	real therm_tke_max(klon,klev) ! Profil de TKE dans les thermiques
	725	real env_tke_max(klon,klev) ! Profil de TKE dans l'environnement
	726
	727	!c---Fermeture statistique
	728	real alp_bl_det(klon) ! ALP d\'terministe du thermique unique
	729	real alp_bl_fluct_m(klon) ! ALP li\'ee aux fluctuations de flux de masse sous-nuageux
	730	real alp_bl_fluct_tke(klon) ! ALP li\'ee aux fluctuations d'\'energie cin\'etique sous-nuageuse
	731	real alp_bl_conv(klon) ! ALP li\'ee \`a grande \'echelle
	732	real alp_bl_stat(klon) ! ALP totale
	733
	734	!ccc fin nrlmd le 10/04/2012
	735
	736	!c Variables locales pour la couche limite (al1):
	737	!c
	738	!cAl1 REAL pblh(klon) ! Hauteur de couche limite
	739	!cAl1 SAVE pblh
	740	!c34EK
	741	!c
	742	!c Variables locales:
	743	!c
	744	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	745	! REAL cdragh(klon) ! drag coefficient pour T and Q
	746	! REAL cdragm(klon) ! drag coefficient pour vent
	747	!cAA
	748	!cAA Pour phytrac
	749	REAL u1(klon) ! vents dans la premiere couche U
	750	REAL v1(klon) ! vents dans la premiere couche V
	751
	752	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule\n
	753	! REAL zxffonte(klon), zxfqcalving(klon),zxfqfonte(klon)
	754
	755	!@$$ LOGICAL offline ! Controle du stockage ds "physique"
	756	!@$$ PARAMETER (offline=.false.)
	757	!@$$ INTEGER physid
	758	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
	759	REAL frac_nucl(klon,klev) ! idem (nucleation)
	760	! RomP >>>
	761	REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
	762	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	763	! REAL beta_prec(klon,klev) ! taux de conv de l'eau cond (utilise)
	764	! RomP <<<
	765	INTEGER :: iii
	766	REAL :: calday
	767
	768	!IM cf FH pour Tiedtke 080604
	769	REAL rain_tiedtke(klon),snow_tiedtke(klon)
	770	!c
	771	!IM 050204 END
	772	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	773	! REAL evap(klon), devap(klon) ! evaporation et sa derivee
	774	REAL devap(klon) ! evaporation et sa derivee
	775
	776	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	777	! REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
	778	REAL dsens(klon) ! chaleur sensible et sa derivee
	779
	780	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	781	! REAL bils(klon) ! bilan de chaleur au sol
	782	! REAL wfbilo(klon,nbsrf) ! bilan d'eau, pour chaque
	783	!C ! type de sous-surface et pondere par la fraction
	784	! REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
	785	!C ! type de sous-surface et pondere par la fraction
	786	REAL slab_wfbils(klon) ! bilan de chaleur au sol pour le cas de slab, sur les points d'ocean
	787
	788	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	789	! REAL fder(klon)
	790	! REAL ve(klon) ! integr. verticale du transport meri. de l'energie
	791	! REAL vq(klon) ! integr. verticale du transport meri. de l'eau
	792	! REAL ue(klon) ! integr. verticale du transport zonal de l'energie
	793	! REAL uq(klon) ! integr. verticale du transport zonal de l'eau
	794	!c
	795
	796	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	797	! REAL frugs(klon,nbsrf)
	798	! REAL zxrugs(klon) ! longueur de rugosite
	799	!c
	800	!c Conditions aux limites
	801	!c
	802	!
	803	REAL :: day_since_equinox
	804	! Date de l'equinoxe de printemps
	805	INTEGER, parameter :: mth_eq=3, day_eq=21
	806	REAL :: jD_eq
	807
	808	LOGICAL, parameter :: new_orbit = .true.
	809
	810	!c
	811	INTEGER lmt_pas
	812	SAVE lmt_pas ! frequence de mise a jour
	813	!$OMP THREADPRIVATE(lmt_pas)
	814	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	815	! real zmasse(klon, llm),exner(klon, llm)
	816	real exner(klon, llm)
	817	!C (column-density of mass of air in a cell, in kg m-2)
	818	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
	819
	820	!IM sorties
	821	REAL un_jour
	822	PARAMETER(un_jour=86400.)
	823	!c======================================================================
	824	!c
	825	!c Declaration des procedures appelees
	826	!c
	827	EXTERNAL angle ! calculer angle zenithal du soleil
	828	EXTERNAL alboc ! calculer l'albedo sur ocean
	829	EXTERNAL ajsec ! ajustement sec
	830	EXTERNAL conlmd ! convection (schema LMD)
	831	!cKE43
	832	EXTERNAL conema3 ! convect4.3
	833	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
	834	!cAA
	835	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
	836	!c ! stockage des coefficients necessaires au
	837	!c ! lessivage OFF-LINE et ON-LINE
	838	EXTERNAL hgardfou ! verifier les temperatures
	839	EXTERNAL nuage ! calculer les proprietes radiatives
	840	!CC EXTERNAL o3cm ! initialiser l'ozone
	841	EXTERNAL orbite ! calculer l'orbite terrestre
	842	EXTERNAL phyetat0 ! lire l'etat initial de la physique
	843	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
	844	EXTERNAL suphel ! initialiser certaines constantes
	845	EXTERNAL transp ! transport total de l'eau et de l'energie
	846	EXTERNAL ecribina ! ecrire le fichier binaire global
	847	EXTERNAL ecribins ! ecrire le fichier binaire global
	848	EXTERNAL ecrirega ! ecrire le fichier binaire regional
	849	EXTERNAL ecriregs ! ecrire le fichier binaire regional
	850	!IM
	851	EXTERNAL haut2bas !variables de haut en bas
	852	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
	853	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
	854	!c EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
	855	!c EXTERNAL moyglo_aire !moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
	856	!c !par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
	857	!c
	858	!c Variables locales
	859	!c
	860	REAL rhcl(klon,klev) ! humiditi relative ciel clair
	861	REAL dialiq(klon,klev) ! eau liquide nuageuse
	862	REAL diafra(klon,klev) ! fraction nuageuse
	863	REAL cldliq(klon,klev) ! eau liquide nuageuse
	864	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	865	! REAL cldfra(klon,klev) ! fraction nuageuse
	866	! REAL cldtau(klon,klev) ! epaisseur optique
	867	! REAL cldemi(klon,klev) ! emissivite infrarouge
	868	!c
	869	!CXXX PB
	870	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	871	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
	872	! REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
	873	! REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
	874	! REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
	875	!c
	876	REAL zxfluxt(klon, klev)
	877	REAL zxfluxq(klon, klev)
	878	REAL zxfluxu(klon, klev)
	879	REAL zxfluxv(klon, klev)
	880	!CXXX
	881	!c
	882	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	883	! REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
	884	! REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
	885	!c Le rayonnement n'est pas calcule tous les pas, il faut donc
	886	!c sauvegarder les sorties du rayonnement
	887	!cym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
	888	!cym SAVE sollwdownclr, toplwdown, toplwdownclr
	889	!cym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
	890	!c
	891	INTEGER itaprad
	892	SAVE itaprad
	893	!$OMP THREADPRIVATE(itaprad)
	894	!c
	895	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
	896	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
	897	!c
	898	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	899	! REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
	900	! REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
	901	!c
	902	! REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
	903	REAL zxsnow_dummy(klon)
	904	!c
	905	REAL dist, rmu0(klon), fract(klon)
	906	REAL zdtime, zlongi
	907	!c
	908	CHARACTER*2 str2
	909	CHARACTER*2 iqn
	910	!c
	911	REAL qcheck
	912	REAL z_avant(klon), z_apres(klon), z_factor(klon)
	913	LOGICAL zx_ajustq
	914	!c
	915	REAL za, zb
	916	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
	917	real zqsat(klon,klev)
	918	INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq, iff
	919	REAL t_coup
	920	PARAMETER (t_coup=234.0)
	921	!c
	922	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	923	! REAL zphi(klon,klev)
	924	!cym A voir plus tard !!
	925	!cym REAL zx_relief(iim,jjmp1)
	926	!cym REAL zx_aire(iim,jjmp1)
	927	!c
	928	!c Grandeurs de sorties
	929	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	930	! REAL s_pblh(klon), s_lcl(klon), s_capCL(klon)
	931	! REAL s_oliqCL(klon), s_cteiCL(klon), s_pblt(klon)
	932	! REAL s_therm(klon), s_trmb1(klon), s_trmb2(klon)
	933	REAL s_capCL(klon)
	934	REAL s_oliqCL(klon), s_cteiCL(klon)
	935	REAL s_trmb1(klon), s_trmb2(klon)
	936	REAL s_trmb3(klon)
	937	!cKE43
	938	!c Variables locales pour la convection de K. Emanuel (sb):
	939	!c
	940	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	941	! REAL upwd(klon,klev) ! saturated updraft mass flux
	942	! REAL dnwd(klon,klev) ! saturated downdraft mass flux
	943	! REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
	944	REAL tvp(klon,klev) ! virtual temp of lifted parcel
	945	! REAL plcl(klon) ! Lifting Condensation Level
	946	! REAL plfc(klon) ! Level of Free Convection
	947	! REAL wbeff(klon) ! saturated updraft velocity at LFC
	948	CHARACTER*40 capemaxcels !max(CAPE)
	949
	950	REAL rflag(klon) ! flag fonctionnement de convect
	951	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
	952
	953	!c -- convect43:
	954	INTEGER ntra ! nb traceurs pour convect4.3
	955	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
	956	REAL dtma_con(klon),dtlcl_con(klon)
	957	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
	958	REAL dplcldt(klon), dplcldr(klon)
	959	!c? . condm_con(klon,klev),conda_con(klon,klev),
	960	!c? . mr_con(klon,klev),ep_con(klon,klev)
	961	!c? . ,sadiab(klon,klev),wadiab(klon,klev)
	962	!c --
	963	!c34EK
	964	!c
	965	!c Variables du changement
	966	!c
	967	!c con: convection
	968	!c lsc: condensation a grande echelle (Large-Scale-Condensation)
	969	!c ajs: ajustement sec
	970	!c eva: evaporation de l'eau liquide nuageuse
	971	!c vdf: couche limite (Vertical DiFfusion)
	972	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	973	! REAL rneb(klon,klev)
	974
	975	! tendance nulles
	976	REAL du0(klon,klev),dv0(klon,klev),dq0(klon,klev),dql0(klon,klev)
	977
	978	!c
	979	!*********************************************************
	980	!* declarations
	981
	982	!*********************************************************
	983	!IM 081204 END
	984	!c
	985	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	986	! REAL pmfu(klon,klev), pmfd(klon,klev)
	987	REAL pen_u(klon,klev), pen_d(klon,klev)
	988	REAL pde_u(klon,klev), pde_d(klon,klev)
	989	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
	990	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	991	! REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
	992	! REAL prfl(klon,klev+1), psfl(klon,klev+1)
	993	!c
	994	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	995	! REAL rain_lsc(klon)
	996	! REAL snow_lsc(klon)
	997	!c
	998	REAL ratqsc(klon,klev)
	999	real ratqsbas,ratqshaut,tau_ratqs
	1000	save ratqsbas,ratqshaut,tau_ratqs
	1001	!$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
	1002	real zpt_conv(klon,klev)
	1003
	1004	!c Parametres lies au nouveau schema de nuages (SB, PDF)
	1005	real fact_cldcon
	1006	real facttemps
	1007	logical ok_newmicro
	1008	save ok_newmicro
	1009	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	1010	! real ref_liq(klon,klev), ref_ice(klon,klev)
	1011
	1012	!$OMP THREADPRIVATE(ok_newmicro)
	1013	save fact_cldcon,facttemps
	1014	!$OMP THREADPRIVATE(fact_cldcon,facttemps)
	1015
	1016	integer iflag_cldcon
	1017	save iflag_cldcon
	1018	!$OMP THREADPRIVATE(iflag_cldcon)
	1019	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	1020	! logical ptconv(klon,klev)
	1021	!IM cf. AM 081204 BEG
	1022	logical ptconvth(klon,klev)
	1023	!IM cf. AM 081204 END
	1024	!c
	1025	!c Variables liees a l'ecriture de la bande histoire physique
	1026	!c
	1027	!c======================================================================
	1028	!c
	1029	!IM cf. AM 081204 BEG
	1030	!c declarations pour sortir sur une sous-region
	1031	integer imin_ins,imax_ins,jmin_ins,jmax_ins
	1032	save imin_ins,imax_ins,jmin_ins,jmax_ins
	1033	!$OMP THREADPRIVATE(imin_ins,imax_ins,jmin_ins,jmax_ins)
	1034	!c real lonmin_ins,lonmax_ins,latmin_ins
	1035	!c s ,latmax_ins
	1036	!c data lonmin_ins,lonmax_ins,latmin_ins
	1037	!c s ,latmax_ins/
	1038	!c valeurs initiales s -5.,20.,41.,55./
	1039	!c s 100.,130.,-20.,20./
	1040	!c s -180.,180.,-90.,90./
	1041	!c======================================================================
	1042	!IM cf. AM 081204 END
	1043
	1044	!c
	1045	integer itau_w ! pas de temps ecriture = itap + itau_phy
	1046	!c
	1047	!c
	1048	!c Variables locales pour effectuer les appels en serie
	1049	!c
	1050	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule\n
	1051	! REAL zx_rh(klon,klev)
	1052	!IM RH a 2m (la surface)
	1053	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	1054	! REAL rh2m(klon), qsat2m(klon)
	1055	! REAL tpot(klon), tpote(klon)
	1056	REAL Lheat
	1057
	1058	INTEGER length
	1059	PARAMETER ( length = 100 )
	1060	REAL tabcntr0( length )
	1061	!c
	1062	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
	1063	!IM
	1064	INTEGER ndex2d1(iwmax)
	1065	!c
	1066	!IM AMIP2 BEG
	1067	REAL moyglo, mountor
	1068	!IM 141004 BEG
	1069	REAL zustrdr(klon), zvstrdr(klon)
	1070	REAL zustrli(klon), zvstrli(klon)
	1071	REAL zustrph(klon), zvstrph(klon)
	1072	REAL zustrhi(klon), zvstrhi(klon)
	1073	REAL aam, torsfc
	1074	!IM 141004 END
	1075	!IM 190504 BEG
	1076	INTEGER ij, imp1jmp1
	1077	PARAMETER(imp1jmp1=(iim+1)*jjmp1)
	1078	!cym A voir plus tard
	1079	REAL zx_tmp(imp1jmp1), airedyn(iim+1,jjmp1)
	1080	REAL padyn(iim+1,jjmp1,klev+1)
	1081	REAL dudyn(iim+1,jjmp1,klev)
	1082	REAL rlatdyn(iim+1,jjmp1)
	1083	!IM 190504 END
	1084	LOGICAL ok_msk
	1085	REAL msk(klon)
	1086	!IM
	1087	REAL airetot, pi
	1088	!cym A voir plus tard
	1089	!cym REAL zm_wo(jjmp1, klev)
	1090	!IM AMIP2 END
	1091	!c
	1092	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
	1093	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
	1094	REAL zx_tmp_fi3d1(klon,klev+1) !variable temporaire pour champs 3D (kelvp1)
	1095	REAL(KIND=8) zx_tmp2_fi3d(klon,klev) ! variable temporaire pour champs 3D
	1096	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
	1097	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
	1098	!c
	1099	INTEGER nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
	1100	INTEGER nid_hfnmc, nid_day_seri, nid_ctesGCM
	1101	SAVE nid_day, nid_mth, nid_ins, nid_mthnmc, nid_daynmc
	1102	SAVE nid_hfnmc, nid_day_seri, nid_ctesGCM
	1103	!$OMP THREADPRIVATE(nid_day, nid_mth, nid_ins)
	1104	!$OMP THREADPRIVATE(nid_mthnmc, nid_daynmc, nid_hfnmc)
	1105	!$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
	1106	!c
	1107	!IM 280405 BEG
	1108	INTEGER nid_bilKPins, nid_bilKPave
	1109	SAVE nid_bilKPins, nid_bilKPave
	1110	!$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
	1111	!c
	1112	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
	1113	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
	1114	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
	1115	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
	1116	!c
	1117	INTEGER nhori, nvert, nvert1, nvert3
	1118	REAL zsto, zsto1, zsto2
	1119	REAL zstophy, zstorad, zstohf, zstoday, zstomth, zout
	1120	REAL zcals(napisccp), zcalh(napisccp), zoutj(napisccp)
	1121	REAL zout_isccp(napisccp)
	1122	SAVE zcals, zcalh, zoutj, zout_isccp
	1123	!$OMP THREADPRIVATE(zcals, zcalh, zoutj, zout_isccp)
	1124
	1125	real zjulian
	1126	save zjulian
	1127	!$OMP THREADPRIVATE(zjulian)
	1128
	1129	character*20 modname
	1130	character*80 abort_message
	1131	logical ok_sync
	1132	real date0
	1133	integer idayref
	1134
	1135	!C essai writephys
	1136	integer fid_day, fid_mth, fid_ins
	1137	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
	1138	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
	1139	parameter (prof2d_on = 1, prof3d_on = 2, &
	1140	& prof2d_av = 3, prof3d_av = 4)
	1141	character*30 nom_fichier
	1142	character*40 varname
	1143	character*40 vartitle
	1144	character*20 varunits
	1145	!C Variables liees au bilan d'energie et d'enthalpi
	1146	REAL ztsol(klon)
	1147	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot &
	1148	& , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	1149	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot &
	1150	& , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
	1151	!$OMP THREADPRIVATE(h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot,
	1152	!$OMP+ h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot)
	1153	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
	1154	REAL d_h_vcol_phy
	1155	REAL fs_bound, fq_bound
	1156	SAVE d_h_vcol_phy
	1157	!$OMP THREADPRIVATE(d_h_vcol_phy)
	1158	REAL zero_v(klon)
	1159	CHARACTER*15 ztit
	1160	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
	1161	SAVE ip_ebil
	1162	DATA ip_ebil/0/
	1163	!$OMP THREADPRIVATE(ip_ebil)
	1164	INTEGER if_ebil ! level for energy conserv. dignostics
	1165	SAVE if_ebil
	1166	!$OMP THREADPRIVATE(if_ebil)
	1167	!c+jld ec_conser
	1168	REAL ZRCPD
	1169	!c-jld ec_conser
	1170	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	1171	! REAL t2m(klon,nbsrf) ! temperature a 2m
	1172	REAL q2m(klon,nbsrf) ! humidite a 2m
	1173
	1174	!IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
	1175	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	1176	! REAL zt2m(klon), zq2m(klon) !temp., hum. 2m moyenne s/ 1 maille
	1177	! REAL zustar(klon),zu10m(klon), zv10m(klon) ! u* et vents a 10m moyennes s/1 maille
	1178	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
	1179	CHARACTER*40 tinst, tave, typeval
	1180	REAL cldtaupi(klon,klev) ! Cloud optical thickness for pre-industrial (pi) aerosols
	1181
	1182	! L. Fita, LMD. January 2014. Defined in output_lmdz_NOmodule
	1183	! REAL re(klon, klev) ! Cloud droplet effective radius
	1184	! REAL fl(klon, klev) ! denominator of re
	1185
	1186	REAL re_top(klon), fl_top(klon) ! CDR at top of liquid water clouds
	1187
	1188	! Aerosol optical properties
	1189	CHARACTER*4, DIMENSION(naero_grp) :: rfname
	1190	REAL, DIMENSION(klon) :: aerindex ! POLDER aerosol index
	1191	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass concentration for all soluble aerosols[ug/m3]
	1192	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi ! - " - (pre-industrial value)
	1193	INTEGER :: naero ! aerosol species
	1194
	1195	! Parameters
	1196	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
	1197	LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
	1198	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
	1199	SAVE ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1
	1200	!$OMP THREADPRIVATE(ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1)
	1201	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
	1202	! false : lecture des aerosol dans un fichier
	1203	!$OMP THREADPRIVATE(aerosol_couple)
	1204	INTEGER, SAVE :: flag_aerosol
	1205	!$OMP THREADPRIVATE(flag_aerosol)
	1206	LOGICAL, SAVE :: new_aod
	1207	!$OMP THREADPRIVATE(new_aod)
	1208	!c
	1209	!c--STRAT AEROSOL
	1210	LOGICAL, SAVE :: flag_aerosol_strat
	1211	!$OMP THREADPRIVATE(flag_aerosol_strat)
	1212	!cc-fin STRAT AEROSOL
	1213	!c
	1214	!c Declaration des constantes et des fonctions thermodynamiques
	1215	!c
	1216	LOGICAL,SAVE :: first=.true.
	1217	!$OMP THREADPRIVATE(first)
	1218
	1219	integer iunit
	1220
	1221	integer, save:: read_climoz ! read ozone climatology
	1222	!C (let it keep the default OpenMP shared attribute)
	1223	!C Allowed values are 0, 1 and 2
	1224	!C 0: do not read an ozone climatology
	1225	!C 1: read a single ozone climatology that will be used day and night
	1226	!C 2: read two ozone climatologies, the average day and night
	1227	!C climatology and the daylight climatology
	1228
	1229	integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
	1230	!C (let it keep the default OpenMP shared attribute)
	1231
	1232	real, pointer, save:: press_climoz(:)
	1233	!C (let it keep the default OpenMP shared attribute)
	1234	! edges of pressure intervals for ozone climatologies, in Pa, in strictly
	1235	! ascending order
	1236
	1237	integer, save:: co3i = 0
	1238	! time index in NetCDF file of current ozone fields
	1239	!$OMP THREADPRIVATE(co3i)
	1240
	1241	integer ro3i
	1242	! required time index in NetCDF file for the ozone fields, between 1
	1243	! and 360
	1244
	1245	INTEGER ierr
	1246	#include "YOMCST.h"
	1247	#include "YOETHF.h"
	1248	#include "FCTTRE.h"
	1249	!IM 100106 BEG : pouvoir sortir les ctes de la physique
	1250	#include "conema3.h"
	1251	#include "fisrtilp.h"
	1252	#include "nuage.h"
	1253	#include "compbl.h"
	1254	!IM 100106 END : pouvoir sortir les ctes de la physique
	1255	!c
	1256	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1257	!c Declarations pour Simulateur COSP
	1258	!c============================================================
	1259	real :: mr_ozone(klon,klev)
	1260
	1261	!IM sorties fichier 1D paramLMDZ_phy.nc
	1262	REAL :: zx_tmp_0d(1,1)
	1263	INTEGER, PARAMETER :: np=1
	1264	REAL,dimension(klon_glo) :: rlat_glo
	1265	REAL,dimension(klon_glo) :: rlon_glo
	1266	REAL gbils(1), gevap(1), gevapt(1), glat(1), gnet0(1), gnet(1)
	1267	REAL grain(1), gtsol(1), gt2m(1), gprw(1)
	1268
	1269	!IM stations CFMIP
	1270	INTEGER, SAVE :: nCFMIP
	1271	!$OMP THREADPRIVATE(nCFMIP)
	1272	INTEGER, PARAMETER :: npCFMIP=120
	1273	INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
	1274	REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
	1275	!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
	1276	INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
	1277	REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
	1278	!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
	1279	INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
	1280	!$OMP THREADPRIVATE(iGCM, jGCM)
	1281	logical, dimension(nfiles) :: phys_out_filestations
	1282	logical, parameter :: lNMC=.FALSE.
	1283
	1284	!IM betaCRF
	1285	REAL, SAVE :: pfree, beta_pbl, beta_free
	1286	!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
	1287	REAL, SAVE :: lon1_beta, lon2_beta, lat1_beta, lat2_beta
	1288	!$OMP THREADPRIVATE(lon1_beta, lon2_beta, lat1_beta, lat2_beta)
	1289	LOGICAL, SAVE :: mskocean_beta
	1290	!$OMP THREADPRIVATE(mskocean_beta)
	1291	REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et cldemirad pour evaluer les retros liees aux CRF
	1292	REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique pour radlwsw pour tester "CRF off"
	1293	REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique pour radlwsw pour tester "CRF off"
	1294	REAL, dimension(klon, klev) :: cldemirad ! emissivite pour radlwsw pour tester "CRF off"
	1295	REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse
	1296
	1297	INTEGER :: nbtr_tmp ! Number of tracer inside concvl
	1298	REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
	1299	integer iostat
	1300
[142]	1301	! L. Fita, LMD. Point for checkings...
[1]	1302	INTEGER :: llp
[183]	1303	llp = 449
[1]	1304
	1305	!c======================================================================
	1306	! Gestion calendrier : mise a jour du module phys_cal_mod
	1307	!
	1308	CALL phys_cal_update(jD_cur,jH_cur)
	1309
	1310	!c======================================================================
	1311	! Ecriture eventuelle d'un profil verticale en entree de la physique.
	1312	! Utilise notamment en 1D mais peut etre active egalement en 3D
	1313	! en imposant la valeur de igout.
	1314	!c======================================================================
	1315
	1316	if (prt_level.ge.1) then
	1317	! Lluis
	1318	! igout=klon/2+1/klon
	1319	igout=llp
	1320	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	1321	write(lunout,*) &
	1322	& 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
	1323	write(lunout,*) &
	1324	& nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
	1325
	1326	write(lunout,*) 'paprs, play, phi, u, v, t'
	1327	do k=1,klev
	1328	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
	1329	& u(igout,k),v(igout,k),t(igout,k)
	1330	enddo
	1331	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
	1332	do k=1,klev
	1333	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
	1334	enddo
	1335	endif
	1336
	1337	!c======================================================================
	1338
	1339	!cym => necessaire pour iflag_con != 2
	1340	pmfd(:,:) = 0.
	1341	pen_u(:,:) = 0.
	1342	pen_d(:,:) = 0.
	1343	pde_d(:,:) = 0.
	1344	pde_u(:,:) = 0.
	1345	aam=0.
	1346
	1347	torsfc=0.
	1348
	1349	forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
	1350
	1351	if (first) then
	1352
	1353	!cCR:nvelles variables convection/poches froides
	1354
	1355	print*, '================================================='
	1356	print*, 'Allocation des variables locales et sauvegardees'
	1357
	1358	call phys_local_var_init
	1359
	1360	!c
	1361	pasphys=pdtphys
	1362	!c appel a la lecture du run.def physique
	1363	call conf_phys(ok_journe, ok_mensuel, &
	1364	& ok_instan, ok_hf, &
	1365	& ok_LES, &
	1366	& callstats, &
	1367	& solarlong0,seuil_inversion, &
	1368	& fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
	1369	& iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
	1370	& ok_ade, ok_aie, ok_cdnc, aerosol_couple, &
	1371	& flag_aerosol, flag_aerosol_strat, new_aod, &
	1372	& bl95_b0, bl95_b1, &
	1373	!c nv flags pour la convection et les poches froides &
	1374	& read_climoz, &
	1375	& alp_offset)
	1376
	1377
	1378	!! call phys_state_var_init(read_climoz)
	1379	call phys_output_var_init
	1380
	1381
	1382	!L. Fita, LMD. November 2013
	1383	! Initializing
	1384	! Variables have been already initialilzed from the WRF variables.
	1385	! Using on purpose subroutines from physiq_limit_variables_mod
	1386	!! CALL neige_initialize()
	1387	!! CALL limit_initialize()
	1388	!! CALL vars_limit_init(klon)
[6]	1389
[1]	1390	! L. Fita, LMD January 2014. Initializing output variables
	1391	CALL init_ovars_lmdz_NOmodule(klon,klev,nbsrf)
	1392
	1393	! Allocating restart variables
	1394	! ALLOCATE(qsol(klon))
	1395	! ALLOCATE(fder(klon))
	1396	! ALLOCATE(snow(klon, nbsrf))
	1397	! ALLOCATE(qsurf(klon, nbsrf))
	1398	! ALLOCATE(evap(klon, nbsrf))
	1399	! ALLOCATE(rugos(klon, nbsrf))
	1400	! ALLOCATE(agesno(klon, nbsrf))
	1401	! ALLOCATE(ftsoil(klon, nsoilmx, nbsrf))
	1402	!! ALLOCATE(restart_runoff(klon))
	1403
	1404	! CALL pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst, evap_rst, &
	1405	! rugos_rst, agesno_rst, ftsoil_rst)
	1406	! CALL fonte_neige_init(restart_runoff)
	1407	! qsol=qsol_rst
	1408
	1409	print*, '================================================='
	1410	!c
	1411	dnwd0=0.0
	1412	ftd=0.0
	1413	fqd=0.0
	1414	cin=0.
	1415	!cym Attention pbase pas initialise dans concvl !!!!
	1416	pbase=0
	1417	!IM 180608
	1418
	1419	itau_con=0
	1420	first=.false.
	1421
	1422	endif ! first
	1423
	1424	modname = 'physiq'
	1425	!IM
	1426	IF (ip_ebil_phy.ge.1) THEN
	1427	DO i=1,klon
	1428	zero_v(i)=0.
	1429	END DO
	1430	END IF
	1431	ok_sync=.TRUE.
	1432
	1433	IF (debut) THEN
	1434	CALL suphel ! initialiser constantes et parametres phys.
	1435	ENDIF
	1436
	1437	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
	1438
	1439
	1440	!c======================================================================
	1441	! Gestion calendrier : mise a jour du module phys_cal_mod
	1442	!
	1443	!c CALL phys_cal_update(jD_cur,jH_cur)
	1444
	1445	!c
	1446	!c Si c'est le debut, il faut initialiser plusieurs choses
	1447	!c ********
	1448	!c
	1449
	1450	IF (debut) THEN
	1451	!rv
	1452	!cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
	1453	!cde la convection a partir des caracteristiques du thermique
	1454	wght_th(:,:)=1.
	1455	lalim_conv(:)=1
	1456	!cRC
	1457	ustar(:,:)=0.
	1458	u10m(:,:)=0.
	1459	v10m(:,:)=0.
	1460	rain_con(:)=0.
	1461	snow_con(:)=0.
	1462	topswai(:)=0.
	1463	topswad(:)=0.
	1464	solswai(:)=0.
	1465	solswad(:)=0.
	1466
	1467	wmax_th(:)=0.
	1468	tau_overturning_th(:)=0.
	1469
	1470	IF (type_trac == 'inca') THEN
	1471	! jg : initialisation jusqu'au ces variables sont dans restart
	1472	ccm(:,:,:) = 0.
	1473	tau_aero(:,:,:,:) = 0.
	1474	piz_aero(:,:,:,:) = 0.
	1475	cg_aero(:,:,:,:) = 0.
	1476	END IF
	1477
	1478	rnebcon0(:,:) = 0.0
	1479	clwcon0(:,:) = 0.0
	1480	rnebcon(:,:) = 0.0
	1481	clwcon(:,:) = 0.0
	1482
	1483	!IM
	1484	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
	1485	!c
	1486	print*,'iflag_coupl,iflag_clos,iflag_wake', &
	1487	& iflag_coupl,iflag_clos,iflag_wake
	1488	print*,'CYCLE_DIURNE', cycle_diurne
	1489	!c
	1490	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
	1491	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
	1492	CALL abort_gcm (modname,abort_message,1)
	1493	ENDIF
	1494	!c
	1495	IF(ok_isccp.AND.iflag_con.LE.2) THEN
	1496	abort_message = 'ISCCP-like outputs may be available for KE &
	1497	&(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
	1498	CALL abort_gcm (modname,abort_message,1)
	1499	ENDIF
	1500	!c
	1501	!c Initialiser les compteurs:
	1502	!c
	1503	itap = 0
	1504	itaprad = 0
	1505
	1506	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1507	!! Un petit travail \`a faire ici.
	1508	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1509
	1510	if (iflag_pbl>1) then
	1511	PRINT*, "Using method MELLOR&YAMADA"
	1512	endif
	1513
	1514	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1515	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
	1516	! dyn3d
	1517	! Attention : la version precedente n'etait pas tres propre.
	1518	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
	1519	! pour obtenir le meme resultat.
	1520	dtime=pdtphys
	1521	radpas = NINT( 86400./dtime/nbapp_rad)
	1522	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1523
	1524	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
	1525
	1526	IF (klon_glo==1) THEN
	1527	coefh=0. ; coefm=0. ; pbl_tke=0.
	1528	coefh(:,2,:)=1.e-2 ; coefm(:,2,:)=1.e-2 ; pbl_tke(:,2,:)=1.e-2
	1529	PRINT*,'FH WARNING : lignes a supprimer'
	1530	ENDIF
	1531	!IM begin
	1532	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
	1533	&,ratqs(1,1)
	1534	!IM end
	1535
	1536
	1537
	1538	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1539	!c
	1540	!C on remet le calendrier a zero
	1541	!c
	1542	IF (raz_date .eq. 1) THEN
	1543	itau_phy = 0
	1544	ENDIF
	1545
	1546	!IM cf. AM 081204 BEG
	1547	PRINT*,'cycle_diurne3 =',cycle_diurne
	1548	!IM cf. AM 081204 END
	1549	!c
	1550	CALL printflag( tabcntr0,radpas,ok_journe, &
	1551	& ok_instan, ok_region )
	1552	!c
	1553	IF (ABS(dtime-pdtphys).GT.0.001) THEN
	1554	WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
	1555	& pdtphys
	1556	abort_message='Pas physique n est pas correct '
	1557	! call abort_gcm(modname,abort_message,1)
	1558	dtime=pdtphys
	1559	ENDIF
	1560	IF (nlon .NE. klon) THEN
	1561	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
	1562	& klon
	1563	abort_message='nlon et klon ne sont pas coherents'
	1564	call abort_gcm(modname,abort_message,1)
	1565	ENDIF
	1566	IF (nlev .NE. klev) THEN
	1567	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
	1568	& klev
	1569	abort_message='nlev et klev ne sont pas coherents'
	1570	call abort_gcm(modname,abort_message,1)
	1571	ENDIF
	1572	!c
	1573	IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
	1574	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	1575	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	1576	abort_message='Nbre d appels au rayonnement insuffisant'
	1577	call abort_gcm(modname,abort_message,1)
	1578	ENDIF
	1579	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
	1580	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
	1581	& ok_cvl
	1582	!c
	1583	!cKE43
	1584	!c Initialisation pour la convection de K.E. (sb):
	1585	IF (iflag_con.GE.3) THEN
	1586
	1587	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
	1588	WRITE(lunout,*) &
	1589	& "On va utiliser le melange convectif des traceurs qui"
	1590	WRITE(lunout,*)"est calcule dans convect4.3"
	1591	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
	1592
	1593	DO i = 1, klon
	1594	ema_cbmf(i) = 0.
	1595	ema_pcb(i) = 0.
	1596	ema_pct(i) = 0.
	1597	!c ema_workcbmf(i) = 0.
	1598	ENDDO
	1599	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
	1600	DO i = 1, klon
	1601	ibas_con(i) = 1
	1602	itop_con(i) = 1
	1603	ENDDO
	1604	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
	1605	!c===============================================================================
	1606	!cCR:04.12.07: initialisations poches froides
	1607	!c Controle de ALE et ALP pour la fermeture convective (jyg)
	1608	if (iflag_wake>=1) then
	1609	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
	1610	& ,alp_bl_prescr, ale_bl_prescr)
	1611	!c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
	1612	!c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
	1613	endif
	1614
	1615	do i = 1,klon
	1616	Ale_bl(i)=0.
	1617	Alp_bl(i)=0.
	1618	enddo
	1619
	1620	!c================================================================================
	1621	!IM stations CFMIP
	1622	nCFMIP=npCFMIP
	1623	OPEN(98,file='npCFMIP_param.data',status='old', &
	1624	& form='formatted',iostat=iostat)
	1625	if (iostat == 0) then
	1626	READ(98,*,end=998) nCFMIP
	1627	998 CONTINUE
	1628	CLOSE(98)
	1629	CONTINUE
	1630	IF(nCFMIP.GT.npCFMIP) THEN
	1631	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
	1632	CALL abort
	1633	else
	1634	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
	1635	ENDIF
	1636
	1637	!c
	1638	ALLOCATE(tabCFMIP(nCFMIP))
	1639	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
	1640	ALLOCATE(tabijGCM(nCFMIP))
	1641	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
	1642	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
	1643	!c
	1644	!c lecture des nCFMIP stations CFMIP, de leur numero
	1645	!c et des coordonnees geographiques lonCFMIP, latCFMIP
	1646	!c
	1647	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
	1648	&lonCFMIP, latCFMIP)
	1649	!c
	1650	!c identification des
	1651	!c 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
	1652	!c 2) indices points tabijGCM de la grille physique 1d sur klon points
	1653	!c 3) indices iGCM, jGCM de la grille physique 2d
	1654	!c
	1655	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
	1656	&tabijGCM, lonGCM, latGCM, iGCM, jGCM)
	1657	!c
	1658	else
	1659	ALLOCATE(tabijGCM(0))
	1660	ALLOCATE(lonGCM(0), latGCM(0))
	1661	ALLOCATE(iGCM(0), jGCM(0))
	1662	end if
	1663	else
	1664	ALLOCATE(tabijGCM(0))
	1665	ALLOCATE(lonGCM(0), latGCM(0))
	1666	ALLOCATE(iGCM(0), jGCM(0))
	1667	ENDIF
	1668
	1669	DO i=1,klon
	1670	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1671	ENDDO
	1672
	1673	!c34EK
	1674	IF (ok_orodr) THEN
	1675
	1676	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1677	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
	1678	! justement quand ok_orodr = false.
	1679	! ce rugoro est utilise par la couche limite et fait double emploi
	1680	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
	1681	! DO i=1,klon
	1682	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1683	! ENDDO
	1684	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1685	IF (ok_strato) THEN
	1686	CALL SUGWD_strato(klon,klev,paprs,pplay)
	1687	ELSE
	1688	CALL SUGWD(klon,klev,paprs,pplay)
	1689	ENDIF
	1690
	1691	DO i=1,klon
	1692	zuthe(i)=0.
	1693	zvthe(i)=0.
	1694	if(zstd(i).gt.10.)then
	1695	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	1696	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	1697	endif
	1698	ENDDO
	1699	ENDIF
	1700	!c
	1701	!c
	1702	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
	1703	WRITE(lunout,*)'La frequence de lecture surface est de ', &
	1704	& lmt_pas
	1705	!c
	1706	capemaxcels = 't_max(X)'
	1707	t2mincels = 't_min(X)'
	1708	t2maxcels = 't_max(X)'
	1709	tinst = 'inst(X)'
	1710	tave = 'ave(X)'
	1711	!IM cf. AM 081204 BEG
	1712	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
	1713	!IM cf. AM 081204 END
	1714	!c
	1715	!c=============================================================
	1716	!c Initialisation des sorties
	1717	!c=============================================================
	1718
	1719	#ifdef CPP_IOIPSL
	1720
	1721	!$OMP MASTER
	1722	call phys_output_open(rlon,rlat,nCFMIP,tabijGCM, &
	1723	& iGCM,jGCM,lonGCM,latGCM, &
	1724	& jjmp1,nlevSTD,clevSTD, &
	1725	& nbteta, ctetaSTD, dtime,ok_veget, &
	1726	& type_ocean,iflag_pbl,ok_mensuel,ok_journe, &
	1727	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
	1728	& read_climoz, phys_out_filestations, &
	1729	& new_aod, aerosol_couple, &
	1730	& flag_aerosol_strat )
	1731	!$OMP END MASTER
	1732	!$OMP BARRIER
	1733
	1734	#ifdef histISCCP
	1735	#include "ini_histISCCP.h"
	1736	#endif
	1737
	1738	#ifdef histNMC
	1739	#include "ini_histhfNMC.h"
	1740	#include "ini_histdayNMC.h"
	1741	#include "ini_histmthNMC.h"
	1742	#endif
	1743
	1744	#include "ini_histday_seri.h"
	1745
	1746	#include "ini_paramLMDZ_phy.h"
	1747
	1748	#endif
	1749	ecrit_reg = ecrit_reg * un_jour
	1750	ecrit_tra = ecrit_tra * un_jour
	1751
	1752	!cXXXPB Positionner date0 pour initialisation de ORCHIDEE
	1753	date0 = jD_ref
	1754	WRITE(,) 'physiq date0 : ',date0
	1755	!c
	1756	!c
	1757	!c
	1758	!c Prescrire l'ozone dans l'atmosphere
	1759	!c
	1760	!c
	1761	!cc DO i = 1, klon
	1762	!cc DO k = 1, klev
	1763	!cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
	1764	!cc ENDDO
	1765	!cc ENDDO
	1766	!c
	1767	IF (type_trac == 'inca') THEN
	1768	#ifdef INCA
	1769	CALL VTe(VTphysiq)
	1770	CALL VTb(VTinca)
	1771	! iii = MOD(NINT(xjour),360)
	1772	! calday = REAL(iii) + jH_cur
	1773	calday = REAL(days_elapsed) + jH_cur
	1774	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
	1775
	1776	CALL chemini( &
	1777	& rg, &
	1778	& ra, &
	1779	& airephy, &
	1780	& rlat, &
	1781	& rlon, &
	1782	& presnivs, &
	1783	& calday, &
	1784	& klon, &
	1785	& nqtot, &
	1786	& pdtphys, &
	1787	& annee_ref, &
	1788	& day_ref, &
	1789	& itau_phy)
	1790
	1791	CALL VTe(VTinca)
	1792	CALL VTb(VTphysiq)
	1793	#endif
	1794	END IF
	1795	!c
	1796	!c
	1797	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1798	! Nouvelle initialisation pour le rayonnement RRTM
	1799	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1800
	1801	call iniradia(klon,klev,paprs(1,1:klev+1))
	1802
	1803	!$OMP single
	1804	if (read_climoz >= 1) then
	1805	call open_climoz(ncid_climoz, press_climoz)
	1806	END IF
	1807	!$OMP end single
	1808	!c
	1809	!IM betaCRF
	1810	pfree=70000. !Pa
	1811	beta_pbl=1.
	1812	beta_free=1.
	1813	lon1_beta=-180.
	1814	lon2_beta=+180.
	1815	lat1_beta=90.
	1816	lat2_beta=-90.
	1817	mskocean_beta=.FALSE.
	1818
	1819	OPEN(99,file='beta_crf.data',status='old', &
	1820	& form='formatted',err=9999)
	1821	READ(99,*,end=9998) pfree
	1822	READ(99,*,end=9998) beta_pbl
	1823	READ(99,*,end=9998) beta_free
	1824	READ(99,*,end=9998) lon1_beta
	1825	READ(99,*,end=9998) lon2_beta
	1826	READ(99,*,end=9998) lat1_beta
	1827	READ(99,*,end=9998) lat2_beta
	1828	READ(99,*,end=9998) mskocean_beta
	1829	9998 Continue
	1830	CLOSE(99)
	1831	9999 Continue
	1832	WRITE(,)'pfree=',pfree
	1833	WRITE(,)'beta_pbl=',beta_pbl
	1834	WRITE(,)'beta_free=',beta_free
	1835	WRITE(,)'lon1_beta=',lon1_beta
	1836	WRITE(,)'lon2_beta=',lon2_beta
	1837	WRITE(,)'lat1_beta=',lat1_beta
	1838	WRITE(,)'lat2_beta=',lat2_beta
	1839	WRITE(,)'mskocean_beta=',mskocean_beta
	1840	ENDIF
	1841	!
	1842	! ************** Fin de IF ( debut ) *************
	1843	!
	1844	!
	1845	! Incrementer le compteur de la physique
	1846	!
	1847	itap = itap + 1
	1848	!c
	1849	!
	1850	! Update fraction of the sub-surfaces (pctsrf) and
	1851	! initialize, where a new fraction has appeared, all variables depending
	1852	! on the surface fraction.
	1853	!
	1854	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
	1855	& pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
	1856
	1857	! Update time and other variables in Reprobus
	1858	IF (type_trac == 'repr') THEN
	1859	#ifdef REPROBUS
	1860	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
	1861	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
	1862	CALL Rtime(debut)
	1863	#endif
	1864	END IF
	1865
	1866
	1867	! Tendances bidons pour les processus qui n'affectent pas certaines
	1868	! variables.
	1869	du0(:,:)=0.
	1870	dv0(:,:)=0.
	1871	dq0(:,:)=0.
	1872	dql0(:,:)=0.
	1873	!c
	1874	!c Mettre a zero des variables de sortie (pour securite)
	1875	!c
	1876	DO i = 1, klon
	1877	d_ps(i) = 0.0
	1878	ENDDO
	1879	DO k = 1, klev
	1880	DO i = 1, klon
	1881	d_t(i,k) = 0.0
	1882	d_u(i,k) = 0.0
	1883	d_v(i,k) = 0.0
	1884	ENDDO
	1885	ENDDO
	1886	DO iq = 1, nqtot
	1887	DO k = 1, klev
	1888	DO i = 1, klon
	1889	d_qx(i,k,iq) = 0.0
	1890	ENDDO
	1891	ENDDO
	1892	ENDDO
	1893	da(:,:)=0.
	1894	mp(:,:)=0.
	1895	phi(:,:,:)=0.
	1896	! RomP >>>
	1897	phi2(:,:,:)=0.
	1898	beta_prec_fisrt(:,:)=0.
	1899	beta_prec(:,:)=0.
	1900	epmlmMm(:,:,:)=0.
	1901	eplaMm(:,:)=0.
	1902	d1a(:,:)=0.
	1903	dam(:,:)=0.
	1904	pmflxr=0.
	1905	pmflxs=0.
	1906	! RomP <<<
	1907
	1908	!c
	1909	!c Ne pas affecter les valeurs entrees de u, v, h, et q
	1910	!c
	1911	DO k = 1, klev
	1912	DO i = 1, klon
	1913	t_seri(i,k) = t(i,k)
	1914	u_seri(i,k) = u(i,k)
	1915	v_seri(i,k) = v(i,k)
	1916	q_seri(i,k) = qx(i,k,ivap)
	1917	ql_seri(i,k) = qx(i,k,iliq)
	1918	qs_seri(i,k) = 0.
	1919	ENDDO
	1920	ENDDO
	1921	tke0(:,:)=pbl_tke(:,:,is_ave)
	1922	IF (nqtot.GE.3) THEN
	1923	DO iq = 3, nqtot
	1924	DO k = 1, klev
	1925	DO i = 1, klon
	1926	tr_seri(i,k,iq-2) = qx(i,k,iq)
	1927	ENDDO
	1928	ENDDO
	1929	ENDDO
	1930	ELSE
	1931	DO k = 1, klev
	1932	DO i = 1, klon
	1933	tr_seri(i,k,1) = 0.0
	1934	ENDDO
	1935	ENDDO
	1936	ENDIF
	1937	!C
	1938	DO i = 1, klon
	1939	ztsol(i) = 0.
	1940	ENDDO
	1941	DO nsrf = 1, nbsrf
	1942	DO i = 1, klon
	1943	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
	1944	ENDDO
	1945	ENDDO
	1946	!IM
	1947	IF (ip_ebil_phy.ge.1) THEN
	1948	ztit='after dynamic'
	1949	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
	1950	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	1951	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1952	!C Comme les tendances de la physique sont ajoute dans la dynamique,
	1953	!C on devrait avoir que la variation d'entalpie par la dynamique
	1954	!C est egale a la variation de la physique au pas de temps precedent.
	1955	!C Donc la somme de ces 2 variations devrait etre nulle.
	1956	call diagphy(airephy,ztit,ip_ebil_phy &
	1957	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	1958	& , zero_v, zero_v, zero_v, ztsol &
	1959	& , d_h_vcol+d_h_vcol_phy, d_qt, 0. &
	1960	& , fs_bound, fq_bound )
	1961	END IF
[142]	1962	PRINT *,' Lluis 1 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	1963
	1964	!c Diagnostiquer la tendance dynamique
	1965	!c
	1966	IF (ancien_ok) THEN
	1967	DO k = 1, klev
	1968	DO i = 1, klon
	1969	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
	1970	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
	1971	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	1972	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
	1973	ENDDO
	1974	ENDDO
	1975	!!! RomP >>> td dyn traceur
	1976	IF (nqtot.GE.3) THEN
	1977	DO iq = 3, nqtot
	1978	DO k = 1, klev
	1979	DO i = 1, klon
	1980	d_tr_dyn(i,k,iq-2)= &
	1981	& (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
	1982	! iiq=niadv(iq)
	1983	! print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
	1984	ENDDO
	1985	ENDDO
	1986	ENDDO
	1987	ENDIF
	1988	!!! RomP <<<
	1989	ELSE
	1990	DO k = 1, klev
	1991	DO i = 1, klon
	1992	d_u_dyn(i,k) = 0.0
	1993	d_v_dyn(i,k) = 0.0
	1994	d_t_dyn(i,k) = 0.0
	1995	d_q_dyn(i,k) = 0.0
	1996	ENDDO
	1997	ENDDO
	1998	!!! RomP >>> td dyn traceur
	1999	IF (nqtot.GE.3) THEN
	2000	DO iq = 3, nqtot
	2001	DO k = 1, klev
	2002	DO i = 1, klon
	2003	d_tr_dyn(i,k,iq-2)= 0.0
	2004	ENDDO
	2005	ENDDO
	2006	ENDDO
	2007	ENDIF
	2008	!!! RomP <<<
	2009	ancien_ok = .TRUE.
	2010	ENDIF
	2011	!c
	2012	!c Ajouter le geopotentiel du sol:
	2013	!c
	2014	DO k = 1, klev
	2015	DO i = 1, klon
	2016	zphi(i,k) = pphi(i,k) + pphis(i)
	2017	ENDDO
	2018	ENDDO
	2019	!c
	2020	!c Verifier les temperatures
	2021	!c
	2022	!IM BEG
	2023	IF (check) THEN
	2024	amn=MIN(ftsol(1,is_ter),1000.)
	2025	amx=MAX(ftsol(1,is_ter),-1000.)
	2026	DO i=2, klon
	2027	amn=MIN(ftsol(i,is_ter),amn)
	2028	amx=MAX(ftsol(i,is_ter),amx)
	2029	ENDDO
	2030	!c
	2031	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
	2032	ENDIF !(check) THEN
	2033
	2034	!IM END
	2035	!c
	2036
	2037	CALL hgardfou(t_seri,ftsol,'debutphy')
	2038	!c
	2039	!IM BEG
	2040	IF (check) THEN
	2041	amn=MIN(ftsol(1,is_ter),1000.)
	2042	amx=MAX(ftsol(1,is_ter),-1000.)
	2043	DO i=2, klon
	2044	amn=MIN(ftsol(i,is_ter),amn)
	2045	amx=MAX(ftsol(i,is_ter),amx)
	2046	ENDDO
	2047	!c
	2048	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
	2049	ENDIF !(check) THEN
	2050
	2051	!IM END
	2052	!c
	2053	!c Mettre en action les conditions aux limites (albedo, sst, etc.).
	2054	!c Prescrire l'ozone et calculer l'albedo sur l'ocean.
	2055	!c
	2056	if (read_climoz >= 1) then
	2057	!C Ozone from a file
	2058	! Update required ozone index:
	2059	ro3i = int((days_elapsed + jh_cur - jh_1jan) &
	2060	& / ioget_year_len(year_cur) * 360.) + 1
	2061	if (ro3i == 361) ro3i = 360
	2062	!C (This should never occur, except perhaps because of roundup
	2063	!C error. See documentation.)
	2064	if (ro3i /= co3i) then
	2065	!C Update ozone field:
	2066	if (read_climoz == 1) then
	2067	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
	2068	& press_in_edg=press_climoz, paprs=paprs, v3=wo)
	2069	else
	2070	!C read_climoz == 2
	2071	call regr_pr_av(ncid_climoz, &
	2072	& (/"tro3 ", "tro3_daylight"/), &
	2073	& julien=ro3i, press_in_edg=press_climoz, paprs=paprs, &
	2074	& v3=wo)
	2075	end if
	2076	! Convert from mole fraction of ozone to column density of ozone in a
	2077	! cell, in kDU:
	2078	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) &
	2079	& * rmo3 / rmd * zmasse / dobson_u / 1e3
	2080	!C (By regridding ozone values for LMDZ only once every 360th of
	2081	!C year, we have already neglected the variation of pressure in one
	2082	!C 360th of year. So do not recompute "wo" at each time step even if
	2083	!C "zmasse" changes a little.)
	2084	co3i = ro3i
	2085	end if
	2086	elseif (MOD(itap-1,lmt_pas) == 0) THEN
	2087	!C Once per day, update ozone from Royer:
	2088	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
	2089	ENDIF
	2090
	2091	!c
	2092	!c Re-evaporer l'eau liquide nuageuse
	2093	!c
	2094	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
	2095	DO i = 1, klon
	2096	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	2097	!c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	2098	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	2099	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
	2100	zb = MAX(0.0,ql_seri(i,k))
	2101	za = - MAX(0.0,ql_seri(i,k)) &
	2102	& * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
	2103	t_seri(i,k) = t_seri(i,k) + za
	2104	q_seri(i,k) = q_seri(i,k) + zb
	2105	ql_seri(i,k) = 0.0
	2106	d_t_eva(i,k) = za
	2107	d_q_eva(i,k) = zb
	2108	ENDDO
	2109	ENDDO
	2110	!IM
	2111	IF (ip_ebil_phy.ge.2) THEN
	2112	ztit='after reevap'
	2113	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
	2114	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2115	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2116	call diagphy(airephy,ztit,ip_ebil_phy &
	2117	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	2118	& , zero_v, zero_v, zero_v, ztsol &
	2119	& , d_h_vcol, d_qt, d_ec &
	2120	& , fs_bound, fq_bound )
	2121	!C
	2122	END IF
[142]	2123	PRINT *,' Lluis 2 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	2124	!c
	2125	!c=========================================================================
	2126	! Calculs de l'orbite.
	2127	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
	2128	! doit donc etre plac\'e avant radlwsw et pbl_surface
	2129
	2130	!!! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2131	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
	2132	day_since_equinox = (jD_cur + jH_cur) - jD_eq
	2133	!
	2134	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
	2135	! solarlong0
	2136	if (solarlong0<-999.) then
	2137	if (new_orbit) then
	2138	! calcul selon la routine utilisee pour les planetes
	2139	call solarlong(day_since_equinox, zlongi, dist)
	2140	else
	2141	! calcul selon la routine utilisee pour l'AR4
	2142	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
	2143	endif
	2144	else
	2145	zlongi=solarlong0 ! longitude solaire vraie
	2146	dist=1. ! distance au soleil / moyenne
	2147	endif
	2148	if(prt_level.ge.1) &
[6]	2149	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist,&
	2150	' jD_cur: ',jD_cur,' jH_cur: ',jH_cur,' days_elapsed: ', &
	2151	REAL(days_elapsed+1)
[1]	2152
	2153
	2154	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2155	! Calcul de l'ensoleillement :
	2156	! ============================
	2157	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
	2158	! l'annee a partir d'une formule analytique.
	2159	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
	2160	! non nul aux poles.
	2161	IF (abs(solarlong0-1000.)<1.e-4) then
	2162	call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
	2163	ELSE
	2164	! Avec ou sans cycle diurne
	2165	IF (cycle_diurne) THEN
	2166	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
	2167	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
	2168	ELSE
	2169	CALL angle(zlongi, rlat, fract, rmu0)
	2170	ENDIF
	2171	ENDIF
	2172
	2173	if (mydebug) then
	2174	call writefield_phy('u_seri',u_seri,llm)
	2175	call writefield_phy('v_seri',v_seri,llm)
	2176	call writefield_phy('t_seri',t_seri,llm)
	2177	call writefield_phy('q_seri',q_seri,llm)
	2178	endif
	2179
	2180	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2181	!c Appel au pbl_surface : Planetary Boudary Layer et Surface
	2182	!c Cela implique tous les interactions des sous-surfaces et la partie diffusion
	2183	!c turbulent du couche limit.
	2184	!c
	2185	!c Certains varibales de sorties de pbl_surface sont utiliser que pour
	2186	!c ecriture des fihiers hist_XXXX.nc, ces sont :
	2187	!c qsol, zq2m, s_pblh, s_lcl,
	2188	!c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
	2189	!c s_therm, s_trmb1, s_trmb2, s_trmb3,
	2190	!c zxrugs, zu10m, zv10m, fder,
	2191	!c zxqsurf, rh2m, zxfluxu, zxfluxv,
	2192	!c frugs, agesno, fsollw, fsolsw,
	2193	!c d_ts, fevap, fluxlat, t2m,
	2194	!c wfbils, wfbilo, fluxt, fluxu, fluxv,
	2195	!c
	2196	!c Certains ne sont pas utiliser du tout :
	2197	!c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
	2198	!c
	2199
	2200	!c Calcul de l'humidite de saturation au niveau du sol
	2201
	2202
	2203
	2204	if (iflag_pbl/=0) then
	2205
	2206	CALL pbl_surface( &
	2207	& dtime, date0, itap, days_elapsed+1, &
	2208	& debut, lafin, &
	2209	& rlon, rlat, rugoro, rmu0, &
	2210	& rain_fall, snow_fall, solsw, sollw, &
	2211	& t_seri, q_seri, u_seri, v_seri, &
	2212	& pplay, paprs, pctsrf, &
	2213	& ftsol, falb1, falb2, ustar, u10m, v10m, &
	2214	& sollwdown, cdragh, cdragm, u1, v1, &
	2215	& albsol1, albsol2, sens, evap, &
	2216	& zxtsol, zxfluxlat, zt2m, qsat2m, &
	2217	& d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
	2218	& coefh, coefm, slab_wfbils, &
	2219	& qsol, zq2m, s_pblh, s_lcl, &
	2220	& s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
	2221	& s_therm, s_trmb1, s_trmb2, s_trmb3, &
	2222	& zxrugs, zustar, zu10m, zv10m, fder, &
	2223	& zxqsurf, rh2m, zxfluxu, zxfluxv, &
	2224	& frugs, agesno, fsollw, fsolsw, &
	2225	& d_ts, fevap, fluxlat, t2m, &
	2226	& wfbils, wfbilo, fluxt, fluxu, fluxv, &
	2227	& dsens, devap, zxsnow, &
	2228	& zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
	2229
	2230	!-----------------------------------------------------------------------------------------
	2231	! ajout des tendances de la diffusion turbulente
	2232	CALL add_phys_tend &
	2233	& (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,'vdf')
	2234	!-----------------------------------------------------------------------------------------
	2235	if (mydebug) then
	2236	call writefield_phy('u_seri',u_seri,llm)
	2237	call writefield_phy('v_seri',v_seri,llm)
	2238	call writefield_phy('t_seri',t_seri,llm)
	2239	call writefield_phy('q_seri',q_seri,llm)
	2240	endif
	2241
	2242	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
	2243	& t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
	2244
	2245
	2246	IF (ip_ebil_phy.ge.2) THEN
	2247	ztit='after surface_main'
	2248	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	2249	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2250	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2251	call diagphy(airephy,ztit,ip_ebil_phy &
	2252	& , zero_v, zero_v, zero_v, zero_v, sens &
	2253	& , evap , zero_v, zero_v, ztsol &
	2254	& , d_h_vcol, d_qt, d_ec &
	2255	& , fs_bound, fq_bound )
	2256	END IF
	2257
	2258	ENDIF
[142]	2259	PRINT *,' Lluis 3 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	2260	!c =================================================================== c
	2261	!c Calcul de Qsat
	2262
	2263	DO k = 1, klev
	2264	DO i = 1, klon
	2265	zx_t = t_seri(i,k)
	2266	IF (thermcep) THEN
	2267	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	2268	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	2269	zx_qs = MIN(0.5,zx_qs)
	2270	zcor = 1./(1.-retv*zx_qs)
	2271	zx_qs = zx_qs*zcor
	2272	ELSE
	2273	IF (zx_t.LT.t_coup) THEN
	2274	zx_qs = qsats(zx_t)/pplay(i,k)
	2275	ELSE
	2276	zx_qs = qsatl(zx_t)/pplay(i,k)
	2277	ENDIF
	2278	ENDIF
	2279	zqsat(i,k)=zx_qs
	2280	ENDDO
	2281	ENDDO
	2282
	2283	if (prt_level.ge.1) then
	2284	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
	2285	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
	2286	endif
	2287	!c
	2288	!c Appeler la convection (au choix)
	2289	!c
	2290	DO k = 1, klev
	2291	DO i = 1, klon
	2292	conv_q(i,k) = d_q_dyn(i,k) &
	2293	& + d_q_vdf(i,k)/dtime
	2294	conv_t(i,k) = d_t_dyn(i,k) &
	2295	& + d_t_vdf(i,k)/dtime
	2296	ENDDO
	2297	ENDDO
	2298	IF (check) THEN
	2299	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2300	WRITE(lunout,*) "avantcon=", za
	2301	ENDIF
	2302	zx_ajustq = .FALSE.
	2303	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
	2304	IF (zx_ajustq) THEN
	2305	DO i = 1, klon
	2306	z_avant(i) = 0.0
	2307	ENDDO
	2308	DO k = 1, klev
	2309	DO i = 1, klon
	2310	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
	2311	& *(paprs(i,k)-paprs(i,k+1))/RG
	2312	ENDDO
	2313	ENDDO
	2314	ENDIF
	2315
	2316	!c Calcule de vitesse verticale a partir de flux de masse verticale
	2317	DO k = 1, klev
	2318	DO i = 1, klon
	2319	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
	2320	END DO
	2321	END DO
	2322	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
	2323	& omega(igout, :)
	2324
	2325	IF (iflag_con.EQ.1) THEN
	2326	abort_message ='reactiver le call conlmd dans physiq.F'
	2327	CALL abort_gcm (modname,abort_message,1)
	2328	!c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
	2329	!c . d_t_con, d_q_con,
	2330	!c . rain_con, snow_con, ibas_con, itop_con)
	2331	ELSE IF (iflag_con.EQ.2) THEN
	2332	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
	2333	& conv_t, conv_q, -evap, omega, &
	2334	& d_t_con, d_q_con, rain_con, snow_con, &
	2335	& pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
	2336	& kcbot, kctop, kdtop, pmflxr, pmflxs)
	2337	d_u_con = 0.
	2338	d_v_con = 0.
	2339
	2340	WHERE (rain_con < 0.) rain_con = 0.
	2341	WHERE (snow_con < 0.) snow_con = 0.
	2342	DO i = 1, klon
	2343	ibas_con(i) = klev+1 - kcbot(i)
	2344	itop_con(i) = klev+1 - kctop(i)
	2345	ENDDO
	2346	ELSE IF (iflag_con.GE.3) THEN
	2347	!c nb of tracers for the KE convection:
	2348	!c MAF la partie traceurs est faite dans phytrac
	2349	!c on met ntra=1 pour limiter les appels mais on peut
	2350	!c supprimer les calculs / ftra.
	2351	ntra = 1
	2352
	2353	!c=====================================================================================
	2354	!cajout pour la parametrisation des poches froides:
	2355	!ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
	2356	do k=1,klev
	2357	do i=1,klon
	2358	if (iflag_wake>=1) then
	2359	t_wake(i,k) = t_seri(i,k) &
	2360	& +(1-wake_s(i))*wake_deltat(i,k)
	2361	q_wake(i,k) = q_seri(i,k) &
	2362	& +(1-wake_s(i))*wake_deltaq(i,k)
	2363	t_undi(i,k) = t_seri(i,k) &
	2364	& -wake_s(i)*wake_deltat(i,k)
	2365	q_undi(i,k) = q_seri(i,k) &
	2366	& -wake_s(i)*wake_deltaq(i,k)
	2367	else
	2368	t_wake(i,k) = t_seri(i,k)
	2369	q_wake(i,k) = q_seri(i,k)
	2370	t_undi(i,k) = t_seri(i,k)
	2371	q_undi(i,k) = q_seri(i,k)
	2372	endif
	2373	enddo
	2374	enddo
	2375
	2376	!cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
	2377	!cc-- pour le soulevement des particules dans le modele convectif
	2378	!c
	2379	do i = 1,klon
	2380	ALE(i) = 0.
	2381	ALP(i) = 0.
	2382	enddo
	2383	!c
	2384	!ccalcul de ale_wake et alp_wake
	2385	if (iflag_wake>=1) then
	2386	if (itap .le. it_wape_prescr) then
	2387	do i = 1,klon
	2388	ale_wake(i) = wape_prescr
	2389	alp_wake(i) = fip_prescr
	2390	enddo
	2391	else
	2392	do i = 1,klon
	2393	!cjyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
	2394	!ccc ale_wake(i) = 0.5wake_cstar(i)*2
	2395	ale_wake(i) = wake_pe(i)
	2396	alp_wake(i) = wake_fip(i)
	2397	enddo
	2398	endif
	2399	else
	2400	do i = 1,klon
	2401	ale_wake(i) = 0.
	2402	alp_wake(i) = 0.
	2403	enddo
	2404	endif
	2405	!ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
	2406	!cdans le thermique sinon
	2407	if (iflag_coupl.eq.0) then
	2408	if (debut.and.prt_level.gt.9) &
	2409	& WRITE(lunout,*)'ALE et ALP imposes'
	2410	do i = 1,klon
	2411	!con ne couple que ale
	2412	!c ALE(i) = max(ale_wake(i),Ale_bl(i))
	2413	ALE(i) = max(ale_wake(i),ale_bl_prescr)
	2414	!con ne couple que alp
	2415	!c ALP(i) = alp_wake(i) + Alp_bl(i)
	2416	ALP(i) = alp_wake(i) + alp_bl_prescr
	2417	enddo
	2418	else
	2419	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
	2420	! do i = 1,klon
	2421	! ALE(i) = max(ale_wake(i),Ale_bl(i))
	2422	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
	2423	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
	2424	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
	2425	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
	2426	! enddo
	2427
	2428	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2429	! Modif FH 2010/04/27. Sans doute temporaire.
	2430	! Deux options pour le alp_offset : constant si >?? 0 ou proportionnel ??a
	2431	! w si <0
	2432	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2433	do i = 1,klon
	2434	ALE(i) = max(ale_wake(i),Ale_bl(i))
	2435	!ccc nrlmd le 10/04/2012----------Stochastic triggering--------------
	2436	if (iflag_trig_bl.ge.1) then
	2437	ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
	2438	endif
	2439	!ccc fin nrlmd le 10/04/2012
	2440	if (alp_offset>=0.) then
	2441	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
	2442	else
	2443	ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
	2444	if (alp(i)<0.) then
	2445	print*,'ALP ',alp(i),alp_wake(i) &
	2446	& ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
	2447	endif
	2448	endif
	2449	enddo
	2450	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2451
	2452	endif
	2453	do i=1,klon
	2454	if (alp(i)>alp_max) then
	2455	IF(prt_level>9)WRITE(lunout,*) &
	2456	& 'WARNING SUPER ALP (seuil=',alp_max, &
	2457	& '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
	2458	alp(i)=alp_max
	2459	endif
	2460	if (ale(i)>ale_max) then
	2461	IF(prt_level>9)WRITE(lunout,*) &
	2462	& 'WARNING SUPER ALE (seuil=',ale_max, &
	2463	& '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
	2464	ale(i)=ale_max
	2465	endif
	2466	enddo
	2467
	2468	!cfin calcul ale et alp
	2469	!c=================================================================================================
	2470
	2471	!c sb, oct02:
	2472	!c Schema de convection modularise et vectorise:
	2473	!c (driver commun aux versions 3 et 4)
	2474	!c
	2475	IF (ok_cvl) THEN ! new driver for convectL
	2476
	2477	IF (type_trac == 'repr') THEN
	2478	nbtr_tmp=ntra
	2479	ELSE
	2480	nbtr_tmp=nbtr
	2481	END IF
	2482	CALL concvl (iflag_con,iflag_clos, &
	2483	& dtime,paprs,pplay,t_undi,q_undi, &
	2484	& t_wake,q_wake,wake_s, &
	2485	& u_seri,v_seri,tr_seri,nbtr_tmp, &
	2486	& ALE,ALP, &
	2487	& ema_work1,ema_work2, &
	2488	& d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
	2489	& rain_con, snow_con, ibas_con, itop_con, sigd, &
	2490	& ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
	2491	& Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
	2492	& pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
	2493	! RomP >>>
	2494	!! . pmflxr,pmflxs,da,phi,mp,
	2495	!! . ftd,fqd,lalim_conv,wght_th) &
	2496	& pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
	2497	& ftd,fqd,lalim_conv,wght_th, &
	2498	& ev, ep,epmlmMm,eplaMm, &
	2499	& wdtrainA,wdtrainM)
	2500	! RomP <<<
	2501
	2502	!IM begin
	2503	!c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
	2504	!c .dnwd0(1,1),ftd(1,1),fqd(1,1)
	2505	!IM end
	2506	!IM cf. FH
	2507	clwcon0=qcondc
	2508	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
	2509
	2510	do i = 1, klon
	2511	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
	2512	enddo
	2513
	2514	ELSE ! ok_cvl
	2515
	2516	!c MAF conema3 ne contient pas les traceurs
	2517	CALL conema3 (dtime, &
	2518	& paprs,pplay,t_seri,q_seri, &
	2519	& u_seri,v_seri,tr_seri,ntra, &
	2520	& ema_work1,ema_work2, &
	2521	& d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
	2522	& rain_con, snow_con, ibas_con, itop_con, &
	2523	& upwd,dnwd,dnwd0,bas,top, &
	2524	& Ma,cape,tvp,rflag, &
	2525	& pbase &
	2526	& ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
	2527	& ,clwcon0)
	2528
	2529	ENDIF ! ok_cvl
	2530
	2531	!c
	2532	!c Correction precip
	2533	rain_con = rain_con * cvl_corr
	2534	snow_con = snow_con * cvl_corr
	2535	!c
	2536
	2537	IF (.NOT. ok_gust) THEN
	2538	do i = 1, klon
	2539	wd(i)=0.0
	2540	enddo
	2541	ENDIF
	2542
	2543	!c =================================================================== c
	2544	!c Calcul des proprietes des nuages convectifs
	2545	!c
	2546
	2547	!c calcul des proprietes des nuages convectifs
	2548	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
	2549	call clouds_gno &
	2550	& (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
	2551
	2552	!c =================================================================== c
	2553
	2554	DO i = 1, klon
	2555	itop_con(i) = min(max(itop_con(i),1),klev)
	2556	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
	2557	ENDDO
	2558
	2559	DO i = 1, klon
	2560	ema_pcb(i) = paprs(i,ibas_con(i))
	2561	ENDDO
	2562	DO i = 1, klon
	2563	! L'idicage de itop_con peut cacher un pb potentiel
	2564	! FH sous la dictee de JYG, CR
	2565	ema_pct(i) = paprs(i,itop_con(i)+1)
	2566
	2567	if (itop_con(i).gt.klev-3) then
	2568	if(prt_level >= 9) then
	2569	write(lunout,*)'La convection monte trop haut '
	2570	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
	2571	endif
	2572	endif
	2573	ENDDO
	2574	ELSE IF (iflag_con.eq.0) THEN
	2575	write(lunout,*) 'On n appelle pas la convection'
	2576	clwcon0=0.
	2577	rnebcon0=0.
	2578	d_t_con=0.
	2579	d_q_con=0.
	2580	d_u_con=0.
	2581	d_v_con=0.
	2582	rain_con=0.
	2583	snow_con=0.
	2584	bas=1
	2585	top=1
	2586	ELSE
	2587	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
	2588	CALL abort
	2589	ENDIF
	2590
	2591	!c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
	2592	!c . d_u_con, d_v_con)
	2593
	2594	!-----------------------------------------------------------------------------------------
	2595	! ajout des tendances de la diffusion turbulente
	2596	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
	2597	!-----------------------------------------------------------------------------------------
	2598
	2599	if (mydebug) then
	2600	call writefield_phy('u_seri',u_seri,llm)
	2601	call writefield_phy('v_seri',v_seri,llm)
	2602	call writefield_phy('t_seri',t_seri,llm)
	2603	call writefield_phy('q_seri',q_seri,llm)
	2604	endif
	2605
	2606	!IM
	2607	IF (ip_ebil_phy.ge.2) THEN
	2608	ztit='after convect'
	2609	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	2610	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2611	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2612	call diagphy(airephy,ztit,ip_ebil_phy &
	2613	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	2614	& , zero_v, rain_con, snow_con, ztsol &
	2615	& , d_h_vcol, d_qt, d_ec &
	2616	& , fs_bound, fq_bound )
	2617	END IF
	2618	!C
[142]	2619	PRINT *,' Lluis 4 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	2620	IF (check) THEN
	2621	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2622	WRITE(lunout,*)"aprescon=", za
	2623	zx_t = 0.0
	2624	za = 0.0
	2625	DO i = 1, klon
	2626	za = za + airephy(i)/REAL(klon)
	2627	zx_t = zx_t + (rain_con(i)+ &
	2628	& snow_con(i))*airephy(i)/REAL(klon)
	2629	ENDDO
	2630	zx_t = zx_t/za*dtime
	2631	WRITE(lunout,*)"Precip=", zx_t
	2632	ENDIF
	2633	IF (zx_ajustq) THEN
	2634	DO i = 1, klon
	2635	z_apres(i) = 0.0
	2636	ENDDO
	2637	DO k = 1, klev
	2638	DO i = 1, klon
	2639	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
	2640	& *(paprs(i,k)-paprs(i,k+1))/RG
	2641	ENDDO
	2642	ENDDO
	2643	DO i = 1, klon
	2644	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
	2645	& /z_apres(i)
	2646	ENDDO
	2647	DO k = 1, klev
	2648	DO i = 1, klon
	2649	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
	2650	& z_factor(i).LT.(1.0-1.0E-08)) THEN
	2651	q_seri(i,k) = q_seri(i,k) * z_factor(i)
	2652	ENDIF
	2653	ENDDO
	2654	ENDDO
	2655	ENDIF
	2656	zx_ajustq=.FALSE.
	2657
	2658	!c
	2659	!c=============================================================================
	2660	!cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
	2661	!cpour la couche limite diffuse pour l instant
	2662	!c
	2663	if (iflag_wake>=1) then
	2664	DO k=1,klev
	2665	DO i=1,klon
	2666	dt_dwn(i,k) = ftd(i,k)
	2667	wdt_PBL(i,k) = 0.
	2668	dq_dwn(i,k) = fqd(i,k)
	2669	wdq_PBL(i,k) = 0.
	2670	M_dwn(i,k) = dnwd0(i,k)
	2671	M_up(i,k) = upwd(i,k)
	2672	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
	2673	udt_PBL(i,k) = 0.
	2674	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
	2675	udq_PBL(i,k) = 0.
	2676	ENDDO
	2677	ENDDO
	2678
	2679	if (iflag_wake==2) then
	2680	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
	2681	DO k = 1,klev
	2682	dt_dwn(:,k)= dt_dwn(:,k)+ &
	2683	& ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
	2684	dq_dwn(:,k)= dq_dwn(:,k)+ &
	2685	& ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
	2686	ENDDO
	2687	endif
	2688	!c
	2689	!ccalcul caracteristiques de la poche froide
	2690	call calWAKE (paprs,pplay,dtime &
	2691	& ,t_seri,q_seri,omega &
	2692	& ,dt_dwn,dq_dwn,M_dwn,M_up &
	2693	& ,dt_a,dq_a,sigd &
	2694	& ,wdt_PBL,wdq_PBL &
	2695	& ,udt_PBL,udq_PBL &
	2696	& ,wake_deltat,wake_deltaq,wake_dth &
	2697	& ,wake_h,wake_s,wake_dens &
	2698	& ,wake_pe,wake_fip,wake_gfl &
	2699	& ,dt_wake,dq_wake &
	2700	& ,wake_k, t_undi,q_undi &
	2701	& ,wake_omgbdth,wake_dp_omgb &
	2702	& ,wake_dtKE,wake_dqKE &
	2703	& ,wake_dtPBL,wake_dqPBL &
	2704	& ,wake_omg,wake_dp_deltomg &
	2705	& ,wake_spread,wake_Cstar,wake_d_deltat_gw &
	2706	& ,wake_ddeltat,wake_ddeltaq)
	2707	!c
	2708	!-----------------------------------------------------------------------------------------
	2709	! ajout des tendances des poches froides
	2710	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
	2711	! coherent avec les autres d_t_...
	2712	d_t_wake(:,:)=dt_wake(:,:)*dtime
	2713	d_q_wake(:,:)=dq_wake(:,:)*dtime
	2714	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
	2715	!-----------------------------------------------------------------------------------------
	2716
	2717	endif
	2718	!c
	2719	!c===================================================================
	2720	!cJYG
	2721	IF (ip_ebil_phy.ge.2) THEN
	2722	ztit='after wake'
	2723	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	2724	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2725	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2726	call diagphy(airephy,ztit,ip_ebil_phy &
	2727	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	2728	& , zero_v, zero_v, zero_v, ztsol &
	2729	& , d_h_vcol, d_qt, d_ec &
	2730	& , fs_bound, fq_bound )
	2731	END IF
[142]	2732	PRINT *,' Lluis 5 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	2733	!c print*,'apres callwake iflag_cldcon=', iflag_cldcon
	2734	!c
	2735	!c===================================================================
	2736	!c Convection seche (thermiques ou ajustement)
	2737	!c===================================================================
	2738	!c
	2739	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
	2740	& ,seuil_inversion,weak_inversion,dthmin)
	2741
	2742
	2743
	2744	d_t_ajsb(:,:)=0.
	2745	d_q_ajsb(:,:)=0.
	2746	d_t_ajs(:,:)=0.
	2747	d_u_ajs(:,:)=0.
	2748	d_v_ajs(:,:)=0.
	2749	d_q_ajs(:,:)=0.
	2750	clwcon0th(:,:)=0.
	2751	!c
	2752	!c fm_therm(:,:)=0.
	2753	!c entr_therm(:,:)=0.
	2754	!c detr_therm(:,:)=0.
	2755	!c
[142]	2756
[1]	2757	IF(prt_level>9)WRITE(lunout,*) &
	2758	& 'AVANT LA CONVECTION SECHE , iflag_thermals=' &
	2759	& ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
	2760	if(iflag_thermals.lt.0) then
	2761	!c Rien
	2762	!c ====
	2763	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
	2764
	2765
	2766	else
	2767
	2768	!c Thermiques
	2769	!c ==========
	2770	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
	2771	& ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
	2772
	2773
	2774	!ccc nrlmd le 10/04/2012
	2775	DO k=1,klev+1
	2776	DO i=1,klon
	2777	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
	2778	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
	2779	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
	2780	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
	2781	ENDDO
	2782	ENDDO
	2783	!ccc fin nrlmd le 10/04/2012
	2784
	2785	if (iflag_thermals>=1) then
	2786	call calltherm(pdtphys &
	2787	& ,pplay,paprs,pphi,weak_inversion &
	2788	& ,u_seri,v_seri,t_seri,q_seri,zqsat,debut &
	2789	& ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
	2790	& ,fm_therm,entr_therm,detr_therm &
	2791	& ,zqasc,clwcon0th,lmax_th,ratqscth &
	2792	& ,ratqsdiff,zqsatth &
	2793	!con rajoute ale et alp, et les caracteristiques de la couche alim &
	2794	& ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
	2795	& ,ztv,zpspsk,ztla,zthl &
	2796	!ccc nrlmd le 10/04/2012 &
	2797	& ,pbl_tke_input,pctsrf,omega,airephy &
	2798	& ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
	2799	& ,n2,s2,ale_bl_stat &
	2800	& ,therm_tke_max,env_tke_max &
	2801	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
	2802	& ,alp_bl_conv,alp_bl_stat &
	2803	!ccc fin nrlmd le 10/04/2012 &
	2804	& ,zqla,ztva )
	2805
	2806	!ccc nrlmd le 10/04/2012
	2807	!c-----------Stochastic triggering-----------
	2808	if (iflag_trig_bl.ge.1) then
	2809	!c
	2810	IF (prt_level .GE. 10) THEN
	2811	print *,'cin, ale_bl_stat, alp_bl_stat ', &
	2812	& cin, ale_bl_stat, alp_bl_stat
	2813	ENDIF
	2814
	2815	!c----Initialisations
	2816	do i=1,klon
	2817	proba_notrig(i)=1.
	2818	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
	2819	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
	2820	tau_trig(i)=tau_trig_shallow
	2821	else
	2822	tau_trig(i)=tau_trig_deep
	2823	endif
	2824	enddo
	2825	!c
	2826	IF (prt_level .GE. 10) THEN
	2827	print *,'random_notrig, tau_trig ', &
	2828	& random_notrig, tau_trig
	2829	print *,'s_trig,s2,n2 ', &
	2830	& s_trig,s2,n2
	2831	ENDIF
	2832
	2833	!c----Tirage al\'eatoire et calcul de ale_bl_trig
	2834	do i=1,klon
	2835	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
	2836	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
	2837	& (n2(i)*dtime/tau_trig(i))
	2838	!c print *, 'proba_notrig(i) ',proba_notrig(i)
	2839	if (random_notrig(i) .ge. proba_notrig(i)) then
	2840	ale_bl_trig(i)=ale_bl_stat(i)
	2841	else
	2842	ale_bl_trig(i)=0.
	2843	endif
	2844	else
	2845	proba_notrig(i)=1.
	2846	random_notrig(i)=0.
	2847	ale_bl_trig(i)=0.
	2848	endif
	2849	enddo
	2850	!c
	2851	IF (prt_level .GE. 10) THEN
	2852	print *,'proba_notrig, ale_bl_trig ', &
	2853	& proba_notrig, ale_bl_trig
	2854	ENDIF
	2855
	2856	endif !(iflag_trig_bl)
	2857
	2858	!c-----------Statistical closure-----------
	2859	if (iflag_clos_bl.ge.1) then
	2860
	2861	do i=1,klon
	2862	alp_bl(i)=alp_bl_stat(i)
	2863	enddo
	2864
	2865	else
	2866
	2867	alp_bl_stat(:)=0.
	2868
	2869	endif !(iflag_clos_bl)
	2870
	2871	IF (prt_level .GE. 10) THEN
	2872	print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
	2873	ENDIF
	2874
	2875	!ccc fin nrlmd le 10/04/2012
	2876
	2877	! ----------------------------------------------------------------------
	2878	! Transport de la TKE par les panaches thermiques.
	2879	! FH : 2010/02/01
	2880	! if (iflag_pbl.eq.10) then
	2881	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
	2882	! s rg,paprs,pbl_tke)
	2883	! endif
	2884	! ----------------------------------------------------------------------
	2885	!IM/FH: 2011/02/23
	2886	! Couplage Thermiques/Emanuel seulement si T<0
	2887	if (iflag_coupl==2) then
	2888	print*,'Couplage Thermiques/Emanuel seulement si T<0'
	2889	do i=1,klon
	2890	if (t_seri(i,lmax_th(i))>273.) then
	2891	Ale_bl(i)=0.
	2892	endif
	2893	enddo
	2894	endif
	2895
	2896	do i=1,klon
	2897	zmax_th(i)=pphi(i,lmax_th(i))/rg
	2898	enddo
	2899
	2900	endif
	2901
	2902
	2903	!c Ajustement sec
	2904	!c ==============
	2905
	2906	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
	2907	! a partir du sommet des thermiques.
	2908	! Dans le cas contraire, on demarre au niveau 1.
	2909
	2910	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
	2911
	2912	if(iflag_thermals.eq.0) then
	2913	IF(prt_level>9)WRITE(lunout,*)'ajsec'
	2914	limbas(:)=1
	2915	else
	2916	limbas(:)=lmax_th(:)
	2917	endif
	2918
	2919	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
	2920	! pour des test de convergence numerique.
	2921	! Le nouveau ajsec est a priori mieux, meme pour le cas
	2922	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
	2923	! non nulles numeriquement pour des mailles non concernees.
	2924
	2925	if (iflag_thermals.eq.0) then
	2926	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
	2927	& , d_t_ajsb, d_q_ajsb)
	2928	else
	2929	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
	2930	& , d_t_ajsb, d_q_ajsb)
	2931	endif
	2932
	2933	!-----------------------------------------------------------------------------------------
	2934	! ajout des tendances de l'ajustement sec ou des thermiques
	2935	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
	2936	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
	2937	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
	2938
	2939	!-----------------------------------------------------------------------------------------
	2940
	2941	endif
	2942
	2943	endif
[142]	2944
[1]	2945	!c
	2946	!c===================================================================
	2947	!IM
	2948	IF (ip_ebil_phy.ge.2) THEN
	2949	ztit='after dry_adjust'
	2950	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	2951	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2952	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2953	call diagphy(airephy,ztit,ip_ebil_phy &
	2954	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	2955	& , zero_v, zero_v, zero_v, ztsol &
	2956	& , d_h_vcol, d_qt, d_ec &
	2957	& , fs_bound, fq_bound )
	2958	END IF
[142]	2959	PRINT *,' Lluis 6 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	2960
	2961	!c-------------------------------------------------------------------------
	2962	! Computation of ratqs, the width (normalized) of the subrid scale
	2963	! water distribution
	2964	CALL calcratqs(klon,klev,prt_level,lunout, &
	2965	& iflag_ratqs,iflag_con,iflag_cldcon,pdtphys, &
	2966	& ratqsbas,ratqshaut,tau_ratqs,fact_cldcon, &
	2967	& ptconv,ptconvth,clwcon0th, rnebcon0th, &
	2968	& paprs,pplay,q_seri,zqsat,fm_therm, &
	2969	& ratqs,ratqsc)
	2970
	2971	!c
	2972	!c Appeler le processus de condensation a grande echelle
	2973	!c et le processus de precipitation
	2974	!c-------------------------------------------------------------------------
	2975	IF (prt_level .GE.10) THEN
	2976	print *,' ->fisrtilp '
	2977	ENDIF
	2978	!c-------------------------------------------------------------------------
	2979	CALL fisrtilp(dtime,paprs,pplay, &
	2980	& t_seri, q_seri,ptconv,ratqs, &
	2981	& d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
	2982	& rain_lsc, snow_lsc, &
	2983	& pfrac_impa, pfrac_nucl, pfrac_1nucl, &
	2984	& frac_impa, frac_nucl, beta_prec_fisrt, &
	2985	& prfl, psfl, rhcl, &
	2986	& zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
	2987
	2988	WHERE (rain_lsc < 0) rain_lsc = 0.
	2989	WHERE (snow_lsc < 0) snow_lsc = 0.
	2990	!-----------------------------------------------------------------------------------------
	2991	! ajout des tendances de la diffusion turbulente
	2992	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
	2993	!-----------------------------------------------------------------------------------------
	2994	DO k = 1, klev
	2995	DO i = 1, klon
	2996	cldfra(i,k) = rneb(i,k)
	2997	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
	2998	ENDDO
	2999	ENDDO
	3000	IF (check) THEN
	3001	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	3002	WRITE(lunout,*)"apresilp=", za
	3003	zx_t = 0.0
	3004	za = 0.0
	3005	DO i = 1, klon
	3006	za = za + airephy(i)/REAL(klon)
	3007	zx_t = zx_t + (rain_lsc(i) &
	3008	& + snow_lsc(i))*airephy(i)/REAL(klon)
	3009	ENDDO
	3010	zx_t = zx_t/za*dtime
	3011	WRITE(lunout,*)"Precip=", zx_t
	3012	ENDIF
	3013	!IM
	3014	IF (ip_ebil_phy.ge.2) THEN
	3015	ztit='after fisrt'
	3016	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	3017	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3018	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3019	call diagphy(airephy,ztit,ip_ebil_phy &
	3020	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	3021	& , zero_v, rain_lsc, snow_lsc, ztsol &
	3022	& , d_h_vcol, d_qt, d_ec &
	3023	& , fs_bound, fq_bound )
	3024	END IF
[142]	3025	PRINT *,' Lluis 7 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	3026	if (mydebug) then
	3027	call writefield_phy('u_seri',u_seri,llm)
	3028	call writefield_phy('v_seri',v_seri,llm)
	3029	call writefield_phy('t_seri',t_seri,llm)
	3030	call writefield_phy('q_seri',q_seri,llm)
	3031	endif
	3032	!c
	3033	!c-------------------------------------------------------------------
	3034	!c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
	3035	!c-------------------------------------------------------------------
	3036
	3037	!c 1. NUAGES CONVECTIFS
	3038	!c
	3039	!IM cf FH
	3040	!c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
	3041	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
	3042	snow_tiedtke=0.
	3043	!c print*,'avant calcul de la pseudo precip '
	3044	!c print*,'iflag_cldcon',iflag_cldcon
	3045	if (iflag_cldcon.eq.-1) then
	3046	rain_tiedtke=rain_con
	3047	else
	3048	!c print*,'calcul de la pseudo precip '
	3049	rain_tiedtke=0.
	3050	!c print*,'calcul de la pseudo precip 0'
	3051	do k=1,klev
	3052	do i=1,klon
	3053	if (d_q_con(i,k).lt.0.) then
	3054	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
	3055	& *(paprs(i,k)-paprs(i,k+1))/rg
	3056	endif
	3057	enddo
	3058	enddo
	3059	endif
	3060	!c
	3061	!c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
	3062	!c
	3063
	3064	!c Nuages diagnostiques pour Tiedtke
	3065	CALL diagcld1(paprs,pplay, &
	3066	!IM cf FH . rain_con,snow_con,ibas_con,itop_con, &
	3067	& rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
	3068	& diafra,dialiq)
	3069	DO k = 1, klev
	3070	DO i = 1, klon
	3071	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	3072	cldliq(i,k) = dialiq(i,k)
	3073	cldfra(i,k) = diafra(i,k)
	3074	ENDIF
	3075	ENDDO
	3076	ENDDO
	3077
	3078	ELSE IF (iflag_cldcon.ge.3) THEN
	3079	!c On prend pour les nuages convectifs le max du calcul de la
	3080	!c convection et du calcul du pas de temps precedent diminue d'un facteur
	3081	!c facttemps
	3082	facteur = pdtphys *facttemps
	3083	do k=1,klev
	3084	do i=1,klon
	3085	rnebcon(i,k)=rnebcon(i,k)*facteur
	3086	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k)) &
	3087	& then
	3088	rnebcon(i,k)=rnebcon0(i,k)
	3089	clwcon(i,k)=clwcon0(i,k)
	3090	endif
	3091	enddo
	3092	enddo
	3093
	3094	!c
	3095	!cjq - introduce the aerosol direct and first indirect radiative forcings
	3096	!cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
	3097	IF (flag_aerosol .gt. 0) THEN
	3098	IF (.NOT. aerosol_couple) &
	3099	& CALL readaerosol_optic( &
	3100	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
	3101	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3102	& mass_solu_aero, mass_solu_aero_pi, &
	3103	& tau_aero, piz_aero, cg_aero, &
	3104	& tausum_aero, tau3d_aero)
	3105	ELSE
	3106	tausum_aero(:,:,:) = 0.
	3107	tau_aero(:,:,:,:) = 0.
	3108	piz_aero(:,:,:,:) = 0.
	3109	cg_aero(:,:,:,:) = 0.
	3110	ENDIF
	3111	!c
	3112	!c--STRAT AEROSOL
	3113	!c--updates tausum_aero,tau_aero,piz_aero,cg_aero
	3114	IF (flag_aerosol_strat) THEN
	3115	PRINT *,'appel a readaerosolstrat', mth_cur
	3116	CALL readaerosolstrato(debut)
	3117	ENDIF
	3118	!c--fin STRAT AEROSOL
	3119
	3120	!IM calcul nuages par le simulateur ISCCP
	3121	!c
	3122	#ifdef histISCCP
	3123	IF (ok_isccp) THEN
	3124	!c
	3125	!IM lecture invtau, tautab des fichiers formattes
	3126	!c
	3127	IF (debut) THEN
	3128	!$OMP MASTER
	3129	!c
	3130	open(99,file='tautab.formatted', FORM='FORMATTED')
	3131	read(99,'(f30.20)') tautab_omp
	3132	close(99)
	3133	!c
	3134	open(99,file='invtau.formatted',form='FORMATTED')
	3135	read(99,'(i10)') invtau_omp
	3136
	3137	!c print*,'calcul_simulISCCP invtau_omp',invtau_omp
	3138	!c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
	3139
	3140	close(99)
	3141	!$OMP END MASTER
	3142	!$OMP BARRIER
	3143	tautab=tautab_omp
	3144	invtau=invtau_omp
	3145	!c
	3146	ENDIF !debut
	3147	!c
	3148	!IM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
	3149	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
	3150	#include "calcul_simulISCCP.h"
	3151	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
	3152	ENDIF !ok_isccp
	3153	#endif
	3154
	3155	!c On prend la somme des fractions nuageuses et des contenus en eau
	3156
	3157	if (iflag_cldcon>=5) then
	3158
	3159	do k=1,klev
	3160	ptconvth(:,k)=fm_therm(:,k+1)>0.
	3161	enddo
	3162
	3163	if (iflag_coupl==4) then
	3164
	3165	! Dans le cas iflag_coupl==4, on prend la somme des convertures
	3166	! convectives et lsc dans la partie des thermiques
	3167	! Le controle par iflag_coupl est peut etre provisoire.
	3168	do k=1,klev
	3169	do i=1,klon
	3170	if (ptconv(i,k).and.ptconvth(i,k)) then
	3171	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
	3172	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
	3173	else if (ptconv(i,k)) then
	3174	cldfra(i,k)=rnebcon(i,k)
	3175	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
	3176	endif
	3177	enddo
	3178	enddo
	3179
	3180	else if (iflag_coupl==5) then
	3181	do k=1,klev
	3182	do i=1,klon
	3183	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
	3184	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
	3185	enddo
	3186	enddo
	3187
	3188	else
	3189
	3190	! Si on est sur un point touche par la convection profonde et pas
	3191	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
	3192	! de la convection profonde.
	3193
	3194	!IM/FH: 2011/02/23
	3195	! definition des points sur lesquels ls thermiques sont actifs
	3196
	3197	do k=1,klev
	3198	do i=1,klon
	3199	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
	3200	cldfra(i,k)=rnebcon(i,k)
	3201	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
	3202	endif
	3203	enddo
	3204	enddo
	3205
	3206	endif
	3207
	3208	else
	3209
	3210	! Ancienne version
	3211	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	3212	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
	3213	endif
	3214
	3215	ENDIF
	3216
	3217	! plulsc(:)=0.
	3218	! do k=1,klev,-1
	3219	! do i=1,klon
	3220	! zzz=prfl(:,k)+psfl(:,k)
	3221	! if (.not.ptconvth.zzz.gt.0.)
	3222	! enddo prfl, psfl,
	3223	! enddo
	3224	!c
	3225	!c 2. NUAGES STARTIFORMES
	3226	!c
	3227	IF (ok_stratus) THEN
	3228	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
	3229	DO k = 1, klev
	3230	DO i = 1, klon
	3231	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	3232	cldliq(i,k) = dialiq(i,k)
	3233	cldfra(i,k) = diafra(i,k)
	3234	ENDIF
	3235	ENDDO
	3236	ENDDO
	3237	ENDIF
	3238	!c
	3239	!c Precipitation totale
	3240	!c
	3241	DO i = 1, klon
	3242	rain_fall(i) = rain_con(i) + rain_lsc(i)
	3243	snow_fall(i) = snow_con(i) + snow_lsc(i)
	3244	ENDDO
	3245	!IM
	3246	IF (ip_ebil_phy.ge.2) THEN
	3247	ztit="after diagcld"
	3248	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	3249	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3250	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3251	call diagphy(airephy,ztit,ip_ebil_phy &
	3252	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	3253	& , zero_v, zero_v, zero_v, ztsol &
	3254	& , d_h_vcol, d_qt, d_ec &
	3255	& , fs_bound, fq_bound )
	3256	END IF
[142]	3257	PRINT *,' Lluis 8 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	3258	!c
	3259	!c Calculer l'humidite relative pour diagnostique
	3260	!c
	3261	DO k = 1, klev
	3262	DO i = 1, klon
	3263	zx_t = t_seri(i,k)
	3264	IF (thermcep) THEN
	3265	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	3266	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	3267	zx_qs = MIN(0.5,zx_qs)
	3268	zcor = 1./(1.-retv*zx_qs)
	3269	zx_qs = zx_qs*zcor
	3270	ELSE
	3271	IF (zx_t.LT.t_coup) THEN
	3272	zx_qs = qsats(zx_t)/pplay(i,k)
	3273	ELSE
	3274	zx_qs = qsatl(zx_t)/pplay(i,k)
	3275	ENDIF
	3276	ENDIF
	3277	zx_rh(i,k) = q_seri(i,k)/zx_qs
	3278	zqsat(i,k)=zx_qs
	3279	ENDDO
	3280	ENDDO
	3281
	3282	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
	3283	!c equivalente a 2m (tpote) pour diagnostique
	3284	!c
	3285	DO i = 1, klon
	3286	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
	3287	IF (thermcep) THEN
	3288	IF(zt2m(i).LT.RTT) then
	3289	Lheat=RLSTT
	3290	ELSE
	3291	Lheat=RLVTT
	3292	ENDIF
	3293	ELSE
	3294	IF (zt2m(i).LT.RTT) THEN
	3295	Lheat=RLSTT
	3296	ELSE
	3297	Lheat=RLVTT
	3298	ENDIF
	3299	ENDIF
	3300	tpote(i) = tpot(i)* &
	3301	& EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
	3302	ENDDO
	3303
	3304	IF (type_trac == 'inca') THEN
	3305	#ifdef INCA
	3306	CALL VTe(VTphysiq)
	3307	CALL VTb(VTinca)
	3308	calday = REAL(days_elapsed + 1) + jH_cur
	3309
	3310	call chemtime(itap+itau_phy-1, date0, dtime)
	3311	IF (config_inca == 'aero') THEN
	3312	CALL AEROSOL_METEO_CALC( &
	3313	& calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
	3314	& prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
	3315	END IF
	3316
	3317	zxsnow_dummy(:) = 0.0
	3318
	3319	CALL chemhook_begin (calday, &
	3320	& days_elapsed+1, &
	3321	& jH_cur, &
	3322	& pctsrf(1,1), &
	3323	& rlat, &
	3324	& rlon, &
	3325	& airephy, &
	3326	& paprs, &
	3327	& pplay, &
	3328	& coefh(:,:,is_ave), &
	3329	& pphi, &
	3330	& t_seri, &
	3331	& u, &
	3332	& v, &
	3333	& wo(:, :, 1), &
	3334	& q_seri, &
	3335	& zxtsol, &
	3336	& zxsnow_dummy, &
	3337	& solsw, &
	3338	& albsol1, &
	3339	& rain_fall, &
	3340	& snow_fall, &
	3341	& itop_con, &
	3342	& ibas_con, &
	3343	& cldfra, &
	3344	& iim, &
	3345	& jjm, &
	3346	& tr_seri, &
	3347	& ftsol, &
	3348	& paprs, &
	3349	& cdragh, &
	3350	& cdragm, &
	3351	& pctsrf, &
	3352	& pdtphys, &
	3353	& itap)
	3354
	3355	CALL VTe(VTinca)
	3356	CALL VTb(VTphysiq)
	3357	#endif
	3358	END IF !type_trac = inca
	3359	!c
	3360	!c Calculer les parametres optiques des nuages et quelques
	3361	!c parametres pour diagnostiques:
	3362	!c
	3363
	3364	IF (aerosol_couple) THEN
	3365	mass_solu_aero(:,:) = ccm(:,:,1)
	3366	mass_solu_aero_pi(:,:) = ccm(:,:,2)
	3367	END IF
	3368
	3369	if (ok_newmicro) then
	3370	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
	3371	& paprs, pplay, t_seri, cldliq, cldfra, &
	3372	& cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
	3373	& flwp, fiwp, flwc, fiwc, &
	3374	& mass_solu_aero, mass_solu_aero_pi, &
	3375	& cldtaupi, re, fl, ref_liq, ref_ice)
	3376	else
	3377	CALL nuage (paprs, pplay, &
	3378	& t_seri, cldliq, cldfra, cldtau, cldemi, &
	3379	& cldh, cldl, cldm, cldt, cldq, &
	3380	& ok_aie, &
	3381	& mass_solu_aero, mass_solu_aero_pi, &
	3382	& bl95_b0, bl95_b1, &
	3383	& cldtaupi, re, fl)
	3384	endif
	3385	!c
	3386	!IM betaCRF
	3387	!c
	3388	cldtaurad = cldtau
	3389	cldtaupirad = cldtaupi
	3390	cldemirad = cldemi &
	3391	&!c
	3392	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
	3393	&lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
	3394	!c
	3395	!c global
	3396	!c
	3397	DO k=1, klev
	3398	DO i=1, klon
	3399	if (pplay(i,k).GE.pfree) THEN
	3400	beta(i,k) = beta_pbl
	3401	else
	3402	beta(i,k) = beta_free
	3403	endif
	3404	if (mskocean_beta) THEN
	3405	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
	3406	endif
	3407	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
	3408	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
	3409	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
	3410	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
	3411	ENDDO
	3412	ENDDO
	3413	!c
	3414	else
	3415	!c
	3416	!c regional
	3417	!c
	3418	DO k=1, klev
	3419	DO i=1,klon
	3420	!c
	3421	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
	3422	& rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
	3423	if (pplay(i,k).GE.pfree) THEN
	3424	beta(i,k) = beta_pbl
	3425	else
	3426	beta(i,k) = beta_free
	3427	endif
	3428	if (mskocean_beta) THEN
	3429	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
	3430	endif
	3431	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
	3432	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
	3433	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
	3434	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
	3435	endif
	3436	!c
	3437	ENDDO
	3438	ENDDO
	3439	!c
	3440	endif
	3441	!c
	3442	!c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
	3443	!c
	3444	IF (MOD(itaprad,radpas).EQ.0) THEN
	3445
	3446	DO i = 1, klon
	3447	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce) &
	3448	& + falb1(i,is_lic) * pctsrf(i,is_lic) &
	3449	& + falb1(i,is_ter) * pctsrf(i,is_ter) &
	3450	& + falb1(i,is_sic) * pctsrf(i,is_sic)
	3451	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce) &
	3452	& + falb2(i,is_lic) * pctsrf(i,is_lic) &
	3453	& + falb2(i,is_ter) * pctsrf(i,is_ter) &
	3454	& + falb2(i,is_sic) * pctsrf(i,is_sic)
	3455	ENDDO
	3456
	3457	if (mydebug) then
	3458	call writefield_phy('u_seri',u_seri,llm)
	3459	call writefield_phy('v_seri',v_seri,llm)
	3460	call writefield_phy('t_seri',t_seri,llm)
	3461	call writefield_phy('q_seri',q_seri,llm)
	3462	endif
	3463
	3464	IF (aerosol_couple) THEN
	3465	#ifdef INCA
	3466	CALL radlwsw_inca &
	3467	& (kdlon,kflev,dist, rmu0, fract, solaire, &
	3468	& paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
	3469	& wo(:, :, 1), &
	3470	& cldfrarad, cldemirad, cldtaurad, &
	3471	& heat,heat0,cool,cool0,radsol,albpla, &
	3472	& topsw,toplw,solsw,sollw, &
	3473	& sollwdown, &
	3474	& topsw0,toplw0,solsw0,sollw0, &
	3475	& lwdn0, lwdn, lwup0, lwup, &
	3476	& swdn0, swdn, swup0, swup, &
	3477	& ok_ade, ok_aie, &
	3478	& tau_aero, piz_aero, cg_aero, &
	3479	& topswad_aero, solswad_aero, &
	3480	& topswad0_aero, solswad0_aero, &
	3481	& topsw_aero, topsw0_aero, &
	3482	& solsw_aero, solsw0_aero, &
	3483	& cldtaupirad, &
	3484	& topswai_aero, solswai_aero)
	3485
	3486	#endif
	3487	ELSE
	3488	!c
	3489	!IM calcul radiatif pour le cas actuel
	3490	!c
	3491	RCO2 = RCO2_act
	3492	RCH4 = RCH4_act
	3493	RN2O = RN2O_act
	3494	RCFC11 = RCFC11_act
	3495	RCFC12 = RCFC12_act
	3496	!c
	3497	IF (prt_level .GE.10) THEN
	3498	print *,' ->radlwsw, number 1 '
	3499	ENDIF
	3500	!c
	3501	CALL radlwsw &
	3502	& (dist, rmu0, fract, &
	3503	& paprs, pplay,zxtsol,albsol1, albsol2, &
	3504	& t_seri,q_seri,wo, &
	3505	& cldfrarad, cldemirad, cldtaurad, &
	3506	& ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
	3507	& flag_aerosol_strat, &
	3508	& tau_aero, piz_aero, cg_aero, &
	3509	& cldtaupirad,new_aod, &
	3510	& zqsat, flwc, fiwc, &
	3511	& heat,heat0,cool,cool0,radsol,albpla, &
	3512	& topsw,toplw,solsw,sollw, &
	3513	& sollwdown, &
	3514	& topsw0,toplw0,solsw0,sollw0, &
	3515	& lwdn0, lwdn, lwup0, lwup, &
	3516	& swdn0, swdn, swup0, swup, &
	3517	& topswad_aero, solswad_aero, &
	3518	& topswai_aero, solswai_aero, &
	3519	& topswad0_aero, solswad0_aero, &
	3520	& topsw_aero, topsw0_aero, &
	3521	& solsw_aero, solsw0_aero, &
	3522	& topswcf_aero, solswcf_aero)
	3523
	3524	!c
	3525	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
	3526	!IM des taux doit etre different du taux actuel
	3527	!IM Par defaut on a les taux perturbes egaux aux taux actuels
	3528	!c
	3529
	3530	if (ok_4xCO2atm) then
	3531	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
	3532	&RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
	3533	&RCFC12_per.NE.RCFC12_act) THEN
	3534	!c
	3535	RCO2 = RCO2_per
	3536	RCH4 = RCH4_per
	3537	RN2O = RN2O_per
	3538	RCFC11 = RCFC11_per
	3539	RCFC12 = RCFC12_per
	3540	!c
	3541	IF (prt_level .GE.10) THEN
	3542	print *,' ->radlwsw, number 2 '
	3543	ENDIF
	3544	!c
	3545	CALL radlwsw &
	3546	& (dist, rmu0, fract, &
	3547	& paprs, pplay,zxtsol,albsol1, albsol2, &
	3548	& t_seri,q_seri,wo, &
	3549	& cldfra, cldemi, cldtau, &
	3550	& ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
	3551	& flag_aerosol_strat, &
	3552	& tau_aero, piz_aero, cg_aero, &
	3553	& cldtaupi,new_aod, &
	3554	& zqsat, flwc, fiwc, &
	3555	& heatp,heat0p,coolp,cool0p,radsolp,albplap, &
	3556	& topswp,toplwp,solswp,sollwp, &
	3557	& sollwdownp, &
	3558	& topsw0p,toplw0p,solsw0p,sollw0p, &
	3559	& lwdn0p, lwdnp, lwup0p, lwupp, &
	3560	& swdn0p, swdnp, swup0p, swupp, &
	3561	& topswad_aerop, solswad_aerop, &
	3562	& topswai_aerop, solswai_aerop, &
	3563	& topswad0_aerop, solswad0_aerop, &
	3564	& topsw_aerop, topsw0_aerop, &
	3565	& solsw_aerop, solsw0_aerop, &
	3566	& topswcf_aerop, solswcf_aerop)
	3567	endif
	3568	endif
	3569	!c
	3570	ENDIF ! aerosol_couple
	3571	itaprad = 0
	3572	ENDIF ! MOD(itaprad,radpas)
	3573	itaprad = itaprad + 1
	3574
	3575	IF (iflag_radia.eq.0) THEN
	3576	IF (prt_level.ge.9) THEN
	3577	PRINT *,'--------------------------------------------------'
	3578	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
	3579	PRINT *,'>>>> heat et cool mis a zero '
	3580	PRINT *,'--------------------------------------------------'
	3581	END IF
	3582	heat=0.
	3583	cool=0.
	3584	sollw=0. ! MPL 01032011
	3585	solsw=0.
	3586	radsol=0.
	3587	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
	3588	swup0=0.
	3589	swdn=0.
	3590	swdn0=0.
	3591	lwup=0.
	3592	lwup0=0.
	3593	lwdn=0.
	3594	lwdn0=0.
	3595	END IF
	3596
	3597	!c
	3598	!c Ajouter la tendance des rayonnements (tous les pas)
	3599	!c
	3600	DO k = 1, klev
	3601	DO i = 1, klon
	3602	t_seri(i,k) = t_seri(i,k) &
	3603	& + (heat(i,k)-cool(i,k)) * dtime/RDAY
	3604	ENDDO
	3605	ENDDO
	3606	!c
	3607	if (mydebug) then
	3608	call writefield_phy('u_seri',u_seri,llm)
	3609	call writefield_phy('v_seri',v_seri,llm)
	3610	call writefield_phy('t_seri',t_seri,llm)
	3611	call writefield_phy('q_seri',q_seri,llm)
	3612	endif
	3613
	3614	!IM
	3615	IF (ip_ebil_phy.ge.2) THEN
	3616	ztit='after rad'
	3617	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	3618	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3619	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3620	call diagphy(airephy,ztit,ip_ebil_phy &
	3621	& , topsw, toplw, solsw, sollw, zero_v &
	3622	& , zero_v, zero_v, zero_v, ztsol &
	3623	& , d_h_vcol, d_qt, d_ec &
	3624	& , fs_bound, fq_bound )
	3625	END IF
[142]	3626	PRINT *,' Lluis 9 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	3627	!c
	3628	!c
	3629	!c Calculer l'hydrologie de la surface
	3630	!c
	3631	!c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
	3632	!c . agesno, ftsol,fqsurf,fsnow, ruis)
	3633	!c
	3634
	3635	!c
	3636	!c Calculer le bilan du sol et la derive de temperature (couplage)
	3637	!c
	3638	DO i = 1, klon
	3639	!c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
	3640	!c a la demande de JLD
	3641	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
	3642	ENDDO
	3643	!c
	3644	!cmoddeblott(jan95)
	3645	!c Appeler le programme de parametrisation de l'orographie
	3646	!c a l'echelle sous-maille:
	3647	!c
	3648	IF (prt_level .GE.10) THEN
	3649	print *,' call orography ? ', ok_orodr
	3650	ENDIF
	3651	!c
	3652	IF (ok_orodr) THEN
	3653	!c
	3654	!c selection des points pour lesquels le shema est actif:
	3655	igwd=0
	3656	DO i=1,klon
	3657	itest(i)=0
	3658	!c IF ((zstd(i).gt.10.0)) THEN
	3659	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	3660	itest(i)=1
	3661	igwd=igwd+1
	3662	idx(igwd)=i
	3663	ENDIF
	3664	ENDDO
	3665	!c igwdim=MAX(1,igwd)
	3666	!c
	3667	IF (ok_strato) THEN
	3668
	3669	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
	3670	& zmea,zstd, zsig, zgam, zthe,zpic,zval, &
	3671	& igwd,idx,itest, &
	3672	& t_seri, u_seri, v_seri, &
	3673	& zulow, zvlow, zustrdr, zvstrdr, &
	3674	& d_t_oro, d_u_oro, d_v_oro)
	3675
	3676	ELSE
	3677	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
	3678	& zmea,zstd, zsig, zgam, zthe,zpic,zval, &
	3679	& igwd,idx,itest, &
	3680	& t_seri, u_seri, v_seri, &
	3681	& zulow, zvlow, zustrdr, zvstrdr, &
	3682	& d_t_oro, d_u_oro, d_v_oro)
	3683	ENDIF
	3684	!c
	3685	!c ajout des tendances
	3686	!-----------------------------------------------------------------------------------------
	3687	! ajout des tendances de la trainee de l'orographie
	3688	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
	3689	!-----------------------------------------------------------------------------------------
	3690	!c
	3691	ENDIF ! fin de test sur ok_orodr
	3692	!c
	3693	if (mydebug) then
	3694	call writefield_phy('u_seri',u_seri,llm)
	3695	call writefield_phy('v_seri',v_seri,llm)
	3696	call writefield_phy('t_seri',t_seri,llm)
	3697	call writefield_phy('q_seri',q_seri,llm)
	3698	endif
	3699
	3700	IF (ok_orolf) THEN
	3701	!c
	3702	!c selection des points pour lesquels le shema est actif:
	3703	igwd=0
	3704	DO i=1,klon
	3705	itest(i)=0
	3706	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	3707	itest(i)=1
	3708	igwd=igwd+1
	3709	idx(igwd)=i
	3710	ENDIF
	3711	ENDDO
	3712	!c igwdim=MAX(1,igwd)
	3713	!c
	3714	IF (ok_strato) THEN
	3715
	3716	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
	3717	& rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
	3718	& igwd,idx,itest, &
	3719	& t_seri, u_seri, v_seri, &
	3720	& zulow, zvlow, zustrli, zvstrli, &
	3721	& d_t_lif, d_u_lif, d_v_lif )
	3722
	3723	ELSE
	3724	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
	3725	& rlat,zmea,zstd,zpic, &
	3726	& itest, &
	3727	& t_seri, u_seri, v_seri, &
	3728	& zulow, zvlow, zustrli, zvstrli, &
	3729	& d_t_lif, d_u_lif, d_v_lif)
	3730	ENDIF
	3731	!c
	3732	!-----------------------------------------------------------------------------------------
	3733	! ajout des tendances de la portance de l'orographie
	3734	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
	3735	!-----------------------------------------------------------------------------------------
	3736	!c
	3737	ENDIF ! fin de test sur ok_orolf
	3738	!C HINES GWD PARAMETRIZATION
	3739
	3740	IF (ok_hines) then
	3741
	3742	CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
	3743	& rlat,t_seri,u_seri,v_seri, &
	3744	& zustrhi,zvstrhi, &
	3745	& d_t_hin, d_u_hin, d_v_hin)
	3746	!c
	3747	!c ajout des tendances
	3748	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
	3749
	3750	ENDIF
	3751	!c
	3752
	3753	!c
	3754	!IM cf. FLott BEG
	3755	!C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
	3756
	3757	if (mydebug) then
	3758	call writefield_phy('u_seri',u_seri,llm)
	3759	call writefield_phy('v_seri',v_seri,llm)
	3760	call writefield_phy('t_seri',t_seri,llm)
	3761	call writefield_phy('q_seri',q_seri,llm)
	3762	endif
	3763
	3764	DO i = 1, klon
	3765	zustrph(i)=0.
	3766	zvstrph(i)=0.
	3767	ENDDO
	3768	DO k = 1, klev
	3769	DO i = 1, klon
	3770	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
	3771	& (paprs(i,k)-paprs(i,k+1))/rg
	3772	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
	3773	& (paprs(i,k)-paprs(i,k+1))/rg
	3774	ENDDO
	3775	ENDDO
	3776	!c
	3777	!IM calcul composantes axiales du moment angulaire et couple des montagnes
	3778	!c
	3779	! IF (is_sequential .and. ok_orodr) THEN
	3780	! CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
	3781	! & ra,rg,romega, &
	3782	! & rlat,rlon,pphis, &
	3783	! & zustrdr,zustrli,zustrph, &
	3784	! & zvstrdr,zvstrli,zvstrph, &
	3785	! & paprs,u,v, &
	3786	! & aam, torsfc)
	3787	! ENDIF
	3788	!IM cf. FLott END
	3789	!IM
	3790	IF (ip_ebil_phy.ge.2) THEN
	3791	ztit='after orography'
	3792	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	3793	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3794	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3795	call diagphy(airephy,ztit,ip_ebil_phy &
	3796	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	3797	& , zero_v, zero_v, zero_v, ztsol &
	3798	& , d_h_vcol, d_qt, d_ec &
	3799	& , fs_bound, fq_bound )
	3800	END IF
[142]	3801	PRINT *,' Lluis 10 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	3802	!c
	3803	!c
	3804	!====================================================================
	3805	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
	3806	!====================================================================
	3807	! Abderrahmane 24.08.09
	3808
	3809	IF (ok_cosp) THEN
	3810	! adeclarer
	3811	#ifdef CPP_COSP
	3812	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
	3813
	3814	print*,'freq_cosp',freq_cosp
	3815	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
	3816	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
	3817	! s ref_liq,ref_ice
	3818	call phys_cosp(itap,dtime,freq_cosp, &
	3819	& ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
	3820	& ecrit_mth,ecrit_day,ecrit_hf, &
	3821	& klon,klev,rlon,rlat,presnivs,overlap, &
	3822	& ref_liq,ref_ice, &
	3823	& pctsrf(:,is_ter)+pctsrf(:,is_lic), &
	3824	& zu10m,zv10m,pphis, &
	3825	& zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
	3826	& qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
	3827	& prfl(:,1:klev),psfl(:,1:klev), &
	3828	& pmflxr(:,1:klev),pmflxs(:,1:klev), &
	3829	& mr_ozone,cldtau, cldemi)
	3830
	3831	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
	3832	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
	3833	! M clMISR,
	3834	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
	3835	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
	3836
	3837	ENDIF
	3838
	3839	#endif
	3840	ENDIF !ok_cosp
	3841	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3842	!cAA
	3843	!cAA Installation de l'interface online-offline pour traceurs
	3844	!cAA
	3845	!c====================================================================
	3846	!c Calcul des tendances traceurs
	3847	!c====================================================================
	3848	!C
	3849
	3850	IF (type_trac=='repr') THEN
	3851	sh_in(:,:) = q_seri(:,:)
	3852	ELSE
	3853	sh_in(:,:) = qx(:,:,ivap)
	3854	END IF
	3855
	3856	call phytrac ( &
	3857	& itap, days_elapsed+1, jH_cur, debut, &
	3858	& lafin, dtime, u, v, t, &
	3859	& paprs, pplay, pmfu, pmfd, &
	3860	& pen_u, pde_u, pen_d, pde_d, &
	3861	& cdragh, coefh(:,:,is_ave), fm_therm, entr_therm, &
	3862	& u1, v1, ftsol, pctsrf, &
	3863	& ustar, u10m, v10m, &
	3864	& rlat, rlon, &
	3865	& frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
	3866	& presnivs, pphis, pphi, albsol1, &
	3867	& sh_in, rhcl, cldfra, rneb, &
	3868	& diafra, cldliq, itop_con, ibas_con, &
	3869	& pmflxr, pmflxs, prfl, psfl, &
	3870	& da, phi, mp, upwd, &
	3871	& phi2, d1a, dam, sij, &
	3872	& wdtrainA, wdtrainM, sigd, clw,elij, &
	3873	& ev, ep, epmlmMm, eplaMm, &
	3874	& dnwd, aerosol_couple, flxmass_w, &
	3875	& tau_aero, piz_aero, cg_aero, ccm, &
	3876	& rfname, &
	3877	& d_tr_dyn, &
	3878	& tr_seri)
	3879
	3880	IF (offline) THEN
	3881
	3882	IF (prt_level.ge.9) &
	3883	& print*,'Attention on met a 0 les thermiques pour phystoke'
	3884	call phystokenc ( &
	3885	& nlon,klev,pdtphys,rlon,rlat, &
	3886	& t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
	3887	& fm_therm,entr_therm, &
	3888	& cdragh,coefh(:,:,is_ave),u1,v1,ftsol,pctsrf, &
	3889	& frac_impa, frac_nucl, &
	3890	& pphis,airephy,dtime,itap, &
	3891	& qx(:,:,ivap),da,phi,mp,upwd,dnwd)
	3892
	3893
	3894	ENDIF
	3895
	3896	!c
	3897	!c Calculer le transport de l'eau et de l'energie (diagnostique)
	3898	!c
	3899	CALL transp (paprs,zxtsol, &
	3900	& t_seri, q_seri, u_seri, v_seri, zphi, &
	3901	& ve, vq, ue, uq)
	3902	!c
	3903	!IM global posePB BEG
	3904	IF(1.EQ.0) THEN
	3905	!c
	3906	CALL transp_lay (paprs,zxtsol, &
	3907	& t_seri, q_seri, u_seri, v_seri, zphi, &
	3908	& ve_lay, vq_lay, ue_lay, uq_lay)
	3909	!c
	3910	ENDIF !(1.EQ.0) THEN
	3911	!IM global posePB END
	3912	!c Accumuler les variables a stocker dans les fichiers histoire:
	3913	!c
	3914
	3915	!================================================================
	3916	! Conversion of kinetic and potential energy into heat, for
	3917	! parameterisation of subgrid-scale motions
	3918	!================================================================
	3919
	3920	d_t_ec(:,:)=0.
	3921	forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
	3922	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
	3923	& u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
	3924	& zmasse,exner,d_t_ec)
	3925	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
	3926
	3927	!IM
	3928	IF (ip_ebil_phy.ge.1) THEN
	3929	ztit='after physic'
	3930	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
	3931	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3932	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3933	!C Comme les tendances de la physique sont ajoute dans la dynamique,
	3934	!C on devrait avoir que la variation d'entalpie par la dynamique
	3935	!C est egale a la variation de la physique au pas de temps precedent.
	3936	!C Donc la somme de ces 2 variations devrait etre nulle.
	3937
[142]	3938	PRINT *,' LLuis sending: d_h_vcol: ',d_h_vcol,' d_qt ',d_qt,' d_ec ',d_ec
	3939
[1]	3940	call diagphy(airephy,ztit,ip_ebil_phy &
	3941	& , topsw, toplw, solsw, sollw, sens &
	3942	& , evap, rain_fall, snow_fall, ztsol &
	3943	& , d_h_vcol, d_qt, d_ec &
	3944	& , fs_bound, fq_bound )
	3945	!C
	3946	d_h_vcol_phy=d_h_vcol
	3947	!C
	3948	END IF
[142]	3949	PRINT *,' Lluis 11 d_h_vcol: ',d_h_vcol,' d_h_vcol_phy: ',d_h_vcol_phy
[1]	3950
	3951	!C
	3952	!c=======================================================================
	3953	!c SORTIES
	3954	!c=======================================================================
	3955
	3956	!IM Interpolation sur les niveaux de pression du NMC
	3957	!c -------------------------------------------------
	3958	!c
	3959	#include "calcul_STDlev.h"
	3960	!c
	3961	!c slp sea level pressure
	3962	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
	3963	!c
	3964	!ccc prw = eau precipitable
	3965	DO i = 1, klon
	3966	prw(i) = 0.
	3967	DO k = 1, klev
	3968	prw(i) = prw(i) + &
	3969	& q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
	3970	ENDDO
	3971	ENDDO
	3972	!c
	3973	!IM initialisation + calculs divers diag AMIP2
	3974	!c
	3975	#include "calcul_divers.h"
	3976	!c
	3977	IF (type_trac == 'inca') THEN
	3978	#ifdef INCA
	3979	CALL VTe(VTphysiq)
	3980	CALL VTb(VTinca)
	3981
	3982	CALL chemhook_end ( &
	3983	& dtime, &
	3984	& pplay, &
	3985	& t_seri, &
	3986	& tr_seri, &
	3987	& nbtr, &
	3988	& paprs, &
	3989	& q_seri, &
	3990	& airephy, &
	3991	& pphi, &
	3992	& pphis, &
	3993	& zx_rh)
	3994
	3995	CALL VTe(VTinca)
	3996	CALL VTb(VTphysiq)
	3997	#endif
	3998	END IF
	3999
	4000
	4001	!c
	4002	!c Convertir les incrementations en tendances
	4003	!c
	4004	IF (prt_level .GE.10) THEN
	4005	print *,'Convertir les incrementations en tendances '
	4006	ENDIF
	4007	!c
	4008	if (mydebug) then
	4009	call writefield_phy('u_seri',u_seri,llm)
	4010	call writefield_phy('v_seri',v_seri,llm)
	4011	call writefield_phy('t_seri',t_seri,llm)
	4012	call writefield_phy('q_seri',q_seri,llm)
	4013	endif
	4014
	4015	DO k = 1, klev
	4016	DO i = 1, klon
	4017	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	4018	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	4019	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
	4020	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
	4021	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
	4022	ENDDO
	4023	ENDDO
	4024	!c
	4025	IF (nqtot.GE.3) THEN
	4026	DO iq = 3, nqtot
	4027	DO k = 1, klev
	4028	DO i = 1, klon
	4029	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
	4030	ENDDO
	4031	ENDDO
	4032	ENDDO
	4033	ENDIF
	4034	!c
	4035	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
	4036	!IM global posePB#include "write_bilKP_ins.h"
	4037	!IM global posePB#include "write_bilKP_ave.h"
	4038	!c
	4039
	4040	!c Sauvegarder les valeurs de t et q a la fin de la physique:
	4041	!c
	4042	DO k = 1, klev
	4043	DO i = 1, klon
	4044	u_ancien(i,k) = u_seri(i,k)
	4045	v_ancien(i,k) = v_seri(i,k)
	4046	t_ancien(i,k) = t_seri(i,k)
	4047	q_ancien(i,k) = q_seri(i,k)
	4048	ENDDO
	4049	ENDDO
	4050
	4051	!!! RomP >>>
	4052	IF (nqtot.GE.3) THEN
	4053	DO iq = 3, nqtot
	4054	DO k = 1, klev
	4055	DO i = 1, klon
	4056	tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
	4057	ENDDO
	4058	ENDDO
	4059	ENDDO
	4060	ENDIF
	4061	!!! RomP <<<
	4062	!==========================================================================
	4063	! Sorties des tendances pour un point particulier
	4064	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
	4065	! pour le debug
	4066	! La valeur de igout est attribuee plus haut dans le programme
	4067	!==========================================================================
	4068
	4069	if (prt_level.ge.1) then
	4070	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	4071	write(lunout,*) &
	4072	& 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
	4073	write(lunout,*) &
	4074	& nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
	4075	& pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
	4076	& pctsrf(igout,is_sic)
	4077	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
	4078	do k=1,klev
	4079	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
	4080	& d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
	4081	& d_t_eva(igout,k)
	4082	enddo
	4083	write(lunout,*) 'cool,heat'
	4084	do k=1,klev
	4085	write(lunout,*) cool(igout,k),heat(igout,k)
	4086	enddo
	4087
	4088	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
	4089	do k=1,klev
	4090	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
	4091	& d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
	4092	enddo
	4093
	4094	write(lunout,*) 'd_ps ',d_ps(igout)
	4095	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
	4096	do k=1,klev,klon
	4097	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
	4098	& d_qx(igout,k,1),d_qx(igout,k,2)
	4099	enddo
	4100	endif
	4101
	4102	!==========================================================================
	4103
	4104	!c============================================================
	4105	!c Calcul de la temperature potentielle
	4106	!c============================================================
	4107	DO k = 1, klev
	4108	DO i = 1, klon
	4109	!cJYG/IM theta en debut du pas de temps
	4110	!cJYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
	4111	!cJYG/IM theta en fin de pas de temps de physique
	4112	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
	4113	!c thetal: 2 lignes suivantes a decommenter si vous avez les fichiers MPL 20130625
	4114	!c fth_fonctions.F90 et parkind1.F90
	4115	!c sinon thetal=theta
	4116	!c thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
	4117	!c : ql_seri(i,k))
	4118	thetal(i,k)=theta(i,k)
	4119	ENDDO
	4120	ENDDO
	4121	!c
	4122
	4123	!c 22.03.04 BEG
	4124	!c=============================================================
	4125	!c Ecriture des sorties
	4126	!c=============================================================
	4127	#ifdef CPP_IOIPSL
	4128
	4129	!c Recupere des varibles calcule dans differents modules
	4130	!c pour ecriture dans histxxx.nc
	4131
	4132	! Get some variables from module fonte_neige_mod
	4133	!L. Fita, LMD. November 2013
	4134	!! It is not working after removing starphy.nc and limit.nc?
	4135	CALL fonte_neige_get_vars(pctsrf, &
	4136	& zxfqcalving, zxfqfonte, zxffonte)
	4137
	4138
	4139
	4140	!c=============================================================
	4141	! Separation entre thermiques et non thermiques dans les sorties
	4142	! de fisrtilp
	4143	!c=============================================================
	4144
	4145	if (iflag_thermals>=1) then
	4146	d_t_lscth=0.
	4147	d_t_lscst=0.
	4148	d_q_lscth=0.
	4149	d_q_lscst=0.
	4150	do k=1,klev
	4151	do i=1,klon
	4152	if (ptconvth(i,k)) then
	4153	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
	4154	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
	4155	else
	4156	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
	4157	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
	4158	endif
	4159	enddo
	4160	enddo
	4161
	4162	do i=1,klon
	4163	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
	4164	plul_th(i)=prfl(i,1)+psfl(i,1)
	4165	enddo
	4166	endif
	4167
	4168	#include "phys_output_write_new.h"
	4169
	4170
	4171	#ifdef histISCCP
	4172	#include "write_histISCCP.h"
	4173	#endif
	4174
	4175	#ifdef histNMC
	4176	#include "write_histhfNMC.h"
	4177	#include "write_histdayNMC.h"
	4178	#include "write_histmthNMC.h"
	4179	#endif
	4180
	4181	#include "write_histday_seri.h"
	4182
	4183	#include "write_paramLMDZ_phy.h"
	4184
	4185	#endif
	4186
	4187	!c 22.03.04 END
	4188	!c
	4189	!c====================================================================
	4190	!c Si c'est la fin, il faut conserver l'etat de redemarrage
	4191	!c====================================================================
	4192	!c
	4193
	4194	!c -----------------------------------------------------------------
	4195	!c WSTATS: Saving statistics
	4196	!c -----------------------------------------------------------------
	4197	!c ("stats" stores and accumulates 8 key variables in file "stats.nc"
	4198	!c which can later be used to make the statistic files of the run:
	4199	!c "stats") only possible in 3D runs !
	4200
	4201
	4202	IF (callstats) THEN
	4203
	4204	call wstats(klon,o_psol%name,"Surface pressure","Pa" &
	4205	& ,2,paprs(:,1))
	4206	call wstats(klon,o_tsol%name,"Surface temperature","K", &
	4207	& 2,zxtsol)
	4208	zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
	4209	call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
	4210	& "kg/(s*m2)",2,zx_tmp_fi2d)
	4211	zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
	4212	call wstats(klon,o_plul%name,"Large-scale Precip", &
	4213	& "kg/(s*m2)",2,zx_tmp_fi2d)
	4214	zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
	4215	call wstats(klon,o_pluc%name,"Convective Precip", &
	4216	& "kg/(s*m2)",2,zx_tmp_fi2d)
	4217	call wstats(klon,o_sols%name,"Solar rad. at surf.", &
	4218	& "W/m2",2,solsw)
	4219	call wstats(klon,o_soll%name,"IR rad. at surf.", &
	4220	& "W/m2",2,sollw)
	4221	zx_tmp_fi2d(:) = topsw(:)-toplw(:)
	4222	call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
	4223	& "W/m2",2,zx_tmp_fi2d)
	4224
	4225
	4226
	4227	call wstats(klon,o_temp%name,"Air temperature","K", &
	4228	& 3,t_seri)
	4229	call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
	4230	& 3,u_seri)
	4231	call wstats(klon,o_vitv%name,"Meridional wind", &
	4232	& "m.s-1",3,v_seri)
	4233	call wstats(klon,o_vitw%name,"Vertical wind", &
	4234	& "m.s-1",3,omega)
	4235	call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
	4236	& 3,q_seri)
	4237
	4238
	4239
	4240	IF(lafin) THEN
	4241	write (,) "Writing stats..."
	4242	call mkstats(ierr)
	4243	ENDIF
	4244
	4245	ENDIF !if callstats
	4246
	4247	IF (lafin) THEN
	4248	itau_phy = itau_phy + itap
	4249	CALL phyredem ("restartphy.nc")
	4250	! open(97,form="unformatted",file="finbin")
	4251	! write(97) u_seri,v_seri,t_seri,q_seri
	4252	! close(97)
	4253	!$OMP MASTER
	4254	if (read_climoz >= 1) then
	4255	if (is_mpi_root) then
	4256	call nf95_close(ncid_climoz)
	4257	end if
	4258	deallocate(press_climoz) ! pointer
	4259	end if
	4260	!$OMP END MASTER
	4261	ENDIF
	4262
	4263	! first=.false.
	4264
	4265	RETURN
	4266	END SUBROUTINE physiq
	4267
	4268	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
	4269	IMPLICIT none
	4270	!c
	4271	!c Calculer et imprimer l'eau totale. A utiliser pour verifier
	4272	!c la conservation de l'eau
	4273	!c
	4274	#include "YOMCST.h"
	4275	INTEGER klon,klev
	4276	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
	4277	REAL aire(klon)
	4278	REAL qtotal, zx, qcheck
	4279	INTEGER i, k
	4280	!c
	4281	zx = 0.0
	4282	DO i = 1, klon
	4283	zx = zx + aire(i)
	4284	ENDDO
	4285	qtotal = 0.0
	4286	DO k = 1, klev
	4287	DO i = 1, klon
	4288	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
	4289	& *(paprs(i,k)-paprs(i,k+1))/RG
	4290	ENDDO
	4291	ENDDO
	4292	!c
	4293	qcheck = qtotal/zx
	4294	!c
	4295	RETURN
	4296	END
	4297
	4298	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	4299	IMPLICIT none
	4300	!c
	4301	!c Tranformer une variable de la grille physique a
	4302	!c la grille d'ecriture
	4303	!c
	4304	INTEGER nfield,nlon,iim,jjmp1, jjm
	4305	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
	4306	!c
	4307	INTEGER i, n, ig
	4308	!c
	4309	jjm = jjmp1 - 1
	4310	DO n = 1, nfield
	4311	DO i=1,iim
	4312	ecrit(i,n) = fi(1,n)
	4313	ecrit(i+jjm*iim,n) = fi(nlon,n)
	4314	ENDDO
	4315	DO ig = 1, nlon - 2
	4316	ecrit(iim+ig,n) = fi(1+ig,n)
	4317	ENDDO
	4318	ENDDO
	4319	RETURN
	4320	END SUBROUTINE gr_fi_ecrit

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

Download in other formats:

Original Format