Context Navigation

physiq_mod.F90 @ 4881

Last change on this file since 4881 was 4881, checked in by evignon, 7 weeks ago
extraction plus propre de la dissipation de TKE
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 212.9 KB

Line
1	!
2	! $Id: physiq_mod.F90 4881 2024-03-28 17:06:40Z evignon $
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, pbl_eps, l_mix, wprime,&
285	t_therm, q_therm, u_therm, v_therm, &
286	cdragh_x, cdragh_w, &
287	cdragm_x, cdragm_w, &
288	kh, kh_x, kh_w, &
289	!
290	wake_k, &
291	alp_wake, &
292	wake_h, wake_omg, &
293	! tendencies of delta T and delta q:
294	d_deltat_wk, d_deltaq_wk, & ! due to wakes
295	d_deltat_wk_gw, d_deltaq_wk_gw, & ! due to wake induced gravity waves
296	d_deltat_vdf, d_deltaq_vdf, & ! due to vertical diffusion
297	d_deltat_the, d_deltaq_the, & ! due to thermals
298	d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
299	! tendencies of wake fractional area and wake number per unit area:
300	d_s_wk, d_s_a_wk, d_dens_wk, d_dens_a_wk, & ! due to wakes
301	!!! d_s_vdf, d_dens_a_vdf, d_dens_vdf, & ! due to vertical diffusion
302	!!! d_s_the, d_dens_a_the, d_dens_the, & ! due to thermals
303	!
304	ptconv, ratqsc, &
305	wbeff, convoccur, zmax_th, &
306	sens, flwp, fiwp, &
307	alp_bl_conv,alp_bl_det, &
308	alp_bl_fluct_m,alp_bl_fluct_tke, &
309	alp_bl_stat, n2, s2, strig, zcong, zlcl_th, &
310	proba_notrig, random_notrig, &
311	!! cv_gen, & !moved to phys_state_var_mod
312	!
313	dnwd0, &
314	omega, &
315	epmax_diag, &
316	! Deep convective variables used in phytrac
317	pmflxr, pmflxs, &
318	wdtrainA, wdtrainS, wdtrainM, &
319	upwd, dnwd, &
320	ep, &
321	da, mp, &
322	phi, &
323	wght_cvfd, &
324	phi2, &
325	d1a, dam, &
326	ev, &
327	elij, &
328	qtaa, &
329	clw, &
330	epmlmMm, eplaMm, &
331	sij, &
332	!
333	rneblsvol, &
334	pfraclr,pfracld, &
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) :: 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, pbl_eps, &
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	l_mix_ave(:,:)=0.
3852	wprime_ave(:,:)=0.
3853
3854	DO nsrf = 1, nbsrf
3855	DO i = 1, klon
3856	l_mix_ave(i,:) = l_mix_ave(i,:) + l_mix(i,:,nsrf)*pctsrf(i,nsrf)
3857	wprime_ave(i,:) = wprime_ave(i,:) + wprime(i,:,nsrf)*pctsrf(i,nsrf)
3858	ENDDO
3859	ENDDO
3860
3861	CALL ratqs_main(klon,klev,nbsrf,prt_level,lunout, &
3862	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3863	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
3864	pctsrf,s_pblh,zstd, &
3865	tau_ratqs,fact_cldcon,wake_s, wake_deltaq, &
3866	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3867	paprs,pplay,t_seri,q_seri, &
3868	qtc_cv, sigt_cv,detrain_cv,fm_cv,fqd,fqcomp,sigd,zqsat, &
3869	omega,pbl_tke(:,:,is_ave),pbl_eps(:,:,is_ave),l_mix_ave,wprime_ave, &
3870	t2m,q2m,fm_therm,entr_therm,detr_therm,cell_area, &
3871	ratqs,ratqsc,ratqs_inter_)
3872
3873	!
3874	! Appeler le processus de condensation a grande echelle
3875	! et le processus de precipitation
3876	!-------------------------------------------------------------------------
3877	IF (prt_level .GE.10) THEN
3878	print *,'itap, ->fisrtilp ',itap
3879	ENDIF
3880	!
3881
3882	picefra(:,:)=0.
3883
3884	IF (ok_new_lscp) THEN
3885
3886	!--mise à jour de flight_m et flight_h2o dans leur module
3887	IF (ok_plane_h2o .OR. ok_plane_contrail) THEN
3888	CALL airplane(debut,pphis,pplay,paprs,t_seri)
3889	ENDIF
3890
3891	CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, &
3892	t_seri, q_seri,ptconv,ratqs, &
3893	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, &
3894	pfraclr,pfracld, &
3895	radocond, picefra, rain_lsc, snow_lsc, &
3896	frac_impa, frac_nucl, beta_prec_fisrt, &
3897	prfl, psfl, rhcl, &
3898	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3899	iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, temp_cltop, &
3900	qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
3901	Tcontr, qcontr, qcontr2, fcontrN, fcontrP , &
3902	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
3903	qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, dqrcol, dqrmelt, &
3904	dqrfreez, dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez)
3905
3906
3907	ELSE
3908
3909	CALL fisrtilp(klon,klev,phys_tstep,paprs,pplay, &
3910	t_seri, q_seri,ptconv,ratqs, &
3911	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, radocond, &
3912	rain_lsc, snow_lsc, &
3913	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3914	frac_impa, frac_nucl, beta_prec_fisrt, &
3915	prfl, psfl, rhcl, &
3916	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3917	iflag_ice_thermo, &
3918	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
3919
3920	ENDIF
3921	!
3922	WHERE (rain_lsc < 0) rain_lsc = 0.
3923	WHERE (snow_lsc < 0) snow_lsc = 0.
3924
3925	!+JLD
3926	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3927	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
3928	! & ,rain_lsc,snow_lsc
3929	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3930	!-JLD
3931	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,dqbs0,paprs, &
3932	'lsc',abortphy,flag_inhib_tend,itap,0)
3933	CALL prt_enerbil('lsc',itap)
3934	rain_num(:)=0.
3935	DO k = 1, klev
3936	DO i = 1, klon
3937	IF (ql_seri(i,k)>oliqmax) THEN
3938	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3939	ql_seri(i,k)=oliqmax
3940	ENDIF
3941	ENDDO
3942	ENDDO
3943	IF (nqo >= 3) THEN
3944	DO k = 1, klev
3945	DO i = 1, klon
3946	IF (qs_seri(i,k)>oicemax) THEN
3947	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3948	qs_seri(i,k)=oicemax
3949	ENDIF
3950	ENDDO
3951	ENDDO
3952	ENDIF
3953
3954
3955	!---------------------------------------------------------------------------
3956	DO k = 1, klev
3957	DO i = 1, klon
3958	cldfra(i,k) = rneb(i,k)
3959	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3960	!EV: en effet etrange, j'ajouterais aussi qs_seri
3961	! plus largement, je nettoierais (enleverrais) ces lignes
3962	IF (.NOT.new_oliq) radocond(i,k) = ql_seri(i,k)
3963	ENDDO
3964	ENDDO
3965
3966
3967	! Option to activate the radiative effect of blowing snow (ok_rad_bs)
3968	! makes sense only if the new large scale condensation scheme is active
3969	! with the ok_icefra_lscp flag active as well
3970
3971	IF (ok_bs .AND. ok_rad_bs) THEN
3972	IF (ok_new_lscp .AND. ok_icefra_lscp) THEN
3973	DO k=1,klev
3974	DO i=1,klon
3975	radocond(i,k)=radocond(i,k)+qbs_seri(i,k)
3976	picefra(i,k)=(radocond(i,k)*picefra(i,k)+qbs_seri(i,k))/(radocond(i,k))
3977	qbsfra=min(qbs_seri(i,k)/qbst_bs,1.0)
3978	cldfra(i,k)=max(cldfra(i,k),qbsfra)
3979	ENDDO
3980	ENDDO
3981	ELSE
3982	WRITE(lunout,*)"PAY ATTENTION, you try to activate the radiative effect of blowing snow"
3983	WRITE(lunout,*)"with ok_new_lscp=false and/or ok_icefra_lscp=false"
3984	abort_message='inconsistency in cloud phase for blowing snow'
3985	CALL abort_physic(modname,abort_message,1)
3986	ENDIF
3987
3988	ENDIF
3989
3990	IF (check) THEN
3991	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3992	WRITE(lunout,*)"apresilp=", za
3993	zx_t = 0.0
3994	za = 0.0
3995	DO i = 1, klon
3996	za = za + cell_area(i)/REAL(klon)
3997	zx_t = zx_t + (rain_lsc(i) &
3998	+ snow_lsc(i))*cell_area(i)/REAL(klon)
3999	ENDDO
4000	zx_t = zx_t/za*phys_tstep
4001	WRITE(lunout,*)"Precip=", zx_t
4002	ENDIF
4003
4004	IF (mydebug) THEN
4005	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4006	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4007	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4008	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4009	ENDIF
4010
4011	!
4012	!-------------------------------------------------------------------
4013	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
4014	!-------------------------------------------------------------------
4015
4016	! 1. NUAGES CONVECTIFS
4017	!
4018	!IM cf FH
4019	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
4020	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
4021	snow_tiedtke=0.
4022	! print*,'avant calcul de la pseudo precip '
4023	! print*,'iflag_cld_th',iflag_cld_th
4024	IF (iflag_cld_th.eq.-1) THEN
4025	rain_tiedtke=rain_con
4026	ELSE
4027	! print*,'calcul de la pseudo precip '
4028	rain_tiedtke=0.
4029	! print*,'calcul de la pseudo precip 0'
4030	DO k=1,klev
4031	DO i=1,klon
4032	IF (d_q_con(i,k).lt.0.) THEN
4033	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
4034	*(paprs(i,k)-paprs(i,k+1))/rg
4035	ENDIF
4036	ENDDO
4037	ENDDO
4038	ENDIF
4039	!
4040	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
4041	!
4042
4043	! Nuages diagnostiques pour Tiedtke
4044	CALL diagcld1(paprs,pplay, &
4045	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
4046	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
4047	diafra,dialiq)
4048	DO k = 1, klev
4049	DO i = 1, klon
4050	IF (diafra(i,k).GT.cldfra(i,k)) THEN
4051	radocond(i,k) = dialiq(i,k)
4052	cldfra(i,k) = diafra(i,k)
4053	ENDIF
4054	ENDDO
4055	ENDDO
4056
4057	ELSE IF (iflag_cld_th.ge.3) THEN
4058	! On prend pour les nuages convectifs le max du calcul de la
4059	! convection et du calcul du pas de temps precedent diminue d'un facteur
4060	! facttemps
4061	facteur = pdtphys *facttemps
4062	DO k=1,klev
4063	DO i=1,klon
4064	rnebcon(i,k)=rnebcon(i,k)*facteur
4065	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
4066	rnebcon(i,k)=rnebcon0(i,k)
4067	clwcon(i,k)=clwcon0(i,k)
4068	ENDIF
4069	ENDDO
4070	ENDDO
4071
4072	! On prend la somme des fractions nuageuses et des contenus en eau
4073
4074	IF (iflag_cld_th>=5) THEN
4075
4076	DO k=1,klev
4077	ptconvth(:,k)=fm_therm(:,k+1)>0.
4078	ENDDO
4079
4080	IF (iflag_coupl==4) THEN
4081
4082	! Dans le cas iflag_coupl==4, on prend la somme des convertures
4083	! convectives et lsc dans la partie des thermiques
4084	! Le controle par iflag_coupl est peut etre provisoire.
4085	DO k=1,klev
4086	DO i=1,klon
4087	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
4088	radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
4089	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
4090	ELSE IF (ptconv(i,k)) THEN
4091	cldfra(i,k)=rnebcon(i,k)
4092	radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
4093	ENDIF
4094	ENDDO
4095	ENDDO
4096
4097	ELSE IF (iflag_coupl==5) THEN
4098	DO k=1,klev
4099	DO i=1,klon
4100	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
4101	radocond(i,k)=radocond(i,k)+rnebcon(i,k)*clwcon(i,k)
4102	ENDDO
4103	ENDDO
4104
4105	ELSE
4106
4107	! Si on est sur un point touche par la convection
4108	! profonde et pas par les thermiques, on prend la
4109	! couverture nuageuse et l'eau nuageuse de la convection
4110	! profonde.
4111
4112	!IM/FH: 2011/02/23
4113	! definition des points sur lesquels ls thermiques sont actifs
4114
4115	DO k=1,klev
4116	DO i=1,klon
4117	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
4118	cldfra(i,k)=rnebcon(i,k)
4119	radocond(i,k)=rnebcon(i,k)*clwcon(i,k)
4120	ENDIF
4121	ENDDO
4122	ENDDO
4123
4124	ENDIF
4125
4126	ELSE
4127
4128	! Ancienne version
4129	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
4130	radocond(:,:)=radocond(:,:)+rnebcon(:,:)*clwcon(:,:)
4131	ENDIF
4132
4133	ENDIF
4134
4135	! plulsc(:)=0.
4136	! do k=1,klev,-1
4137	! do i=1,klon
4138	! zzz=prfl(:,k)+psfl(:,k)
4139	! if (.not.ptconvth.zzz.gt.0.)
4140	! enddo prfl, psfl,
4141	! enddo
4142	!
4143	! 2. NUAGES STARTIFORMES
4144	!
4145	IF (ok_stratus) THEN
4146	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
4147	DO k = 1, klev
4148	DO i = 1, klon
4149	IF (diafra(i,k).GT.cldfra(i,k)) THEN
4150	radocond(i,k) = dialiq(i,k)
4151	cldfra(i,k) = diafra(i,k)
4152	ENDIF
4153	ENDDO
4154	ENDDO
4155	ENDIF
4156	!
4157	! Precipitation totale
4158	!
4159	DO i = 1, klon
4160	rain_fall(i) = rain_con(i) + rain_lsc(i)
4161	snow_fall(i) = snow_con(i) + snow_lsc(i)
4162	ENDDO
4163	!
4164	! Calculer l'humidite relative pour diagnostique
4165	!
4166	DO k = 1, klev
4167	DO i = 1, klon
4168	zx_t = t_seri(i,k)
4169	IF (thermcep) THEN
4170	!! if (iflag_ice_thermo.eq.0) then !jyg
4171	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
4172	!! else !jyg
4173	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
4174	!! endif !jyg
4175	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
4176	zx_qs = MIN(0.5,zx_qs)
4177	zcor = 1./(1.-retv*zx_qs)
4178	zx_qs = zx_qs*zcor
4179	ELSE
4180	!! IF (zx_t.LT.t_coup) THEN !jyg
4181	IF (zx_t.LT.rtt) THEN !jyg
4182	zx_qs = qsats(zx_t)/pplay(i,k)
4183	ELSE
4184	zx_qs = qsatl(zx_t)/pplay(i,k)
4185	ENDIF
4186	ENDIF
4187	zx_rh(i,k) = q_seri(i,k)/zx_qs
4188	IF (iflag_ice_thermo .GT. 0) THEN
4189	zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
4190	zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
4191	ENDIF
4192	zqsat(i,k)=zx_qs
4193	ENDDO
4194	ENDDO
4195
4196	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
4197	! equivalente a 2m (tpote) pour diagnostique
4198	!
4199	DO i = 1, klon
4200	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
4201	IF (thermcep) THEN
4202	IF(zt2m(i).LT.RTT) then
4203	Lheat=RLSTT
4204	ELSE
4205	Lheat=RLVTT
4206	ENDIF
4207	ELSE
4208	IF (zt2m(i).LT.RTT) THEN
4209	Lheat=RLSTT
4210	ELSE
4211	Lheat=RLVTT
4212	ENDIF
4213	ENDIF
4214	tpote(i) = tpot(i)* &
4215	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
4216	ENDDO
4217
4218	IF (ANY(type_trac == ['inca','inco'])) THEN ! ModThL
4219	#ifdef INCA
4220	CALL VTe(VTphysiq)
4221	CALL VTb(VTinca)
4222	calday = REAL(days_elapsed + 1) + jH_cur
4223
4224	CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
4225	CALL AEROSOL_METEO_CALC( &
4226	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
4227	prfl,psfl,pctsrf,cell_area, &
4228	latitude_deg,longitude_deg,u10m,v10m)
4229
4230	zxsnow_dummy(:) = 0.0
4231
4232	CALL chemhook_begin (calday, &
4233	days_elapsed+1, &
4234	jH_cur, &
4235	pctsrf(1,1), &
4236	latitude_deg, &
4237	longitude_deg, &
4238	cell_area, &
4239	paprs, &
4240	pplay, &
4241	coefh(1:klon,1:klev,is_ave), &
4242	pphi, &
4243	t_seri, &
4244	u, &
4245	v, &
4246	rot, &
4247	wo(:, :, 1), &
4248	q_seri, &
4249	zxtsol, &
4250	zt2m, &
4251	zxsnow_dummy, &
4252	solsw, &
4253	albsol1, &
4254	rain_fall, &
4255	snow_fall, &
4256	itop_con, &
4257	ibas_con, &
4258	cldfra, &
4259	nbp_lon, &
4260	nbp_lat-1, &
4261	tr_seri(:,:,1+nqCO2:nbtr), &
4262	ftsol, &
4263	paprs, &
4264	cdragh, &
4265	cdragm, &
4266	pctsrf, &
4267	pdtphys, &
4268	itap)
4269
4270	CALL VTe(VTinca)
4271	CALL VTb(VTphysiq)
4272	#endif
4273	ENDIF !type_trac = inca or inco
4274	IF (type_trac == 'repr') THEN
4275	#ifdef REPROBUS
4276	!CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap)
4277	CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap)
4278	#endif
4279	ENDIF
4280
4281	!
4282	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
4283	!
4284	IF (MOD(itaprad,radpas).EQ.0) THEN
4285
4286	!
4287	!jq - introduce the aerosol direct and first indirect radiative forcings
4288	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
4289	IF (flag_aerosol .GT. 0) THEN
4290	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4291	IF (.NOT. aerosol_couple) THEN
4292	!
4293	CALL readaerosol_optic( &
4294	debut, flag_aerosol, itap, jD_cur-jD_ref, &
4295	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4296	mass_solu_aero, mass_solu_aero_pi, &
4297	tau_aero, piz_aero, cg_aero, &
4298	tausum_aero, tau3d_aero)
4299	ENDIF
4300	ELSE IF (iflag_rrtm .EQ.1) THEN ! RRTM radiation
4301	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
4302	abort_message='config_inca=aero et rrtm=1 impossible'
4303	CALL abort_physic(modname,abort_message,1)
4304	ELSE
4305	!
4306	#ifdef CPP_RRTM
4307	IF (NSW.EQ.6) THEN
4308	!--new aerosol properties SW and LW
4309	!
4310	#ifdef CPP_Dust
4311	!--SPL aerosol model
4312	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
4313	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
4314	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4315	tausum_aero, tau3d_aero)
4316	#else
4317	!--climatologies or INCA aerosols
4318	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
4319	flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
4320	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4321	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
4322	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4323	tausum_aero, drytausum_aero, tau3d_aero)
4324	#endif
4325
4326	IF (flag_aerosol .EQ. 7) THEN
4327	CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg, &
4328	tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
4329	ENDIF
4330
4331	!
4332	ELSE IF (NSW.EQ.2) THEN
4333	!--for now we use the old aerosol properties
4334	!
4335	CALL readaerosol_optic( &
4336	debut, flag_aerosol, itap, jD_cur-jD_ref, &
4337	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4338	mass_solu_aero, mass_solu_aero_pi, &
4339	tau_aero, piz_aero, cg_aero, &
4340	tausum_aero, tau3d_aero)
4341	!
4342	!--natural aerosols
4343	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
4344	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
4345	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
4346	!--all aerosols
4347	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
4348	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
4349	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
4350	!
4351	!--no LW optics
4352	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4353	!
4354	ELSE
4355	abort_message='Only NSW=2 or 6 are possible with ' &
4356	// 'aerosols and iflag_rrtm=1'
4357	CALL abort_physic(modname,abort_message,1)
4358	ENDIF
4359	#else
4360	abort_message='You should compile with -rrtm if running ' &
4361	// 'with iflag_rrtm=1'
4362	CALL abort_physic(modname,abort_message,1)
4363	#endif
4364	!
4365	ENDIF
4366	ELSE IF (iflag_rrtm .EQ.2) THEN ! ecrad RADIATION
4367	#ifdef CPP_ECRAD
4368	!--climatologies or INCA aerosols
4369	CALL readaerosol_optic_ecrad( debut, aerosol_couple, ok_alw, ok_volcan, &
4370	flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
4371	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
4372	tr_seri, mass_solu_aero, mass_solu_aero_pi)
4373	#else
4374	abort_message='You should compile with -rad ecrad if running with iflag_rrtm=2'
4375	CALL abort_physic(modname,abort_message,1)
4376	#endif
4377	ENDIF
4378
4379	ELSE !--flag_aerosol = 0
4380	tausum_aero(:,:,:) = 0.
4381	drytausum_aero(:,:) = 0.
4382	mass_solu_aero(:,:) = 0.
4383	mass_solu_aero_pi(:,:) = 0.
4384	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
4385	tau_aero(:,:,:,:) = 1.e-15
4386	piz_aero(:,:,:,:) = 1.
4387	cg_aero(:,:,:,:) = 0.
4388	ELSE
4389	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
4390	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
4391	piz_aero_sw_rrtm(:,:,:,:) = 1.0
4392	cg_aero_sw_rrtm(:,:,:,:) = 0.0
4393	ENDIF
4394	ENDIF
4395	!
4396	!--WMO criterion to determine tropopause
4397	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
4398	!
4399	!--STRAT AEROSOL
4400	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
4401	IF (flag_aerosol_strat.GT.0) THEN
4402	IF (prt_level .GE.10) THEN
4403	PRINT *,'appel a readaerosolstrat', mth_cur
4404	ENDIF
4405	IF (iflag_rrtm.EQ.0) THEN
4406	IF (flag_aerosol_strat.EQ.1) THEN
4407	CALL readaerosolstrato(debut)
4408	ELSE
4409	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
4410	CALL abort_physic(modname,abort_message,1)
4411	ENDIF
4412	ELSE
4413	#ifdef CPP_RRTM
4414	#ifndef CPP_StratAer
4415	!--prescribed strat aerosols
4416	!--only in the case of non-interactive strat aerosols
4417	IF (flag_aerosol_strat.EQ.1) THEN
4418	CALL readaerosolstrato1_rrtm(debut)
4419	ELSEIF (flag_aerosol_strat.EQ.2) THEN
4420	CALL readaerosolstrato2_rrtm(debut, ok_volcan)
4421	ELSE
4422	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
4423	CALL abort_physic(modname,abort_message,1)
4424	ENDIF
4425	#endif
4426	#else
4427	abort_message='You should compile with -rrtm if running ' &
4428	// 'with iflag_rrtm=1'
4429	CALL abort_physic(modname,abort_message,1)
4430	#endif
4431	ENDIF
4432	ELSE
4433	tausum_aero(:,:,id_STRAT_phy) = 0.
4434	ENDIF
4435	!
4436	#ifdef CPP_RRTM
4437	#ifdef CPP_StratAer
4438	!--compute stratospheric mask
4439	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
4440	!--interactive strat aerosols
4441	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
4442	#endif
4443	#endif
4444	!--fin STRAT AEROSOL
4445	!
4446
4447	! Calculer les parametres optiques des nuages et quelques
4448	! parametres pour diagnostiques:
4449	!
4450	IF (aerosol_couple.AND.config_inca=='aero') THEN
4451	mass_solu_aero(:,:) = ccm(:,:,1)
4452	mass_solu_aero_pi(:,:) = ccm(:,:,2)
4453	ENDIF
4454
4455	!Rajout appel a interface calcul proprietes optiques des nuages
4456	CALL call_cloud_optics_prop(klon, klev, ok_newmicro, &
4457	paprs, pplay, t_seri, radocond, picefra, cldfra, &
4458	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
4459	flwp, fiwp, flwc, fiwc, ok_aie, &
4460	mass_solu_aero, mass_solu_aero_pi, &
4461	cldtaupi, distcltop, temp_cltop, re, fl, ref_liq, ref_ice, &
4462	ref_liq_pi, ref_ice_pi, scdnc, cldncl, reffclwtop, lcc, reffclws, &
4463	reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra, &
4464	zfice, dNovrN, ptconv, rnebcon, clwcon)
4465
4466	!
4467	!IM betaCRF
4468	!
4469	cldtaurad = cldtau
4470	cldtaupirad = cldtaupi
4471	cldemirad = cldemi
4472	cldfrarad = cldfra
4473
4474	!
4475	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
4476	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
4477	!
4478	! global
4479	!
4480	!IM 251017 begin
4481	! print*,'physiq betaCRF global zdtime=',zdtime
4482	!IM 251017 end
4483	DO k=1, klev
4484	DO i=1, klon
4485	IF (pplay(i,k).GE.pfree) THEN
4486	beta(i,k) = beta_pbl
4487	ELSE
4488	beta(i,k) = beta_free
4489	ENDIF
4490	IF (mskocean_beta) THEN
4491	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4492	ENDIF
4493	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4494	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4495	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4496	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4497	ENDDO
4498	ENDDO
4499	!
4500	ELSE
4501	!
4502	! regional
4503	!
4504	DO k=1, klev
4505	DO i=1,klon
4506	!
4507	IF (longitude_deg(i).ge.lon1_beta.AND. &
4508	longitude_deg(i).le.lon2_beta.AND. &
4509	latitude_deg(i).le.lat1_beta.AND. &
4510	latitude_deg(i).ge.lat2_beta) THEN
4511	IF (pplay(i,k).GE.pfree) THEN
4512	beta(i,k) = beta_pbl
4513	ELSE
4514	beta(i,k) = beta_free
4515	ENDIF
4516	IF (mskocean_beta) THEN
4517	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4518	ENDIF
4519	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4520	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4521	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4522	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4523	ENDIF
4524	!
4525	ENDDO
4526	ENDDO
4527	!
4528	ENDIF
4529
4530	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
4531	IF (ok_chlorophyll) THEN
4532	print*,"-- reading chlorophyll"
4533	CALL readchlorophyll(debut)
4534	ENDIF
4535
4536	!--if ok_suntime_rrtm we use ancillay data for RSUN
4537	!--previous values are therefore overwritten
4538	!--this is needed for CMIP6 runs
4539	!--and only possible for new radiation scheme
4540	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
4541	#ifdef CPP_RRTM
4542	CALL read_rsun_rrtm(debut)
4543	#endif
4544	ENDIF
4545
4546	IF (mydebug) THEN
4547	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4548	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4549	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4550	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4551	ENDIF
4552
4553	!
4554	!sonia : If Iflag_radia >=2, pertubation of some variables
4555	!input to radiation (DICE)
4556	!
4557	IF (iflag_radia .ge. 2) THEN
4558	zsav_tsol (:) = zxtsol(:)
4559	CALL perturb_radlwsw(zxtsol,iflag_radia)
4560	ENDIF
4561
4562	IF (aerosol_couple.AND.config_inca=='aero') THEN
4563	#ifdef INCA
4564	CALL radlwsw_inca &
4565	(chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
4566	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
4567	size(wo,3), wo, &
4568	cldfrarad, cldemirad, cldtaurad, &
4569	heat,heat0,cool,cool0,albpla, &
4570	topsw,toplw,solsw,sollw, &
4571	sollwdown, &
4572	topsw0,toplw0,solsw0,sollw0, &
4573	lwdn0, lwdn, lwup0, lwup, &
4574	swdn0, swdn, swup0, swup, &
4575	ok_ade, ok_aie, &
4576	tau_aero, piz_aero, cg_aero, &
4577	topswad_aero, solswad_aero, &
4578	topswad0_aero, solswad0_aero, &
4579	topsw_aero, topsw0_aero, &
4580	solsw_aero, solsw0_aero, &
4581	cldtaupirad, &
4582	topswai_aero, solswai_aero)
4583	#endif
4584	ELSE
4585	!
4586	!IM calcul radiatif pour le cas actuel
4587	!
4588	RCO2 = RCO2_act
4589	RCH4 = RCH4_act
4590	RN2O = RN2O_act
4591	RCFC11 = RCFC11_act
4592	RCFC12 = RCFC12_act
4593	!
4594	!--interactive CO2 in ppm from carbon cycle
4595	IF (carbon_cycle_rad) RCO2=RCO2_glo
4596	!
4597	IF (prt_level .GE.10) THEN
4598	print *,' ->radlwsw, number 1 '
4599	ENDIF
4600	!
4601	! AI namelist utilise pour l appel principal de radlwsw (ecrad)
4602	namelist_ecrad_file='namelist_ecrad'
4603	!
4604	CALL radlwsw &
4605	(debut, dist, rmu0, fract, &
4606	!albedo SB >>>
4607	! paprs, pplay,zxtsol,albsol1, albsol2, &
4608	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4609	!albedo SB <<<
4610	t_seri,q_seri,wo, &
4611	cldfrarad, cldemirad, cldtaurad, &
4612	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4613	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4614	tau_aero, piz_aero, cg_aero, &
4615	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4616	! Rajoute par OB pour RRTM
4617	tau_aero_lw_rrtm, &
4618	cldtaupirad, &
4619	! zqsat, flwcrad, fiwcrad, &
4620	zqsat, flwc, fiwc, &
4621	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4622	namelist_ecrad_file, &
4623	heat,heat0,cool,cool0,albpla, &
4624	heat_volc,cool_volc, &
4625	topsw,toplw,solsw,solswfdiff,sollw, &
4626	sollwdown, &
4627	topsw0,toplw0,solsw0,sollw0, &
4628	lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, &
4629	swdnc0, swdn0, swdn, swupc0, swup0, swup, &
4630	topswad_aero, solswad_aero, &
4631	topswai_aero, solswai_aero, &
4632	topswad0_aero, solswad0_aero, &
4633	topsw_aero, topsw0_aero, &
4634	solsw_aero, solsw0_aero, &
4635	topswcf_aero, solswcf_aero, &
4636	!-C. Kleinschmitt for LW diagnostics
4637	toplwad_aero, sollwad_aero,&
4638	toplwai_aero, sollwai_aero, &
4639	toplwad0_aero, sollwad0_aero,&
4640	!-end
4641	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4642	ZSWFT0_i, ZFSDN0, ZFSUP0, &
4643	cloud_cover_sw)
4644
4645	!lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
4646	!schemes
4647	toplw = toplw + betalwoff * (toplw0 - toplw)
4648	sollw = sollw + betalwoff * (sollw0 - sollw)
4649	lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
4650	lwup = lwup + betalwoff * (lwup0 - lwup)
4651	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
4652	sollwdown(:))
4653	cool = cool + betalwoff * (cool0 - cool)
4654
4655	IF (.NOT. using_xios) THEN
4656	!
4657	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
4658	!IM des taux doit etre different du taux actuel
4659	!IM Par defaut on a les taux perturbes egaux aux taux actuels
4660	!
4661	IF (RCO2_per.NE.RCO2_act.OR. &
4662	RCH4_per.NE.RCH4_act.OR. &
4663	RN2O_per.NE.RN2O_act.OR. &
4664	RCFC11_per.NE.RCFC11_act.OR. &
4665	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
4666	ENDIF
4667	!
4668	IF (ok_4xCO2atm) THEN
4669	!
4670	RCO2 = RCO2_per
4671	RCH4 = RCH4_per
4672	RN2O = RN2O_per
4673	RCFC11 = RCFC11_per
4674	RCFC12 = RCFC12_per
4675	!
4676	IF (prt_level .GE.10) THEN
4677	print *,' ->radlwsw, number 2 '
4678	ENDIF
4679	!
4680	namelist_ecrad_file='namelist_ecrad'
4681	!
4682	CALL radlwsw &
4683	(debut, dist, rmu0, fract, &
4684	!albedo SB >>>
4685	! paprs, pplay,zxtsol,albsol1, albsol2, &
4686	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4687	!albedo SB <<<
4688	t_seri,q_seri,wo, &
4689	cldfrarad, cldemirad, cldtaurad, &
4690	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4691	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4692	tau_aero, piz_aero, cg_aero, &
4693	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4694	! Rajoute par OB pour RRTM
4695	tau_aero_lw_rrtm, &
4696	cldtaupi, &
4697	! zqsat, flwcrad, fiwcrad, &
4698	zqsat, flwc, fiwc, &
4699	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4700	namelist_ecrad_file, &
4701	heatp,heat0p,coolp,cool0p,albplap, &
4702	heat_volc,cool_volc, &
4703	topswp,toplwp,solswp,solswfdiffp,sollwp, &
4704	sollwdownp, &
4705	topsw0p,toplw0p,solsw0p,sollw0p, &
4706	lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp, &
4707	swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
4708	topswad_aerop, solswad_aerop, &
4709	topswai_aerop, solswai_aerop, &
4710	topswad0_aerop, solswad0_aerop, &
4711	topsw_aerop, topsw0_aerop, &
4712	solsw_aerop, solsw0_aerop, &
4713	topswcf_aerop, solswcf_aerop, &
4714	!-C. Kleinschmitt for LW diagnostics
4715	toplwad_aerop, sollwad_aerop,&
4716	toplwai_aerop, sollwai_aerop, &
4717	toplwad0_aerop, sollwad0_aerop,&
4718	!-end
4719	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4720	ZSWFT0_i, ZFSDN0, ZFSUP0, &
4721	cloud_cover_sw)
4722	ENDIF !ok_4xCO2atm
4723
4724	! A.I aout 2023
4725	! Effet 3D des nuages Ecrad
4726	! a passer : nom du ficher namelist et cles ok_3Deffect
4727	! a declarer comme iflag_rrtm et a lire dans physiq.def
4728	#ifdef CPP_ECRAD
4729	IF (ok_3Deffect) then
4730	! print*,'ok_3Deffect = ',ok_3Deffect
4731	namelist_ecrad_file='namelist_ecrad_s2'
4732	CALL radlwsw &
4733	(debut, dist, rmu0, fract, &
4734	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4735	t_seri,q_seri,wo, &
4736	cldfrarad, cldemirad, cldtaurad, &
4737	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, flag_volc_surfstrat, &
4738	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, &
4739	tau_aero, piz_aero, cg_aero, &
4740	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4741	tau_aero_lw_rrtm, &
4742	cldtaupi, &
4743	zqsat, flwc, fiwc, &
4744	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4745	namelist_ecrad_file, &
4746	! A modifier
4747	heat_s2,heat0_s2,cool_s2,cool0_s2,albpla_s2, &
4748	heat_volc,cool_volc, &
4749	topsw_s2,toplw_s2,solsw_s2,solswfdiff_s2,sollw_s2, &
4750	sollwdown_s2, &
4751	topsw0_s2,toplw0_s2,solsw0_s2,sollw0_s2, &
4752	lwdnc0_s2, lwdn0_s2, lwdn_s2, lwupc0_s2, lwup0_s2, lwup_s2, &
4753	swdnc0_s2, swdn0_s2, swdn_s2, swupc0_s2, swup0_s2, swup_s2, &
4754	topswad_aero_s2, solswad_aero_s2, &
4755	topswai_aero_s2, solswai_aero_s2, &
4756	topswad0_aero_s2, solswad0_aero_s2, &
4757	topsw_aero_s2, topsw0_aero_s2, &
4758	solsw_aero_s2, solsw0_aero_s2, &
4759	topswcf_aero_s2, solswcf_aero_s2, &
4760	!-C. Kleinschmitt for LW diagnostics
4761	toplwad_aero_s2, sollwad_aero_s2,&
4762	toplwai_aero_s2, sollwai_aero_s2, &
4763	toplwad0_aero_s2, sollwad0_aero_s2,&
4764	!-end
4765	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4766	ZSWFT0_i, ZFSDN0, ZFSUP0, &
4767	cloud_cover_sw_s2)
4768	ENDIF ! ok_3Deffect
4769	#endif
4770
4771	ENDIF ! aerosol_couple
4772	itaprad = 0
4773	!
4774	! If Iflag_radia >=2, reset pertubed variables
4775	!
4776	IF (iflag_radia .ge. 2) THEN
4777	zxtsol(:) = zsav_tsol (:)
4778	ENDIF
4779	ENDIF ! MOD(itaprad,radpas)
4780	itaprad = itaprad + 1
4781
4782	IF (iflag_radia.eq.0) THEN
4783	IF (prt_level.ge.9) THEN
4784	PRINT *,'--------------------------------------------------'
4785	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4786	PRINT *,'>>>> heat et cool mis a zero '
4787	PRINT *,'--------------------------------------------------'
4788	ENDIF
4789	heat=0.
4790	cool=0.
4791	sollw=0. ! MPL 01032011
4792	solsw=0.
4793	radsol=0.
4794	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4795	swup0=0.
4796	lwup=0.
4797	lwup0=0.
4798	lwdn=0.
4799	lwdn0=0.
4800	ENDIF
4801
4802	!
4803	! Calculer radsol a l'exterieur de radlwsw
4804	! pour prendre en compte le cycle diurne
4805	! recode par Olivier Boucher en sept 2015
4806	!
4807	radsol=solsw*swradcorr+sollw
4808
4809	IF (ok_4xCO2atm) THEN
4810	radsolp=solswp*swradcorr+sollwp
4811	ENDIF
4812
4813	!
4814	! Ajouter la tendance des rayonnements (tous les pas)
4815	! avec une correction pour le cycle diurne dans le SW
4816	!
4817
4818	DO k=1, klev
4819	d_t_swr(:,k)=swradcorr(:)heat(:,k)phys_tstep/RDAY
4820	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)phys_tstep/RDAY
4821	d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
4822	d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
4823	ENDDO
4824
4825	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,dqbs0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4826	CALL prt_enerbil('SW',itap)
4827	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,dqbs0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4828	CALL prt_enerbil('LW',itap)
4829
4830	!
4831	IF (mydebug) THEN
4832	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4833	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4834	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4835	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4836	ENDIF
4837
4838	! Calculer l'hydrologie de la surface
4839	!
4840	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4841	! . agesno, ftsol,fqsurf,fsnow, ruis)
4842	!
4843
4844	!
4845	! Calculer le bilan du sol et la derive de temperature (couplage)
4846	!
4847	DO i = 1, klon
4848	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4849	! a la demande de JLD
4850	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4851	ENDDO
4852	!
4853	!moddeblott(jan95)
4854	! Appeler le programme de parametrisation de l'orographie
4855	! a l'echelle sous-maille:
4856	!
4857	IF (prt_level .GE.10) THEN
4858	print *,' call orography ? ', ok_orodr
4859	ENDIF
4860	!
4861	IF (ok_orodr) THEN
4862	!
4863	! selection des points pour lesquels le shema est actif:
4864	igwd=0
4865	DO i=1,klon
4866	itest(i)=0
4867	zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
4868	!zrel_oro: relative mountain height wrt relief explained by mean slope
4869	! -> condition on zrel_oro can deactivate the drag on tilted planar terrains
4870	! such as ice sheets (work by V. Wiener)
4871	! zpmm_orodr_t and zstd_orodr_t are activation thresholds set by F. Lott to
4872	! earn computation time but they are not physical.
4873	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
4874	itest(i)=1
4875	igwd=igwd+1
4876	idx(igwd)=i
4877	ENDIF
4878	ENDDO
4879	! igwdim=MAX(1,igwd)
4880	!
4881	IF (ok_strato) THEN
4882
4883	CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
4884	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4885	igwd,idx,itest, &
4886	t_seri, u_seri, v_seri, &
4887	zulow, zvlow, zustrdr, zvstrdr, &
4888	d_t_oro, d_u_oro, d_v_oro)
4889
4890	ELSE
4891	CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
4892	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4893	igwd,idx,itest, &
4894	t_seri, u_seri, v_seri, &
4895	zulow, zvlow, zustrdr, zvstrdr, &
4896	d_t_oro, d_u_oro, d_v_oro)
4897	ENDIF
4898	!
4899	! ajout des tendances
4900	!-----------------------------------------------------------------------
4901	! ajout des tendances de la trainee de l'orographie
4902	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,dqbs0,paprs,'oro', &
4903	abortphy,flag_inhib_tend,itap,0)
4904	CALL prt_enerbil('oro',itap)
4905	!----------------------------------------------------------------------
4906	!
4907	ENDIF ! fin de test sur ok_orodr
4908	!
4909	IF (mydebug) THEN
4910	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4911	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4912	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4913	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4914	ENDIF
4915
4916	IF (ok_orolf) THEN
4917	!
4918	! selection des points pour lesquels le shema est actif:
4919	igwd=0
4920	DO i=1,klon
4921	itest(i)=0
4922	!zrel_oro: relative mountain height wrt relief explained by mean slope
4923	! -> condition on zrel_oro can deactivate the lifting on tilted planar terrains
4924	! such as ice sheets (work by V. Wiener)
4925	zrel_oro(i)=zstd(i)/(max(zsig(i),1.E-8)*sqrt(cell_area(i)))
4926	IF (((zpic(i)-zmea(i)).GT.zpmm_orolf_t).AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
4927	itest(i)=1
4928	igwd=igwd+1
4929	idx(igwd)=i
4930	ENDIF
4931	ENDDO
4932	! igwdim=MAX(1,igwd)
4933	!
4934	IF (ok_strato) THEN
4935
4936	CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
4937	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4938	igwd,idx,itest, &
4939	t_seri, u_seri, v_seri, &
4940	zulow, zvlow, zustrli, zvstrli, &
4941	d_t_lif, d_u_lif, d_v_lif )
4942
4943	ELSE
4944	CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
4945	latitude_deg,zmea,zstd,zpic, &
4946	itest, &
4947	t_seri, u_seri, v_seri, &
4948	zulow, zvlow, zustrli, zvstrli, &
4949	d_t_lif, d_u_lif, d_v_lif)
4950	ENDIF
4951
4952	! ajout des tendances de la portance de l'orographie
4953	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, dqbs0,paprs, &
4954	'lif', abortphy,flag_inhib_tend,itap,0)
4955	CALL prt_enerbil('lif',itap)
4956	ENDIF ! fin de test sur ok_orolf
4957
4958	IF (ok_hines) then
4959	! HINES GWD PARAMETRIZATION
4960	east_gwstress=0.
4961	west_gwstress=0.
4962	du_gwd_hines=0.
4963	dv_gwd_hines=0.
4964	CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
4965	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4966	du_gwd_hines, dv_gwd_hines)
4967	zustr_gwd_hines=0.
4968	zvstr_gwd_hines=0.
4969	DO k = 1, klev
4970	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
4971	* (paprs(:, k)-paprs(:, k+1))/rg
4972	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
4973	* (paprs(:, k)-paprs(:, k+1))/rg
4974	ENDDO
4975
4976	d_t_hin(:, :)=0.
4977	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4978	dqi0, dqbs0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4979	CALL prt_enerbil('hin',itap)
4980	ENDIF
4981
4982	IF (.not. ok_hines .and. ok_gwd_rando) then
4983	! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
4984	CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
4985	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4986	dv_gwd_front, east_gwstress, west_gwstress)
4987	zustr_gwd_front=0.
4988	zvstr_gwd_front=0.
4989	DO k = 1, klev
4990	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
4991	* (paprs(:, k)-paprs(:, k+1))/rg
4992	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
4993	* (paprs(:, k)-paprs(:, k+1))/rg
4994	ENDDO
4995
4996	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, dqbs0, &
4997	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4998	CALL prt_enerbil('front_gwd_rando',itap)
4999	ENDIF
5000
5001	IF (ok_gwd_rando) THEN
5002	CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
5003	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
5004	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
5005	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, dqbs0, &
5006	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
5007	CALL prt_enerbil('flott_gwd_rando',itap)
5008	zustr_gwd_rando=0.
5009	zvstr_gwd_rando=0.
5010	DO k = 1, klev
5011	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
5012	* (paprs(:, k)-paprs(:, k+1))/rg
5013	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
5014	* (paprs(:, k)-paprs(:, k+1))/rg
5015	ENDDO
5016	ENDIF
5017
5018	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
5019
5020	IF (mydebug) THEN
5021	CALL writefield_phy('u_seri',u_seri,nbp_lev)
5022	CALL writefield_phy('v_seri',v_seri,nbp_lev)
5023	CALL writefield_phy('t_seri',t_seri,nbp_lev)
5024	CALL writefield_phy('q_seri',q_seri,nbp_lev)
5025	ENDIF
5026
5027	DO i = 1, klon
5028	zustrph(i)=0.
5029	zvstrph(i)=0.
5030	ENDDO
5031	DO k = 1, klev
5032	DO i = 1, klon
5033	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
5034	(paprs(i,k)-paprs(i,k+1))/rg
5035	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
5036	(paprs(i,k)-paprs(i,k+1))/rg
5037	ENDDO
5038	ENDDO
5039	!
5040	!IM calcul composantes axiales du moment angulaire et couple des montagnes
5041	!
5042	IF (is_sequential .and. ok_orodr) THEN
5043	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
5044	ra,rg,romega, &
5045	latitude_deg,longitude_deg,pphis, &
5046	zustrdr,zustrli,zustrph, &
5047	zvstrdr,zvstrli,zvstrph, &
5048	paprs,u,v, &
5049	aam, torsfc)
5050	ENDIF
5051	!IM cf. FLott END
5052	!DC Calcul de la tendance due au methane
5053	IF (ok_qch4) THEN
5054	! d_q_ch4: H2O source from CH4 in MMR/s (mass mixing ratio/s or kg H2O/kg air/s)
5055	#ifdef CPP_StratAer
5056	CALL stratH2O_methox(debut,paprs,d_q_ch4)
5057	#else
5058	! ECMWF routine METHOX
5059	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
5060	#endif
5061	! add humidity tendency due to methane
5062	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
5063	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, dqbs0, paprs, &
5064	'q_ch4', abortphy,flag_inhib_tend,itap,0)
5065	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep ! update with H2O conserv done in add_phys_tend
5066	ENDIF
5067	!
5068	!
5069	#ifdef CPP_StratAer
5070	IF (ok_qemiss) THEN
5071	flh2o=1
5072	IF(flag_verbose_strataer) THEN
5073	print *,'IN physiq_mod: ok_qemiss =yes (',ok_qemiss,'), flh2o=',flh2o
5074	print *,'IN physiq_mod: flag_emit=',flag_emit,', nErupt=',nErupt
5075	print *,'IN physiq_mod: nAerErupt=',nAerErupt
5076	ENDIF
5077
5078	SELECT CASE(flag_emit)
5079	CASE(1) ! emission volc H2O in LMDZ
5080	DO ieru=1, nErupt
5081	IF (year_cur==year_emit_vol(ieru).AND.&
5082	mth_cur==mth_emit_vol(ieru).AND.&
5083	day_cur>=day_emit_vol(ieru).AND.&
5084	day_cur<(day_emit_vol(ieru)+injdur)) THEN
5085
5086	IF(flag_verbose_strataer) print *,'IN physiq_mod: date=',year_cur,mth_cur,day_cur
5087	! initialisation of q tendency emission
5088	d_q_emiss(:,:)=0.
5089	! daily injection mass emission - NL
5090	m_H2O_emiss_vol_daily = m_H2O_emiss_vol(ieru)/(REAL(injdur)&
5091	*REAL(ponde_lonlat_vol(ieru)))
5092	!
5093	CALL STRATEMIT(pdtphys,pdtphys,latitude_deg,longitude_deg,t_seri,&
5094	pplay,paprs,tr_seri,&
5095	m_H2O_emiss_vol_daily,&
5096	xlat_min_vol(ieru),xlat_max_vol(ieru),&
5097	xlon_min_vol(ieru),xlon_max_vol(ieru),&
5098	altemiss_vol(ieru),sigma_alt_vol(ieru),1,1.,&
5099	nAerErupt+1,0)
5100
5101	IF(flag_verbose_strataer) print *,'IN physiq_mod: min max d_q_emiss=',&
5102	minval(d_q_emiss),maxval(d_q_emiss)
5103
5104	CALL add_phys_tend(du0, dv0, dt0, d_q_emiss, dql0, dqi0, dqbs0, paprs, &
5105	'q_emiss',abortphy,flag_inhib_tend,itap,0)
5106	IF (abortphy==1) Print*,'ERROR ABORT TEND EMISS'
5107	ENDIF
5108	ENDDO
5109	flh2o=0
5110	END SELECT ! emission scenario (flag_emit)
5111	ENDIF
5112	#endif
5113
5114	!===============================================================
5115	! Additional tendency of TKE due to orography
5116	!===============================================================
5117	!
5118	! Inititialization
5119	!------------------
5120
5121	addtkeoro=0
5122	CALL getin_p('addtkeoro',addtkeoro)
5123
5124	IF (prt_level.ge.5) &
5125	print*,'addtkeoro', addtkeoro
5126
5127	alphatkeoro=1.
5128	CALL getin_p('alphatkeoro',alphatkeoro)
5129	alphatkeoro=min(max(0.,alphatkeoro),1.)
5130
5131	smallscales_tkeoro=.FALSE.
5132	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
5133
5134
5135	dtadd(:,:)=0.
5136	duadd(:,:)=0.
5137	dvadd(:,:)=0.
5138
5139	! Choices for addtkeoro:
5140	! ** 0 no TKE tendency from orography
5141	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
5142	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
5143	!
5144
5145	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
5146	! -------------------------------------------
5147
5148
5149	! selection des points pour lesquels le schema est actif:
5150
5151
5152	IF (addtkeoro .EQ. 1 ) THEN
5153
5154	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
5155	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
5156
5157	ELSE IF (addtkeoro .EQ. 2) THEN
5158
5159	IF (smallscales_tkeoro) THEN
5160	igwd=0
5161	DO i=1,klon
5162	itest(i)=0
5163	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
5164	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
5165	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
5166	IF ((zstd(i).GT.1.0) .AND.(zrel_oro(i).LE.zrel_oro_t)) THEN
5167	itest(i)=1
5168	igwd=igwd+1
5169	idx(igwd)=i
5170	ENDIF
5171	ENDDO
5172
5173	ELSE
5174
5175	igwd=0
5176	DO i=1,klon
5177	itest(i)=0
5178	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
5179	itest(i)=1
5180	igwd=igwd+1
5181	idx(igwd)=i
5182	ENDIF
5183	ENDDO
5184
5185	ENDIF
5186
5187	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
5188	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
5189	igwd,idx,itest, &
5190	t_seri, u_seri, v_seri, &
5191	zulow, zvlow, zustrdr, zvstrdr, &
5192	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
5193
5194	zustrdr(:)=0.
5195	zvstrdr(:)=0.
5196	zulow(:)=0.
5197	zvlow(:)=0.
5198
5199	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
5200	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
5201	ENDIF
5202
5203
5204	! TKE update from subgrid temperature and wind tendencies
5205	!----------------------------------------------------------
5206	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
5207
5208
5209	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
5210	!
5211	! Prevent pbl_tke_w from becoming negative
5212	wake_delta_pbl_tke(:,:,:) = max(wake_delta_pbl_tke(:,:,:), -pbl_tke(:,:,:))
5213	!
5214
5215	ENDIF
5216	! -----
5217	!===============================================================
5218
5219
5220	!====================================================================
5221	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
5222	!====================================================================
5223	! Abderrahmane 24.08.09
5224
5225	IF (ok_cosp) THEN
5226	! adeclarer
5227	#ifdef CPP_COSP
5228	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5229
5230	IF (prt_level .GE.10) THEN
5231	print*,'freq_cosp',freq_cosp
5232	ENDIF
5233	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
5234	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
5235	! s ref_liq,ref_ice
5236	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
5237	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
5238	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
5239	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
5240	JrNt,ref_liq,ref_ice, &
5241	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
5242	zu10m,zv10m,pphis, &
5243	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
5244	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
5245	prfl(:,1:klev),psfl(:,1:klev), &
5246	pmflxr(:,1:klev),pmflxs(:,1:klev), &
5247	mr_ozone,cldtau, cldemi)
5248
5249	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
5250	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
5251	! M clMISR,
5252	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
5253	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
5254
5255	ENDIF
5256	#endif
5257
5258	#ifdef CPP_COSP2
5259	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5260
5261	IF (prt_level .GE.10) THEN
5262	print*,'freq_cosp',freq_cosp
5263	ENDIF
5264	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
5265	print*,'Dans physiq.F avant appel '
5266	! s ref_liq,ref_ice
5267	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
5268	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
5269	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
5270	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
5271	JrNt,ref_liq,ref_ice, &
5272	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
5273	zu10m,zv10m,pphis, &
5274	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
5275	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
5276	prfl(:,1:klev),psfl(:,1:klev), &
5277	pmflxr(:,1:klev),pmflxs(:,1:klev), &
5278	mr_ozone,cldtau, cldemi)
5279	ENDIF
5280	#endif
5281
5282	#ifdef CPP_COSPV2
5283	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5284	! IF (MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
5285
5286	IF (prt_level .GE.10) THEN
5287	print*,'freq_cosp',freq_cosp
5288	ENDIF
5289	DO k = 1, klev
5290	DO i = 1, klon
5291	phicosp(i,k) = pphi(i,k) + pphis(i)
5292	ENDDO
5293	ENDDO
5294	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
5295	print*,'Dans physiq.F avant appel '
5296	! s ref_liq,ref_ice
5297	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
5298	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
5299	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
5300	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
5301	JrNt,ref_liq,ref_ice, &
5302	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
5303	zu10m,zv10m,pphis, &
5304	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
5305	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
5306	prfl(:,1:klev),psfl(:,1:klev), &
5307	pmflxr(:,1:klev),pmflxs(:,1:klev), &
5308	mr_ozone,cldtau, cldemi)
5309	ENDIF
5310	#endif
5311
5312	ENDIF !ok_cosp
5313
5314
5315	! Marine
5316
5317	IF (ok_airs) then
5318
5319	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
5320	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
5321	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
5322	& map_prop_hc,map_prop_hist,&
5323	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
5324	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
5325	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
5326	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
5327	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
5328	& map_ntot,map_hc,map_hist,&
5329	& map_Cb,map_ThCi,map_Anv,&
5330	& alt_tropo )
5331	ENDIF
5332
5333	ENDIF ! ok_airs
5334
5335
5336	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5337	!AA
5338	!AA Installation de l'interface online-offline pour traceurs
5339	!AA
5340	!====================================================================
5341	! Calcul des tendances traceurs
5342	!====================================================================
5343	!
5344
5345	IF (type_trac == 'repr') THEN
5346	!MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
5347	!MM dans Reprobus
5348	sh_in(:,:) = q_seri(:,:)
5349	#ifdef REPROBUS
5350	d_q_rep(:,:) = 0.
5351	d_ql_rep(:,:) = 0.
5352	d_qi_rep(:,:) = 0.
5353	#endif
5354	ELSE
5355	sh_in(:,:) = qx(:,:,ivap)
5356	IF (nqo >= 3) THEN
5357	ch_in(:,:) = qx(:,:,iliq) + qx(:,:,isol)
5358	ELSE
5359	ch_in(:,:) = qx(:,:,iliq)
5360	ENDIF
5361	ENDIF
5362
5363	#ifdef CPP_Dust
5364	! Avec SPLA, iflag_phytrac est forcé =1
5365	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
5366	pdtphys,ftsol, & ! I
5367	t,q_seri,paprs,pplay,RHcl, & ! I
5368	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
5369	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
5370	u_seri, v_seri, latitude_deg, longitude_deg, &
5371	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
5372	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
5373	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
5374	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
5375	fm_therm, entr_therm, rneb, & ! I
5376	beta_prec_fisrt,beta_prec, & !I
5377	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
5378	d_tr_dyn,tr_seri)
5379
5380	#else
5381	IF (iflag_phytrac == 1 ) THEN
5382	CALL phytrac ( &
5383	itap, days_elapsed+1, jH_cur, debut, &
5384	lafin, phys_tstep, u, v, t, &
5385	paprs, pplay, pmfu, pmfd, &
5386	pen_u, pde_u, pen_d, pde_d, &
5387	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
5388	u1, v1, ftsol, pctsrf, &
5389	zustar, zu10m, zv10m, &
5390	wstar(:,is_ave), ale_bl, ale_wake, &
5391	latitude_deg, longitude_deg, &
5392	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
5393	presnivs, pphis, pphi, albsol1, &
5394	sh_in, ch_in, rhcl, cldfra, rneb, &
5395	diafra, radocond, itop_con, ibas_con, &
5396	pmflxr, pmflxs, prfl, psfl, &
5397	da, phi, mp, upwd, &
5398	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
5399	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
5400	ev, ep, epmlmMm, eplaMm, & !<<RomP
5401	dnwd, aerosol_couple, flxmass_w, &
5402	tau_aero, piz_aero, cg_aero, ccm, &
5403	rfname, &
5404	d_tr_dyn, & !<<RomP
5405	tr_seri, init_source)
5406	#ifdef REPROBUS
5407
5408
5409	print*,'avt add phys rep',abortphy
5410
5411	CALL add_phys_tend &
5412	(du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,dqbs0,paprs,&
5413	'rep',abortphy,flag_inhib_tend,itap,0)
5414	IF (abortphy==1) Print*,'ERROR ABORT REP'
5415
5416	print*,'apr add phys rep',abortphy
5417
5418	#endif
5419	ENDIF ! (iflag_phytrac=1)
5420
5421	#endif
5422	!ENDIF ! (iflag_phytrac=1)
5423
5424	IF (offline) THEN
5425
5426	IF (prt_level.ge.9) &
5427	print*,'Attention on met a 0 les thermiques pour phystoke'
5428	CALL phystokenc ( &
5429	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
5430	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
5431	fm_therm,entr_therm, &
5432	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
5433	frac_impa, frac_nucl, &
5434	pphis,cell_area,phys_tstep,itap, &
5435	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
5436
5437
5438	ENDIF
5439
5440	!
5441	! Calculer le transport de l'eau et de l'energie (diagnostique)
5442	!
5443	CALL transp (paprs,zxtsol, t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
5444	ue, ve, uq, vq, uwat, vwat)
5445	!
5446	!IM global posePB BEG
5447	IF(1.EQ.0) THEN
5448	!
5449	CALL transp_lay (paprs,zxtsol, t_seri, q_seri, u_seri, v_seri, zphi, &
5450	ve_lay, vq_lay, ue_lay, uq_lay)
5451	!
5452	ENDIF !(1.EQ.0) THEN
5453	!IM global posePB END
5454	!
5455	! Accumuler les variables a stocker dans les fichiers histoire:
5456	!
5457
5458	!================================================================
5459	! Conversion of kinetic and potential energy into heat, for
5460	! parameterisation of subgrid-scale motions
5461	!================================================================
5462
5463	d_t_ec(:,:)=0.
5464	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
5465	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx,ivap,iliq,isol, &
5466	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
5467	zmasse,exner,d_t_ec)
5468	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
5469
5470	!=======================================================================
5471	! SORTIES
5472	!=======================================================================
5473	!
5474	!IM initialisation + calculs divers diag AMIP2
5475	!
5476	include "calcul_divers.h"
5477	!
5478	!IM Interpolation sur les niveaux de pression du NMC
5479	! -------------------------------------------------
5480	!
5481	include "calcul_STDlev.h"
5482	!
5483	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
5484	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
5485	!
5486	!cc prw = eau precipitable
5487	! prlw = colonne eau liquide
5488	! prlw = colonne eau solide
5489	! prbsw = colonne neige soufflee
5490	prw(:) = 0.
5491	prlw(:) = 0.
5492	prsw(:) = 0.
5493	prbsw(:) = 0.
5494	DO k = 1, klev
5495	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
5496	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
5497	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
5498	prbsw(:)= prbsw(:) + qbs_seri(:,k)*zmasse(:,k)
5499	ENDDO
5500	!
5501	IF (ANY(type_trac == ['inca','inco'])) THEN
5502	#ifdef INCA
5503	CALL VTe(VTphysiq)
5504	CALL VTb(VTinca)
5505
5506	CALL chemhook_end ( &
5507	phys_tstep, &
5508	pplay, &
5509	t_seri, &
5510	tr_seri(:,:,1+nqCO2:nbtr), &
5511	nbtr, &
5512	paprs, &
5513	q_seri, &
5514	cell_area, &
5515	pphi, &
5516	pphis, &
5517	zx_rh, &
5518	aps, bps, ap, bp, lafin)
5519
5520	CALL VTe(VTinca)
5521	CALL VTb(VTphysiq)
5522	#endif
5523	ENDIF
5524
5525	IF (type_trac == 'repr') THEN
5526	#ifdef REPROBUS
5527	CALL coord_hyb_rep(paprs, pplay, aps, bps, ap, bp, cell_area)
5528	#endif
5529	ENDIF
5530
5531	!
5532	! Convertir les incrementations en tendances
5533	!
5534	IF (prt_level .GE.10) THEN
5535	print *,'Convertir les incrementations en tendances '
5536	ENDIF
5537	!
5538	IF (mydebug) THEN
5539	CALL writefield_phy('u_seri',u_seri,nbp_lev)
5540	CALL writefield_phy('v_seri',v_seri,nbp_lev)
5541	CALL writefield_phy('t_seri',t_seri,nbp_lev)
5542	CALL writefield_phy('q_seri',q_seri,nbp_lev)
5543	ENDIF
5544
5545	DO k = 1, klev
5546	DO i = 1, klon
5547	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
5548	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
5549	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
5550	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
5551	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
5552	!CR: on ajoute le contenu en glace
5553	IF (nqo >= 3) THEN
5554	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
5555	ENDIF
5556	!--ice_sursat: nqo=4, on ajoute rneb
5557	IF (nqo.ge.4 .and. ok_ice_sursat) THEN
5558	d_qx(i,k,irneb) = ( rneb_seri(i,k) - qx(i,k,irneb) ) / phys_tstep
5559	ENDIF
5560
5561	IF (nqo.ge.4 .and. ok_bs) THEN
5562	d_qx(i,k,ibs) = ( qbs_seri(i,k) - qx(i,k,ibs) ) / phys_tstep
5563	ENDIF
5564
5565	ENDDO
5566	ENDDO
5567	!
5568	! DC: All iterations are cycled if nqtot==nqo, so no nqtot>nqo condition required
5569	itr = 0
5570	DO iq = 1, nqtot
5571	IF(.NOT.tracers(iq)%isInPhysics) CYCLE
5572	itr = itr+1
5573	DO k = 1, klev
5574	DO i = 1, klon
5575	d_qx(i,k,iq) = ( tr_seri(i,k,itr) - qx(i,k,iq) ) / phys_tstep
5576	ENDDO
5577	ENDDO
5578	ENDDO
5579	!
5580	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
5581	!IM global posePB include "write_bilKP_ins.h"
5582	!IM global posePB include "write_bilKP_ave.h"
5583	!
5584
5585	!--OB mass fixer
5586	!--profile is corrected to force mass conservation of water
5587	IF (mass_fixer) THEN
5588	qql2(:)=0.0
5589	DO k = 1, klev
5590	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
5591	ENDDO
5592
5593	#ifdef CPP_StratAer
5594	IF (ok_qemiss) THEN
5595	DO k = 1, klev
5596	qql1(:) = qql1(:)+d_q_emiss(:,k)*zmasse(:,k)
5597	ENDDO
5598	ENDIF
5599	#endif
5600	IF (ok_qch4) THEN
5601	DO k = 1, klev
5602	qql1(:) = qql1(:)+d_q_ch4_dtime(:,k)*zmasse(:,k)
5603	ENDDO
5604	ENDIF
5605
5606	DO i = 1, klon
5607	!--compute ratio of what q+ql should be with conservation to what it is
5608	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
5609	DO k = 1, klev
5610	q_seri(i,k) =q_seri(i,k)*corrqql
5611	ql_seri(i,k)=ql_seri(i,k)*corrqql
5612	ENDDO
5613	ENDDO
5614	ENDIF
5615	!--fin mass fixer
5616
5617	! Sauvegarder les valeurs de t et q a la fin de la physique:
5618	!
5619	u_ancien(:,:) = u_seri(:,:)
5620	v_ancien(:,:) = v_seri(:,:)
5621	t_ancien(:,:) = t_seri(:,:)
5622	q_ancien(:,:) = q_seri(:,:)
5623	ql_ancien(:,:) = ql_seri(:,:)
5624	qs_ancien(:,:) = qs_seri(:,:)
5625	qbs_ancien(:,:) = qbs_seri(:,:)
5626	rneb_ancien(:,:) = rneb_seri(:,:)
5627	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
5628	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
5629	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
5630	CALL water_int(klon,klev,qbs_ancien,zmasse,prbsw_ancien)
5631	! !! RomP >>>
5632	IF (nqtot > nqo) tr_ancien(:,:,:) = tr_seri(:,:,:)
5633	! !! RomP <<<
5634	!==========================================================================
5635	! Sorties des tendances pour un point particulier
5636	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
5637	! pour le debug
5638	! La valeur de igout est attribuee plus haut dans le programme
5639	!==========================================================================
5640
5641	IF (prt_level.ge.1) THEN
5642	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
5643	write(lunout,*) &
5644	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
5645	write(lunout,*) &
5646	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
5647	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
5648	pctsrf(igout,is_sic)
5649	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
5650	DO k=1,klev
5651	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
5652	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
5653	d_t_eva(igout,k)
5654	ENDDO
5655	write(lunout,*) 'cool,heat'
5656	DO k=1,klev
5657	write(lunout,*) cool(igout,k),heat(igout,k)
5658	ENDDO
5659
5660	!jyg< (En attendant de statuer sur le sort de d_t_oli)
5661	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5662	!jyg! do k=1,klev
5663	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
5664	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5665	!jyg! enddo
5666	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5667	DO k=1,klev
5668	write(lunout,*) d_t_vdf(igout,k), &
5669	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5670	ENDDO
5671	!>jyg
5672
5673	write(lunout,*) 'd_ps ',d_ps(igout)
5674	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
5675	DO k=1,klev
5676	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
5677	d_qx(igout,k,1),d_qx(igout,k,2)
5678	ENDDO
5679	ENDIF
5680
5681	!============================================================
5682	! Calcul de la temperature potentielle
5683	!============================================================
5684	DO k = 1, klev
5685	DO i = 1, klon
5686	!JYG/IM theta en debut du pas de temps
5687	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
5688	!JYG/IM theta en fin de pas de temps de physique
5689	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
5690	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
5691	! MPL 20130625
5692	! fth_fonctions.F90 et parkind1.F90
5693	! sinon thetal=theta
5694	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
5695	! : ql_seri(i,k))
5696	thetal(i,k)=theta(i,k)
5697	ENDDO
5698	ENDDO
5699	!
5700
5701	! 22.03.04 BEG
5702	!=============================================================
5703	! Ecriture des sorties
5704	!=============================================================
5705	#ifdef CPP_IOIPSL
5706
5707	! Recupere des varibles calcule dans differents modules
5708	! pour ecriture dans histxxx.nc
5709
5710	! Get some variables from module fonte_neige_mod
5711	CALL fonte_neige_get_vars(pctsrf, &
5712	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
5713
5714
5715	!=============================================================
5716	! Separation entre thermiques et non thermiques dans les sorties
5717	! de fisrtilp
5718	!=============================================================
5719
5720	IF (iflag_thermals>=1) THEN
5721	d_t_lscth=0.
5722	d_t_lscst=0.
5723	d_q_lscth=0.
5724	d_q_lscst=0.
5725	DO k=1,klev
5726	DO i=1,klon
5727	IF (ptconvth(i,k)) THEN
5728	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5729	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5730	ELSE
5731	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5732	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5733	ENDIF
5734	ENDDO
5735	ENDDO
5736
5737	DO i=1,klon
5738	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
5739	plul_th(i)=prfl(i,1)+psfl(i,1)
5740	ENDDO
5741	ENDIF
5742
5743	!On effectue les sorties:
5744
5745	#ifdef CPP_Dust
5746	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
5747	pplay, lmax_th, aerosol_couple, &
5748	ok_ade, ok_aie, ivap, ok_sync, &
5749	ptconv, read_climoz, clevSTD, &
5750	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
5751	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5752	#else
5753	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
5754	pplay, lmax_th, aerosol_couple, &
5755	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, ibs, &
5756	ok_sync, ptconv, read_climoz, clevSTD, &
5757	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
5758	flag_aerosol, flag_aerosol_strat, ok_cdnc,t, u1, v1)
5759	#endif
5760
5761	#ifndef CPP_XIOS
5762	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
5763	#endif
5764
5765	#endif
5766	! Petit appelle de sorties pour accompagner le travail sur phyex
5767	if ( iflag_physiq == 1 ) then
5768	call output_physiqex(debut,jD_eq,pdtphys,presnivs,paprs,u,v,t,qx,cldfra,0.t,0.t,0.*t,pbl_tke,theta)
5769	endif
5770
5771	!====================================================================
5772	! Arret du modele apres hgardfou en cas de detection d'un
5773	! plantage par hgardfou
5774	!====================================================================
5775
5776	IF (abortphy==1) THEN
5777	abort_message ='Plantage hgardfou'
5778	CALL abort_physic (modname,abort_message,1)
5779	ENDIF
5780
5781	! 22.03.04 END
5782	!
5783	!====================================================================
5784	! Si c'est la fin, il faut conserver l'etat de redemarrage
5785	!====================================================================
5786	!
5787
5788	! Disabling calls to the prt_alerte function
5789	alert_first_call = .FALSE.
5790
5791
5792	IF (lafin) THEN
5793	itau_phy = itau_phy + itap
5794	CALL phyredem ("restartphy.nc")
5795	! open(97,form="unformatted",file="finbin")
5796	! write(97) u_seri,v_seri,t_seri,q_seri
5797	! close(97)
5798
5799	IF (is_omp_master) THEN
5800
5801	IF (read_climoz >= 1) THEN
5802	IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5803	DEALLOCATE(press_edg_climoz)
5804	DEALLOCATE(press_cen_climoz)
5805	ENDIF
5806
5807	ENDIF
5808
5809	IF (using_xios) THEN
5810
5811	#ifdef INCA
5812	IF (type_trac == 'inca') THEN
5813	IF (is_omp_master .AND. grid_type==unstructured) THEN
5814	CALL finalize_inca
5815	ENDIF
5816	ENDIF
5817	#endif
5818
5819	IF (is_omp_master .and. grid_type==unstructured) CALL xios_context_finalize
5820	ENDIF
5821
5822	WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
5823
5824	ENDIF
5825
5826	! first=.false.
5827
5828	END SUBROUTINE physiq
5829
5830	END MODULE physiq_mod

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

Download in other formats:

Original Format