Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

physiq_mod.F90

Last change on this file was 4669, checked in by Laurent Fairhead, 15 months ago
Merged with trunk revision 4586 corresponding to june 2023 testing
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 199.5 KB

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

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

Download in other formats:

Original Format