Context Navigation

physiq.F90 @ 6

Last change on this file since 6 was 6, checked in by lfita, 10 years ago

Fixing problem with the solar radiation:

1.- Correct values for the julian and the hour (0,1) to the LMDZ physiq

Adding NaN checking for 'PSFC' and 'T_2'

File size: 166.0 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
	1301	! Lluis
	1302	INTEGER :: llp
	1303	llp = 644
	1304
	1305	!c======================================================================
	1306	! Gestion calendrier : mise a jour du module phys_cal_mod
	1307	!
	1308	CALL phys_cal_update(jD_cur,jH_cur)
[6]	1309	PRINT *,' Lluis physiq: jD_cur: ',jD_cur, ' jH_cur: ',jH_cur, &
	1310	' days_elapsed: ',days_elapsed
[1]	1311
	1312	!c======================================================================
	1313	! Ecriture eventuelle d'un profil verticale en entree de la physique.
	1314	! Utilise notamment en 1D mais peut etre active egalement en 3D
	1315	! en imposant la valeur de igout.
	1316	!c======================================================================
	1317
	1318	if (prt_level.ge.1) then
	1319	! Lluis
	1320	! igout=klon/2+1/klon
	1321	igout=llp
	1322	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	1323	write(lunout,*) &
	1324	& 'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
	1325	write(lunout,*) &
	1326	& nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
	1327
	1328	write(lunout,*) 'paprs, play, phi, u, v, t'
	1329	do k=1,klev
	1330	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
	1331	& u(igout,k),v(igout,k),t(igout,k)
	1332	enddo
	1333	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
	1334	do k=1,klev
	1335	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
	1336	enddo
	1337	endif
	1338
	1339	!c======================================================================
	1340
	1341	!cym => necessaire pour iflag_con != 2
	1342	pmfd(:,:) = 0.
	1343	pen_u(:,:) = 0.
	1344	pen_d(:,:) = 0.
	1345	pde_d(:,:) = 0.
	1346	pde_u(:,:) = 0.
	1347	aam=0.
	1348
	1349	torsfc=0.
	1350
	1351	forall (k=1: llm) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
	1352
	1353	if (first) then
	1354
	1355	!cCR:nvelles variables convection/poches froides
	1356
	1357	print*, '================================================='
	1358	print*, 'Allocation des variables locales et sauvegardees'
	1359
	1360	call phys_local_var_init
	1361
	1362	!c
	1363	pasphys=pdtphys
	1364	!c appel a la lecture du run.def physique
	1365	call conf_phys(ok_journe, ok_mensuel, &
	1366	& ok_instan, ok_hf, &
	1367	& ok_LES, &
	1368	& callstats, &
	1369	& solarlong0,seuil_inversion, &
	1370	& fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
	1371	& iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
	1372	& ok_ade, ok_aie, ok_cdnc, aerosol_couple, &
	1373	& flag_aerosol, flag_aerosol_strat, new_aod, &
	1374	& bl95_b0, bl95_b1, &
	1375	!c nv flags pour la convection et les poches froides &
	1376	& read_climoz, &
	1377	& alp_offset)
	1378
	1379
	1380	!! call phys_state_var_init(read_climoz)
	1381	call phys_output_var_init
	1382
	1383
	1384	!L. Fita, LMD. November 2013
	1385	! Initializing
	1386	! Variables have been already initialilzed from the WRF variables.
	1387	! Using on purpose subroutines from physiq_limit_variables_mod
	1388	!! CALL neige_initialize()
	1389	!! CALL limit_initialize()
	1390	!! CALL vars_limit_init(klon)
