Context Navigation

physiq_mod.F90 @ 4804

Last change on this file since 4804 was 4803, checked in by evignon, 17 months ago
implementation sous flag des premiers changements concernant le traitement des precipitations grande echelle dans le cadre de l'atelier nuages Audran, Lea, Niels, Gwendal et Etienne
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 212.4 KB

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

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

Download in other formats:

Original Format