Context Navigation

physiq_mod.F90 @ 4692

Last change on this file since 4692 was 4692, checked in by Laurent Fairhead, 13 months ago
Continuing work from poihl workshop on replayisation of routines: removed most of the modules (discussion needed for rest), declared calling arguments with INTENTs and created the ini routine
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 210.3 KB

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

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

Download in other formats:

Original Format