Context Navigation

physiq_mod.F90 @ 4872

Last change on this file since 4872 was 4872, checked in by idelkadi, 8 weeks ago
Initialization of input fields for the 1st call in LMDZ to the COSP simulator (version 2)
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.0 KB

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

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

Download in other formats:

Original Format