Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

physiq.F @ 3990

Last change on this file since 3990 was 1634, checked in by Laurent Fairhead, 12 years ago

Version du code LMDZ utilisée dans la configuration IPSLCM5B
de l'exercice CMIP5/AR5

Actual version of the code used for the IPSLCM5B configuration
of the CMIP5/AR5 exercise

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

File size: 136.2 KB

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

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

Download in other formats:

Original Format