Context Navigation

physiq_mod.F90 @ 5418

Last change on this file since 5418 was 4575, checked in by Laurent Fairhead, 18 months ago
Sortir les profils de flux sensible et d'eau dus à la diffusion turbulente (déja calculés auparavant) FC
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: 199.0 KB

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

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

Download in other formats:

Original Format