Context Navigation

physiq_mod.F90 @ 5306

Last change on this file since 5306 was 4758, checked in by idelkadi, 11 months ago
Update of Ecrad in LMDZ version 1.6.1 https://github.com/lguez/ecrad.git Continuation of work on the implementation of the double call of Ecrad in LMDZ
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.3 KB

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

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

Download in other formats:

Original Format