Context Navigation

physiq_mod.F90 @ 5097

Last change on this file since 5097 was 5091, checked in by abarral, 2 months ago
Move lmdz_netcdf_format.F90 -> lmdz_cppkeys_wrapper.F90 to handle other CPP keys Replace all (except wrapper) use of CPP_PHYS by fortran logical Refactor makelmdz_fcm (put blocks into functions, use modern bash)
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 212.6 KB

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

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

Download in other formats:

Original Format