Context Navigation

physiq_mod.F90 @ 4141

Last change on this file since 4141 was 4141, checked in by Laurent Fairhead, 2 years ago
Outillage de physiq_mod.F90 en directives openacc pour prendre en compte le portage de ajsec.F90
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 196.3 KB

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

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

Download in other formats:

Original Format