Context Navigation

physiq_mod.F90 @ 5054

Last change on this file since 5054 was 2879, checked in by Laurent Fairhead, 7 years ago
Inclusion of r2877 corrections in CM6.0.10-LR branch
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: 156.4 KB

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

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

Download in other formats:

Original Format