[6]	1391
[1]	1392	! L. Fita, LMD January 2014. Initializing output variables
	1393	CALL init_ovars_lmdz_NOmodule(klon,klev,nbsrf)
	1394
	1395	! Allocating restart variables
	1396	! ALLOCATE(qsol(klon))
	1397	! ALLOCATE(fder(klon))
	1398	! ALLOCATE(snow(klon, nbsrf))
	1399	! ALLOCATE(qsurf(klon, nbsrf))
	1400	! ALLOCATE(evap(klon, nbsrf))
	1401	! ALLOCATE(rugos(klon, nbsrf))
	1402	! ALLOCATE(agesno(klon, nbsrf))
	1403	! ALLOCATE(ftsoil(klon, nsoilmx, nbsrf))
	1404	!! ALLOCATE(restart_runoff(klon))
	1405
	1406	! CALL pbl_surface_init(qsol_rst, fder_rst, snow_rst, qsurf_rst, evap_rst, &
	1407	! rugos_rst, agesno_rst, ftsoil_rst)
	1408	! CALL fonte_neige_init(restart_runoff)
	1409	! qsol=qsol_rst
	1410
	1411	print*, '================================================='
	1412	!c
	1413	dnwd0=0.0
	1414	ftd=0.0
	1415	fqd=0.0
	1416	cin=0.
	1417	!cym Attention pbase pas initialise dans concvl !!!!
	1418	pbase=0
	1419	!IM 180608
	1420
	1421	itau_con=0
	1422	first=.false.
	1423
	1424	endif ! first
	1425
	1426	modname = 'physiq'
	1427	!IM
	1428	IF (ip_ebil_phy.ge.1) THEN
	1429	DO i=1,klon
	1430	zero_v(i)=0.
	1431	END DO
	1432	END IF
	1433	ok_sync=.TRUE.
	1434
	1435	IF (debut) THEN
	1436	CALL suphel ! initialiser constantes et parametres phys.
	1437	ENDIF
	1438
	1439	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
	1440
	1441
	1442	!c======================================================================
	1443	! Gestion calendrier : mise a jour du module phys_cal_mod
	1444	!
	1445	!c CALL phys_cal_update(jD_cur,jH_cur)
	1446
	1447	!c
	1448	!c Si c'est le debut, il faut initialiser plusieurs choses
	1449	!c ********
	1450	!c
	1451
	1452	IF (debut) THEN
	1453	!rv
	1454	!cCRinitialisation de wght_th et lalim_conv pour la definition de la couche alimentation
	1455	!cde la convection a partir des caracteristiques du thermique
	1456	wght_th(:,:)=1.
	1457	lalim_conv(:)=1
	1458	!cRC
	1459	ustar(:,:)=0.
	1460	u10m(:,:)=0.
	1461	v10m(:,:)=0.
	1462	rain_con(:)=0.
	1463	snow_con(:)=0.
	1464	topswai(:)=0.
	1465	topswad(:)=0.
	1466	solswai(:)=0.
	1467	solswad(:)=0.
	1468
	1469	wmax_th(:)=0.
	1470	tau_overturning_th(:)=0.
	1471
	1472	IF (type_trac == 'inca') THEN
	1473	! jg : initialisation jusqu'au ces variables sont dans restart
	1474	ccm(:,:,:) = 0.
	1475	tau_aero(:,:,:,:) = 0.
	1476	piz_aero(:,:,:,:) = 0.
	1477	cg_aero(:,:,:,:) = 0.
	1478	END IF
	1479
	1480	rnebcon0(:,:) = 0.0
	1481	clwcon0(:,:) = 0.0
	1482	rnebcon(:,:) = 0.0
	1483	clwcon(:,:) = 0.0
	1484
	1485	!IM
	1486	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
	1487	!c
	1488	print*,'iflag_coupl,iflag_clos,iflag_wake', &
	1489	& iflag_coupl,iflag_clos,iflag_wake
	1490	print*,'CYCLE_DIURNE', cycle_diurne
	1491	!c
	1492	IF (iflag_con.EQ.2.AND.iflag_cldcon.GT.-1) THEN
	1493	abort_message = 'Tiedtke needs iflag_cldcon=-2 or -1'
	1494	CALL abort_gcm (modname,abort_message,1)
	1495	ENDIF
	1496	!c
	1497	IF(ok_isccp.AND.iflag_con.LE.2) THEN
	1498	abort_message = 'ISCCP-like outputs may be available for KE &
	1499	&(iflag_con >= 3); for Tiedtke (iflag_con=-2) put ok_isccp=n'
	1500	CALL abort_gcm (modname,abort_message,1)
	1501	ENDIF
	1502	!c
	1503	!c Initialiser les compteurs:
	1504	!c
	1505	itap = 0
	1506	itaprad = 0
	1507
	1508	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1509	!! Un petit travail \`a faire ici.
	1510	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1511
	1512	if (iflag_pbl>1) then
	1513	PRINT*, "Using method MELLOR&YAMADA"
	1514	endif
	1515
	1516	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1517	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
	1518	! dyn3d
	1519	! Attention : la version precedente n'etait pas tres propre.
	1520	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
	1521	! pour obtenir le meme resultat.
	1522	dtime=pdtphys
	1523	radpas = NINT( 86400./dtime/nbapp_rad)
	1524	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1525
	1526	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
	1527
	1528	IF (klon_glo==1) THEN
	1529	coefh=0. ; coefm=0. ; pbl_tke=0.
	1530	coefh(:,2,:)=1.e-2 ; coefm(:,2,:)=1.e-2 ; pbl_tke(:,2,:)=1.e-2
	1531	PRINT*,'FH WARNING : lignes a supprimer'
	1532	ENDIF
	1533	!IM begin
	1534	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
	1535	&,ratqs(1,1)
	1536	!IM end
	1537
	1538
	1539
	1540	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1541	!c
	1542	!C on remet le calendrier a zero
	1543	!c
	1544	IF (raz_date .eq. 1) THEN
	1545	itau_phy = 0
	1546	ENDIF
	1547
	1548	!IM cf. AM 081204 BEG
	1549	PRINT*,'cycle_diurne3 =',cycle_diurne
	1550	!IM cf. AM 081204 END
	1551	!c
	1552	CALL printflag( tabcntr0,radpas,ok_journe, &
	1553	& ok_instan, ok_region )
	1554	!c
	1555	IF (ABS(dtime-pdtphys).GT.0.001) THEN
	1556	WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
	1557	& pdtphys
	1558	abort_message='Pas physique n est pas correct '
	1559	! call abort_gcm(modname,abort_message,1)
	1560	dtime=pdtphys
	1561	ENDIF
	1562	IF (nlon .NE. klon) THEN
	1563	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
	1564	& klon
	1565	abort_message='nlon et klon ne sont pas coherents'
	1566	call abort_gcm(modname,abort_message,1)
	1567	ENDIF
	1568	IF (nlev .NE. klev) THEN
	1569	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
	1570	& klev
	1571	abort_message='nlev et klev ne sont pas coherents'
	1572	call abort_gcm(modname,abort_message,1)
	1573	ENDIF
	1574	!c
	1575	IF (dtime*REAL(radpas).GT.21600..AND.cycle_diurne) THEN
	1576	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	1577	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	1578	abort_message='Nbre d appels au rayonnement insuffisant'
	1579	call abort_gcm(modname,abort_message,1)
	1580	ENDIF
	1581	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
	1582	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
	1583	& ok_cvl
	1584	!c
	1585	!cKE43
	1586	!c Initialisation pour la convection de K.E. (sb):
	1587	IF (iflag_con.GE.3) THEN
	1588
	1589	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
	1590	WRITE(lunout,*) &
	1591	& "On va utiliser le melange convectif des traceurs qui"
	1592	WRITE(lunout,*)"est calcule dans convect4.3"
	1593	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
	1594
	1595	DO i = 1, klon
	1596	ema_cbmf(i) = 0.
	1597	ema_pcb(i) = 0.
	1598	ema_pct(i) = 0.
	1599	!c ema_workcbmf(i) = 0.
	1600	ENDDO
	1601	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
	1602	DO i = 1, klon
	1603	ibas_con(i) = 1
	1604	itop_con(i) = 1
	1605	ENDDO
	1606	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
	1607	!c===============================================================================
	1608	!cCR:04.12.07: initialisations poches froides
	1609	!c Controle de ALE et ALP pour la fermeture convective (jyg)
	1610	if (iflag_wake>=1) then
	1611	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
	1612	& ,alp_bl_prescr, ale_bl_prescr)
	1613	!c 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
	1614	!c print*,'apres ini_wake iflag_cldcon=', iflag_cldcon
	1615	endif
	1616
	1617	do i = 1,klon
	1618	Ale_bl(i)=0.
	1619	Alp_bl(i)=0.
	1620	enddo
	1621
	1622	!c================================================================================
	1623	!IM stations CFMIP
	1624	nCFMIP=npCFMIP
	1625	OPEN(98,file='npCFMIP_param.data',status='old', &
	1626	& form='formatted',iostat=iostat)
	1627	if (iostat == 0) then
	1628	READ(98,*,end=998) nCFMIP
	1629	998 CONTINUE
	1630	CLOSE(98)
	1631	CONTINUE
	1632	IF(nCFMIP.GT.npCFMIP) THEN
	1633	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
	1634	CALL abort
	1635	else
	1636	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
	1637	ENDIF
	1638
	1639	!c
	1640	ALLOCATE(tabCFMIP(nCFMIP))
	1641	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
	1642	ALLOCATE(tabijGCM(nCFMIP))
	1643	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
	1644	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
	1645	!c
	1646	!c lecture des nCFMIP stations CFMIP, de leur numero
	1647	!c et des coordonnees geographiques lonCFMIP, latCFMIP
	1648	!c
	1649	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
	1650	&lonCFMIP, latCFMIP)
	1651	!c
	1652	!c identification des
	1653	!c 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
	1654	!c 2) indices points tabijGCM de la grille physique 1d sur klon points
	1655	!c 3) indices iGCM, jGCM de la grille physique 2d
	1656	!c
	1657	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
	1658	&tabijGCM, lonGCM, latGCM, iGCM, jGCM)
	1659	!c
	1660	else
	1661	ALLOCATE(tabijGCM(0))
	1662	ALLOCATE(lonGCM(0), latGCM(0))
	1663	ALLOCATE(iGCM(0), jGCM(0))
	1664	end if
	1665	else
	1666	ALLOCATE(tabijGCM(0))
	1667	ALLOCATE(lonGCM(0), latGCM(0))
	1668	ALLOCATE(iGCM(0), jGCM(0))
	1669	ENDIF
	1670
	1671	DO i=1,klon
	1672	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1673	ENDDO
	1674
	1675	!c34EK
	1676	IF (ok_orodr) THEN
	1677
	1678	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1679	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
	1680	! justement quand ok_orodr = false.
	1681	! ce rugoro est utilise par la couche limite et fait double emploi
	1682	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
	1683	! DO i=1,klon
	1684	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1685	! ENDDO
	1686	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1687	IF (ok_strato) THEN
	1688	CALL SUGWD_strato(klon,klev,paprs,pplay)
	1689	ELSE
	1690	CALL SUGWD(klon,klev,paprs,pplay)
	1691	ENDIF
	1692
	1693	DO i=1,klon
	1694	zuthe(i)=0.
	1695	zvthe(i)=0.
	1696	if(zstd(i).gt.10.)then
	1697	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	1698	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	1699	endif
	1700	ENDDO
	1701	ENDIF
	1702	!c
	1703	!c
	1704	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
	1705	WRITE(lunout,*)'La frequence de lecture surface est de ', &
	1706	& lmt_pas
	1707	!c
	1708	capemaxcels = 't_max(X)'
	1709	t2mincels = 't_min(X)'
	1710	t2maxcels = 't_max(X)'
	1711	tinst = 'inst(X)'
	1712	tave = 'ave(X)'
	1713	!IM cf. AM 081204 BEG
	1714	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
	1715	!IM cf. AM 081204 END
	1716	!c
	1717	!c=============================================================
	1718	!c Initialisation des sorties
	1719	!c=============================================================
	1720
	1721	#ifdef CPP_IOIPSL
	1722
	1723	!$OMP MASTER
	1724	call phys_output_open(rlon,rlat,nCFMIP,tabijGCM, &
	1725	& iGCM,jGCM,lonGCM,latGCM, &
	1726	& jjmp1,nlevSTD,clevSTD, &
	1727	& nbteta, ctetaSTD, dtime,ok_veget, &
	1728	& type_ocean,iflag_pbl,ok_mensuel,ok_journe, &
	1729	& ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
	1730	& read_climoz, phys_out_filestations, &
	1731	& new_aod, aerosol_couple, &
	1732	& flag_aerosol_strat )
	1733	!$OMP END MASTER
	1734	!$OMP BARRIER
	1735
	1736	#ifdef histISCCP
	1737	#include "ini_histISCCP.h"
	1738	#endif
	1739
	1740	#ifdef histNMC
	1741	#include "ini_histhfNMC.h"
	1742	#include "ini_histdayNMC.h"
	1743	#include "ini_histmthNMC.h"
	1744	#endif
	1745
	1746	#include "ini_histday_seri.h"
	1747
	1748	#include "ini_paramLMDZ_phy.h"
	1749
	1750	#endif
	1751	ecrit_reg = ecrit_reg * un_jour
	1752	ecrit_tra = ecrit_tra * un_jour
	1753
	1754	!cXXXPB Positionner date0 pour initialisation de ORCHIDEE
	1755	date0 = jD_ref
	1756	WRITE(,) 'physiq date0 : ',date0
	1757	!c
	1758	!c
	1759	!c
	1760	!c Prescrire l'ozone dans l'atmosphere
	1761	!c
	1762	!c
	1763	!cc DO i = 1, klon
	1764	!cc DO k = 1, klev
	1765	!cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
	1766	!cc ENDDO
	1767	!cc ENDDO
	1768	!c
	1769	IF (type_trac == 'inca') THEN
	1770	#ifdef INCA
	1771	CALL VTe(VTphysiq)
	1772	CALL VTb(VTinca)
	1773	! iii = MOD(NINT(xjour),360)
	1774	! calday = REAL(iii) + jH_cur
	1775	calday = REAL(days_elapsed) + jH_cur
	1776	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
	1777
	1778	CALL chemini( &
	1779	& rg, &
	1780	& ra, &
	1781	& airephy, &
	1782	& rlat, &
	1783	& rlon, &
	1784	& presnivs, &
	1785	& calday, &
	1786	& klon, &
	1787	& nqtot, &
	1788	& pdtphys, &
	1789	& annee_ref, &
	1790	& day_ref, &
	1791	& itau_phy)
	1792
	1793	CALL VTe(VTinca)
	1794	CALL VTb(VTphysiq)
	1795	#endif
	1796	END IF
	1797	!c
	1798	!c
	1799	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1800	! Nouvelle initialisation pour le rayonnement RRTM
	1801	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1802
	1803	call iniradia(klon,klev,paprs(1,1:klev+1))
	1804
	1805	!$OMP single
	1806	if (read_climoz >= 1) then
	1807	call open_climoz(ncid_climoz, press_climoz)
	1808	END IF
	1809	!$OMP end single
	1810	!c
	1811	!IM betaCRF
	1812	pfree=70000. !Pa
	1813	beta_pbl=1.
	1814	beta_free=1.
	1815	lon1_beta=-180.
	1816	lon2_beta=+180.
	1817	lat1_beta=90.
	1818	lat2_beta=-90.
	1819	mskocean_beta=.FALSE.
	1820
	1821	OPEN(99,file='beta_crf.data',status='old', &
	1822	& form='formatted',err=9999)
	1823	READ(99,*,end=9998) pfree
	1824	READ(99,*,end=9998) beta_pbl
	1825	READ(99,*,end=9998) beta_free
	1826	READ(99,*,end=9998) lon1_beta
	1827	READ(99,*,end=9998) lon2_beta
	1828	READ(99,*,end=9998) lat1_beta
	1829	READ(99,*,end=9998) lat2_beta
	1830	READ(99,*,end=9998) mskocean_beta
	1831	9998 Continue
	1832	CLOSE(99)
	1833	9999 Continue
	1834	WRITE(,)'pfree=',pfree
	1835	WRITE(,)'beta_pbl=',beta_pbl
	1836	WRITE(,)'beta_free=',beta_free
	1837	WRITE(,)'lon1_beta=',lon1_beta
	1838	WRITE(,)'lon2_beta=',lon2_beta
	1839	WRITE(,)'lat1_beta=',lat1_beta
	1840	WRITE(,)'lat2_beta=',lat2_beta
	1841	WRITE(,)'mskocean_beta=',mskocean_beta
	1842	ENDIF
	1843	!
	1844	! ************** Fin de IF ( debut ) *************
	1845	!
	1846	!
	1847	! Incrementer le compteur de la physique
	1848	!
	1849	itap = itap + 1
	1850	!c
	1851	!
	1852	! Update fraction of the sub-surfaces (pctsrf) and
	1853	! initialize, where a new fraction has appeared, all variables depending
	1854	! on the surface fraction.
	1855	!
	1856	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
	1857	& pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
	1858
	1859	! Update time and other variables in Reprobus
	1860	IF (type_trac == 'repr') THEN
	1861	#ifdef REPROBUS
	1862	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
	1863	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
	1864	CALL Rtime(debut)
	1865	#endif
	1866	END IF
	1867
	1868
	1869	! Tendances bidons pour les processus qui n'affectent pas certaines
	1870	! variables.
	1871	du0(:,:)=0.
	1872	dv0(:,:)=0.
	1873	dq0(:,:)=0.
	1874	dql0(:,:)=0.
	1875	!c
	1876	!c Mettre a zero des variables de sortie (pour securite)
	1877	!c
	1878	DO i = 1, klon
	1879	d_ps(i) = 0.0
	1880	ENDDO
	1881	DO k = 1, klev
	1882	DO i = 1, klon
	1883	d_t(i,k) = 0.0
	1884	d_u(i,k) = 0.0
	1885	d_v(i,k) = 0.0
	1886	ENDDO
	1887	ENDDO
	1888	DO iq = 1, nqtot
	1889	DO k = 1, klev
	1890	DO i = 1, klon
	1891	d_qx(i,k,iq) = 0.0
	1892	ENDDO
	1893	ENDDO
	1894	ENDDO
	1895	da(:,:)=0.
	1896	mp(:,:)=0.
	1897	phi(:,:,:)=0.
	1898	! RomP >>>
	1899	phi2(:,:,:)=0.
	1900	beta_prec_fisrt(:,:)=0.
	1901	beta_prec(:,:)=0.
	1902	epmlmMm(:,:,:)=0.
	1903	eplaMm(:,:)=0.
	1904	d1a(:,:)=0.
	1905	dam(:,:)=0.
	1906	pmflxr=0.
	1907	pmflxs=0.
	1908	! RomP <<<
	1909
	1910	!c
	1911	!c Ne pas affecter les valeurs entrees de u, v, h, et q
	1912	!c
	1913	DO k = 1, klev
	1914	DO i = 1, klon
	1915	t_seri(i,k) = t(i,k)
	1916	u_seri(i,k) = u(i,k)
	1917	v_seri(i,k) = v(i,k)
	1918	q_seri(i,k) = qx(i,k,ivap)
	1919	ql_seri(i,k) = qx(i,k,iliq)
	1920	qs_seri(i,k) = 0.
	1921	ENDDO
	1922	ENDDO
	1923	tke0(:,:)=pbl_tke(:,:,is_ave)
	1924	IF (nqtot.GE.3) THEN
	1925	DO iq = 3, nqtot
	1926	DO k = 1, klev
	1927	DO i = 1, klon
	1928	tr_seri(i,k,iq-2) = qx(i,k,iq)
	1929	ENDDO
	1930	ENDDO
	1931	ENDDO
	1932	ELSE
	1933	DO k = 1, klev
	1934	DO i = 1, klon
	1935	tr_seri(i,k,1) = 0.0
	1936	ENDDO
	1937	ENDDO
	1938	ENDIF
	1939	!C
	1940	DO i = 1, klon
	1941	ztsol(i) = 0.
	1942	ENDDO
	1943	DO nsrf = 1, nbsrf
	1944	DO i = 1, klon
	1945	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
	1946	ENDDO
	1947	ENDDO
	1948	!IM
	1949	IF (ip_ebil_phy.ge.1) THEN
	1950	ztit='after dynamic'
	1951	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
	1952	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	1953	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1954	!C Comme les tendances de la physique sont ajoute dans la dynamique,
	1955	!C on devrait avoir que la variation d'entalpie par la dynamique
	1956	!C est egale a la variation de la physique au pas de temps precedent.
	1957	!C Donc la somme de ces 2 variations devrait etre nulle.
	1958	call diagphy(airephy,ztit,ip_ebil_phy &
	1959	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	1960	& , zero_v, zero_v, zero_v, ztsol &
	1961	& , d_h_vcol+d_h_vcol_phy, d_qt, 0. &
	1962	& , fs_bound, fq_bound )
	1963	END IF
	1964
	1965	!c Diagnostiquer la tendance dynamique
	1966	!c
	1967	IF (ancien_ok) THEN
	1968	DO k = 1, klev
	1969	DO i = 1, klon
	1970	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
	1971	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
	1972	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
	1973	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
	1974	ENDDO
	1975	ENDDO
	1976	!!! RomP >>> td dyn traceur
	1977	IF (nqtot.GE.3) THEN
	1978	DO iq = 3, nqtot
	1979	DO k = 1, klev
	1980	DO i = 1, klon
	1981	d_tr_dyn(i,k,iq-2)= &
	1982	& (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
	1983	! iiq=niadv(iq)
	1984	! print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
	1985	ENDDO
	1986	ENDDO
	1987	ENDDO
	1988	ENDIF
	1989	!!! RomP <<<
	1990	ELSE
	1991	DO k = 1, klev
	1992	DO i = 1, klon
	1993	d_u_dyn(i,k) = 0.0
	1994	d_v_dyn(i,k) = 0.0
	1995	d_t_dyn(i,k) = 0.0
	1996	d_q_dyn(i,k) = 0.0
	1997	ENDDO
	1998	ENDDO
	1999	!!! RomP >>> td dyn traceur
	2000	IF (nqtot.GE.3) THEN
	2001	DO iq = 3, nqtot
	2002	DO k = 1, klev
	2003	DO i = 1, klon
	2004	d_tr_dyn(i,k,iq-2)= 0.0
	2005	ENDDO
	2006	ENDDO
	2007	ENDDO
	2008	ENDIF
	2009	!!! RomP <<<
	2010	ancien_ok = .TRUE.
	2011	ENDIF
	2012	!c
	2013	!c Ajouter le geopotentiel du sol:
	2014	!c
	2015	DO k = 1, klev
	2016	DO i = 1, klon
	2017	zphi(i,k) = pphi(i,k) + pphis(i)
	2018	ENDDO
	2019	ENDDO
	2020	!c
	2021	!c Verifier les temperatures
	2022	!c
	2023	!IM BEG
	2024	IF (check) THEN
	2025	amn=MIN(ftsol(1,is_ter),1000.)
	2026	amx=MAX(ftsol(1,is_ter),-1000.)
	2027	DO i=2, klon
	2028	amn=MIN(ftsol(i,is_ter),amn)
	2029	amx=MAX(ftsol(i,is_ter),amx)
	2030	ENDDO
	2031	!c
	2032	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
	2033	ENDIF !(check) THEN
	2034
	2035	!IM END
	2036	!c
	2037
	2038	CALL hgardfou(t_seri,ftsol,'debutphy')
	2039	!c
	2040	!IM BEG
	2041	IF (check) THEN
	2042	amn=MIN(ftsol(1,is_ter),1000.)
	2043	amx=MAX(ftsol(1,is_ter),-1000.)
	2044	DO i=2, klon
	2045	amn=MIN(ftsol(i,is_ter),amn)
	2046	amx=MAX(ftsol(i,is_ter),amx)
	2047	ENDDO
	2048	!c
	2049	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
	2050	ENDIF !(check) THEN
	2051
	2052	!IM END
	2053	!c
	2054	!c Mettre en action les conditions aux limites (albedo, sst, etc.).
	2055	!c Prescrire l'ozone et calculer l'albedo sur l'ocean.
	2056	!c
	2057	if (read_climoz >= 1) then
	2058	!C Ozone from a file
	2059	! Update required ozone index:
	2060	ro3i = int((days_elapsed + jh_cur - jh_1jan) &
	2061	& / ioget_year_len(year_cur) * 360.) + 1
	2062	if (ro3i == 361) ro3i = 360
	2063	!C (This should never occur, except perhaps because of roundup
	2064	!C error. See documentation.)
	2065	if (ro3i /= co3i) then
	2066	!C Update ozone field:
	2067	if (read_climoz == 1) then
	2068	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
	2069	& press_in_edg=press_climoz, paprs=paprs, v3=wo)
	2070	else
	2071	!C read_climoz == 2
	2072	call regr_pr_av(ncid_climoz, &
	2073	& (/"tro3 ", "tro3_daylight"/), &
	2074	& julien=ro3i, press_in_edg=press_climoz, paprs=paprs, &
	2075	& v3=wo)
	2076	end if
	2077	! Convert from mole fraction of ozone to column density of ozone in a
	2078	! cell, in kDU:
	2079	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) &
	2080	& * rmo3 / rmd * zmasse / dobson_u / 1e3
	2081	!C (By regridding ozone values for LMDZ only once every 360th of
	2082	!C year, we have already neglected the variation of pressure in one
	2083	!C 360th of year. So do not recompute "wo" at each time step even if
	2084	!C "zmasse" changes a little.)
	2085	co3i = ro3i
	2086	end if
	2087	elseif (MOD(itap-1,lmt_pas) == 0) THEN
	2088	!C Once per day, update ozone from Royer:
	2089	wo(:, :, 1) = ozonecm(rlat, paprs, rjour=real(days_elapsed+1))
	2090	ENDIF
	2091
	2092	!c
	2093	!c Re-evaporer l'eau liquide nuageuse
	2094	!c
	2095	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
	2096	DO i = 1, klon
	2097	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	2098	!c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	2099	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	2100	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
	2101	zb = MAX(0.0,ql_seri(i,k))
	2102	za = - MAX(0.0,ql_seri(i,k)) &
	2103	& * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
	2104	t_seri(i,k) = t_seri(i,k) + za
	2105	q_seri(i,k) = q_seri(i,k) + zb
	2106	ql_seri(i,k) = 0.0
	2107	d_t_eva(i,k) = za
	2108	d_q_eva(i,k) = zb
	2109	ENDDO
	2110	ENDDO
	2111	!IM
	2112	IF (ip_ebil_phy.ge.2) THEN
	2113	ztit='after reevap'
	2114	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
	2115	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2116	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2117	call diagphy(airephy,ztit,ip_ebil_phy &
	2118	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	2119	& , zero_v, zero_v, zero_v, ztsol &
	2120	& , d_h_vcol, d_qt, d_ec &
	2121	& , fs_bound, fq_bound )
	2122	!C
	2123	END IF
	2124
	2125	!c
	2126	!c=========================================================================
	2127	! Calculs de l'orbite.
	2128	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
	2129	! doit donc etre plac\'e avant radlwsw et pbl_surface
	2130
	2131	!!! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2132	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
	2133	day_since_equinox = (jD_cur + jH_cur) - jD_eq
	2134	!
	2135	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
	2136	! solarlong0
	2137	if (solarlong0<-999.) then
	2138	if (new_orbit) then
	2139	! calcul selon la routine utilisee pour les planetes
	2140	call solarlong(day_since_equinox, zlongi, dist)
	2141	else
	2142	! calcul selon la routine utilisee pour l'AR4
	2143	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
	2144	endif
	2145	else
	2146	zlongi=solarlong0 ! longitude solaire vraie
	2147	dist=1. ! distance au soleil / moyenne
	2148	endif
	2149	if(prt_level.ge.1) &
[6]	2150	& write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist,&
	2151	' jD_cur: ',jD_cur,' jH_cur: ',jH_cur,' days_elapsed: ', &
	2152	REAL(days_elapsed+1)
[1]	2153
	2154
	2155	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2156	! Calcul de l'ensoleillement :
	2157	! ============================
	2158	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
	2159	! l'annee a partir d'une formule analytique.
	2160	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
	2161	! non nul aux poles.
	2162	IF (abs(solarlong0-1000.)<1.e-4) then
[6]	2163	PRINT *,' Lluis calling zenang_an'
[1]	2164	call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
	2165	ELSE
	2166	! Avec ou sans cycle diurne
	2167	IF (cycle_diurne) THEN
	2168	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
	2169	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
[6]	2170	PRINT *,' Lluis: radpas: ',radpas,' fract: ',fract(llp), &
	2171	' zdtime: ',zdtime,' jH_cur: ',jH_cur,' rmu0: ',rmu0(llp),&
	2172	' dtime: ',dtime,' nbapp_rad: ',nbapp_rad
	2173	PRINT *,' Lluis @rlat@: ',rlat
	2174	PRINT *,' Lluis @rlon@: ',rlon
	2175	PRINT *,' Lluis @rmu0@: ',rmu0
	2176	PRINT *,' Lluis @fract@: ',fract
[1]	2177	ELSE
	2178	CALL angle(zlongi, rlat, fract, rmu0)
	2179	ENDIF
	2180	ENDIF
[6]	2181	PRINT *,' Lluis zenang_an rlat: ',rlat(llp), ' rlon:',rlon(llp), &
	2182	' rmu0: ',rmu0(llp),' fract: ',fract(llp),' jH_cur: ',jH_cur, &
	2183	' cycle_diurne: ',cycle_diurne
[1]	2184
	2185	if (mydebug) then
	2186	call writefield_phy('u_seri',u_seri,llm)
	2187	call writefield_phy('v_seri',v_seri,llm)
	2188	call writefield_phy('t_seri',t_seri,llm)
	2189	call writefield_phy('q_seri',q_seri,llm)
	2190	endif
	2191
	2192	!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2193	!c Appel au pbl_surface : Planetary Boudary Layer et Surface
	2194	!c Cela implique tous les interactions des sous-surfaces et la partie diffusion
	2195	!c turbulent du couche limit.
	2196	!c
	2197	!c Certains varibales de sorties de pbl_surface sont utiliser que pour
	2198	!c ecriture des fihiers hist_XXXX.nc, ces sont :
	2199	!c qsol, zq2m, s_pblh, s_lcl,
	2200	!c s_capCL, s_oliqCL, s_cteiCL,s_pblT,
	2201	!c s_therm, s_trmb1, s_trmb2, s_trmb3,
	2202	!c zxrugs, zu10m, zv10m, fder,
	2203	!c zxqsurf, rh2m, zxfluxu, zxfluxv,
	2204	!c frugs, agesno, fsollw, fsolsw,
	2205	!c d_ts, fevap, fluxlat, t2m,
	2206	!c wfbils, wfbilo, fluxt, fluxu, fluxv,
	2207	!c
	2208	!c Certains ne sont pas utiliser du tout :
	2209	!c dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
	2210	!c
	2211
	2212	!c Calcul de l'humidite de saturation au niveau du sol
	2213
	2214
	2215
	2216	if (iflag_pbl/=0) then
	2217
	2218	PRINT *,' Lluis before pbl_surface qsol: ',qsol(llp), &
[6]	2219	' rmu0: ',rmu0(llp),' jH_cur: ',jH_cur
[1]	2220
	2221	CALL pbl_surface( &
	2222	& dtime, date0, itap, days_elapsed+1, &
	2223	& debut, lafin, &
	2224	& rlon, rlat, rugoro, rmu0, &
	2225	& rain_fall, snow_fall, solsw, sollw, &
	2226	& t_seri, q_seri, u_seri, v_seri, &
	2227	& pplay, paprs, pctsrf, &
	2228	& ftsol, falb1, falb2, ustar, u10m, v10m, &
	2229	& sollwdown, cdragh, cdragm, u1, v1, &
	2230	& albsol1, albsol2, sens, evap, &
	2231	& zxtsol, zxfluxlat, zt2m, qsat2m, &
	2232	& d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
	2233	& coefh, coefm, slab_wfbils, &
	2234	& qsol, zq2m, s_pblh, s_lcl, &
	2235	& s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
	2236	& s_therm, s_trmb1, s_trmb2, s_trmb3, &
	2237	& zxrugs, zustar, zu10m, zv10m, fder, &
	2238	& zxqsurf, rh2m, zxfluxu, zxfluxv, &
	2239	& frugs, agesno, fsollw, fsolsw, &
	2240	& d_ts, fevap, fluxlat, t2m, &
	2241	& wfbils, wfbilo, fluxt, fluxu, fluxv, &
	2242	& dsens, devap, zxsnow, &
	2243	& zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
	2244
[6]	2245	PRINT *,' Lluis after pbl_surface qsol: ',qsol(llp), &
	2246	' rmu0: ',rmu0(llp),' jH_cur: ',jH_cur
[1]	2247
	2248	!-----------------------------------------------------------------------------------------
	2249	! ajout des tendances de la diffusion turbulente
	2250	CALL add_phys_tend &
	2251	& (d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,'vdf')
	2252	!-----------------------------------------------------------------------------------------
	2253	if (mydebug) then
	2254	call writefield_phy('u_seri',u_seri,llm)
	2255	call writefield_phy('v_seri',v_seri,llm)
	2256	call writefield_phy('t_seri',t_seri,llm)
	2257	call writefield_phy('q_seri',q_seri,llm)
	2258	endif
	2259
	2260	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
	2261	& t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
	2262
	2263
	2264	IF (ip_ebil_phy.ge.2) THEN
	2265	ztit='after surface_main'
	2266	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	2267	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2268	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2269	call diagphy(airephy,ztit,ip_ebil_phy &
	2270	& , zero_v, zero_v, zero_v, zero_v, sens &
	2271	& , evap , zero_v, zero_v, ztsol &
	2272	& , d_h_vcol, d_qt, d_ec &
	2273	& , fs_bound, fq_bound )
	2274	END IF
	2275
	2276	ENDIF
	2277	!c =================================================================== c
	2278	!c Calcul de Qsat
	2279
	2280	DO k = 1, klev
	2281	DO i = 1, klon
	2282	zx_t = t_seri(i,k)
	2283	IF (thermcep) THEN
	2284	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	2285	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	2286	zx_qs = MIN(0.5,zx_qs)
	2287	zcor = 1./(1.-retv*zx_qs)
	2288	zx_qs = zx_qs*zcor
	2289	ELSE
	2290	IF (zx_t.LT.t_coup) THEN
	2291	zx_qs = qsats(zx_t)/pplay(i,k)
	2292	ELSE
	2293	zx_qs = qsatl(zx_t)/pplay(i,k)
	2294	ENDIF
	2295	ENDIF
	2296	zqsat(i,k)=zx_qs
	2297	ENDDO
	2298	ENDDO
	2299
	2300	if (prt_level.ge.1) then
	2301	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
	2302	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
	2303	endif
	2304	!c
	2305	!c Appeler la convection (au choix)
	2306	!c
	2307	DO k = 1, klev
	2308	DO i = 1, klon
	2309	conv_q(i,k) = d_q_dyn(i,k) &
	2310	& + d_q_vdf(i,k)/dtime
	2311	conv_t(i,k) = d_t_dyn(i,k) &
	2312	& + d_t_vdf(i,k)/dtime
	2313	ENDDO
	2314	ENDDO
	2315	IF (check) THEN
	2316	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2317	WRITE(lunout,*) "avantcon=", za
	2318	ENDIF
	2319	zx_ajustq = .FALSE.
	2320	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
	2321	IF (zx_ajustq) THEN
	2322	DO i = 1, klon
	2323	z_avant(i) = 0.0
	2324	ENDDO
	2325	DO k = 1, klev
	2326	DO i = 1, klon
	2327	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
	2328	& *(paprs(i,k)-paprs(i,k+1))/RG
	2329	ENDDO
	2330	ENDDO
	2331	ENDIF
	2332
	2333	!c Calcule de vitesse verticale a partir de flux de masse verticale
	2334	DO k = 1, klev
	2335	DO i = 1, klon
	2336	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
	2337	END DO
	2338	END DO
	2339	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
	2340	& omega(igout, :)
	2341
	2342	IF (iflag_con.EQ.1) THEN
	2343	abort_message ='reactiver le call conlmd dans physiq.F'
	2344	CALL abort_gcm (modname,abort_message,1)
	2345	!c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
	2346	!c . d_t_con, d_q_con,
	2347	!c . rain_con, snow_con, ibas_con, itop_con)
	2348	ELSE IF (iflag_con.EQ.2) THEN
	2349	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
	2350	& conv_t, conv_q, -evap, omega, &
	2351	& d_t_con, d_q_con, rain_con, snow_con, &
	2352	& pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
	2353	& kcbot, kctop, kdtop, pmflxr, pmflxs)
	2354	d_u_con = 0.
	2355	d_v_con = 0.
	2356
	2357	WHERE (rain_con < 0.) rain_con = 0.
	2358	WHERE (snow_con < 0.) snow_con = 0.
	2359	DO i = 1, klon
	2360	ibas_con(i) = klev+1 - kcbot(i)
	2361	itop_con(i) = klev+1 - kctop(i)
	2362	ENDDO
	2363	ELSE IF (iflag_con.GE.3) THEN
	2364	!c nb of tracers for the KE convection:
	2365	!c MAF la partie traceurs est faite dans phytrac
	2366	!c on met ntra=1 pour limiter les appels mais on peut
	2367	!c supprimer les calculs / ftra.
	2368	ntra = 1
	2369
	2370	!c=====================================================================================
	2371	!cajout pour la parametrisation des poches froides:
	2372	!ccalcul de t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
	2373	do k=1,klev
	2374	do i=1,klon
	2375	if (iflag_wake>=1) then
	2376	t_wake(i,k) = t_seri(i,k) &
	2377	& +(1-wake_s(i))*wake_deltat(i,k)
	2378	q_wake(i,k) = q_seri(i,k) &
	2379	& +(1-wake_s(i))*wake_deltaq(i,k)
	2380	t_undi(i,k) = t_seri(i,k) &
	2381	& -wake_s(i)*wake_deltat(i,k)
	2382	q_undi(i,k) = q_seri(i,k) &
	2383	& -wake_s(i)*wake_deltaq(i,k)
	2384	else
	2385	t_wake(i,k) = t_seri(i,k)
	2386	q_wake(i,k) = q_seri(i,k)
	2387	t_undi(i,k) = t_seri(i,k)
	2388	q_undi(i,k) = q_seri(i,k)
	2389	endif
	2390	enddo
	2391	enddo
	2392
	2393	!cc-- Calcul de l'energie disponible ALE (J/kg) et de la puissance disponible ALP (W/m2)
	2394	!cc-- pour le soulevement des particules dans le modele convectif
	2395	!c
	2396	do i = 1,klon
	2397	ALE(i) = 0.
	2398	ALP(i) = 0.
	2399	enddo
	2400	!c
	2401	!ccalcul de ale_wake et alp_wake
	2402	if (iflag_wake>=1) then
	2403	if (itap .le. it_wape_prescr) then
	2404	do i = 1,klon
	2405	ale_wake(i) = wape_prescr
	2406	alp_wake(i) = fip_prescr
	2407	enddo
	2408	else
	2409	do i = 1,klon
	2410	!cjyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
	2411	!ccc ale_wake(i) = 0.5wake_cstar(i)*2
	2412	ale_wake(i) = wake_pe(i)
	2413	alp_wake(i) = wake_fip(i)
	2414	enddo
	2415	endif
	2416	else
	2417	do i = 1,klon
	2418	ale_wake(i) = 0.
	2419	alp_wake(i) = 0.
	2420	enddo
	2421	endif
	2422	!ccombinaison avec ale et alp de couche limite: constantes si pas de couplage, valeurs calculees
	2423	!cdans le thermique sinon
	2424	if (iflag_coupl.eq.0) then
	2425	if (debut.and.prt_level.gt.9) &
	2426	& WRITE(lunout,*)'ALE et ALP imposes'
	2427	do i = 1,klon
	2428	!con ne couple que ale
	2429	!c ALE(i) = max(ale_wake(i),Ale_bl(i))
	2430	ALE(i) = max(ale_wake(i),ale_bl_prescr)
	2431	!con ne couple que alp
	2432	!c ALP(i) = alp_wake(i) + Alp_bl(i)
	2433	ALP(i) = alp_wake(i) + alp_bl_prescr
	2434	enddo
	2435	else
	2436	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
	2437	! do i = 1,klon
	2438	! ALE(i) = max(ale_wake(i),Ale_bl(i))
	2439	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
	2440	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
	2441	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
	2442	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
	2443	! enddo
	2444
	2445	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2446	! Modif FH 2010/04/27. Sans doute temporaire.
	2447	! Deux options pour le alp_offset : constant si >?? 0 ou proportionnel ??a
	2448	! w si <0
	2449	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2450	do i = 1,klon
	2451	ALE(i) = max(ale_wake(i),Ale_bl(i))
	2452	!ccc nrlmd le 10/04/2012----------Stochastic triggering--------------
	2453	if (iflag_trig_bl.ge.1) then
	2454	ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
	2455	endif
	2456	!ccc fin nrlmd le 10/04/2012
	2457	if (alp_offset>=0.) then
	2458	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
	2459	else
	2460	ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
	2461	if (alp(i)<0.) then
	2462	print*,'ALP ',alp(i),alp_wake(i) &
	2463	& ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
	2464	endif
	2465	endif
	2466	enddo
	2467	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2468
	2469	endif
	2470	do i=1,klon
	2471	if (alp(i)>alp_max) then
	2472	IF(prt_level>9)WRITE(lunout,*) &
	2473	& 'WARNING SUPER ALP (seuil=',alp_max, &
	2474	& '): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
	2475	alp(i)=alp_max
	2476	endif
	2477	if (ale(i)>ale_max) then
	2478	IF(prt_level>9)WRITE(lunout,*) &
	2479	& 'WARNING SUPER ALE (seuil=',ale_max, &
	2480	& '): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
	2481	ale(i)=ale_max
	2482	endif
	2483	enddo
	2484
	2485	!cfin calcul ale et alp
	2486	!c=================================================================================================
	2487
	2488
	2489	!c sb, oct02:
	2490	!c Schema de convection modularise et vectorise:
	2491	!c (driver commun aux versions 3 et 4)
	2492	!c
	2493	IF (ok_cvl) THEN ! new driver for convectL
	2494
	2495	IF (type_trac == 'repr') THEN
	2496	nbtr_tmp=ntra
	2497	ELSE
	2498	nbtr_tmp=nbtr
	2499	END IF
	2500	CALL concvl (iflag_con,iflag_clos, &
	2501	& dtime,paprs,pplay,t_undi,q_undi, &
	2502	& t_wake,q_wake,wake_s, &
	2503	& u_seri,v_seri,tr_seri,nbtr_tmp, &
	2504	& ALE,ALP, &
	2505	& ema_work1,ema_work2, &
	2506	& d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
	2507	& rain_con, snow_con, ibas_con, itop_con, sigd, &
	2508	& ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
	2509	& Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
	2510	& pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
	2511	! RomP >>>
	2512	!! . pmflxr,pmflxs,da,phi,mp,
	2513	!! . ftd,fqd,lalim_conv,wght_th) &
	2514	& pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
	2515	& ftd,fqd,lalim_conv,wght_th, &
	2516	& ev, ep,epmlmMm,eplaMm, &
	2517	& wdtrainA,wdtrainM)
	2518	! RomP <<<
	2519
	2520	!IM begin
	2521	!c print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
	2522	!c .dnwd0(1,1),ftd(1,1),fqd(1,1)
	2523	!IM end
	2524	!IM cf. FH
	2525	clwcon0=qcondc
	2526	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
	2527
	2528	do i = 1, klon
	2529	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
	2530	enddo
	2531
	2532	ELSE ! ok_cvl
	2533
	2534	!c MAF conema3 ne contient pas les traceurs
	2535	CALL conema3 (dtime, &
	2536	& paprs,pplay,t_seri,q_seri, &
	2537	& u_seri,v_seri,tr_seri,ntra, &
	2538	& ema_work1,ema_work2, &
	2539	& d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
	2540	& rain_con, snow_con, ibas_con, itop_con, &
	2541	& upwd,dnwd,dnwd0,bas,top, &
	2542	& Ma,cape,tvp,rflag, &
	2543	& pbase &
	2544	& ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
	2545	& ,clwcon0)
	2546
	2547	ENDIF ! ok_cvl
	2548
	2549	!c
	2550	!c Correction precip
	2551	rain_con = rain_con * cvl_corr
	2552	snow_con = snow_con * cvl_corr
	2553	!c
	2554
	2555	IF (.NOT. ok_gust) THEN
	2556	do i = 1, klon
	2557	wd(i)=0.0
	2558	enddo
	2559	ENDIF
	2560
	2561	!c =================================================================== c
	2562	!c Calcul des proprietes des nuages convectifs
	2563	!c
	2564
	2565	!c calcul des proprietes des nuages convectifs
	2566	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
	2567	call clouds_gno &
	2568	& (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
	2569
	2570	!c =================================================================== c
	2571
	2572	DO i = 1, klon
	2573	itop_con(i) = min(max(itop_con(i),1),klev)
	2574	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
	2575	ENDDO
	2576
	2577	DO i = 1, klon
	2578	ema_pcb(i) = paprs(i,ibas_con(i))
	2579	ENDDO
	2580	DO i = 1, klon
	2581	! L'idicage de itop_con peut cacher un pb potentiel
	2582	! FH sous la dictee de JYG, CR
	2583	ema_pct(i) = paprs(i,itop_con(i)+1)
	2584
	2585	if (itop_con(i).gt.klev-3) then
	2586	if(prt_level >= 9) then
	2587	write(lunout,*)'La convection monte trop haut '
	2588	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
	2589	endif
	2590	endif
	2591	ENDDO
	2592	ELSE IF (iflag_con.eq.0) THEN
	2593	write(lunout,*) 'On n appelle pas la convection'
	2594	clwcon0=0.
	2595	rnebcon0=0.
	2596	d_t_con=0.
	2597	d_q_con=0.
	2598	d_u_con=0.
	2599	d_v_con=0.
	2600	rain_con=0.
	2601	snow_con=0.
	2602	bas=1
	2603	top=1
	2604	ELSE
	2605	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
	2606	CALL abort
	2607	ENDIF
	2608
	2609	!c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
	2610	!c . d_u_con, d_v_con)
	2611
	2612	!-----------------------------------------------------------------------------------------
	2613	! ajout des tendances de la diffusion turbulente
	2614	CALL add_phys_tend(d_u_con,d_v_con,d_t_con,d_q_con,dql0,'con')
	2615	!-----------------------------------------------------------------------------------------
	2616
	2617	if (mydebug) then
	2618	call writefield_phy('u_seri',u_seri,llm)
	2619	call writefield_phy('v_seri',v_seri,llm)
	2620	call writefield_phy('t_seri',t_seri,llm)
	2621	call writefield_phy('q_seri',q_seri,llm)
	2622	endif
	2623
	2624	!IM
	2625	IF (ip_ebil_phy.ge.2) THEN
	2626	ztit='after convect'
	2627	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	2628	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2629	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2630	call diagphy(airephy,ztit,ip_ebil_phy &
	2631	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	2632	& , zero_v, rain_con, snow_con, ztsol &
	2633	& , d_h_vcol, d_qt, d_ec &
	2634	& , fs_bound, fq_bound )
	2635	END IF
	2636	!C
	2637	IF (check) THEN
	2638	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	2639	WRITE(lunout,*)"aprescon=", za
	2640	zx_t = 0.0
	2641	za = 0.0
	2642	DO i = 1, klon
	2643	za = za + airephy(i)/REAL(klon)
	2644	zx_t = zx_t + (rain_con(i)+ &
	2645	& snow_con(i))*airephy(i)/REAL(klon)
	2646	ENDDO
	2647	zx_t = zx_t/za*dtime
	2648	WRITE(lunout,*)"Precip=", zx_t
	2649	ENDIF
	2650	IF (zx_ajustq) THEN
	2651	DO i = 1, klon
	2652	z_apres(i) = 0.0
	2653	ENDDO
	2654	DO k = 1, klev
	2655	DO i = 1, klon
	2656	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
	2657	& *(paprs(i,k)-paprs(i,k+1))/RG
	2658	ENDDO
	2659	ENDDO
	2660	DO i = 1, klon
	2661	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
	2662	& /z_apres(i)
	2663	ENDDO
	2664	DO k = 1, klev
	2665	DO i = 1, klon
	2666	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
	2667	& z_factor(i).LT.(1.0-1.0E-08)) THEN
	2668	q_seri(i,k) = q_seri(i,k) * z_factor(i)
	2669	ENDIF
	2670	ENDDO
	2671	ENDDO
	2672	ENDIF
	2673	zx_ajustq=.FALSE.
	2674
	2675	!c
	2676	!c=============================================================================
	2677	!cRR:Evolution de la poche froide: on ne fait pas de separation wake/env
	2678	!cpour la couche limite diffuse pour l instant
	2679	!c
	2680	if (iflag_wake>=1) then
	2681	DO k=1,klev
	2682	DO i=1,klon
	2683	dt_dwn(i,k) = ftd(i,k)
	2684	wdt_PBL(i,k) = 0.
	2685	dq_dwn(i,k) = fqd(i,k)
	2686	wdq_PBL(i,k) = 0.
	2687	M_dwn(i,k) = dnwd0(i,k)
	2688	M_up(i,k) = upwd(i,k)
	2689	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
	2690	udt_PBL(i,k) = 0.
	2691	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
	2692	udq_PBL(i,k) = 0.
	2693	ENDDO
	2694	ENDDO
	2695
	2696	if (iflag_wake==2) then
	2697	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
	2698	DO k = 1,klev
	2699	dt_dwn(:,k)= dt_dwn(:,k)+ &
	2700	& ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
	2701	dq_dwn(:,k)= dq_dwn(:,k)+ &
	2702	& ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
	2703	ENDDO
	2704	endif
	2705	!c
	2706	!ccalcul caracteristiques de la poche froide
	2707	call calWAKE (paprs,pplay,dtime &
	2708	& ,t_seri,q_seri,omega &
	2709	& ,dt_dwn,dq_dwn,M_dwn,M_up &
	2710	& ,dt_a,dq_a,sigd &
	2711	& ,wdt_PBL,wdq_PBL &
	2712	& ,udt_PBL,udq_PBL &
	2713	& ,wake_deltat,wake_deltaq,wake_dth &
	2714	& ,wake_h,wake_s,wake_dens &
	2715	& ,wake_pe,wake_fip,wake_gfl &
	2716	& ,dt_wake,dq_wake &
	2717	& ,wake_k, t_undi,q_undi &
	2718	& ,wake_omgbdth,wake_dp_omgb &
	2719	& ,wake_dtKE,wake_dqKE &
	2720	& ,wake_dtPBL,wake_dqPBL &
	2721	& ,wake_omg,wake_dp_deltomg &
	2722	& ,wake_spread,wake_Cstar,wake_d_deltat_gw &
	2723	& ,wake_ddeltat,wake_ddeltaq)
	2724	!c
	2725	!-----------------------------------------------------------------------------------------
	2726	! ajout des tendances des poches froides
	2727	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
	2728	! coherent avec les autres d_t_...
	2729	d_t_wake(:,:)=dt_wake(:,:)*dtime
	2730	d_q_wake(:,:)=dq_wake(:,:)*dtime
	2731	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,'wake')
	2732	!-----------------------------------------------------------------------------------------
	2733
	2734	endif
	2735	!c
	2736	!c===================================================================
	2737	!cJYG
	2738	IF (ip_ebil_phy.ge.2) THEN
	2739	ztit='after wake'
	2740	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	2741	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2742	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2743	call diagphy(airephy,ztit,ip_ebil_phy &
	2744	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	2745	& , zero_v, zero_v, zero_v, ztsol &
	2746	& , d_h_vcol, d_qt, d_ec &
	2747	& , fs_bound, fq_bound )
	2748	END IF
	2749
	2750	!c print*,'apres callwake iflag_cldcon=', iflag_cldcon
	2751	!c
	2752	!c===================================================================
	2753	!c Convection seche (thermiques ou ajustement)
	2754	!c===================================================================
	2755	!c
	2756	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
	2757	& ,seuil_inversion,weak_inversion,dthmin)
	2758
	2759
	2760
	2761	d_t_ajsb(:,:)=0.
	2762	d_q_ajsb(:,:)=0.
	2763	d_t_ajs(:,:)=0.
	2764	d_u_ajs(:,:)=0.
	2765	d_v_ajs(:,:)=0.
	2766	d_q_ajs(:,:)=0.
	2767	clwcon0th(:,:)=0.
	2768	!c
	2769	!c fm_therm(:,:)=0.
	2770	!c entr_therm(:,:)=0.
	2771	!c detr_therm(:,:)=0.
	2772	!c
	2773	IF(prt_level>9)WRITE(lunout,*) &
	2774	& 'AVANT LA CONVECTION SECHE , iflag_thermals=' &
	2775	& ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
	2776	if(iflag_thermals.lt.0) then
	2777	!c Rien
	2778	!c ====
	2779	IF(prt_level>9)WRITE(lunout,*)'pas de convection'
	2780
	2781
	2782	else
	2783
	2784	!c Thermiques
	2785	!c ==========
	2786	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
	2787	& ,iflag_thermals,' nsplit_thermals=',nsplit_thermals
	2788
	2789
	2790	!ccc nrlmd le 10/04/2012
	2791	DO k=1,klev+1
	2792	DO i=1,klon
	2793	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
	2794	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
	2795	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
	2796	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
	2797	ENDDO
	2798	ENDDO
	2799	!ccc fin nrlmd le 10/04/2012
	2800
	2801	if (iflag_thermals>=1) then
	2802	call calltherm(pdtphys &
	2803	& ,pplay,paprs,pphi,weak_inversion &
	2804	& ,u_seri,v_seri,t_seri,q_seri,zqsat,debut &
	2805	& ,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
	2806	& ,fm_therm,entr_therm,detr_therm &
	2807	& ,zqasc,clwcon0th,lmax_th,ratqscth &
	2808	& ,ratqsdiff,zqsatth &
	2809	!con rajoute ale et alp, et les caracteristiques de la couche alim &
	2810	& ,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
	2811	& ,ztv,zpspsk,ztla,zthl &
	2812	!ccc nrlmd le 10/04/2012 &
	2813	& ,pbl_tke_input,pctsrf,omega,airephy &
	2814	& ,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
	2815	& ,n2,s2,ale_bl_stat &
	2816	& ,therm_tke_max,env_tke_max &
	2817	& ,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
	2818	& ,alp_bl_conv,alp_bl_stat &
	2819	!ccc fin nrlmd le 10/04/2012 &
	2820	& ,zqla,ztva )
	2821
	2822	!ccc nrlmd le 10/04/2012
	2823	!c-----------Stochastic triggering-----------
	2824	if (iflag_trig_bl.ge.1) then
	2825	!c
	2826	IF (prt_level .GE. 10) THEN
	2827	print *,'cin, ale_bl_stat, alp_bl_stat ', &
	2828	& cin, ale_bl_stat, alp_bl_stat
	2829	ENDIF
	2830
	2831	!c----Initialisations
	2832	do i=1,klon
	2833	proba_notrig(i)=1.
	2834	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
	2835	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
	2836	tau_trig(i)=tau_trig_shallow
	2837	else
	2838	tau_trig(i)=tau_trig_deep
	2839	endif
	2840	enddo
	2841	!c
	2842	IF (prt_level .GE. 10) THEN
	2843	print *,'random_notrig, tau_trig ', &
	2844	& random_notrig, tau_trig
	2845	print *,'s_trig,s2,n2 ', &
	2846	& s_trig,s2,n2
	2847	ENDIF
	2848
	2849	!c----Tirage al\'eatoire et calcul de ale_bl_trig
	2850	do i=1,klon
	2851	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
	2852	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
	2853	& (n2(i)*dtime/tau_trig(i))
	2854	!c print *, 'proba_notrig(i) ',proba_notrig(i)
	2855	if (random_notrig(i) .ge. proba_notrig(i)) then
	2856	ale_bl_trig(i)=ale_bl_stat(i)
	2857	else
	2858	ale_bl_trig(i)=0.
	2859	endif
	2860	else
	2861	proba_notrig(i)=1.
	2862	random_notrig(i)=0.
	2863	ale_bl_trig(i)=0.
	2864	endif
	2865	enddo
	2866	!c
	2867	IF (prt_level .GE. 10) THEN
	2868	print *,'proba_notrig, ale_bl_trig ', &
	2869	& proba_notrig, ale_bl_trig
	2870	ENDIF
	2871
	2872	endif !(iflag_trig_bl)
	2873
	2874	!c-----------Statistical closure-----------
	2875	if (iflag_clos_bl.ge.1) then
	2876
	2877	do i=1,klon
	2878	alp_bl(i)=alp_bl_stat(i)
	2879	enddo
	2880
	2881	else
	2882
	2883	alp_bl_stat(:)=0.
	2884
	2885	endif !(iflag_clos_bl)
	2886
	2887	IF (prt_level .GE. 10) THEN
	2888	print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
	2889	ENDIF
	2890
	2891	!ccc fin nrlmd le 10/04/2012
	2892
	2893	! ----------------------------------------------------------------------
	2894	! Transport de la TKE par les panaches thermiques.
	2895	! FH : 2010/02/01
	2896	! if (iflag_pbl.eq.10) then
	2897	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
	2898	! s rg,paprs,pbl_tke)
	2899	! endif
	2900	! ----------------------------------------------------------------------
	2901	!IM/FH: 2011/02/23
	2902	! Couplage Thermiques/Emanuel seulement si T<0
	2903	if (iflag_coupl==2) then
	2904	print*,'Couplage Thermiques/Emanuel seulement si T<0'
	2905	do i=1,klon
	2906	if (t_seri(i,lmax_th(i))>273.) then
	2907	Ale_bl(i)=0.
	2908	endif
	2909	enddo
	2910	endif
	2911
	2912	do i=1,klon
	2913	zmax_th(i)=pphi(i,lmax_th(i))/rg
	2914	enddo
	2915
	2916	endif
	2917
	2918
	2919	!c Ajustement sec
	2920	!c ==============
	2921
	2922	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
	2923	! a partir du sommet des thermiques.
	2924	! Dans le cas contraire, on demarre au niveau 1.
	2925
	2926	if (iflag_thermals.ge.13.or.iflag_thermals.eq.0) then
	2927
	2928	if(iflag_thermals.eq.0) then
	2929	IF(prt_level>9)WRITE(lunout,*)'ajsec'
	2930	limbas(:)=1
	2931	else
	2932	limbas(:)=lmax_th(:)
	2933	endif
	2934
	2935	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
	2936	! pour des test de convergence numerique.
	2937	! Le nouveau ajsec est a priori mieux, meme pour le cas
	2938	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
	2939	! non nulles numeriquement pour des mailles non concernees.
	2940
	2941	if (iflag_thermals.eq.0) then
	2942	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
	2943	& , d_t_ajsb, d_q_ajsb)
	2944	else
	2945	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
	2946	& , d_t_ajsb, d_q_ajsb)
	2947	endif
	2948
	2949	!-----------------------------------------------------------------------------------------
	2950	! ajout des tendances de l'ajustement sec ou des thermiques
	2951	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,'ajsb')
	2952	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
	2953	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
	2954
	2955	!-----------------------------------------------------------------------------------------
	2956
	2957	endif
	2958
	2959	endif
	2960	!c
	2961	!c===================================================================
	2962	!IM
	2963	IF (ip_ebil_phy.ge.2) THEN
	2964	ztit='after dry_adjust'
	2965	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	2966	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2967	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2968	call diagphy(airephy,ztit,ip_ebil_phy &
	2969	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	2970	& , zero_v, zero_v, zero_v, ztsol &
	2971	& , d_h_vcol, d_qt, d_ec &
	2972	& , fs_bound, fq_bound )
	2973	END IF
	2974
	2975
	2976	!c-------------------------------------------------------------------------
	2977	! Computation of ratqs, the width (normalized) of the subrid scale
	2978	! water distribution
	2979	CALL calcratqs(klon,klev,prt_level,lunout, &
	2980	& iflag_ratqs,iflag_con,iflag_cldcon,pdtphys, &
	2981	& ratqsbas,ratqshaut,tau_ratqs,fact_cldcon, &
	2982	& ptconv,ptconvth,clwcon0th, rnebcon0th, &
	2983	& paprs,pplay,q_seri,zqsat,fm_therm, &
	2984	& ratqs,ratqsc)
	2985
	2986
	2987	!c
	2988	!c Appeler le processus de condensation a grande echelle
	2989	!c et le processus de precipitation
	2990	!c-------------------------------------------------------------------------
	2991	IF (prt_level .GE.10) THEN
	2992	print *,' ->fisrtilp '
	2993	ENDIF
	2994	!c-------------------------------------------------------------------------
	2995	CALL fisrtilp(dtime,paprs,pplay, &
	2996	& t_seri, q_seri,ptconv,ratqs, &
	2997	& d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq, &
	2998	& rain_lsc, snow_lsc, &
	2999	& pfrac_impa, pfrac_nucl, pfrac_1nucl, &
	3000	& frac_impa, frac_nucl, beta_prec_fisrt, &
	3001	& prfl, psfl, rhcl, &
	3002	& zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon )
	3003
	3004	WHERE (rain_lsc < 0) rain_lsc = 0.
	3005	WHERE (snow_lsc < 0) snow_lsc = 0.
	3006	!-----------------------------------------------------------------------------------------
	3007	! ajout des tendances de la diffusion turbulente
	3008	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,'lsc')
	3009	!-----------------------------------------------------------------------------------------
	3010	DO k = 1, klev
	3011	DO i = 1, klon
	3012	cldfra(i,k) = rneb(i,k)
	3013	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
	3014	ENDDO
	3015	ENDDO
	3016	IF (check) THEN
	3017	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
	3018	WRITE(lunout,*)"apresilp=", za
	3019	zx_t = 0.0
	3020	za = 0.0
	3021	DO i = 1, klon
	3022	za = za + airephy(i)/REAL(klon)
	3023	zx_t = zx_t + (rain_lsc(i) &
	3024	& + snow_lsc(i))*airephy(i)/REAL(klon)
	3025	ENDDO
	3026	zx_t = zx_t/za*dtime
	3027	WRITE(lunout,*)"Precip=", zx_t
	3028	ENDIF
	3029	!IM
	3030	IF (ip_ebil_phy.ge.2) THEN
	3031	ztit='after fisrt'
	3032	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	3033	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3034	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3035	call diagphy(airephy,ztit,ip_ebil_phy &
	3036	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	3037	& , zero_v, rain_lsc, snow_lsc, ztsol &
	3038	& , d_h_vcol, d_qt, d_ec &
	3039	& , fs_bound, fq_bound )
	3040	END IF
	3041	if (mydebug) then
	3042	call writefield_phy('u_seri',u_seri,llm)
	3043	call writefield_phy('v_seri',v_seri,llm)
	3044	call writefield_phy('t_seri',t_seri,llm)
	3045	call writefield_phy('q_seri',q_seri,llm)
	3046	endif
	3047
	3048	!c
	3049	!c-------------------------------------------------------------------
	3050	!c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
	3051	!c-------------------------------------------------------------------
	3052
	3053	!c 1. NUAGES CONVECTIFS
	3054	!c
	3055	!IM cf FH
	3056	!c IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
	3057	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
	3058	snow_tiedtke=0.
	3059	!c print*,'avant calcul de la pseudo precip '
	3060	!c print*,'iflag_cldcon',iflag_cldcon
	3061	if (iflag_cldcon.eq.-1) then
	3062	rain_tiedtke=rain_con
	3063	else
	3064	!c print*,'calcul de la pseudo precip '
	3065	rain_tiedtke=0.
	3066	!c print*,'calcul de la pseudo precip 0'
	3067	do k=1,klev
	3068	do i=1,klon
	3069	if (d_q_con(i,k).lt.0.) then
	3070	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
	3071	& *(paprs(i,k)-paprs(i,k+1))/rg
	3072	endif
	3073	enddo
	3074	enddo
	3075	endif
	3076	!c
	3077	!c call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
	3078	!c
	3079
	3080	!c Nuages diagnostiques pour Tiedtke
	3081	CALL diagcld1(paprs,pplay, &
	3082	!IM cf FH . rain_con,snow_con,ibas_con,itop_con, &
	3083	& rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
	3084	& diafra,dialiq)
	3085	DO k = 1, klev
	3086	DO i = 1, klon
	3087	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	3088	cldliq(i,k) = dialiq(i,k)
	3089	cldfra(i,k) = diafra(i,k)
	3090	ENDIF
	3091	ENDDO
	3092	ENDDO
	3093
	3094	ELSE IF (iflag_cldcon.ge.3) THEN
	3095	!c On prend pour les nuages convectifs le max du calcul de la
	3096	!c convection et du calcul du pas de temps precedent diminue d'un facteur
	3097	!c facttemps
	3098	facteur = pdtphys *facttemps
	3099	do k=1,klev
	3100	do i=1,klon
	3101	rnebcon(i,k)=rnebcon(i,k)*facteur
	3102	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k)) &
	3103	& then
	3104	rnebcon(i,k)=rnebcon0(i,k)
	3105	clwcon(i,k)=clwcon0(i,k)
	3106	endif
	3107	enddo
	3108	enddo
	3109
	3110	!c
	3111	!cjq - introduce the aerosol direct and first indirect radiative forcings
	3112	!cjq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
	3113	IF (flag_aerosol .gt. 0) THEN
	3114	IF (.NOT. aerosol_couple) &
	3115	& CALL readaerosol_optic( &
	3116	& debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
	3117	& pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3118	& mass_solu_aero, mass_solu_aero_pi, &
	3119	& tau_aero, piz_aero, cg_aero, &
	3120	& tausum_aero, tau3d_aero)
	3121	ELSE
	3122	tausum_aero(:,:,:) = 0.
	3123	tau_aero(:,:,:,:) = 0.
	3124	piz_aero(:,:,:,:) = 0.
	3125	cg_aero(:,:,:,:) = 0.
	3126	ENDIF
	3127	!c
	3128	!c--STRAT AEROSOL
	3129	!c--updates tausum_aero,tau_aero,piz_aero,cg_aero
	3130	IF (flag_aerosol_strat) THEN
	3131	PRINT *,'appel a readaerosolstrat', mth_cur
	3132	CALL readaerosolstrato(debut)
	3133	ENDIF
	3134	!c--fin STRAT AEROSOL
	3135
	3136	!IM calcul nuages par le simulateur ISCCP
	3137	!c
	3138	#ifdef histISCCP
	3139	IF (ok_isccp) THEN
	3140	!c
	3141	!IM lecture invtau, tautab des fichiers formattes
	3142	!c
	3143	IF (debut) THEN
	3144	!$OMP MASTER
	3145	!c
	3146	open(99,file='tautab.formatted', FORM='FORMATTED')
	3147	read(99,'(f30.20)') tautab_omp
	3148	close(99)
	3149	!c
	3150	open(99,file='invtau.formatted',form='FORMATTED')
	3151	read(99,'(i10)') invtau_omp
	3152
	3153	!c print*,'calcul_simulISCCP invtau_omp',invtau_omp
	3154	!c write(6,'(a,8i10)') 'invtau_omp',(invtau_omp(i),i=1,100)
	3155
	3156	close(99)
	3157	!$OMP END MASTER
	3158	!$OMP BARRIER
	3159	tautab=tautab_omp
	3160	invtau=invtau_omp
	3161	!c
	3162	ENDIF !debut
	3163	!c
	3164	!IM appel simulateur toutes les NINT(freq_ISCCP/dtime) heures
	3165	IF (MOD(itap,NINT(freq_ISCCP/dtime)).EQ.0) THEN
	3166	#include "calcul_simulISCCP.h"
	3167	ENDIF !(MOD(itap,NINT(freq_ISCCP/dtime))
	3168	ENDIF !ok_isccp
	3169	#endif
	3170
	3171	!c On prend la somme des fractions nuageuses et des contenus en eau
	3172
	3173	if (iflag_cldcon>=5) then
	3174
	3175	do k=1,klev
	3176	ptconvth(:,k)=fm_therm(:,k+1)>0.
	3177	enddo
	3178
	3179	if (iflag_coupl==4) then
	3180
	3181	! Dans le cas iflag_coupl==4, on prend la somme des convertures
	3182	! convectives et lsc dans la partie des thermiques
	3183	! Le controle par iflag_coupl est peut etre provisoire.
	3184	do k=1,klev
	3185	do i=1,klon
	3186	if (ptconv(i,k).and.ptconvth(i,k)) then
	3187	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
	3188	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
	3189	else if (ptconv(i,k)) then
	3190	cldfra(i,k)=rnebcon(i,k)
	3191	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
	3192	endif
	3193	enddo
	3194	enddo
	3195
	3196	else if (iflag_coupl==5) then
	3197	do k=1,klev
	3198	do i=1,klon
	3199	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
	3200	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
	3201	enddo
	3202	enddo
	3203
	3204	else
	3205
	3206	! Si on est sur un point touche par la convection profonde et pas
	3207	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
	3208	! de la convection profonde.
	3209
	3210	!IM/FH: 2011/02/23
	3211	! definition des points sur lesquels ls thermiques sont actifs
	3212
	3213	do k=1,klev
	3214	do i=1,klon
	3215	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
	3216	cldfra(i,k)=rnebcon(i,k)
	3217	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
	3218	endif
	3219	enddo
	3220	enddo
	3221
	3222	endif
	3223
	3224	else
	3225
	3226	! Ancienne version
	3227	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	3228	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
	3229	endif
	3230
	3231	ENDIF
	3232
	3233	! plulsc(:)=0.
	3234	! do k=1,klev,-1
	3235	! do i=1,klon
	3236	! zzz=prfl(:,k)+psfl(:,k)
	3237	! if (.not.ptconvth.zzz.gt.0.)
	3238	! enddo prfl, psfl,
	3239	! enddo
	3240	!c
	3241	!c 2. NUAGES STARTIFORMES
	3242	!c
	3243	IF (ok_stratus) THEN
	3244	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
	3245	DO k = 1, klev
	3246	DO i = 1, klon
	3247	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	3248	cldliq(i,k) = dialiq(i,k)
	3249	cldfra(i,k) = diafra(i,k)
	3250	ENDIF
	3251	ENDDO
	3252	ENDDO
	3253	ENDIF
	3254	!c
	3255	!c Precipitation totale
	3256	!c
	3257	DO i = 1, klon
	3258	rain_fall(i) = rain_con(i) + rain_lsc(i)
	3259	snow_fall(i) = snow_con(i) + snow_lsc(i)
	3260	ENDDO
	3261	!IM
	3262	IF (ip_ebil_phy.ge.2) THEN
	3263	ztit="after diagcld"
	3264	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	3265	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3266	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3267	call diagphy(airephy,ztit,ip_ebil_phy &
	3268	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	3269	& , zero_v, zero_v, zero_v, ztsol &
	3270	& , d_h_vcol, d_qt, d_ec &
	3271	& , fs_bound, fq_bound )
	3272	END IF
	3273	!c
	3274	!c Calculer l'humidite relative pour diagnostique
	3275	!c
	3276	DO k = 1, klev
	3277	DO i = 1, klon
	3278	zx_t = t_seri(i,k)
	3279	IF (thermcep) THEN
	3280	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	3281	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	3282	zx_qs = MIN(0.5,zx_qs)
	3283	zcor = 1./(1.-retv*zx_qs)
	3284	zx_qs = zx_qs*zcor
	3285	ELSE
	3286	IF (zx_t.LT.t_coup) THEN
	3287	zx_qs = qsats(zx_t)/pplay(i,k)
	3288	ELSE
	3289	zx_qs = qsatl(zx_t)/pplay(i,k)
	3290	ENDIF
	3291	ENDIF
	3292	zx_rh(i,k) = q_seri(i,k)/zx_qs
	3293	zqsat(i,k)=zx_qs
	3294	ENDDO
	3295	ENDDO
	3296
	3297	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
	3298	!c equivalente a 2m (tpote) pour diagnostique
	3299	!c
	3300	DO i = 1, klon
	3301	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
	3302	IF (thermcep) THEN
	3303	IF(zt2m(i).LT.RTT) then
	3304	Lheat=RLSTT
	3305	ELSE
	3306	Lheat=RLVTT
	3307	ENDIF
	3308	ELSE
	3309	IF (zt2m(i).LT.RTT) THEN
	3310	Lheat=RLSTT
	3311	ELSE
	3312	Lheat=RLVTT
	3313	ENDIF
	3314	ENDIF
	3315	tpote(i) = tpot(i)* &
	3316	& EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
	3317	ENDDO
	3318
	3319	IF (type_trac == 'inca') THEN
	3320	#ifdef INCA
	3321	CALL VTe(VTphysiq)
	3322	CALL VTb(VTinca)
	3323	calday = REAL(days_elapsed + 1) + jH_cur
	3324
	3325	call chemtime(itap+itau_phy-1, date0, dtime)
	3326	IF (config_inca == 'aero') THEN
	3327	CALL AEROSOL_METEO_CALC( &
	3328	& calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
	3329	& prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
	3330	END IF
	3331
	3332	zxsnow_dummy(:) = 0.0
	3333
	3334	CALL chemhook_begin (calday, &
	3335	& days_elapsed+1, &
	3336	& jH_cur, &
	3337	& pctsrf(1,1), &
	3338	& rlat, &
	3339	& rlon, &
	3340	& airephy, &
	3341	& paprs, &
	3342	& pplay, &
	3343	& coefh(:,:,is_ave), &
	3344	& pphi, &
	3345	& t_seri, &
	3346	& u, &
	3347	& v, &
	3348	& wo(:, :, 1), &
	3349	& q_seri, &
	3350	& zxtsol, &
	3351	& zxsnow_dummy, &
	3352	& solsw, &
	3353	& albsol1, &
	3354	& rain_fall, &
	3355	& snow_fall, &
	3356	& itop_con, &
	3357	& ibas_con, &
	3358	& cldfra, &
	3359	& iim, &
	3360	& jjm, &
	3361	& tr_seri, &
	3362	& ftsol, &
	3363	& paprs, &
	3364	& cdragh, &
	3365	& cdragm, &
	3366	& pctsrf, &
	3367	& pdtphys, &
	3368	& itap)
	3369
	3370	CALL VTe(VTinca)
	3371	CALL VTb(VTphysiq)
	3372	#endif
	3373	END IF !type_trac = inca
	3374	!c
	3375	!c Calculer les parametres optiques des nuages et quelques
	3376	!c parametres pour diagnostiques:
	3377	!c
	3378
	3379	IF (aerosol_couple) THEN
	3380	mass_solu_aero(:,:) = ccm(:,:,1)
	3381	mass_solu_aero_pi(:,:) = ccm(:,:,2)
	3382	END IF
	3383
	3384	if (ok_newmicro) then
	3385	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
	3386	& paprs, pplay, t_seri, cldliq, cldfra, &
	3387	& cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
	3388	& flwp, fiwp, flwc, fiwc, &
	3389	& mass_solu_aero, mass_solu_aero_pi, &
	3390	& cldtaupi, re, fl, ref_liq, ref_ice)
	3391	else
	3392	CALL nuage (paprs, pplay, &
	3393	& t_seri, cldliq, cldfra, cldtau, cldemi, &
	3394	& cldh, cldl, cldm, cldt, cldq, &
	3395	& ok_aie, &
	3396	& mass_solu_aero, mass_solu_aero_pi, &
	3397	& bl95_b0, bl95_b1, &
	3398	& cldtaupi, re, fl)
	3399	endif
	3400	!c
	3401	!IM betaCRF
	3402	!c
	3403	cldtaurad = cldtau
	3404	cldtaupirad = cldtaupi
	3405	cldemirad = cldemi &
	3406	&!c
	3407	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
	3408	&lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
	3409	!c
	3410	!c global
	3411	!c
	3412	DO k=1, klev
	3413	DO i=1, klon
	3414	if (pplay(i,k).GE.pfree) THEN
	3415	beta(i,k) = beta_pbl
	3416	else
	3417	beta(i,k) = beta_free
	3418	endif
	3419	if (mskocean_beta) THEN
	3420	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
	3421	endif
	3422	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
	3423	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
	3424	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
	3425	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
	3426	ENDDO
	3427	ENDDO
	3428	!c
	3429	else
	3430	!c
	3431	!c regional
	3432	!c
	3433	DO k=1, klev
	3434	DO i=1,klon
	3435	!c
	3436	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
	3437	& rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
	3438	if (pplay(i,k).GE.pfree) THEN
	3439	beta(i,k) = beta_pbl
	3440	else
	3441	beta(i,k) = beta_free
	3442	endif
	3443	if (mskocean_beta) THEN
	3444	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
	3445	endif
	3446	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
	3447	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
	3448	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
	3449	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
	3450	endif
	3451	!c
	3452	ENDDO
	3453	ENDDO
	3454	!c
	3455	endif
	3456	!c
	3457	!c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
	3458	!c
	3459	IF (MOD(itaprad,radpas).EQ.0) THEN
	3460
	3461	DO i = 1, klon
	3462	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce) &
	3463	& + falb1(i,is_lic) * pctsrf(i,is_lic) &
	3464	& + falb1(i,is_ter) * pctsrf(i,is_ter) &
	3465	& + falb1(i,is_sic) * pctsrf(i,is_sic)
	3466	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce) &
	3467	& + falb2(i,is_lic) * pctsrf(i,is_lic) &
	3468	& + falb2(i,is_ter) * pctsrf(i,is_ter) &
	3469	& + falb2(i,is_sic) * pctsrf(i,is_sic)
	3470	ENDDO
	3471
	3472	if (mydebug) then
	3473	call writefield_phy('u_seri',u_seri,llm)
	3474	call writefield_phy('v_seri',v_seri,llm)
	3475	call writefield_phy('t_seri',t_seri,llm)
	3476	call writefield_phy('q_seri',q_seri,llm)
	3477	endif
	3478
	3479	IF (aerosol_couple) THEN
	3480	#ifdef INCA
	3481	CALL radlwsw_inca &
	3482	& (kdlon,kflev,dist, rmu0, fract, solaire, &
	3483	& paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
	3484	& wo(:, :, 1), &
	3485	& cldfrarad, cldemirad, cldtaurad, &
	3486	& heat,heat0,cool,cool0,radsol,albpla, &
	3487	& topsw,toplw,solsw,sollw, &
	3488	& sollwdown, &
	3489	& topsw0,toplw0,solsw0,sollw0, &
	3490	& lwdn0, lwdn, lwup0, lwup, &
	3491	& swdn0, swdn, swup0, swup, &
	3492	& ok_ade, ok_aie, &
	3493	& tau_aero, piz_aero, cg_aero, &
	3494	& topswad_aero, solswad_aero, &
	3495	& topswad0_aero, solswad0_aero, &
	3496	& topsw_aero, topsw0_aero, &
	3497	& solsw_aero, solsw0_aero, &
	3498	& cldtaupirad, &
	3499	& topswai_aero, solswai_aero)
	3500
	3501	#endif
	3502	ELSE
	3503	!c
	3504	!IM calcul radiatif pour le cas actuel
	3505	!c
	3506	RCO2 = RCO2_act
	3507	RCH4 = RCH4_act
	3508	RN2O = RN2O_act
	3509	RCFC11 = RCFC11_act
	3510	RCFC12 = RCFC12_act
	3511	!c
	3512	IF (prt_level .GE.10) THEN
	3513	print *,' ->radlwsw, number 1 '
	3514	ENDIF
	3515	!c
[6]	3516	PRINT *,' Lluis before calling radlwsw rmu0: ',rmu0(llp), &
	3517	' dist: ',dist,' frac: ',fract(llp),' jH_cur: ',jH_cur
[1]	3518	CALL radlwsw &
	3519	& (dist, rmu0, fract, &
	3520	& paprs, pplay,zxtsol,albsol1, albsol2, &
	3521	& t_seri,q_seri,wo, &
	3522	& cldfrarad, cldemirad, cldtaurad, &
	3523	& ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
	3524	& flag_aerosol_strat, &
	3525	& tau_aero, piz_aero, cg_aero, &
	3526	& cldtaupirad,new_aod, &
	3527	& zqsat, flwc, fiwc, &
	3528	& heat,heat0,cool,cool0,radsol,albpla, &
	3529	& topsw,toplw,solsw,sollw, &
	3530	& sollwdown, &
	3531	& topsw0,toplw0,solsw0,sollw0, &
	3532	& lwdn0, lwdn, lwup0, lwup, &
	3533	& swdn0, swdn, swup0, swup, &
	3534	& topswad_aero, solswad_aero, &
	3535	& topswai_aero, solswai_aero, &
	3536	& topswad0_aero, solswad0_aero, &
	3537	& topsw_aero, topsw0_aero, &
	3538	& solsw_aero, solsw0_aero, &
	3539	& topswcf_aero, solswcf_aero)
[6]	3540	PRINT *,' Lluis after calling radlwsw rmu0: ',rmu0(llp), &
	3541	' dist: ',dist,' frac: ',fract(llp),' jH_cur: ',jH_cur
[1]	3542
	3543	!c
	3544	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
	3545	!IM des taux doit etre different du taux actuel
	3546	!IM Par defaut on a les taux perturbes egaux aux taux actuels
	3547	!c
	3548
	3549	if (ok_4xCO2atm) then
	3550	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
	3551	&RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
	3552	&RCFC12_per.NE.RCFC12_act) THEN
	3553	!c
	3554	RCO2 = RCO2_per
	3555	RCH4 = RCH4_per
	3556	RN2O = RN2O_per
	3557	RCFC11 = RCFC11_per
	3558	RCFC12 = RCFC12_per
	3559	!c
	3560	IF (prt_level .GE.10) THEN
	3561	print *,' ->radlwsw, number 2 '
	3562	ENDIF
	3563	!c
	3564	CALL radlwsw &
	3565	& (dist, rmu0, fract, &
	3566	& paprs, pplay,zxtsol,albsol1, albsol2, &
	3567	& t_seri,q_seri,wo, &
	3568	& cldfra, cldemi, cldtau, &
	3569	& ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
	3570	& flag_aerosol_strat, &
	3571	& tau_aero, piz_aero, cg_aero, &
	3572	& cldtaupi,new_aod, &
	3573	& zqsat, flwc, fiwc, &
	3574	& heatp,heat0p,coolp,cool0p,radsolp,albplap, &
	3575	& topswp,toplwp,solswp,sollwp, &
	3576	& sollwdownp, &
	3577	& topsw0p,toplw0p,solsw0p,sollw0p, &
	3578	& lwdn0p, lwdnp, lwup0p, lwupp, &
	3579	& swdn0p, swdnp, swup0p, swupp, &
	3580	& topswad_aerop, solswad_aerop, &
	3581	& topswai_aerop, solswai_aerop, &
	3582	& topswad0_aerop, solswad0_aerop, &
	3583	& topsw_aerop, topsw0_aerop, &
	3584	& solsw_aerop, solsw0_aerop, &
	3585	& topswcf_aerop, solswcf_aerop)
	3586	endif
	3587	endif
	3588	!c
	3589	ENDIF ! aerosol_couple
	3590	itaprad = 0
	3591	ENDIF ! MOD(itaprad,radpas)
	3592	itaprad = itaprad + 1
	3593
	3594	IF (iflag_radia.eq.0) THEN
	3595	IF (prt_level.ge.9) THEN
	3596	PRINT *,'--------------------------------------------------'
	3597	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
	3598	PRINT *,'>>>> heat et cool mis a zero '
	3599	PRINT *,'--------------------------------------------------'
	3600	END IF
	3601	heat=0.
	3602	cool=0.
	3603	sollw=0. ! MPL 01032011
	3604	solsw=0.
	3605	radsol=0.
	3606	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
	3607	swup0=0.
	3608	swdn=0.
	3609	swdn0=0.
	3610	lwup=0.
	3611	lwup0=0.
	3612	lwdn=0.
	3613	lwdn0=0.
	3614	END IF
	3615
	3616	!c
	3617	!c Ajouter la tendance des rayonnements (tous les pas)
	3618	!c
	3619	DO k = 1, klev
	3620	DO i = 1, klon
	3621	t_seri(i,k) = t_seri(i,k) &
	3622	& + (heat(i,k)-cool(i,k)) * dtime/RDAY
	3623	ENDDO
	3624	ENDDO
	3625	!c
	3626	if (mydebug) then
	3627	call writefield_phy('u_seri',u_seri,llm)
	3628	call writefield_phy('v_seri',v_seri,llm)
	3629	call writefield_phy('t_seri',t_seri,llm)
	3630	call writefield_phy('q_seri',q_seri,llm)
	3631	endif
	3632
	3633	!IM
	3634	IF (ip_ebil_phy.ge.2) THEN
	3635	ztit='after rad'
	3636	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	3637	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3638	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3639	call diagphy(airephy,ztit,ip_ebil_phy &
	3640	& , topsw, toplw, solsw, sollw, zero_v &
	3641	& , zero_v, zero_v, zero_v, ztsol &
	3642	& , d_h_vcol, d_qt, d_ec &
	3643	& , fs_bound, fq_bound )
	3644	END IF
	3645	!c
	3646	!c
	3647	!c Calculer l'hydrologie de la surface
	3648	!c
	3649	!c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
	3650	!c . agesno, ftsol,fqsurf,fsnow, ruis)
	3651	!c
	3652
	3653	!c
	3654	!c Calculer le bilan du sol et la derive de temperature (couplage)
	3655	!c
	3656	DO i = 1, klon
	3657	!c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
	3658	!c a la demande de JLD
	3659	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
	3660	ENDDO
	3661	!c
	3662	!cmoddeblott(jan95)
	3663	!c Appeler le programme de parametrisation de l'orographie
	3664	!c a l'echelle sous-maille:
	3665	!c
	3666	IF (prt_level .GE.10) THEN
	3667	print *,' call orography ? ', ok_orodr
	3668	ENDIF
	3669	!c
	3670	IF (ok_orodr) THEN
	3671	!c
	3672	!c selection des points pour lesquels le shema est actif:
	3673	igwd=0
	3674	DO i=1,klon
	3675	itest(i)=0
	3676	!c IF ((zstd(i).gt.10.0)) THEN
	3677	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	3678	itest(i)=1
	3679	igwd=igwd+1
	3680	idx(igwd)=i
	3681	ENDIF
	3682	ENDDO
	3683	!c igwdim=MAX(1,igwd)
	3684	!c
	3685	IF (ok_strato) THEN
	3686
	3687	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
	3688	& zmea,zstd, zsig, zgam, zthe,zpic,zval, &
	3689	& igwd,idx,itest, &
	3690	& t_seri, u_seri, v_seri, &
	3691	& zulow, zvlow, zustrdr, zvstrdr, &
	3692	& d_t_oro, d_u_oro, d_v_oro)
	3693
	3694	ELSE
	3695	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
	3696	& zmea,zstd, zsig, zgam, zthe,zpic,zval, &
	3697	& igwd,idx,itest, &
	3698	& t_seri, u_seri, v_seri, &
	3699	& zulow, zvlow, zustrdr, zvstrdr, &
	3700	& d_t_oro, d_u_oro, d_v_oro)
	3701	ENDIF
	3702	!c
	3703	!c ajout des tendances
	3704	!-----------------------------------------------------------------------------------------
	3705	! ajout des tendances de la trainee de l'orographie
	3706	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,'oro')
	3707	!-----------------------------------------------------------------------------------------
	3708	!c
	3709	ENDIF ! fin de test sur ok_orodr
	3710	!c
	3711	if (mydebug) then
	3712	call writefield_phy('u_seri',u_seri,llm)
	3713	call writefield_phy('v_seri',v_seri,llm)
	3714	call writefield_phy('t_seri',t_seri,llm)
	3715	call writefield_phy('q_seri',q_seri,llm)
	3716	endif
	3717
	3718	IF (ok_orolf) THEN
	3719	!c
	3720	!c selection des points pour lesquels le shema est actif:
	3721	igwd=0
	3722	DO i=1,klon
	3723	itest(i)=0
	3724	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	3725	itest(i)=1
	3726	igwd=igwd+1
	3727	idx(igwd)=i
	3728	ENDIF
	3729	ENDDO
	3730	!c igwdim=MAX(1,igwd)
	3731	!c
	3732	IF (ok_strato) THEN
	3733
	3734	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
	3735	& rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
	3736	& igwd,idx,itest, &
	3737	& t_seri, u_seri, v_seri, &
	3738	& zulow, zvlow, zustrli, zvstrli, &
	3739	& d_t_lif, d_u_lif, d_v_lif )
	3740
	3741	ELSE
	3742	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
	3743	& rlat,zmea,zstd,zpic, &
	3744	& itest, &
	3745	& t_seri, u_seri, v_seri, &
	3746	& zulow, zvlow, zustrli, zvstrli, &
	3747	& d_t_lif, d_u_lif, d_v_lif)
	3748	ENDIF
	3749	!c
	3750	!-----------------------------------------------------------------------------------------
	3751	! ajout des tendances de la portance de l'orographie
	3752	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,'lif')
	3753	!-----------------------------------------------------------------------------------------
	3754	!c
	3755	ENDIF ! fin de test sur ok_orolf
	3756	!C HINES GWD PARAMETRIZATION
	3757
	3758	IF (ok_hines) then
	3759
	3760	CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
	3761	& rlat,t_seri,u_seri,v_seri, &
	3762	& zustrhi,zvstrhi, &
	3763	& d_t_hin, d_u_hin, d_v_hin)
	3764	!c
	3765	!c ajout des tendances
	3766	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,'hin')
	3767
	3768	ENDIF
	3769	!c
	3770
	3771	!c
	3772	!IM cf. FLott BEG
	3773	!C STRESS NECESSAIRES: TOUTE LA PHYSIQUE
	3774
	3775	if (mydebug) then
	3776	call writefield_phy('u_seri',u_seri,llm)
	3777	call writefield_phy('v_seri',v_seri,llm)
	3778	call writefield_phy('t_seri',t_seri,llm)
	3779	call writefield_phy('q_seri',q_seri,llm)
	3780	endif
	3781
	3782	DO i = 1, klon
	3783	zustrph(i)=0.
	3784	zvstrph(i)=0.
	3785	ENDDO
	3786	DO k = 1, klev
	3787	DO i = 1, klon
	3788	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
	3789	& (paprs(i,k)-paprs(i,k+1))/rg
	3790	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
	3791	& (paprs(i,k)-paprs(i,k+1))/rg
	3792	ENDDO
	3793	ENDDO
	3794	!c
	3795	!IM calcul composantes axiales du moment angulaire et couple des montagnes
	3796	!c
	3797	! IF (is_sequential .and. ok_orodr) THEN
	3798	! CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
	3799	! & ra,rg,romega, &
	3800	! & rlat,rlon,pphis, &
	3801	! & zustrdr,zustrli,zustrph, &
	3802	! & zvstrdr,zvstrli,zvstrph, &
	3803	! & paprs,u,v, &
	3804	! & aam, torsfc)
	3805	! ENDIF
	3806	!IM cf. FLott END
	3807	!IM
	3808	IF (ip_ebil_phy.ge.2) THEN
	3809	ztit='after orography'
	3810	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
	3811	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3812	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3813	call diagphy(airephy,ztit,ip_ebil_phy &
	3814	& , zero_v, zero_v, zero_v, zero_v, zero_v &
	3815	& , zero_v, zero_v, zero_v, ztsol &
	3816	& , d_h_vcol, d_qt, d_ec &
	3817	& , fs_bound, fq_bound )
	3818	END IF
	3819	!c
	3820	!c
	3821	!====================================================================
	3822	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
	3823	!====================================================================
	3824	! Abderrahmane 24.08.09
	3825
	3826	IF (ok_cosp) THEN
	3827	! adeclarer
	3828	#ifdef CPP_COSP
	3829	IF (MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
	3830
	3831	print*,'freq_cosp',freq_cosp
	3832	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
	3833	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
	3834	! s ref_liq,ref_ice
	3835	call phys_cosp(itap,dtime,freq_cosp, &
	3836	& ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
	3837	& ecrit_mth,ecrit_day,ecrit_hf, &
	3838	& klon,klev,rlon,rlat,presnivs,overlap, &
	3839	& ref_liq,ref_ice, &
	3840	& pctsrf(:,is_ter)+pctsrf(:,is_lic), &
	3841	& zu10m,zv10m,pphis, &
	3842	& zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
	3843	& qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
	3844	& prfl(:,1:klev),psfl(:,1:klev), &
	3845	& pmflxr(:,1:klev),pmflxs(:,1:klev), &
	3846	& mr_ozone,cldtau, cldemi)
	3847
	3848	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
	3849	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
	3850	! M clMISR,
	3851	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
	3852	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
	3853
	3854	ENDIF
	3855
	3856	#endif
	3857	ENDIF !ok_cosp
	3858	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	3859	!cAA
	3860	!cAA Installation de l'interface online-offline pour traceurs
	3861	!cAA
	3862	!c====================================================================
	3863	!c Calcul des tendances traceurs
	3864	!c====================================================================
	3865	!C
	3866
	3867	IF (type_trac=='repr') THEN
	3868	sh_in(:,:) = q_seri(:,:)
	3869	ELSE
	3870	sh_in(:,:) = qx(:,:,ivap)
	3871	END IF
	3872
	3873	call phytrac ( &
	3874	& itap, days_elapsed+1, jH_cur, debut, &
	3875	& lafin, dtime, u, v, t, &
	3876	& paprs, pplay, pmfu, pmfd, &
	3877	& pen_u, pde_u, pen_d, pde_d, &
	3878	& cdragh, coefh(:,:,is_ave), fm_therm, entr_therm, &
	3879	& u1, v1, ftsol, pctsrf, &
	3880	& ustar, u10m, v10m, &
	3881	& rlat, rlon, &
	3882	& frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
	3883	& presnivs, pphis, pphi, albsol1, &
	3884	& sh_in, rhcl, cldfra, rneb, &
	3885	& diafra, cldliq, itop_con, ibas_con, &
	3886	& pmflxr, pmflxs, prfl, psfl, &
	3887	& da, phi, mp, upwd, &
	3888	& phi2, d1a, dam, sij, &
	3889	& wdtrainA, wdtrainM, sigd, clw,elij, &
	3890	& ev, ep, epmlmMm, eplaMm, &
	3891	& dnwd, aerosol_couple, flxmass_w, &
	3892	& tau_aero, piz_aero, cg_aero, ccm, &
	3893	& rfname, &
	3894	& d_tr_dyn, &
	3895	& tr_seri)
	3896
	3897	IF (offline) THEN
	3898
	3899	IF (prt_level.ge.9) &
	3900	& print*,'Attention on met a 0 les thermiques pour phystoke'
	3901	call phystokenc ( &
	3902	& nlon,klev,pdtphys,rlon,rlat, &
	3903	& t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
	3904	& fm_therm,entr_therm, &
	3905	& cdragh,coefh(:,:,is_ave),u1,v1,ftsol,pctsrf, &
	3906	& frac_impa, frac_nucl, &
	3907	& pphis,airephy,dtime,itap, &
	3908	& qx(:,:,ivap),da,phi,mp,upwd,dnwd)
	3909
	3910
	3911	ENDIF
	3912
	3913	!c
	3914	!c Calculer le transport de l'eau et de l'energie (diagnostique)
	3915	!c
	3916	CALL transp (paprs,zxtsol, &
	3917	& t_seri, q_seri, u_seri, v_seri, zphi, &
	3918	& ve, vq, ue, uq)
	3919	!c
	3920	!IM global posePB BEG
	3921	IF(1.EQ.0) THEN
	3922	!c
	3923	CALL transp_lay (paprs,zxtsol, &
	3924	& t_seri, q_seri, u_seri, v_seri, zphi, &
	3925	& ve_lay, vq_lay, ue_lay, uq_lay)
	3926	!c
	3927	ENDIF !(1.EQ.0) THEN
	3928	!IM global posePB END
	3929	!c Accumuler les variables a stocker dans les fichiers histoire:
	3930	!c
	3931
	3932	!================================================================
	3933	! Conversion of kinetic and potential energy into heat, for
	3934	! parameterisation of subgrid-scale motions
	3935	!================================================================
	3936
	3937	d_t_ec(:,:)=0.
	3938	forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
	3939	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
	3940	& u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
	3941	& zmasse,exner,d_t_ec)
	3942	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
	3943
	3944	!IM
	3945	IF (ip_ebil_phy.ge.1) THEN
	3946	ztit='after physic'
	3947	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
	3948	& , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3949	& , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3950	!C Comme les tendances de la physique sont ajoute dans la dynamique,
	3951	!C on devrait avoir que la variation d'entalpie par la dynamique
	3952	!C est egale a la variation de la physique au pas de temps precedent.
	3953	!C Donc la somme de ces 2 variations devrait etre nulle.
	3954
	3955	call diagphy(airephy,ztit,ip_ebil_phy &
	3956	& , topsw, toplw, solsw, sollw, sens &
	3957	& , evap, rain_fall, snow_fall, ztsol &
	3958	& , d_h_vcol, d_qt, d_ec &
	3959	& , fs_bound, fq_bound )
	3960	!C
	3961	d_h_vcol_phy=d_h_vcol
	3962	!C
	3963	END IF
	3964	PRINT *,'Lluis Reaching the SORTIES point'
	3965
	3966	!C
	3967	!c=======================================================================
	3968	!c SORTIES
	3969	!c=======================================================================
	3970
	3971	!IM Interpolation sur les niveaux de pression du NMC
	3972	!c -------------------------------------------------
	3973	!c
	3974	#include "calcul_STDlev.h"
	3975	!c
	3976	!c slp sea level pressure
	3977	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
	3978	!c
	3979	!ccc prw = eau precipitable
	3980	DO i = 1, klon
	3981	prw(i) = 0.
	3982	DO k = 1, klev
	3983	prw(i) = prw(i) + &
	3984	& q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
	3985	ENDDO
	3986	ENDDO
	3987	!c
	3988	!IM initialisation + calculs divers diag AMIP2
	3989	!c
	3990	#include "calcul_divers.h"
	3991	!c
	3992	IF (type_trac == 'inca') THEN
	3993	#ifdef INCA
	3994	CALL VTe(VTphysiq)
	3995	CALL VTb(VTinca)
	3996
	3997	CALL chemhook_end ( &
	3998	& dtime, &
	3999	& pplay, &
	4000	& t_seri, &
	4001	& tr_seri, &
	4002	& nbtr, &
	4003	& paprs, &
	4004	& q_seri, &
	4005	& airephy, &
	4006	& pphi, &
	4007	& pphis, &
	4008	& zx_rh)
	4009
	4010	CALL VTe(VTinca)
	4011	CALL VTb(VTphysiq)
	4012	#endif
	4013	END IF
	4014
	4015
	4016	!c
	4017	!c Convertir les incrementations en tendances
	4018	!c
	4019	IF (prt_level .GE.10) THEN
	4020	print *,'Convertir les incrementations en tendances '
	4021	ENDIF
	4022	!c
	4023	if (mydebug) then
	4024	call writefield_phy('u_seri',u_seri,llm)
	4025	call writefield_phy('v_seri',v_seri,llm)
	4026	call writefield_phy('t_seri',t_seri,llm)
	4027	call writefield_phy('q_seri',q_seri,llm)
	4028	endif
	4029
	4030	DO k = 1, klev
	4031	DO i = 1, klon
	4032	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	4033	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	4034	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
	4035	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
	4036	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
	4037	ENDDO
	4038	ENDDO
	4039	!c
	4040	IF (nqtot.GE.3) THEN
	4041	DO iq = 3, nqtot
	4042	DO k = 1, klev
	4043	DO i = 1, klon
	4044	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
	4045	ENDDO
	4046	ENDDO
	4047	ENDDO
	4048	ENDIF
	4049	!c
	4050	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
	4051	!IM global posePB#include "write_bilKP_ins.h"
	4052	!IM global posePB#include "write_bilKP_ave.h"
	4053	!c
	4054
	4055	!c Sauvegarder les valeurs de t et q a la fin de la physique:
	4056	!c
	4057	DO k = 1, klev
	4058	DO i = 1, klon
	4059	u_ancien(i,k) = u_seri(i,k)
	4060	v_ancien(i,k) = v_seri(i,k)
	4061	t_ancien(i,k) = t_seri(i,k)
	4062	q_ancien(i,k) = q_seri(i,k)
	4063	ENDDO
	4064	ENDDO
	4065
	4066	!!! RomP >>>
	4067	IF (nqtot.GE.3) THEN
	4068	DO iq = 3, nqtot
	4069	DO k = 1, klev
	4070	DO i = 1, klon
	4071	tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
	4072	ENDDO
	4073	ENDDO
	4074	ENDDO
	4075	ENDIF
	4076	!!! RomP <<<
	4077	!==========================================================================
	4078	! Sorties des tendances pour un point particulier
	4079	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
	4080	! pour le debug
	4081	! La valeur de igout est attribuee plus haut dans le programme
	4082	!==========================================================================
	4083
	4084	if (prt_level.ge.1) then
	4085	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	4086	write(lunout,*) &
	4087	& 'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
	4088	write(lunout,*) &
	4089	& nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
	4090	& pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
	4091	& pctsrf(igout,is_sic)
	4092	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
	4093	do k=1,klev
	4094	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
	4095	& d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
	4096	& d_t_eva(igout,k)
	4097	enddo
	4098	write(lunout,*) 'cool,heat'
	4099	do k=1,klev
	4100	write(lunout,*) cool(igout,k),heat(igout,k)
	4101	enddo
	4102
	4103	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
	4104	do k=1,klev
	4105	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
	4106	& d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
	4107	enddo
	4108
	4109	write(lunout,*) 'd_ps ',d_ps(igout)
	4110	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
	4111	do k=1,klev,klon
	4112	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
	4113	& d_qx(igout,k,1),d_qx(igout,k,2)
	4114	enddo
	4115	endif
	4116
	4117	!==========================================================================
	4118
	4119	!c============================================================
	4120	!c Calcul de la temperature potentielle
	4121	!c============================================================
	4122	DO k = 1, klev
	4123	DO i = 1, klon
	4124	!cJYG/IM theta en debut du pas de temps
	4125	!cJYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
	4126	!cJYG/IM theta en fin de pas de temps de physique
	4127	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
	4128	!c thetal: 2 lignes suivantes a decommenter si vous avez les fichiers MPL 20130625
	4129	!c fth_fonctions.F90 et parkind1.F90
	4130	!c sinon thetal=theta
	4131	!c thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
	4132	!c : ql_seri(i,k))
	4133	thetal(i,k)=theta(i,k)
	4134	ENDDO
	4135	ENDDO
	4136	!c
	4137
	4138	!c 22.03.04 BEG
	4139	!c=============================================================
	4140	!c Ecriture des sorties
	4141	!c=============================================================
	4142	#ifdef CPP_IOIPSL
	4143
	4144	!c Recupere des varibles calcule dans differents modules
	4145	!c pour ecriture dans histxxx.nc
	4146
	4147	! Get some variables from module fonte_neige_mod
	4148	!L. Fita, LMD. November 2013
	4149	!! It is not working after removing starphy.nc and limit.nc?
	4150	CALL fonte_neige_get_vars(pctsrf, &
	4151	& zxfqcalving, zxfqfonte, zxffonte)
	4152
	4153
	4154
	4155	!c=============================================================
	4156	! Separation entre thermiques et non thermiques dans les sorties
	4157	! de fisrtilp
	4158	!c=============================================================
	4159
	4160	if (iflag_thermals>=1) then
	4161	d_t_lscth=0.
	4162	d_t_lscst=0.
	4163	d_q_lscth=0.
	4164	d_q_lscst=0.
	4165	do k=1,klev
	4166	do i=1,klon
	4167	if (ptconvth(i,k)) then
	4168	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
	4169	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
	4170	else
	4171	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
	4172	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
	4173	endif
	4174	enddo
	4175	enddo
	4176
	4177	do i=1,klon
	4178	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
	4179	plul_th(i)=prfl(i,1)+psfl(i,1)
	4180	enddo
	4181	endif
	4182
	4183	PRINT *,' Lluis WRITING outputs! itap: ',itap,' itau_phy: ',itau_phy, &
	4184	' itau_w: ',itau_w
	4185	PRINT *,' Lluis writting outputs qsol: ',qsol(llp), &
	4186	' ftsol: ',ftsol(llp,:)
	4187
	4188	#include "phys_output_write_new.h"
	4189
	4190
	4191	#ifdef histISCCP
	4192	#include "write_histISCCP.h"
	4193	#endif
	4194
	4195	PRINT *,' Lluis WRITING histfiles!'
	4196
	4197	#ifdef histNMC
	4198	#include "write_histhfNMC.h"
	4199	#include "write_histdayNMC.h"
	4200	#include "write_histmthNMC.h"
	4201	#endif
	4202
	4203	#include "write_histday_seri.h"
	4204
	4205	#include "write_paramLMDZ_phy.h"
	4206
	4207	#endif
	4208
	4209	!c 22.03.04 END
	4210	!c
	4211	!c====================================================================
	4212	!c Si c'est la fin, il faut conserver l'etat de redemarrage
	4213	!c====================================================================
	4214	!c
	4215
	4216	!c -----------------------------------------------------------------
	4217	!c WSTATS: Saving statistics
	4218	!c -----------------------------------------------------------------
	4219	!c ("stats" stores and accumulates 8 key variables in file "stats.nc"
	4220	!c which can later be used to make the statistic files of the run:
	4221	!c "stats") only possible in 3D runs !
	4222
	4223
	4224	IF (callstats) THEN
	4225
	4226	call wstats(klon,o_psol%name,"Surface pressure","Pa" &
	4227	& ,2,paprs(:,1))
	4228	call wstats(klon,o_tsol%name,"Surface temperature","K", &
	4229	& 2,zxtsol)
	4230	zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
	4231	call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
	4232	& "kg/(s*m2)",2,zx_tmp_fi2d)
	4233	zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
	4234	call wstats(klon,o_plul%name,"Large-scale Precip", &
	4235	& "kg/(s*m2)",2,zx_tmp_fi2d)
	4236	zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
	4237	call wstats(klon,o_pluc%name,"Convective Precip", &
	4238	& "kg/(s*m2)",2,zx_tmp_fi2d)
	4239	call wstats(klon,o_sols%name,"Solar rad. at surf.", &
	4240	& "W/m2",2,solsw)
	4241	call wstats(klon,o_soll%name,"IR rad. at surf.", &
	4242	& "W/m2",2,sollw)
	4243	zx_tmp_fi2d(:) = topsw(:)-toplw(:)
	4244	call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
	4245	& "W/m2",2,zx_tmp_fi2d)
	4246
	4247
	4248
	4249	call wstats(klon,o_temp%name,"Air temperature","K", &
	4250	& 3,t_seri)
	4251	call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
	4252	& 3,u_seri)
	4253	call wstats(klon,o_vitv%name,"Meridional wind", &
	4254	& "m.s-1",3,v_seri)
	4255	call wstats(klon,o_vitw%name,"Vertical wind", &
	4256	& "m.s-1",3,omega)
	4257	call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
	4258	& 3,q_seri)
	4259
	4260
	4261
	4262	IF(lafin) THEN
	4263	write (,) "Writing stats..."
	4264	call mkstats(ierr)
	4265	ENDIF
	4266
	4267	ENDIF !if callstats
	4268
	4269	IF (lafin) THEN
	4270	itau_phy = itau_phy + itap
	4271	CALL phyredem ("restartphy.nc")
	4272	! open(97,form="unformatted",file="finbin")
	4273	! write(97) u_seri,v_seri,t_seri,q_seri
	4274	! close(97)
	4275	!$OMP MASTER
	4276	if (read_climoz >= 1) then
	4277	if (is_mpi_root) then
	4278	call nf95_close(ncid_climoz)
	4279	end if
	4280	deallocate(press_climoz) ! pointer
	4281	end if
	4282	!$OMP END MASTER
	4283	ENDIF
	4284
	4285	! first=.false.
	4286
	4287	! Lluis
	4288	PRINT *,' Lluis: ',klev,' UBOUNDS: ',UBOUND(t_seri), UBOUND(u_seri), &
	4289	UBOUND(d_q_con), UBOUND(d_t_con)
	4290	PRINT *,' Lluis llp ',llp,' itap: ',itap,' zlev t_seri u_seri d_q_con d_t_con_____'
	4291	DO i=1,klev
	4292	PRINT *,i,t_seri(llp,i), u_seri(llp,i), d_q_con(llp,i), d_t_con(llp,i)
	4293	END DO
	4294
	4295
	4296	RETURN
	4297	END SUBROUTINE physiq
	4298
	4299	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
	4300	IMPLICIT none
	4301	!c
	4302	!c Calculer et imprimer l'eau totale. A utiliser pour verifier
	4303	!c la conservation de l'eau
	4304	!c
	4305	#include "YOMCST.h"
	4306	INTEGER klon,klev
	4307	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
	4308	REAL aire(klon)
	4309	REAL qtotal, zx, qcheck
	4310	INTEGER i, k
	4311	!c
	4312	zx = 0.0
	4313	DO i = 1, klon
	4314	zx = zx + aire(i)
	4315	ENDDO
	4316	qtotal = 0.0
	4317	DO k = 1, klev
	4318	DO i = 1, klon
	4319	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
	4320	& *(paprs(i,k)-paprs(i,k+1))/RG
	4321	ENDDO
	4322	ENDDO
	4323	!c
	4324	qcheck = qtotal/zx
	4325	!c
	4326	RETURN
	4327	END
	4328
	4329	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
	4330	IMPLICIT none
	4331	!c
	4332	!c Tranformer une variable de la grille physique a
	4333	!c la grille d'ecriture
	4334	!c
	4335	INTEGER nfield,nlon,iim,jjmp1, jjm
	4336	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
	4337	!c
	4338	INTEGER i, n, ig
	4339	!c
	4340	jjm = jjmp1 - 1
	4341	DO n = 1, nfield
	4342	DO i=1,iim
	4343	ecrit(i,n) = fi(1,n)
	4344	ecrit(i+jjm*iim,n) = fi(nlon,n)
	4345	ENDDO
	4346	DO ig = 1, nlon - 2
	4347	ecrit(iim+ig,n) = fi(1+ig,n)
	4348	ENDDO
	4349	ENDDO
	4350	RETURN
	4351	END SUBROUTINE gr_fi_ecrit

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

Download in other formats:

Original Format