Context Navigation

physiq.F @ 1864

Last change on this file since 1864 was 1864, checked in by Laurent Fairhead, 11 years ago

Création d'une nouvelle testing:

merge des modifications du trunk entre r1796 et r1860

New testing version

merged modifications between r1796 and r1860 from the trunk

i.e.
svn merge -r1796:1860 http://svn.lmd.jussieu.fr/LMDZ/LMDZ5/trunk

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

File size: 133.1 KB

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

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

Download in other formats:

Original Format