Context Navigation

physiq_mod.F90 @ 5456

Last change on this file since 5456 was 3629, checked in by acozic, 5 years ago

Add new grid, new axis and new variables for cmip protocole and dr2xml

field_group id="coord_hyb"
grid_ref="klevp1_bnds"
grid_ref="klev_bnds"
domain id="greordered"
axis id="axis_lat"
axis id="bnds"
axis id="klevp1"
axis id="klev"

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 161.4 KB

Rev	Line
[2418]	1	!
[1279]	2	! $Id: physiq_mod.F90 3629 2020-02-10 08:54:26Z aborella $
[2418]	3	!
[1862]	4	!#define IO_DEBUG
[2418]	5	MODULE physiq_mod
[766]	6
[2469]	7	IMPLICIT NONE
[2418]	8
	9	CONTAINS
	10
[2469]	11	SUBROUTINE physiq (nlon,nlev, &
	12	debut,lafin,pdtphys_, &
	13	paprs,pplay,pphi,pphis,presnivs, &
	14	u,v,rot,t,qx, &
	15	flxmass_w, &
	16	d_u, d_v, d_t, d_qx, d_ps)
[524]	17
[2469]	18	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
	19	histwrite, ju2ymds, ymds2ju, getin
	20	USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
	21	USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
	22	year_cur, mth_cur,jD_cur, jH_cur, jD_ref
	23	USE write_field_phy
	24	USE dimphy
	25	USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac
	26	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo
	27	USE mod_phys_lmdz_para
	28	USE iophy
	29	USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
	30	USE phystokenc_mod, ONLY: offline, phystokenc
	31	USE time_phylmdz_mod, only: raz_date, day_step_phy, update_time
	32	USE vampir
	33	USE pbl_surface_mod, ONLY : pbl_surface
	34	USE change_srf_frac_mod
	35	USE surface_data, ONLY : type_ocean, ok_veget, ok_snow
[2606]	36	USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
	37	! [Variables internes non sauvegardees de la physique]
	38	! Variables locales pour effectuer les appels en serie
	39	t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri, &
	40	! Dynamic tendencies (diagnostics)
	41	d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn, &
	42	d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
	43	! Physic tendencies
	44	d_t_con,d_q_con,d_u_con,d_v_con, &
	45	d_tr, & !! to be removed?? (jyg)
	46	d_t_wake,d_q_wake, &
	47	d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
	48	d_t_ajsb,d_q_ajsb, &
	49	d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
	50	d_t_ajs_w,d_q_ajs_w, &
	51	d_t_ajs_x,d_q_ajs_x, &
	52	!
	53	d_t_eva,d_q_eva, &
	54	d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
	55	d_t_lscst,d_q_lscst, &
	56	d_t_lscth,d_q_lscth, &
	57	plul_st,plul_th, &
	58	!
	59	d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
	60	d_t_vdf_w,d_q_vdf_w, &
	61	d_t_vdf_x,d_q_vdf_x, &
	62	d_ts, &
	63	!
	64	d_t_oli,d_u_oli,d_v_oli, &
	65	d_t_oro,d_u_oro,d_v_oro, &
	66	d_t_lif,d_u_lif,d_v_lif, &
	67	d_t_ec, &
	68	!
	69	du_gwd_hines,dv_gwd_hines,d_t_hin, &
	70	dv_gwd_rando,dv_gwd_front, &
	71	east_gwstress,west_gwstress, &
	72	d_q_ch4, &
	73	! Special RRTM
	74	ZLWFT0_i,ZSWFT0_i,ZFLDN0, &
	75	ZFLUP0,ZFSDN0,ZFSUP0, &
	76	!
	77	topswad_aero,solswad_aero, &
	78	topswai_aero,solswai_aero, &
	79	topswad0_aero,solswad0_aero, &
	80	!LW additional
	81	toplwad_aero,sollwad_aero, &
	82	toplwai_aero,sollwai_aero, &
	83	toplwad0_aero,sollwad0_aero, &
	84	!
	85	topsw_aero,solsw_aero, &
	86	topsw0_aero,solsw0_aero, &
	87	topswcf_aero,solswcf_aero, &
	88	tausum_aero,tau3d_aero, &
	89	!
	90	!variables CFMIP2/CMIP5
	91	topswad_aerop, solswad_aerop, &
	92	topswai_aerop, solswai_aerop, &
	93	topswad0_aerop, solswad0_aerop, &
	94	topsw_aerop, topsw0_aerop, &
	95	solsw_aerop, solsw0_aerop, &
	96	topswcf_aerop, solswcf_aerop, &
	97	!LW diagnostics
	98	toplwad_aerop, sollwad_aerop, &
	99	toplwai_aerop, sollwai_aerop, &
	100	toplwad0_aerop, sollwad0_aerop, &
	101	!
	102	ptstar, pt0, slp, &
	103	!
	104	bils, &
	105	!
	106	cldh, cldl,cldm, cldq, cldt, &
	107	JrNt, &
	108	dthmin, evap, fder, plcl, plfc, &
	109	prw, prlw, prsw, &
	110	s_lcl, s_pblh, s_pblt, s_therm, &
	111	cdragm, cdragh, &
	112	zustar, zu10m, zv10m, rh2m, qsat2m, &
	113	zq2m, zt2m, weak_inversion, &
	114	zt2m_min_mon, zt2m_max_mon, & ! pour calcul_divers.h
	115	t2m_min_mon, t2m_max_mon, & ! pour calcul_divers.h
	116	!
	117	s_pblh_x, s_pblh_w, &
	118	s_lcl_x, s_lcl_w, &
	119	!
	120	slab_wfbils, tpot, tpote, &
	121	ue, uq, ve, vq, zxffonte, &
	122	zxfqcalving, zxfluxlat, &
	123	zxrunofflic, &
	124	zxtsol, snow_lsc, zxfqfonte, zxqsurf, &
	125	rain_lsc, rain_num, &
	126	!
	127	sens_x, sens_w, &
	128	zxfluxlat_x, zxfluxlat_w, &
	129	!
	130	dtvdf_x, dtvdf_w, &
	131	dqvdf_x, dqvdf_w, &
	132	undi_tke, wake_tke, &
	133	pbl_tke_input, &
	134	t_therm, q_therm, u_therm, v_therm, &
	135	cdragh_x, cdragh_w, &
	136	cdragm_x, cdragm_w, &
	137	kh, kh_x, kh_w, &
	138	!
	139	ale_wake, alp_wake, &
	140	wake_h, wbeff, zmax_th, &
	141	sens, flwp, fiwp, &
	142	ale_bl_stat,alp_bl_conv,alp_bl_det, &
	143	alp_bl_fluct_m,alp_bl_fluct_tke, &
	144	alp_bl_stat, n2, s2, &
	145	proba_notrig, random_notrig, &
	146	!
	147	dnwd, dnwd0, &
	148	upwd, omega, &
	149	epmax_diag, &
	150	ep, &
	151	cldemi, &
	152	cldfra, cldtau, fiwc, &
	153	fl, re, flwc, &
	154	ref_liq, ref_ice, theta, &
	155	ref_liq_pi, ref_ice_pi, &
	156	zphi, wake_omg, zx_rh, &
	157	pmfd, pmfu, &
	158	!
	159	t2m, fluxlat, &
	160	fsollw, evap_pot, &
	161	fsolsw, wfbils, wfbilo, &
	162
	163	pmflxr, pmflxs, prfl, &
	164	psfl, fraca, Vprecip, &
	165	zw2, &
	166
	167	fluxu, fluxv, &
	168	fluxt, &
	169
	170	uwriteSTD, vwriteSTD, & !pour calcul_STDlev.h
	171	wwriteSTD, phiwriteSTD, & !pour calcul_STDlev.h
	172	qwriteSTD, twriteSTD, rhwriteSTD, & !pour calcul_STDlev.h
	173
	174	wdtrainA, wdtrainM, &
	175	beta_prec, &
	176	rneb, &
	177	zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
	178	!
[2469]	179	USE phys_state_var_mod ! Variables sauvegardees de la physique
	180	USE phys_output_var_mod ! Variables pour les ecritures des sorties
	181	USE phys_output_write_mod
	182	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
	183	USE phys_output_mod
	184	USE phys_output_ctrlout_mod
	185	USE iophy
	186	use open_climoz_m, only: open_climoz ! ozone climatology from a file
	187	use regr_pr_av_m, only: regr_pr_av
	188	use netcdf95, only: nf95_close
	189	!IM for NMC files
	190	! use netcdf, only: nf90_fill_real
	191	use netcdf
	192	use mod_phys_lmdz_mpi_data, only: is_mpi_root
	193	USE aero_mod
	194	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
	195	use conf_phys_m, only: conf_phys
	196	use radlwsw_m, only: radlwsw
	197	use phyaqua_mod, only: zenang_an
	198	USE time_phylmdz_mod, only: day_step_phy, annee_ref, day_ref, itau_phy, &
	199	start_time, pdtphys, day_ini
	200	USE tracinca_mod, ONLY: config_inca
[2271]	201	#ifdef CPP_XIOS
[2469]	202	USE wxios, ONLY: missing_val, missing_val_omp
	203	USE xios, ONLY: xios_get_field_attr
[2271]	204	#endif
[1565]	205	#ifdef REPROBUS
[2469]	206	USE CHEM_REP, ONLY : Init_chem_rep_xjour
[1565]	207	#endif
[2469]	208	USE indice_sol_mod
	209	USE phytrac_mod, ONLY : phytrac
[782]	210
[2009]	211	#ifdef CPP_RRTM
[2517]	212	USE YOERAD, ONLY : NRADLP
[2524]	213	USE YOESW, ONLY : RSUN
[2009]	214	#endif
[2469]	215	USE ioipsl_getin_p_mod, ONLY : getin_p
[2003]	216
[2590]	217	USE paramLMDZ_phy_mod
[2294]	218
[2611]	219	USE cmp_seri_mod
	220
[2469]	221	!IM stations CFMIP
	222	USE CFMIP_point_locations
	223	use FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
	224	use ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
[3629]	225	USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
	226
	227
[2469]	228	IMPLICIT none
	229	!>======================================================================
	230	!!
	231	!! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
	232	!!
	233	!! Objet: Moniteur general de la physique du modele
	234	!!AA Modifications quant aux traceurs :
	235	!!AA - uniformisation des parametrisations ds phytrac
	236	!!AA - stockage des moyennes des champs necessaires
	237	!!AA en mode traceur off-line
	238	!!======================================================================
	239	!! CLEFS CPP POUR LES IO
	240	!! =====================
[1352]	241	#define histNMC
[2469]	242	!!======================================================================
	243	!! modif ( P. Le Van , 12/10/98 )
	244	!!
	245	!! Arguments:
	246	!!
	247	!! nlon----input-I-nombre de points horizontaux
	248	!! nlev----input-I-nombre de couches verticales, doit etre egale a klev
	249	!! debut---input-L-variable logique indiquant le premier passage
	250	!! lafin---input-L-variable logique indiquant le dernier passage
	251	!! jD_cur -R-jour courant a l'appel de la physique (jour julien)
	252	!! jH_cur -R-heure courante a l'appel de la physique (jour julien)
	253	!! pdtphys-input-R-pas d'integration pour la physique (seconde)
	254	!! paprs---input-R-pression pour chaque inter-couche (en Pa)
	255	!! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
	256	!! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
	257	!! pphis---input-R-geopotentiel du sol
	258	!! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
	259	!! u-------input-R-vitesse dans la direction X (de O a E) en m/s
	260	!! v-------input-R-vitesse Y (de S a N) en m/s
	261	!! t-------input-R-temperature (K)
	262	!! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
	263	!! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
	264	!! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
[2496]	265	!! d_ql_dyn-input-R-tendance dynamique pour "ql" (kg/kg/s)
	266	!! d_qs_dyn-input-R-tendance dynamique pour "qs" (kg/kg/s)
[2469]	267	!! flxmass_w -input-R- flux de masse verticale
	268	!! d_u-----output-R-tendance physique de "u" (m/s/s)
	269	!! d_v-----output-R-tendance physique de "v" (m/s/s)
	270	!! d_t-----output-R-tendance physique de "t" (K/s)
	271	!! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
	272	!! d_ps----output-R-tendance physique de la pression au sol
	273	!!======================================================================
	274	integer jjmp1
	275	! parameter (jjmp1=jjm+1-1/jjm) ! => (jjmp1=nbp_lat-1/(nbp_lat-1))
	276	! integer iip1
	277	! parameter (iip1=iim+1)
[782]	278
[2469]	279	include "regdim.h"
	280	include "dimsoil.h"
	281	include "clesphys.h"
	282	include "thermcell.h"
	283	!======================================================================
	284	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
	285	PARAMETER (ok_cvl=.TRUE.)
	286	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
	287	PARAMETER (ok_gust=.FALSE.)
	288	integer iflag_radia ! active ou non le rayonnement (MPL)
	289	save iflag_radia
	290	!$OMP THREADPRIVATE(iflag_radia)
	291	!======================================================================
	292	LOGICAL check ! Verifier la conservation du modele en eau
	293	PARAMETER (check=.FALSE.)
	294	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
	295	PARAMETER (ok_stratus=.FALSE.)
	296	!======================================================================
	297	REAL amn, amx
	298	INTEGER igout
	299	!======================================================================
	300	! Clef controlant l'activation du cycle diurne:
	301	! en attente du codage des cles par Fred
	302	INTEGER iflag_cycle_diurne
	303	PARAMETER (iflag_cycle_diurne=1)
	304	!======================================================================
	305	! Modele thermique du sol, a activer pour le cycle diurne:
	306	!cc LOGICAL soil_model
	307	!cc PARAMETER (soil_model=.FALSE.)
	308	!======================================================================
	309	! Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
	310	! le calcul du rayonnement est celle apres la precipitation des nuages.
	311	! Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
	312	! la condensation et la precipitation. Cette cle augmente les impacts
	313	! radiatifs des nuages.
	314	!cc LOGICAL new_oliq
	315	!cc PARAMETER (new_oliq=.FALSE.)
	316	!======================================================================
	317	! Clefs controlant deux parametrisations de l'orographie:
	318	!c LOGICAL ok_orodr
	319	!cc PARAMETER (ok_orodr=.FALSE.)
	320	!cc LOGICAL ok_orolf
	321	!cc PARAMETER (ok_orolf=.FALSE.)
	322	!======================================================================
	323	LOGICAL ok_journe ! sortir le fichier journalier
	324	save ok_journe
	325	!$OMP THREADPRIVATE(ok_journe)
	326	!
	327	LOGICAL ok_mensuel ! sortir le fichier mensuel
	328	save ok_mensuel
	329	!$OMP THREADPRIVATE(ok_mensuel)
	330	!
	331	LOGICAL ok_instan ! sortir le fichier instantane
	332	save ok_instan
	333	!$OMP THREADPRIVATE(ok_instan)
	334	!
	335	LOGICAL ok_LES ! sortir le fichier LES
	336	save ok_LES
	337	!$OMP THREADPRIVATE(ok_LES)
	338	!
	339	LOGICAL callstats ! sortir le fichier stats
	340	save callstats
	341	!$OMP THREADPRIVATE(callstats)
	342	!
	343	LOGICAL ok_region ! sortir le fichier regional
	344	PARAMETER (ok_region=.FALSE.)
	345	!======================================================================
	346	real seuil_inversion
	347	save seuil_inversion
	348	!$OMP THREADPRIVATE(seuil_inversion)
	349	integer iflag_ratqs
	350	save iflag_ratqs
	351	!$OMP THREADPRIVATE(iflag_ratqs)
	352	real facteur
[1507]	353
[2469]	354	REAL wmax_th(klon)
	355	REAL tau_overturning_th(klon)
[878]	356
[2469]	357	integer lmax_th(klon)
	358	integer limbas(klon)
	359	real ratqscth(klon,klev)
	360	real ratqsdiff(klon,klev)
	361	real zqsatth(klon,klev)
[878]	362
[2469]	363	!======================================================================
	364	!
	365	INTEGER ivap ! indice de traceurs pour vapeur d'eau
	366	PARAMETER (ivap=1)
	367	INTEGER iliq ! indice de traceurs pour eau liquide
	368	PARAMETER (iliq=2)
	369	!CR: on ajoute la phase glace
	370	INTEGER isol ! indice de traceurs pour eau glace
	371	PARAMETER (isol=3)
	372	!
	373	!
	374	! Variables argument:
	375	!
	376	INTEGER nlon
	377	INTEGER nlev
	378	REAL,INTENT(IN) :: pdtphys_
	379	! NB: pdtphys to be used in physics is in time_phylmdz_mod
	380	LOGICAL debut, lafin
	381	REAL paprs(klon,klev+1)
	382	REAL pplay(klon,klev)
	383	REAL pphi(klon,klev)
	384	REAL pphis(klon)
	385	REAL presnivs(klev)
	386	REAL znivsig(klev)
	387	real pir
[719]	388
[2469]	389	REAL u(klon,klev)
	390	REAL v(klon,klev)
[2333]	391
[2469]	392	REAL, intent(in):: rot(klon, klev)
	393	! relative vorticity, in s-1, needed for frontal waves
[2333]	394
[2469]	395	REAL t(klon,klev),thetal(klon,klev)
	396	! thetal: ligne suivante a decommenter si vous avez les fichiers
	397	! MPL 20130625
	398	! fth_fonctions.F90 et parkind1.F90
	399	! sinon thetal=theta
	400	! REAL fth_thetae,fth_thetav,fth_thetal
	401	REAL qx(klon,klev,nqtot)
	402	REAL flxmass_w(klon,klev)
	403	REAL d_u(klon,klev)
	404	REAL d_v(klon,klev)
	405	REAL d_t(klon,klev)
	406	REAL d_qx(klon,klev,nqtot)
	407	REAL d_ps(klon)
	408	! Variables pour le transport convectif
	409	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
	410	real wght_cvfd(klon,klev)
[2271]	411	#ifndef CPP_XIOS
[2469]	412	REAL, SAVE :: missing_val
[2271]	413	#endif
[2469]	414	! Variables pour le lessivage convectif
	415	! RomP >>>
	416	real phi2(klon,klev,klev)
	417	real d1a(klon,klev),dam(klon,klev)
[2481]	418	real ev(klon,klev)
[2469]	419	real clw(klon,klev),elij(klon,klev,klev)
	420	real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
	421	! RomP <<<
	422	!IM definition dynamique o_trac dans phys_output_open
	423	! type(ctrl_out) :: o_trac(nqtot)
[524]	424
[2469]	425	! variables a une pression donnee
	426	!
	427	include "declare_STDlev.h"
	428	!
	429	!
	430	include "radopt.h"
	431	!
	432	!
	433	INTEGER debug
	434	INTEGER n
	435	!ym INTEGER npoints
	436	!ym PARAMETER(npoints=klon)
	437	!
	438	INTEGER nregISCtot
	439	PARAMETER(nregISCtot=1)
	440	!
	441	! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties
	442	! sur 1 region rectangulaire y compris pour 1 point
	443	! imin_debut : indice minimum de i; nbpti : nombre de points en
	444	! direction i (longitude)
	445	! jmin_debut : indice minimum de j; nbptj : nombre de points en
	446	! direction j (latitude)
	447	INTEGER imin_debut, nbpti
	448	INTEGER jmin_debut, nbptj
	449	!IM: region='3d' <==> sorties en global
	450	CHARACTER*3 region
	451	PARAMETER(region='3d')
	452	logical ok_hf
	453	!
	454	save ok_hf
	455	!$OMP THREADPRIVATE(ok_hf)
[524]	456
[2469]	457	INTEGER,PARAMETER :: longcles=20
	458	REAL,SAVE :: clesphy0(longcles)
	459	!$OMP THREADPRIVATE(clesphy0)
	460	!
	461	! Variables propres a la physique
	462	INTEGER itap
	463	SAVE itap ! compteur pour la physique
	464	!$OMP THREADPRIVATE(itap)
[2235]	465
[2469]	466	INTEGER, SAVE :: abortphy=0 ! Reprere si on doit arreter en fin de phys
	467	!$OMP THREADPRIVATE(abortphy)
	468	!
	469	REAL,save :: solarlong0
	470	!$OMP THREADPRIVATE(solarlong0)
[987]	471
[2469]	472	!
	473	! Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
	474	!
	475	!IM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
	476	REAL zulow(klon),zvlow(klon)
	477	!
	478	INTEGER igwd,idx(klon),itest(klon)
	479	!
	480	! REAL,allocatable,save :: run_off_lic_0(:)
	481	! !$OMP THREADPRIVATE(run_off_lic_0)
	482	!ym SAVE run_off_lic_0
	483	!KE43
	484	! Variables liees a la convection de K. Emanuel (sb):
	485	!
	486	REAL bas, top ! cloud base and top levels
	487	SAVE bas
	488	SAVE top
	489	!$OMP THREADPRIVATE(bas, top)
	490	!------------------------------------------------------------------
	491	! Upmost level reached by deep convection and related variable (jyg)
	492	!
	493	INTEGER izero
	494	INTEGER k_upper_cv
	495	!------------------------------------------------------------------
	496	!
	497	!==========================================================================
	498	!CR04.12.07: on ajoute les nouvelles variables du nouveau schema
	499	!de convection avec poches froides
	500	! Variables li\'ees \`a la poche froide (jyg)
[879]	501
[2469]	502	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
	503	!
	504	REAL wape_prescr, fip_prescr
	505	INTEGER it_wape_prescr
	506	SAVE wape_prescr, fip_prescr, it_wape_prescr
	507	!$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
	508	!
	509	! variables supplementaires de concvl
	510	REAL Tconv(klon,klev)
	511	REAL sij(klon,klev,klev)
[970]	512
[2469]	513	real, save :: alp_bl_prescr=0.
	514	real, save :: ale_bl_prescr=0.
[979]	515
[2469]	516	real, save :: wake_s_min_lsp=0.1
[1516]	517
[2469]	518	!$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
	519	!$OMP THREADPRIVATE(wake_s_min_lsp)
[970]	520
[1516]	521
[2469]	522	real ok_wk_lsp(klon)
[1516]	523
[2469]	524	!RC
	525	! Variables li\'ees \`a la poche froide (jyg et rr)
	526	! Version diagnostique pour l'instant : pas de r\'etroaction sur
	527	! la convection
[879]	528
[2469]	529	REAL t_wake(klon,klev),q_wake(klon,klev) ! wake pour la convection
[879]	530
[2469]	531	REAL wake_dth(klon,klev) ! wake : temp pot difference
[879]	532
[2469]	533	REAL wake_d_deltat_gw(klon,klev)! wake : delta T tendency due to
	534	! Gravity Wave (/s)
	535	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta
	536	! transported by LS omega
	537	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of
	538	! large scale omega
	539	REAL wake_dtKE(klon,klev) ! Wake : differential heating
	540	! (wake - unpertubed) CONV
	541	REAL wake_dqKE(klon,klev) ! Wake : differential moistening
	542	! (wake - unpertubed) CONV
	543	REAL wake_dtPBL(klon,klev) ! Wake : differential heating
	544	! (wake - unpertubed) PBL
	545	REAL wake_dqPBL(klon,klev) ! Wake : differential moistening
	546	! (wake - unpertubed) PBL
	547	REAL wake_ddeltat(klon,klev),wake_ddeltaq(klon,klev)
	548	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
	549	REAL wake_spread(klon,klev) ! spreading term in wake_delt
	550	!
	551	!pourquoi y'a pas de save??
	552	!
	553	INTEGER wake_k(klon) ! Wake sommet
	554	!
	555	REAL t_undi(klon,klev) ! temperature moyenne dans la zone
	556	! non perturbee
	557	REAL q_undi(klon,klev) ! humidite moyenne dans la zone
	558	! non perturbee
	559	!
	560	!jyg<
	561	!cc REAL wake_pe(klon) ! Wake potential energy - WAPE
	562	!>jyg
