Context Navigation

physiq_mod.F90 @ 5229

Last change on this file since 5229 was 4485, checked in by evignon, 20 months ago
premier commit pour l'ajout de la neige soufflee sur la nouvelle branche
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: 196.2 KB

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

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

Download in other formats:

Original Format