Context Navigation

physiq_mod.F90 @ 5449

Last change on this file since 5449 was 5226, checked in by abarral, 4 months ago
Merge r5208 r5209
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 212.4 KB

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

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

Download in other formats:

Original Format