Context Navigation

physiq_mod.F90 @ 4137

Last change on this file since 4137 was 4137, checked in by Ehouarn Millour, 2 years ago
Minor follow-up to r4135. EM
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 4137 2022-04-22 16:01:52Z emillour $
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
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	! FH : if ok_sync=.true. , the time axis is written at each time step
1675	! in the output files. Only at the end in the opposite case
1676	ok_sync=.FALSE.
1677	CALL getin_p('ok_sync',ok_sync)
1678	!$OMP MASTER
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)
1689	!$OMP END MASTER
1690	!$OMP BARRIER
1691
1692	freq_outNMC(1) = ecrit_files(7)
1693	freq_outNMC(2) = ecrit_files(8)
1694	freq_outNMC(3) = ecrit_files(9)
1695	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1696	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1697	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1698
1699	#ifndef CPP_XIOS
1700	CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
1701	#endif
1702
1703	#endif
1704	ecrit_reg = ecrit_reg * un_jour
1705	ecrit_tra = ecrit_tra * un_jour
1706
1707	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
1708	date0 = jD_ref
1709	WRITE(,) 'physiq date0 : ',date0
1710	!
1711
1712	! CALL create_climoz(read_climoz)
1713	IF (.NOT. create_etat0_limit) CALL init_aero_fromfile(flag_aerosol) !! initialise aero from file for XIOS interpolation (unstructured_grid)
1714	IF (.NOT. create_etat0_limit) CALL init_readaerosolstrato(flag_aerosol_strat) !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
1715
1716	#ifdef CPP_COSP
1717	IF (ok_cosp) THEN
1718	! DO k = 1, klev
1719	! DO i = 1, klon
1720	! phicosp(i,k) = pphi(i,k) + pphis(i)
1721	! ENDDO
1722	! ENDDO
1723	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
1724	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1725	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1726	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1727	JrNt,ref_liq,ref_ice, &
1728	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1729	zu10m,zv10m,pphis, &
1730	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1731	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1732	prfl(:,1:klev),psfl(:,1:klev), &
1733	pmflxr(:,1:klev),pmflxs(:,1:klev), &
1734	mr_ozone,cldtau, cldemi)
1735	ENDIF
1736	#endif
1737
1738	#ifdef CPP_COSP2
1739	IF (ok_cosp) THEN
1740	! DO k = 1, klev
1741	! DO i = 1, klon
1742	! phicosp(i,k) = pphi(i,k) + pphis(i)
1743	! ENDDO
1744	! ENDDO
1745	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
1746	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1747	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1748	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1749	JrNt,ref_liq,ref_ice, &
1750	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1751	zu10m,zv10m,pphis, &
1752	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1753	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1754	prfl(:,1:klev),psfl(:,1:klev), &
1755	pmflxr(:,1:klev),pmflxs(:,1:klev), &
1756	mr_ozone,cldtau, cldemi)
1757	ENDIF
1758	#endif
1759
1760	#ifdef CPP_COSPV2
1761	IF (ok_cosp) THEN
1762	DO k = 1, klev
1763	DO i = 1, klon
1764	phicosp(i,k) = pphi(i,k) + pphis(i)
1765	ENDDO
1766	ENDDO
1767	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
1768	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1769	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1770	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1771	JrNt,ref_liq,ref_ice, &
1772	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1773	zu10m,zv10m,pphis, &
1774	phicosp,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1775	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1776	prfl(:,1:klev),psfl(:,1:klev), &
1777	pmflxr(:,1:klev),pmflxs(:,1:klev), &
1778	mr_ozone,cldtau, cldemi)
1779	ENDIF
1780	#endif
1781
1782	!
1783	!
1784	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1785	! Nouvelle initialisation pour le rayonnement RRTM
1786	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1787
1788	CALL iniradia(klon,klev,paprs(1,1:klev+1))
1789
1790	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1791	CALL wake_ini(rg,rd,rv,prt_level)
1792	CALL thermcell_ini(iflag_thermals,prt_level,tau_thermals,lunout, &
1793	& RG,RD,RCPD,RKAPPA,RLVTT,RETV)
1794
1795	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1796
1797	!
1798	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1799	! Initialisation des champs dans phytrac* qui sont utilises par phys_output_write*
1800	!
1801	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1802
1803	#ifdef CPP_Dust
1804	! Quand on utilise SPLA, on force iflag_phytrac=1
1805	CALL phytracr_spl_out_init()
1806	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
1807	pplay, lmax_th, aerosol_couple, &
1808	ok_ade, ok_aie, ivap, ok_sync, &
1809	ptconv, read_climoz, clevSTD, &
1810	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
1811	flag_aerosol, flag_aerosol_strat, ok_cdnc)
1812	#else
1813	! phys_output_write écrit des variables traceurs seulement si iflag_phytrac == 1
1814	! donc seulement dans ce cas on doit appeler phytrac_init()
1815	IF (iflag_phytrac == 1 ) THEN
1816	CALL phytrac_init()
1817	ENDIF
1818	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
1819	pplay, lmax_th, aerosol_couple, &
1820	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ok_sync,&
1821	ptconv, read_climoz, clevSTD, &
1822	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
1823	flag_aerosol, flag_aerosol_strat, ok_cdnc)
1824	#endif
1825
1826
1827	#ifdef CPP_XIOS
1828	IF (is_omp_master) CALL xios_update_calendar(1)
1829	#endif
1830	IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
1831	CALL create_etat0_limit_unstruct
1832	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1833
1834	!jyg<
1835	IF (iflag_pbl<=1) THEN
1836	! No TKE for Standard Physics
1837	pbl_tke(:,:,:)=0.
1838
1839	ELSE IF (klon_glo==1) THEN
1840	pbl_tke(:,:,is_ave) = 0.
1841	DO nsrf=1,nbsrf
1842	DO k = 1,klev+1
1843	pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
1844	+pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
1845	ENDDO
1846	ENDDO
1847	ELSE
1848	pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
1849	!>jyg
1850	ENDIF
1851	!IM begin
1852	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1853	,ratqs(1,1)
1854	!IM end
1855
1856
1857	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1858	!
1859	! on remet le calendrier a zero
1860	!
1861	IF (raz_date .eq. 1) THEN
1862	itau_phy = 0
1863	ENDIF
1864
1865	! IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
1866	! WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
1867	! pdtphys
1868	! abort_message='Pas physique n est pas correct '
1869	! ! call abort_physic(modname,abort_message,1)
1870	! phys_tstep=pdtphys
1871	! ENDIF
1872	IF (nlon .NE. klon) THEN
1873	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1874	klon
1875	abort_message='nlon et klon ne sont pas coherents'
1876	CALL abort_physic(modname,abort_message,1)
1877	ENDIF
1878	IF (nlev .NE. klev) THEN
1879	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1880	klev
1881	abort_message='nlev et klev ne sont pas coherents'
1882	CALL abort_physic(modname,abort_message,1)
1883	ENDIF
1884	!
1885	IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
1886	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1887	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1888	abort_message='Nbre d appels au rayonnement insuffisant'
1889	CALL abort_physic(modname,abort_message,1)
1890	ENDIF
1891
1892	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1893	! Initialisation pour la convection de K.E. et pour les poches froides
1894	!
1895	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1896
1897	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1898	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", ok_cvl
1899	!
1900	!KE43
1901	! Initialisation pour la convection de K.E. (sb):
1902	IF (iflag_con.GE.3) THEN
1903
1904	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1905	WRITE(lunout,*) &
1906	"On va utiliser le melange convectif des traceurs qui"
1907	WRITE(lunout,*)"est calcule dans convect4.3"
1908	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1909
1910	DO i = 1, klon
1911	ema_cbmf(i) = 0.
1912	ema_pcb(i) = 0.
1913	ema_pct(i) = 0.
1914	! ema_workcbmf(i) = 0.
1915	ENDDO
1916	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1917	DO i = 1, klon
1918	ibas_con(i) = 1
1919	itop_con(i) = 1
1920	ENDDO
1921	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1922	!================================================================
1923	!CR:04.12.07: initialisations poches froides
1924	! Controle de ALE et ALP pour la fermeture convective (jyg)
1925	IF (iflag_wake>=1) THEN
1926	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1927	,alp_bl_prescr, ale_bl_prescr)
1928	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1929	! print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
1930	!
1931	! Initialize tendencies of wake state variables (for some flag values
1932	! they are not computed).
1933	d_deltat_wk(:,:) = 0.
1934	d_deltaq_wk(:,:) = 0.
1935	d_deltat_wk_gw(:,:) = 0.
1936	d_deltaq_wk_gw(:,:) = 0.
1937	d_deltat_vdf(:,:) = 0.
1938	d_deltaq_vdf(:,:) = 0.
1939	d_deltat_the(:,:) = 0.
1940	d_deltaq_the(:,:) = 0.
1941	d_deltat_ajs_cv(:,:) = 0.
1942	d_deltaq_ajs_cv(:,:) = 0.
1943	d_s_wk(:) = 0.
1944	d_dens_wk(:) = 0.
1945	ENDIF ! (iflag_wake>=1)
1946
1947	! do i = 1,klon
1948	! Ale_bl(i)=0.
1949	! Alp_bl(i)=0.
1950	! enddo
1951
1952	!ELSE
1953	! ALLOCATE(tabijGCM(0))
1954	! ALLOCATE(lonGCM(0), latGCM(0))
1955	! ALLOCATE(iGCM(0), jGCM(0))
1956	ENDIF ! (iflag_con.GE.3)
1957	!
1958	DO i=1,klon
1959	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1960	ENDDO
1961
1962	!34EK
1963	IF (ok_orodr) THEN
1964
1965	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1966	! FH sans doute a enlever de finitivement ou, si on le
1967	! garde, l'activer justement quand ok_orodr = false.
1968	! ce rugoro est utilise par la couche limite et fait double emploi
1969	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1970	! DO i=1,klon
1971	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1972	! ENDDO
1973	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1974	IF (ok_strato) THEN
1975	CALL SUGWD_strato(klon,klev,paprs,pplay)
1976	ELSE
1977	CALL SUGWD(klon,klev,paprs,pplay)
1978	ENDIF
1979
1980	DO i=1,klon
1981	zuthe(i)=0.
1982	zvthe(i)=0.
1983	IF (zstd(i).gt.10.) THEN
1984	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1985	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1986	ENDIF
1987	ENDDO
1988	ENDIF
1989	!
1990	!
1991	lmt_pas = NINT(86400./phys_tstep * 1.0) ! tous les jours
1992	WRITE(lunout,*)'La frequence de lecture surface est de ', &
1993	lmt_pas
1994	!
1995	capemaxcels = 't_max(X)'
1996	t2mincels = 't_min(X)'
1997	t2maxcels = 't_max(X)'
1998	tinst = 'inst(X)'
1999	tave = 'ave(X)'
2000	!IM cf. AM 081204 BEG
2001	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
2002	!IM cf. AM 081204 END
2003	!
2004	!=============================================================
2005	! Initialisation des sorties
2006	!=============================================================
2007
2008	#ifdef CPP_XIOS
2009	! Get "missing_val" value from XML files (from temperature variable)
2010	!$OMP MASTER
2011	CALL xios_get_field_attr("temp",default_value=missing_val_omp)
2012	!$OMP END MASTER
2013	!$OMP BARRIER
2014	missing_val=missing_val_omp
2015	#endif
2016
2017	#ifdef CPP_XIOS
2018	! Need to put this initialisation after phyetat0 as in the coupled model the XIOS context is only
2019	! initialised at that moment
2020	! Get "missing_val" value from XML files (from temperature variable)
2021	!$OMP MASTER
2022	CALL xios_get_field_attr("temp",default_value=missing_val_omp)
2023	!$OMP END MASTER
2024	!$OMP BARRIER
2025	missing_val=missing_val_omp
2026	!
2027	! Now we activate some double radiation call flags only if some
2028	! diagnostics are requested, otherwise there is no point in doing this
2029	IF (is_master) THEN
2030	!--setting up swaero_diag to TRUE in XIOS case
2031	IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
2032	xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
2033	xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR. &
2034	(iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
2035	xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0")))) &
2036	!!!--for now these fields are not in the XML files so they are omitted
2037	!!! xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
2038	swaero_diag=.TRUE.
2039
2040	!--setting up swaerofree_diag to TRUE in XIOS case
2041	IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
2042	xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR. &
2043	xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
2044	xios_field_is_active("LWupTOAcleanclr")) &
2045	swaerofree_diag=.TRUE.
2046
2047	!--setting up dryaod_diag to TRUE in XIOS case
2048	DO naero = 1, naero_tot-1
2049	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
2050	ENDDO
2051	!
2052	!--setting up ok_4xCO2atm to TRUE in XIOS case
2053	IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
2054	xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
2055	xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
2056	xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
2057	xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
2058	xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
2059	ok_4xCO2atm=.TRUE.
2060	ENDIF
2061	!$OMP BARRIER
2062	CALL bcast(swaero_diag)
2063	CALL bcast(swaerofree_diag)
2064	CALL bcast(dryaod_diag)
2065	CALL bcast(ok_4xCO2atm)
2066	#endif
2067	!
2068	CALL printflag( tabcntr0,radpas,ok_journe, &
2069	ok_instan, ok_region )
2070	!
2071	!
2072	! Prescrire l'ozone dans l'atmosphere
2073	!
2074	!c DO i = 1, klon
2075	!c DO k = 1, klev
2076	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
2077	!c ENDDO
2078	!c ENDDO
2079	!
2080	IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN ! ModThL
2081	#ifdef INCA
2082	CALL VTe(VTphysiq)
2083	CALL VTb(VTinca)
2084	calday = REAL(days_elapsed) + jH_cur
2085	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
2086
2087	call init_const_lmdz( &
2088	ndays, nbsrf, is_oce,is_sic, is_ter,is_lic, calend, &
2089	config_inca)
2090
2091	CALL init_inca_geometry( &
2092	longitude, latitude, &
2093	boundslon, boundslat, &
2094	dx, dy, cell_area, ind_cell_glo)
2095
2096
2097	CALL chemini( pplay, &
2098	nbp_lon, nbp_lat, &
2099	latitude_deg, &
2100	longitude_deg, &
2101	presnivs, &
2102	calday, &
2103	klon, &
2104	nqtot, &
2105	nqo+nqCO2, &
2106	pdtphys, &
2107	annee_ref, &
2108	year_cur, &
2109	day_ref, &
2110	day_ini, &
2111	start_time, &
2112	itau_phy, &
2113	date0, &
2114	io_lon, &
2115	io_lat, &
2116	chemistry_couple, &
2117	init_source, &
2118	init_tauinca, &
2119	init_pizinca, &
2120	init_cginca, &
2121	init_ccminca)
2122
2123
2124	! initialisation des variables depuis le restart de inca
2125	ccm(:,:,:) = init_ccminca
2126	tau_aero(:,:,:,:) = init_tauinca
2127	piz_aero(:,:,:,:) = init_pizinca
2128	cg_aero(:,:,:,:) = init_cginca
2129	!
2130
2131
2132	CALL VTe(VTinca)
2133	CALL VTb(VTphysiq)
2134	#endif
2135	ENDIF
2136	!
2137	IF (type_trac == 'repr') THEN
2138	#ifdef REPROBUS
2139	CALL chemini_rep( &
2140	presnivs, &
2141	pdtphys, &
2142	annee_ref, &
2143	day_ref, &
2144	day_ini, &
2145	start_time, &
2146	itau_phy, &
2147	io_lon, &
2148	io_lat)
2149	#endif
2150	ENDIF
2151
2152	!$omp single
2153	IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, &
2154	press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
2155	!$omp end single
2156	!
2157	!IM betaCRF
2158	pfree=70000. !Pa
2159	beta_pbl=1.
2160	beta_free=1.
2161	lon1_beta=-180.
2162	lon2_beta=+180.
2163	lat1_beta=90.
2164	lat2_beta=-90.
2165	mskocean_beta=.FALSE.
2166
2167	!albedo SB >>>
2168	SELECT CASE(nsw)
2169	CASE(2)
2170	SFRWL(1)=0.45538747
2171	SFRWL(2)=0.54461211
2172	CASE(4)
2173	SFRWL(1)=0.45538747
2174	SFRWL(2)=0.32870591
2175	SFRWL(3)=0.18568763
2176	SFRWL(4)=3.02191470E-02
2177	CASE(6)
2178	SFRWL(1)=1.28432794E-03
2179	SFRWL(2)=0.12304168
2180	SFRWL(3)=0.33106142
2181	SFRWL(4)=0.32870591
2182	SFRWL(5)=0.18568763
2183	SFRWL(6)=3.02191470E-02
2184	END SELECT
2185	!albedo SB <<<
2186
2187	OPEN(99,file='beta_crf.data',status='old', &
2188	form='formatted',err=9999)
2189	READ(99,*,end=9998) pfree
2190	READ(99,*,end=9998) beta_pbl
2191	READ(99,*,end=9998) beta_free
2192	READ(99,*,end=9998) lon1_beta
2193	READ(99,*,end=9998) lon2_beta
2194	READ(99,*,end=9998) lat1_beta
2195	READ(99,*,end=9998) lat2_beta
2196	READ(99,*,end=9998) mskocean_beta
2197	9998 Continue
2198	CLOSE(99)
2199	9999 Continue
2200	WRITE(,)'pfree=',pfree
2201	WRITE(,)'beta_pbl=',beta_pbl
2202	WRITE(,)'beta_free=',beta_free
2203	WRITE(,)'lon1_beta=',lon1_beta
2204	WRITE(,)'lon2_beta=',lon2_beta
2205	WRITE(,)'lat1_beta=',lat1_beta
2206	WRITE(,)'lat2_beta=',lat2_beta
2207	WRITE(,)'mskocean_beta=',mskocean_beta
2208
2209	!lwoff=y : offset LW CRE for radiation code and other schemes
2210	!lwoff=y : betalwoff=1.
2211	betalwoff=0.
2212	IF (ok_lwoff) THEN
2213	betalwoff=1.
2214	ENDIF
2215	WRITE(,)'ok_lwoff=',ok_lwoff
2216	!
2217	!lwoff=y to begin only sollw and sollwdown are set up to CS values
2218	sollw = sollw + betalwoff * (sollw0 - sollw)
2219	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
2220	sollwdown(:))
2221
2222
2223
2224	ENDIF
2225	!
2226	! ************** Fin de IF ( debut ) *************
2227	!
2228	!
2229	! Incrementer le compteur de la physique
2230	!
2231	itap = itap + 1
2232	IF (is_master .OR. prt_level > 9) THEN
2233	IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
2234	WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
2235	WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
2236	100 FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
2237	ENDIF
2238	ENDIF
2239	!
2240	!
2241	! Update fraction of the sub-surfaces (pctsrf) and
2242	! initialize, where a new fraction has appeared, all variables depending
2243	! on the surface fraction.
2244	!
2245	CALL change_srf_frac(itap, phys_tstep, days_elapsed+1, &
2246	pctsrf, fevap, z0m, z0h, agesno, &
2247	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
2248
2249	! Update time and other variables in Reprobus
2250	IF (type_trac == 'repr') THEN
2251	#ifdef REPROBUS
2252	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
2253	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
2254	CALL Rtime(debut)
2255	#endif
2256	ENDIF
2257
2258	! Tendances bidons pour les processus qui n'affectent pas certaines
2259	! variables.
2260	du0(:,:)=0.
2261	dv0(:,:)=0.
2262	dt0 = 0.
2263	dq0(:,:)=0.
2264	dql0(:,:)=0.
2265	dqi0(:,:)=0.
2266	dsig0(:) = 0.
2267	ddens0(:) = 0.
2268	wkoccur1(:)=1
2269	!
2270	! Mettre a zero des variables de sortie (pour securite)
2271	!
2272	DO i = 1, klon
2273	d_ps(i) = 0.0
2274	ENDDO
2275	DO k = 1, klev
2276	DO i = 1, klon
2277	d_t(i,k) = 0.0
2278	d_u(i,k) = 0.0
2279	d_v(i,k) = 0.0
2280	ENDDO
2281	ENDDO
2282	DO iq = 1, nqtot
2283	DO k = 1, klev
2284	DO i = 1, klon
2285	d_qx(i,k,iq) = 0.0
2286	ENDDO
2287	ENDDO
2288	ENDDO
2289	beta_prec_fisrt(:,:)=0.
2290	beta_prec(:,:)=0.
2291	!
2292	! Output variables from the convective scheme should not be set to 0
2293	! since convection is not always called at every time step.
2294	IF (ok_bug_cv_trac) THEN
2295	da(:,:)=0.
2296	mp(:,:)=0.
2297	phi(:,:,:)=0.
2298	! RomP >>>
2299	phi2(:,:,:)=0.
2300	epmlmMm(:,:,:)=0.
2301	eplaMm(:,:)=0.
2302	d1a(:,:)=0.
2303	dam(:,:)=0.
2304	pmflxr(:,:)=0.
2305	pmflxs(:,:)=0.
2306	! RomP <<<
2307	ENDIF
2308	!
2309	! Ne pas affecter les valeurs entrees de u, v, h, et q
2310	!
2311	DO k = 1, klev
2312	DO i = 1, klon
2313	t_seri(i,k) = t(i,k)
2314	u_seri(i,k) = u(i,k)
2315	v_seri(i,k) = v(i,k)
2316	q_seri(i,k) = qx(i,k,ivap)
2317	ql_seri(i,k) = qx(i,k,iliq)
2318	!CR: ATTENTION, on rajoute la variable glace
2319	IF (nqo.EQ.2) THEN !--vapour and liquid only
2320	qs_seri(i,k) = 0.
2321	rneb_seri(i,k) = 0.
2322	ELSE IF (nqo.EQ.3) THEN !--vapour, liquid and ice
2323	qs_seri(i,k) = qx(i,k,isol)
2324	rneb_seri(i,k) = 0.
2325	ELSE IF (nqo.EQ.4) THEN !--vapour, liquid, ice and rneb
2326	qs_seri(i,k) = qx(i,k,isol)
2327	rneb_seri(i,k) = qx(i,k,irneb)
2328	ENDIF
2329	ENDDO
2330	ENDDO
2331	!
2332	!--OB mass fixer
2333	IF (mass_fixer) THEN
2334	!--store initial water burden
2335	qql1(:)=0.0
2336	DO k = 1, klev
2337	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
2338	ENDDO
2339	ENDIF
2340	!--fin mass fixer
2341
2342	tke0(:,:)=pbl_tke(:,:,is_ave)
2343	IF (nqtot > nqo) THEN
2344	! water isotopes are not included in tr_seri
2345	itr = 0
2346	DO iq = 1, nqtot
2347	IF(.NOT.tracers(iq)%isInPhysics) CYCLE
2348	itr = itr+1
2349	DO k = 1, klev
2350	DO i = 1, klon
2351	tr_seri(i,k,itr) = qx(i,k,iq)
2352	ENDDO
2353	ENDDO
2354	ENDDO
2355	ELSE
2356	! DC: make sure the final "1" index was meant for 1st H2O phase (vapor) !!!
2357	tr_seri(:,:,strIdx(tracers(:)%name,addPhase('H2O','g'))) = 0.0
2358	ENDIF
2359	!
2360	! Temporary solutions adressing ticket #104 and the non initialisation of tr_ancien
2361	! LF
2362	IF (debut) THEN
2363	WRITE(lunout,*)' WARNING: tr_ancien initialised to tr_seri'
2364	itr = 0
2365	do iq = 1, nqtot
2366	IF(.NOT.tracers(iq)%isInPhysics) CYCLE
2367	itr = itr+1
2368	tr_ancien(:,:,itr)=tr_seri(:,:,itr)
2369	enddo
2370	ENDIF
2371	!
2372	DO i = 1, klon
2373	ztsol(i) = 0.
2374	ENDDO
2375	DO nsrf = 1, nbsrf
2376	DO i = 1, klon
2377	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2378	ENDDO
2379	ENDDO
2380	! Initialize variables used for diagnostic purpose
2381	IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
2382
2383	! Diagnostiquer la tendance dynamique
2384	!
2385	IF (ancien_ok) THEN
2386	!
2387	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
2388	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
2389	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
2390	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
2391	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
2392	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
2393	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
2394	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
2395	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
2396	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
2397	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
2398	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
2399	! !! RomP >>> td dyn traceur
2400	IF (nqtot > nqo) d_tr_dyn(:,:,:)=(tr_seri(:,:,:)-tr_ancien(:,:,:))/phys_tstep
2401	! !! RomP <<<
2402	!!d_rneb_dyn(:,:)=(rneb_seri(:,:)-rneb_ancien(:,:))/phys_tstep
2403	d_rneb_dyn(:,:)=0.0
2404	ELSE
2405	d_u_dyn(:,:) = 0.0
2406	d_v_dyn(:,:) = 0.0
2407	d_t_dyn(:,:) = 0.0
2408	d_q_dyn(:,:) = 0.0
2409	d_ql_dyn(:,:) = 0.0
2410	d_qs_dyn(:,:) = 0.0
2411	d_q_dyn2d(:) = 0.0
2412	d_ql_dyn2d(:) = 0.0
2413	d_qs_dyn2d(:) = 0.0
2414	! !! RomP >>> td dyn traceur
2415	IF (nqtot > nqo) d_tr_dyn(:,:,:)= 0.0
2416	! !! RomP <<<
2417	d_rneb_dyn(:,:)=0.0
2418	ancien_ok = .TRUE.
2419	ENDIF
2420	!
2421	! Ajouter le geopotentiel du sol:
2422	!
2423	DO k = 1, klev
2424	DO i = 1, klon
2425	zphi(i,k) = pphi(i,k) + pphis(i)
2426	ENDDO
2427	ENDDO
2428	!
2429	! Verifier les temperatures
2430	!
2431	!IM BEG
2432	IF (check) THEN
2433	amn=MIN(ftsol(1,is_ter),1000.)
2434	amx=MAX(ftsol(1,is_ter),-1000.)
2435	DO i=2, klon
2436	amn=MIN(ftsol(i,is_ter),amn)
2437	amx=MAX(ftsol(i,is_ter),amx)
2438	ENDDO
2439	!
2440	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2441	ENDIF !(check) THEN
2442	!IM END
2443	!
2444	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
2445	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
2446
2447	!
2448	!IM BEG
2449	IF (check) THEN
2450	amn=MIN(ftsol(1,is_ter),1000.)
2451	amx=MAX(ftsol(1,is_ter),-1000.)
2452	DO i=2, klon
2453	amn=MIN(ftsol(i,is_ter),amn)
2454	amx=MAX(ftsol(i,is_ter),amx)
2455	ENDDO
2456	!
2457	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2458	ENDIF !(check) THEN
2459	!IM END
2460	!
2461	! Mettre en action les conditions aux limites (albedo, sst, etc.).
2462	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
2463	!
2464	! Update ozone if day change
2465	IF (MOD(itap-1,lmt_pas) == 0) THEN
2466	IF (read_climoz <= 0) THEN
2467	! Once per day, update ozone from Royer:
2468	IF (solarlong0<-999.) then
2469	! Generic case with evolvoing season
2470	zzz=real(days_elapsed+1)
2471	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
2472	! Particular case with annual mean insolation
2473	zzz=real(90) ! could be revisited
2474	IF (read_climoz/=-1) THEN
2475	abort_message ='read_climoz=-1 is recommended when ' &
2476	// 'solarlong0=1000.'
2477	CALL abort_physic (modname,abort_message,1)
2478	ENDIF
2479	ELSE
2480	! Case where the season is imposed with solarlong0
2481	zzz=real(90) ! could be revisited
2482	ENDIF
2483
2484	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
2485	#ifdef REPROBUS
2486	ptrop=dyn_tropopause(t_seri, ztsol, paprs, pplay, rot)/100.
2487	DO i = 1, klon
2488	Z1=t_seri(i,itroprep(i)+1)
2489	Z2=t_seri(i,itroprep(i))
2490	fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2491	B=Z2-fac*alog(pplay(i,itroprep(i)))
2492	ttrop(i)= fac*alog(ptrop(i))+B
2493	!
2494	Z1= 1.e-3 * ( pphi(i,itroprep(i)+1)+pphis(i) ) / gravit
2495	Z2= 1.e-3 * ( pphi(i,itroprep(i)) +pphis(i) ) / gravit
2496	fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2497	B=Z2-fac*alog(pplay(i,itroprep(i)))
2498	ztrop(i)=fac*alog(ptrop(i))+B
2499	ENDDO
2500	#endif
2501	ELSE
2502	!--- ro3i = elapsed days number since current year 1st january, 0h
2503	ro3i=days_elapsed+jh_cur-jh_1jan
2504	!--- scaling for old style files (360 records)
2505	IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
2506	IF(adjust_tropopause) THEN
2507	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2508	ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), &
2509	time_climoz , longitude_deg, latitude_deg, &
2510	dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
2511	ELSE
2512	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2513	ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), &
2514	time_climoz )
2515	ENDIF
2516	! Convert from mole fraction of ozone to column density of ozone in a
2517	! cell, in kDU:
2518	FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2519	* zmasse / dobson_u / 1e3
2520	! (By regridding ozone values for LMDZ only once a day, we
2521	! have already neglected the variation of pressure in one
2522	! day. So do not recompute "wo" at each time step even if
2523	! "zmasse" changes a little.)
2524	ENDIF
2525	ENDIF
2526	!
2527	! Re-evaporer l'eau liquide nuageuse
2528	!
2529	CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2530	& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
2531
2532	CALL add_phys_tend &
2533	(du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
2534	'eva',abortphy,flag_inhib_tend,itap,0)
2535	CALL prt_enerbil('eva',itap)
2536
2537	!=========================================================================
2538	! Calculs de l'orbite.
2539	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2540	! doit donc etre plac\'e avant radlwsw et pbl_surface
2541
2542	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2543	CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2544	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2545	!
2546	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2547	! solarlong0
2548	IF (solarlong0<-999.) THEN
2549	IF (new_orbit) THEN
2550	! calcul selon la routine utilisee pour les planetes
2551	CALL solarlong(day_since_equinox, zlongi, dist)
2552	ELSE
2553	! calcul selon la routine utilisee pour l'AR4
2554	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2555	ENDIF
2556	ELSE
2557	zlongi=solarlong0 ! longitude solaire vraie
2558	dist=1. ! distance au soleil / moyenne
2559	ENDIF
2560
2561	IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2562
2563
2564	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2565	! Calcul de l'ensoleillement :
2566	! ============================
2567	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2568	! l'annee a partir d'une formule analytique.
2569	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2570	! non nul aux poles.
2571	IF (abs(solarlong0-1000.)<1.e-4) THEN
2572	CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
2573	latitude_deg,longitude_deg,rmu0,fract)
2574	swradcorr(:) = 1.0
2575	JrNt(:) = 1.0
2576	zrmu0(:) = rmu0(:)
2577	ELSE
2578	! recode par Olivier Boucher en sept 2015
2579	SELECT CASE (iflag_cycle_diurne)
2580	CASE(0)
2581	! Sans cycle diurne
2582	CALL angle(zlongi, latitude_deg, fract, rmu0)
2583	swradcorr = 1.0
2584	JrNt = 1.0
2585	zrmu0 = rmu0
2586	CASE(1)
2587	! Avec cycle diurne sans application des poids
2588	! bit comparable a l ancienne formulation cycle_diurne=true
2589	! on integre entre gmtime et gmtime+radpas
2590	zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
2591	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2592	latitude_deg,longitude_deg,rmu0,fract)
2593	zrmu0 = rmu0
2594	swradcorr = 1.0
2595	! Calcul du flag jour-nuit
2596	JrNt = 0.0
2597	WHERE (fract.GT.0.0) JrNt = 1.0
2598	CASE(2)
2599	! Avec cycle diurne sans application des poids
2600	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2601	! Comme cette routine est appele a tous les pas de temps de
2602	! la physique meme si le rayonnement n'est pas appele je
2603	! remonte en arriere les radpas-1 pas de temps
2604	! suivant. Petite ruse avec MOD pour prendre en compte le
2605	! premier pas de temps de la physique pendant lequel
2606	! itaprad=0
2607	zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)
2608	zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
2609	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2610	latitude_deg,longitude_deg,rmu0,fract)
2611	!
2612	! Calcul des poids
2613	!
2614	zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
2615	zdtime2=0.0 !--pas de temps de la physique qui se termine
2616	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2617	latitude_deg,longitude_deg,zrmu0,zfract)
2618	swradcorr = 0.0
2619	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2620	swradcorr=zfract/fract*zrmu0/rmu0
2621	! Calcul du flag jour-nuit
2622	JrNt = 0.0
2623	WHERE (zfract.GT.0.0) JrNt = 1.0
2624	END SELECT
2625	ENDIF
2626	sza_o = ACOS (rmu0) *180./pi
2627
2628	IF (mydebug) THEN
2629	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2630	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2631	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2632	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2633	ENDIF
2634
2635	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2636	! Appel au pbl_surface : Planetary Boudary Layer et Surface
2637	! Cela implique tous les interactions des sous-surfaces et la
2638	! partie diffusion turbulent du couche limit.
2639	!
2640	! Certains varibales de sorties de pbl_surface sont utiliser que pour
2641	! ecriture des fihiers hist_XXXX.nc, ces sont :
2642	! qsol, zq2m, s_pblh, s_lcl,
2643	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2644	! s_therm, s_trmb1, s_trmb2, s_trmb3,
2645	! zu10m, zv10m, fder,
2646	! zxqsurf, delta_qsurf,
2647	! rh2m, zxfluxu, zxfluxv,
2648	! frugs, agesno, fsollw, fsolsw,
2649	! d_ts, fevap, fluxlat, t2m,
2650	! wfbils, wfbilo, fluxt, fluxu, fluxv,
2651	!
2652	! Certains ne sont pas utiliser du tout :
2653	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2654	!
2655
2656	! Calcul de l'humidite de saturation au niveau du sol
2657
2658
2659
2660	IF (iflag_pbl/=0) THEN
2661
2662	!jyg+nrlmd<
2663	!!jyg IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2664	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
2665	print *,'debut du splitting de la PBL, wake_s = ', wake_s(:)
2666	print *,'debut du splitting de la PBL, wake_deltat = ', wake_deltat(:,1)
2667	print *,'debut du splitting de la PBL, wake_deltaq = ', wake_deltaq(:,1)
2668	ENDIF
2669	! !!
2670	!>jyg+nrlmd
2671	!
2672	!-------gustiness calculation-------!
2673	!ym : Warning gustiness non inialized for iflag_gusts=2 & iflag_gusts=3
2674	gustiness=0 !ym missing init
2675
2676	IF (iflag_gusts==0) THEN
2677	gustiness(1:klon)=0
2678	ELSE IF (iflag_gusts==1) THEN
2679	gustiness(1:klon)=f_gust_blale_bl(1:klon)+f_gust_wkale_wake(1:klon)
2680	ELSE IF (iflag_gusts==2) THEN
2681	gustiness(1:klon)=f_gust_blale_bl_stat(1:klon)+f_gust_wkale_wake(1:klon)
2682	! ELSE IF (iflag_gusts==2) THEN
2683	! do i = 1, klon
2684	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
2685	! *ale_wake(i) !! need to make sigma_wk accessible here
2686	! enddo
2687	! ELSE IF (iflag_gusts==3) THEN
2688	! do i = 1, klon
2689	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
2690	! enddo
2691	ENDIF
2692
2693	CALL pbl_surface( &
2694	phys_tstep, date0, itap, days_elapsed+1, &
2695	debut, lafin, &
2696	longitude_deg, latitude_deg, rugoro, zrmu0, &
2697	sollwdown, cldt, &
2698	rain_fall, snow_fall, solsw, solswfdiff, sollw, &
2699	gustiness, &
2700	t_seri, q_seri, u_seri, v_seri, &
2701	!nrlmd+jyg<
2702	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2703	!>nrlmd+jyg
2704	pplay, paprs, pctsrf, &
2705	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2706	!albedo SB <<<
2707	cdragh, cdragm, u1, v1, &
2708	beta_aridity, &
2709	!albedo SB >>>
2710	! albsol1, albsol2, sens, evap, &
2711	albsol_dir, albsol_dif, sens, evap, &
2712	!albedo SB <<<
2713	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
2714	zxtsol, zxfluxlat, zt2m, qsat2m, zn2mout, &
2715	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2716	!nrlmd<
2717	!jyg<
2718	d_t_vdf_w, d_q_vdf_w, &
2719	d_t_vdf_x, d_q_vdf_x, &
2720	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2721	!>jyg
2722	delta_tsurf,wake_dens, &
2723	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2724	kh,kh_x,kh_w, &
2725	!>nrlmd
2726	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2727	slab_wfbils, &
2728	qsol, zq2m, s_pblh, s_lcl, &
2729	!jyg<
2730	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2731	!>jyg
2732	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2733	s_therm, s_trmb1, s_trmb2, s_trmb3, &
2734	zustar, zu10m, zv10m, fder, &
2735	zxqsurf, delta_qsurf, rh2m, zxfluxu, zxfluxv, &
2736	z0m, z0h, agesno, fsollw, fsolsw, &
2737	d_ts, fevap, fluxlat, t2m, &
2738	wfbils, wfbilo, wfevap, wfrain, wfsnow, &
2739	fluxt, fluxu, fluxv, &
2740	dsens, devap, zxsnow, &
2741	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
2742	!nrlmd+jyg<
2743	wake_delta_pbl_TKE, &
2744	!>nrlmd+jyg
2745	treedrg )
2746	!FC
2747	!
2748	! Add turbulent diffusion tendency to the wake difference variables
2749	!!jyg IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2750	IF (mod(iflag_pbl_split,10) .NE. 0) THEN
2751	!jyg<
2752	d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2753	d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2754	CALL add_wake_tend &
2755	(d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
2756	ELSE
2757	d_deltat_vdf(:,:) = 0.
2758	d_deltaq_vdf(:,:) = 0.
2759	!>jyg
2760	ENDIF
2761
2762	!---------------------------------------------------------------------
2763	! ajout des tendances de la diffusion turbulente
2764	IF (klon_glo==1) THEN
2765	CALL add_pbl_tend &
2766	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2767	'vdf',abortphy,flag_inhib_tend,itap)
2768	ELSE
2769	CALL add_phys_tend &
2770	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2771	'vdf',abortphy,flag_inhib_tend,itap,0)
2772	ENDIF
2773	CALL prt_enerbil('vdf',itap)
2774	!--------------------------------------------------------------------
2775
2776	IF (mydebug) THEN
2777	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2778	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2779	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2780	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2781	ENDIF
2782
2783	!albedo SB >>>
2784	albsol1=0.
2785	albsol2=0.
2786	falb1=0.
2787	falb2=0.
2788	SELECT CASE(nsw)
2789	CASE(2)
2790	albsol1=albsol_dir(:,1)
2791	albsol2=albsol_dir(:,2)
2792	falb1=falb_dir(:,1,:)
2793	falb2=falb_dir(:,2,:)
2794	CASE(4)
2795	albsol1=albsol_dir(:,1)
2796	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3) &
2797	+albsol_dir(:,4)*SFRWL(4)
2798	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2799	falb1=falb_dir(:,1,:)
2800	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3) &
2801	+falb_dir(:,4,:)*SFRWL(4)
2802	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2803	CASE(6)
2804	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2) &
2805	+albsol_dir(:,3)*SFRWL(3)
2806	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2807	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5) &
2808	+albsol_dir(:,6)*SFRWL(6)
2809	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2810	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2) &
2811	+falb_dir(:,3,:)*SFRWL(3)
2812	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2813	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5) &
2814	+falb_dir(:,6,:)*SFRWL(6)
2815	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2816	END SELECt
2817	!albedo SB <<<
2818
2819
2820	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2821	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2822
2823	ENDIF
2824	! =================================================================== c
2825	! Calcul de Qsat
2826
2827	DO k = 1, klev
2828	DO i = 1, klon
2829	zx_t = t_seri(i,k)
2830	IF (thermcep) THEN
2831	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2832	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2833	zx_qs = MIN(0.5,zx_qs)
2834	zcor = 1./(1.-retv*zx_qs)
2835	zx_qs = zx_qs*zcor
2836	ELSE
2837	!! IF (zx_t.LT.t_coup) THEN !jyg
2838	IF (zx_t.LT.rtt) THEN !jyg
2839	zx_qs = qsats(zx_t)/pplay(i,k)
2840	ELSE
2841	zx_qs = qsatl(zx_t)/pplay(i,k)
2842	ENDIF
2843	ENDIF
2844	zqsat(i,k)=zx_qs
2845	ENDDO
2846	ENDDO
2847
2848	IF (prt_level.ge.1) THEN
2849	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2850	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2851	ENDIF
2852	!
2853	! Appeler la convection (au choix)
2854	!
2855	DO k = 1, klev
2856	DO i = 1, klon
2857	conv_q(i,k) = d_q_dyn(i,k) &
2858	+ d_q_vdf(i,k)/phys_tstep
2859	conv_t(i,k) = d_t_dyn(i,k) &
2860	+ d_t_vdf(i,k)/phys_tstep
2861	ENDDO
2862	ENDDO
2863	IF (check) THEN
2864	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2865	WRITE(lunout,*) "avantcon=", za
2866	ENDIF
2867	zx_ajustq = .FALSE.
2868	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2869	IF (zx_ajustq) THEN
2870	DO i = 1, klon
2871	z_avant(i) = 0.0
2872	ENDDO
2873	DO k = 1, klev
2874	DO i = 1, klon
2875	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2876	*(paprs(i,k)-paprs(i,k+1))/RG
2877	ENDDO
2878	ENDDO
2879	ENDIF
2880
2881	! Calcule de vitesse verticale a partir de flux de masse verticale
2882	DO k = 1, klev
2883	DO i = 1, klon
2884	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
2885	ENDDO
2886	ENDDO
2887
2888	IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2889	omega(igout, :)
2890	!
2891	! Appel de la convection tous les "cvpas"
2892	!
2893	!!jyg IF (MOD(itapcv,cvpas).EQ.0) THEN
2894	!! print *,' physiq : itapcv, cvpas, itap-1, cvpas_0 ', &
2895	!! itapcv, cvpas, itap-1, cvpas_0
2896	IF (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itap-1,cvpas_0).EQ.0) THEN
2897
2898	!
2899	! Mettre a zero des variables de sortie (pour securite)
2900	!
2901	pmflxr(:,:) = 0.
2902	pmflxs(:,:) = 0.
2903	wdtrainA(:,:) = 0.
2904	wdtrainS(:,:) = 0.
2905	wdtrainM(:,:) = 0.
2906	upwd(:,:) = 0.
2907	dnwd(:,:) = 0.
2908	ep(:,:) = 0.
2909	da(:,:)=0.
2910	mp(:,:)=0.
2911	wght_cvfd(:,:)=0.
2912	phi(:,:,:)=0.
2913	phi2(:,:,:)=0.
2914	epmlmMm(:,:,:)=0.
2915	eplaMm(:,:)=0.
2916	d1a(:,:)=0.
2917	dam(:,:)=0.
2918	elij(:,:,:)=0.
2919	ev(:,:)=0.
2920	qtaa(:,:)=0.
2921	clw(:,:)=0.
2922	sij(:,:,:)=0.
2923	!
2924	IF (iflag_con.EQ.1) THEN
2925	abort_message ='reactiver le call conlmd dans physiq.F'
2926	CALL abort_physic (modname,abort_message,1)
2927	! CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
2928	! . d_t_con, d_q_con,
2929	! . rain_con, snow_con, ibas_con, itop_con)
2930	ELSE IF (iflag_con.EQ.2) THEN
2931	CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
2932	conv_t, conv_q, -evap, omega, &
2933	d_t_con, d_q_con, rain_con, snow_con, &
2934	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2935	kcbot, kctop, kdtop, pmflxr, pmflxs)
2936	d_u_con = 0.
2937	d_v_con = 0.
2938
2939	WHERE (rain_con < 0.) rain_con = 0.
2940	WHERE (snow_con < 0.) snow_con = 0.
2941	DO i = 1, klon
2942	ibas_con(i) = klev+1 - kcbot(i)
2943	itop_con(i) = klev+1 - kctop(i)
2944	ENDDO
2945	ELSE IF (iflag_con.GE.3) THEN
2946	! nb of tracers for the KE convection:
2947	! MAF la partie traceurs est faite dans phytrac
2948	! on met ntra=1 pour limiter les appels mais on peut
2949	! supprimer les calculs / ftra.
2950	ntra = 1
2951
2952	!=======================================================================
2953	!ajout pour la parametrisation des poches froides: calcul de
2954	!t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
2955	IF (iflag_wake>=1) THEN
2956	DO k=1,klev
2957	DO i=1,klon
2958	t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
2959	q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
2960	t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2961	q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2962	ENDDO
2963	ENDDO
2964	ELSE
2965	t_w(:,:) = t_seri(:,:)
2966	q_w(:,:) = q_seri(:,:)
2967	t_x(:,:) = t_seri(:,:)
2968	q_x(:,:) = q_seri(:,:)
2969	ENDIF
2970	!
2971	!jyg<
2972	! Perform dry adiabatic adjustment on wake profile
2973	! The corresponding tendencies are added to the convective tendencies
2974	! after the call to the convective scheme.
2975	IF (iflag_wake>=1) then
2976	IF (iflag_adjwk >= 1) THEN
2977	limbas(:) = 1
2978	CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
2979	d_t_adjwk, d_q_adjwk)
2980	!
2981	DO k=1,klev
2982	DO i=1,klon
2983	IF (wake_s(i) .GT. 1.e-3) THEN
2984	t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
2985	q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
2986	d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
2987	d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
2988	ELSE
2989	d_deltat_ajs_cv(i,k) = 0.
2990	d_deltaq_ajs_cv(i,k) = 0.
2991	ENDIF
2992	ENDDO
2993	ENDDO
2994	IF (iflag_adjwk == 2) THEN
2995	CALL add_wake_tend &
2996	(d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
2997	ENDIF ! (iflag_adjwk == 2)
2998	ENDIF ! (iflag_adjwk >= 1)
2999	ENDIF ! (iflag_wake>=1)
3000	!>jyg
3001	!
3002
3003	!! print *,'physiq. q_w(1,k), q_x(1,k) ', &
3004	!! (k, q_w(1,k), q_x(1,k),k=1,25)
3005
3006	!jyg<
3007	CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri, &
3008	alp_offset, it_wape_prescr, wape_prescr, fip_prescr, &
3009	ale_bl_prescr, alp_bl_prescr, &
3010	wake_pe, wake_fip, &
3011	Ale_bl, Ale_bl_trig, Alp_bl, &
3012	Ale, Alp , Ale_wake, Alp_wake)
3013	!>jyg
3014	!
3015	! sb, oct02:
3016	! Schema de convection modularise et vectorise:
3017	! (driver commun aux versions 3 et 4)
3018	!
3019	IF (ok_cvl) THEN ! new driver for convectL
3020	!
3021	!jyg<
3022	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3023	! Calculate the upmost level of deep convection loops: k_upper_cv
3024	! (near 22 km)
3025	k_upper_cv = klev
3026	!izero = klon/2+1/klon
3027	!DO k = klev,1,-1
3028	! IF (pphi(izero,k) > 22.e4) k_upper_cv = k
3029	!ENDDO
3030	! FH : nouveau calcul base sur un profil global sans quoi
3031	! le modele etait sensible au decoupage de domaines
3032	DO k = klev,1,-1
3033	IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
3034	ENDDO
3035	IF (prt_level .ge. 5) THEN
3036	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
3037	k_upper_cv
3038	ENDIF
3039	!
3040	!>jyg
3041	IF (type_trac == 'repr') THEN
3042	nbtr_tmp=ntra
3043	ELSE
3044	nbtr_tmp=nbtr
3045	ENDIF
3046	!jyg iflag_con est dans clesphys
3047	!c CALL concvl (iflag_con,iflag_clos,
3048	CALL concvl (iflag_clos, &
3049	phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
3050	t_w,q_w,wake_s, &
3051	u_seri,v_seri,tr_seri,nbtr_tmp, &
3052	ALE,ALP, &
3053	sig1,w01, &
3054	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
3055	rain_con, snow_con, ibas_con, itop_con, sigd, &
3056	ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
3057	Ma,mipsh,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
3058	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
3059	! RomP >>>
3060	!! . pmflxr,pmflxs,da,phi,mp,
3061	!! . ftd,fqd,lalim_conv,wght_th)
3062	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,qtaa,clw,elij, &
3063	ftd,fqd,lalim_conv,wght_th, &
3064	ev, ep,epmlmMm,eplaMm, &
3065	wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
3066	tau_cld_cv,coefw_cld_cv,epmax_diag)
3067
3068	! RomP <<<
3069
3070	!IM begin
3071	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
3072	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
3073	!IM end
3074	!IM cf. FH
3075	clwcon0=qcondc
3076	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
3077	!
3078	!jyg<
3079	! If convective tendencies are too large, then call convection
3080	! every time step
3081	cvpas = cvpas_0
3082	DO k=1,k_upper_cv
3083	DO i=1,klon
3084	IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
3085	d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
3086	dtcon_multistep_max = 3.
3087	dqcon_multistep_max = 0.02
3088	ENDIF
3089	ENDDO
3090	ENDDO
3091	!
3092	DO k=1,k_upper_cv
3093	DO i=1,klon
3094	!! IF (abs(d_t_con(i,k)) > 0.24 .OR. &
3095	!! abs(d_q_con(i,k)) > 2.e-2) THEN
3096	IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
3097	abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
3098	cvpas = 1
3099	!! print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
3100	!! i,k,d_t_con(i,k),d_q_con(i,k)
3101	ENDIF
3102	ENDDO
3103	ENDDO
3104	!!! Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
3105	!!! call bcast(cvpas)
3106	!!! ------------------------------------------------------------
3107	!>jyg
3108	!
3109	DO i = 1, klon
3110	IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
3111	ENDDO
3112	!
3113	!jyg<
3114	! Add the tendency due to the dry adjustment of the wake profile
3115	IF (iflag_wake>=1) THEN
3116	IF (iflag_adjwk == 2) THEN
3117	DO k=1,klev
3118	DO i=1,klon
3119	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
3120	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
3121	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
3122	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
3123	ENDDO
3124	ENDDO
3125	ENDIF ! (iflag_adjwk = 2)
3126	ENDIF ! (iflag_wake>=1)
3127	!>jyg
3128	!
3129	ELSE ! ok_cvl
3130
3131	! MAF conema3 ne contient pas les traceurs
3132	CALL conema3 (phys_tstep, &
3133	paprs,pplay,t_seri,q_seri, &
3134	u_seri,v_seri,tr_seri,ntra, &
3135	sig1,w01, &
3136	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
3137	rain_con, snow_con, ibas_con, itop_con, &
3138	upwd,dnwd,dnwd0,bas,top, &
3139	Ma,cape,tvp,rflag, &
3140	pbase &
3141	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
3142	,clwcon0)
3143
3144	ENDIF ! ok_cvl
3145
3146	!
3147	! Correction precip
3148	rain_con = rain_con * cvl_corr
3149	snow_con = snow_con * cvl_corr
3150	!
3151
3152	IF (.NOT. ok_gust) THEN
3153	do i = 1, klon
3154	wd(i)=0.0
3155	enddo
3156	ENDIF
3157
3158	! =================================================================== c
3159	! Calcul des proprietes des nuages convectifs
3160	!
3161
3162	! calcul des proprietes des nuages convectifs
3163	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
3164	IF (iflag_cld_cv == 0) THEN
3165	CALL clouds_gno &
3166	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
3167	ELSE
3168	CALL clouds_bigauss &
3169	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
3170	ENDIF
3171
3172
3173	! =================================================================== c
3174
3175	DO i = 1, klon
3176	itop_con(i) = min(max(itop_con(i),1),klev)
3177	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
3178	ENDDO
3179
3180	DO i = 1, klon
3181	! C Risi modif: pour éviter pb de dépassement d'indice dans les cas
3182	! où i n'est pas un point convectif et donc ibas_con(i)=0
3183	! c'est un pb indépendant des isotopes
3184	if (ibas_con(i) > 0) then
3185	ema_pcb(i) = paprs(i,ibas_con(i))
3186	else
3187	ema_pcb(i) = 0.0
3188	endif
3189	ENDDO
3190	DO i = 1, klon
3191	! L'idicage de itop_con peut cacher un pb potentiel
3192	! FH sous la dictee de JYG, CR
3193	ema_pct(i) = paprs(i,itop_con(i)+1)
3194
3195	IF (itop_con(i).gt.klev-3) THEN
3196	IF (prt_level >= 9) THEN
3197	write(lunout,*)'La convection monte trop haut '
3198	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
3199	ENDIF
3200	ENDIF
3201	ENDDO
3202	ELSE IF (iflag_con.eq.0) THEN
3203	write(lunout,*) 'On n appelle pas la convection'
3204	clwcon0=0.
3205	rnebcon0=0.
3206	d_t_con=0.
3207	d_q_con=0.
3208	d_u_con=0.
3209	d_v_con=0.
3210	rain_con=0.
3211	snow_con=0.
3212	bas=1
3213	top=1
3214	ELSE
3215	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
3216	CALL abort_physic("physiq", "", 1)
3217	ENDIF
3218
3219	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
3220	! . d_u_con, d_v_con)
3221
3222	!jyg Reinitialize proba_notrig and itapcv when convection has been called
3223	proba_notrig(:) = 1.
3224	itapcv = 0
3225	ENDIF ! (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
3226	!
3227	itapcv = itapcv+1
3228	!
3229	! Compter les steps ou cvpas=1
3230	IF (cvpas == 1) THEN
3231	Ncvpaseq1 = Ncvpaseq1+1
3232	ENDIF
3233	IF (mod(itap,1000) == 0) THEN
3234	print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
3235	ENDIF
3236
3237	!!!jyg Appel diagnostique a add_phys_tend pour tester la conservation de
3238	!!! l'energie dans les courants satures.
3239	!! d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
3240	!! d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
3241	!! dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
3242	!! CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat, &
3243	!! dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
3244	!! itap, 1)
3245	!! call prt_enerbil('convection_sat',itap)
3246	!!
3247	!!
3248	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
3249	'convection',abortphy,flag_inhib_tend,itap,0)
3250	CALL prt_enerbil('convection',itap)
3251
3252	!-------------------------------------------------------------------------
3253
3254	IF (mydebug) THEN
3255	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3256	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3257	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3258	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3259	ENDIF
3260
3261	IF (check) THEN
3262	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3263	WRITE(lunout,*)"aprescon=", za
3264	zx_t = 0.0
3265	za = 0.0
3266	DO i = 1, klon
3267	za = za + cell_area(i)/REAL(klon)
3268	zx_t = zx_t + (rain_con(i)+ &
3269	snow_con(i))*cell_area(i)/REAL(klon)
3270	ENDDO
3271	zx_t = zx_t/za*phys_tstep
3272	WRITE(lunout,*)"Precip=", zx_t
3273	ENDIF
3274	IF (zx_ajustq) THEN
3275	DO i = 1, klon
3276	z_apres(i) = 0.0
3277	ENDDO
3278	DO k = 1, klev
3279	DO i = 1, klon
3280	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
3281	*(paprs(i,k)-paprs(i,k+1))/RG
3282	ENDDO
3283	ENDDO
3284	DO i = 1, klon
3285	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*phys_tstep) &
3286	/z_apres(i)
3287	ENDDO
3288	DO k = 1, klev
3289	DO i = 1, klon
3290	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
3291	z_factor(i).LT.(1.0-1.0E-08)) THEN
3292	q_seri(i,k) = q_seri(i,k) * z_factor(i)
3293	ENDIF
3294	ENDDO
3295	ENDDO
3296	ENDIF
3297	zx_ajustq=.FALSE.
3298
3299	!
3300	!==========================================================================
3301	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
3302	!pour la couche limite diffuse pour l instant
3303	!
3304	!
3305	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
3306	! il faut rajouter cette tendance calcul\'ee hors des poches
3307	! froides
3308	!
3309	IF (iflag_wake>=1) THEN
3310	!
3311	!
3312	! Call wakes every "wkpas" step
3313	!
3314	IF (MOD(itapwk,wkpas).EQ.0) THEN
3315	!
3316	DO k=1,klev
3317	DO i=1,klon
3318	dt_dwn(i,k) = ftd(i,k)
3319	dq_dwn(i,k) = fqd(i,k)
3320	M_dwn(i,k) = dnwd0(i,k)
3321	M_up(i,k) = upwd(i,k)
3322	dt_a(i,k) = d_t_con(i,k)/phys_tstep - ftd(i,k)
3323	dq_a(i,k) = d_q_con(i,k)/phys_tstep - fqd(i,k)
3324	ENDDO
3325	ENDDO
3326
3327	IF (iflag_wake==2) THEN
3328	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3329	DO k = 1,klev
3330	dt_dwn(:,k)= dt_dwn(:,k)+ &
3331	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
3332	dq_dwn(:,k)= dq_dwn(:,k)+ &
3333	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
3334	ENDDO
3335	ELSEIF (iflag_wake==3) THEN
3336	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3337	DO k = 1,klev
3338	DO i=1,klon
3339	IF (rneb(i,k)==0.) THEN
3340	! On ne tient compte des tendances qu'en dehors des
3341	! nuages (c'est-\`a-dire a priri dans une region ou
3342	! l'eau se reevapore).
3343	dt_dwn(i,k)= dt_dwn(i,k)+ &
3344	ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
3345	dq_dwn(i,k)= dq_dwn(i,k)+ &
3346	ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
3347	ENDIF
3348	ENDDO
3349	ENDDO
3350	ENDIF
3351
3352	!
3353	!calcul caracteristiques de la poche froide
3354	CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
3355	t_seri, q_seri, omega, &
3356	dt_dwn, dq_dwn, M_dwn, M_up, &
3357	dt_a, dq_a, cv_gen, &
3358	sigd, cin, &
3359	wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens, &
3360	wake_dth, wake_h, &
3361	!! wake_pe, wake_fip, wake_gfl, &
3362	wake_pe, wake_fip_0, wake_gfl, & !! jyg
3363	d_t_wake, d_q_wake, &
3364	wake_k, t_x, q_x, &
3365	wake_omgbdth, wake_dp_omgb, &
3366	wake_dtKE, wake_dqKE, &
3367	wake_omg, wake_dp_deltomg, &
3368	wake_spread, wake_Cstar, d_deltat_wk_gw, &
3369	d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk)
3370	!
3371	!jyg Reinitialize itapwk when wakes have been called
3372	itapwk = 0
3373	ENDIF ! (MOD(itapwk,wkpas).EQ.0)
3374	!
3375	itapwk = itapwk+1
3376	!
3377	!-----------------------------------------------------------------------
3378	! ajout des tendances des poches froides
3379	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
3380	abortphy,flag_inhib_tend,itap,0)
3381	CALL prt_enerbil('wake',itap)
3382	!------------------------------------------------------------------------
3383
3384	! Increment Wake state variables
3385	IF (iflag_wake_tend .GT. 0.) THEN
3386
3387	CALL add_wake_tend &
3388	(d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, &
3389	'wake', abortphy)
3390	CALL prt_enerbil('wake',itap)
3391	ENDIF ! (iflag_wake_tend .GT. 0.)
3392	!
3393	IF (prt_level .GE. 10) THEN
3394	print *,' physiq, after calwake, wake_s: ',wake_s(:)
3395	print *,' physiq, after calwake, wake_deltat: ',wake_deltat(:,1)
3396	print *,' physiq, after calwake, wake_deltaq: ',wake_deltaq(:,1)
3397	ENDIF
3398
3399	IF (iflag_alp_wk_cond .GT. 0.) THEN
3400
3401	CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
3402	wake_fip)
3403	ELSE
3404	wake_fip(:) = wake_fip_0(:)
3405	ENDIF ! (iflag_alp_wk_cond .GT. 0.)
3406
3407	ENDIF ! (iflag_wake>=1)
3408	!
3409	!===================================================================
3410	! Convection seche (thermiques ou ajustement)
3411	!===================================================================
3412	!
3413	CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
3414	,seuil_inversion,weak_inversion,dthmin)
3415
3416
3417
3418	d_t_ajsb(:,:)=0.
3419	d_q_ajsb(:,:)=0.
3420	d_t_ajs(:,:)=0.
3421	d_u_ajs(:,:)=0.
3422	d_v_ajs(:,:)=0.
3423	d_q_ajs(:,:)=0.
3424	clwcon0th(:,:)=0.
3425	!
3426	! fm_therm(:,:)=0.
3427	! entr_therm(:,:)=0.
3428	! detr_therm(:,:)=0.
3429	!
3430	IF (prt_level>9) WRITE(lunout,*) &
3431	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
3432	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3433	IF (iflag_thermals<0) THEN
3434	! Rien
3435	! ====
3436	IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
3437
3438
3439	ELSE
3440
3441	! Thermiques
3442	! ==========
3443	IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
3444	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3445
3446
3447	!cc nrlmd le 10/04/2012
3448	DO k=1,klev+1
3449	DO i=1,klon
3450	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
3451	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
3452	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
3453	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
3454	ENDDO
3455	ENDDO
3456	!cc fin nrlmd le 10/04/2012
3457
3458	IF (iflag_thermals>=1) THEN
3459	!jyg<
3460	!! IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3461	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3462	! Appel des thermiques avec les profils exterieurs aux poches
3463	DO k=1,klev
3464	DO i=1,klon
3465	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3466	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
3467	u_therm(i,k) = u_seri(i,k)
3468	v_therm(i,k) = v_seri(i,k)
3469	ENDDO
3470	ENDDO
3471	ELSE
3472	! Appel des thermiques avec les profils moyens
3473	DO k=1,klev
3474	DO i=1,klon
3475	t_therm(i,k) = t_seri(i,k)
3476	q_therm(i,k) = q_seri(i,k)
3477	u_therm(i,k) = u_seri(i,k)
3478	v_therm(i,k) = v_seri(i,k)
3479	ENDDO
3480	ENDDO
3481	ENDIF
3482	!>jyg
3483	CALL calltherm(pdtphys &
3484	,pplay,paprs,pphi,weak_inversion &
3485	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
3486	,u_therm,v_therm,t_therm,q_therm,zqsat,debut & !jyg
3487	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
3488	,fm_therm,entr_therm,detr_therm &
3489	,zqasc,clwcon0th,lmax_th,ratqscth &
3490	,ratqsdiff,zqsatth &
3491	!on rajoute ale et alp, et les
3492	!caracteristiques de la couche alim
3493	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
3494	,ztv,zpspsk,ztla,zthl &
3495	!cc nrlmd le 10/04/2012
3496	,pbl_tke_input,pctsrf,omega,cell_area &
3497	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
3498	,n2,s2,ale_bl_stat &
3499	,therm_tke_max,env_tke_max &
3500	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
3501	,alp_bl_conv,alp_bl_stat &
3502	!cc fin nrlmd le 10/04/2012
3503	,zqla,ztva )
3504	!
3505	!jyg<
3506	!!jyg IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3507	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3508	! Si les thermiques ne sont presents que hors des
3509	! poches, la tendance moyenne associ\'ee doit etre
3510	! multipliee par la fraction surfacique qu'ils couvrent.
3511	DO k=1,klev
3512	DO i=1,klon
3513	!
3514	d_deltat_the(i,k) = - d_t_ajs(i,k)
3515	d_deltaq_the(i,k) = - d_q_ajs(i,k)
3516	!
3517	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
3518	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
3519	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
3520	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
3521	!
3522	ENDDO
3523	ENDDO
3524	!
3525	IF (ok_bug_split_th) THEN
3526	CALL add_wake_tend &
3527	(d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
3528	ELSE
3529	CALL add_wake_tend &
3530	(d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
3531	ENDIF
3532	CALL prt_enerbil('the',itap)
3533	!
3534	ENDIF ! (mod(iflag_pbl_split/10,10) .GE. 1)
3535	!
3536	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
3537	dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
3538	CALL prt_enerbil('thermals',itap)
3539	!
3540	!
3541	CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area, &
3542	cin, s2, n2, &
3543	ale_bl_trig, ale_bl_stat, ale_bl, &
3544	alp_bl, alp_bl_stat, &
3545	proba_notrig, random_notrig, cv_gen)
3546	!>jyg
3547
3548	! ------------------------------------------------------------------
3549	! Transport de la TKE par les panaches thermiques.
3550	! FH : 2010/02/01
3551	! if (iflag_pbl.eq.10) then
3552	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3553	! s rg,paprs,pbl_tke)
3554	! endif
3555	! -------------------------------------------------------------------
3556
3557	DO i=1,klon
3558	! zmax_th(i)=pphi(i,lmax_th(i))/rg
3559	!CR:04/05/12:correction calcul zmax
3560	zmax_th(i)=zmax0(i)
3561	ENDDO
3562
3563	ENDIF
3564
3565	! Ajustement sec
3566	! ==============
3567
3568	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3569	! a partir du sommet des thermiques.
3570	! Dans le cas contraire, on demarre au niveau 1.
3571
3572	IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
3573
3574	IF (iflag_thermals.eq.0) THEN
3575	IF (prt_level>9) WRITE(lunout,*)'ajsec'
3576	limbas(:)=1
3577	ELSE
3578	limbas(:)=lmax_th(:)
3579	ENDIF
3580
3581	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3582	! pour des test de convergence numerique.
3583	! Le nouveau ajsec est a priori mieux, meme pour le cas
3584	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3585	! non nulles numeriquement pour des mailles non concernees.
3586
3587	IF (iflag_thermals==0) THEN
3588	! Calling adjustment alone (but not the thermal plume model)
3589	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3590	, d_t_ajsb, d_q_ajsb)
3591	ELSE IF (iflag_thermals>0) THEN
3592	! Calling adjustment above the top of thermal plumes
3593	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3594	, d_t_ajsb, d_q_ajsb)
3595	ENDIF
3596
3597	!--------------------------------------------------------------------
3598	! ajout des tendances de l'ajustement sec ou des thermiques
3599	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
3600	'ajsb',abortphy,flag_inhib_tend,itap,0)
3601	CALL prt_enerbil('ajsb',itap)
3602	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3603	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3604
3605	!---------------------------------------------------------------------
3606
3607	ENDIF
3608
3609	ENDIF
3610	!
3611	!===================================================================
3612	! Computation of ratqs, the width (normalized) of the subrid scale
3613	! water distribution
3614
3615	tke_dissip_ave(:,:)=0.
3616	l_mix_ave(:,:)=0.
3617	wprime_ave(:,:)=0.
3618
3619	DO nsrf = 1, nbsrf
3620	DO i = 1, klon
3621	tke_dissip_ave(i,:) = tke_dissip_ave(i,:) + tke_dissip(i,:,nsrf)*pctsrf(i,nsrf)
3622	l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
3623	wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
3624	ENDDO
3625	ENDDO
3626
3627	CALL calcratqs(klon,klev,prt_level,lunout, &
3628	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3629	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
3630	tau_ratqs,fact_cldcon,wake_s, wake_deltaq, &
3631	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3632	paprs,pplay,t_seri,q_seri, qtc_cv, sigt_cv, zqsat, &
3633	pbl_tke(:,:,is_ave),tke_dissip_ave,l_mix_ave,wprime_ave,t2m,q2m,fm_therm, &
3634	ratqs,ratqsc,ratqs_inter)
3635
3636	!
3637	! Appeler le processus de condensation a grande echelle
3638	! et le processus de precipitation
3639	!-------------------------------------------------------------------------
3640	IF (prt_level .GE.10) THEN
3641	print *,'itap, ->fisrtilp ',itap
3642	ENDIF
3643	!
3644
3645	picefra(:,:)=0.
3646
3647	IF (ok_new_lscp) THEN
3648
3649	!--mise à jour de flight_m et flight_h2o dans leur module
3650	IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
3651	CALL airplane(debut,pphis,pplay,paprs,t_seri)
3652	ENDIF
3653
3654	CALL lscp(phys_tstep,missing_val,paprs,pplay, &
3655	t_seri, q_seri,ptconv,ratqs, &
3656	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneb_seri, &
3657	cldliq, picefra, rain_lsc, snow_lsc, &
3658	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3659	frac_impa, frac_nucl, beta_prec_fisrt, &
3660	prfl, psfl, rhcl, &
3661	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3662	iflag_ice_thermo, ok_ice_sursat)
3663
3664	ELSE
3665
3666	CALL fisrtilp(phys_tstep,paprs,pplay, &
3667	t_seri, q_seri,ptconv,ratqs, &
3668	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3669	rain_lsc, snow_lsc, &
3670	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3671	frac_impa, frac_nucl, beta_prec_fisrt, &
3672	prfl, psfl, rhcl, &
3673	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3674	iflag_ice_thermo)
3675
3676	ENDIF
3677	!
3678	WHERE (rain_lsc < 0) rain_lsc = 0.
3679	WHERE (snow_lsc < 0) snow_lsc = 0.
3680
3681	!+JLD
3682	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3683	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
3684	! & ,rain_lsc,snow_lsc
3685	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3686	!-JLD
3687	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
3688	'lsc',abortphy,flag_inhib_tend,itap,0)
3689	CALL prt_enerbil('lsc',itap)
3690	rain_num(:)=0.
3691	DO k = 1, klev
3692	DO i = 1, klon
3693	IF (ql_seri(i,k)>oliqmax) THEN
3694	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3695	ql_seri(i,k)=oliqmax
3696	ENDIF
3697	ENDDO
3698	ENDDO
3699	IF (nqo >= 3) THEN
3700	DO k = 1, klev
3701	DO i = 1, klon
3702	IF (qs_seri(i,k)>oicemax) THEN
3703	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3704	qs_seri(i,k)=oicemax
3705	ENDIF
3706	ENDDO
3707	ENDDO
3708	ENDIF
3709
3710	!---------------------------------------------------------------------------
3711	DO k = 1, klev
3712	DO i = 1, klon
3713	cldfra(i,k) = rneb(i,k)
3714	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3715	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3716	ENDDO
3717	ENDDO
3718	IF (check) THEN
3719	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3720	WRITE(lunout,*)"apresilp=", za
3721	zx_t = 0.0
3722	za = 0.0
3723	DO i = 1, klon
3724	za = za + cell_area(i)/REAL(klon)
3725	zx_t = zx_t + (rain_lsc(i) &
3726	+ snow_lsc(i))*cell_area(i)/REAL(klon)
3727	ENDDO
3728	zx_t = zx_t/za*phys_tstep
3729	WRITE(lunout,*)"Precip=", zx_t
3730	ENDIF
3731
3732	IF (mydebug) THEN
3733	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3734	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3735	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3736	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3737	ENDIF
3738
3739	!
3740	!-------------------------------------------------------------------
3741	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3742	!-------------------------------------------------------------------
3743
3744	! 1. NUAGES CONVECTIFS
3745	!
3746	!IM cf FH
3747	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3748	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3749	snow_tiedtke=0.
3750	! print*,'avant calcul de la pseudo precip '
3751	! print*,'iflag_cld_th',iflag_cld_th
3752	IF (iflag_cld_th.eq.-1) THEN
3753	rain_tiedtke=rain_con
3754	ELSE
3755	! print*,'calcul de la pseudo precip '
3756	rain_tiedtke=0.
3757	! print*,'calcul de la pseudo precip 0'
3758	DO k=1,klev
3759	DO i=1,klon
3760	IF (d_q_con(i,k).lt.0.) THEN
3761	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3762	*(paprs(i,k)-paprs(i,k+1))/rg
3763	ENDIF
3764	ENDDO
3765	ENDDO
3766	ENDIF
3767	!
3768	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3769	!
3770
3771	! Nuages diagnostiques pour Tiedtke
3772	CALL diagcld1(paprs,pplay, &
3773	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3774	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3775	diafra,dialiq)
3776	DO k = 1, klev
3777	DO i = 1, klon
3778	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3779	cldliq(i,k) = dialiq(i,k)
3780	cldfra(i,k) = diafra(i,k)
3781	ENDIF
3782	ENDDO
3783	ENDDO
3784
3785	ELSE IF (iflag_cld_th.ge.3) THEN
3786	! On prend pour les nuages convectifs le max du calcul de la
3787	! convection et du calcul du pas de temps precedent diminue d'un facteur
3788	! facttemps
3789	facteur = pdtphys *facttemps
3790	DO k=1,klev
3791	DO i=1,klon
3792	rnebcon(i,k)=rnebcon(i,k)*facteur
3793	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
3794	rnebcon(i,k)=rnebcon0(i,k)
3795	clwcon(i,k)=clwcon0(i,k)
3796	ENDIF
3797	ENDDO
3798	ENDDO
3799
3800	! On prend la somme des fractions nuageuses et des contenus en eau
3801
3802	IF (iflag_cld_th>=5) THEN
3803
3804	DO k=1,klev
3805	ptconvth(:,k)=fm_therm(:,k+1)>0.
3806	ENDDO
3807
3808	IF (iflag_coupl==4) THEN
3809
3810	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3811	! convectives et lsc dans la partie des thermiques
3812	! Le controle par iflag_coupl est peut etre provisoire.
3813	DO k=1,klev
3814	DO i=1,klon
3815	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
3816	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3817	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3818	ELSE IF (ptconv(i,k)) THEN
3819	cldfra(i,k)=rnebcon(i,k)
3820	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3821	ENDIF
3822	ENDDO
3823	ENDDO
3824
3825	ELSE IF (iflag_coupl==5) THEN
3826	DO k=1,klev
3827	DO i=1,klon
3828	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3829	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3830	ENDDO
3831	ENDDO
3832
3833	ELSE
3834
3835	! Si on est sur un point touche par la convection
3836	! profonde et pas par les thermiques, on prend la
3837	! couverture nuageuse et l'eau nuageuse de la convection
3838	! profonde.
3839
3840	!IM/FH: 2011/02/23
3841	! definition des points sur lesquels ls thermiques sont actifs
3842
3843	DO k=1,klev
3844	DO i=1,klon
3845	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
3846	cldfra(i,k)=rnebcon(i,k)
3847	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3848	ENDIF
3849	ENDDO
3850	ENDDO
3851
3852	ENDIF
3853
3854	ELSE
3855
3856	! Ancienne version
3857	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3858	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3859	ENDIF
3860
3861	ENDIF
3862
3863	! plulsc(:)=0.
3864	! do k=1,klev,-1
3865	! do i=1,klon
3866	! zzz=prfl(:,k)+psfl(:,k)
3867	! if (.not.ptconvth.zzz.gt.0.)
3868	! enddo prfl, psfl,
3869	! enddo
3870	!
3871	! 2. NUAGES STARTIFORMES
3872	!
3873	IF (ok_stratus) THEN
3874	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3875	DO k = 1, klev
3876	DO i = 1, klon
3877	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3878	cldliq(i,k) = dialiq(i,k)
3879	cldfra(i,k) = diafra(i,k)
3880	ENDIF
3881	ENDDO
3882	ENDDO
3883	ENDIF
3884	!
3885	! Precipitation totale
3886	!
3887	DO i = 1, klon
3888	rain_fall(i) = rain_con(i) + rain_lsc(i)
3889	snow_fall(i) = snow_con(i) + snow_lsc(i)
3890	ENDDO
3891	!
3892	! Calculer l'humidite relative pour diagnostique
3893	!
3894	DO k = 1, klev
3895	DO i = 1, klon
3896	zx_t = t_seri(i,k)
3897	IF (thermcep) THEN
3898	!! if (iflag_ice_thermo.eq.0) then !jyg
3899	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3900	!! else !jyg
3901	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
3902	!! endif !jyg
3903	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3904	zx_qs = MIN(0.5,zx_qs)
3905	zcor = 1./(1.-retv*zx_qs)
3906	zx_qs = zx_qs*zcor
3907	ELSE
3908	!! IF (zx_t.LT.t_coup) THEN !jyg
3909	IF (zx_t.LT.rtt) THEN !jyg
3910	zx_qs = qsats(zx_t)/pplay(i,k)
3911	ELSE
3912	zx_qs = qsatl(zx_t)/pplay(i,k)
3913	ENDIF
3914	ENDIF
3915	zx_rh(i,k) = q_seri(i,k)/zx_qs
3916	IF (iflag_ice_thermo .GT. 0) THEN
3917	zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
3918	zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
3919	ENDIF
3920	zqsat(i,k)=zx_qs
3921	ENDDO
3922	ENDDO
3923
3924	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3925	! equivalente a 2m (tpote) pour diagnostique
3926	!
3927	DO i = 1, klon
3928	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3929	IF (thermcep) THEN
3930	IF(zt2m(i).LT.RTT) then
3931	Lheat=RLSTT
3932	ELSE
3933	Lheat=RLVTT
3934	ENDIF
3935	ELSE
3936	IF (zt2m(i).LT.RTT) THEN
3937	Lheat=RLSTT
3938	ELSE
3939	Lheat=RLVTT
3940	ENDIF
3941	ENDIF
3942	tpote(i) = tpot(i)* &
3943	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3944	ENDDO
3945
3946	IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN ! ModThL
3947	#ifdef INCA
3948	CALL VTe(VTphysiq)
3949	CALL VTb(VTinca)
3950	calday = REAL(days_elapsed + 1) + jH_cur
3951
3952	CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
3953	CALL AEROSOL_METEO_CALC( &
3954	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3955	prfl,psfl,pctsrf,cell_area, &
3956	latitude_deg,longitude_deg,u10m,v10m)
3957
3958	zxsnow_dummy(:) = 0.0
3959
3960	CALL chemhook_begin (calday, &
3961	days_elapsed+1, &
3962	jH_cur, &
3963	pctsrf(1,1), &
3964	latitude_deg, &
3965	longitude_deg, &
3966	cell_area, &
3967	paprs, &
3968	pplay, &
3969	coefh(1:klon,1:klev,is_ave), &
3970	pphi, &
3971	t_seri, &
3972	u, &
3973	v, &
3974	rot, &
3975	wo(:, :, 1), &
3976	q_seri, &
3977	zxtsol, &
3978	zt2m, &
3979	zxsnow_dummy, &
3980	solsw, &
3981	albsol1, &
3982	rain_fall, &
3983	snow_fall, &
3984	itop_con, &
3985	ibas_con, &
3986	cldfra, &
3987	nbp_lon, &
3988	nbp_lat-1, &
3989	tr_seri(:,:,1+nqCO2:nbtr), &
3990	ftsol, &
3991	paprs, &
3992	cdragh, &
3993	cdragm, &
3994	pctsrf, &
3995	pdtphys, &
3996	itap)
3997
3998	CALL VTe(VTinca)
3999	CALL VTb(VTphysiq)
4000	#endif
4001	ENDIF !type_trac = inca or inco
4002	IF (type_trac == 'repr') THEN
4003	#ifdef REPROBUS
4004	!CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap)
4005	CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap)
4006	#endif
4007	ENDIF
4008
4009	!
4010	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
4011	!
4012	IF (MOD(itaprad,radpas).EQ.0) THEN
4013
4014	!
4015	!jq - introduce the aerosol direct and first indirect radiative forcings
4016	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
4017	IF (flag_aerosol .GT. 0) THEN
4018	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4019	IF (.NOT. aerosol_couple) THEN
4020	!
4021	CALL readaerosol_optic( &
4022	debut, flag_aerosol, itap, jD_cur-jD_ref, &
4023	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4024	mass_solu_aero, mass_solu_aero_pi, &
4025	tau_aero, piz_aero, cg_aero, &
4026	tausum_aero, tau3d_aero)
4027	ENDIF
4028	ELSE ! RRTM radiation
4029	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
4030	abort_message='config_inca=aero et rrtm=1 impossible'
4031	CALL abort_physic(modname,abort_message,1)
4032	ELSE
4033	!
4034	#ifdef CPP_RRTM
4035	IF (NSW.EQ.6) THEN
4036	!--new aerosol properties SW and LW
4037	!
4038	#ifdef CPP_Dust
4039	!--SPL aerosol model
4040	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
4041	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
4042	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4043	tausum_aero, tau3d_aero)
4044	#else
4045	!--climatologies or INCA aerosols
4046	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
4047	flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
4048	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4049	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
4050	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4051	tausum_aero, drytausum_aero, tau3d_aero)
4052	#endif
4053
4054	IF (flag_aerosol .EQ. 7) THEN
4055	CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg, &
4056	tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
4057	ENDIF
4058
4059	!
4060	ELSE IF (NSW.EQ.2) THEN
4061	!--for now we use the old aerosol properties
4062	!
4063	CALL readaerosol_optic( &
4064	debut, flag_aerosol, itap, jD_cur-jD_ref, &
4065	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4066	mass_solu_aero, mass_solu_aero_pi, &
4067	tau_aero, piz_aero, cg_aero, &
4068	tausum_aero, tau3d_aero)
4069	!
4070	!--natural aerosols
4071	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
4072	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
4073	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
4074	!--all aerosols
4075	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
4076	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
4077	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
4078	!
4079	!--no LW optics
4080	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4081	!
4082	ELSE
4083	abort_message='Only NSW=2 or 6 are possible with ' &
4084	// 'aerosols and iflag_rrtm=1'
4085	CALL abort_physic(modname,abort_message,1)
4086	ENDIF
4087	#else
4088	abort_message='You should compile with -rrtm if running ' &
4089	// 'with iflag_rrtm=1'
4090	CALL abort_physic(modname,abort_message,1)
4091	#endif
4092	!
4093	ENDIF
4094	ENDIF
4095	ELSE !--flag_aerosol = 0
4096	tausum_aero(:,:,:) = 0.
4097	drytausum_aero(:,:) = 0.
4098	mass_solu_aero(:,:) = 0.
4099	mass_solu_aero_pi(:,:) = 0.
4100	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4101	tau_aero(:,:,:,:) = 1.e-15
4102	piz_aero(:,:,:,:) = 1.
4103	cg_aero(:,:,:,:) = 0.
4104	ELSE
4105	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
4106	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4107	piz_aero_sw_rrtm(:,:,:,:) = 1.0
4108	cg_aero_sw_rrtm(:,:,:,:) = 0.0
4109	ENDIF
4110	ENDIF
4111	!
4112	!--WMO criterion to determine tropopause
4113	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
4114	!
4115	!--STRAT AEROSOL
4116	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
4117	IF (flag_aerosol_strat.GT.0) THEN
4118	IF (prt_level .GE.10) THEN
4119	PRINT *,'appel a readaerosolstrat', mth_cur
4120	ENDIF
4121	IF (iflag_rrtm.EQ.0) THEN
4122	IF (flag_aerosol_strat.EQ.1) THEN
4123	CALL readaerosolstrato(debut)
4124	ELSE
4125	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
4126	CALL abort_physic(modname,abort_message,1)
4127	ENDIF
4128	ELSE
4129	#ifdef CPP_RRTM
4130	#ifndef CPP_StratAer
4131	!--prescribed strat aerosols
4132	!--only in the case of non-interactive strat aerosols
4133	IF (flag_aerosol_strat.EQ.1) THEN
4134	CALL readaerosolstrato1_rrtm(debut)
4135	ELSEIF (flag_aerosol_strat.EQ.2) THEN
4136	CALL readaerosolstrato2_rrtm(debut, ok_volcan)
4137	ELSE
4138	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
4139	CALL abort_physic(modname,abort_message,1)
4140	ENDIF
4141	#endif
4142	#else
4143	abort_message='You should compile with -rrtm if running ' &
4144	// 'with iflag_rrtm=1'
4145	CALL abort_physic(modname,abort_message,1)
4146	#endif
4147	ENDIF
4148	ELSE
4149	tausum_aero(:,:,id_STRAT_phy) = 0.
4150	ENDIF
4151	!
4152	#ifdef CPP_RRTM
4153	#ifdef CPP_StratAer
4154	!--compute stratospheric mask
4155	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
4156	!--interactive strat aerosols
4157	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
4158	#endif
4159	#endif
4160	!--fin STRAT AEROSOL
4161	!
4162
4163	! Calculer les parametres optiques des nuages et quelques
4164	! parametres pour diagnostiques:
4165	!
4166	IF (aerosol_couple.AND.config_inca=='aero') THEN
4167	mass_solu_aero(:,:) = ccm(:,:,1)
4168	mass_solu_aero_pi(:,:) = ccm(:,:,2)
4169	ENDIF
4170
4171	IF (ok_newmicro) then
4172	! AI IF (iflag_rrtm.NE.0) THEN
4173	IF (iflag_rrtm.EQ.1) THEN
4174	#ifdef CPP_RRTM
4175	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
4176	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
4177	// 'pour ok_cdnc'
4178	CALL abort_physic(modname,abort_message,1)
4179	ENDIF
4180	#else
4181
4182	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
4183	CALL abort_physic(modname,abort_message,1)
4184	#endif
4185	ENDIF
4186	CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
4187	paprs, pplay, t_seri, cldliq, picefra, cldfra, &
4188	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
4189	flwp, fiwp, flwc, fiwc, &
4190	mass_solu_aero, mass_solu_aero_pi, &
4191	cldtaupi, latitude_deg, re, fl, ref_liq, ref_ice, &
4192	ref_liq_pi, ref_ice_pi)
4193	ELSE
4194	CALL nuage (paprs, pplay, &
4195	t_seri, cldliq, picefra, cldfra, cldtau, cldemi, &
4196	cldh, cldl, cldm, cldt, cldq, &
4197	ok_aie, &
4198	mass_solu_aero, mass_solu_aero_pi, &
4199	bl95_b0, bl95_b1, &
4200	cldtaupi, re, fl)
4201	ENDIF
4202	!
4203	!IM betaCRF
4204	!
4205	cldtaurad = cldtau
4206	cldtaupirad = cldtaupi
4207	cldemirad = cldemi
4208	cldfrarad = cldfra
4209
4210	!
4211	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
4212	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
4213	!
4214	! global
4215	!
4216	!IM 251017 begin
4217	! print*,'physiq betaCRF global zdtime=',zdtime
4218	!IM 251017 end
4219	DO k=1, klev
4220	DO i=1, klon
4221	IF (pplay(i,k).GE.pfree) THEN
4222	beta(i,k) = beta_pbl
4223	ELSE
4224	beta(i,k) = beta_free
4225	ENDIF
4226	IF (mskocean_beta) THEN
4227	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4228	ENDIF
4229	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4230	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4231	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4232	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4233	ENDDO
4234	ENDDO
4235	!
4236	ELSE
4237	!
4238	! regional
4239	!
4240	DO k=1, klev
4241	DO i=1,klon
4242	!
4243	IF (longitude_deg(i).ge.lon1_beta.AND. &
4244	longitude_deg(i).le.lon2_beta.AND. &
4245	latitude_deg(i).le.lat1_beta.AND. &
4246	latitude_deg(i).ge.lat2_beta) THEN
4247	IF (pplay(i,k).GE.pfree) THEN
4248	beta(i,k) = beta_pbl
4249	ELSE
4250	beta(i,k) = beta_free
4251	ENDIF
4252	IF (mskocean_beta) THEN
4253	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4254	ENDIF
4255	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4256	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4257	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4258	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4259	ENDIF
4260	!
4261	ENDDO
4262	ENDDO
4263	!
4264	ENDIF
4265
4266	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
4267	IF (ok_chlorophyll) THEN
4268	print*,"-- reading chlorophyll"
4269	CALL readchlorophyll(debut)
4270	ENDIF
4271
4272	!--if ok_suntime_rrtm we use ancillay data for RSUN
4273	!--previous values are therefore overwritten
4274	!--this is needed for CMIP6 runs
4275	!--and only possible for new radiation scheme
4276	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
4277	#ifdef CPP_RRTM
4278	CALL read_rsun_rrtm(debut)
4279	#endif
4280	ENDIF
4281
4282	IF (mydebug) THEN
4283	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4284	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4285	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4286	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4287	ENDIF
4288
4289	!
4290	!sonia : If Iflag_radia >=2, pertubation of some variables
4291	!input to radiation (DICE)
4292	!
4293	IF (iflag_radia .ge. 2) THEN
4294	zsav_tsol (:) = zxtsol(:)
4295	CALL perturb_radlwsw(zxtsol,iflag_radia)
4296	ENDIF
4297
4298	IF (aerosol_couple.AND.config_inca=='aero') THEN
4299	#ifdef INCA
4300	CALL radlwsw_inca &
4301	(chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
4302	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
4303	size(wo,3), wo, &
4304	cldfrarad, cldemirad, cldtaurad, &
4305	heat,heat0,cool,cool0,albpla, &
4306	topsw,toplw,solsw,sollw, &
4307	sollwdown, &
4308	topsw0,toplw0,solsw0,sollw0, &
4309	lwdn0, lwdn, lwup0, lwup, &
4310	swdn0, swdn, swup0, swup, &
4311	ok_ade, ok_aie, &
4312	tau_aero, piz_aero, cg_aero, &
4313	topswad_aero, solswad_aero, &
4314	topswad0_aero, solswad0_aero, &
4315	topsw_aero, topsw0_aero, &
4316	solsw_aero, solsw0_aero, &
4317	cldtaupirad, &
4318	topswai_aero, solswai_aero)
4319	#endif
4320	ELSE
4321	!
4322	!IM calcul radiatif pour le cas actuel
4323	!
4324	RCO2 = RCO2_act
4325	RCH4 = RCH4_act
4326	RN2O = RN2O_act
4327	RCFC11 = RCFC11_act
4328	RCFC12 = RCFC12_act
4329	!
4330	!--interactive CO2 in ppm from carbon cycle
4331	IF (carbon_cycle_rad.AND..NOT.debut) THEN
4332	RCO2=RCO2_glo
4333	ENDIF
4334	!
4335	IF (prt_level .GE.10) THEN
4336	print *,' ->radlwsw, number 1 '
4337	ENDIF
4338	!
4339	CALL radlwsw &
4340	(dist, rmu0, fract, &
4341	!albedo SB >>>
4342	! paprs, pplay,zxtsol,albsol1, albsol2, &
4343	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4344	!albedo SB <<<
4345	t_seri,q_seri,wo, &
4346	cldfrarad, cldemirad, cldtaurad, &
4347	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4348	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4349	tau_aero, piz_aero, cg_aero, &
4350	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4351	! Rajoute par OB pour RRTM
4352	tau_aero_lw_rrtm, &
4353	cldtaupirad, &
4354	! zqsat, flwcrad, fiwcrad, &
4355	zqsat, flwc, fiwc, &
4356	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4357	heat,heat0,cool,cool0,albpla, &
4358	heat_volc,cool_volc, &
4359	topsw,toplw,solsw,solswfdiff,sollw, &
4360	sollwdown, &
4361	topsw0,toplw0,solsw0,sollw0, &
4362	lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, &
4363	swdnc0, swdn0, swdn, swupc0, swup0, swup, &
4364	topswad_aero, solswad_aero, &
4365	topswai_aero, solswai_aero, &
4366	topswad0_aero, solswad0_aero, &
4367	topsw_aero, topsw0_aero, &
4368	solsw_aero, solsw0_aero, &
4369	topswcf_aero, solswcf_aero, &
4370	!-C. Kleinschmitt for LW diagnostics
4371	toplwad_aero, sollwad_aero,&
4372	toplwai_aero, sollwai_aero, &
4373	toplwad0_aero, sollwad0_aero,&
4374	!-end
4375	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4376	ZSWFT0_i, ZFSDN0, ZFSUP0)
4377
4378	!lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
4379	!schemes
4380	toplw = toplw + betalwoff * (toplw0 - toplw)
4381	sollw = sollw + betalwoff * (sollw0 - sollw)
4382	lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
4383	lwup = lwup + betalwoff * (lwup0 - lwup)
4384	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
4385	sollwdown(:))
4386	cool = cool + betalwoff * (cool0 - cool)
4387
4388	#ifndef CPP_XIOS
4389	!--OB 30/05/2016 modified 21/10/2016
4390	!--here we return swaero_diag and dryaod_diag to FALSE
4391	!--and histdef will switch it back to TRUE if necessary
4392	!--this is necessary to get the right swaero at first step
4393	!--but only in the case of no XIOS as XIOS is covered elsewhere
4394	IF (debut) swaerofree_diag = .FALSE.
4395	IF (debut) swaero_diag = .FALSE.
4396	IF (debut) dryaod_diag = .FALSE.
4397	!--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
4398	!--as for swaero_diag, see above
4399	IF (debut) ok_4xCO2atm = .FALSE.
4400
4401	!
4402	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
4403	!IM des taux doit etre different du taux actuel
4404	!IM Par defaut on a les taux perturbes egaux aux taux actuels
4405	!
4406	IF (RCO2_per.NE.RCO2_act.OR. &
4407	RCH4_per.NE.RCH4_act.OR. &
4408	RN2O_per.NE.RN2O_act.OR. &
4409	RCFC11_per.NE.RCFC11_act.OR. &
4410	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
4411	#endif
4412	!
4413	IF (ok_4xCO2atm) THEN
4414	!
4415	RCO2 = RCO2_per
4416	RCH4 = RCH4_per
4417	RN2O = RN2O_per
4418	RCFC11 = RCFC11_per
4419	RCFC12 = RCFC12_per
4420	!
4421	IF (prt_level .GE.10) THEN
4422	print *,' ->radlwsw, number 2 '
4423	ENDIF
4424	!
4425	CALL radlwsw &
4426	(dist, rmu0, fract, &
4427	!albedo SB >>>
4428	! paprs, pplay,zxtsol,albsol1, albsol2, &
4429	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4430	!albedo SB <<<
4431	t_seri,q_seri,wo, &
4432	cldfrarad, cldemirad, cldtaurad, &
4433	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4434	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4435	tau_aero, piz_aero, cg_aero, &
4436	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4437	! Rajoute par OB pour RRTM
4438	tau_aero_lw_rrtm, &
4439	cldtaupi, &
4440	! zqsat, flwcrad, fiwcrad, &
4441	zqsat, flwc, fiwc, &
4442	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4443	heatp,heat0p,coolp,cool0p,albplap, &
4444	heat_volc,cool_volc, &
4445	topswp,toplwp,solswp,solswfdiffp,sollwp, &
4446	sollwdownp, &
4447	topsw0p,toplw0p,solsw0p,sollw0p, &
4448	lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp, &
4449	swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
4450	topswad_aerop, solswad_aerop, &
4451	topswai_aerop, solswai_aerop, &
4452	topswad0_aerop, solswad0_aerop, &
4453	topsw_aerop, topsw0_aerop, &
4454	solsw_aerop, solsw0_aerop, &
4455	topswcf_aerop, solswcf_aerop, &
4456	!-C. Kleinschmitt for LW diagnostics
4457	toplwad_aerop, sollwad_aerop,&
4458	toplwai_aerop, sollwai_aerop, &
4459	toplwad0_aerop, sollwad0_aerop,&
4460	!-end
4461	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4462	ZSWFT0_i, ZFSDN0, ZFSUP0)
4463	ENDIF !ok_4xCO2atm
4464	ENDIF ! aerosol_couple
4465	itaprad = 0
4466	!
4467	! If Iflag_radia >=2, reset pertubed variables
4468	!
4469	IF (iflag_radia .ge. 2) THEN
4470	zxtsol(:) = zsav_tsol (:)
4471	ENDIF
4472	ENDIF ! MOD(itaprad,radpas)
4473	itaprad = itaprad + 1
4474
4475	IF (iflag_radia.eq.0) THEN
4476	IF (prt_level.ge.9) THEN
4477	PRINT *,'--------------------------------------------------'
4478	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4479	PRINT *,'>>>> heat et cool mis a zero '
4480	PRINT *,'--------------------------------------------------'
4481	ENDIF
4482	heat=0.
4483	cool=0.
4484	sollw=0. ! MPL 01032011
4485	solsw=0.
4486	radsol=0.
4487	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4488	swup0=0.
4489	lwup=0.
4490	lwup0=0.
4491	lwdn=0.
4492	lwdn0=0.
4493	ENDIF
4494
4495	!
4496	! Calculer radsol a l'exterieur de radlwsw
4497	! pour prendre en compte le cycle diurne
4498	! recode par Olivier Boucher en sept 2015
4499	!
4500	radsol=solsw*swradcorr+sollw
4501
4502	IF (ok_4xCO2atm) THEN
4503	radsolp=solswp*swradcorr+sollwp
4504	ENDIF
4505
4506	!
4507	! Ajouter la tendance des rayonnements (tous les pas)
4508	! avec une correction pour le cycle diurne dans le SW
4509	!
4510
4511	DO k=1, klev
4512	d_t_swr(:,k)=swradcorr(:)heat(:,k)phys_tstep/RDAY
4513	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)phys_tstep/RDAY
4514	d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
4515	d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
4516	ENDDO
4517
4518	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4519	CALL prt_enerbil('SW',itap)
4520	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4521	CALL prt_enerbil('LW',itap)
4522
4523	!
4524	IF (mydebug) THEN
4525	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4526	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4527	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4528	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4529	ENDIF
4530
4531	! Calculer l'hydrologie de la surface
4532	!
4533	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4534	! . agesno, ftsol,fqsurf,fsnow, ruis)
4535	!
4536
4537	!
4538	! Calculer le bilan du sol et la derive de temperature (couplage)
4539	!
4540	DO i = 1, klon
4541	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4542	! a la demande de JLD
4543	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4544	ENDDO
4545	!
4546	!moddeblott(jan95)
4547	! Appeler le programme de parametrisation de l'orographie
4548	! a l'echelle sous-maille:
4549	!
4550	IF (prt_level .GE.10) THEN
4551	print *,' call orography ? ', ok_orodr
4552	ENDIF
4553	!
4554	IF (ok_orodr) THEN
4555	!
4556	! selection des points pour lesquels le shema est actif:
4557	igwd=0
4558	DO i=1,klon
4559	itest(i)=0
4560	! IF ((zstd(i).gt.10.0)) THEN
4561	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4562	itest(i)=1
4563	igwd=igwd+1
4564	idx(igwd)=i
4565	ENDIF
4566	ENDDO
4567	! igwdim=MAX(1,igwd)
4568	!
4569	IF (ok_strato) THEN
4570
4571	CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
4572	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4573	igwd,idx,itest, &
4574	t_seri, u_seri, v_seri, &
4575	zulow, zvlow, zustrdr, zvstrdr, &
4576	d_t_oro, d_u_oro, d_v_oro)
4577
4578	ELSE
4579	CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
4580	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4581	igwd,idx,itest, &
4582	t_seri, u_seri, v_seri, &
4583	zulow, zvlow, zustrdr, zvstrdr, &
4584	d_t_oro, d_u_oro, d_v_oro)
4585	ENDIF
4586	!
4587	! ajout des tendances
4588	!-----------------------------------------------------------------------
4589	! ajout des tendances de la trainee de l'orographie
4590	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
4591	abortphy,flag_inhib_tend,itap,0)
4592	CALL prt_enerbil('oro',itap)
4593	!----------------------------------------------------------------------
4594	!
4595	ENDIF ! fin de test sur ok_orodr
4596	!
4597	IF (mydebug) THEN
4598	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4599	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4600	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4601	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4602	ENDIF
4603
4604	IF (ok_orolf) THEN
4605	!
4606	! selection des points pour lesquels le shema est actif:
4607	igwd=0
4608	DO i=1,klon
4609	itest(i)=0
4610	IF ((zpic(i)-zmea(i)).GT.100.) THEN
4611	itest(i)=1
4612	igwd=igwd+1
4613	idx(igwd)=i
4614	ENDIF
4615	ENDDO
4616	! igwdim=MAX(1,igwd)
4617	!
4618	IF (ok_strato) THEN
4619
4620	CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
4621	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4622	igwd,idx,itest, &
4623	t_seri, u_seri, v_seri, &
4624	zulow, zvlow, zustrli, zvstrli, &
4625	d_t_lif, d_u_lif, d_v_lif )
4626
4627	ELSE
4628	CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
4629	latitude_deg,zmea,zstd,zpic, &
4630	itest, &
4631	t_seri, u_seri, v_seri, &
4632	zulow, zvlow, zustrli, zvstrli, &
4633	d_t_lif, d_u_lif, d_v_lif)
4634	ENDIF
4635
4636	! ajout des tendances de la portance de l'orographie
4637	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
4638	'lif', abortphy,flag_inhib_tend,itap,0)
4639	CALL prt_enerbil('lif',itap)
4640	ENDIF ! fin de test sur ok_orolf
4641
4642	IF (ok_hines) then
4643	! HINES GWD PARAMETRIZATION
4644	east_gwstress=0.
4645	west_gwstress=0.
4646	du_gwd_hines=0.
4647	dv_gwd_hines=0.
4648	CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
4649	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4650	du_gwd_hines, dv_gwd_hines)
4651	zustr_gwd_hines=0.
4652	zvstr_gwd_hines=0.
4653	DO k = 1, klev
4654	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
4655	* (paprs(:, k)-paprs(:, k+1))/rg
4656	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
4657	* (paprs(:, k)-paprs(:, k+1))/rg
4658	ENDDO
4659
4660	d_t_hin(:, :)=0.
4661	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4662	dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4663	CALL prt_enerbil('hin',itap)
4664	ENDIF
4665
4666	IF (.not. ok_hines .and. ok_gwd_rando) then
4667	! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
4668	CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
4669	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4670	dv_gwd_front, east_gwstress, west_gwstress)
4671	zustr_gwd_front=0.
4672	zvstr_gwd_front=0.
4673	DO k = 1, klev
4674	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
4675	* (paprs(:, k)-paprs(:, k+1))/rg
4676	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
4677	* (paprs(:, k)-paprs(:, k+1))/rg
4678	ENDDO
4679
4680	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
4681	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4682	CALL prt_enerbil('front_gwd_rando',itap)
4683	ENDIF
4684
4685	IF (ok_gwd_rando) THEN
4686	CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
4687	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4688	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4689	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
4690	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4691	CALL prt_enerbil('flott_gwd_rando',itap)
4692	zustr_gwd_rando=0.
4693	zvstr_gwd_rando=0.
4694	DO k = 1, klev
4695	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
4696	* (paprs(:, k)-paprs(:, k+1))/rg
4697	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
4698	* (paprs(:, k)-paprs(:, k+1))/rg
4699	ENDDO
4700	ENDIF
4701
4702	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4703
4704	IF (mydebug) THEN
4705	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4706	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4707	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4708	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4709	ENDIF
4710
4711	DO i = 1, klon
4712	zustrph(i)=0.
4713	zvstrph(i)=0.
4714	ENDDO
4715	DO k = 1, klev
4716	DO i = 1, klon
4717	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
4718	(paprs(i,k)-paprs(i,k+1))/rg
4719	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
4720	(paprs(i,k)-paprs(i,k+1))/rg
4721	ENDDO
4722	ENDDO
4723	!
4724	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4725	!
4726	IF (is_sequential .and. ok_orodr) THEN
4727	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4728	ra,rg,romega, &
4729	latitude_deg,longitude_deg,pphis, &
4730	zustrdr,zustrli,zustrph, &
4731	zvstrdr,zvstrli,zvstrph, &
4732	paprs,u,v, &
4733	aam, torsfc)
4734	ENDIF
4735	!IM cf. FLott END
4736	!DC Calcul de la tendance due au methane
4737	IF (ok_qch4) THEN
4738	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4739	! ajout de la tendance d'humidite due au methane
4740	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
4741	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
4742	'q_ch4', abortphy,flag_inhib_tend,itap,0)
4743	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
4744	ENDIF
4745	!
4746	!
4747
4748	!===============================================================
4749	! Additional tendency of TKE due to orography
4750	!===============================================================
4751	!
4752	! Inititialization
4753	!------------------
4754
4755	addtkeoro=0
4756	CALL getin_p('addtkeoro',addtkeoro)
4757
4758	IF (prt_level.ge.5) &
4759	print*,'addtkeoro', addtkeoro
4760
4761	alphatkeoro=1.
4762	CALL getin_p('alphatkeoro',alphatkeoro)
4763	alphatkeoro=min(max(0.,alphatkeoro),1.)
4764
4765	smallscales_tkeoro=.FALSE.
4766	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4767
4768
4769	dtadd(:,:)=0.
4770	duadd(:,:)=0.
4771	dvadd(:,:)=0.
4772
4773	! Choices for addtkeoro:
4774	! ** 0 no TKE tendency from orography
4775	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4776	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4777	!
4778
4779	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4780	! -------------------------------------------
4781
4782
4783	! selection des points pour lesquels le schema est actif:
4784
4785
4786	IF (addtkeoro .EQ. 1 ) THEN
4787
4788	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4789	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4790
4791	ELSE IF (addtkeoro .EQ. 2) THEN
4792
4793	IF (smallscales_tkeoro) THEN
4794	igwd=0
4795	DO i=1,klon
4796	itest(i)=0
4797	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4798	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4799	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4800	IF (zstd(i).GT.1.0) THEN
4801	itest(i)=1
4802	igwd=igwd+1
4803	idx(igwd)=i
4804	ENDIF
4805	ENDDO
4806
4807	ELSE
4808
4809	igwd=0
4810	DO i=1,klon
4811	itest(i)=0
4812	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4813	itest(i)=1
4814	igwd=igwd+1
4815	idx(igwd)=i
4816	ENDIF
4817	ENDDO
4818
4819	ENDIF
4820
4821	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
4822	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4823	igwd,idx,itest, &
4824	t_seri, u_seri, v_seri, &
4825	zulow, zvlow, zustrdr, zvstrdr, &
4826	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4827
4828	zustrdr(:)=0.
4829	zvstrdr(:)=0.
4830	zulow(:)=0.
4831	zvlow(:)=0.
4832
4833	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4834	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4835	ENDIF
4836
4837
4838	! TKE update from subgrid temperature and wind tendencies
4839	!----------------------------------------------------------
4840	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4841
4842
4843	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
4844	!
4845	! Prevent pbl_tke_w from becoming negative
4846	wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
4847	!
4848
4849	ENDIF
4850	! -----
4851	!===============================================================
4852
4853
4854	!====================================================================
4855	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4856	!====================================================================
4857	! Abderrahmane 24.08.09
4858
4859	IF (ok_cosp) THEN
4860	! adeclarer
4861	#ifdef CPP_COSP
4862	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4863
4864	IF (prt_level .GE.10) THEN
4865	print*,'freq_cosp',freq_cosp
4866	ENDIF
4867	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4868	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4869	! s ref_liq,ref_ice
4870	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
4871	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4872	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4873	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4874	JrNt,ref_liq,ref_ice, &
4875	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4876	zu10m,zv10m,pphis, &
4877	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4878	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4879	prfl(:,1:klev),psfl(:,1:klev), &
4880	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4881	mr_ozone,cldtau, cldemi)
4882
4883	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4884	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4885	! M clMISR,
4886	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4887	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4888
4889	ENDIF
4890	#endif
4891
4892	#ifdef CPP_COSP2
4893	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4894
4895	IF (prt_level .GE.10) THEN
4896	print*,'freq_cosp',freq_cosp
4897	ENDIF
4898	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4899	print*,'Dans physiq.F avant appel '
4900	! s ref_liq,ref_ice
4901	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
4902	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4903	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4904	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4905	JrNt,ref_liq,ref_ice, &
4906	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4907	zu10m,zv10m,pphis, &
4908	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4909	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4910	prfl(:,1:klev),psfl(:,1:klev), &
4911	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4912	mr_ozone,cldtau, cldemi)
4913	ENDIF
4914	#endif
4915
4916	#ifdef CPP_COSPV2
4917	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4918	! IF (MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4919
4920	IF (prt_level .GE.10) THEN
4921	print*,'freq_cosp',freq_cosp
4922	ENDIF
4923	DO k = 1, klev
4924	DO i = 1, klon
4925	phicosp(i,k) = pphi(i,k) + pphis(i)
4926	ENDDO
4927	ENDDO
4928	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4929	print*,'Dans physiq.F avant appel '
4930	! s ref_liq,ref_ice
4931	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
4932	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4933	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4934	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4935	JrNt,ref_liq,ref_ice, &
4936	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4937	zu10m,zv10m,pphis, &
4938	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4939	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4940	prfl(:,1:klev),psfl(:,1:klev), &
4941	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4942	mr_ozone,cldtau, cldemi)
4943	ENDIF
4944	#endif
4945
4946	ENDIF !ok_cosp
4947
4948
4949	! Marine
4950
4951	IF (ok_airs) then
4952
4953	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
4954	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4955	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4956	& map_prop_hc,map_prop_hist,&
4957	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4958	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4959	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4960	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4961	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4962	& map_ntot,map_hc,map_hist,&
4963	& map_Cb,map_ThCi,map_Anv,&
4964	& alt_tropo )
4965	ENDIF
4966
4967	ENDIF ! ok_airs
4968
4969
4970	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4971	!AA
4972	!AA Installation de l'interface online-offline pour traceurs
4973	!AA
4974	!====================================================================
4975	! Calcul des tendances traceurs
4976	!====================================================================
4977	!
4978
4979	IF (type_trac=='repr') THEN
4980	!MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
4981	!MM dans Reprobus
4982	sh_in(:,:) = q_seri(:,:)
4983	#ifdef REPROBUS
4984	d_q_rep(:,:) = 0.
4985	d_ql_rep(:,:) = 0.
4986	d_qi_rep(:,:) = 0.
4987	#endif
4988	ELSE
4989	sh_in(:,:) = qx(:,:,ivap)
4990	IF (nqo >= 3) THEN
4991	ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
4992	ELSE
4993	ch_in(:,:) = qx(:,:,iliq)
4994	ENDIF
4995	ENDIF
4996
4997	#ifdef CPP_Dust
4998	! Avec SPLA, iflag_phytrac est forcé =1
4999	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
5000	pdtphys,ftsol, & ! I
5001	t,q_seri,paprs,pplay,RHcl, & ! I
5002	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
5003	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
5004	u_seri, v_seri, latitude_deg, longitude_deg, &
5005	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
5006	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
5007	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
5008	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
5009	fm_therm, entr_therm, rneb, & ! I
5010	beta_prec_fisrt,beta_prec, & !I
5011	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
5012	d_tr_dyn,tr_seri)
5013
5014	#else
5015	IF (iflag_phytrac == 1 ) THEN
5016	CALL phytrac ( &
5017	itap, days_elapsed+1, jH_cur, debut, &
5018	lafin, phys_tstep, u, v, t, &
5019	paprs, pplay, pmfu, pmfd, &
5020	pen_u, pde_u, pen_d, pde_d, &
5021	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
5022	u1, v1, ftsol, pctsrf, &
5023	zustar, zu10m, zv10m, &
5024	wstar(:,is_ave), ale_bl, ale_wake, &
5025	latitude_deg, longitude_deg, &
5026	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
5027	presnivs, pphis, pphi, albsol1, &
5028	sh_in, ch_in, rhcl, cldfra, rneb, &
5029	diafra, cldliq, itop_con, ibas_con, &
5030	pmflxr, pmflxs, prfl, psfl, &
5031	da, phi, mp, upwd, &
5032	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
5033	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
5034	ev, ep, epmlmMm, eplaMm, & !<<RomP
5035	dnwd, aerosol_couple, flxmass_w, &
5036	tau_aero, piz_aero, cg_aero, ccm, &
5037	rfname, &
5038	d_tr_dyn, & !<<RomP
5039	tr_seri, init_source)
5040	#ifdef REPROBUS
5041
5042
5043	print*,'avt add phys rep',abortphy
5044
5045	CALL add_phys_tend &
5046	(du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,paprs,&
5047	'rep',abortphy,flag_inhib_tend,itap,0)
5048	IF (abortphy==1) Print*,'ERROR ABORT REP'
5049
5050	print*,'apr add phys rep',abortphy
5051
5052	#endif
5053	ENDIF ! (iflag_phytrac=1)
5054
5055	#endif
5056	!ENDIF ! (iflag_phytrac=1)
5057
5058	IF (offline) THEN
5059
5060	IF (prt_level.ge.9) &
5061	print*,'Attention on met a 0 les thermiques pour phystoke'
5062	CALL phystokenc ( &
5063	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
5064	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
5065	fm_therm,entr_therm, &
5066	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
5067	frac_impa, frac_nucl, &
5068	pphis,cell_area,phys_tstep,itap, &
5069	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
5070
5071
5072	ENDIF
5073
5074	!
5075	! Calculer le transport de l'eau et de l'energie (diagnostique)
5076	!
5077	CALL transp (paprs,zxtsol, &
5078	t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
5079	ve, vq, ue, uq, vwat, uwat)
5080	!
5081	!IM global posePB BEG
5082	IF(1.EQ.0) THEN
5083	!
5084	CALL transp_lay (paprs,zxtsol, &
5085	t_seri, q_seri, u_seri, v_seri, zphi, &
5086	ve_lay, vq_lay, ue_lay, uq_lay)
5087	!
5088	ENDIF !(1.EQ.0) THEN
5089	!IM global posePB END
5090	! Accumuler les variables a stocker dans les fichiers histoire:
5091	!
5092
5093	!================================================================
5094	! Conversion of kinetic and potential energy into heat, for
5095	! parameterisation of subgrid-scale motions
5096	!================================================================
5097
5098	d_t_ec(:,:)=0.
5099	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
5100	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
5101	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
5102	zmasse,exner,d_t_ec)
5103	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
5104
5105	!=======================================================================
5106	! SORTIES
5107	!=======================================================================
5108	!
5109	!IM initialisation + calculs divers diag AMIP2
5110	!
5111	include "calcul_divers.h"
5112	!
5113	!IM Interpolation sur les niveaux de pression du NMC
5114	! -------------------------------------------------
5115	!
5116	include "calcul_STDlev.h"
5117	!
5118	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
5119	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
5120	!
5121	!cc prw = eau precipitable
5122	! prlw = colonne eau liquide
5123	! prlw = colonne eau solide
5124	prw(:) = 0.
5125	prlw(:) = 0.
5126	prsw(:) = 0.
5127	DO k = 1, klev
5128	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
5129	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
5130	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
5131	ENDDO
5132	!
5133	IF (type_trac == 'inca' .OR. type_trac == 'inco') THEN
5134	#ifdef INCA
5135	CALL VTe(VTphysiq)
5136	CALL VTb(VTinca)
5137
5138	CALL chemhook_end ( &
5139	phys_tstep, &
5140	pplay, &
5141	t_seri, &
5142	tr_seri(:,:,1+nqCO2:nbtr), &
5143	nbtr, &
5144	paprs, &
5145	q_seri, &
5146	cell_area, &
5147	pphi, &
5148	pphis, &
5149	zx_rh, &
5150	aps, bps, ap, bp, lafin)
5151
5152	CALL VTe(VTinca)
5153	CALL VTb(VTphysiq)
5154	#endif
5155	ENDIF
5156
5157
5158	!
5159	! Convertir les incrementations en tendances
5160	!
5161	IF (prt_level .GE.10) THEN
5162	print *,'Convertir les incrementations en tendances '
5163	ENDIF
5164	!
5165	IF (mydebug) THEN
5166	CALL writefield_phy('u_seri',u_seri,nbp_lev)
5167	CALL writefield_phy('v_seri',v_seri,nbp_lev)
5168	CALL writefield_phy('t_seri',t_seri,nbp_lev)
5169	CALL writefield_phy('q_seri',q_seri,nbp_lev)
5170	ENDIF
5171
5172	DO k = 1, klev
5173	DO i = 1, klon
5174	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
5175	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
5176	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
5177	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
5178	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
5179	!CR: on ajoute le contenu en glace
5180	IF (nqo >= 3) THEN
5181	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
5182	ENDIF
5183	!--ice_sursat: nqo=4, on ajoute rneb
5184	IF (nqo == 4) THEN
5185	d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
5186	ENDIF
5187	ENDDO
5188	ENDDO
5189	!
5190	IF (nqtot > nqo) THEN
5191	itr = 0
5192	DO iq = 1, nqtot
5193	IF(.NOT.tracers(iq)%isInPhysics) CYCLE
5194	itr = itr+1
5195	DO k = 1, klev
5196	DO i = 1, klon
5197	d_qx(i,k,iq) = ( tr_seri(i,k,itr) - qx(i,k,iq) ) / phys_tstep
5198	ENDDO
5199	ENDDO
5200	ENDDO
5201	ENDIF
5202	!
5203	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
5204	!IM global posePB include "write_bilKP_ins.h"
5205	!IM global posePB include "write_bilKP_ave.h"
5206	!
5207
5208	!--OB mass fixer
5209	!--profile is corrected to force mass conservation of water
5210	IF (mass_fixer) THEN
5211	qql2(:)=0.0
5212	DO k = 1, klev
5213	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
5214	ENDDO
5215	DO i = 1, klon
5216	!--compute ratio of what q+ql should be with conservation to what it is
5217	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
5218	DO k = 1, klev
5219	q_seri(i,k) =q_seri(i,k)*corrqql
5220	ql_seri(i,k)=ql_seri(i,k)*corrqql
5221	ENDDO
5222	ENDDO
5223	ENDIF
5224	!--fin mass fixer
5225
5226	! Sauvegarder les valeurs de t et q a la fin de la physique:
5227	!
5228	u_ancien(:,:) = u_seri(:,:)
5229	v_ancien(:,:) = v_seri(:,:)
5230	t_ancien(:,:) = t_seri(:,:)
5231	q_ancien(:,:) = q_seri(:,:)
5232	ql_ancien(:,:) = ql_seri(:,:)
5233	qs_ancien(:,:) = qs_seri(:,:)
5234	rneb_ancien(:,:) = rneb_seri(:,:)
5235	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
5236	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
5237	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
5238	! !! RomP >>>
5239	IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
5240	! !! RomP <<<
5241	!==========================================================================
5242	! Sorties des tendances pour un point particulier
5243	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
5244	! pour le debug
5245	! La valeur de igout est attribuee plus haut dans le programme
5246	!==========================================================================
5247
5248	IF (prt_level.ge.1) THEN
5249	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
5250	write(lunout,*) &
5251	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
5252	write(lunout,*) &
5253	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
5254	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
5255	pctsrf(igout,is_sic)
5256	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
5257	DO k=1,klev
5258	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
5259	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
5260	d_t_eva(igout,k)
5261	ENDDO
5262	write(lunout,*) 'cool,heat'
5263	DO k=1,klev
5264	write(lunout,*) cool(igout,k),heat(igout,k)
5265	ENDDO
5266
5267	!jyg< (En attendant de statuer sur le sort de d_t_oli)
5268	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5269	!jyg! do k=1,klev
5270	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
5271	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5272	!jyg! enddo
5273	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5274	DO k=1,klev
5275	write(lunout,*) d_t_vdf(igout,k), &
5276	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5277	ENDDO
5278	!>jyg
5279
5280	write(lunout,*) 'd_ps ',d_ps(igout)
5281	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
5282	DO k=1,klev
5283	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
5284	d_qx(igout,k,1),d_qx(igout,k,2)
5285	ENDDO
5286	ENDIF
5287
5288	!============================================================
5289	! Calcul de la temperature potentielle
5290	!============================================================
5291	DO k = 1, klev
5292	DO i = 1, klon
5293	!JYG/IM theta en debut du pas de temps
5294	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
5295	!JYG/IM theta en fin de pas de temps de physique
5296	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
5297	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
5298	! MPL 20130625
5299	! fth_fonctions.F90 et parkind1.F90
5300	! sinon thetal=theta
5301	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
5302	! : ql_seri(i,k))
5303	thetal(i,k)=theta(i,k)
5304	ENDDO
5305	ENDDO
5306	!
5307
5308	! 22.03.04 BEG
5309	!=============================================================
5310	! Ecriture des sorties
5311	!=============================================================
5312	#ifdef CPP_IOIPSL
5313
5314	! Recupere des varibles calcule dans differents modules
5315	! pour ecriture dans histxxx.nc
5316
5317	! Get some variables from module fonte_neige_mod
5318	CALL fonte_neige_get_vars(pctsrf, &
5319	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
5320
5321
5322	!=============================================================
5323	! Separation entre thermiques et non thermiques dans les sorties
5324	! de fisrtilp
5325	!=============================================================
5326
5327	IF (iflag_thermals>=1) THEN
5328	d_t_lscth=0.
5329	d_t_lscst=0.
5330	d_q_lscth=0.
5331	d_q_lscst=0.
5332	DO k=1,klev
5333	DO i=1,klon
5334	IF (ptconvth(i,k)) THEN
5335	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5336	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5337	ELSE
5338	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5339	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5340	ENDIF
5341	ENDDO
5342	ENDDO
5343
5344	DO i=1,klon
5345	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
5346	plul_th(i)=prfl(i,1)+psfl(i,1)
5347	ENDDO
5348	ENDIF
5349
5350	!On effectue les sorties:
5351
5352	#ifdef CPP_Dust
5353	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
5354	pplay, lmax_th, aerosol_couple, &
5355	ok_ade, ok_aie, ivap, ok_sync, &
5356	ptconv, read_climoz, clevSTD, &
5357	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
5358	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5359	#else
5360	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
5361	pplay, lmax_th, aerosol_couple, &
5362	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, &
5363	ok_sync, ptconv, read_climoz, clevSTD, &
5364	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
5365	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5366	#endif
5367
5368	#ifndef CPP_XIOS
5369	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
5370	#endif
5371
5372	#endif
5373
5374	!====================================================================
5375	! Arret du modele apres hgardfou en cas de detection d'un
5376	! plantage par hgardfou
5377	!====================================================================
5378
5379	IF (abortphy==1) THEN
5380	abort_message ='Plantage hgardfou'
5381	CALL abort_physic (modname,abort_message,1)
5382	ENDIF
5383
5384	! 22.03.04 END
5385	!
5386	!====================================================================
5387	! Si c'est la fin, il faut conserver l'etat de redemarrage
5388	!====================================================================
5389	!
5390
5391	! Disabling calls to the prt_alerte function
5392	alert_first_call = .FALSE.
5393
5394	IF (lafin) THEN
5395	itau_phy = itau_phy + itap
5396	CALL phyredem ("restartphy.nc")
5397	! open(97,form="unformatted",file="finbin")
5398	! write(97) u_seri,v_seri,t_seri,q_seri
5399	! close(97)
5400
5401	IF (is_omp_master) THEN
5402
5403	IF (read_climoz >= 1) THEN
5404	IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5405	DEALLOCATE(press_edg_climoz) ! pointer
5406	DEALLOCATE(press_cen_climoz) ! pointer
5407	ENDIF
5408
5409	ENDIF
5410	#ifdef CPP_XIOS
5411	IF (is_omp_master) CALL xios_context_finalize
5412
5413	#ifdef INCA
5414	if (type_trac == 'inca' ) then
5415	IF (is_omp_master .and. grid_type==unstructured) THEN
5416	CALL finalize_inca
5417	ENDIF
5418	endif
5419	#endif
5420
5421	#endif
5422	WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
5423	ENDIF
5424
5425	! first=.false.
5426
5427	END SUBROUTINE physiq
5428
5429	END MODULE physiq_mod

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

Download in other formats:

Original Format