Context Navigation

physiq_mod.F90 @ 4547

Last change on this file since 4547 was 4547, checked in by fhourdin, 12 months ago
Details pour travail sur PHHYEX
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: 198.9 KB

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

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

Download in other formats:

Original Format