Context Navigation

physiq_mod.F90 @ 5194

Last change on this file since 5194 was 5194, checked in by abarral, 5 weeks ago
Merge r5189, r5191 from trunk
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: 210.8 KB

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

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

Download in other formats:

Original Format