[879]	563
[2469]	564	REAL wake_gfl(klon) ! Gust Front Length
	565	REAL wake_dens(klon)
	566	!
	567	!
	568	REAL dt_dwn(klon,klev)
	569	REAL dq_dwn(klon,klev)
	570	REAL wdt_PBL(klon,klev)
	571	REAL udt_PBL(klon,klev)
	572	REAL wdq_PBL(klon,klev)
	573	REAL udq_PBL(klon,klev)
	574	REAL M_dwn(klon,klev)
	575	REAL M_up(klon,klev)
	576	REAL dt_a(klon,klev)
	577	REAL dq_a(klon,klev)
	578	REAL d_t_adjwk(klon,klev) !jyg
	579	REAL d_q_adjwk(klon,klev) !jyg
	580	LOGICAL,SAVE :: ok_adjwk=.FALSE.
	581	!$OMP THREADPRIVATE(ok_adjwk)
[2613]	582	REAL,SAVE :: oliqmax=999.
	583	!$OMP THREADPRIVATE(oliqmax)
[2469]	584	REAL, SAVE :: alp_offset
	585	!$OMP THREADPRIVATE(alp_offset)
[1403]	586
[2469]	587	! !!
	588	!=================================================================
	589	! PROVISOIRE : DECOUPLAGE PBL/WAKE
	590	! --------------------------------
	591	REAL wake_deltat_sav(klon,klev)
	592	REAL wake_deltaq_sav(klon,klev)
	593	!=================================================================
[2159]	594
[2469]	595	!
	596	!RR:fin declarations poches froides
	597	!==========================================================================
[1032]	598
[2469]	599	REAL ztv(klon,klev),ztva(klon,klev)
	600	REAL zpspsk(klon,klev)
	601	REAL ztla(klon,klev),zqla(klon,klev)
	602	REAL zthl(klon,klev)
[1638]	603
[2469]	604	!cc nrlmd le 10/04/2012
[1638]	605
[2469]	606	!--------Stochastic Boundary Layer Triggering: ALE_BL--------
	607	!---Propri\'et\'es du thermiques au LCL
	608	real zlcl_th(klon) ! Altitude du LCL calcul\'e
	609	! continument (pcon dans
	610	! thermcell_main.F90)
	611	real fraca0(klon) ! Fraction des thermiques au LCL
	612	real w0(klon) ! Vitesse des thermiques au LCL
	613	real w_conv(klon) ! Vitesse verticale de grande \'echelle au LCL
	614	real tke0(klon,klev+1) ! TKE au d\'ebut du pas de temps
	615	real therm_tke_max0(klon) ! TKE dans les thermiques au LCL
	616	real env_tke_max0(klon) ! TKE dans l'environnement au LCL
[1638]	617
[2469]	618	!---D\'eclenchement stochastique
	619	integer :: tau_trig(klon)
[1638]	620
[2469]	621	REAL,SAVE :: random_notrig_max=1.
	622	!$OMP THREADPRIVATE(random_notrig_max)
[2294]	623
[2469]	624	!--------Statistical Boundary Layer Closure: ALP_BL--------
	625	!---Profils de TKE dans et hors du thermique
	626	real therm_tke_max(klon,klev) ! Profil de TKE dans les thermiques
	627	real env_tke_max(klon,klev) ! Profil de TKE dans l'environnement
[1638]	628
	629
[2469]	630	!cc fin nrlmd le 10/04/2012
[782]	631
[2469]	632	! Variables locales pour la couche limite (al1):
	633	!
	634	!Al1 REAL pblh(klon) ! Hauteur de couche limite
	635	!Al1 SAVE pblh
	636	!34EK
	637	!
	638	! Variables locales:
	639	!
	640	!AA
	641	!AA Pour phytrac
	642	REAL u1(klon) ! vents dans la premiere couche U
	643	REAL v1(klon) ! vents dans la premiere couche V
[524]	644
[2469]	645	!@$$ LOGICAL offline ! Controle du stockage ds "physique"
	646	!@$$ PARAMETER (offline=.false.)
	647	!@$$ INTEGER physid
	648	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
	649	REAL frac_nucl(klon,klev) ! idem (nucleation)
	650	! RomP >>>
	651	REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
	652	! RomP <<<
[2068]	653
[2469]	654	REAL :: calday
[782]	655
[2469]	656	!IM cf FH pour Tiedtke 080604
	657	REAL rain_tiedtke(klon),snow_tiedtke(klon)
	658	!
	659	!IM 050204 END
	660	REAL devap(klon) ! evaporation et sa derivee
	661	REAL dsens(klon) ! chaleur sensible et sa derivee
[1279]	662
[2469]	663	!
	664	! Conditions aux limites
	665	!
	666	!
	667	REAL :: day_since_equinox
	668	! Date de l'equinoxe de printemps
	669	INTEGER, parameter :: mth_eq=3, day_eq=21
	670	REAL :: jD_eq
[1279]	671
[2469]	672	LOGICAL, parameter :: new_orbit = .true.
[524]	673
[2469]	674	!
	675	INTEGER lmt_pas
	676	SAVE lmt_pas ! frequence de mise a jour
	677	!$OMP THREADPRIVATE(lmt_pas)
	678	real zmasse(klon, nbp_lev),exner(klon, nbp_lev)
	679	! (column-density of mass of air in a cell, in kg m-2)
	680	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
[1797]	681
[2469]	682	!IM sorties
	683	REAL un_jour
	684	PARAMETER(un_jour=86400.)
	685	INTEGER itapm1 !pas de temps de la physique du(es) mois precedents
	686	SAVE itapm1 !mis a jour le dernier pas de temps du mois en cours
	687	!$OMP THREADPRIVATE(itapm1)
	688	!======================================================================
	689	!
	690	! Declaration des procedures appelees
	691	!
	692	EXTERNAL angle ! calculer angle zenithal du soleil
	693	EXTERNAL alboc ! calculer l'albedo sur ocean
	694	EXTERNAL ajsec ! ajustement sec
	695	EXTERNAL conlmd ! convection (schema LMD)
	696	!KE43
	697	EXTERNAL conema3 ! convect4.3
	698	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
	699	!AA
	700	! JBM (3/14) fisrtilp_tr not loaded
	701	! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
	702	! ! stockage des coefficients necessaires au
	703	! ! lessivage OFF-LINE et ON-LINE
	704	EXTERNAL hgardfou ! verifier les temperatures
	705	EXTERNAL nuage ! calculer les proprietes radiatives
	706	!C EXTERNAL o3cm ! initialiser l'ozone
	707	EXTERNAL orbite ! calculer l'orbite terrestre
	708	EXTERNAL phyetat0 ! lire l'etat initial de la physique
	709	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
	710	EXTERNAL suphel ! initialiser certaines constantes
	711	EXTERNAL transp ! transport total de l'eau et de l'energie
	712	!IM
	713	EXTERNAL haut2bas !variables de haut en bas
	714	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
	715	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
	716	! EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
	717	! EXTERNAL moyglo_aire
	718	! moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
	719	! par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
	720	!
	721	!
	722	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	723	! Local variables
	724	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	725	!
	726	REAL rhcl(klon,klev) ! humiditi relative ciel clair
	727	REAL dialiq(klon,klev) ! eau liquide nuageuse
	728	REAL diafra(klon,klev) ! fraction nuageuse
	729	REAL cldliq(klon,klev) ! eau liquide nuageuse
	730	!
	731	!XXX PB
	732	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
	733	!
	734	REAL zxfluxt(klon, klev)
	735	REAL zxfluxq(klon, klev)
	736	REAL zxfluxu(klon, klev)
	737	REAL zxfluxv(klon, klev)
[1797]	738
[2469]	739	! Le rayonnement n'est pas calcule tous les pas, il faut donc
	740	! sauvegarder les sorties du rayonnement
	741	!ym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
	742	!ym SAVE sollwdownclr, toplwdown, toplwdownclr
	743	!ym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
	744	!
	745	INTEGER itaprad
	746	SAVE itaprad
	747	!$OMP THREADPRIVATE(itaprad)
	748	!
	749	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
	750	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
[1797]	751
[2469]	752	!
	753	! REAL zxsnow(klon)
	754	REAL zxsnow_dummy(klon)
	755	REAL zsav_tsol(klon)
	756	!
	757	REAL dist, rmu0(klon), fract(klon)
	758	REAL zrmu0(klon), zfract(klon)
	759	REAL zdtime, zdtime1, zdtime2, zlongi
	760	!
	761	REAL qcheck
	762	REAL z_avant(klon), z_apres(klon), z_factor(klon)
	763	LOGICAL zx_ajustq
	764	!
	765	REAL za, zb
	766	REAL zx_t, zx_qs, zdelta, zcor, zlvdcp, zlsdcp
	767	real zqsat(klon,klev)
	768	!
	769	INTEGER i, k, iq, ig, j, nsrf, ll, l, iiq
	770	!
	771	REAL t_coup
	772	PARAMETER (t_coup=234.0)
[1797]	773
[2469]	774	!ym A voir plus tard !!
	775	!ym REAL zx_relief(iim,jjmp1)
	776	!ym REAL zx_aire(iim,jjmp1)
	777	!
	778	! Grandeurs de sorties
	779	REAL s_capCL(klon)
	780	REAL s_oliqCL(klon), s_cteiCL(klon)
	781	REAL s_trmb1(klon), s_trmb2(klon)
	782	REAL s_trmb3(klon)
	783	!KE43
	784	! Variables locales pour la convection de K. Emanuel (sb):
[524]	785
[2469]	786	REAL tvp(klon,klev) ! virtual temp of lifted parcel
	787	CHARACTER*40 capemaxcels !max(CAPE)
[1412]	788
[2469]	789	REAL rflag(klon) ! flag fonctionnement de convect
	790	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
[904]	791
[2469]	792	! -- convect43:
	793	INTEGER ntra ! nb traceurs pour convect4.3
	794	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
	795	REAL dplcldt(klon), dplcldr(klon)
	796	!? . condm_con(klon,klev),conda_con(klon,klev),
	797	!? . mr_con(klon,klev),ep_con(klon,klev)
	798	!? . ,sadiab(klon,klev),wadiab(klon,klev)
	799	! --
	800	!34EK
	801	!
	802	! Variables du changement
	803	!
	804	! con: convection
	805	! lsc: condensation a grande echelle (Large-Scale-Condensation)
	806	! ajs: ajustement sec
	807	! eva: evaporation de l'eau liquide nuageuse
	808	! vdf: couche limite (Vertical DiFfusion)
[2611]	809	!
[2469]	810	! tendance nulles
	811	REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
[2611]	812	!
	813	! Flag pour pouvoir ne pas ajouter les tendances.
	814	! Par defaut, les tendances doivente etre ajoutees et
	815	! flag_inhib_tend = 0
	816	! flag_inhib_tend > 0 : tendances non ajoutees, avec un nombre
	817	! croissant de print quand la valeur du flag augmente
	818	!!! attention, ce flag doit etre change avec prudence !!!
	819	INTEGER :: flag_inhib_tend = 0 ! 0 is the default value
	820	!! INTEGER :: flag_inhib_tend = 2
[524]	821
[2469]	822	!
	823	!********************************************************
	824	! declarations
[524]	825
[2469]	826	!********************************************************
	827	!IM 081204 END
	828	!
	829	REAL pen_u(klon,klev), pen_d(klon,klev)
	830	REAL pde_u(klon,klev), pde_d(klon,klev)
	831	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
	832	!
	833	REAL ratqsc(klon,klev)
	834	real ratqsbas,ratqshaut,tau_ratqs
	835	save ratqsbas,ratqshaut,tau_ratqs
	836	!$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
[2534]	837	REAL, SAVE :: ratqsp0=50000., ratqsdp=20000.
	838	!$OMP THREADPRIVATE(ratqsp0, ratqsdp)
[644]	839
[2469]	840	! Parametres lies au nouveau schema de nuages (SB, PDF)
	841	real fact_cldcon
	842	real facttemps
	843	logical ok_newmicro
	844	save ok_newmicro
	845	!$OMP THREADPRIVATE(ok_newmicro)
	846	!real ref_liq_pi(klon,klev), ref_ice_pi(klon,klev)
	847	save fact_cldcon,facttemps
	848	!$OMP THREADPRIVATE(fact_cldcon,facttemps)
[524]	849
[2469]	850	integer iflag_cld_th
	851	save iflag_cld_th
	852	!$OMP THREADPRIVATE(iflag_cld_th)
	853	logical ptconv(klon,klev)
	854	!IM cf. AM 081204 BEG
	855	logical ptconvth(klon,klev)
	856	!IM cf. AM 081204 END
	857	!
	858	! Variables liees a l'ecriture de la bande histoire physique
	859	!
	860	!======================================================================
	861	!
[2068]	862
[2469]	863	!
	864	integer itau_w ! pas de temps ecriture = itap + itau_phy
	865	!
	866	!
	867	! Variables locales pour effectuer les appels en serie
	868	!
	869	!IM RH a 2m (la surface)
	870	REAL Lheat
[524]	871
[2469]	872	INTEGER length
	873	PARAMETER ( length = 100 )
	874	REAL tabcntr0( length )
	875	!
	876	INTEGER ndex2d(nbp_lon*nbp_lat)
	877	!IM
	878	!
	879	!IM AMIP2 BEG
	880	REAL moyglo, mountor
	881	!IM 141004 BEG
	882	REAL zustrdr(klon), zvstrdr(klon)
	883	REAL zustrli(klon), zvstrli(klon)
	884	REAL zustrph(klon), zvstrph(klon)
	885	REAL aam, torsfc
	886	!IM 141004 END
	887	!IM 190504 BEG
	888	INTEGER ij
	889	! INTEGER imp1jmp1
	890	! PARAMETER(imp1jmp1=(iim+1)*jjmp1)
	891	!ym A voir plus tard
	892	! REAL zx_tmp((nbp_lon+1)*nbp_lat)
	893	! REAL airedyn(nbp_lon+1,nbp_lat)
	894	!IM 190504 END
	895	LOGICAL ok_msk
	896	REAL msk(klon)
	897	!IM
	898	REAL airetot, pi
	899	!ym A voir plus tard
	900	!ym REAL zm_wo(jjmp1, klev)
	901	!IM AMIP2 END
	902	!
	903	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
	904	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
	905	REAL zx_tmp_2d(nbp_lon,nbp_lat)
	906	REAL zx_lon(nbp_lon,nbp_lat)
	907	REAL zx_lat(nbp_lon,nbp_lat)
	908	!
	909	INTEGER nid_day_seri, nid_ctesGCM
	910	SAVE nid_day_seri, nid_ctesGCM
	911	!$OMP THREADPRIVATE(nid_day_seri,nid_ctesGCM)
	912	!
	913	!IM 280405 BEG
	914	! INTEGER nid_bilKPins, nid_bilKPave
	915	! SAVE nid_bilKPins, nid_bilKPave
	916	! !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
	917	!
	918	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
	919	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
	920	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
	921	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
	922	!
[2590]	923	REAL zjulian
	924	SAVE zjulian
	925	!$OMP THREADPRIVATE(zjulian)
	926
[2469]	927	INTEGER nhori, nvert
	928	REAL zsto
	929	REAL zstophy, zout
[2068]	930
[2469]	931	character*20 modname
	932	character*80 abort_message
	933	logical, save :: ok_sync, ok_sync_omp
	934	!$OMP THREADPRIVATE(ok_sync)
	935	real date0
	936	integer idayref
[524]	937
[2469]	938	! essai writephys
	939	integer fid_day, fid_mth, fid_ins
	940	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
	941	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
	942	parameter (prof2d_on = 1, prof3d_on = 2, &
	943	prof2d_av = 3, prof3d_av = 4)
	944	! Variables liees au bilan d'energie et d'enthalpi
	945	REAL ztsol(klon)
	946	REAL d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec
	947	REAL d_h_vcol_phy
	948	REAL fs_bound, fq_bound
	949	SAVE d_h_vcol_phy
	950	!$OMP THREADPRIVATE(d_h_vcol_phy)
	951	REAL zero_v(klon)
	952	CHARACTER*40 ztit
	953	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
	954	SAVE ip_ebil
	955	DATA ip_ebil/0/
	956	!$OMP THREADPRIVATE(ip_ebil)
	957	INTEGER if_ebil ! level for energy conserv. dignostics
	958	SAVE if_ebil
	959	!$OMP THREADPRIVATE(if_ebil)
	960	REAL q2m(klon,nbsrf) ! humidite a 2m
[524]	961
[2469]	962	!IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
	963	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
	964	CHARACTER*40 tinst, tave, typeval
	965	REAL cldtaupi(klon,klev) ! Cloud optical thickness for
	966	! pre-industrial (pi) aerosols
[524]	967
[959]	968
[2469]	969	! Aerosol optical properties
	970	CHARACTER*4, DIMENSION(naero_grp) :: rfname
	971	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass
	972	! concentration
	973	! for all soluble
	974	! aerosols[ug/m3]
	975	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi
	976	! - " - (pre-industrial value)
[1279]	977
[2469]	978	! Parameters
	979	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
	980	LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
	981	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
	982	SAVE ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1
	983	!$OMP THREADPRIVATE(ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1)
	984	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
	985	! false : lecture des aerosol dans un fichier
	986	!$OMP THREADPRIVATE(aerosol_couple)
[3340]	987	LOGICAL, SAVE :: chemistry_couple ! true : use INCA chemistry O3 / false : use offline chemistry O3
	988	!$OMP THREADPRIVATE(chemistry_couple)
[2469]	989	INTEGER, SAVE :: flag_aerosol
	990	!$OMP THREADPRIVATE(flag_aerosol)
	991	LOGICAL, SAVE :: new_aod
	992	!$OMP THREADPRIVATE(new_aod)
	993	!
	994	!--STRAT AEROSOL
[2530]	995	INTEGER, SAVE :: flag_aerosol_strat
[2469]	996	!$OMP THREADPRIVATE(flag_aerosol_strat)
	997	!c-fin STRAT AEROSOL
	998	!
	999	! Declaration des constantes et des fonctions thermodynamiques
	1000	!
	1001	LOGICAL,SAVE :: first=.true.
	1002	!$OMP THREADPRIVATE(first)
[1279]	1003
[2469]	1004	integer, save:: read_climoz ! read ozone climatology
	1005	! (let it keep the default OpenMP shared attribute)
	1006	! Allowed values are 0, 1 and 2
	1007	! 0: do not read an ozone climatology
	1008	! 1: read a single ozone climatology that will be used day and night
	1009	! 2: read two ozone climatologies, the average day and night
	1010	! climatology and the daylight climatology
[1279]	1011
[2469]	1012	integer, save:: ncid_climoz ! NetCDF file containing ozone climatologies
	1013	! (let it keep the default OpenMP shared attribute)
[1279]	1014
[2469]	1015	real, pointer, save:: press_climoz(:)
	1016	! (let it keep the default OpenMP shared attribute)
	1017	! edges of pressure intervals for ozone climatologies, in Pa, in strictly
	1018	! ascending order
[1279]	1019
[2469]	1020	integer, save:: co3i = 0
	1021	! time index in NetCDF file of current ozone fields
	1022	!$OMP THREADPRIVATE(co3i)
[1863]	1023
[2469]	1024	integer ro3i
	1025	! required time index in NetCDF file for the ozone fields, between 1
	1026	! and 360
[1863]	1027
[2469]	1028	INTEGER ierr
	1029	include "YOMCST.h"
	1030	include "YOETHF.h"
	1031	include "FCTTRE.h"
	1032	!IM 100106 BEG : pouvoir sortir les ctes de la physique
	1033	include "conema3.h"
	1034	include "fisrtilp.h"
	1035	include "nuage.h"
	1036	include "compbl.h"
	1037	!IM 100106 END : pouvoir sortir les ctes de la physique
	1038	!
	1039	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1040	! Declarations pour Simulateur COSP
	1041	!============================================================
	1042	real :: mr_ozone(klon,klev)
[1538]	1043
[2469]	1044	!IM stations CFMIP
	1045	INTEGER, SAVE :: nCFMIP
	1046	!$OMP THREADPRIVATE(nCFMIP)
	1047	INTEGER, PARAMETER :: npCFMIP=120
	1048	INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
	1049	REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
	1050	!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
	1051	INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
	1052	REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
	1053	!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
	1054	INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
	1055	!$OMP THREADPRIVATE(iGCM, jGCM)
	1056	logical, dimension(nfiles) :: phys_out_filestations
	1057	logical, parameter :: lNMC=.FALSE.
[1539]	1058
[2469]	1059	!IM betaCRF
	1060	REAL, SAVE :: pfree, beta_pbl, beta_free
	1061	!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
	1062	REAL, SAVE :: lon1_beta, lon2_beta, lat1_beta, lat2_beta
	1063	!$OMP THREADPRIVATE(lon1_beta, lon2_beta, lat1_beta, lat2_beta)
	1064	LOGICAL, SAVE :: mskocean_beta
	1065	!$OMP THREADPRIVATE(mskocean_beta)
	1066	REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et
	1067	! cldemirad pour evaluer les
	1068	! retros liees aux CRF
	1069	REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique
	1070	! pour radlwsw pour
	1071	! tester "CRF off"
	1072	REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique
	1073	! pour radlwsw pour
	1074	! tester "CRF off"
	1075	REAL, dimension(klon, klev) :: cldemirad ! emissivite pour
	1076	! radlwsw pour tester
	1077	! "CRF off"
	1078	REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse
[1735]	1079
[3433]	1080	#ifdef INCA
	1081	! set de variables utilisees pour l'initialisation des valeurs provenant de INCA
	1082	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_tauinca
	1083	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_pizinca
	1084	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_cginca
	1085	REAL, DIMENSION(klon,klev,nbands) :: init_ccminca
	1086	#endif
	1087	REAL, DIMENSION(klon,nbtr) :: init_source
	1088
[2469]	1089	INTEGER :: nbtr_tmp ! Number of tracer inside concvl
	1090	REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
[2907]	1091	REAL, dimension(klon,klev) :: ch_in ! Condensed humidity entering in phytrac (eau liquide)
[2469]	1092	integer iostat
[1539]	1093
[2469]	1094	REAL zzz
	1095	!albedo SB >>>
	1096	real,dimension(6),save :: SFRWL
	1097	!albedo SB <<<
