Context Navigation

physiq.F @ 5017

Last change on this file since 5017 was 1399, checked in by idelkadi, 14 years ago

Modifications relatives a l'inclusion des nouveaux nuages d'Arnaud Jam

calltherm.F90 : remontee de variables suplementaires.

Passage de ztv, zpspsk, zthla, zthl en argument

fisrtilp.F : le programme de nuages avec appel a cloudth

Ajout des arguments : zqta, fraca, ztv, zpspsk, ztla, ztlh,
iflag_cldcon

physiq.F : appel modifie a calltherm et fisrtilp
thermcell_main.F90 : changement dans les flag_thermals_ed
thermcell_plume.F90 : version modifiee entrainement et detrainement

(on retrouve la precendente pour iflag_thermals_ed=10)

cloudth.F90 : PDF d'eau bi-gaussiennes

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

File size: 123.6 KB

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

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

Download in other formats:

Original Format