Context Navigation

physiq_mod.F90 @ 4722

Last change on this file since 4722 was 4722, checked in by Laurent Fairhead, 8 months ago
Modification by O. Torres to the cdrag routines to include different bulk formulae to calculate cdrag coefficients over ocean as well as an iteration of that calculation. The iteration is controlled by flag ok_cdrag_iter which if set to FALSE by default to converge with previous results. The choice of bulk formulae is set with the choix_bulk parameter The number of iterations to run is set with nit_bulk OT, PB, CD, LF
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 210.8 KB

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

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

Download in other formats:

Original Format