Context Navigation

physiq_mod.F90 @ 3865

Last change on this file since 3865 was 3865, checked in by lmdz-users, 3 years ago
Modifications from Thibaut to create an ESM with interactive CO2 + INCA aerosols
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: 188.8 KB

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

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

Download in other formats:

Original Format