[1955]	1098
[2485]	1099	!--OB variables for mass fixer (hard coded for now)
[2477]	1100	logical, parameter :: mass_fixer=.false.
[2476]	1101	real qql1(klon),qql2(klon),zdz,corrqql
	1102
[2469]	1103	! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
	1104	jjmp1=nbp_lat
[2344]	1105
[2469]	1106	!======================================================================
	1107	! Gestion calendrier : mise a jour du module phys_cal_mod
	1108	!
	1109	pdtphys=pdtphys_
	1110	CALL update_time(pdtphys)
[1355]	1111
[2469]	1112	!======================================================================
	1113	! Ecriture eventuelle d'un profil verticale en entree de la physique.
	1114	! Utilise notamment en 1D mais peut etre active egalement en 3D
	1115	! en imposant la valeur de igout.
	1116	!======================================================================d
	1117	if (prt_level.ge.1) then
	1118	igout=klon/2+1/klon
	1119	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	1120	write(lunout,*) 'igout, lat, lon ',igout, latitude_deg(igout), &
	1121	longitude_deg(igout)
	1122	write(lunout,*) &
	1123	'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
	1124	write(lunout,*) &
	1125	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
[879]	1126
[2469]	1127	write(lunout,*) 'paprs, play, phi, u, v, t'
	1128	do k=1,klev
	1129	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
	1130	u(igout,k),v(igout,k),t(igout,k)
	1131	enddo
	1132	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
	1133	do k=1,klev
	1134	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
	1135	enddo
	1136	endif
[879]	1137
[2469]	1138	!======================================================================
[879]	1139
[2469]	1140	if (first) then
[1403]	1141
[2469]	1142	!CR:nvelles variables convection/poches froides
	1143
	1144	print*, '================================================='
	1145	print*, 'Allocation des variables locales et sauvegardees'
	1146	call phys_local_var_init
	1147	!
	1148	pasphys=pdtphys
	1149	! appel a la lecture du run.def physique
	1150	call conf_phys(ok_journe, ok_mensuel, &
	1151	ok_instan, ok_hf, &
	1152	ok_LES, &
	1153	callstats, &
	1154	solarlong0,seuil_inversion, &
	1155	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
	1156	iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
[3340]	1157	ok_ade, ok_aie, ok_cdnc, aerosol_couple, chemistry_couple, &
[2469]	1158	flag_aerosol, flag_aerosol_strat, new_aod, &
	1159	bl95_b0, bl95_b1, &
	1160	! nv flags pour la convection et les
	1161	! poches froides
	1162	read_climoz, &
	1163	alp_offset)
	1164	call phys_state_var_init(read_climoz)
	1165	call phys_output_var_init
	1166	print*, '================================================='
	1167	!
	1168	!CR: check sur le nb de traceurs de l eau
	1169	if ((iflag_ice_thermo.gt.0).and.(nqo==2)) then
	1170	WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
	1171	'(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
[2224]	1172	STOP
[2469]	1173	endif
[2224]	1174
[2469]	1175	dnwd0=0.0
	1176	ftd=0.0
	1177	fqd=0.0
	1178	cin=0.
	1179	!ym Attention pbase pas initialise dans concvl !!!!
	1180	pbase=0
	1181	!IM 180608
[904]	1182
[2469]	1183	itau_con=0
	1184	first=.false.
[1797]	1185
[2469]	1186	endif ! first
[1797]	1187
[2469]	1188	!ym => necessaire pour iflag_con != 2
	1189	pmfd(:,:) = 0.
	1190	pen_u(:,:) = 0.
	1191	pen_d(:,:) = 0.
	1192	pde_d(:,:) = 0.
	1193	pde_u(:,:) = 0.
	1194	aam=0.
	1195	d_t_adjwk(:,:)=0
	1196	d_q_adjwk(:,:)=0
[1797]	1197
[2469]	1198	alp_bl_conv(:)=0.
[2245]	1199
[2469]	1200	torsfc=0.
	1201	forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
[1797]	1202
[2469]	1203	modname = 'physiq'
	1204	!IM
	1205	IF (ip_ebil_phy.ge.1) THEN
	1206	DO i=1,klon
	1207	zero_v(i)=0.
	1208	END DO
	1209	END IF
[644]	1210
[2469]	1211	IF (debut) THEN
	1212	CALL suphel ! initialiser constantes et parametres phys.
	1213	CALL getin_p('random_notrig_max',random_notrig_max)
	1214	CALL getin_p('ok_adjwk',ok_adjwk)
[2613]	1215	CALL getin_p('oliqmax',oliqmax)
[2534]	1216	CALL getin_p('ratqsp0',ratqsp0)
	1217	CALL getin_p('ratqsdp',ratqsdp)
[2469]	1218	ENDIF
[878]	1219
[2469]	1220	if(prt_level.ge.1) print*,'CONVERGENCE PHYSIQUE THERM 1 '
[1279]	1221
[959]	1222
[2469]	1223	!======================================================================
	1224	! Gestion calendrier : mise a jour du module phys_cal_mod
	1225	!
	1226	! CALL phys_cal_update(jD_cur,jH_cur)
[1279]	1227
[2469]	1228	!
	1229	! Si c'est le debut, il faut initialiser plusieurs choses
	1230	! ********
	1231	!
	1232	IF (debut) THEN
	1233	!rv CRinitialisation de wght_th et lalim_conv pour la
	1234	!definition de la couche alimentation de la convection a partir
	1235	!des caracteristiques du thermique
	1236	wght_th(:,:)=1.
	1237	lalim_conv(:)=1
	1238	!RC
	1239	ustar(:,:)=0.
[2569]	1240	! u10m(:,:)=0.
	1241	! v10m(:,:)=0.
[2469]	1242	rain_con(:)=0.
	1243	snow_con(:)=0.
	1244	topswai(:)=0.
	1245	topswad(:)=0.
	1246	solswai(:)=0.
	1247	solswad(:)=0.
[959]	1248
[2469]	1249	wmax_th(:)=0.
	1250	tau_overturning_th(:)=0.
[645]	1251
[2469]	1252	IF (type_trac == 'inca') THEN
	1253	! jg : initialisation jusqu'au ces variables sont dans restart
	1254	ccm(:,:,:) = 0.
	1255	tau_aero(:,:,:,:) = 0.
	1256	piz_aero(:,:,:,:) = 0.
	1257	cg_aero(:,:,:,:) = 0.
[2372]	1258
[2469]	1259	config_inca='none' ! default
	1260	CALL getin_p('config_inca',config_inca)
[2372]	1261
[2469]	1262	ELSE
	1263	config_inca='none' ! default
	1264	END IF
[782]	1265
[2469]	1266	IF (aerosol_couple .AND. (config_inca /= "aero" &
	1267	.AND. config_inca /= "aeNP ")) THEN
	1268	abort_message &
	1269	= 'if aerosol_couple is activated, config_inca need to be ' &
	1270	// 'aero or aeNP'
	1271	CALL abort_physic (modname,abort_message,1)
	1272	ENDIF
[2443]	1273
	1274
[1863]	1275
[2469]	1276	rnebcon0(:,:) = 0.0
	1277	clwcon0(:,:) = 0.0
	1278	rnebcon(:,:) = 0.0
	1279	clwcon(:,:) = 0.0
[1863]	1280
[2469]	1281	!IM
	1282	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
	1283	!
	1284	print*,'iflag_coupl,iflag_clos,iflag_wake', &
	1285	iflag_coupl,iflag_clos,iflag_wake
	1286	print*,'iflag_CYCLE_DIURNE', iflag_cycle_diurne
	1287	!
	1288	IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
	1289	abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
	1290	CALL abort_physic (modname,abort_message,1)
	1291	ENDIF
	1292	!
	1293	!
	1294	! Initialiser les compteurs:
	1295	!
	1296	itap = 0
	1297	itaprad = 0
