Context Navigation

physiq_mod.F90 @ 4562

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

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

Download in other formats:

Original Format