Context Navigation

physiq_mod.F90 @ 4768

Last change on this file since 4768 was 4755, checked in by lebasn, 8 months ago
Methox: Update comments to correct units + bugfix on output variable
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: 211.6 KB

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

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

Download in other formats:

Original Format