Context Navigation

physiq_mod.F90 @ 4715

Last change on this file since 4715 was 4715, checked in by Laurent Fairhead, 8 months ago
Final (hopefully) commit from the newmicro replayisation workshop. The final USE statements that were still included in the cloud_optics_prop routine were moved to the call_cloud_optics_prop routine that sets up the interface between LMDZ and the parametrization. LF for LR, MCD, AI, EV, JBM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 210.3 KB

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

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

Download in other formats:

Original Format