Context Navigation

physiq_mod.F90 @ 4445

Last change on this file since 4445 was 4442, checked in by Laurent Fairhead, 17 months ago
GPU modifications to calling routines
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: 197.2 KB

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

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

Download in other formats:

Original Format