Context Navigation

physiq.F @ 1509

Last change on this file since 1509 was 1507, checked in by idelkadi, 13 years ago

Introduction de nouvelles facon de calucler les ratqs pour firstilp quand iflag_cldcon>=7 : ratqs devient pronos

tic avec une relaxation vers un ratqs issu de la convection (ratqsc) quand celle-ci est active et une relaxation ve
rs ratqss dans le cas contraire.
Les constantes de temps sont a 3h en dur dans le code (en attendant mieux)
On met un coeff iflag_cldcon/100 devant ratqsc
On plafone le ratqs a 0.5 (la aussi en dur dans le code en attendant de stabiliser une facon de faire).

Introduction de nouveaux diagnostics concernant la separation entre la contribution des thermiques et le reste d

ans fisrtilp (dq/tlscth, dq/tlscst,plulth, plulst)
Correction de la sortie de la hauteur des thermiques zmax_th

Introduction d'une option iflag_wake=2 pour prendre en compte les precipitations issues de fisrtipl dans le forc

age des wake.

Correction d'un bug sur l'identification des points affectes par les thermiques (ptconvth(:,k)=fm_therm(:,k+1)>0

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

File size: 128.5 KB

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

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

Download in other formats:

Original Format