Context Navigation

physiq.F @ 1267

Last change on this file since 1267 was 1266, checked in by Laurent Fairhead, 15 years ago
Correction d'un bug dans aeropt_5wv.F90 : une variable était utilisée sans être calculée Ajout de tests sur l'écriture de variables dans les fichiers de sortie Changement des arguments d'appels à INCA ACo
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 117.4 KB

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

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

Download in other formats:

Original Format