Context Navigation

physiq_mod.F90 @ 2994

Last change on this file since 2994 was 2994, checked in by oboucher, 7 years ago
Addition to previous commit. Call of stratomask to diagnose tropopause level according to WMO definition Previous commit aimed at diagnosing tropospheric and stratospheric O3 columns in DU.
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: 162.0 KB

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

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

Download in other formats:

Original Format