Context Navigation

physiq_mod.F90 @ 4870

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

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

Download in other formats:

Original Format