Context Navigation

physiq_mod.F90 @ 4135

Last change on this file since 4135 was 4135, checked in by Laurent Fairhead, 2 years ago
Cleaning up the USEs in physiq_mod to prepare inclusion of acc data regions and got rid of the last getin in physiq_mod.F90 :-)
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 `Id`
File size: 195.3 KB

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

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

Download in other formats:

Original Format