Context Navigation

physiq_mod.F90 @ 4150

Last change on this file since 4150 was 4150, checked in by Ehouarn Millour, 2 years ago
Fix previously introduced bug in the code reorganization for OpenACC. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 196.2 KB

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

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

Download in other formats:

Original Format