Context Navigation

physiq_mod.F90 @ 5157

Last change on this file since 5157 was 5153, checked in by abarral, 7 weeks ago
Revert FCTTRE to INCLUDE to assess impact of inlining
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: 211.4 KB

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

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

Download in other formats:

Original Format