Context Navigation

physiq_mod.F90 @ 2900

Last change on this file since 2900 was 2897, checked in by fhourdin, 7 years ago
Introduction d'une possible prise en compte de génération de TKE par les ondes de relief. Etienne Vignon et FH
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 `Author Date Id Revision`
File size: 160.9 KB

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

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

Download in other formats:

Original Format