Context Navigation

physiq_mod.F90 @ 5133

Last change on this file since 5133 was 5133, checked in by abarral, 5 months ago
Fix 1D, rrtm & ecrad compilation
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.0 KB

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

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

Download in other formats:

Original Format