Context Navigation

physiq_mod.F90 @ 5134

Last change on this file since 5134 was 5131, checked in by acozic, 5 months ago
update to write mass flow files in physiq - work begin with rev[4608]
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: 213.1 KB

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

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

Download in other formats:

Original Format