[878]	1298
[2469]	1299	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1300	!! Un petit travail \`a faire ici.
	1301	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[878]	1302
[2469]	1303	if (iflag_pbl>1) then
	1304	PRINT*, "Using method MELLOR&YAMADA"
	1305	endif
[956]	1306
[2469]	1307	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1308	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
	1309	! phylmd plutot que dyn3d
	1310	! Attention : la version precedente n'etait pas tres propre.
	1311	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
	1312	! pour obtenir le meme resultat.
	1313	dtime=pdtphys
	1314	IF (MOD(INT(86400./dtime),nbapp_rad).EQ.0) THEN
	1315	radpas = NINT( 86400./dtime/nbapp_rad)
	1316	ELSE
	1317	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
	1318	'multiple de nbapp_rad'
	1319	WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
	1320	'mais 1+1<>2'
	1321	abort_message='nbre de pas de temps physique n est pas multiple ' &
	1322	// 'de nbapp_rad'
	1323	call abort_physic(modname,abort_message,1)
	1324	ENDIF
	1325	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[524]	1326
[2469]	1327	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
[2565]	1328	!jyg<
[2469]	1329	IF (klon_glo==1) THEN
[2565]	1330	pbl_tke(:,:,is_ave) = 0.
	1331	DO nsrf=1,nbsrf
	1332	DO k = 1,klev+1
	1333	pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
	1334	+pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
	1335	ENDDO
	1336	ENDDO
	1337	!>jyg
[2469]	1338	ENDIF
	1339	!IM begin
	1340	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
	1341	,ratqs(1,1)
	1342	!IM end
[878]	1343
	1344
[524]	1345
[2469]	1346	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1347	!
	1348	! on remet le calendrier a zero
	1349	!
	1350	IF (raz_date .eq. 1) THEN
	1351	itau_phy = 0
	1352	ENDIF
[524]	1353
[2469]	1354	CALL printflag( tabcntr0,radpas,ok_journe, &
	1355	ok_instan, ok_region )
	1356	!
	1357	IF (ABS(dtime-pdtphys).GT.0.001) THEN
	1358	WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
	1359	pdtphys
	1360	abort_message='Pas physique n est pas correct '
	1361	! call abort_physic(modname,abort_message,1)
	1362	dtime=pdtphys
	1363	ENDIF
	1364	IF (nlon .NE. klon) THEN
	1365	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
	1366	klon
	1367	abort_message='nlon et klon ne sont pas coherents'
	1368	call abort_physic(modname,abort_message,1)
	1369	ENDIF
	1370	IF (nlev .NE. klev) THEN
	1371	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
	1372	klev
	1373	abort_message='nlev et klev ne sont pas coherents'
	1374	call abort_physic(modname,abort_message,1)
	1375	ENDIF
	1376	!
	1377	IF (dtime*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
	1378	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
	1379	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
	1380	abort_message='Nbre d appels au rayonnement insuffisant'
	1381	call abort_physic(modname,abort_message,1)
	1382	ENDIF
	1383	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
	1384	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
	1385	ok_cvl
	1386	!
	1387	!KE43
	1388	! Initialisation pour la convection de K.E. (sb):
	1389	IF (iflag_con.GE.3) THEN
[524]	1390
[2469]	1391	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
	1392	WRITE(lunout,*) &
	1393	"On va utiliser le melange convectif des traceurs qui"
	1394	WRITE(lunout,*)"est calcule dans convect4.3"
	1395	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
[524]	1396
[2469]	1397	DO i = 1, klon
	1398	ema_cbmf(i) = 0.
	1399	ema_pcb(i) = 0.
	1400	ema_pct(i) = 0.
	1401	! ema_workcbmf(i) = 0.
	1402	ENDDO
	1403	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
	1404	DO i = 1, klon
	1405	ibas_con(i) = 1
	1406	itop_con(i) = 1
	1407	ENDDO
	1408	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
	1409	!================================================================
	1410	!CR:04.12.07: initialisations poches froides
	1411	! Controle de ALE et ALP pour la fermeture convective (jyg)
	1412	if (iflag_wake>=1) then
	1413	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
	1414	,alp_bl_prescr, ale_bl_prescr)
	1415	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
	1416	! print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
	1417	endif
[973]	1418
[2469]	1419	! do i = 1,klon
	1420	! Ale_bl(i)=0.
	1421	! Alp_bl(i)=0.
	1422	! enddo
[1638]	1423
[2469]	1424	!===================================================================
	1425	!IM stations CFMIP
	1426	nCFMIP=npCFMIP
	1427	OPEN(98,file='npCFMIP_param.data',status='old', &
	1428	form='formatted',iostat=iostat)
	1429	if (iostat == 0) then
	1430	READ(98,*,end=998) nCFMIP
	1431	998 CONTINUE
	1432	CLOSE(98)
	1433	CONTINUE
	1434	IF(nCFMIP.GT.npCFMIP) THEN
	1435	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
	1436	call abort_physic("physiq", "", 1)
	1437	else
	1438	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
	1439	ENDIF
[1279]	1440
[2469]	1441	!
	1442	ALLOCATE(tabCFMIP(nCFMIP))
	1443	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
	1444	ALLOCATE(tabijGCM(nCFMIP))
	1445	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
	1446	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
	1447	!
	1448	! lecture des nCFMIP stations CFMIP, de leur numero
	1449	! et des coordonnees geographiques lonCFMIP, latCFMIP
	1450	!
	1451	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
	1452	lonCFMIP, latCFMIP)
	1453	!
	1454	! identification des
	1455	! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
	1456	! grille de LMDZ
	1457	! 2) indices points tabijGCM de la grille physique 1d sur
	1458	! klon points
	1459	! 3) indices iGCM, jGCM de la grille physique 2d
	1460	!
	1461	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
	1462	tabijGCM, lonGCM, latGCM, iGCM, jGCM)
	1463	!
	1464	else
	1465	ALLOCATE(tabijGCM(0))
	1466	ALLOCATE(lonGCM(0), latGCM(0))
	1467	ALLOCATE(iGCM(0), jGCM(0))
	1468	end if
	1469	else
	1470	ALLOCATE(tabijGCM(0))
	1471	ALLOCATE(lonGCM(0), latGCM(0))
	1472	ALLOCATE(iGCM(0), jGCM(0))
	1473	ENDIF
[878]	1474
[2469]	1475	DO i=1,klon
	1476	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1477	ENDDO
[1863]	1478
[2469]	1479	!34EK
	1480	IF (ok_orodr) THEN
[524]	1481
[2469]	1482	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1483	! FH sans doute a enlever de finitivement ou, si on le
	1484	! garde, l'activer justement quand ok_orodr = false.
	1485	! ce rugoro est utilise par la couche limite et fait double emploi
	1486	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
	1487	! DO i=1,klon
	1488	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
	1489	! ENDDO
	1490	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1491	IF (ok_strato) THEN
	1492	CALL SUGWD_strato(klon,klev,paprs,pplay)
	1493	ELSE
	1494	CALL SUGWD(klon,klev,paprs,pplay)
	1495	ENDIF
[1863]	1496
[2469]	1497	DO i=1,klon
	1498	zuthe(i)=0.
	1499	zvthe(i)=0.
	1500	if(zstd(i).gt.10.)then
	1501	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
	1502	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
	1503	endif
	1504	ENDDO
	1505	ENDIF
	1506	!
	1507	!
	1508	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
	1509	WRITE(lunout,*)'La frequence de lecture surface est de ', &
	1510	lmt_pas
	1511	!
	1512	capemaxcels = 't_max(X)'
	1513	t2mincels = 't_min(X)'
	1514	t2maxcels = 't_max(X)'
	1515	tinst = 'inst(X)'
	1516	tave = 'ave(X)'
	1517	!IM cf. AM 081204 BEG
	1518	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
	1519	!IM cf. AM 081204 END
	1520	!
	1521	!=============================================================
	1522	! Initialisation des sorties
	1523	!=============================================================
	1524
[524]	1525	#ifdef CPP_IOIPSL
	1526
[2469]	1527	!$OMP MASTER
	1528	! FH : if ok_sync=.true. , the time axis is written at each time step
	1529	! in the output files. Only at the end in the opposite case
	1530	ok_sync_omp=.false.
	1531	CALL getin('ok_sync',ok_sync_omp)
	1532	call phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
	1533	iGCM,jGCM,lonGCM,latGCM, &
	1534	jjmp1,nlevSTD,clevSTD,rlevSTD, dtime,ok_veget, &
	1535	type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
	1536	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
	1537	read_climoz, phys_out_filestations, &
	1538	new_aod, aerosol_couple, &
	1539	flag_aerosol_strat, pdtphys, paprs, pphis, &
	1540	pplay, lmax_th, ptconv, ptconvth, ivap, &
	1541	d_t, qx, d_qx, zmasse, ok_sync_omp)
	1542	!$OMP END MASTER
	1543	!$OMP BARRIER
	1544	ok_sync=ok_sync_omp
[909]	1545
[2469]	1546	freq_outNMC(1) = ecrit_files(7)
	1547	freq_outNMC(2) = ecrit_files(8)
	1548	freq_outNMC(3) = ecrit_files(9)
	1549	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
	1550	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
	1551	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
[524]	1552
[2469]	1553	include "ini_histday_seri.h"
[524]	1554
[2590]	1555	CALL ini_paramLMDZ_phy(dtime,nid_ctesGCM)
[524]	1556
[644]	1557	#endif
[2469]	1558	ecrit_reg = ecrit_reg * un_jour
	1559	ecrit_tra = ecrit_tra * un_jour
[1863]	1560
[2469]	1561	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
	1562	date0 = jD_ref
	1563	WRITE(,) 'physiq date0 : ',date0
	1564	!
	1565	!
	1566	!
	1567	! Prescrire l'ozone dans l'atmosphere
	1568	!
	1569	!
	1570	!c DO i = 1, klon
	1571	!c DO k = 1, klev
	1572	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
	1573	!c ENDDO
	1574	!c ENDDO
	1575	!
	1576	IF (type_trac == 'inca') THEN
[524]	1577	#ifdef INCA
[2469]	1578	CALL VTe(VTphysiq)
	1579	CALL VTb(VTinca)
	1580	calday = REAL(days_elapsed) + jH_cur
	1581	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
[959]	1582
[2469]	1583	CALL chemini( &
	1584	rg, &
	1585	ra, &
	1586	cell_area, &
	1587	latitude_deg, &
	1588	longitude_deg, &
	1589	presnivs, &
	1590	calday, &
	1591	klon, &
	1592	nqtot, &
[2566]	1593	nqo, &
[2469]	1594	pdtphys, &
	1595	annee_ref, &
[2914]	1596	year_cur, &
[2469]	1597	day_ref, &
	1598	day_ini, &
	1599	start_time, &
	1600	itau_phy, &
[2907]	1601	date0, &
[2469]	1602	io_lon, &
[3340]	1603	io_lat, &
[3433]	1604	chemistry_couple, &
	1605	init_source, &
	1606	init_tauinca, &
	1607	init_pizinca, &
	1608	init_cginca, &
	1609	init_ccminca)
[959]	1610
[3433]	1611	! initialisation des variables depuis le restart de inca
	1612	ccm(:,:,:) = init_ccminca
	1613	tau_aero(:,:,:,:) = init_tauinca
	1614	piz_aero(:,:,:,:) = init_pizinca
	1615	cg_aero(:,:,:,:) = init_cginca
	1616
[2469]	1617	CALL VTe(VTinca)
	1618	CALL VTb(VTphysiq)
[524]	1619	#endif
[2469]	1620	END IF
	1621	!
	1622	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1623	! Nouvelle initialisation pour le rayonnement RRTM
	1624	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[998]	1625
[2469]	1626	call iniradia(klon,klev,paprs(1,1:klev+1))
[998]	1627
[2469]	1628	!$omp single
	1629	if (read_climoz >= 1) then
	1630	call open_climoz(ncid_climoz, press_climoz)
	1631	END IF
	1632	!$omp end single
	1633	!
	1634	!IM betaCRF
	1635	pfree=70000. !Pa
	1636	beta_pbl=1.
	1637	beta_free=1.
	1638	lon1_beta=-180.
	1639	lon2_beta=+180.
	1640	lat1_beta=90.
	1641	lat2_beta=-90.
	1642	mskocean_beta=.FALSE.
[1539]	1643
[2469]	1644	!albedo SB >>>
	1645	select case(nsw)
	1646	case(2)
	1647	SFRWL(1)=0.45538747
	1648	SFRWL(2)=0.54461211
	1649	case(4)
	1650	SFRWL(1)=0.45538747
	1651	SFRWL(2)=0.32870591
	1652	SFRWL(3)=0.18568763
	1653	SFRWL(4)=3.02191470E-02
	1654	case(6)
	1655	SFRWL(1)=1.28432794E-03
	1656	SFRWL(2)=0.12304168
	1657	SFRWL(3)=0.33106142
	1658	SFRWL(4)=0.32870591
	1659	SFRWL(5)=0.18568763
	1660	SFRWL(6)=3.02191470E-02
	1661	end select
[2227]	1662
	1663
[2469]	1664	!albedo SB <<<
[2227]	1665
[2469]	1666	OPEN(99,file='beta_crf.data',status='old', &
	1667	form='formatted',err=9999)
	1668	READ(99,*,end=9998) pfree
	1669	READ(99,*,end=9998) beta_pbl
	1670	READ(99,*,end=9998) beta_free
	1671	READ(99,*,end=9998) lon1_beta
	1672	READ(99,*,end=9998) lon2_beta
	1673	READ(99,*,end=9998) lat1_beta
	1674	READ(99,*,end=9998) lat2_beta
	1675	READ(99,*,end=9998) mskocean_beta
	1676	9998 Continue
	1677	CLOSE(99)
	1678	9999 Continue
	1679	WRITE(,)'pfree=',pfree
	1680	WRITE(,)'beta_pbl=',beta_pbl
	1681	WRITE(,)'beta_free=',beta_free
	1682	WRITE(,)'lon1_beta=',lon1_beta
	1683	WRITE(,)'lon2_beta=',lon2_beta
	1684	WRITE(,)'lat1_beta=',lat1_beta
	1685	WRITE(,)'lat2_beta=',lat2_beta
	1686	WRITE(,)'mskocean_beta=',mskocean_beta
	1687	ENDIF
	1688	!
	1689	! ************** Fin de IF ( debut ) *************
	1690	!
	1691	!
	1692	! Incrementer le compteur de la physique
	1693	!
	1694	itap = itap + 1
	1695	!
	1696	!
	1697	! Update fraction of the sub-surfaces (pctsrf) and
	1698	! initialize, where a new fraction has appeared, all variables depending
	1699	! on the surface fraction.
	1700	!
	1701	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
	1702	pctsrf, fevap, z0m, z0h, agesno, &
	1703	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
[996]	1704
[2469]	1705	! Update time and other variables in Reprobus
	1706	IF (type_trac == 'repr') THEN
[1565]	1707	#ifdef REPROBUS
[2469]	1708	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
	1709	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
	1710	CALL Rtime(debut)
[1565]	1711	#endif
[2469]	1712	END IF
[1565]	1713
	1714
[2469]	1715	! Tendances bidons pour les processus qui n'affectent pas certaines
	1716	! variables.
	1717	du0(:,:)=0.
	1718	dv0(:,:)=0.
	1719	dt0 = 0.
	1720	dq0(:,:)=0.
	1721	dql0(:,:)=0.
	1722	dqi0(:,:)=0.
	1723	!
	1724	! Mettre a zero des variables de sortie (pour securite)
	1725	!
	1726	DO i = 1, klon
	1727	d_ps(i) = 0.0
	1728	ENDDO
	1729	DO k = 1, klev
	1730	DO i = 1, klon
	1731	d_t(i,k) = 0.0
	1732	d_u(i,k) = 0.0
	1733	d_v(i,k) = 0.0
	1734	ENDDO
	1735	ENDDO
	1736	DO iq = 1, nqtot
	1737	DO k = 1, klev
	1738	DO i = 1, klon
	1739	d_qx(i,k,iq) = 0.0
	1740	ENDDO
	1741	ENDDO
	1742	ENDDO
	1743	da(:,:)=0.
	1744	mp(:,:)=0.
	1745	phi(:,:,:)=0.
	1746	! RomP >>>
	1747	phi2(:,:,:)=0.
	1748	beta_prec_fisrt(:,:)=0.
	1749	beta_prec(:,:)=0.
	1750	epmlmMm(:,:,:)=0.
	1751	eplaMm(:,:)=0.
	1752	d1a(:,:)=0.
	1753	dam(:,:)=0.
	1754	pmflxr=0.
	1755	pmflxs=0.
	1756	! RomP <<<
[1742]	1757
[2469]	1758	!
	1759	! Ne pas affecter les valeurs entrees de u, v, h, et q
	1760	!
	1761	DO k = 1, klev
	1762	DO i = 1, klon
	1763	t_seri(i,k) = t(i,k)
	1764	u_seri(i,k) = u(i,k)
	1765	v_seri(i,k) = v(i,k)
	1766	q_seri(i,k) = qx(i,k,ivap)
	1767	ql_seri(i,k) = qx(i,k,iliq)
	1768	!CR: ATTENTION, on rajoute la variable glace
	1769	if (nqo.eq.2) then
	1770	qs_seri(i,k) = 0.
	1771	else if (nqo.eq.3) then
	1772	qs_seri(i,k) = qx(i,k,isol)
	1773	endif
	1774	ENDDO
	1775	ENDDO
[2476]	1776	!
	1777	!--OB mass fixer
	1778	IF (mass_fixer) THEN
	1779	!--store initial water burden
	1780	qql1(:)=0.0
[2499]	1781	DO k = 1, klev
	1782	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
[2476]	1783	ENDDO
	1784	ENDIF
	1785	!--fin mass fixer
	1786
[2469]	1787	tke0(:,:)=pbl_tke(:,:,is_ave)
	1788	!CR:Nombre de traceurs de l'eau: nqo
	1789	! IF (nqtot.GE.3) THEN
	1790	IF (nqtot.GE.(nqo+1)) THEN
	1791	! DO iq = 3, nqtot
	1792	DO iq = nqo+1, nqtot
	1793	DO k = 1, klev
	1794	DO i = 1, klon
	1795	! tr_seri(i,k,iq-2) = qx(i,k,iq)
	1796	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
	1797	ENDDO
	1798	ENDDO
	1799	ENDDO
	1800	ELSE
	1801	DO k = 1, klev
	1802	DO i = 1, klon
	1803	tr_seri(i,k,1) = 0.0
	1804	ENDDO
	1805	ENDDO
	1806	ENDIF
	1807	!
	1808	DO i = 1, klon
	1809	ztsol(i) = 0.
	1810	ENDDO
	1811	DO nsrf = 1, nbsrf
	1812	DO i = 1, klon
	1813	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
	1814	ENDDO
	1815	ENDDO
[2611]	1816	! Initialize variables used for diagnostic purpose
	1817	if (flag_inhib_tend .ne. 0) call init_cmp_seri
[2469]	1818	!IM
	1819	IF (ip_ebil_phy.ge.1) THEN
	1820	ztit='after dynamic'
	1821	CALL diagetpq(cell_area,ztit,ip_ebil_phy,1,1,dtime &
	1822	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	1823	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	1824	! Comme les tendances de la physique sont ajoute dans la dynamique,
	1825	! on devrait avoir que la variation d'entalpie par la dynamique
	1826	! est egale a la variation de la physique au pas de temps precedent.
	1827	! Donc la somme de ces 2 variations devrait etre nulle.
	1828	call diagphy(cell_area,ztit,ip_ebil_phy &
	1829	, zero_v, zero_v, zero_v, zero_v, zero_v &
	1830	, zero_v, zero_v, zero_v, ztsol &
	1831	, d_h_vcol+d_h_vcol_phy, d_qt, 0. &
	1832	, fs_bound, fq_bound )
	1833	END IF
[524]	1834
[2469]	1835	! Diagnostiquer la tendance dynamique
	1836	!
	1837	IF (ancien_ok) THEN
[2499]	1838	!
	1839	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/dtime
	1840	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/dtime
	1841	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/dtime
	1842	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/dtime
	1843	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/dtime
	1844	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/dtime
	1845	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
	1846	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/dtime
	1847	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
	1848	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/dtime
	1849	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
	1850	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/dtime
[2469]	1851	! !! RomP >>> td dyn traceur
[2499]	1852	IF (nqtot.GT.nqo) THEN ! jyg
[2469]	1853	DO iq = nqo+1, nqtot ! jyg
[2499]	1854	d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/dtime ! jyg
[2469]	1855	ENDDO
	1856	ENDIF
	1857	! !! RomP <<<
	1858	ELSE
[2499]	1859	d_u_dyn(:,:) = 0.0
	1860	d_v_dyn(:,:) = 0.0
	1861	d_t_dyn(:,:) = 0.0
	1862	d_q_dyn(:,:) = 0.0
	1863	d_ql_dyn(:,:) = 0.0
	1864	d_qs_dyn(:,:) = 0.0
	1865	d_q_dyn2d(:) = 0.0
	1866	d_ql_dyn2d(:) = 0.0
	1867	d_qs_dyn2d(:) = 0.0
[2469]	1868	! !! RomP >>> td dyn traceur
[2499]	1869	IF (nqtot.GT.nqo) THEN ! jyg
	1870	DO iq = nqo+1, nqtot ! jyg
	1871	d_tr_dyn(:,:,iq-nqo)= 0.0 ! jyg
[2469]	1872	ENDDO
	1873	ENDIF
	1874	! !! RomP <<<
	1875	ancien_ok = .TRUE.
	1876	ENDIF
	1877	!
	1878	! Ajouter le geopotentiel du sol:
	1879	!
	1880	DO k = 1, klev
	1881	DO i = 1, klon
	1882	zphi(i,k) = pphi(i,k) + pphis(i)
	1883	ENDDO
	1884	ENDDO
	1885	!
	1886	! Verifier les temperatures
	1887	!
	1888	!IM BEG
	1889	IF (check) THEN
	1890	amn=MIN(ftsol(1,is_ter),1000.)
	1891	amx=MAX(ftsol(1,is_ter),-1000.)
	1892	DO i=2, klon
	1893	amn=MIN(ftsol(i,is_ter),amn)
	1894	amx=MAX(ftsol(i,is_ter),amx)
	1895	ENDDO
	1896	!
	1897	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
	1898	ENDIF !(check) THEN
	1899	!IM END
	1900	!
	1901	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
	1902	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
[2235]	1903
[2469]	1904	!
	1905	!IM BEG
	1906	IF (check) THEN
	1907	amn=MIN(ftsol(1,is_ter),1000.)
	1908	amx=MAX(ftsol(1,is_ter),-1000.)
	1909	DO i=2, klon
	1910	amn=MIN(ftsol(i,is_ter),amn)
	1911	amx=MAX(ftsol(i,is_ter),amx)
	1912	ENDDO
	1913	!
	1914	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
	1915	ENDIF !(check) THEN
	1916	!IM END
	1917	!
	1918	! Mettre en action les conditions aux limites (albedo, sst, etc.).
	1919	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
	1920	!
	1921	if (read_climoz >= 1) then
	1922	! Ozone from a file
	1923	! Update required ozone index:
	1924	ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
	1925	if (ro3i == 361) ro3i = 360
	1926	! (This should never occur, except perhaps because of roundup
	1927	! error. See documentation.)
	1928	if (ro3i /= co3i) then
	1929	! Update ozone field:
	1930	if (read_climoz == 1) then
	1931	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
	1932	press_in_edg=press_climoz, paprs=paprs, v3=wo)
	1933	else
	1934	! read_climoz == 2
	1935	call regr_pr_av(ncid_climoz, (/"tro3 ", &
	1936	"tro3_daylight"/), julien=ro3i, press_in_edg=press_climoz, &
	1937	paprs=paprs, v3=wo)
	1938	end if
	1939	! Convert from mole fraction of ozone to column density of ozone in a
	1940	! cell, in kDU:
	1941	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
	1942	* zmasse / dobson_u / 1e3
	1943	! (By regridding ozone values for LMDZ only once every 360th of
	1944	! year, we have already neglected the variation of pressure in one
	1945	! 360th of year. So do not recompute "wo" at each time step even if
	1946	! "zmasse" changes a little.)
	1947	co3i = ro3i
	1948	end if
	1949	ELSEIF (MOD(itap-1,lmt_pas) == 0) THEN
	1950	! Once per day, update ozone from Royer:
[1955]	1951
[2469]	1952	IF (solarlong0<-999.) then
	1953	! Generic case with evolvoing season
	1954	zzz=real(days_elapsed+1)
	1955	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
	1956	! Particular case with annual mean insolation
	1957	zzz=real(90) ! could be revisited
	1958	IF (read_climoz/=-1) THEN
	1959	abort_message ='read_climoz=-1 is recommended when ' &
	1960	// 'solarlong0=1000.'
	1961	CALL abort_physic (modname,abort_message,1)
	1962	ENDIF
	1963	ELSE
	1964	! Case where the season is imposed with solarlong0
	1965	zzz=real(90) ! could be revisited
	1966	ENDIF
	1967	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
	1968	ENDIF
	1969	!
	1970	! Re-evaporer l'eau liquide nuageuse
	1971	!
	1972	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
	1973	DO i = 1, klon
	1974	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1975	!jyg<
	1976	! Attention : Arnaud a propose des formules completement differentes
	1977	! A verifier !!!
	1978	zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
	1979	IF (iflag_ice_thermo .EQ. 0) THEN
	1980	zlsdcp=zlvdcp
	1981	ENDIF
	1982	!>jyg
[1849]	1983
[2469]	1984	if (iflag_ice_thermo.eq.0) then
	1985	!pas necessaire a priori
[2086]	1986
[2469]	1987	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
	1988	zb = MAX(0.0,ql_seri(i,k))
	1989	za = - MAX(0.0,ql_seri(i,k)) &
	1990	* (zlvdcp(1.-zdelta)+zlsdcpzdelta)
	1991	t_seri(i,k) = t_seri(i,k) + za
	1992	q_seri(i,k) = q_seri(i,k) + zb
	1993	ql_seri(i,k) = 0.0
	1994	d_t_eva(i,k) = za
	1995	d_q_eva(i,k) = zb
[2086]	1996
[2469]	1997	else
[2086]	1998
[2469]	1999	!CR: on r\'e-\'evapore eau liquide et glace
[2086]	2000
[2469]	2001	! zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
	2002	! zb = MAX(0.0,ql_seri(i,k))
	2003	! za = - MAX(0.0,ql_seri(i,k)) &
	2004	! * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
	2005	zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
	2006	za = - MAX(0.0,ql_seri(i,k))*zlvdcp &
	2007	- MAX(0.0,qs_seri(i,k))*zlsdcp
	2008	t_seri(i,k) = t_seri(i,k) + za
	2009	q_seri(i,k) = q_seri(i,k) + zb
	2010	ql_seri(i,k) = 0.0
	2011	!on \'evapore la glace
	2012	qs_seri(i,k) = 0.0
	2013	d_t_eva(i,k) = za
	2014	d_q_eva(i,k) = zb
	2015	endif
[782]	2016
[2469]	2017	ENDDO
	2018	ENDDO
	2019	!IM
	2020	IF (ip_ebil_phy.ge.2) THEN
	2021	ztit='after reevap'
	2022	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,1,dtime &
	2023	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2024	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2025	call diagphy(cell_area,ztit,ip_ebil_phy &
	2026	, zero_v, zero_v, zero_v, zero_v, zero_v &
	2027	, zero_v, zero_v, zero_v, ztsol &
	2028	, d_h_vcol, d_qt, d_ec &
	2029	, fs_bound, fq_bound )
	2030	!
	2031	END IF
[883]	2032
[2469]	2033	!
	2034	!=========================================================================
	2035	! Calculs de l'orbite.
	2036	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
	2037	! doit donc etre plac\'e avant radlwsw et pbl_surface
[883]	2038
[2469]	2039	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2040	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
	2041	day_since_equinox = (jD_cur + jH_cur) - jD_eq
	2042	!
	2043	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
	2044	! solarlong0
	2045	if (solarlong0<-999.) then
	2046	if (new_orbit) then
	2047	! calcul selon la routine utilisee pour les planetes
	2048	call solarlong(day_since_equinox, zlongi, dist)
	2049	else
	2050	! calcul selon la routine utilisee pour l'AR4
	2051	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
	2052	endif
	2053	else
	2054	zlongi=solarlong0 ! longitude solaire vraie
	2055	dist=1. ! distance au soleil / moyenne
	2056	endif
	2057	if(prt_level.ge.1) &
	2058	write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
[1529]	2059
[524]	2060
[2469]	2061	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2062	! Calcul de l'ensoleillement :
	2063	! ============================
	2064	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
	2065	! l'annee a partir d'une formule analytique.
	2066	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
	2067	! non nul aux poles.
	2068	IF (abs(solarlong0-1000.)<1.e-4) then
	2069	call zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
	2070	latitude_deg,longitude_deg,rmu0,fract)
	2071	JrNt = 1.0
	2072	ELSE
	2073	! recode par Olivier Boucher en sept 2015
	2074	SELECT CASE (iflag_cycle_diurne)
	2075	CASE(0)
	2076	! Sans cycle diurne
	2077	CALL angle(zlongi, latitude_deg, fract, rmu0)
	2078	swradcorr = 1.0
	2079	JrNt = 1.0
	2080	zrmu0 = rmu0
	2081	CASE(1)
	2082	! Avec cycle diurne sans application des poids
	2083	! bit comparable a l ancienne formulation cycle_diurne=true
	2084	! on integre entre gmtime et gmtime+radpas
	2085	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
	2086	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
	2087	latitude_deg,longitude_deg,rmu0,fract)
	2088	zrmu0 = rmu0
	2089	swradcorr = 1.0
	2090	! Calcul du flag jour-nuit
	2091	JrNt = 0.0
	2092	WHERE (fract.GT.0.0) JrNt = 1.0
	2093	CASE(2)
	2094	! Avec cycle diurne sans application des poids
	2095	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
	2096	! Comme cette routine est appele a tous les pas de temps de
	2097	! la physique meme si le rayonnement n'est pas appele je
	2098	! remonte en arriere les radpas-1 pas de temps
	2099	! suivant. Petite ruse avec MOD pour prendre en compte le
	2100	! premier pas de temps de la physique pendant lequel
	2101	! itaprad=0
	2102	zdtime1=dtime*REAL(-MOD(itaprad,radpas)-1)
	2103	zdtime2=dtime*REAL(radpas-MOD(itaprad,radpas)-1)
	2104	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
	2105	latitude_deg,longitude_deg,rmu0,fract)
	2106	!
	2107	! Calcul des poids
	2108	!
	2109	zdtime1=-dtime !--on corrige le rayonnement pour representer le
	2110	zdtime2=0.0 !--pas de temps de la physique qui se termine
	2111	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
	2112	latitude_deg,longitude_deg,zrmu0,zfract)
	2113	swradcorr = 0.0
	2114	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
	2115	swradcorr=zfract/fract*zrmu0/rmu0
	2116	! Calcul du flag jour-nuit
	2117	JrNt = 0.0
	2118	WHERE (zfract.GT.0.0) JrNt = 1.0
	2119	END SELECT
	2120	ENDIF
[782]	2121
[2469]	2122	if (mydebug) then
	2123	call writefield_phy('u_seri',u_seri,nbp_lev)
	2124	call writefield_phy('v_seri',v_seri,nbp_lev)
	2125	call writefield_phy('t_seri',t_seri,nbp_lev)
	2126	call writefield_phy('q_seri',q_seri,nbp_lev)
	2127	endif
[883]	2128
[2469]	2129	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	2130	! Appel au pbl_surface : Planetary Boudary Layer et Surface
	2131	! Cela implique tous les interactions des sous-surfaces et la
	2132	! partie diffusion turbulent du couche limit.
	2133	!
	2134	! Certains varibales de sorties de pbl_surface sont utiliser que pour
	2135	! ecriture des fihiers hist_XXXX.nc, ces sont :
	2136	! qsol, zq2m, s_pblh, s_lcl,
	2137	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
	2138	! s_therm, s_trmb1, s_trmb2, s_trmb3,
	2139	! zu10m, zv10m, fder,
	2140	! zxqsurf, rh2m, zxfluxu, zxfluxv,
	2141	! frugs, agesno, fsollw, fsolsw,
	2142	! d_ts, fevap, fluxlat, t2m,
	2143	! wfbils, wfbilo, fluxt, fluxu, fluxv,
	2144	!
	2145	! Certains ne sont pas utiliser du tout :
	2146	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
	2147	!
[1724]	2148
[2469]	2149	! Calcul de l'humidite de saturation au niveau du sol
[1724]	2150
	2151
[996]	2152
[2469]	2153	if (iflag_pbl/=0) then
[2240]	2154
[2469]	2155	!jyg+nrlmd<
	2156	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
	2157	print *,'debut du splitting de la PBL'
	2158	ENDIF
	2159	! !!
	2160	!=================================================================
	2161	! PROVISOIRE : DECOUPLAGE PBL/WAKE
	2162	! --------------------------------
	2163	!
	2164	!! wake_deltat_sav(:,:)=wake_deltat(:,:)
	2165	!! wake_deltaq_sav(:,:)=wake_deltaq(:,:)
	2166	!! wake_deltat(:,:)=0.
	2167	!! wake_deltaq(:,:)=0.
	2168	!=================================================================
	2169	!>jyg+nrlmd
	2170	!
	2171	!-------gustiness calculation-------!
	2172	IF (iflag_gusts==0) THEN
	2173	gustiness(1:klon)=0
	2174	ELSE IF (iflag_gusts==1) THEN
	2175	do i = 1, klon
	2176	gustiness(i)=f_gust_blale_bl(i)+f_gust_wkale_wake(i)
	2177	enddo
	2178	! ELSE IF (iflag_gusts==2) THEN
	2179	! do i = 1, klon
	2180	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
	2181	! *ale_wake(i) !! need to make sigma_wk accessible here
	2182	! enddo
	2183	! ELSE IF (iflag_gusts==3) THEN
	2184	! do i = 1, klon
	2185	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
	2186	! enddo
	2187	ENDIF
[2278]	2188
	2189
[1067]	2190
[2469]	2191	CALL pbl_surface( &
	2192	dtime, date0, itap, days_elapsed+1, &
	2193	debut, lafin, &
	2194	longitude_deg, latitude_deg, rugoro, zrmu0, &
	2195	zsig, sollwdown, pphi, cldt, &
	2196	rain_fall, snow_fall, solsw, sollw, &
	2197	gustiness, &
	2198	t_seri, q_seri, u_seri, v_seri, &
	2199	!nrlmd+jyg<
	2200	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
	2201	!>nrlmd+jyg
	2202	pplay, paprs, pctsrf, &
	2203	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
	2204	!albedo SB <<<
	2205	cdragh, cdragm, u1, v1, &
	2206	!albedo SB >>>
	2207	! albsol1, albsol2, sens, evap, &
	2208	albsol_dir, albsol_dif, sens, evap, &
	2209	!albedo SB <<<
	2210	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
	2211	zxtsol, zxfluxlat, zt2m, qsat2m, &
	2212	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
	2213	!nrlmd<
	2214	!jyg<
	2215	d_t_vdf_w, d_q_vdf_w, &
	2216	d_t_vdf_x, d_q_vdf_x, &
	2217	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
	2218	!>jyg
	2219	delta_tsurf,wake_dens, &
	2220	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
	2221	kh,kh_x,kh_w, &
	2222	!>nrlmd
	2223	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
	2224	slab_wfbils, &
	2225	qsol, zq2m, s_pblh, s_lcl, &
	2226	!jyg<
	2227	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
	2228	!>jyg
	2229	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
	2230	s_therm, s_trmb1, s_trmb2, s_trmb3, &
	2231	zustar, zu10m, zv10m, fder, &
	2232	zxqsurf, rh2m, zxfluxu, zxfluxv, &
	2233	z0m, z0h, agesno, fsollw, fsolsw, &
	2234	d_ts, fevap, fluxlat, t2m, &
	2235	wfbils, wfbilo, fluxt, fluxu, fluxv, &
	2236	dsens, devap, zxsnow, &
	2237	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
	2238	!nrlmd+jyg<
	2239	wake_delta_pbl_TKE &
	2240	!>nrlmd+jyg
	2241	)
	2242	!
	2243	!=================================================================
	2244	! PROVISOIRE : DECOUPLAGE PBL/WAKE
	2245	! --------------------------------
	2246	!
	2247	!! wake_deltat(:,:)=wake_deltat_sav(:,:)
	2248	!! wake_deltaq(:,:)=wake_deltaq_sav(:,:)
	2249	!=================================================================
	2250	!
	2251	! Add turbulent diffusion tendency to the wake difference variables
	2252	IF (mod(iflag_pbl_split,2) .NE. 0) THEN
	2253	wake_deltat(:,:) = wake_deltat(:,:) + (d_t_vdf_w(:,:)-d_t_vdf_x(:,:))
	2254	wake_deltaq(:,:) = wake_deltaq(:,:) + (d_q_vdf_w(:,:)-d_q_vdf_x(:,:))
	2255	ENDIF
[1624]	2256
[766]	2257
[2469]	2258	!---------------------------------------------------------------------
	2259	! ajout des tendances de la diffusion turbulente
	2260	IF (klon_glo==1) THEN
	2261	CALL add_pbl_tend &
	2262	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
[2611]	2263	'vdf',abortphy,flag_inhib_tend)
[2469]	2264	ELSE
	2265	CALL add_phys_tend &
	2266	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
[2611]	2267	'vdf',abortphy,flag_inhib_tend)
[2469]	2268	ENDIF
	2269	!--------------------------------------------------------------------
[766]	2270
[2469]	2271	if (mydebug) then
	2272	call writefield_phy('u_seri',u_seri,nbp_lev)
	2273	call writefield_phy('v_seri',v_seri,nbp_lev)
	2274	call writefield_phy('t_seri',t_seri,nbp_lev)
	2275	call writefield_phy('q_seri',q_seri,nbp_lev)
	2276	endif
[2227]	2277
	2278
[2469]	2279	!albedo SB >>>
	2280	albsol1=0.
	2281	albsol2=0.
	2282	falb1=0.
	2283	falb2=0.
	2284	select case(nsw)
	2285	case(2)
	2286	albsol1=albsol_dir(:,1)
	2287	albsol2=albsol_dir(:,2)
	2288	falb1=falb_dir(:,1,:)
	2289	falb2=falb_dir(:,2,:)
	2290	case(4)
	2291	albsol1=albsol_dir(:,1)
	2292	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3) &
	2293	+albsol_dir(:,4)*SFRWL(4)
	2294	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
	2295	falb1=falb_dir(:,1,:)
	2296	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3) &
	2297	+falb_dir(:,4,:)*SFRWL(4)
	2298	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
	2299	case(6)
	2300	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2) &
	2301	+albsol_dir(:,3)*SFRWL(3)
	2302	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
	2303	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5) &
	2304	+albsol_dir(:,6)*SFRWL(6)
	2305	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
	2306	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2) &
	2307	+falb_dir(:,3,:)*SFRWL(3)
	2308	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
	2309	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5) &
	2310	+falb_dir(:,6,:)*SFRWL(6)
	2311	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
	2312	end select
	2313	!albedo SB <<<
[2227]	2314
[766]	2315
[2469]	2316	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
	2317	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
[1724]	2318
[524]	2319
[2469]	2320	IF (ip_ebil_phy.ge.2) THEN
	2321	ztit='after surface_main'
	2322	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
	2323	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2324	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2325	call diagphy(cell_area,ztit,ip_ebil_phy &
	2326	, zero_v, zero_v, zero_v, zero_v, sens &
	2327	, evap , zero_v, zero_v, ztsol &
	2328	, d_h_vcol, d_qt, d_ec &
	2329	, fs_bound, fq_bound )
	2330	END IF
[881]	2331
[2469]	2332	ENDIF
	2333	! =================================================================== c
	2334	! Calcul de Qsat
[881]	2335
[2469]	2336	DO k = 1, klev
	2337	DO i = 1, klon
	2338	zx_t = t_seri(i,k)
	2339	IF (thermcep) THEN
	2340	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	2341	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	2342	zx_qs = MIN(0.5,zx_qs)
	2343	zcor = 1./(1.-retv*zx_qs)
	2344	zx_qs = zx_qs*zcor
	2345	ELSE
	2346	!! IF (zx_t.LT.t_coup) THEN !jyg
	2347	IF (zx_t.LT.rtt) THEN !jyg
	2348	zx_qs = qsats(zx_t)/pplay(i,k)
	2349	ELSE
	2350	zx_qs = qsatl(zx_t)/pplay(i,k)
	2351	ENDIF
	2352	ENDIF
	2353	zqsat(i,k)=zx_qs
	2354	ENDDO
	2355	ENDDO
[959]	2356
[2469]	2357	if (prt_level.ge.1) then
	2358	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
	2359	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
	2360	endif
	2361	!
	2362	! Appeler la convection (au choix)
	2363	!
	2364	DO k = 1, klev
	2365	DO i = 1, klon
	2366	conv_q(i,k) = d_q_dyn(i,k) &
	2367	+ d_q_vdf(i,k)/dtime
	2368	conv_t(i,k) = d_t_dyn(i,k) &
	2369	+ d_t_vdf(i,k)/dtime
	2370	ENDDO
	2371	ENDDO
	2372	IF (check) THEN
	2373	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
	2374	WRITE(lunout,*) "avantcon=", za
	2375	ENDIF
	2376	zx_ajustq = .FALSE.
	2377	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
	2378	IF (zx_ajustq) THEN
	2379	DO i = 1, klon
	2380	z_avant(i) = 0.0
	2381	ENDDO
	2382	DO k = 1, klev
	2383	DO i = 1, klon
	2384	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
	2385	*(paprs(i,k)-paprs(i,k+1))/RG
	2386	ENDDO
	2387	ENDDO
	2388	ENDIF
[959]	2389
[2469]	2390	! Calcule de vitesse verticale a partir de flux de masse verticale
	2391	DO k = 1, klev
	2392	DO i = 1, klon
	2393	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
	2394	END DO
	2395	END DO
	2396	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
	2397	omega(igout, :)
[1015]	2398
[2469]	2399	IF (iflag_con.EQ.1) THEN
	2400	abort_message ='reactiver le call conlmd dans physiq.F'
	2401	CALL abort_physic (modname,abort_message,1)
	2402	! CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
	2403	! . d_t_con, d_q_con,
	2404	! . rain_con, snow_con, ibas_con, itop_con)
	2405	ELSE IF (iflag_con.EQ.2) THEN
	2406	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
	2407	conv_t, conv_q, -evap, omega, &
	2408	d_t_con, d_q_con, rain_con, snow_con, &
	2409	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
	2410	kcbot, kctop, kdtop, pmflxr, pmflxs)
	2411	d_u_con = 0.
	2412	d_v_con = 0.
[879]	2413
[2469]	2414	WHERE (rain_con < 0.) rain_con = 0.
	2415	WHERE (snow_con < 0.) snow_con = 0.
	2416	DO i = 1, klon
	2417	ibas_con(i) = klev+1 - kcbot(i)
	2418	itop_con(i) = klev+1 - kctop(i)
	2419	ENDDO
	2420	ELSE IF (iflag_con.GE.3) THEN
	2421	! nb of tracers for the KE convection:
	2422	! MAF la partie traceurs est faite dans phytrac
	2423	! on met ntra=1 pour limiter les appels mais on peut
	2424	! supprimer les calculs / ftra.
	2425	ntra = 1
	2426
	2427	!=======================================================================
	2428	!ajout pour la parametrisation des poches froides: calcul de
	2429	!t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
	2430	do k=1,klev
	2431	do i=1,klon
	2432	if (iflag_wake>=1) then
	2433	t_wake(i,k) = t_seri(i,k) &
	2434	+(1-wake_s(i))*wake_deltat(i,k)
	2435	q_wake(i,k) = q_seri(i,k) &
	2436	+(1-wake_s(i))*wake_deltaq(i,k)
	2437	t_undi(i,k) = t_seri(i,k) &
	2438	-wake_s(i)*wake_deltat(i,k)
	2439	q_undi(i,k) = q_seri(i,k) &
	2440	-wake_s(i)*wake_deltaq(i,k)
	2441	else
	2442	t_wake(i,k) = t_seri(i,k)
	2443	q_wake(i,k) = q_seri(i,k)
	2444	t_undi(i,k) = t_seri(i,k)
	2445	q_undi(i,k) = q_seri(i,k)
	2446	endif
	2447	enddo
	2448	enddo
	2449	!
	2450	!jyg<
	2451	! Perform dry adiabatic adjustment on wake profile
	2452	! The corresponding tendencies are added to the convective tendencies
	2453	! after the call to the convective scheme.
	2454	IF (iflag_wake>=1) then
	2455	IF (ok_adjwk) THEN
	2456	limbas(:) = 1
	2457	CALL ajsec(paprs, pplay, t_wake, q_wake, limbas, &
[2309]	2458	d_t_adjwk, d_q_adjwk)
[2318]	2459	ENDIF
[2469]	2460	!
	2461	DO k=1,klev
	2462	DO i=1,klon
	2463	IF (wake_s(i) .GT. 1.e-3) THEN
	2464	t_wake(i,k) = t_wake(i,k) + d_t_adjwk(i,k)
	2465	q_wake(i,k) = q_wake(i,k) + d_q_adjwk(i,k)
	2466	wake_deltat(i,k) = wake_deltat(i,k) + d_t_adjwk(i,k)
	2467	wake_deltaq(i,k) = wake_deltaq(i,k) + d_q_adjwk(i,k)
	2468	ENDIF
	2469	ENDDO
	2470	ENDDO
	2471	ENDIF ! (iflag_wake>=1)
	2472	!>jyg
	2473	!
[2513]	2474	!jyg<
	2475	CALL alpale( debut, itap, dtime, paprs, omega, t_seri, &
	2476	alp_offset, it_wape_prescr, wape_prescr, fip_prescr, &
	2477	ale_bl_prescr, alp_bl_prescr, &
	2478	wake_pe, wake_fip, &
	2479	Ale_bl, Ale_bl_trig, Alp_bl, &
[2554]	2480	Ale, Alp , Ale_wake, Alp_wake)
[2513]	2481	!>jyg
	2482	!
[2469]	2483	! sb, oct02:
	2484	! Schema de convection modularise et vectorise:
	2485	! (driver commun aux versions 3 et 4)
	2486	!
	2487	IF (ok_cvl) THEN ! new driver for convectL
	2488	!
	2489	!jyg<
	2490	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	2491	! Calculate the upmost level of deep convection loops: k_upper_cv
	2492	! (near 22 km)
	2493	izero = klon/2+1/klon
	2494	k_upper_cv = klev
	2495	DO k = klev,1,-1
	2496	IF (pphi(izero,k) > 22.e4) k_upper_cv = k
	2497	ENDDO
	2498	IF (prt_level .ge. 5) THEN
	2499	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
	2500	k_upper_cv
	2501	ENDIF
	2502	!
	2503	!>jyg
	2504	IF (type_trac == 'repr') THEN
	2505	nbtr_tmp=ntra
	2506	ELSE
	2507	nbtr_tmp=nbtr
	2508	END IF
	2509	!jyg iflag_con est dans clesphys
	2510	!c CALL concvl (iflag_con,iflag_clos,
	2511	CALL concvl (iflag_clos, &
	2512	dtime, paprs, pplay, k_upper_cv, t_undi,q_undi, &
	2513	t_wake,q_wake,wake_s, &
	2514	u_seri,v_seri,tr_seri,nbtr_tmp, &
	2515	ALE,ALP, &
	2516	sig1,w01, &
	2517	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
	2518	rain_con, snow_con, ibas_con, itop_con, sigd, &
	2519	ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
	2520	Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
	2521	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
	2522	! RomP >>>
	2523	!! . pmflxr,pmflxs,da,phi,mp,
	2524	!! . ftd,fqd,lalim_conv,wght_th)
	2525	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
	2526	ftd,fqd,lalim_conv,wght_th, &
	2527	ev, ep,epmlmMm,eplaMm, &
	2528	wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
[2481]	2529	tau_cld_cv,coefw_cld_cv,epmax_diag)
[2469]	2530	! RomP <<<
[619]	2531
[2469]	2532	!IM begin
	2533	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
	2534	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
	2535	!IM end
	2536	!IM cf. FH
	2537	clwcon0=qcondc
	2538	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
[524]	2539
[2469]	2540	do i = 1, klon
	2541	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
	2542	enddo
	2543	!
	2544	!jyg<
	2545	! Add the tendency due to the dry adjustment of the wake profile
	2546	IF (iflag_wake>=1) THEN
	2547	DO k=1,klev
	2548	DO i=1,klon
	2549	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
	2550	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
	2551	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
	2552	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
	2553	ENDDO
	2554	ENDDO
	2555	ENDIF
	2556	!>jyg
	2557	!
	2558	ELSE ! ok_cvl
[1412]	2559
[2469]	2560	! MAF conema3 ne contient pas les traceurs
	2561	CALL conema3 (dtime, &
	2562	paprs,pplay,t_seri,q_seri, &
	2563	u_seri,v_seri,tr_seri,ntra, &
	2564	sig1,w01, &
	2565	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
	2566	rain_con, snow_con, ibas_con, itop_con, &
	2567	upwd,dnwd,dnwd0,bas,top, &
	2568	Ma,cape,tvp,rflag, &
	2569	pbase &
	2570	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
	2571	,clwcon0)
[524]	2572
[2469]	2573	ENDIF ! ok_cvl
[524]	2574
[2469]	2575	!
	2576	! Correction precip
	2577	rain_con = rain_con * cvl_corr
	2578	snow_con = snow_con * cvl_corr
	2579	!
[766]	2580
[2469]	2581	IF (.NOT. ok_gust) THEN
	2582	do i = 1, klon
	2583	wd(i)=0.0
	2584	enddo
	2585	ENDIF
[524]	2586
[2469]	2587	! =================================================================== c
	2588	! Calcul des proprietes des nuages convectifs
	2589	!
[524]	2590
[2469]	2591	! calcul des proprietes des nuages convectifs
	2592	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
	2593	IF (iflag_cld_cv == 0) THEN
	2594	call clouds_gno &
	2595	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
	2596	ELSE
	2597	call clouds_bigauss &
	2598	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
	2599	ENDIF
[524]	2600
[2205]	2601
[2469]	2602	! =================================================================== c
[524]	2603
[2469]	2604	DO i = 1, klon
	2605	itop_con(i) = min(max(itop_con(i),1),klev)
	2606	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
	2607	ENDDO
[1428]	2608
[2469]	2609	DO i = 1, klon
	2610	ema_pcb(i) = paprs(i,ibas_con(i))
	2611	ENDDO
	2612	DO i = 1, klon
	2613	! L'idicage de itop_con peut cacher un pb potentiel
	2614	! FH sous la dictee de JYG, CR
	2615	ema_pct(i) = paprs(i,itop_con(i)+1)
[879]	2616
[2469]	2617	if (itop_con(i).gt.klev-3) then
	2618	if(prt_level >= 9) then
	2619	write(lunout,*)'La convection monte trop haut '
	2620	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
	2621	endif
	2622	endif
	2623	ENDDO
	2624	ELSE IF (iflag_con.eq.0) THEN
	2625	write(lunout,*) 'On n appelle pas la convection'
	2626	clwcon0=0.
	2627	rnebcon0=0.
	2628	d_t_con=0.
	2629	d_q_con=0.
	2630	d_u_con=0.
	2631	d_v_con=0.
	2632	rain_con=0.
	2633	snow_con=0.
	2634	bas=1
	2635	top=1
	2636	ELSE
	2637	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
	2638	call abort_physic("physiq", "", 1)
	2639	ENDIF
[524]	2640
[2469]	2641	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
	2642	! . d_u_con, d_v_con)
[524]	2643
[2469]	2644	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
[2611]	2645	'convection',abortphy,flag_inhib_tend)
[2235]	2646
[2469]	2647	!-------------------------------------------------------------------------
[766]	2648
[2469]	2649	if (mydebug) then
	2650	call writefield_phy('u_seri',u_seri,nbp_lev)
	2651	call writefield_phy('v_seri',v_seri,nbp_lev)
	2652	call writefield_phy('t_seri',t_seri,nbp_lev)
	2653	call writefield_phy('q_seri',q_seri,nbp_lev)
	2654	endif
[766]	2655
[2469]	2656	!IM
	2657	IF (ip_ebil_phy.ge.2) THEN
	2658	ztit='after convect'
	2659	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
	2660	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2661	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2662	call diagphy(cell_area,ztit,ip_ebil_phy &
	2663	, zero_v, zero_v, zero_v, zero_v, zero_v &
	2664	, zero_v, rain_con, snow_con, ztsol &
	2665	, d_h_vcol, d_qt, d_ec &
	2666	, fs_bound, fq_bound )
	2667	END IF
	2668	!
	2669	IF (check) THEN
	2670	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
	2671	WRITE(lunout,*)"aprescon=", za
	2672	zx_t = 0.0
	2673	za = 0.0
	2674	DO i = 1, klon
	2675	za = za + cell_area(i)/REAL(klon)
	2676	zx_t = zx_t + (rain_con(i)+ &
	2677	snow_con(i))*cell_area(i)/REAL(klon)
	2678	ENDDO
	2679	zx_t = zx_t/za*dtime
	2680	WRITE(lunout,*)"Precip=", zx_t
	2681	ENDIF
	2682	IF (zx_ajustq) THEN
	2683	DO i = 1, klon
	2684	z_apres(i) = 0.0
	2685	ENDDO
	2686	DO k = 1, klev
	2687	DO i = 1, klon
	2688	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
	2689	*(paprs(i,k)-paprs(i,k+1))/RG
	2690	ENDDO
	2691	ENDDO
	2692	DO i = 1, klon
	2693	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
	2694	/z_apres(i)
	2695	ENDDO
	2696	DO k = 1, klev
	2697	DO i = 1, klon
	2698	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
	2699	z_factor(i).LT.(1.0-1.0E-08)) THEN
	2700	q_seri(i,k) = q_seri(i,k) * z_factor(i)
	2701	ENDIF
	2702	ENDDO
	2703	ENDDO
	2704	ENDIF
	2705	zx_ajustq=.FALSE.
[879]	2706
[2469]	2707	!
	2708	!==========================================================================
	2709	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
	2710	!pour la couche limite diffuse pour l instant
	2711	!
	2712	!
	2713	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
	2714	! il faut rajouter cette tendance calcul\'ee hors des poches
	2715	! froides
	2716	!
	2717	if (iflag_wake>=1) then
[2159]	2718	DO k=1,klev
[2469]	2719	DO i=1,klon
	2720	dt_dwn(i,k) = ftd(i,k)
	2721	dq_dwn(i,k) = fqd(i,k)
	2722	M_dwn(i,k) = dnwd0(i,k)
	2723	M_up(i,k) = upwd(i,k)
	2724	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
	2725	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
	2726	ENDDO
[2159]	2727	ENDDO
[2469]	2728	!nrlmd+jyg<
[2159]	2729	DO k=1,klev
[2469]	2730	DO i=1,klon
	2731	wdt_PBL(i,k) = 0.
	2732	wdq_PBL(i,k) = 0.
	2733	udt_PBL(i,k) = 0.
	2734	udq_PBL(i,k) = 0.
	2735	ENDDO
[2159]	2736	ENDDO
[2469]	2737	!
	2738	IF (mod(iflag_pbl_split,2) .EQ. 1) THEN
	2739	DO k=1,klev
	2740	DO i=1,klon
	2741	wdt_PBL(i,k) = wdt_PBL(i,k) + d_t_vdf_w(i,k)/dtime
	2742	wdq_PBL(i,k) = wdq_PBL(i,k) + d_q_vdf_w(i,k)/dtime
	2743	udt_PBL(i,k) = udt_PBL(i,k) + d_t_vdf_x(i,k)/dtime
	2744	udq_PBL(i,k) = udq_PBL(i,k) + d_q_vdf_x(i,k)/dtime
	2745	!! dt_dwn(i,k) = dt_dwn(i,k) + d_t_vdf_w(i,k)/dtime
	2746	!! dq_dwn(i,k) = dq_dwn(i,k) + d_q_vdf_w(i,k)/dtime
	2747	!! dt_a (i,k) = dt_a(i,k) + d_t_vdf_x(i,k)/dtime
	2748	!! dq_a (i,k) = dq_a(i,k) + d_q_vdf_x(i,k)/dtime
	2749	ENDDO
	2750	ENDDO
	2751	ENDIF
	2752	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
	2753	DO k=1,klev
	2754	DO i=1,klon
	2755	!! dt_dwn(i,k) = dt_dwn(i,k) + 0.
	2756	!! dq_dwn(i,k) = dq_dwn(i,k) + 0.
	2757	!! dt_a(i,k) = dt_a(i,k) + d_t_ajs(i,k)/dtime
	2758	!! dq_a(i,k) = dq_a(i,k) + d_q_ajs(i,k)/dtime
	2759	udt_PBL(i,k) = udt_PBL(i,k) + d_t_ajs(i,k)/dtime
	2760	udq_PBL(i,k) = udq_PBL(i,k) + d_q_ajs(i,k)/dtime
	2761	ENDDO
	2762	ENDDO
	2763	ENDIF
	2764	!>nrlmd+jyg
[1507]	2765
[2469]	2766	IF (iflag_wake==2) THEN
	2767	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
	2768	DO k = 1,klev
	2769	dt_dwn(:,k)= dt_dwn(:,k)+ &
	2770	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
	2771	dq_dwn(:,k)= dq_dwn(:,k)+ &
	2772	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
	2773	ENDDO
	2774	ELSEIF (iflag_wake==3) THEN
	2775	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
	2776	DO k = 1,klev
	2777	DO i=1,klon
	2778	IF (rneb(i,k)==0.) THEN
	2779	! On ne tient compte des tendances qu'en dehors des
	2780	! nuages (c'est-\`a-dire a priri dans une region ou
	2781	! l'eau se reevapore).
	2782	dt_dwn(i,k)= dt_dwn(i,k)+ &
	2783	ok_wk_lsp(i)*d_t_lsc(i,k)/dtime
	2784	dq_dwn(i,k)= dq_dwn(i,k)+ &
	2785	ok_wk_lsp(i)*d_q_lsc(i,k)/dtime
	2786	ENDIF
	2787	ENDDO
	2788	ENDDO
	2789	ENDIF
[2105]	2790
[2469]	2791	!
	2792	!calcul caracteristiques de la poche froide
	2793	call calWAKE (paprs,pplay,dtime &
	2794	,t_seri,q_seri,omega &
	2795	,dt_dwn,dq_dwn,M_dwn,M_up &
	2796	,dt_a,dq_a,sigd &
	2797	,wdt_PBL,wdq_PBL &
	2798	,udt_PBL,udq_PBL &
	2799	,wake_deltat,wake_deltaq,wake_dth &
	2800	,wake_h,wake_s,wake_dens &
	2801	,wake_pe,wake_fip,wake_gfl &
	2802	,dt_wake,dq_wake &
	2803	,wake_k, t_undi,q_undi &
	2804	,wake_omgbdth,wake_dp_omgb &
	2805	,wake_dtKE,wake_dqKE &
	2806	,wake_dtPBL,wake_dqPBL &
	2807	,wake_omg,wake_dp_deltomg &
	2808	,wake_spread,wake_Cstar,wake_d_deltat_gw &
	2809	,wake_ddeltat,wake_ddeltaq)
	2810	!
	2811	!-----------------------------------------------------------------------
	2812	! ajout des tendances des poches froides
	2813	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
	2814	! coherent avec les autres d_t_...
	2815	d_t_wake(:,:)=dt_wake(:,:)*dtime
	2816	d_q_wake(:,:)=dq_wake(:,:)*dtime
	2817	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
[2611]	2818	abortphy,flag_inhib_tend)
[2469]	2819	!------------------------------------------------------------------------
[879]	2820
[2469]	2821	endif ! (iflag_wake>=1)
	2822	!
	2823	!===================================================================
	2824	!JYG
	2825	IF (ip_ebil_phy.ge.2) THEN
	2826	ztit='after wake'
	2827	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
	2828	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	2829	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	2830	call diagphy(cell_area,ztit,ip_ebil_phy &
	2831	, zero_v, zero_v, zero_v, zero_v, zero_v &
	2832	, zero_v, zero_v, zero_v, ztsol &
	2833	, d_h_vcol, d_qt, d_ec &
	2834	, fs_bound, fq_bound )
	2835	END IF
[1403]	2836
[2469]	2837	! print*,'apres callwake iflag_cld_th=', iflag_cld_th
	2838	!
	2839	!===================================================================
	2840	! Convection seche (thermiques ou ajustement)
	2841	!===================================================================
	2842	!
	2843	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
	2844	,seuil_inversion,weak_inversion,dthmin)
[878]	2845
	2846
	2847
[2469]	2848	d_t_ajsb(:,:)=0.
	2849	d_q_ajsb(:,:)=0.
	2850	d_t_ajs(:,:)=0.
	2851	d_u_ajs(:,:)=0.
	2852	d_v_ajs(:,:)=0.
	2853	d_q_ajs(:,:)=0.
	2854	clwcon0th(:,:)=0.
	2855	!
	2856	! fm_therm(:,:)=0.
	2857	! entr_therm(:,:)=0.
	2858	! detr_therm(:,:)=0.
	2859	!
	2860	IF(prt_level>9)WRITE(lunout,*) &
	2861	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
	2862	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
	2863	if(iflag_thermals<0) then
	2864	! Rien
	2865	! ====
	2866	IF(prt_level>9)WRITE(lunout,*)'pas de convection seche'
[541]	2867
[878]	2868
[2469]	2869	else
[878]	2870
[2469]	2871	! Thermiques
	2872	! ==========
	2873	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
	2874	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
[878]	2875
	2876
[2469]	2877	!cc nrlmd le 10/04/2012
	2878	DO k=1,klev+1
	2879	DO i=1,klon
	2880	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
	2881	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
	2882	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
	2883	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
[2159]	2884	ENDDO
[2469]	2885	ENDDO
	2886	!cc fin nrlmd le 10/04/2012
[1403]	2887
[2469]	2888	if (iflag_thermals>=1) then
	2889	!jyg<
	2890	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
	2891	! Appel des thermiques avec les profils exterieurs aux poches
	2892	DO k=1,klev
	2893	DO i=1,klon
	2894	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
	2895	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
[2606]	2896	u_therm(i,k) = u_seri(i,k)
	2897	v_therm(i,k) = v_seri(i,k)
[2469]	2898	ENDDO
	2899	ENDDO
	2900	ELSE
	2901	! Appel des thermiques avec les profils moyens
	2902	DO k=1,klev
	2903	DO i=1,klon
	2904	t_therm(i,k) = t_seri(i,k)
	2905	q_therm(i,k) = q_seri(i,k)
[2606]	2906	u_therm(i,k) = u_seri(i,k)
	2907	v_therm(i,k) = v_seri(i,k)
[2469]	2908	ENDDO
	2909	ENDDO
	2910	ENDIF
	2911	!>jyg
	2912	call calltherm(pdtphys &
	2913	,pplay,paprs,pphi,weak_inversion &
[2606]	2914	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
	2915	,u_therm,v_therm,t_therm,q_therm,zqsat,debut & !jyg
[2469]	2916	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
	2917	,fm_therm,entr_therm,detr_therm &
	2918	,zqasc,clwcon0th,lmax_th,ratqscth &
	2919	,ratqsdiff,zqsatth &
	2920	!on rajoute ale et alp, et les
	2921	!caracteristiques de la couche alim
	2922	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
	2923	,ztv,zpspsk,ztla,zthl &
	2924	!cc nrlmd le 10/04/2012
	2925	,pbl_tke_input,pctsrf,omega,cell_area &
	2926	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
	2927	,n2,s2,ale_bl_stat &
	2928	,therm_tke_max,env_tke_max &
	2929	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
	2930	,alp_bl_conv,alp_bl_stat &
	2931	!cc fin nrlmd le 10/04/2012
	2932	,zqla,ztva )
	2933	!
	2934	!jyg<
	2935	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
	2936	! Si les thermiques ne sont presents que hors des
	2937	! poches, la tendance moyenne associ\'ee doit etre
	2938	! multipliee par la fraction surfacique qu'ils couvrent.
	2939	DO k=1,klev
	2940	DO i=1,klon
	2941	!
	2942	wake_deltat(i,k) = wake_deltat(i,k) - d_t_ajs(i,k)
	2943	wake_deltaq(i,k) = wake_deltaq(i,k) - d_q_ajs(i,k)
	2944	!
[2606]	2945	!!!t_seri(i,k) = t_therm(i,k) + wake_s(i)*wake_deltat(i,k)
	2946	!!!q_seri(i,k) = q_therm(i,k) + wake_s(i)*wake_deltaq(i,k)
	2947	!
[2469]	2948	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
	2949	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
	2950	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
	2951	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
	2952	!
	2953	ENDDO
	2954	ENDDO
[2606]	2955	!!! ELSE
	2956	!!! DO k=1,klev
	2957	!!! DO i=1,klon
	2958	!!! t_seri(i,k) = t_therm(i,k)
	2959	!!! q_seri(i,k) = q_therm(i,k)
	2960	!!! ENDDO
	2961	!!! ENDDO
[2469]	2962	ENDIF
[2606]	2963	!
	2964	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
[2611]	2965	dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend)
[2606]	2966	!
[2469]	2967	!>jyg
[2513]	2968	!jyg<
	2969	!
[2565]	2970	CALL alpale_th( dtime, lmax_th, t_seri, cell_area, &
[2513]	2971	cin, s2, n2, &
	2972	ale_bl_trig, ale_bl_stat, ale_bl, &
[2556]	2973	alp_bl, alp_bl_stat, &
	2974	proba_notrig, random_notrig)
[1638]	2975
[2554]	2976	! ------------------------------------------------------------------
	2977	! Transport de la TKE par les panaches thermiques.
	2978	! FH : 2010/02/01
	2979	! if (iflag_pbl.eq.10) then
	2980	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
	2981	! s rg,paprs,pbl_tke)
	2982	! endif
	2983	! -------------------------------------------------------------------
	2984
[2469]	2985	do i=1,klon
	2986	! zmax_th(i)=pphi(i,lmax_th(i))/rg
	2987	!CR:04/05/12:correction calcul zmax
	2988	zmax_th(i)=zmax0(i)
	2989	enddo
[1507]	2990
[2469]	2991	endif
[878]	2992
	2993
[2469]	2994	! Ajustement sec
	2995	! ==============
[878]	2996
[2469]	2997	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
	2998	! a partir du sommet des thermiques.
	2999	! Dans le cas contraire, on demarre au niveau 1.
[878]	3000
[2469]	3001	if (iflag_thermals>=13.or.iflag_thermals<=0) then
[878]	3002
[2469]	3003	if(iflag_thermals.eq.0) then
	3004	IF(prt_level>9)WRITE(lunout,*)'ajsec'
	3005	limbas(:)=1
	3006	else
	3007	limbas(:)=lmax_th(:)
	3008	endif
[878]	3009
[2469]	3010	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
	3011	! pour des test de convergence numerique.
	3012	! Le nouveau ajsec est a priori mieux, meme pour le cas
	3013	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
	3014	! non nulles numeriquement pour des mailles non concernees.
[878]	3015
[2469]	3016	if (iflag_thermals==0) then
	3017	! Calling adjustment alone (but not the thermal plume model)
	3018	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
	3019	, d_t_ajsb, d_q_ajsb)
	3020	else if (iflag_thermals>0) then
	3021	! Calling adjustment above the top of thermal plumes
	3022	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
	3023	, d_t_ajsb, d_q_ajsb)
	3024	endif
[878]	3025
[2469]	3026	!--------------------------------------------------------------------
	3027	! ajout des tendances de l'ajustement sec ou des thermiques
	3028	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
[2611]	3029	'ajsb',abortphy,flag_inhib_tend)
[2469]	3030	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
	3031	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
[904]	3032
[2469]	3033	!---------------------------------------------------------------------
[878]	3034
[2469]	3035	endif
[524]	3036
[2469]	3037	endif
	3038	!
	3039	!===================================================================
	3040	!IM
	3041	IF (ip_ebil_phy.ge.2) THEN
	3042	ztit='after dry_adjust'
	3043	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
	3044	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3045	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3046	call diagphy(cell_area,ztit,ip_ebil_phy &
	3047	, zero_v, zero_v, zero_v, zero_v, zero_v &
	3048	, zero_v, zero_v, zero_v, ztsol &
	3049	, d_h_vcol, d_qt, d_ec &
	3050	, fs_bound, fq_bound )
	3051	END IF
[524]	3052
	3053
[2469]	3054	!-------------------------------------------------------------------------
	3055	! Computation of ratqs, the width (normalized) of the subrid scale
	3056	! water distribution
	3057	CALL calcratqs(klon,klev,prt_level,lunout, &
	3058	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
[2534]	3059	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
	3060	tau_ratqs,fact_cldcon, &
[2469]	3061	ptconv,ptconvth,clwcon0th, rnebcon0th, &
	3062	paprs,pplay,q_seri,zqsat,fm_therm, &
	3063	ratqs,ratqsc)
[1032]	3064
[2100]	3065
[2469]	3066	!
	3067	! Appeler le processus de condensation a grande echelle
	3068	! et le processus de precipitation
	3069	!-------------------------------------------------------------------------
	3070	IF (prt_level .GE.10) THEN
	3071	print *,'itap, ->fisrtilp ',itap
	3072	ENDIF
	3073	!
	3074	CALL fisrtilp(dtime,paprs,pplay, &
	3075	t_seri, q_seri,ptconv,ratqs, &
	3076	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
	3077	rain_lsc, snow_lsc, &
	3078	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
	3079	frac_impa, frac_nucl, beta_prec_fisrt, &
	3080	prfl, psfl, rhcl, &
	3081	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
	3082	iflag_ice_thermo)
	3083	!
	3084	WHERE (rain_lsc < 0) rain_lsc = 0.
	3085	WHERE (snow_lsc < 0) snow_lsc = 0.
[766]	3086
[2469]	3087	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
[2611]	3088	'lsc',abortphy,flag_inhib_tend)
[2613]	3089	rain_num(:)=0.
	3090	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
	3091	DO i = 1, klon
	3092	IF (ql_seri(i,k)>oliqmax) THEN
	3093	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
	3094	ql_seri(i,k)=oliqmax
	3095	ENDIF
	3096	ENDDO
	3097	ENDDO
	3098
[2524]	3099	!---------------------------------------------------------------------------
[2469]	3100	DO k = 1, klev
	3101	DO i = 1, klon
	3102	cldfra(i,k) = rneb(i,k)
	3103	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
	3104	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
	3105	ENDDO
	3106	ENDDO
	3107	IF (check) THEN
	3108	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
	3109	WRITE(lunout,*)"apresilp=", za
	3110	zx_t = 0.0
	3111	za = 0.0
	3112	DO i = 1, klon
	3113	za = za + cell_area(i)/REAL(klon)
	3114	zx_t = zx_t + (rain_lsc(i) &
	3115	+ snow_lsc(i))*cell_area(i)/REAL(klon)
	3116	ENDDO
	3117	zx_t = zx_t/za*dtime
	3118	WRITE(lunout,*)"Precip=", zx_t
	3119	ENDIF
	3120	!IM
	3121	IF (ip_ebil_phy.ge.2) THEN
	3122	ztit='after fisrt'
	3123	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
	3124	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3125	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3126	call diagphy(cell_area,ztit,ip_ebil_phy &
	3127	, zero_v, zero_v, zero_v, zero_v, zero_v &
	3128	, zero_v, rain_lsc, snow_lsc, ztsol &
	3129	, d_h_vcol, d_qt, d_ec &
	3130	, fs_bound, fq_bound )
	3131	END IF
[766]	3132
[2469]	3133	if (mydebug) then
	3134	call writefield_phy('u_seri',u_seri,nbp_lev)
	3135	call writefield_phy('v_seri',v_seri,nbp_lev)
	3136	call writefield_phy('t_seri',t_seri,nbp_lev)
	3137	call writefield_phy('q_seri',q_seri,nbp_lev)
	3138	endif
[524]	3139
[2469]	3140	!
	3141	!-------------------------------------------------------------------
	3142	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
	3143	!-------------------------------------------------------------------
[524]	3144
[2469]	3145	! 1. NUAGES CONVECTIFS
	3146	!
	3147	!IM cf FH
	3148	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
	3149	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
	3150	snow_tiedtke=0.
	3151	! print*,'avant calcul de la pseudo precip '
	3152	! print*,'iflag_cld_th',iflag_cld_th
	3153	if (iflag_cld_th.eq.-1) then
	3154	rain_tiedtke=rain_con
	3155	else
	3156	! print*,'calcul de la pseudo precip '
	3157	rain_tiedtke=0.
	3158	! print*,'calcul de la pseudo precip 0'
	3159	do k=1,klev
	3160	do i=1,klon
	3161	if (d_q_con(i,k).lt.0.) then
	3162	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
	3163	*(paprs(i,k)-paprs(i,k+1))/rg
	3164	endif
	3165	enddo
	3166	enddo
	3167	endif
	3168	!
	3169	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
	3170	!
[524]	3171
[2469]	3172	! Nuages diagnostiques pour Tiedtke
	3173	CALL diagcld1(paprs,pplay, &
	3174	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
	3175	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
	3176	diafra,dialiq)
	3177	DO k = 1, klev
	3178	DO i = 1, klon
	3179	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	3180	cldliq(i,k) = dialiq(i,k)
	3181	cldfra(i,k) = diafra(i,k)
	3182	ENDIF
	3183	ENDDO
	3184	ENDDO
[524]	3185
[2469]	3186	ELSE IF (iflag_cld_th.ge.3) THEN
	3187	! On prend pour les nuages convectifs le max du calcul de la
	3188	! convection et du calcul du pas de temps precedent diminue d'un facteur
	3189	! facttemps
	3190	facteur = pdtphys *facttemps
	3191	do k=1,klev
	3192	do i=1,klon
	3193	rnebcon(i,k)=rnebcon(i,k)*facteur
	3194	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k)) &
	3195	then
	3196	rnebcon(i,k)=rnebcon0(i,k)
	3197	clwcon(i,k)=clwcon0(i,k)
	3198	endif
	3199	enddo
	3200	enddo
	3201
	3202	! On prend la somme des fractions nuageuses et des contenus en eau
[524]	3203
[2469]	3204	if (iflag_cld_th>=5) then
[1411]	3205
[2469]	3206	do k=1,klev
	3207	ptconvth(:,k)=fm_therm(:,k+1)>0.
	3208	enddo
[1496]	3209
[2469]	3210	if (iflag_coupl==4) then
[1496]	3211
[2469]	3212	! Dans le cas iflag_coupl==4, on prend la somme des convertures
	3213	! convectives et lsc dans la partie des thermiques
	3214	! Le controle par iflag_coupl est peut etre provisoire.
	3215	do k=1,klev
	3216	do i=1,klon
	3217	if (ptconv(i,k).and.ptconvth(i,k)) then
	3218	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
	3219	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
	3220	else if (ptconv(i,k)) then
	3221	cldfra(i,k)=rnebcon(i,k)
	3222	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
	3223	endif
	3224	enddo
	3225	enddo
[1496]	3226
[2469]	3227	else if (iflag_coupl==5) then
	3228	do k=1,klev
	3229	do i=1,klon
	3230	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
	3231	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
	3232	enddo
	3233	enddo
[1525]	3234
[2469]	3235	else
[1525]	3236
[2469]	3237	! Si on est sur un point touche par la convection
	3238	! profonde et pas par les thermiques, on prend la
	3239	! couverture nuageuse et l'eau nuageuse de la convection
	3240	! profonde.
[1411]	3241
[2469]	3242	!IM/FH: 2011/02/23
	3243	! definition des points sur lesquels ls thermiques sont actifs
[1496]	3244
[2469]	3245	do k=1,klev
	3246	do i=1,klon
	3247	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
	3248	cldfra(i,k)=rnebcon(i,k)
	3249	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
	3250	endif
	3251	enddo
	3252	enddo
[1496]	3253
[2469]	3254	endif
[1496]	3255
[2469]	3256	else
[1496]	3257
[2469]	3258	! Ancienne version
	3259	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
	3260	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
	3261	endif
[1411]	3262
[2469]	3263	ENDIF
[1507]	3264
[2469]	3265	! plulsc(:)=0.
	3266	! do k=1,klev,-1
	3267	! do i=1,klon
	3268	! zzz=prfl(:,k)+psfl(:,k)
	3269	! if (.not.ptconvth.zzz.gt.0.)
	3270	! enddo prfl, psfl,
	3271	! enddo
	3272	!
	3273	! 2. NUAGES STARTIFORMES
	3274	!
	3275	IF (ok_stratus) THEN
	3276	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
	3277	DO k = 1, klev
	3278	DO i = 1, klon
	3279	IF (diafra(i,k).GT.cldfra(i,k)) THEN
	3280	cldliq(i,k) = dialiq(i,k)
	3281	cldfra(i,k) = diafra(i,k)
	3282	ENDIF
	3283	ENDDO
	3284	ENDDO
	3285	ENDIF
	3286	!
	3287	! Precipitation totale
	3288	!
	3289	DO i = 1, klon
	3290	rain_fall(i) = rain_con(i) + rain_lsc(i)
	3291	snow_fall(i) = snow_con(i) + snow_lsc(i)
	3292	ENDDO
	3293	!IM
	3294	IF (ip_ebil_phy.ge.2) THEN
	3295	ztit="after diagcld"
	3296	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
	3297	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3298	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3299	call diagphy(cell_area,ztit,ip_ebil_phy &
	3300	, zero_v, zero_v, zero_v, zero_v, zero_v &
	3301	, zero_v, zero_v, zero_v, ztsol &
	3302	, d_h_vcol, d_qt, d_ec &
	3303	, fs_bound, fq_bound )
	3304	END IF
	3305	!
	3306	! Calculer l'humidite relative pour diagnostique
	3307	!
	3308	DO k = 1, klev
	3309	DO i = 1, klon
	3310	zx_t = t_seri(i,k)
	3311	IF (thermcep) THEN
	3312	!! if (iflag_ice_thermo.eq.0) then !jyg
	3313	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
	3314	!! else !jyg
	3315	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
	3316	!! endif !jyg
	3317	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
	3318	zx_qs = MIN(0.5,zx_qs)
	3319	zcor = 1./(1.-retv*zx_qs)
	3320	zx_qs = zx_qs*zcor
	3321	ELSE
	3322	!! IF (zx_t.LT.t_coup) THEN !jyg
	3323	IF (zx_t.LT.rtt) THEN !jyg
	3324	zx_qs = qsats(zx_t)/pplay(i,k)
	3325	ELSE
	3326	zx_qs = qsatl(zx_t)/pplay(i,k)
	3327	ENDIF
	3328	ENDIF
	3329	zx_rh(i,k) = q_seri(i,k)/zx_qs
	3330	zqsat(i,k)=zx_qs
	3331	ENDDO
	3332	ENDDO
[782]	3333
[2469]	3334	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
	3335	! equivalente a 2m (tpote) pour diagnostique
	3336	!
	3337	DO i = 1, klon
	3338	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
	3339	IF (thermcep) THEN
	3340	IF(zt2m(i).LT.RTT) then
	3341	Lheat=RLSTT
	3342	ELSE
	3343	Lheat=RLVTT
	3344	ENDIF
	3345	ELSE
	3346	IF (zt2m(i).LT.RTT) THEN
	3347	Lheat=RLSTT
	3348	ELSE
	3349	Lheat=RLVTT
	3350	ENDIF
	3351	ENDIF
	3352	tpote(i) = tpot(i)* &
	3353	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
	3354	ENDDO
[524]	3355
[2469]	3356	IF (type_trac == 'inca') THEN
[524]	3357	#ifdef INCA
[2469]	3358	CALL VTe(VTphysiq)
	3359	CALL VTb(VTinca)
	3360	calday = REAL(days_elapsed + 1) + jH_cur
[524]	3361
[2469]	3362	call chemtime(itap+itau_phy-1, date0, dtime, itap)
[3606]	3363	CALL AEROSOL_METEO_CALC( &
	3364	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
	3365	prfl,psfl,pctsrf,cell_area, &
	3366	latitude_deg,longitude_deg,u10m,v10m)
[524]	3367
[2469]	3368	zxsnow_dummy(:) = 0.0
[625]	3369
[2469]	3370	CALL chemhook_begin (calday, &
	3371	days_elapsed+1, &
	3372	jH_cur, &
	3373	pctsrf(1,1), &
	3374	latitude_deg, &
	3375	longitude_deg, &
	3376	cell_area, &
	3377	paprs, &
	3378	pplay, &
	3379	coefh(1:klon,1:klev,is_ave), &
	3380	pphi, &
	3381	t_seri, &
	3382	u, &
	3383	v, &
[3611]	3384	rot, &
[2469]	3385	wo(:, :, 1), &
	3386	q_seri, &
	3387	zxtsol, &
[3611]	3388	zt2m, &
[2469]	3389	zxsnow_dummy, &
	3390	solsw, &
	3391	albsol1, &
	3392	rain_fall, &
	3393	snow_fall, &
	3394	itop_con, &
	3395	ibas_con, &
	3396	cldfra, &
	3397	nbp_lon, &
	3398	nbp_lat-1, &
	3399	tr_seri, &
	3400	ftsol, &
	3401	paprs, &
	3402	cdragh, &
	3403	cdragm, &
	3404	pctsrf, &
	3405	pdtphys, &
	3406	itap)
[616]	3407
[2469]	3408	CALL VTe(VTinca)
	3409	CALL VTb(VTphysiq)
[959]	3410	#endif
[2469]	3411	END IF !type_trac = inca
[2618]	3412
	3413
[2469]	3414	!
[2618]	3415	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
	3416	!
	3417	IF (MOD(itaprad,radpas).EQ.0) THEN
[959]	3418
[2618]	3419	!
	3420	!jq - introduce the aerosol direct and first indirect radiative forcings
	3421	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
	3422	IF (flag_aerosol .gt. 0) THEN
	3423	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
	3424	IF (.NOT. aerosol_couple) THEN
	3425	!
	3426	CALL readaerosol_optic( &
	3427	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
	3428	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3429	mass_solu_aero, mass_solu_aero_pi, &
	3430	tau_aero, piz_aero, cg_aero, &
	3431	tausum_aero, tau3d_aero)
	3432	ENDIF
	3433	ELSE ! RRTM radiation
	3434	!
	3435	#ifdef CPP_RRTM
	3436	IF (NSW.EQ.6) THEN
	3437	!--new aerosol properties
	3438	!
	3439	CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
	3440	new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
	3441	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3442	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
	3443	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
	3444	tausum_aero, tau3d_aero)
[955]	3445
[2618]	3446	ELSE IF (NSW.EQ.2) THEN
	3447	!--for now we use the old aerosol properties
	3448	!
	3449	CALL readaerosol_optic( &
	3450	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
	3451	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3452	mass_solu_aero, mass_solu_aero_pi, &
	3453	tau_aero, piz_aero, cg_aero, &
	3454	tausum_aero, tau3d_aero)
	3455	!
	3456	!--natural aerosols
	3457	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
	3458	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
	3459	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
	3460	!--all aerosols
	3461	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
	3462	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
	3463	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
	3464	ELSE
	3465	abort_message='Only NSW=2 or 6 are possible with ' &
	3466	// 'aerosols and iflag_rrtm=1'
	3467	call abort_physic(modname,abort_message,1)
	3468	ENDIF
	3469
	3470	!--call LW optical properties for tropospheric aerosols
	3471	!--only works for INCA aerosol (aerosol_couple = TRUE)
	3472	CALL aeropt_lw_rrtm(aerosol_couple,paprs,tr_seri)
	3473	!
	3474	#else
	3475	abort_message='You should compile with -rrtm if running ' &
	3476	// 'with iflag_rrtm=1'
	3477	call abort_physic(modname,abort_message,1)
	3478	#endif
	3479	!
	3480	ENDIF
	3481	ELSE
	3482	tausum_aero(:,:,:) = 0.
	3483	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
	3484	tau_aero(:,:,:,:) = 1.e-15
	3485	piz_aero(:,:,:,:) = 1.
	3486	cg_aero(:,:,:,:) = 0.
	3487	ELSE
	3488	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
	3489	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
	3490	piz_aero_sw_rrtm(:,:,:,:) = 1.0
	3491	cg_aero_sw_rrtm(:,:,:,:) = 0.0
	3492	ENDIF
	3493	ENDIF
	3494	!
	3495	!--STRAT AEROSOL
	3496	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
	3497	IF (flag_aerosol_strat.GT.0) THEN
	3498	IF (prt_level .GE.10) THEN
	3499	PRINT *,'appel a readaerosolstrat', mth_cur
	3500	ENDIF
	3501	IF (iflag_rrtm.EQ.0) THEN
	3502	IF (flag_aerosol_strat.EQ.1) THEN
	3503	CALL readaerosolstrato(debut)
	3504	ELSE
	3505	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
	3506	CALL abort_physic(modname,abort_message,1)
	3507	ENDIF
	3508	ELSE
[2009]	3509	#ifdef CPP_RRTM
[2618]	3510	IF (flag_aerosol_strat.EQ.1) THEN
	3511	CALL readaerosolstrato1_rrtm(debut)
	3512	ELSEIF (flag_aerosol_strat.EQ.2) THEN
	3513	CALL stratosphere_mask(t_seri, pplay, latitude_deg)
	3514	CALL readaerosolstrato2_rrtm(debut)
	3515	ELSE
	3516	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
	3517	CALL abort_physic(modname,abort_message,1)
	3518	ENDIF
	3519	#else
	3520	abort_message='You should compile with -rrtm if running ' &
	3521	// 'with iflag_rrtm=1'
	3522	CALL abort_physic(modname,abort_message,1)
	3523	#endif
	3524	ENDIF
	3525	ENDIF
	3526	!--fin STRAT AEROSOL
	3527	!
	3528
	3529	! Calculer les parametres optiques des nuages et quelques
	3530	! parametres pour diagnostiques:
	3531	!
	3532	IF (aerosol_couple.AND.config_inca=='aero') THEN
	3533	mass_solu_aero(:,:) = ccm(:,:,1)
	3534	mass_solu_aero_pi(:,:) = ccm(:,:,2)
	3535	END IF
	3536
	3537	IF (ok_newmicro) then
	3538	IF (iflag_rrtm.NE.0) THEN
	3539	#ifdef CPP_RRTM
	3540	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
[2469]	3541	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
	3542	// 'pour ok_cdnc'
[2618]	3543	CALL abort_physic(modname,abort_message,1)
	3544	ENDIF
[2009]	3545	#else
	3546
[2618]	3547	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
	3548	CALL abort_physic(modname,abort_message,1)
[2009]	3549	#endif
[2618]	3550	ENDIF
	3551	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
	3552	paprs, pplay, t_seri, cldliq, cldfra, &
	3553	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
	3554	flwp, fiwp, flwc, fiwc, &
	3555	mass_solu_aero, mass_solu_aero_pi, &
	3556	cldtaupi, re, fl, ref_liq, ref_ice, &
	3557	ref_liq_pi, ref_ice_pi)
	3558	ELSE
	3559	CALL nuage (paprs, pplay, &
	3560	t_seri, cldliq, cldfra, cldtau, cldemi, &
	3561	cldh, cldl, cldm, cldt, cldq, &
	3562	ok_aie, &
	3563	mass_solu_aero, mass_solu_aero_pi, &
	3564	bl95_b0, bl95_b1, &
	3565	cldtaupi, re, fl)
[2469]	3566	ENDIF
	3567	!
[2618]	3568	!IM betaCRF
[2469]	3569	!
[2618]	3570	cldtaurad = cldtau
	3571	cldtaupirad = cldtaupi
	3572	cldemirad = cldemi
	3573	cldfrarad = cldfra
	3574
[2469]	3575	!
[2618]	3576	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
	3577	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
	3578	!
	3579	! global
	3580	!
	3581	DO k=1, klev
	3582	DO i=1, klon
	3583	IF (pplay(i,k).GE.pfree) THEN
[2469]	3584	beta(i,k) = beta_pbl
[2618]	3585	ELSE
[2469]	3586	beta(i,k) = beta_free
[2618]	3587	ENDIF
	3588	IF (mskocean_beta) THEN
[2469]	3589	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
[2618]	3590	ENDIF
[2469]	3591	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
	3592	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
	3593	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
	3594	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
[2618]	3595	ENDDO
	3596	ENDDO
	3597	!
	3598	ELSE
	3599	!
	3600	! regional
	3601	!
	3602	DO k=1, klev
	3603	DO i=1,klon
	3604	!
	3605	IF (longitude_deg(i).ge.lon1_beta.AND. &
	3606	longitude_deg(i).le.lon2_beta.AND. &
	3607	latitude_deg(i).le.lat1_beta.AND. &
	3608	latitude_deg(i).ge.lat2_beta) THEN
	3609	IF (pplay(i,k).GE.pfree) THEN
	3610	beta(i,k) = beta_pbl
	3611	ELSE
	3612	beta(i,k) = beta_free
	3613	ENDIF
	3614	IF (mskocean_beta) THEN
	3615	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
	3616	ENDIF
	3617	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
	3618	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
	3619	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
	3620	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
	3621	ENDIF
[2469]	3622	!
[2618]	3623	ENDDO
[2469]	3624	ENDDO
	3625	!
[2618]	3626	ENDIF
[766]	3627
[2618]	3628	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
	3629	IF (ok_chlorophyll) THEN
[2469]	3630	print*,"-- reading chlorophyll"
[2618]	3631	CALL readchlorophyll(debut)
	3632	ENDIF
[1863]	3633
[2614]	3634	!
	3635	!jq - introduce the aerosol direct and first indirect radiative forcings
	3636	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
	3637	IF (flag_aerosol .gt. 0) THEN
	3638	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
	3639	IF (.NOT. aerosol_couple) THEN
	3640	!
	3641	CALL readaerosol_optic( &
	3642	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
	3643	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3644	mass_solu_aero, mass_solu_aero_pi, &
	3645	tau_aero, piz_aero, cg_aero, &
	3646	tausum_aero, tau3d_aero)
	3647	ENDIF
	3648	ELSE ! RRTM radiation
	3649	!
	3650	#ifdef CPP_RRTM
	3651	IF (NSW.EQ.6) THEN
	3652	!--new aerosol properties
	3653	!
	3654	CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
	3655	new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
	3656	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3657	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
	3658	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
	3659	tausum_aero, tau3d_aero)
	3660
	3661	ELSE IF (NSW.EQ.2) THEN
	3662	!--for now we use the old aerosol properties
	3663	!
	3664	CALL readaerosol_optic( &
	3665	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
	3666	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
	3667	mass_solu_aero, mass_solu_aero_pi, &
	3668	tau_aero, piz_aero, cg_aero, &
	3669	tausum_aero, tau3d_aero)
	3670	!
	3671	!--natural aerosols
	3672	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
	3673	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
	3674	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
	3675	!--all aerosols
	3676	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
	3677	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
	3678	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
	3679	ELSE
	3680	abort_message='Only NSW=2 or 6 are possible with ' &
	3681	// 'aerosols and iflag_rrtm=1'
	3682	call abort_physic(modname,abort_message,1)
	3683	ENDIF
	3684
	3685	!--call LW optical properties for tropospheric aerosols
	3686	!--only works for INCA aerosol (aerosol_couple = TRUE)
	3687	CALL aeropt_lw_rrtm(aerosol_couple,paprs,tr_seri)
	3688	!
	3689	#else
	3690	abort_message='You should compile with -rrtm if running ' &
	3691	// 'with iflag_rrtm=1'
	3692	call abort_physic(modname,abort_message,1)
	3693	#endif
	3694	!
	3695	ENDIF
	3696	ELSE
	3697	tausum_aero(:,:,:) = 0.
	3698	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
	3699	tau_aero(:,:,:,:) = 1.e-15
	3700	piz_aero(:,:,:,:) = 1.
	3701	cg_aero(:,:,:,:) = 0.
	3702	ELSE
	3703	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
	3704	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
	3705	piz_aero_sw_rrtm(:,:,:,:) = 1.0
	3706	cg_aero_sw_rrtm(:,:,:,:) = 0.0
	3707	ENDIF
	3708	ENDIF
	3709	!
	3710	!--STRAT AEROSOL
	3711	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
	3712	IF (flag_aerosol_strat.GT.0) THEN
	3713	IF (prt_level .GE.10) THEN
	3714	PRINT *,'appel a readaerosolstrat', mth_cur
	3715	ENDIF
	3716	IF (iflag_rrtm.EQ.0) THEN
	3717	IF (flag_aerosol_strat.EQ.1) THEN
	3718	CALL readaerosolstrato(debut)
	3719	ELSE
	3720	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
	3721	call abort_physic(modname,abort_message,1)
	3722	ENDIF
	3723	ELSE
	3724	#ifdef CPP_RRTM
	3725	IF (flag_aerosol_strat.EQ.1) THEN
	3726	CALL readaerosolstrato1_rrtm(debut)
	3727	ELSEIF (flag_aerosol_strat.EQ.2) THEN
	3728	CALL stratosphere_mask(t_seri, pplay, latitude_deg)
	3729	CALL readaerosolstrato2_rrtm(debut)
	3730	ELSE
	3731	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
	3732	call abort_physic(modname,abort_message,1)
	3733	ENDIF
	3734	#else
	3735	abort_message='You should compile with -rrtm if running ' &
	3736	// 'with iflag_rrtm=1'
	3737	call abort_physic(modname,abort_message,1)
	3738	#endif
	3739	ENDIF
	3740	ENDIF
	3741	!--fin STRAT AEROSOL
	3742
[2524]	3743	!--if ok_suntime_rrtm we use ancillay data for RSUN
	3744	!--previous values are therefore overwritten
	3745	!--this is needed for CMIP6 runs
	3746	!--and only possible for new radiation scheme
	3747	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
[2525]	3748	#ifdef CPP_RRTM
[2524]	3749	CALL read_rsun_rrtm(debut)
[2525]	3750	#endif
[3576]	3751	ELSE IF (ok_suntime_rrtm) THEN
	3752	CALL read_rsun_volcanoes(debut)
	3753	ENDIF
[2524]	3754
	3755	if (mydebug) then
	3756	call writefield_phy('u_seri',u_seri,nbp_lev)
	3757	call writefield_phy('v_seri',v_seri,nbp_lev)
	3758	call writefield_phy('t_seri',t_seri,nbp_lev)
	3759	call writefield_phy('q_seri',q_seri,nbp_lev)
	3760	endif
	3761
[2469]	3762	!
	3763	!sonia : If Iflag_radia >=2, pertubation of some variables
	3764	!input to radiation (DICE)
	3765	!
	3766	IF (iflag_radia .ge. 2) THEN
	3767	zsav_tsol (:) = zxtsol(:)
	3768	call perturb_radlwsw(zxtsol,iflag_radia)
	3769	ENDIF
[2328]	3770
[3514]	3771	IF (aerosol_couple.AND.iflag_rrtm==0) THEN
[959]	3772	#ifdef INCA
[2469]	3773	CALL radlwsw_inca &
[3340]	3774	(chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
[2469]	3775	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
[2907]	3776	size(wo,3), wo(:, :, 1), &
[2469]	3777	cldfrarad, cldemirad, cldtaurad, &
	3778	heat,heat0,cool,cool0,albpla, &
	3779	topsw,toplw,solsw,sollw, &
	3780	sollwdown, &
	3781	topsw0,toplw0,solsw0,sollw0, &
	3782	lwdn0, lwdn, lwup0, lwup, &
	3783	swdn0, swdn, swup0, swup, &
	3784	ok_ade, ok_aie, &
	3785	tau_aero, piz_aero, cg_aero, &
	3786	topswad_aero, solswad_aero, &
	3787	topswad0_aero, solswad0_aero, &
	3788	topsw_aero, topsw0_aero, &
	3789	solsw_aero, solsw0_aero, &
	3790	cldtaupirad, &
	3791	topswai_aero, solswai_aero)
[955]	3792	#endif
[2469]	3793	ELSE
	3794	!
	3795	!IM calcul radiatif pour le cas actuel
	3796	!
	3797	RCO2 = RCO2_act
	3798	RCH4 = RCH4_act
	3799	RN2O = RN2O_act
	3800	RCFC11 = RCFC11_act
	3801	RCFC12 = RCFC12_act
	3802	!
	3803	IF (prt_level .GE.10) THEN
	3804	print *,' ->radlwsw, number 1 '
	3805	ENDIF
	3806	!
	3807	CALL radlwsw &
	3808	(dist, rmu0, fract, &
	3809	!albedo SB >>>
	3810	! paprs, pplay,zxtsol,albsol1, albsol2, &
	3811	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
	3812	!albedo SB <<<
	3813	t_seri,q_seri,wo, &
	3814	cldfrarad, cldemirad, cldtaurad, &
[2530]	3815	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
[2469]	3816	flag_aerosol_strat, &
	3817	tau_aero, piz_aero, cg_aero, &
	3818	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
	3819	! Rajoute par OB pour RRTM
	3820	tau_aero_lw_rrtm, &
	3821	cldtaupirad,new_aod, &
	3822	zqsat, flwc, fiwc, &
	3823	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
	3824	heat,heat0,cool,cool0,albpla, &
	3825	topsw,toplw,solsw,sollw, &
	3826	sollwdown, &
	3827	topsw0,toplw0,solsw0,sollw0, &
	3828	lwdn0, lwdn, lwup0, lwup, &
	3829	swdn0, swdn, swup0, swup, &
	3830	topswad_aero, solswad_aero, &
	3831	topswai_aero, solswai_aero, &
	3832	topswad0_aero, solswad0_aero, &
	3833	topsw_aero, topsw0_aero, &
	3834	solsw_aero, solsw0_aero, &
	3835	topswcf_aero, solswcf_aero, &
	3836	!-C. Kleinschmitt for LW diagnostics
	3837	toplwad_aero, sollwad_aero,&
	3838	toplwai_aero, sollwai_aero, &
	3839	toplwad0_aero, sollwad0_aero,&
	3840	!-end
	3841	ZLWFT0_i, ZFLDN0, ZFLUP0, &
	3842	ZSWFT0_i, ZFSDN0, ZFSUP0)
[879]	3843
[2529]	3844	!--OB 30/05/2016
	3845	!--here we return swaero_diag to FALSE
	3846	!--and histdef will switch it back to TRUE if necessary
	3847	!--this is necessary to get the right swaero at first step
	3848	IF (debut) swaero_diag = .FALSE.
[2469]	3849	!
	3850	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
	3851	!IM des taux doit etre different du taux actuel
	3852	!IM Par defaut on a les taux perturbes egaux aux taux actuels
	3853	!
	3854	if (ok_4xCO2atm) then
	3855	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
	3856	RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
	3857	RCFC12_per.NE.RCFC12_act) THEN
	3858	!
	3859	RCO2 = RCO2_per
	3860	RCH4 = RCH4_per
	3861	RN2O = RN2O_per
	3862	RCFC11 = RCFC11_per
	3863	RCFC12 = RCFC12_per
	3864	!
	3865	IF (prt_level .GE.10) THEN
	3866	print *,' ->radlwsw, number 2 '
	3867	ENDIF
	3868	!
	3869	CALL radlwsw &
	3870	(dist, rmu0, fract, &
	3871	!albedo SB >>>
	3872	! paprs, pplay,zxtsol,albsol1, albsol2, &
	3873	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
	3874	!albedo SB <<<
	3875	t_seri,q_seri,wo, &
	3876	cldfra, cldemi, cldtau, &
[2530]	3877	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
[2469]	3878	flag_aerosol_strat, &
	3879	tau_aero, piz_aero, cg_aero, &
	3880	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
	3881	! Rajoute par OB pour RRTM
	3882	tau_aero_lw_rrtm, &
	3883	cldtaupi,new_aod, &
	3884	zqsat, flwc, fiwc, &
	3885	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
	3886	heatp,heat0p,coolp,cool0p,albplap, &
	3887	topswp,toplwp,solswp,sollwp, &
	3888	sollwdownp, &
	3889	topsw0p,toplw0p,solsw0p,sollw0p, &
	3890	lwdn0p, lwdnp, lwup0p, lwupp, &
	3891	swdn0p, swdnp, swup0p, swupp, &
	3892	topswad_aerop, solswad_aerop, &
	3893	topswai_aerop, solswai_aerop, &
	3894	topswad0_aerop, solswad0_aerop, &
	3895	topsw_aerop, topsw0_aerop, &
	3896	solsw_aerop, solsw0_aerop, &
	3897	topswcf_aerop, solswcf_aerop, &
	3898	!-C. Kleinschmitt for LW diagnostics
	3899	toplwad_aerop, sollwad_aerop,&
	3900	toplwai_aerop, sollwai_aerop, &
	3901	toplwad0_aerop, sollwad0_aerop,&
	3902	!-end
	3903	ZLWFT0_i, ZFLDN0, ZFLUP0, &
	3904	ZSWFT0_i, ZFSDN0, ZFSUP0)
	3905	endif
	3906	endif
	3907	!
	3908	ENDIF ! aerosol_couple
	3909	itaprad = 0
	3910	!
	3911	! If Iflag_radia >=2, reset pertubed variables
	3912	!
	3913	IF (iflag_radia .ge. 2) THEN
	3914	zxtsol(:) = zsav_tsol (:)
	3915	ENDIF
	3916	ENDIF ! MOD(itaprad,radpas)
	3917	itaprad = itaprad + 1
[879]	3918
[2469]	3919	IF (iflag_radia.eq.0) THEN
	3920	IF (prt_level.ge.9) THEN
	3921	PRINT *,'--------------------------------------------------'
	3922	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
	3923	PRINT *,'>>>> heat et cool mis a zero '
	3924	PRINT *,'--------------------------------------------------'
	3925	END IF
	3926	heat=0.
	3927	cool=0.
	3928	sollw=0. ! MPL 01032011
	3929	solsw=0.
	3930	radsol=0.
	3931	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
	3932	swup0=0.
	3933	lwup=0.
	3934	lwup0=0.
	3935	lwdn=0.
	3936	lwdn0=0.
	3937	END IF
[782]	3938
[2469]	3939	!
	3940	! Calculer radsol a l'exterieur de radlwsw
	3941	! pour prendre en compte le cycle diurne
	3942	! recode par Olivier Boucher en sept 2015
	3943	!
	3944	radsol=solsw*swradcorr+sollw
[2618]	3945
[2469]	3946	if (ok_4xCO2atm) then
	3947	radsolp=solswp*swradcorr+sollwp
	3948	endif
[2359]	3949
[2469]	3950	!
	3951	! Ajouter la tendance des rayonnements (tous les pas)
	3952	! avec une correction pour le cycle diurne dans le SW
	3953	!
[2359]	3954
[2469]	3955	DO k=1, klev
	3956	d_t_swr(:,k)=swradcorr(:)heat(:,k)dtime/RDAY
	3957	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)dtime/RDAY
	3958	d_t_lwr(:,k)=-cool(:,k)*dtime/RDAY
	3959	d_t_lw0(:,k)=-cool0(:,k)*dtime/RDAY
	3960	ENDDO
[2194]	3961
[2611]	3962	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend)
	3963	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend)
[1863]	3964
[2469]	3965	!
	3966	if (mydebug) then
	3967	call writefield_phy('u_seri',u_seri,nbp_lev)
	3968	call writefield_phy('v_seri',v_seri,nbp_lev)
	3969	call writefield_phy('t_seri',t_seri,nbp_lev)
	3970	call writefield_phy('q_seri',q_seri,nbp_lev)
	3971	endif
[1863]	3972
[2469]	3973	!IM
	3974	IF (ip_ebil_phy.ge.2) THEN
	3975	ztit='after rad'
	3976	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
	3977	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	3978	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	3979	call diagphy(cell_area,ztit,ip_ebil_phy &
	3980	, topsw, toplw, solsw, sollw, zero_v &
	3981	, zero_v, zero_v, zero_v, ztsol &
	3982	, d_h_vcol, d_qt, d_ec &
	3983	, fs_bound, fq_bound )
	3984	END IF
	3985	!
	3986	!
	3987	! Calculer l'hydrologie de la surface
	3988	!
	3989	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
	3990	! . agesno, ftsol,fqsurf,fsnow, ruis)
	3991	!
[1001]	3992
[2469]	3993	!
	3994	! Calculer le bilan du sol et la derive de temperature (couplage)
	3995	!
	3996	DO i = 1, klon
	3997	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
	3998	! a la demande de JLD
	3999	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
	4000	ENDDO
	4001	!
	4002	!moddeblott(jan95)
	4003	! Appeler le programme de parametrisation de l'orographie
	4004	! a l'echelle sous-maille:
	4005	!
	4006	IF (prt_level .GE.10) THEN
	4007	print *,' call orography ? ', ok_orodr
	4008	ENDIF
	4009	!
	4010	IF (ok_orodr) THEN
	4011	!
	4012	! selection des points pour lesquels le shema est actif:
	4013	igwd=0
	4014	DO i=1,klon
	4015	itest(i)=0
	4016	! IF ((zstd(i).gt.10.0)) THEN
	4017	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
	4018	itest(i)=1
	4019	igwd=igwd+1
	4020	idx(igwd)=i
	4021	ENDIF
	4022	ENDDO
	4023	! igwdim=MAX(1,igwd)
	4024	!
	4025	IF (ok_strato) THEN
[1863]	4026
[2469]	4027	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
	4028	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
	4029	igwd,idx,itest, &
	4030	t_seri, u_seri, v_seri, &
	4031	zulow, zvlow, zustrdr, zvstrdr, &
	4032	d_t_oro, d_u_oro, d_v_oro)
[1863]	4033
[2469]	4034	ELSE
	4035	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
	4036	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
	4037	igwd,idx,itest, &
	4038	t_seri, u_seri, v_seri, &
	4039	zulow, zvlow, zustrdr, zvstrdr, &
	4040	d_t_oro, d_u_oro, d_v_oro)
	4041	ENDIF
	4042	!
	4043	! ajout des tendances
	4044	!-----------------------------------------------------------------------
	4045	! ajout des tendances de la trainee de l'orographie
	4046	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
[2611]	4047	abortphy,flag_inhib_tend)
[2469]	4048	!----------------------------------------------------------------------
	4049	!
	4050	ENDIF ! fin de test sur ok_orodr
	4051	!
	4052	if (mydebug) then
	4053	call writefield_phy('u_seri',u_seri,nbp_lev)
	4054	call writefield_phy('v_seri',v_seri,nbp_lev)
	4055	call writefield_phy('t_seri',t_seri,nbp_lev)
	4056	call writefield_phy('q_seri',q_seri,nbp_lev)
	4057	endif
[1001]	4058
[2469]	4059	IF (ok_orolf) THEN
	4060	!
	4061	! selection des points pour lesquels le shema est actif:
	4062	igwd=0
	4063	DO i=1,klon
	4064	itest(i)=0
	4065	IF ((zpic(i)-zmea(i)).GT.100.) THEN
	4066	itest(i)=1
	4067	igwd=igwd+1
	4068	idx(igwd)=i
	4069	ENDIF
	4070	ENDDO
	4071	! igwdim=MAX(1,igwd)
	4072	!
	4073	IF (ok_strato) THEN
[1001]	4074
[2469]	4075	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
	4076	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
	4077	igwd,idx,itest, &
	4078	t_seri, u_seri, v_seri, &
	4079	zulow, zvlow, zustrli, zvstrli, &
	4080	d_t_lif, d_u_lif, d_v_lif )
[2333]	4081
[2469]	4082	ELSE
	4083	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
	4084	latitude_deg,zmea,zstd,zpic, &
	4085	itest, &
	4086	t_seri, u_seri, v_seri, &
	4087	zulow, zvlow, zustrli, zvstrli, &
	4088	d_t_lif, d_u_lif, d_v_lif)
	4089	ENDIF
[1638]	4090
[2469]	4091	! ajout des tendances de la portance de l'orographie
	4092	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
[2611]	4093	'lif', abortphy,flag_inhib_tend)
[2469]	4094	ENDIF ! fin de test sur ok_orolf
[1638]	4095
[2469]	4096	IF (ok_hines) then
	4097	! HINES GWD PARAMETRIZATION
	4098	east_gwstress=0.
	4099	west_gwstress=0.
	4100	du_gwd_hines=0.
	4101	dv_gwd_hines=0.
	4102	CALL hines_gwd(klon, klev, dtime, paprs, pplay, latitude_deg, t_seri, &
	4103	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
	4104	du_gwd_hines, dv_gwd_hines)
	4105	zustr_gwd_hines=0.
	4106	zvstr_gwd_hines=0.
	4107	DO k = 1, klev
	4108	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/dtime &
	4109	* (paprs(:, k)-paprs(:, k+1))/rg
	4110	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/dtime &
	4111	* (paprs(:, k)-paprs(:, k+1))/rg
	4112	ENDDO
[1001]	4113
[2469]	4114	d_t_hin(:, :)=0.
	4115	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
[2611]	4116	dqi0, paprs, 'hin', abortphy,flag_inhib_tend)
[2469]	4117	ENDIF
[2333]	4118
[2469]	4119	IF (.not. ok_hines .and. ok_gwd_rando) then
	4120	CALL acama_GWD_rando(DTIME, pplay, latitude_deg, t_seri, u_seri, &
	4121	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
	4122	dv_gwd_front, east_gwstress, west_gwstress)
	4123	zustr_gwd_front=0.
	4124	zvstr_gwd_front=0.
	4125	DO k = 1, klev
	4126	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/dtime &
	4127	* (paprs(:, k)-paprs(:, k+1))/rg
	4128	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/dtime &
	4129	* (paprs(:, k)-paprs(:, k+1))/rg
	4130	ENDDO
[644]	4131
[2469]	4132	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
[2611]	4133	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend)
[2469]	4134	ENDIF
[1938]	4135
[2469]	4136	if (ok_gwd_rando) then
	4137	call FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
	4138	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
	4139	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
	4140	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
[2611]	4141	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend)
[2469]	4142	zustr_gwd_rando=0.
	4143	zvstr_gwd_rando=0.
	4144	DO k = 1, klev
	4145	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/dtime &
	4146	* (paprs(:, k)-paprs(:, k+1))/rg
	4147	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/dtime &
	4148	* (paprs(:, k)-paprs(:, k+1))/rg
	4149	ENDDO
	4150	end if
[766]	4151
[2469]	4152	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
[1279]	4153
[2469]	4154	if (mydebug) then
	4155	call writefield_phy('u_seri',u_seri,nbp_lev)
	4156	call writefield_phy('v_seri',v_seri,nbp_lev)
	4157	call writefield_phy('t_seri',t_seri,nbp_lev)
	4158	call writefield_phy('q_seri',q_seri,nbp_lev)
	4159	endif
[2136]	4160
[2469]	4161	DO i = 1, klon
	4162	zustrph(i)=0.
	4163	zvstrph(i)=0.
	4164	ENDDO
	4165	DO k = 1, klev
	4166	DO i = 1, klon
	4167	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
	4168	(paprs(i,k)-paprs(i,k+1))/rg
	4169	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
	4170	(paprs(i,k)-paprs(i,k+1))/rg
	4171	ENDDO
	4172	ENDDO
	4173	!
	4174	!IM calcul composantes axiales du moment angulaire et couple des montagnes
	4175	!
	4176	IF (is_sequential .and. ok_orodr) THEN
	4177	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
	4178	ra,rg,romega, &
	4179	latitude_deg,longitude_deg,pphis, &
	4180	zustrdr,zustrli,zustrph, &
	4181	zvstrdr,zvstrli,zvstrph, &
	4182	paprs,u,v, &
	4183	aam, torsfc)
	4184	ENDIF
	4185	!IM cf. FLott END
	4186	!IM
	4187	IF (ip_ebil_phy.ge.2) THEN
	4188	ztit='after orography'
	4189	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
	4190	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	4191	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	4192	call diagphy(cell_area,ztit,ip_ebil_phy &
	4193	, zero_v, zero_v, zero_v, zero_v, zero_v &
	4194	, zero_v, zero_v, zero_v, ztsol &
	4195	, d_h_vcol, d_qt, d_ec &
	4196	, fs_bound, fq_bound )
	4197	END IF
[1863]	4198
[2469]	4199	!DC Calcul de la tendance due au methane
	4200	IF(ok_qch4) THEN
	4201	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
	4202	! ajout de la tendance d'humidite due au methane
	4203	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4*dtime, dql0, dqi0, paprs, &
[2611]	4204	'q_ch4', abortphy,flag_inhib_tend)
[2469]	4205	END IF
	4206	!
	4207	!
	4208	!====================================================================
	4209	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
	4210	!====================================================================
	4211	! Abderrahmane 24.08.09
	4212
	4213	IF (ok_cosp) THEN
	4214	! adeclarer
[1279]	4215	#ifdef CPP_COSP
[2469]	4216	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
[1279]	4217
[2469]	4218	IF (prt_level .GE.10) THEN
	4219	print*,'freq_cosp',freq_cosp
	4220	ENDIF
	4221	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
	4222	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
	4223	! s ref_liq,ref_ice
	4224	call phys_cosp(itap,dtime,freq_cosp, &
	4225	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
	4226	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
	4227	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
	4228	JrNt,ref_liq,ref_ice, &
	4229	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
	4230	zu10m,zv10m,pphis, &
	4231	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
	4232	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
	4233	prfl(:,1:klev),psfl(:,1:klev), &
	4234	pmflxr(:,1:klev),pmflxs(:,1:klev), &
	4235	mr_ozone,cldtau, cldemi)
[1412]	4236
[2469]	4237	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
	4238	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
	4239	! M clMISR,
	4240	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
	4241	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
[1279]	4242
[2469]	4243	ENDIF
[1279]	4244
	4245	#endif
[2469]	4246	ENDIF !ok_cosp
[2580]	4247
	4248
	4249	! Marine
	4250
	4251	IF (ok_airs) then
	4252
	4253	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/dtime)).EQ.0) THEN
	4254	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', &
	4255	& ok_airs, freq_airs
	4256	call simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
	4257	& map_prop_hc,map_prop_hist,&
	4258	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
	4259	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
	4260	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
	4261	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
	4262	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
	4263	& map_ntot,map_hc,map_hist,&
	4264	& map_Cb,map_ThCi,map_Anv,&
	4265	& alt_tropo )
	4266	ENDIF
	4267
	4268	ENDIF ! ok_airs
	4269
	4270
[2469]	4271	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	4272	!AA
	4273	!AA Installation de l'interface online-offline pour traceurs
	4274	!AA
	4275	!====================================================================
	4276	! Calcul des tendances traceurs
	4277	!====================================================================
	4278	!
[959]	4279
[2469]	4280	IF (type_trac=='repr') THEN
	4281	sh_in(:,:) = q_seri(:,:)
	4282	ELSE
	4283	sh_in(:,:) = qx(:,:,ivap)
[2907]	4284	ch_in(:,:) = qx(:,:,iliq)
[2469]	4285	END IF
[1565]	4286
[2469]	4287	call phytrac ( &
	4288	itap, days_elapsed+1, jH_cur, debut, &
	4289	lafin, dtime, u, v, t, &
	4290	paprs, pplay, pmfu, pmfd, &
	4291	pen_u, pde_u, pen_d, pde_d, &
	4292	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
	4293	u1, v1, ftsol, pctsrf, &
	4294	zustar, zu10m, zv10m, &
	4295	wstar(:,is_ave), ale_bl, ale_wake, &
	4296	latitude_deg, longitude_deg, &
	4297	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
	4298	presnivs, pphis, pphi, albsol1, &
[2907]	4299	sh_in, ch_in, rhcl, cldfra, rneb, &
[2469]	4300	diafra, cldliq, itop_con, ibas_con, &
	4301	pmflxr, pmflxs, prfl, psfl, &
	4302	da, phi, mp, upwd, &
	4303	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
	4304	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
	4305	ev, ep, epmlmMm, eplaMm, & !<<RomP
	4306	dnwd, aerosol_couple, flxmass_w, &
	4307	tau_aero, piz_aero, cg_aero, ccm, &
	4308	rfname, &
	4309	d_tr_dyn, & !<<RomP
[3433]	4310	tr_seri, init_source)
[524]	4311
[2469]	4312	IF (offline) THEN
[524]	4313
[2469]	4314	IF (prt_level.ge.9) &
	4315	print*,'Attention on met a 0 les thermiques pour phystoke'
	4316	call phystokenc ( &
	4317	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
	4318	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
	4319	fm_therm,entr_therm, &
	4320	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
	4321	frac_impa, frac_nucl, &
	4322	pphis,cell_area,dtime,itap, &
	4323	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
[524]	4324
	4325
[2469]	4326	ENDIF
[524]	4327
[2469]	4328	!
	4329	! Calculer le transport de l'eau et de l'energie (diagnostique)
	4330	!
	4331	CALL transp (paprs,zxtsol, &
	4332	t_seri, q_seri, u_seri, v_seri, zphi, &
	4333	ve, vq, ue, uq)
	4334	!
	4335	!IM global posePB BEG
	4336	IF(1.EQ.0) THEN
	4337	!
	4338	CALL transp_lay (paprs,zxtsol, &
	4339	t_seri, q_seri, u_seri, v_seri, zphi, &
	4340	ve_lay, vq_lay, ue_lay, uq_lay)
	4341	!
	4342	ENDIF !(1.EQ.0) THEN
	4343	!IM global posePB END
	4344	! Accumuler les variables a stocker dans les fichiers histoire:
	4345	!
[1279]	4346
[2469]	4347	!================================================================
	4348	! Conversion of kinetic and potential energy into heat, for
	4349	! parameterisation of subgrid-scale motions
	4350	!================================================================
[1753]	4351
[2469]	4352	d_t_ec(:,:)=0.
	4353	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
	4354	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
	4355	u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
	4356	zmasse,exner,d_t_ec)
	4357	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
[1753]	4358
[2469]	4359	!IM
	4360	IF (ip_ebil_phy.ge.1) THEN
	4361	ztit='after physic'
	4362	CALL diagetpq(cell_area,ztit,ip_ebil_phy,1,1,dtime &
	4363	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
	4364	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
	4365	! Comme les tendances de la physique sont ajoute dans la dynamique,
	4366	! on devrait avoir que la variation d'entalpie par la dynamique
	4367	! est egale a la variation de la physique au pas de temps precedent.
	4368	! Donc la somme de ces 2 variations devrait etre nulle.
[1279]	4369
[2469]	4370	call diagphy(cell_area,ztit,ip_ebil_phy &
	4371	, topsw, toplw, solsw, sollw, sens &
	4372	, evap, rain_fall, snow_fall, ztsol &
	4373	, d_h_vcol, d_qt, d_ec &
	4374	, fs_bound, fq_bound )
	4375	!
	4376	d_h_vcol_phy=d_h_vcol
	4377	!
	4378	END IF
	4379	!
	4380	!=======================================================================
	4381	! SORTIES
	4382	!=======================================================================
	4383	!
	4384	!IM initialisation + calculs divers diag AMIP2
	4385	!
	4386	include "calcul_divers.h"
	4387	!
	4388	!IM Interpolation sur les niveaux de pression du NMC
	4389	! -------------------------------------------------
[2271]	4390	#ifdef CPP_XIOS
[2469]	4391	!$OMP MASTER
	4392	!On recupere la valeur de la missing value donnee dans le xml
	4393	CALL xios_get_field_attr("t850",default_value=missing_val_omp)
	4394	! PRINT *,"ARNAUD value missing ",missing_val_omp
	4395	!$OMP END MASTER
	4396	!$OMP BARRIER
	4397	missing_val=missing_val_omp
[2271]	4398	#endif
	4399	#ifndef CPP_XIOS
[2469]	4400	missing_val=missing_val_nf90
[2271]	4401	#endif
[2469]	4402	!
	4403	include "calcul_STDlev.h"
	4404	!
	4405	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
	4406	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
	4407	!
[2496]	4408	!cc prw = eau precipitable
	4409	! prlw = colonne eau liquide
	4410	! prlw = colonne eau solide
[2499]	4411	prw(:) = 0.
	4412	prlw(:) = 0.
	4413	prsw(:) = 0.
	4414	DO k = 1, klev
	4415	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
	4416	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
	4417	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
[2469]	4418	ENDDO
	4419	!
	4420	IF (type_trac == 'inca') THEN
[655]	4421	#ifdef INCA
[2469]	4422	CALL VTe(VTphysiq)
	4423	CALL VTb(VTinca)
[959]	4424
[2469]	4425	CALL chemhook_end ( &
	4426	dtime, &
	4427	pplay, &
	4428	t_seri, &
	4429	tr_seri, &
	4430	nbtr, &
	4431	paprs, &
	4432	q_seri, &
	4433	cell_area, &
	4434	pphi, &
	4435	pphis, &
[2907]	4436	zx_rh, &
[3629]	4437	aps, bps, ap, bp)
[959]	4438
[2469]	4439	CALL VTe(VTinca)
	4440	CALL VTb(VTphysiq)
[655]	4441	#endif
[2469]	4442	END IF
[655]	4443
[1753]	4444
[2469]	4445	!
	4446	! Convertir les incrementations en tendances
	4447	!
	4448	IF (prt_level .GE.10) THEN
	4449	print *,'Convertir les incrementations en tendances '
	4450	ENDIF
	4451	!
	4452	if (mydebug) then
	4453	call writefield_phy('u_seri',u_seri,nbp_lev)
	4454	call writefield_phy('v_seri',v_seri,nbp_lev)
	4455	call writefield_phy('t_seri',t_seri,nbp_lev)
	4456	call writefield_phy('q_seri',q_seri,nbp_lev)
	4457	endif
[766]	4458
[2469]	4459	DO k = 1, klev
	4460	DO i = 1, klon
	4461	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
	4462	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
	4463	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
	4464	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
	4465	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
	4466	!CR: on ajoute le contenu en glace
	4467	if (nqo.eq.3) then
	4468	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
	4469	endif
	4470	ENDDO
	4471	ENDDO
	4472	!
	4473	!CR: nb de traceurs eau: nqo
	4474	! IF (nqtot.GE.3) THEN
	4475	IF (nqtot.GE.(nqo+1)) THEN
	4476	! DO iq = 3, nqtot
	4477	DO iq = nqo+1, nqtot
	4478	DO k = 1, klev
	4479	DO i = 1, klon
	4480	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
	4481	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
	4482	ENDDO
	4483	ENDDO
	4484	ENDDO
	4485	ENDIF
	4486	!
	4487	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
	4488	!IM global posePB include "write_bilKP_ins.h"
	4489	!IM global posePB include "write_bilKP_ave.h"
	4490	!
[1412]	4491
[2489]	4492	!--OB mass fixer
	4493	!--profile is corrected to force mass conservation of water
	4494	IF (mass_fixer) THEN
	4495	qql2(:)=0.0
[2499]	4496	DO k = 1, klev
	4497	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
[2489]	4498	ENDDO
	4499	DO i = 1, klon
	4500	!--compute ratio of what q+ql should be with conservation to what it is
	4501	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
	4502	DO k = 1, klev
	4503	q_seri(i,k) =q_seri(i,k)*corrqql
	4504	ql_seri(i,k)=ql_seri(i,k)*corrqql
	4505	ENDDO
	4506	ENDDO
	4507	ENDIF
	4508	!--fin mass fixer
	4509
[2469]	4510	! Sauvegarder les valeurs de t et q a la fin de la physique:
	4511	!
[2499]	4512	u_ancien(:,:) = u_seri(:,:)
	4513	v_ancien(:,:) = v_seri(:,:)
	4514	t_ancien(:,:) = t_seri(:,:)
	4515	q_ancien(:,:) = q_seri(:,:)
	4516	ql_ancien(:,:) = ql_seri(:,:)
	4517	qs_ancien(:,:) = qs_seri(:,:)
	4518	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
	4519	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
	4520	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
[2469]	4521	! !! RomP >>>
	4522	!CR: nb de traceurs eau: nqo
[2499]	4523	IF (nqtot.GT.nqo) THEN
[2469]	4524	DO iq = nqo+1, nqtot
[2499]	4525	tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
[2469]	4526	ENDDO
	4527	ENDIF
	4528	! !! RomP <<<
	4529	!==========================================================================
	4530	! Sorties des tendances pour un point particulier
	4531	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
	4532	! pour le debug
	4533	! La valeur de igout est attribuee plus haut dans le programme
	4534	!==========================================================================
[879]	4535
[2469]	4536	if (prt_level.ge.1) then
	4537	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
	4538	write(lunout,*) &
	4539	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
	4540	write(lunout,*) &
	4541	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
	4542	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
	4543	pctsrf(igout,is_sic)
	4544	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
	4545	do k=1,klev
	4546	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
	4547	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
	4548	d_t_eva(igout,k)
	4549	enddo
	4550	write(lunout,*) 'cool,heat'
	4551	do k=1,klev
	4552	write(lunout,*) cool(igout,k),heat(igout,k)
	4553	enddo
[879]	4554
[2469]	4555	!jyg< (En attendant de statuer sur le sort de d_t_oli)
	4556	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
	4557	!jyg! do k=1,klev
	4558	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
	4559	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
	4560	!jyg! enddo
	4561	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
	4562	do k=1,klev
	4563	write(lunout,*) d_t_vdf(igout,k), &
	4564	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
	4565	enddo
	4566	!>jyg
[879]	4567
[2469]	4568	write(lunout,*) 'd_ps ',d_ps(igout)
	4569	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
	4570	do k=1,klev
	4571	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
	4572	d_qx(igout,k,1),d_qx(igout,k,2)
	4573	enddo
	4574	endif
[879]	4575
[2469]	4576	!==========================================================================
[879]	4577
[2469]	4578	!============================================================
	4579	! Calcul de la temperature potentielle
	4580	!============================================================
	4581	DO k = 1, klev
	4582	DO i = 1, klon
	4583	!JYG/IM theta en debut du pas de temps
	4584	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
	4585	!JYG/IM theta en fin de pas de temps de physique
	4586	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
	4587	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
	4588	! MPL 20130625
	4589	! fth_fonctions.F90 et parkind1.F90
	4590	! sinon thetal=theta
	4591	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
	4592	! : ql_seri(i,k))
	4593	thetal(i,k)=theta(i,k)
	4594	ENDDO
	4595	ENDDO
	4596	!
[879]	4597
[2469]	4598	! 22.03.04 BEG
	4599	!=============================================================
	4600	! Ecriture des sorties
	4601	!=============================================================
[524]	4602	#ifdef CPP_IOIPSL
	4603
[2469]	4604	! Recupere des varibles calcule dans differents modules
	4605	! pour ecriture dans histxxx.nc
[782]	4606
[2469]	4607	! Get some variables from module fonte_neige_mod
	4608	CALL fonte_neige_get_vars(pctsrf, &
[2517]	4609	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
[782]	4610
[1507]	4611
[2469]	4612	!=============================================================
	4613	! Separation entre thermiques et non thermiques dans les sorties
	4614	! de fisrtilp
	4615	!=============================================================
[1507]	4616
[2469]	4617	if (iflag_thermals>=1) then
	4618	d_t_lscth=0.
	4619	d_t_lscst=0.
	4620	d_q_lscth=0.
	4621	d_q_lscst=0.
	4622	do k=1,klev
	4623	do i=1,klon
	4624	if (ptconvth(i,k)) then
	4625	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
	4626	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
	4627	else
	4628	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
	4629	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
	4630	endif
	4631	enddo
	4632	enddo
[1507]	4633
[2469]	4634	do i=1,klon
	4635	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
	4636	plul_th(i)=prfl(i,1)+psfl(i,1)
	4637	enddo
	4638	endif
[909]	4639
[1791]	4640
[2469]	4641	!On effectue les sorties:
[1791]	4642
[2469]	4643	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
	4644	pplay, lmax_th, aerosol_couple, &
[2496]	4645	ok_ade, ok_aie, ivap, iliq, isol, new_aod, &
	4646	ok_sync, ptconv, read_climoz, clevSTD, &
[2469]	4647	ptconvth, d_t, qx, d_qx, zmasse, &
	4648	flag_aerosol, flag_aerosol_strat, ok_cdnc)
[1791]	4649
[2469]	4650	include "write_histday_seri.h"
[687]	4651
[2590]	4652	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
[687]	4653
[524]	4654	#endif
	4655
[2235]	4656
[2469]	4657	!====================================================================
	4658	! Arret du modele apres hgardfou en cas de detection d'un
	4659	! plantage par hgardfou
	4660	!====================================================================
[2235]	4661
	4662	IF (abortphy==1) THEN
	4663	abort_message ='Plantage hgardfou'
[2311]	4664	CALL abort_physic (modname,abort_message,1)
[2235]	4665	ENDIF
	4666
[2469]	4667	! 22.03.04 END
	4668	!
	4669	!====================================================================
	4670	! Si c'est la fin, il faut conserver l'etat de redemarrage
	4671	!====================================================================
	4672	!
[782]	4673
[2469]	4674	IF (lafin) THEN
	4675	itau_phy = itau_phy + itap
	4676	CALL phyredem ("restartphy.nc")
	4677	! open(97,form="unformatted",file="finbin")
	4678	! write(97) u_seri,v_seri,t_seri,q_seri
	4679	! close(97)
	4680	!$OMP MASTER
	4681	if (read_climoz >= 1) then
	4682	if (is_mpi_root) then
	4683	call nf95_close(ncid_climoz)
	4684	end if
	4685	deallocate(press_climoz) ! pointer
	4686	end if
	4687	!$OMP END MASTER
	4688	ENDIF
[1863]	4689
[2469]	4690	! first=.false.
[1863]	4691
[2418]	4692
[2469]	4693	END SUBROUTINE physiq
[2418]	4694
	4695	END MODULE physiq_mod

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

Download in other formats:

Original Format