Context Navigation

physiq_mod.F90 @ 5087

Last change on this file since 5087 was 5087, checked in by abarral, 4 months ago
remove fixed-form ^{\s+& remaining in .f90,.F90}
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 212.6 KB

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

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

Download in other formats:

Original Format