Context Navigation

physiq_mod.F90 @ 2863

Last change on this file since 2863 was 2863, checked in by Ehouarn Millour, 7 years ago
Added missing declaration and closing bracket introduced in rev 2854. EM
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: 157.0 KB

Line
1	!
2	! $Id: physiq_mod.F90 2863 2017-04-24 08:14:02Z emillour $
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	drytausum_aero, &
95	!
96	!variables CFMIP2/CMIP5
97	topswad_aerop, solswad_aerop, &
98	topswai_aerop, solswai_aerop, &
99	topswad0_aerop, solswad0_aerop, &
100	topsw_aerop, topsw0_aerop, &
101	solsw_aerop, solsw0_aerop, &
102	topswcf_aerop, solswcf_aerop, &
103	!LW diagnostics
104	toplwad_aerop, sollwad_aerop, &
105	toplwai_aerop, sollwai_aerop, &
106	toplwad0_aerop, sollwad0_aerop, &
107	!
108	ptstar, pt0, slp, &
109	!
110	bils, &
111	!
112	cldh, cldl,cldm, cldq, cldt, &
113	JrNt, &
114	dthmin, evap, fder, plcl, plfc, &
115	prw, prlw, prsw, &
116	s_lcl, s_pblh, s_pblt, s_therm, &
117	cdragm, cdragh, &
118	zustar, zu10m, zv10m, rh2m, qsat2m, &
119	zq2m, zt2m, weak_inversion, &
120	zt2m_min_mon, zt2m_max_mon, & ! pour calcul_divers.h
121	t2m_min_mon, t2m_max_mon, & ! pour calcul_divers.h
122	!
123	s_pblh_x, s_pblh_w, &
124	s_lcl_x, s_lcl_w, &
125	!
126	slab_wfbils, tpot, tpote, &
127	ue, uq, ve, vq, zxffonte, &
128	zxfqcalving, zxfluxlat, &
129	zxrunofflic, &
130	zxtsol, snow_lsc, zxfqfonte, zxqsurf, &
131	rain_lsc, rain_num, &
132	!
133	sens_x, sens_w, &
134	zxfluxlat_x, zxfluxlat_w, &
135	!
136	dtvdf_x, dtvdf_w, &
137	dqvdf_x, dqvdf_w, &
138	pbl_tke_input, &
139	t_therm, q_therm, u_therm, v_therm, &
140	cdragh_x, cdragh_w, &
141	cdragm_x, cdragm_w, &
142	kh, kh_x, kh_w, &
143	!
144	wake_k, &
145	ale_wake, alp_wake, &
146	wake_h, wake_omg, &
147	! tendencies of delta T and delta q:
148	d_deltat_wk, d_deltaq_wk, & ! due to wakes
149	d_deltat_wk_gw, d_deltaq_wk_gw, & ! due to wake induced gravity waves
150	d_deltat_vdf, d_deltaq_vdf, & ! due to vertical diffusion
151	d_deltat_the, d_deltaq_the, & ! due to thermals
152	d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
153	! tendencies of wake fractional area and wake number per unit area:
154	d_s_wk, d_dens_wk, & ! due to wakes
155	!!! d_s_vdf, d_dens_vdf, & ! due to vertical diffusion
156	!!! d_s_the, d_dens_the, & ! due to thermals
157	!
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	logical ptconv(klon,klev)
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	INTEGER :: naero
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
1584	!--setting up dryaod_diag to TRUE in XIOS case
1585	DO naero = 1, naero_tot-1
1586	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
1587	ENDDO
1588	#endif
1589
1590	#ifdef CPP_IOIPSL
1591
1592	!$OMP MASTER
1593	! FH : if ok_sync=.true. , the time axis is written at each time step
1594	! in the output files. Only at the end in the opposite case
1595	ok_sync_omp=.false.
1596	CALL getin('ok_sync',ok_sync_omp)
1597	CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
1598	iGCM,jGCM,lonGCM,latGCM, &
1599	jjmp1,nlevSTD,clevSTD,rlevSTD, dtime,ok_veget, &
1600	type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1601	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1602	read_climoz, phys_out_filestations, &
1603	new_aod, aerosol_couple, &
1604	flag_aerosol_strat, pdtphys, paprs, pphis, &
1605	pplay, lmax_th, ptconv, ptconvth, ivap, &
1606	d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
1607	!$OMP END MASTER
1608	!$OMP BARRIER
1609	ok_sync=ok_sync_omp
1610
1611	freq_outNMC(1) = ecrit_files(7)
1612	freq_outNMC(2) = ecrit_files(8)
1613	freq_outNMC(3) = ecrit_files(9)
1614	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1615	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1616	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1617
1618	#ifndef CPP_XIOS
1619	CALL ini_paramLMDZ_phy(dtime,nid_ctesGCM)
1620	#endif
1621
1622	#endif
1623	ecrit_reg = ecrit_reg * un_jour
1624	ecrit_tra = ecrit_tra * un_jour
1625
1626	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
1627	date0 = jD_ref
1628	WRITE(,) 'physiq date0 : ',date0
1629	!
1630	!
1631	!
1632	! Prescrire l'ozone dans l'atmosphere
1633	!
1634	!
1635	!c DO i = 1, klon
1636	!c DO k = 1, klev
1637	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1638	!c ENDDO
1639	!c ENDDO
1640	!
1641	IF (type_trac == 'inca') THEN
1642	#ifdef INCA
1643	CALL VTe(VTphysiq)
1644	CALL VTb(VTinca)
1645	calday = REAL(days_elapsed) + jH_cur
1646	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1647
1648	CALL chemini( &
1649	rg, &
1650	ra, &
1651	cell_area, &
1652	latitude_deg, &
1653	longitude_deg, &
1654	presnivs, &
1655	calday, &
1656	klon, &
1657	nqtot, &
1658	nqo, &
1659	pdtphys, &
1660	annee_ref, &
1661	day_ref, &
1662	day_ini, &
1663	start_time, &
1664	itau_phy, &
1665	date0, &
1666	io_lon, &
1667	io_lat)
1668
1669	CALL VTe(VTinca)
1670	CALL VTb(VTphysiq)
1671	#endif
1672	ENDIF
1673	!
1674	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1675	! Nouvelle initialisation pour le rayonnement RRTM
1676	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1677
1678	CALL iniradia(klon,klev,paprs(1,1:klev+1))
1679
1680	!$omp single
1681	IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, &
1682	press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
1683	!$omp end single
1684	!
1685	!IM betaCRF
1686	pfree=70000. !Pa
1687	beta_pbl=1.
1688	beta_free=1.
1689	lon1_beta=-180.
1690	lon2_beta=+180.
1691	lat1_beta=90.
1692	lat2_beta=-90.
1693	mskocean_beta=.FALSE.
1694
1695	!albedo SB >>>
1696	select case(nsw)
1697	case(2)
1698	SFRWL(1)=0.45538747
1699	SFRWL(2)=0.54461211
1700	case(4)
1701	SFRWL(1)=0.45538747
1702	SFRWL(2)=0.32870591
1703	SFRWL(3)=0.18568763
1704	SFRWL(4)=3.02191470E-02
1705	case(6)
1706	SFRWL(1)=1.28432794E-03
1707	SFRWL(2)=0.12304168
1708	SFRWL(3)=0.33106142
1709	SFRWL(4)=0.32870591
1710	SFRWL(5)=0.18568763
1711	SFRWL(6)=3.02191470E-02
1712	end select
1713
1714
1715	!albedo SB <<<
1716
1717	OPEN(99,file='beta_crf.data',status='old', &
1718	form='formatted',err=9999)
1719	READ(99,*,end=9998) pfree
1720	READ(99,*,end=9998) beta_pbl
1721	READ(99,*,end=9998) beta_free
1722	READ(99,*,end=9998) lon1_beta
1723	READ(99,*,end=9998) lon2_beta
1724	READ(99,*,end=9998) lat1_beta
1725	READ(99,*,end=9998) lat2_beta
1726	READ(99,*,end=9998) mskocean_beta
1727	9998 Continue
1728	CLOSE(99)
1729	9999 Continue
1730	WRITE(,)'pfree=',pfree
1731	WRITE(,)'beta_pbl=',beta_pbl
1732	WRITE(,)'beta_free=',beta_free
1733	WRITE(,)'lon1_beta=',lon1_beta
1734	WRITE(,)'lon2_beta=',lon2_beta
1735	WRITE(,)'lat1_beta=',lat1_beta
1736	WRITE(,)'lat2_beta=',lat2_beta
1737	WRITE(,)'mskocean_beta=',mskocean_beta
1738	ENDIF
1739	!
1740	! ************** Fin de IF ( debut ) *************
1741	!
1742	!
1743	! Incrementer le compteur de la physique
1744	!
1745	itap = itap + 1
1746	IF (is_master .OR. prt_level > 9) THEN
1747	IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
1748	WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
1749	WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
1750	100 FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
1751	ENDIF
1752	ENDIF
1753	!
1754	!
1755	! Update fraction of the sub-surfaces (pctsrf) and
1756	! initialize, where a new fraction has appeared, all variables depending
1757	! on the surface fraction.
1758	!
1759	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
1760	pctsrf, fevap, z0m, z0h, agesno, &
1761	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
1762
1763	! Update time and other variables in Reprobus
1764	IF (type_trac == 'repr') THEN
1765	#ifdef REPROBUS
1766	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1767	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1768	CALL Rtime(debut)
1769	#endif
1770	ENDIF
1771
1772
1773	! Tendances bidons pour les processus qui n'affectent pas certaines
1774	! variables.
1775	du0(:,:)=0.
1776	dv0(:,:)=0.
1777	dt0 = 0.
1778	dq0(:,:)=0.
1779	dql0(:,:)=0.
1780	dqi0(:,:)=0.
1781	dsig0(:) = 0.
1782	ddens0(:) = 0.
1783	wkoccur1(:)=1
1784	!
1785	! Mettre a zero des variables de sortie (pour securite)
1786	!
1787	DO i = 1, klon
1788	d_ps(i) = 0.0
1789	ENDDO
1790	DO k = 1, klev
1791	DO i = 1, klon
1792	d_t(i,k) = 0.0
1793	d_u(i,k) = 0.0
1794	d_v(i,k) = 0.0
1795	ENDDO
1796	ENDDO
1797	DO iq = 1, nqtot
1798	DO k = 1, klev
1799	DO i = 1, klon
1800	d_qx(i,k,iq) = 0.0
1801	ENDDO
1802	ENDDO
1803	ENDDO
1804	da(:,:)=0.
1805	mp(:,:)=0.
1806	phi(:,:,:)=0.
1807	! RomP >>>
1808	phi2(:,:,:)=0.
1809	beta_prec_fisrt(:,:)=0.
1810	beta_prec(:,:)=0.
1811	epmlmMm(:,:,:)=0.
1812	eplaMm(:,:)=0.
1813	d1a(:,:)=0.
1814	dam(:,:)=0.
1815	pmflxr=0.
1816	pmflxs=0.
1817	! RomP <<<
1818
1819	!
1820	! Ne pas affecter les valeurs entrees de u, v, h, et q
1821	!
1822	DO k = 1, klev
1823	DO i = 1, klon
1824	t_seri(i,k) = t(i,k)
1825	u_seri(i,k) = u(i,k)
1826	v_seri(i,k) = v(i,k)
1827	q_seri(i,k) = qx(i,k,ivap)
1828	ql_seri(i,k) = qx(i,k,iliq)
1829	!CR: ATTENTION, on rajoute la variable glace
1830	IF (nqo.eq.2) THEN
1831	qs_seri(i,k) = 0.
1832	ELSE IF (nqo.eq.3) THEN
1833	qs_seri(i,k) = qx(i,k,isol)
1834	ENDIF
1835	ENDDO
1836	ENDDO
1837	!
1838	!--OB mass fixer
1839	IF (mass_fixer) THEN
1840	!--store initial water burden
1841	qql1(:)=0.0
1842	DO k = 1, klev
1843	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
1844	ENDDO
1845	ENDIF
1846	!--fin mass fixer
1847
1848	tke0(:,:)=pbl_tke(:,:,is_ave)
1849	!CR:Nombre de traceurs de l'eau: nqo
1850	! IF (nqtot.GE.3) THEN
1851	IF (nqtot.GE.(nqo+1)) THEN
1852	! DO iq = 3, nqtot
1853	DO iq = nqo+1, nqtot
1854	DO k = 1, klev
1855	DO i = 1, klon
1856	! tr_seri(i,k,iq-2) = qx(i,k,iq)
1857	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
1858	ENDDO
1859	ENDDO
1860	ENDDO
1861	ELSE
1862	DO k = 1, klev
1863	DO i = 1, klon
1864	tr_seri(i,k,1) = 0.0
1865	ENDDO
1866	ENDDO
1867	ENDIF
1868	!
1869	DO i = 1, klon
1870	ztsol(i) = 0.
1871	ENDDO
1872	DO nsrf = 1, nbsrf
1873	DO i = 1, klon
1874	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1875	ENDDO
1876	ENDDO
1877	! Initialize variables used for diagnostic purpose
1878	IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
1879
1880	! Diagnostiquer la tendance dynamique
1881	!
1882	IF (ancien_ok) THEN
1883	!
1884	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/dtime
1885	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/dtime
1886	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/dtime
1887	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/dtime
1888	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/dtime
1889	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/dtime
1890	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
1891	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/dtime
1892	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
1893	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/dtime
1894	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
1895	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/dtime
1896	! !! RomP >>> td dyn traceur
1897	IF (nqtot.GT.nqo) THEN ! jyg
1898	DO iq = nqo+1, nqtot ! jyg
1899	d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/dtime ! jyg
1900	ENDDO
1901	ENDIF
1902	! !! RomP <<<
1903	ELSE
1904	d_u_dyn(:,:) = 0.0
1905	d_v_dyn(:,:) = 0.0
1906	d_t_dyn(:,:) = 0.0
1907	d_q_dyn(:,:) = 0.0
1908	d_ql_dyn(:,:) = 0.0
1909	d_qs_dyn(:,:) = 0.0
1910	d_q_dyn2d(:) = 0.0
1911	d_ql_dyn2d(:) = 0.0
1912	d_qs_dyn2d(:) = 0.0
1913	! !! RomP >>> td dyn traceur
1914	IF (nqtot.GT.nqo) THEN ! jyg
1915	DO iq = nqo+1, nqtot ! jyg
1916	d_tr_dyn(:,:,iq-nqo)= 0.0 ! jyg
1917	ENDDO
1918	ENDIF
1919	! !! RomP <<<
1920	ancien_ok = .TRUE.
1921	ENDIF
1922	!
1923	! Ajouter le geopotentiel du sol:
1924	!
1925	DO k = 1, klev
1926	DO i = 1, klon
1927	zphi(i,k) = pphi(i,k) + pphis(i)
1928	ENDDO
1929	ENDDO
1930	!
1931	! Verifier les temperatures
1932	!
1933	!IM BEG
1934	IF (check) THEN
1935	amn=MIN(ftsol(1,is_ter),1000.)
1936	amx=MAX(ftsol(1,is_ter),-1000.)
1937	DO i=2, klon
1938	amn=MIN(ftsol(i,is_ter),amn)
1939	amx=MAX(ftsol(i,is_ter),amx)
1940	ENDDO
1941	!
1942	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1943	ENDIF !(check) THEN
1944	!IM END
1945	!
1946	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
1947	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
1948
1949	!
1950	!IM BEG
1951	IF (check) THEN
1952	amn=MIN(ftsol(1,is_ter),1000.)
1953	amx=MAX(ftsol(1,is_ter),-1000.)
1954	DO i=2, klon
1955	amn=MIN(ftsol(i,is_ter),amn)
1956	amx=MAX(ftsol(i,is_ter),amx)
1957	ENDDO
1958	!
1959	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1960	ENDIF !(check) THEN
1961	!IM END
1962	!
1963	! Mettre en action les conditions aux limites (albedo, sst, etc.).
1964	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
1965	!
1966	! Update ozone if day change
1967	IF (MOD(itap-1,lmt_pas) == 0) THEN
1968	IF (read_climoz <= 0) THEN
1969	! Once per day, update ozone from Royer:
1970	IF (solarlong0<-999.) then
1971	! Generic case with evolvoing season
1972	zzz=real(days_elapsed+1)
1973	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
1974	! Particular case with annual mean insolation
1975	zzz=real(90) ! could be revisited
1976	IF (read_climoz/=-1) THEN
1977	abort_message ='read_climoz=-1 is recommended when ' &
1978	// 'solarlong0=1000.'
1979	CALL abort_physic (modname,abort_message,1)
1980	ENDIF
1981	ELSE
1982	! Case where the season is imposed with solarlong0
1983	zzz=real(90) ! could be revisited
1984	ENDIF
1985
1986	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
1987	ELSE
1988	!--- ro3i = elapsed days number since current year 1st january, 0h
1989	ro3i=days_elapsed+jh_cur-jh_1jan
1990	!--- scaling for old style files (360 records)
1991	IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
1992	IF(adjust_tropopause) THEN
1993	CALL dyn_tropopause(t_seri, ztsol, paprs, pplay, presnivs, rot, &
1994	pr_tropopause)
1995	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
1996	ro3i, press_edg_climoz, paprs, wo, time_climoz, &
1997	latitude_deg, press_cen_climoz, pr_tropopause)
1998	ELSE
1999	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2000	ro3i, press_edg_climoz, paprs, wo, time_climoz)
2001	END IF
2002	! Convert from mole fraction of ozone to column density of ozone in a
2003	! cell, in kDU:
2004	FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2005	* zmasse / dobson_u / 1e3
2006	! (By regridding ozone values for LMDZ only once a day, we
2007	! have already neglected the variation of pressure in one
2008	! day. So do not recompute "wo" at each time step even if
2009	! "zmasse" changes a little.)
2010	ENDIF
2011	ENDIF
2012	!
2013	! Re-evaporer l'eau liquide nuageuse
2014	!
2015	CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2016	& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
2017
2018	CALL add_phys_tend &
2019	(du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
2020	'eva',abortphy,flag_inhib_tend,itap,0)
2021	call prt_enerbil('eva',itap)
2022
2023	!=========================================================================
2024	! Calculs de l'orbite.
2025	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2026	! doit donc etre plac\'e avant radlwsw et pbl_surface
2027
2028	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2029	CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2030	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2031	!
2032	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2033	! solarlong0
2034	IF (solarlong0<-999.) THEN
2035	IF (new_orbit) THEN
2036	! calcul selon la routine utilisee pour les planetes
2037	CALL solarlong(day_since_equinox, zlongi, dist)
2038	ELSE
2039	! calcul selon la routine utilisee pour l'AR4
2040	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2041	ENDIF
2042	ELSE
2043	zlongi=solarlong0 ! longitude solaire vraie
2044	dist=1. ! distance au soleil / moyenne
2045	ENDIF
2046
2047	IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2048
2049
2050	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2051	! Calcul de l'ensoleillement :
2052	! ============================
2053	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2054	! l'annee a partir d'une formule analytique.
2055	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2056	! non nul aux poles.
2057	IF (abs(solarlong0-1000.)<1.e-4) THEN
2058	CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
2059	latitude_deg,longitude_deg,rmu0,fract)
2060	JrNt = 1.0
2061	ELSE
2062	! recode par Olivier Boucher en sept 2015
2063	SELECT CASE (iflag_cycle_diurne)
2064	CASE(0)
2065	! Sans cycle diurne
2066	CALL angle(zlongi, latitude_deg, fract, rmu0)
2067	swradcorr = 1.0
2068	JrNt = 1.0
2069	zrmu0 = rmu0
2070	CASE(1)
2071	! Avec cycle diurne sans application des poids
2072	! bit comparable a l ancienne formulation cycle_diurne=true
2073	! on integre entre gmtime et gmtime+radpas
2074	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
2075	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2076	latitude_deg,longitude_deg,rmu0,fract)
2077	zrmu0 = rmu0
2078	swradcorr = 1.0
2079	! Calcul du flag jour-nuit
2080	JrNt = 0.0
2081	WHERE (fract.GT.0.0) JrNt = 1.0
2082	CASE(2)
2083	! Avec cycle diurne sans application des poids
2084	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2085	! Comme cette routine est appele a tous les pas de temps de
2086	! la physique meme si le rayonnement n'est pas appele je
2087	! remonte en arriere les radpas-1 pas de temps
2088	! suivant. Petite ruse avec MOD pour prendre en compte le
2089	! premier pas de temps de la physique pendant lequel
2090	! itaprad=0
2091	zdtime1=dtime*REAL(-MOD(itaprad,radpas)-1)
2092	zdtime2=dtime*REAL(radpas-MOD(itaprad,radpas)-1)
2093	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2094	latitude_deg,longitude_deg,rmu0,fract)
2095	!
2096	! Calcul des poids
2097	!
2098	zdtime1=-dtime !--on corrige le rayonnement pour representer le
2099	zdtime2=0.0 !--pas de temps de la physique qui se termine
2100	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2101	latitude_deg,longitude_deg,zrmu0,zfract)
2102	swradcorr = 0.0
2103	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2104	swradcorr=zfract/fract*zrmu0/rmu0
2105	! Calcul du flag jour-nuit
2106	JrNt = 0.0
2107	WHERE (zfract.GT.0.0) JrNt = 1.0
2108	END SELECT
2109	ENDIF
2110
2111	IF (mydebug) THEN
2112	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2113	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2114	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2115	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2116	ENDIF
2117
2118	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2119	! Appel au pbl_surface : Planetary Boudary Layer et Surface
2120	! Cela implique tous les interactions des sous-surfaces et la
2121	! partie diffusion turbulent du couche limit.
2122	!
2123	! Certains varibales de sorties de pbl_surface sont utiliser que pour
2124	! ecriture des fihiers hist_XXXX.nc, ces sont :
2125	! qsol, zq2m, s_pblh, s_lcl,
2126	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2127	! s_therm, s_trmb1, s_trmb2, s_trmb3,
2128	! zu10m, zv10m, fder,
2129	! zxqsurf, rh2m, zxfluxu, zxfluxv,
2130	! frugs, agesno, fsollw, fsolsw,
2131	! d_ts, fevap, fluxlat, t2m,
2132	! wfbils, wfbilo, fluxt, fluxu, fluxv,
2133	!
2134	! Certains ne sont pas utiliser du tout :
2135	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2136	!
2137
2138	! Calcul de l'humidite de saturation au niveau du sol
2139
2140
2141
2142	IF (iflag_pbl/=0) THEN
2143
2144	!jyg+nrlmd<
2145	!!jyg IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2146	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
2147	print *,'debut du splitting de la PBL'
2148	ENDIF
2149	! !!
2150	!>jyg+nrlmd
2151	!
2152	!-------gustiness calculation-------!
2153	IF (iflag_gusts==0) THEN
2154	gustiness(1:klon)=0
2155	ELSE IF (iflag_gusts==1) THEN
2156	do i = 1, klon
2157	gustiness(i)=f_gust_blale_bl(i)+f_gust_wkale_wake(i)
2158	enddo
2159	! ELSE IF (iflag_gusts==2) THEN
2160	! do i = 1, klon
2161	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
2162	! *ale_wake(i) !! need to make sigma_wk accessible here
2163	! enddo
2164	! ELSE IF (iflag_gusts==3) THEN
2165	! do i = 1, klon
2166	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
2167	! enddo
2168	ENDIF
2169
2170
2171
2172	CALL pbl_surface( &
2173	dtime, date0, itap, days_elapsed+1, &
2174	debut, lafin, &
2175	longitude_deg, latitude_deg, rugoro, zrmu0, &
2176	zsig, sollwdown, pphi, cldt, &
2177	rain_fall, snow_fall, solsw, sollw, &
2178	gustiness, &
2179	t_seri, q_seri, u_seri, v_seri, &
2180	!nrlmd+jyg<
2181	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2182	!>nrlmd+jyg
2183	pplay, paprs, pctsrf, &
2184	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2185	!albedo SB <<<
2186	cdragh, cdragm, u1, v1, &
2187	!albedo SB >>>
2188	! albsol1, albsol2, sens, evap, &
2189	albsol_dir, albsol_dif, sens, evap, &
2190	!albedo SB <<<
2191	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
2192	zxtsol, zxfluxlat, zt2m, qsat2m, &
2193	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2194	!nrlmd<
2195	!jyg<
2196	d_t_vdf_w, d_q_vdf_w, &
2197	d_t_vdf_x, d_q_vdf_x, &
2198	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2199	!>jyg
2200	delta_tsurf,wake_dens, &
2201	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2202	kh,kh_x,kh_w, &
2203	!>nrlmd
2204	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2205	slab_wfbils, &
2206	qsol, zq2m, s_pblh, s_lcl, &
2207	!jyg<
2208	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2209	!>jyg
2210	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2211	s_therm, s_trmb1, s_trmb2, s_trmb3, &
2212	zustar, zu10m, zv10m, fder, &
2213	zxqsurf, rh2m, zxfluxu, zxfluxv, &
2214	z0m, z0h, agesno, fsollw, fsolsw, &
2215	d_ts, fevap, fluxlat, t2m, &
2216	wfbils, wfbilo, wfevap, wfrain, wfsnow, &
2217	fluxt, fluxu, fluxv, &
2218	dsens, devap, zxsnow, &
2219	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
2220	!nrlmd+jyg<
2221	wake_delta_pbl_TKE &
2222	!>nrlmd+jyg
2223	)
2224	!
2225	! Add turbulent diffusion tendency to the wake difference variables
2226	!!jyg IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2227	IF (mod(iflag_pbl_split,10) .NE. 0) THEN
2228	!jyg<
2229	d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2230	d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2231	CALL add_wake_tend &
2232	(d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, wkoccur1, 'vdf', abortphy)
2233	ELSE
2234	d_deltat_vdf(:,:) = 0.
2235	d_deltaq_vdf(:,:) = 0.
2236	!>jyg
2237	ENDIF
2238
2239
2240	!---------------------------------------------------------------------
2241	! ajout des tendances de la diffusion turbulente
2242	IF (klon_glo==1) THEN
2243	CALL add_pbl_tend &
2244	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2245	'vdf',abortphy,flag_inhib_tend,itap)
2246	ELSE
2247	CALL add_phys_tend &
2248	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2249	'vdf',abortphy,flag_inhib_tend,itap,0)
2250	ENDIF
2251	call prt_enerbil('vdf',itap)
2252	!--------------------------------------------------------------------
2253
2254	IF (mydebug) THEN
2255	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2256	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2257	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2258	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2259	ENDIF
2260
2261	!albedo SB >>>
2262	albsol1=0.
2263	albsol2=0.
2264	falb1=0.
2265	falb2=0.
2266	SELECT CASE(nsw)
2267	CASE(2)
2268	albsol1=albsol_dir(:,1)
2269	albsol2=albsol_dir(:,2)
2270	falb1=falb_dir(:,1,:)
2271	falb2=falb_dir(:,2,:)
2272	CASE(4)
2273	albsol1=albsol_dir(:,1)
2274	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3) &
2275	+albsol_dir(:,4)*SFRWL(4)
2276	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2277	falb1=falb_dir(:,1,:)
2278	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3) &
2279	+falb_dir(:,4,:)*SFRWL(4)
2280	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2281	CASE(6)
2282	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2) &
2283	+albsol_dir(:,3)*SFRWL(3)
2284	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2285	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5) &
2286	+albsol_dir(:,6)*SFRWL(6)
2287	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2288	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2) &
2289	+falb_dir(:,3,:)*SFRWL(3)
2290	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2291	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5) &
2292	+falb_dir(:,6,:)*SFRWL(6)
2293	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2294	END SELECt
2295	!albedo SB <<<
2296
2297
2298	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2299	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2300
2301	ENDIF
2302	! =================================================================== c
2303	! Calcul de Qsat
2304
2305	DO k = 1, klev
2306	DO i = 1, klon
2307	zx_t = t_seri(i,k)
2308	IF (thermcep) THEN
2309	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2310	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2311	zx_qs = MIN(0.5,zx_qs)
2312	zcor = 1./(1.-retv*zx_qs)
2313	zx_qs = zx_qs*zcor
2314	ELSE
2315	!! IF (zx_t.LT.t_coup) THEN !jyg
2316	IF (zx_t.LT.rtt) THEN !jyg
2317	zx_qs = qsats(zx_t)/pplay(i,k)
2318	ELSE
2319	zx_qs = qsatl(zx_t)/pplay(i,k)
2320	ENDIF
2321	ENDIF
2322	zqsat(i,k)=zx_qs
2323	ENDDO
2324	ENDDO
2325
2326	IF (prt_level.ge.1) THEN
2327	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2328	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2329	ENDIF
2330	!
2331	! Appeler la convection (au choix)
2332	!
2333	DO k = 1, klev
2334	DO i = 1, klon
2335	conv_q(i,k) = d_q_dyn(i,k) &
2336	+ d_q_vdf(i,k)/dtime
2337	conv_t(i,k) = d_t_dyn(i,k) &
2338	+ d_t_vdf(i,k)/dtime
2339	ENDDO
2340	ENDDO
2341	IF (check) THEN
2342	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2343	WRITE(lunout,*) "avantcon=", za
2344	ENDIF
2345	zx_ajustq = .FALSE.
2346	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2347	IF (zx_ajustq) THEN
2348	DO i = 1, klon
2349	z_avant(i) = 0.0
2350	ENDDO
2351	DO k = 1, klev
2352	DO i = 1, klon
2353	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2354	*(paprs(i,k)-paprs(i,k+1))/RG
2355	ENDDO
2356	ENDDO
2357	ENDIF
2358
2359	! Calcule de vitesse verticale a partir de flux de masse verticale
2360	DO k = 1, klev
2361	DO i = 1, klon
2362	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
2363	ENDDO
2364	ENDDO
2365
2366	IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2367	omega(igout, :)
2368	!
2369	! Appel de la convection tous les "cvpas"
2370	!
2371	IF (MOD(itapcv,cvpas).EQ.0) THEN
2372
2373	IF (iflag_con.EQ.1) THEN
2374	abort_message ='reactiver le call conlmd dans physiq.F'
2375	CALL abort_physic (modname,abort_message,1)
2376	! CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2377	! . d_t_con, d_q_con,
2378	! . rain_con, snow_con, ibas_con, itop_con)
2379	ELSE IF (iflag_con.EQ.2) THEN
2380	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
2381	conv_t, conv_q, -evap, omega, &
2382	d_t_con, d_q_con, rain_con, snow_con, &
2383	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2384	kcbot, kctop, kdtop, pmflxr, pmflxs)
2385	d_u_con = 0.
2386	d_v_con = 0.
2387
2388	WHERE (rain_con < 0.) rain_con = 0.
2389	WHERE (snow_con < 0.) snow_con = 0.
2390	DO i = 1, klon
2391	ibas_con(i) = klev+1 - kcbot(i)
2392	itop_con(i) = klev+1 - kctop(i)
2393	ENDDO
2394	ELSE IF (iflag_con.GE.3) THEN
2395	! nb of tracers for the KE convection:
2396	! MAF la partie traceurs est faite dans phytrac
2397	! on met ntra=1 pour limiter les appels mais on peut
2398	! supprimer les calculs / ftra.
2399	ntra = 1
2400
2401	!=======================================================================
2402	!ajout pour la parametrisation des poches froides: calcul de
2403	!t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
2404	IF (iflag_wake>=1) THEN
2405	DO k=1,klev
2406	DO i=1,klon
2407	t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
2408	q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
2409	t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2410	q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2411	ENDDO
2412	ENDDO
2413	ELSE
2414	t_w(:,:) = t_seri(:,:)
2415	q_w(:,:) = q_seri(:,:)
2416	t_x(:,:) = t_seri(:,:)
2417	q_x(:,:) = q_seri(:,:)
2418	ENDIF
2419	!
2420	!jyg<
2421	! Perform dry adiabatic adjustment on wake profile
2422	! The corresponding tendencies are added to the convective tendencies
2423	! after the call to the convective scheme.
2424	IF (iflag_wake>=1) then
2425	IF (ok_adjwk) THEN
2426	limbas(:) = 1
2427	CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
2428	d_t_adjwk, d_q_adjwk)
2429	!
2430	DO k=1,klev
2431	DO i=1,klon
2432	IF (wake_s(i) .GT. 1.e-3) THEN
2433	t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
2434	q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
2435	d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
2436	d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
2437	ELSE
2438	d_deltat_ajs_cv(i,k) = 0.
2439	d_deltaq_ajs_cv(i,k) = 0.
2440	ENDIF
2441	ENDDO
2442	ENDDO
2443	CALL add_wake_tend &
2444	(d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, wkoccur1, 'ajs_cv', abortphy)
2445	ENDIF ! (ok_adjwk)
2446	ENDIF ! (iflag_wake>=1)
2447	!>jyg
2448	!
2449
2450	!! print *,'physiq. q_w(1,k), q_x(1,k) ', &
2451	!! (k, q_w(1,k), q_x(1,k),k=1,25)
2452
2453	!jyg<
2454	CALL alpale( debut, itap, dtime, paprs, omega, t_seri, &
2455	alp_offset, it_wape_prescr, wape_prescr, fip_prescr, &
2456	ale_bl_prescr, alp_bl_prescr, &
2457	wake_pe, wake_fip, &
2458	Ale_bl, Ale_bl_trig, Alp_bl, &
2459	Ale, Alp , Ale_wake, Alp_wake)
2460	!>jyg
2461	!
2462	! sb, oct02:
2463	! Schema de convection modularise et vectorise:
2464	! (driver commun aux versions 3 et 4)
2465	!
2466	IF (ok_cvl) THEN ! new driver for convectL
2467	!
2468	!jyg<
2469	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2470	! Calculate the upmost level of deep convection loops: k_upper_cv
2471	! (near 22 km)
2472	izero = klon/2+1/klon
2473	k_upper_cv = klev
2474	DO k = klev,1,-1
2475	IF (pphi(izero,k) > 22.e4) k_upper_cv = k
2476	ENDDO
2477	IF (prt_level .ge. 5) THEN
2478	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
2479	k_upper_cv
2480	ENDIF
2481	!
2482	!>jyg
2483	IF (type_trac == 'repr') THEN
2484	nbtr_tmp=ntra
2485	ELSE
2486	nbtr_tmp=nbtr
2487	ENDIF
2488	!jyg iflag_con est dans clesphys
2489	!c CALL concvl (iflag_con,iflag_clos,
2490	CALL concvl (iflag_clos, &
2491	dtime, paprs, pplay, k_upper_cv, t_x,q_x, &
2492	t_w,q_w,wake_s, &
2493	u_seri,v_seri,tr_seri,nbtr_tmp, &
2494	ALE,ALP, &
2495	sig1,w01, &
2496	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2497	rain_con, snow_con, ibas_con, itop_con, sigd, &
2498	ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
2499	Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2500	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2501	! RomP >>>
2502	!! . pmflxr,pmflxs,da,phi,mp,
2503	!! . ftd,fqd,lalim_conv,wght_th)
2504	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2505	ftd,fqd,lalim_conv,wght_th, &
2506	ev, ep,epmlmMm,eplaMm, &
2507	wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2508	tau_cld_cv,coefw_cld_cv,epmax_diag)
2509
2510	! RomP <<<
2511
2512	!IM begin
2513	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2514	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2515	!IM end
2516	!IM cf. FH
2517	clwcon0=qcondc
2518	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2519
2520	DO i = 1, klon
2521	IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2522	ENDDO
2523	!
2524	!jyg<
2525	! Add the tendency due to the dry adjustment of the wake profile
2526	IF (iflag_wake>=1) THEN
2527	DO k=1,klev
2528	DO i=1,klon
2529	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
2530	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
2531	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
2532	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
2533	ENDDO
2534	ENDDO
2535	ENDIF
2536	!>jyg
2537	!
2538	ELSE ! ok_cvl
2539
2540	! MAF conema3 ne contient pas les traceurs
2541	CALL conema3 (dtime, &
2542	paprs,pplay,t_seri,q_seri, &
2543	u_seri,v_seri,tr_seri,ntra, &
2544	sig1,w01, &
2545	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2546	rain_con, snow_con, ibas_con, itop_con, &
2547	upwd,dnwd,dnwd0,bas,top, &
2548	Ma,cape,tvp,rflag, &
2549	pbase &
2550	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2551	,clwcon0)
2552
2553	ENDIF ! ok_cvl
2554
2555	!
2556	! Correction precip
2557	rain_con = rain_con * cvl_corr
2558	snow_con = snow_con * cvl_corr
2559	!
2560
2561	IF (.NOT. ok_gust) THEN
2562	do i = 1, klon
2563	wd(i)=0.0
2564	enddo
2565	ENDIF
2566
2567	! =================================================================== c
2568	! Calcul des proprietes des nuages convectifs
2569	!
2570
2571	! calcul des proprietes des nuages convectifs
2572	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2573	IF (iflag_cld_cv == 0) THEN
2574	CALL clouds_gno &
2575	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2576	ELSE
2577	CALL clouds_bigauss &
2578	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
2579	ENDIF
2580
2581
2582	! =================================================================== c
2583
2584	DO i = 1, klon
2585	itop_con(i) = min(max(itop_con(i),1),klev)
2586	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2587	ENDDO
2588
2589	DO i = 1, klon
2590	ema_pcb(i) = paprs(i,ibas_con(i))
2591	ENDDO
2592	DO i = 1, klon
2593	! L'idicage de itop_con peut cacher un pb potentiel
2594	! FH sous la dictee de JYG, CR
2595	ema_pct(i) = paprs(i,itop_con(i)+1)
2596
2597	IF (itop_con(i).gt.klev-3) THEN
2598	IF (prt_level >= 9) THEN
2599	write(lunout,*)'La convection monte trop haut '
2600	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2601	ENDIF
2602	ENDIF
2603	ENDDO
2604	ELSE IF (iflag_con.eq.0) THEN
2605	write(lunout,*) 'On n appelle pas la convection'
2606	clwcon0=0.
2607	rnebcon0=0.
2608	d_t_con=0.
2609	d_q_con=0.
2610	d_u_con=0.
2611	d_v_con=0.
2612	rain_con=0.
2613	snow_con=0.
2614	bas=1
2615	top=1
2616	ELSE
2617	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2618	CALL abort_physic("physiq", "", 1)
2619	ENDIF
2620
2621	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2622	! . d_u_con, d_v_con)
2623
2624	!jyg Reinitialize proba_notrig and itapcv when convection has been called
2625	proba_notrig(:) = 1.
2626	itapcv = 0
2627	ENDIF ! (MOD(itapcv,cvpas).EQ.0)
2628	!
2629	itapcv = itapcv+1
2630
2631	!!!jyg Appel diagnostique a add_phys_tend pour tester la conservation de
2632	!!! l'energie dans les courants satures.
2633	!! d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
2634	!! d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
2635	!! dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
2636	!! CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat, &
2637	!! dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
2638	!! itap, 1)
2639	!! call prt_enerbil('convection_sat',itap)
2640	!!
2641	!!
2642	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
2643	'convection',abortphy,flag_inhib_tend,itap,0)
2644	call prt_enerbil('convection',itap)
2645
2646	!-------------------------------------------------------------------------
2647
2648	IF (mydebug) THEN
2649	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2650	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2651	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2652	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2653	ENDIF
2654
2655	IF (check) THEN
2656	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2657	WRITE(lunout,*)"aprescon=", za
2658	zx_t = 0.0
2659	za = 0.0
2660	DO i = 1, klon
2661	za = za + cell_area(i)/REAL(klon)
2662	zx_t = zx_t + (rain_con(i)+ &
2663	snow_con(i))*cell_area(i)/REAL(klon)
2664	ENDDO
2665	zx_t = zx_t/za*dtime
2666	WRITE(lunout,*)"Precip=", zx_t
2667	ENDIF
2668	IF (zx_ajustq) THEN
2669	DO i = 1, klon
2670	z_apres(i) = 0.0
2671	ENDDO
2672	DO k = 1, klev
2673	DO i = 1, klon
2674	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2675	*(paprs(i,k)-paprs(i,k+1))/RG
2676	ENDDO
2677	ENDDO
2678	DO i = 1, klon
2679	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2680	/z_apres(i)
2681	ENDDO
2682	DO k = 1, klev
2683	DO i = 1, klon
2684	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2685	z_factor(i).LT.(1.0-1.0E-08)) THEN
2686	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2687	ENDIF
2688	ENDDO
2689	ENDDO
2690	ENDIF
2691	zx_ajustq=.FALSE.
2692
2693	!
2694	!==========================================================================
2695	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
2696	!pour la couche limite diffuse pour l instant
2697	!
2698	!
2699	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
2700	! il faut rajouter cette tendance calcul\'ee hors des poches
2701	! froides
2702	!
2703	IF (iflag_wake>=1) THEN
2704	!
2705	!
2706	! Call wakes every "wkpas" step
2707	!
2708	IF (MOD(itapwk,wkpas).EQ.0) THEN
2709	!
2710	DO k=1,klev
2711	DO i=1,klon
2712	dt_dwn(i,k) = ftd(i,k)
2713	dq_dwn(i,k) = fqd(i,k)
2714	M_dwn(i,k) = dnwd0(i,k)
2715	M_up(i,k) = upwd(i,k)
2716	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2717	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2718	ENDDO
2719	ENDDO
2720
2721	IF (iflag_wake==2) THEN
2722	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2723	DO k = 1,klev
2724	dt_dwn(:,k)= dt_dwn(:,k)+ &
2725	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2726	dq_dwn(:,k)= dq_dwn(:,k)+ &
2727	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2728	ENDDO
2729	ELSEIF (iflag_wake==3) THEN
2730	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2731	DO k = 1,klev
2732	DO i=1,klon
2733	IF (rneb(i,k)==0.) THEN
2734	! On ne tient compte des tendances qu'en dehors des
2735	! nuages (c'est-\`a-dire a priri dans une region ou
2736	! l'eau se reevapore).
2737	dt_dwn(i,k)= dt_dwn(i,k)+ &
2738	ok_wk_lsp(i)*d_t_lsc(i,k)/dtime
2739	dq_dwn(i,k)= dq_dwn(i,k)+ &
2740	ok_wk_lsp(i)*d_q_lsc(i,k)/dtime
2741	ENDIF
2742	ENDDO
2743	ENDDO
2744	ENDIF
2745
2746	!
2747	!calcul caracteristiques de la poche froide
2748	CALL calWAKE (iflag_wake_tend, paprs, pplay, dtime, &
2749	t_seri, q_seri, omega, &
2750	dt_dwn, dq_dwn, M_dwn, M_up, &
2751	dt_a, dq_a, &
2752	sigd, &
2753	wake_deltat, wake_deltaq, wake_s, wake_dens, &
2754	wake_dth, wake_h, &
2755	wake_pe, wake_fip, wake_gfl, &
2756	d_t_wake, d_q_wake, &
2757	wake_k, t_x, q_x, &
2758	wake_omgbdth, wake_dp_omgb, &
2759	wake_dtKE, wake_dqKE, &
2760	wake_omg, wake_dp_deltomg, &
2761	wake_spread, wake_Cstar, d_deltat_wk_gw, &
2762	d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk)
2763	!
2764	!jyg Reinitialize itapwk when wakes have been called
2765	itapwk = 0
2766	ENDIF ! (MOD(itapwk,wkpas).EQ.0)
2767	!
2768	itapwk = itapwk+1
2769	!
2770	!-----------------------------------------------------------------------
2771	! ajout des tendances des poches froides
2772	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
2773	abortphy,flag_inhib_tend,itap,0)
2774	call prt_enerbil('wake',itap)
2775	!------------------------------------------------------------------------
2776
2777	! Increment Wake state variables
2778	IF (iflag_wake_tend .GT. 0.) THEN
2779
2780	CALL add_wake_tend &
2781	(d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk, wake_k, &
2782	'wake', abortphy)
2783	call prt_enerbil('wake',itap)
2784	ENDIF ! (iflag_wake_tend .GT. 0.)
2785
2786	ENDIF ! (iflag_wake>=1)
2787	!
2788	!===================================================================
2789	! Convection seche (thermiques ou ajustement)
2790	!===================================================================
2791	!
2792	CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
2793	,seuil_inversion,weak_inversion,dthmin)
2794
2795
2796
2797	d_t_ajsb(:,:)=0.
2798	d_q_ajsb(:,:)=0.
2799	d_t_ajs(:,:)=0.
2800	d_u_ajs(:,:)=0.
2801	d_v_ajs(:,:)=0.
2802	d_q_ajs(:,:)=0.
2803	clwcon0th(:,:)=0.
2804	!
2805	! fm_therm(:,:)=0.
2806	! entr_therm(:,:)=0.
2807	! detr_therm(:,:)=0.
2808	!
2809	IF (prt_level>9) WRITE(lunout,*) &
2810	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
2811	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2812	IF (iflag_thermals<0) THEN
2813	! Rien
2814	! ====
2815	IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
2816
2817
2818	ELSE
2819
2820	! Thermiques
2821	! ==========
2822	IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
2823	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2824
2825
2826	!cc nrlmd le 10/04/2012
2827	DO k=1,klev+1
2828	DO i=1,klon
2829	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2830	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2831	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2832	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2833	ENDDO
2834	ENDDO
2835	!cc fin nrlmd le 10/04/2012
2836
2837	IF (iflag_thermals>=1) THEN
2838	!jyg<
2839	!! IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2840	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
2841	! Appel des thermiques avec les profils exterieurs aux poches
2842	DO k=1,klev
2843	DO i=1,klon
2844	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2845	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2846	u_therm(i,k) = u_seri(i,k)
2847	v_therm(i,k) = v_seri(i,k)
2848	ENDDO
2849	ENDDO
2850	ELSE
2851	! Appel des thermiques avec les profils moyens
2852	DO k=1,klev
2853	DO i=1,klon
2854	t_therm(i,k) = t_seri(i,k)
2855	q_therm(i,k) = q_seri(i,k)
2856	u_therm(i,k) = u_seri(i,k)
2857	v_therm(i,k) = v_seri(i,k)
2858	ENDDO
2859	ENDDO
2860	ENDIF
2861	!>jyg
2862	CALL calltherm(pdtphys &
2863	,pplay,paprs,pphi,weak_inversion &
2864	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
2865	,u_therm,v_therm,t_therm,q_therm,zqsat,debut & !jyg
2866	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
2867	,fm_therm,entr_therm,detr_therm &
2868	,zqasc,clwcon0th,lmax_th,ratqscth &
2869	,ratqsdiff,zqsatth &
2870	!on rajoute ale et alp, et les
2871	!caracteristiques de la couche alim
2872	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
2873	,ztv,zpspsk,ztla,zthl &
2874	!cc nrlmd le 10/04/2012
2875	,pbl_tke_input,pctsrf,omega,cell_area &
2876	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
2877	,n2,s2,ale_bl_stat &
2878	,therm_tke_max,env_tke_max &
2879	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
2880	,alp_bl_conv,alp_bl_stat &
2881	!cc fin nrlmd le 10/04/2012
2882	,zqla,ztva )
2883	!
2884	!jyg<
2885	!!jyg IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2886	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
2887	! Si les thermiques ne sont presents que hors des
2888	! poches, la tendance moyenne associ\'ee doit etre
2889	! multipliee par la fraction surfacique qu'ils couvrent.
2890	DO k=1,klev
2891	DO i=1,klon
2892	!
2893	d_deltat_the(i,k) = - d_t_ajs(i,k)
2894	d_deltaq_the(i,k) = - d_q_ajs(i,k)
2895	!
2896	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
2897	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
2898	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
2899	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
2900	!
2901	ENDDO
2902	ENDDO
2903	!
2904	CALL add_wake_tend &
2905	(d_deltat_the, d_deltaq_the, dsig0, ddens0, wkoccur1, 'the', abortphy)
2906	call prt_enerbil('the',itap)
2907	!
2908	ENDIF ! (mod(iflag_pbl_split/10,10) .GE. 1)
2909	!
2910	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
2911	dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
2912	call prt_enerbil('thermals',itap)
2913	!
2914	!
2915	CALL alpale_th( dtime, lmax_th, t_seri, cell_area, &
2916	cin, s2, n2, &
2917	ale_bl_trig, ale_bl_stat, ale_bl, &
2918	alp_bl, alp_bl_stat, &
2919	proba_notrig, random_notrig)
2920	!>jyg
2921
2922	! ------------------------------------------------------------------
2923	! Transport de la TKE par les panaches thermiques.
2924	! FH : 2010/02/01
2925	! if (iflag_pbl.eq.10) then
2926	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2927	! s rg,paprs,pbl_tke)
2928	! endif
2929	! -------------------------------------------------------------------
2930
2931	DO i=1,klon
2932	! zmax_th(i)=pphi(i,lmax_th(i))/rg
2933	!CR:04/05/12:correction calcul zmax
2934	zmax_th(i)=zmax0(i)
2935	ENDDO
2936
2937	ENDIF
2938
2939	! Ajustement sec
2940	! ==============
2941
2942	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
2943	! a partir du sommet des thermiques.
2944	! Dans le cas contraire, on demarre au niveau 1.
2945
2946	IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
2947
2948	IF (iflag_thermals.eq.0) THEN
2949	IF (prt_level>9) WRITE(lunout,*)'ajsec'
2950	limbas(:)=1
2951	ELSE
2952	limbas(:)=lmax_th(:)
2953	ENDIF
2954
2955	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2956	! pour des test de convergence numerique.
2957	! Le nouveau ajsec est a priori mieux, meme pour le cas
2958	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2959	! non nulles numeriquement pour des mailles non concernees.
2960
2961	IF (iflag_thermals==0) THEN
2962	! Calling adjustment alone (but not the thermal plume model)
2963	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
2964	, d_t_ajsb, d_q_ajsb)
2965	ELSE IF (iflag_thermals>0) THEN
2966	! Calling adjustment above the top of thermal plumes
2967	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
2968	, d_t_ajsb, d_q_ajsb)
2969	ENDIF
2970
2971	!--------------------------------------------------------------------
2972	! ajout des tendances de l'ajustement sec ou des thermiques
2973	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
2974	'ajsb',abortphy,flag_inhib_tend,itap,0)
2975	call prt_enerbil('ajsb',itap)
2976	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2977	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2978
2979	!---------------------------------------------------------------------
2980
2981	ENDIF
2982
2983	ENDIF
2984	!
2985	!===================================================================
2986	! Computation of ratqs, the width (normalized) of the subrid scale
2987	! water distribution
2988	CALL calcratqs(klon,klev,prt_level,lunout, &
2989	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
2990	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
2991	tau_ratqs,fact_cldcon, &
2992	ptconv,ptconvth,clwcon0th, rnebcon0th, &
2993	paprs,pplay,q_seri,zqsat,fm_therm, &
2994	ratqs,ratqsc)
2995
2996
2997	!
2998	! Appeler le processus de condensation a grande echelle
2999	! et le processus de precipitation
3000	!-------------------------------------------------------------------------
3001	IF (prt_level .GE.10) THEN
3002	print *,'itap, ->fisrtilp ',itap
3003	ENDIF
3004	!
3005	CALL fisrtilp(dtime,paprs,pplay, &
3006	t_seri, q_seri,ptconv,ratqs, &
3007	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3008	rain_lsc, snow_lsc, &
3009	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3010	frac_impa, frac_nucl, beta_prec_fisrt, &
3011	prfl, psfl, rhcl, &
3012	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3013	iflag_ice_thermo)
3014	!
3015	WHERE (rain_lsc < 0) rain_lsc = 0.
3016	WHERE (snow_lsc < 0) snow_lsc = 0.
3017
3018	!+JLD
3019	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3020	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
3021	! & ,rain_lsc,snow_lsc
3022	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3023	!-JLD
3024	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
3025	'lsc',abortphy,flag_inhib_tend,itap,0)
3026	call prt_enerbil('lsc',itap)
3027	rain_num(:)=0.
3028	DO k = 1, klev
3029	DO i = 1, klon
3030	IF (ql_seri(i,k)>oliqmax) THEN
3031	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3032	ql_seri(i,k)=oliqmax
3033	ENDIF
3034	ENDDO
3035	ENDDO
3036	IF (nqo==3) THEN
3037	DO k = 1, klev
3038	DO i = 1, klon
3039	IF (qs_seri(i,k)>oicemax) THEN
3040	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3041	qs_seri(i,k)=oicemax
3042	ENDIF
3043	ENDDO
3044	ENDDO
3045	ENDIF
3046
3047	!---------------------------------------------------------------------------
3048	DO k = 1, klev
3049	DO i = 1, klon
3050	cldfra(i,k) = rneb(i,k)
3051	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3052	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3053	ENDDO
3054	ENDDO
3055	IF (check) THEN
3056	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3057	WRITE(lunout,*)"apresilp=", za
3058	zx_t = 0.0
3059	za = 0.0
3060	DO i = 1, klon
3061	za = za + cell_area(i)/REAL(klon)
3062	zx_t = zx_t + (rain_lsc(i) &
3063	+ snow_lsc(i))*cell_area(i)/REAL(klon)
3064	ENDDO
3065	zx_t = zx_t/za*dtime
3066	WRITE(lunout,*)"Precip=", zx_t
3067	ENDIF
3068
3069	IF (mydebug) THEN
3070	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3071	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3072	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3073	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3074	ENDIF
3075
3076	!
3077	!-------------------------------------------------------------------
3078	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3079	!-------------------------------------------------------------------
3080
3081	! 1. NUAGES CONVECTIFS
3082	!
3083	!IM cf FH
3084	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3085	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3086	snow_tiedtke=0.
3087	! print*,'avant calcul de la pseudo precip '
3088	! print*,'iflag_cld_th',iflag_cld_th
3089	IF (iflag_cld_th.eq.-1) THEN
3090	rain_tiedtke=rain_con
3091	ELSE
3092	! print*,'calcul de la pseudo precip '
3093	rain_tiedtke=0.
3094	! print*,'calcul de la pseudo precip 0'
3095	DO k=1,klev
3096	DO i=1,klon
3097	IF (d_q_con(i,k).lt.0.) THEN
3098	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3099	*(paprs(i,k)-paprs(i,k+1))/rg
3100	ENDIF
3101	ENDDO
3102	ENDDO
3103	ENDIF
3104	!
3105	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3106	!
3107
3108	! Nuages diagnostiques pour Tiedtke
3109	CALL diagcld1(paprs,pplay, &
3110	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3111	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3112	diafra,dialiq)
3113	DO k = 1, klev
3114	DO i = 1, klon
3115	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3116	cldliq(i,k) = dialiq(i,k)
3117	cldfra(i,k) = diafra(i,k)
3118	ENDIF
3119	ENDDO
3120	ENDDO
3121
3122	ELSE IF (iflag_cld_th.ge.3) THEN
3123	! On prend pour les nuages convectifs le max du calcul de la
3124	! convection et du calcul du pas de temps precedent diminue d'un facteur
3125	! facttemps
3126	facteur = pdtphys *facttemps
3127	DO k=1,klev
3128	DO i=1,klon
3129	rnebcon(i,k)=rnebcon(i,k)*facteur
3130	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
3131	rnebcon(i,k)=rnebcon0(i,k)
3132	clwcon(i,k)=clwcon0(i,k)
3133	ENDIF
3134	ENDDO
3135	ENDDO
3136
3137	! On prend la somme des fractions nuageuses et des contenus en eau
3138
3139	IF (iflag_cld_th>=5) THEN
3140
3141	DO k=1,klev
3142	ptconvth(:,k)=fm_therm(:,k+1)>0.
3143	ENDDO
3144
3145	IF (iflag_coupl==4) THEN
3146
3147	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3148	! convectives et lsc dans la partie des thermiques
3149	! Le controle par iflag_coupl est peut etre provisoire.
3150	DO k=1,klev
3151	DO i=1,klon
3152	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
3153	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3154	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3155	ELSE IF (ptconv(i,k)) THEN
3156	cldfra(i,k)=rnebcon(i,k)
3157	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3158	ENDIF
3159	ENDDO
3160	ENDDO
3161
3162	ELSE IF (iflag_coupl==5) THEN
3163	DO k=1,klev
3164	DO i=1,klon
3165	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3166	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3167	ENDDO
3168	ENDDO
3169
3170	ELSE
3171
3172	! Si on est sur un point touche par la convection
3173	! profonde et pas par les thermiques, on prend la
3174	! couverture nuageuse et l'eau nuageuse de la convection
3175	! profonde.
3176
3177	!IM/FH: 2011/02/23
3178	! definition des points sur lesquels ls thermiques sont actifs
3179
3180	DO k=1,klev
3181	DO i=1,klon
3182	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
3183	cldfra(i,k)=rnebcon(i,k)
3184	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3185	ENDIF
3186	ENDDO
3187	ENDDO
3188
3189	ENDIF
3190
3191	ELSE
3192
3193	! Ancienne version
3194	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3195	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3196	ENDIF
3197
3198	ENDIF
3199
3200	! plulsc(:)=0.
3201	! do k=1,klev,-1
3202	! do i=1,klon
3203	! zzz=prfl(:,k)+psfl(:,k)
3204	! if (.not.ptconvth.zzz.gt.0.)
3205	! enddo prfl, psfl,
3206	! enddo
3207	!
3208	! 2. NUAGES STARTIFORMES
3209	!
3210	IF (ok_stratus) THEN
3211	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3212	DO k = 1, klev
3213	DO i = 1, klon
3214	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3215	cldliq(i,k) = dialiq(i,k)
3216	cldfra(i,k) = diafra(i,k)
3217	ENDIF
3218	ENDDO
3219	ENDDO
3220	ENDIF
3221	!
3222	! Precipitation totale
3223	!
3224	DO i = 1, klon
3225	rain_fall(i) = rain_con(i) + rain_lsc(i)
3226	snow_fall(i) = snow_con(i) + snow_lsc(i)
3227	ENDDO
3228	!
3229	! Calculer l'humidite relative pour diagnostique
3230	!
3231	DO k = 1, klev
3232	DO i = 1, klon
3233	zx_t = t_seri(i,k)
3234	IF (thermcep) THEN
3235	!! if (iflag_ice_thermo.eq.0) then !jyg
3236	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3237	!! else !jyg
3238	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
3239	!! endif !jyg
3240	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3241	zx_qs = MIN(0.5,zx_qs)
3242	zcor = 1./(1.-retv*zx_qs)
3243	zx_qs = zx_qs*zcor
3244	ELSE
3245	!! IF (zx_t.LT.t_coup) THEN !jyg
3246	IF (zx_t.LT.rtt) THEN !jyg
3247	zx_qs = qsats(zx_t)/pplay(i,k)
3248	ELSE
3249	zx_qs = qsatl(zx_t)/pplay(i,k)
3250	ENDIF
3251	ENDIF
3252	zx_rh(i,k) = q_seri(i,k)/zx_qs
3253	zqsat(i,k)=zx_qs
3254	ENDDO
3255	ENDDO
3256
3257	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3258	! equivalente a 2m (tpote) pour diagnostique
3259	!
3260	DO i = 1, klon
3261	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3262	IF (thermcep) THEN
3263	IF(zt2m(i).LT.RTT) then
3264	Lheat=RLSTT
3265	ELSE
3266	Lheat=RLVTT
3267	ENDIF
3268	ELSE
3269	IF (zt2m(i).LT.RTT) THEN
3270	Lheat=RLSTT
3271	ELSE
3272	Lheat=RLVTT
3273	ENDIF
3274	ENDIF
3275	tpote(i) = tpot(i)* &
3276	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3277	ENDDO
3278
3279	IF (type_trac == 'inca') THEN
3280	#ifdef INCA
3281	CALL VTe(VTphysiq)
3282	CALL VTb(VTinca)
3283	calday = REAL(days_elapsed + 1) + jH_cur
3284
3285	CALL chemtime(itap+itau_phy-1, date0, dtime, itap)
3286	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
3287	CALL AEROSOL_METEO_CALC( &
3288	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3289	prfl,psfl,pctsrf,cell_area, &
3290	latitude_deg,longitude_deg,u10m,v10m)
3291	ENDIF
3292
3293	zxsnow_dummy(:) = 0.0
3294
3295	CALL chemhook_begin (calday, &
3296	days_elapsed+1, &
3297	jH_cur, &
3298	pctsrf(1,1), &
3299	latitude_deg, &
3300	longitude_deg, &
3301	cell_area, &
3302	paprs, &
3303	pplay, &
3304	coefh(1:klon,1:klev,is_ave), &
3305	pphi, &
3306	t_seri, &
3307	u, &
3308	v, &
3309	wo(:, :, 1), &
3310	q_seri, &
3311	zxtsol, &
3312	zxsnow_dummy, &
3313	solsw, &
3314	albsol1, &
3315	rain_fall, &
3316	snow_fall, &
3317	itop_con, &
3318	ibas_con, &
3319	cldfra, &
3320	nbp_lon, &
3321	nbp_lat-1, &
3322	tr_seri, &
3323	ftsol, &
3324	paprs, &
3325	cdragh, &
3326	cdragm, &
3327	pctsrf, &
3328	pdtphys, &
3329	itap)
3330
3331	CALL VTe(VTinca)
3332	CALL VTb(VTphysiq)
3333	#endif
3334	ENDIF !type_trac = inca
3335
3336
3337	!
3338	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3339	!
3340	IF (MOD(itaprad,radpas).EQ.0) THEN
3341
3342	!
3343	!jq - introduce the aerosol direct and first indirect radiative forcings
3344	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3345	IF (flag_aerosol .GT. 0) THEN
3346	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3347	IF (.NOT. aerosol_couple) THEN
3348	!
3349	CALL readaerosol_optic( &
3350	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3351	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3352	mass_solu_aero, mass_solu_aero_pi, &
3353	tau_aero, piz_aero, cg_aero, &
3354	tausum_aero, tau3d_aero)
3355	ENDIF
3356	ELSE ! RRTM radiation
3357	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3358	abort_message='config_inca=aero et rrtm=1 impossible'
3359	CALL abort_physic(modname,abort_message,1)
3360	ELSE
3361	!
3362	#ifdef CPP_RRTM
3363	IF (NSW.EQ.6) THEN
3364	!--new aerosol properties SW and LW
3365	!
3366	#ifdef CPP_Dust
3367	!--SPL aerosol model
3368	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
3369	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3370	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3371	tausum_aero, tau3d_aero)
3372	#else
3373	!--climatologies or INCA aerosols
3374	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, &
3375	new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
3376	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3377	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3378	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3379	tausum_aero, drytausum_aero, tau3d_aero)
3380	#endif
3381	!
3382	ELSE IF (NSW.EQ.2) THEN
3383	!--for now we use the old aerosol properties
3384	!
3385	CALL readaerosol_optic( &
3386	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3387	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3388	mass_solu_aero, mass_solu_aero_pi, &
3389	tau_aero, piz_aero, cg_aero, &
3390	tausum_aero, tau3d_aero)
3391	!
3392	!--natural aerosols
3393	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
3394	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
3395	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
3396	!--all aerosols
3397	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
3398	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
3399	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
3400	!
3401	!--no LW optics
3402	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3403	!
3404	ELSE
3405	abort_message='Only NSW=2 or 6 are possible with ' &
3406	// 'aerosols and iflag_rrtm=1'
3407	CALL abort_physic(modname,abort_message,1)
3408	ENDIF
3409	#else
3410	abort_message='You should compile with -rrtm if running ' &
3411	// 'with iflag_rrtm=1'
3412	CALL abort_physic(modname,abort_message,1)
3413	#endif
3414	!
3415	ENDIF
3416	ENDIF
3417	ELSE !--flag_aerosol = 0
3418	tausum_aero(:,:,:) = 0.
3419	drytausum_aero(:,:) = 0.
3420	mass_solu_aero(:,:) = 0.
3421	mass_solu_aero_pi(:,:) = 0.
3422	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3423	tau_aero(:,:,:,:) = 1.e-15
3424	piz_aero(:,:,:,:) = 1.
3425	cg_aero(:,:,:,:) = 0.
3426	ELSE
3427	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
3428	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3429	piz_aero_sw_rrtm(:,:,:,:) = 1.0
3430	cg_aero_sw_rrtm(:,:,:,:) = 0.0
3431	ENDIF
3432	ENDIF
3433	!
3434	!--STRAT AEROSOL
3435	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
3436	IF (flag_aerosol_strat.GT.0) THEN
3437	IF (prt_level .GE.10) THEN
3438	PRINT *,'appel a readaerosolstrat', mth_cur
3439	ENDIF
3440	IF (iflag_rrtm.EQ.0) THEN
3441	IF (flag_aerosol_strat.EQ.1) THEN
3442	CALL readaerosolstrato(debut)
3443	ELSE
3444	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
3445	CALL abort_physic(modname,abort_message,1)
3446	ENDIF
3447	ELSE
3448	#ifdef CPP_RRTM
3449	#ifndef CPP_StratAer
3450	!--prescribed strat aerosols
3451	!--only in the case of non-interactive strat aerosols
3452	IF (flag_aerosol_strat.EQ.1) THEN
3453	CALL readaerosolstrato1_rrtm(debut)
3454	ELSEIF (flag_aerosol_strat.EQ.2) THEN
3455	CALL stratosphere_mask(t_seri, pplay, latitude_deg)
3456	CALL readaerosolstrato2_rrtm(debut)
3457	ELSE
3458	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
3459	CALL abort_physic(modname,abort_message,1)
3460	ENDIF
3461	#endif
3462	#else
3463	abort_message='You should compile with -rrtm if running ' &
3464	// 'with iflag_rrtm=1'
3465	CALL abort_physic(modname,abort_message,1)
3466	#endif
3467	ENDIF
3468	ENDIF
3469	!
3470	#ifdef CPP_RRTM
3471	#ifdef CPP_StratAer
3472	!--compute stratospheric mask
3473	CALL stratosphere_mask(t_seri, pplay, latitude_deg)
3474	!--interactive strat aerosols
3475	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
3476	#endif
3477	#endif
3478	!--fin STRAT AEROSOL
3479	!
3480
3481	! Calculer les parametres optiques des nuages et quelques
3482	! parametres pour diagnostiques:
3483	!
3484	IF (aerosol_couple.AND.config_inca=='aero') THEN
3485	mass_solu_aero(:,:) = ccm(:,:,1)
3486	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3487	ENDIF
3488
3489	IF (ok_newmicro) then
3490	IF (iflag_rrtm.NE.0) THEN
3491	#ifdef CPP_RRTM
3492	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3493	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
3494	// 'pour ok_cdnc'
3495	CALL abort_physic(modname,abort_message,1)
3496	ENDIF
3497	#else
3498
3499	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
3500	CALL abort_physic(modname,abort_message,1)
3501	#endif
3502	ENDIF
3503	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3504	paprs, pplay, t_seri, cldliq, cldfra, &
3505	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3506	flwp, fiwp, flwc, fiwc, &
3507	mass_solu_aero, mass_solu_aero_pi, &
3508	cldtaupi, re, fl, ref_liq, ref_ice, &
3509	ref_liq_pi, ref_ice_pi)
3510	ELSE
3511	CALL nuage (paprs, pplay, &
3512	t_seri, cldliq, cldfra, cldtau, cldemi, &
3513	cldh, cldl, cldm, cldt, cldq, &
3514	ok_aie, &
3515	mass_solu_aero, mass_solu_aero_pi, &
3516	bl95_b0, bl95_b1, &
3517	cldtaupi, re, fl)
3518	ENDIF
3519	!
3520	!IM betaCRF
3521	!
3522	cldtaurad = cldtau
3523	cldtaupirad = cldtaupi
3524	cldemirad = cldemi
3525	cldfrarad = cldfra
3526
3527	!
3528	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3529	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3530	!
3531	! global
3532	!
3533	DO k=1, klev
3534	DO i=1, klon
3535	IF (pplay(i,k).GE.pfree) THEN
3536	beta(i,k) = beta_pbl
3537	ELSE
3538	beta(i,k) = beta_free
3539	ENDIF
3540	IF (mskocean_beta) THEN
3541	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3542	ENDIF
3543	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3544	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3545	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3546	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3547	ENDDO
3548	ENDDO
3549	!
3550	ELSE
3551	!
3552	! regional
3553	!
3554	DO k=1, klev
3555	DO i=1,klon
3556	!
3557	IF (longitude_deg(i).ge.lon1_beta.AND. &
3558	longitude_deg(i).le.lon2_beta.AND. &
3559	latitude_deg(i).le.lat1_beta.AND. &
3560	latitude_deg(i).ge.lat2_beta) THEN
3561	IF (pplay(i,k).GE.pfree) THEN
3562	beta(i,k) = beta_pbl
3563	ELSE
3564	beta(i,k) = beta_free
3565	ENDIF
3566	IF (mskocean_beta) THEN
3567	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3568	ENDIF
3569	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3570	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3571	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3572	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3573	ENDIF
3574	!
3575	ENDDO
3576	ENDDO
3577	!
3578	ENDIF
3579
3580	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
3581	IF (ok_chlorophyll) THEN
3582	print*,"-- reading chlorophyll"
3583	CALL readchlorophyll(debut)
3584	ENDIF
3585
3586	!--if ok_suntime_rrtm we use ancillay data for RSUN
3587	!--previous values are therefore overwritten
3588	!--this is needed for CMIP6 runs
3589	!--and only possible for new radiation scheme
3590	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
3591	#ifdef CPP_RRTM
3592	CALL read_rsun_rrtm(debut)
3593	#endif
3594	ENDIF
3595
3596	IF (mydebug) THEN
3597	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3598	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3599	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3600	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3601	ENDIF
3602
3603	!
3604	!sonia : If Iflag_radia >=2, pertubation of some variables
3605	!input to radiation (DICE)
3606	!
3607	IF (iflag_radia .ge. 2) THEN
3608	zsav_tsol (:) = zxtsol(:)
3609	CALL perturb_radlwsw(zxtsol,iflag_radia)
3610	ENDIF
3611
3612	IF (aerosol_couple.AND.config_inca=='aero') THEN
3613	#ifdef INCA
3614	CALL radlwsw_inca &
3615	(kdlon,kflev,dist, rmu0, fract, solaire, &
3616	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3617	size(wo,3), wo, &
3618	cldfrarad, cldemirad, cldtaurad, &
3619	heat,heat0,cool,cool0,albpla, &
3620	topsw,toplw,solsw,sollw, &
3621	sollwdown, &
3622	topsw0,toplw0,solsw0,sollw0, &
3623	lwdn0, lwdn, lwup0, lwup, &
3624	swdn0, swdn, swup0, swup, &
3625	ok_ade, ok_aie, &
3626	tau_aero, piz_aero, cg_aero, &
3627	topswad_aero, solswad_aero, &
3628	topswad0_aero, solswad0_aero, &
3629	topsw_aero, topsw0_aero, &
3630	solsw_aero, solsw0_aero, &
3631	cldtaupirad, &
3632	topswai_aero, solswai_aero)
3633	#endif
3634	ELSE
3635	!
3636	!IM calcul radiatif pour le cas actuel
3637	!
3638	RCO2 = RCO2_act
3639	RCH4 = RCH4_act
3640	RN2O = RN2O_act
3641	RCFC11 = RCFC11_act
3642	RCFC12 = RCFC12_act
3643	!
3644	IF (prt_level .GE.10) THEN
3645	print *,' ->radlwsw, number 1 '
3646	ENDIF
3647	!
3648	CALL radlwsw &
3649	(dist, rmu0, fract, &
3650	!albedo SB >>>
3651	! paprs, pplay,zxtsol,albsol1, albsol2, &
3652	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3653	!albedo SB <<<
3654	t_seri,q_seri,wo, &
3655	cldfrarad, cldemirad, cldtaurad, &
3656	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
3657	flag_aerosol_strat, &
3658	tau_aero, piz_aero, cg_aero, &
3659	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3660	! Rajoute par OB pour RRTM
3661	tau_aero_lw_rrtm, &
3662	cldtaupirad,new_aod, &
3663	zqsat, flwc, fiwc, &
3664	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3665	heat,heat0,cool,cool0,albpla, &
3666	topsw,toplw,solsw,sollw, &
3667	sollwdown, &
3668	topsw0,toplw0,solsw0,sollw0, &
3669	lwdn0, lwdn, lwup0, lwup, &
3670	swdn0, swdn, swup0, swup, &
3671	topswad_aero, solswad_aero, &
3672	topswai_aero, solswai_aero, &
3673	topswad0_aero, solswad0_aero, &
3674	topsw_aero, topsw0_aero, &
3675	solsw_aero, solsw0_aero, &
3676	topswcf_aero, solswcf_aero, &
3677	!-C. Kleinschmitt for LW diagnostics
3678	toplwad_aero, sollwad_aero,&
3679	toplwai_aero, sollwai_aero, &
3680	toplwad0_aero, sollwad0_aero,&
3681	!-end
3682	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3683	ZSWFT0_i, ZFSDN0, ZFSUP0)
3684
3685	#ifndef CPP_XIOS
3686	!--OB 30/05/2016 modified 21/10/2016
3687	!--here we return swaero_diag and dryaod_diag to FALSE
3688	!--and histdef will switch it back to TRUE if necessary
3689	!--this is necessary to get the right swaero at first step
3690	!--but only in the case of no XIOS as XIOS is covered elsewhere
3691	IF (debut) swaero_diag = .FALSE.
3692	IF (debut) dryaod_diag = .FALSE.
3693	#endif
3694	!
3695	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3696	!IM des taux doit etre different du taux actuel
3697	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3698	!
3699	IF (ok_4xCO2atm) THEN
3700	IF (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3701	RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3702	RCFC12_per.NE.RCFC12_act) THEN
3703	!
3704	RCO2 = RCO2_per
3705	RCH4 = RCH4_per
3706	RN2O = RN2O_per
3707	RCFC11 = RCFC11_per
3708	RCFC12 = RCFC12_per
3709	!
3710	IF (prt_level .GE.10) THEN
3711	print *,' ->radlwsw, number 2 '
3712	ENDIF
3713	!
3714	CALL radlwsw &
3715	(dist, rmu0, fract, &
3716	!albedo SB >>>
3717	! paprs, pplay,zxtsol,albsol1, albsol2, &
3718	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3719	!albedo SB <<<
3720	t_seri,q_seri,wo, &
3721	cldfrarad, cldemirad, cldtaurad, &
3722	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
3723	flag_aerosol_strat, &
3724	tau_aero, piz_aero, cg_aero, &
3725	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3726	! Rajoute par OB pour RRTM
3727	tau_aero_lw_rrtm, &
3728	cldtaupi,new_aod, &
3729	zqsat, flwc, fiwc, &
3730	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3731	heatp,heat0p,coolp,cool0p,albplap, &
3732	topswp,toplwp,solswp,sollwp, &
3733	sollwdownp, &
3734	topsw0p,toplw0p,solsw0p,sollw0p, &
3735	lwdn0p, lwdnp, lwup0p, lwupp, &
3736	swdn0p, swdnp, swup0p, swupp, &
3737	topswad_aerop, solswad_aerop, &
3738	topswai_aerop, solswai_aerop, &
3739	topswad0_aerop, solswad0_aerop, &
3740	topsw_aerop, topsw0_aerop, &
3741	solsw_aerop, solsw0_aerop, &
3742	topswcf_aerop, solswcf_aerop, &
3743	!-C. Kleinschmitt for LW diagnostics
3744	toplwad_aerop, sollwad_aerop,&
3745	toplwai_aerop, sollwai_aerop, &
3746	toplwad0_aerop, sollwad0_aerop,&
3747	!-end
3748	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3749	ZSWFT0_i, ZFSDN0, ZFSUP0)
3750	endif
3751	endif
3752	!
3753	ENDIF ! aerosol_couple
3754	itaprad = 0
3755	!
3756	! If Iflag_radia >=2, reset pertubed variables
3757	!
3758	IF (iflag_radia .ge. 2) THEN
3759	zxtsol(:) = zsav_tsol (:)
3760	ENDIF
3761	ENDIF ! MOD(itaprad,radpas)
3762	itaprad = itaprad + 1
3763
3764	IF (iflag_radia.eq.0) THEN
3765	IF (prt_level.ge.9) THEN
3766	PRINT *,'--------------------------------------------------'
3767	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3768	PRINT *,'>>>> heat et cool mis a zero '
3769	PRINT *,'--------------------------------------------------'
3770	ENDIF
3771	heat=0.
3772	cool=0.
3773	sollw=0. ! MPL 01032011
3774	solsw=0.
3775	radsol=0.
3776	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3777	swup0=0.
3778	lwup=0.
3779	lwup0=0.
3780	lwdn=0.
3781	lwdn0=0.
3782	ENDIF
3783
3784	!
3785	! Calculer radsol a l'exterieur de radlwsw
3786	! pour prendre en compte le cycle diurne
3787	! recode par Olivier Boucher en sept 2015
3788	!
3789	radsol=solsw*swradcorr+sollw
3790
3791	IF (ok_4xCO2atm) THEN
3792	radsolp=solswp*swradcorr+sollwp
3793	ENDIF
3794
3795	!
3796	! Ajouter la tendance des rayonnements (tous les pas)
3797	! avec une correction pour le cycle diurne dans le SW
3798	!
3799
3800	DO k=1, klev
3801	d_t_swr(:,k)=swradcorr(:)heat(:,k)dtime/RDAY
3802	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)dtime/RDAY
3803	d_t_lwr(:,k)=-cool(:,k)*dtime/RDAY
3804	d_t_lw0(:,k)=-cool0(:,k)*dtime/RDAY
3805	ENDDO
3806
3807	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
3808	call prt_enerbil('SW',itap)
3809	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
3810	call prt_enerbil('LW',itap)
3811
3812	!
3813	IF (mydebug) THEN
3814	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3815	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3816	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3817	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3818	ENDIF
3819
3820	! Calculer l'hydrologie de la surface
3821	!
3822	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3823	! . agesno, ftsol,fqsurf,fsnow, ruis)
3824	!
3825
3826	!
3827	! Calculer le bilan du sol et la derive de temperature (couplage)
3828	!
3829	DO i = 1, klon
3830	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3831	! a la demande de JLD
3832	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3833	ENDDO
3834	!
3835	!moddeblott(jan95)
3836	! Appeler le programme de parametrisation de l'orographie
3837	! a l'echelle sous-maille:
3838	!
3839	IF (prt_level .GE.10) THEN
3840	print *,' call orography ? ', ok_orodr
3841	ENDIF
3842	!
3843	IF (ok_orodr) THEN
3844	!
3845	! selection des points pour lesquels le shema est actif:
3846	igwd=0
3847	DO i=1,klon
3848	itest(i)=0
3849	! IF ((zstd(i).gt.10.0)) THEN
3850	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3851	itest(i)=1
3852	igwd=igwd+1
3853	idx(igwd)=i
3854	ENDIF
3855	ENDDO
3856	! igwdim=MAX(1,igwd)
3857	!
3858	IF (ok_strato) THEN
3859
3860	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
3861	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3862	igwd,idx,itest, &
3863	t_seri, u_seri, v_seri, &
3864	zulow, zvlow, zustrdr, zvstrdr, &
3865	d_t_oro, d_u_oro, d_v_oro)
3866
3867	ELSE
3868	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
3869	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3870	igwd,idx,itest, &
3871	t_seri, u_seri, v_seri, &
3872	zulow, zvlow, zustrdr, zvstrdr, &
3873	d_t_oro, d_u_oro, d_v_oro)
3874	ENDIF
3875	!
3876	! ajout des tendances
3877	!-----------------------------------------------------------------------
3878	! ajout des tendances de la trainee de l'orographie
3879	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
3880	abortphy,flag_inhib_tend,itap,0)
3881	call prt_enerbil('oro',itap)
3882	!----------------------------------------------------------------------
3883	!
3884	ENDIF ! fin de test sur ok_orodr
3885	!
3886	IF (mydebug) THEN
3887	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3888	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3889	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3890	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3891	ENDIF
3892
3893	IF (ok_orolf) THEN
3894	!
3895	! selection des points pour lesquels le shema est actif:
3896	igwd=0
3897	DO i=1,klon
3898	itest(i)=0
3899	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3900	itest(i)=1
3901	igwd=igwd+1
3902	idx(igwd)=i
3903	ENDIF
3904	ENDDO
3905	! igwdim=MAX(1,igwd)
3906	!
3907	IF (ok_strato) THEN
3908
3909	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
3910	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3911	igwd,idx,itest, &
3912	t_seri, u_seri, v_seri, &
3913	zulow, zvlow, zustrli, zvstrli, &
3914	d_t_lif, d_u_lif, d_v_lif )
3915
3916	ELSE
3917	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
3918	latitude_deg,zmea,zstd,zpic, &
3919	itest, &
3920	t_seri, u_seri, v_seri, &
3921	zulow, zvlow, zustrli, zvstrli, &
3922	d_t_lif, d_u_lif, d_v_lif)
3923	ENDIF
3924
3925	! ajout des tendances de la portance de l'orographie
3926	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
3927	'lif', abortphy,flag_inhib_tend,itap,0)
3928	call prt_enerbil('lif',itap)
3929	ENDIF ! fin de test sur ok_orolf
3930
3931	IF (ok_hines) then
3932	! HINES GWD PARAMETRIZATION
3933	east_gwstress=0.
3934	west_gwstress=0.
3935	du_gwd_hines=0.
3936	dv_gwd_hines=0.
3937	CALL hines_gwd(klon, klev, dtime, paprs, pplay, latitude_deg, t_seri, &
3938	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
3939	du_gwd_hines, dv_gwd_hines)
3940	zustr_gwd_hines=0.
3941	zvstr_gwd_hines=0.
3942	DO k = 1, klev
3943	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/dtime &
3944	* (paprs(:, k)-paprs(:, k+1))/rg
3945	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/dtime &
3946	* (paprs(:, k)-paprs(:, k+1))/rg
3947	ENDDO
3948
3949	d_t_hin(:, :)=0.
3950	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
3951	dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
3952	call prt_enerbil('hin',itap)
3953	ENDIF
3954
3955	IF (.not. ok_hines .and. ok_gwd_rando) then
3956	CALL acama_GWD_rando(DTIME, pplay, latitude_deg, t_seri, u_seri, &
3957	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
3958	dv_gwd_front, east_gwstress, west_gwstress)
3959	zustr_gwd_front=0.
3960	zvstr_gwd_front=0.
3961	DO k = 1, klev
3962	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/dtime &
3963	* (paprs(:, k)-paprs(:, k+1))/rg
3964	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/dtime &
3965	* (paprs(:, k)-paprs(:, k+1))/rg
3966	ENDDO
3967
3968	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
3969	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
3970	call prt_enerbil('front_gwd_rando',itap)
3971	ENDIF
3972
3973	IF (ok_gwd_rando) THEN
3974	CALL FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
3975	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
3976	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
3977	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
3978	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
3979	call prt_enerbil('flott_gwd_rando',itap)
3980	zustr_gwd_rando=0.
3981	zvstr_gwd_rando=0.
3982	DO k = 1, klev
3983	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/dtime &
3984	* (paprs(:, k)-paprs(:, k+1))/rg
3985	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/dtime &
3986	* (paprs(:, k)-paprs(:, k+1))/rg
3987	ENDDO
3988	ENDIF
3989
3990	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3991
3992	IF (mydebug) THEN
3993	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3994	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3995	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3996	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3997	ENDIF
3998
3999	DO i = 1, klon
4000	zustrph(i)=0.
4001	zvstrph(i)=0.
4002	ENDDO
4003	DO k = 1, klev
4004	DO i = 1, klon
4005	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
4006	(paprs(i,k)-paprs(i,k+1))/rg
4007	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
4008	(paprs(i,k)-paprs(i,k+1))/rg
4009	ENDDO
4010	ENDDO
4011	!
4012	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4013	!
4014	IF (is_sequential .and. ok_orodr) THEN
4015	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4016	ra,rg,romega, &
4017	latitude_deg,longitude_deg,pphis, &
4018	zustrdr,zustrli,zustrph, &
4019	zvstrdr,zvstrli,zvstrph, &
4020	paprs,u,v, &
4021	aam, torsfc)
4022	ENDIF
4023	!IM cf. FLott END
4024	!DC Calcul de la tendance due au methane
4025	IF(ok_qch4) THEN
4026	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4027	! ajout de la tendance d'humidite due au methane
4028	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*dtime
4029	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
4030	'q_ch4', abortphy,flag_inhib_tend,itap,0)
4031	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/dtime
4032	ENDIF
4033	!
4034	!
4035	!====================================================================
4036	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4037	!====================================================================
4038	! Abderrahmane 24.08.09
4039
4040	IF (ok_cosp) THEN
4041	! adeclarer
4042	#ifdef CPP_COSP
4043	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
4044
4045	IF (prt_level .GE.10) THEN
4046	print*,'freq_cosp',freq_cosp
4047	ENDIF
4048	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4049	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4050	! s ref_liq,ref_ice
4051	CALL phys_cosp(itap,dtime,freq_cosp, &
4052	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4053	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4054	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4055	JrNt,ref_liq,ref_ice, &
4056	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4057	zu10m,zv10m,pphis, &
4058	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4059	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4060	prfl(:,1:klev),psfl(:,1:klev), &
4061	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4062	mr_ozone,cldtau, cldemi)
4063
4064	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4065	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4066	! M clMISR,
4067	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4068	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4069
4070	ENDIF
4071
4072	#endif
4073	ENDIF !ok_cosp
4074
4075
4076	! Marine
4077
4078	IF (ok_airs) then
4079
4080	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/dtime)).EQ.0) THEN
4081	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4082	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4083	& map_prop_hc,map_prop_hist,&
4084	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4085	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4086	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4087	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4088	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4089	& map_ntot,map_hc,map_hist,&
4090	& map_Cb,map_ThCi,map_Anv,&
4091	& alt_tropo )
4092	ENDIF
4093
4094	ENDIF ! ok_airs
4095
4096
4097	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4098	!AA
4099	!AA Installation de l'interface online-offline pour traceurs
4100	!AA
4101	!====================================================================
4102	! Calcul des tendances traceurs
4103	!====================================================================
4104	!
4105
4106	IF (type_trac=='repr') THEN
4107	sh_in(:,:) = q_seri(:,:)
4108	ELSE
4109	sh_in(:,:) = qx(:,:,ivap)
4110	ch_in(:,:) = qx(:,:,iliq)
4111	ENDIF
4112
4113	#ifdef CPP_Dust
4114	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
4115	pdtphys,ftsol, & ! I
4116	t,q_seri,paprs,pplay,RHcl, & ! I
4117	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
4118	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
4119	u_seri, v_seri, latitude_deg, longitude_deg, &
4120	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
4121	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
4122	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
4123	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
4124	fm_therm, entr_therm, rneb, & ! I
4125	beta_prec_fisrt,beta_prec, & !I
4126	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
4127	d_tr_dyn,tr_seri)
4128
4129	#else
4130
4131	CALL phytrac ( &
4132	itap, days_elapsed+1, jH_cur, debut, &
4133	lafin, dtime, u, v, t, &
4134	paprs, pplay, pmfu, pmfd, &
4135	pen_u, pde_u, pen_d, pde_d, &
4136	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
4137	u1, v1, ftsol, pctsrf, &
4138	zustar, zu10m, zv10m, &
4139	wstar(:,is_ave), ale_bl, ale_wake, &
4140	latitude_deg, longitude_deg, &
4141	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4142	presnivs, pphis, pphi, albsol1, &
4143	sh_in, ch_in, rhcl, cldfra, rneb, &
4144	diafra, cldliq, itop_con, ibas_con, &
4145	pmflxr, pmflxs, prfl, psfl, &
4146	da, phi, mp, upwd, &
4147	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
4148	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
4149	ev, ep, epmlmMm, eplaMm, & !<<RomP
4150	dnwd, aerosol_couple, flxmass_w, &
4151	tau_aero, piz_aero, cg_aero, ccm, &
4152	rfname, &
4153	d_tr_dyn, & !<<RomP
4154	tr_seri)
4155	#endif
4156
4157	IF (offline) THEN
4158
4159	IF (prt_level.ge.9) &
4160	print*,'Attention on met a 0 les thermiques pour phystoke'
4161	CALL phystokenc ( &
4162	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4163	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4164	fm_therm,entr_therm, &
4165	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4166	frac_impa, frac_nucl, &
4167	pphis,cell_area,dtime,itap, &
4168	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4169
4170
4171	ENDIF
4172
4173	!
4174	! Calculer le transport de l'eau et de l'energie (diagnostique)
4175	!
4176	CALL transp (paprs,zxtsol, &
4177	t_seri, q_seri, u_seri, v_seri, zphi, &
4178	ve, vq, ue, uq)
4179	!
4180	!IM global posePB BEG
4181	IF(1.EQ.0) THEN
4182	!
4183	CALL transp_lay (paprs,zxtsol, &
4184	t_seri, q_seri, u_seri, v_seri, zphi, &
4185	ve_lay, vq_lay, ue_lay, uq_lay)
4186	!
4187	ENDIF !(1.EQ.0) THEN
4188	!IM global posePB END
4189	! Accumuler les variables a stocker dans les fichiers histoire:
4190	!
4191
4192	!================================================================
4193	! Conversion of kinetic and potential energy into heat, for
4194	! parameterisation of subgrid-scale motions
4195	!================================================================
4196
4197	d_t_ec(:,:)=0.
4198	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4199	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
4200	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4201	zmasse,exner,d_t_ec)
4202	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4203
4204	!=======================================================================
4205	! SORTIES
4206	!=======================================================================
4207	!
4208	!IM initialisation + calculs divers diag AMIP2
4209	!
4210	include "calcul_divers.h"
4211	!
4212	!IM Interpolation sur les niveaux de pression du NMC
4213	! -------------------------------------------------
4214	#ifdef CPP_XIOS
4215	!$OMP MASTER
4216	!On recupere la valeur de la missing value donnee dans le xml
4217	CALL xios_get_field_attr("t850",default_value=missing_val_omp)
4218	! PRINT *,"ARNAUD value missing ",missing_val_omp
4219	!$OMP END MASTER
4220	!$OMP BARRIER
4221	missing_val=missing_val_omp
4222	#endif
4223	#ifndef CPP_XIOS
4224	missing_val=missing_val_nf90
4225	#endif
4226	!
4227	include "calcul_STDlev.h"
4228	!
4229	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
4230	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
4231	!
4232	!cc prw = eau precipitable
4233	! prlw = colonne eau liquide
4234	! prlw = colonne eau solide
4235	prw(:) = 0.
4236	prlw(:) = 0.
4237	prsw(:) = 0.
4238	DO k = 1, klev
4239	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
4240	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
4241	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
4242	ENDDO
4243	!
4244	IF (type_trac == 'inca') THEN
4245	#ifdef INCA
4246	CALL VTe(VTphysiq)
4247	CALL VTb(VTinca)
4248
4249	CALL chemhook_end ( &
4250	dtime, &
4251	pplay, &
4252	t_seri, &
4253	tr_seri, &
4254	nbtr, &
4255	paprs, &
4256	q_seri, &
4257	cell_area, &
4258	pphi, &
4259	pphis, &
4260	zx_rh, &
4261	aps, bps)
4262
4263	CALL VTe(VTinca)
4264	CALL VTb(VTphysiq)
4265	#endif
4266	ENDIF
4267
4268
4269	!
4270	! Convertir les incrementations en tendances
4271	!
4272	IF (prt_level .GE.10) THEN
4273	print *,'Convertir les incrementations en tendances '
4274	ENDIF
4275	!
4276	IF (mydebug) THEN
4277	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4278	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4279	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4280	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4281	ENDIF
4282
4283	DO k = 1, klev
4284	DO i = 1, klon
4285	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
4286	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
4287	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
4288	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
4289	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
4290	!CR: on ajoute le contenu en glace
4291	IF (nqo.eq.3) THEN
4292	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
4293	ENDIF
4294	ENDDO
4295	ENDDO
4296	!
4297	!CR: nb de traceurs eau: nqo
4298	! IF (nqtot.GE.3) THEN
4299	IF (nqtot.GE.(nqo+1)) THEN
4300	! DO iq = 3, nqtot
4301	DO iq = nqo+1, nqtot
4302	DO k = 1, klev
4303	DO i = 1, klon
4304	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
4305	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
4306	ENDDO
4307	ENDDO
4308	ENDDO
4309	ENDIF
4310	!
4311	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4312	!IM global posePB include "write_bilKP_ins.h"
4313	!IM global posePB include "write_bilKP_ave.h"
4314	!
4315
4316	!--OB mass fixer
4317	!--profile is corrected to force mass conservation of water
4318	IF (mass_fixer) THEN
4319	qql2(:)=0.0
4320	DO k = 1, klev
4321	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
4322	ENDDO
4323	DO i = 1, klon
4324	!--compute ratio of what q+ql should be with conservation to what it is
4325	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
4326	DO k = 1, klev
4327	q_seri(i,k) =q_seri(i,k)*corrqql
4328	ql_seri(i,k)=ql_seri(i,k)*corrqql
4329	ENDDO
4330	ENDDO
4331	ENDIF
4332	!--fin mass fixer
4333
4334	! Sauvegarder les valeurs de t et q a la fin de la physique:
4335	!
4336	u_ancien(:,:) = u_seri(:,:)
4337	v_ancien(:,:) = v_seri(:,:)
4338	t_ancien(:,:) = t_seri(:,:)
4339	q_ancien(:,:) = q_seri(:,:)
4340	ql_ancien(:,:) = ql_seri(:,:)
4341	qs_ancien(:,:) = qs_seri(:,:)
4342	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
4343	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
4344	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
4345	! !! RomP >>>
4346	!CR: nb de traceurs eau: nqo
4347	IF (nqtot.GT.nqo) THEN
4348	DO iq = nqo+1, nqtot
4349	tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
4350	ENDDO
4351	ENDIF
4352	! !! RomP <<<
4353	!==========================================================================
4354	! Sorties des tendances pour un point particulier
4355	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4356	! pour le debug
4357	! La valeur de igout est attribuee plus haut dans le programme
4358	!==========================================================================
4359
4360	IF (prt_level.ge.1) THEN
4361	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4362	write(lunout,*) &
4363	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4364	write(lunout,*) &
4365	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4366	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4367	pctsrf(igout,is_sic)
4368	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4369	DO k=1,klev
4370	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4371	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4372	d_t_eva(igout,k)
4373	ENDDO
4374	write(lunout,*) 'cool,heat'
4375	DO k=1,klev
4376	write(lunout,*) cool(igout,k),heat(igout,k)
4377	ENDDO
4378
4379	!jyg< (En attendant de statuer sur le sort de d_t_oli)
4380	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4381	!jyg! do k=1,klev
4382	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4383	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4384	!jyg! enddo
4385	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4386	DO k=1,klev
4387	write(lunout,*) d_t_vdf(igout,k), &
4388	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4389	ENDDO
4390	!>jyg
4391
4392	write(lunout,*) 'd_ps ',d_ps(igout)
4393	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4394	DO k=1,klev
4395	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4396	d_qx(igout,k,1),d_qx(igout,k,2)
4397	ENDDO
4398	ENDIF
4399
4400	!============================================================
4401	! Calcul de la temperature potentielle
4402	!============================================================
4403	DO k = 1, klev
4404	DO i = 1, klon
4405	!JYG/IM theta en debut du pas de temps
4406	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4407	!JYG/IM theta en fin de pas de temps de physique
4408	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4409	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
4410	! MPL 20130625
4411	! fth_fonctions.F90 et parkind1.F90
4412	! sinon thetal=theta
4413	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4414	! : ql_seri(i,k))
4415	thetal(i,k)=theta(i,k)
4416	ENDDO
4417	ENDDO
4418	!
4419
4420	! 22.03.04 BEG
4421	!=============================================================
4422	! Ecriture des sorties
4423	!=============================================================
4424	#ifdef CPP_IOIPSL
4425
4426	! Recupere des varibles calcule dans differents modules
4427	! pour ecriture dans histxxx.nc
4428
4429	! Get some variables from module fonte_neige_mod
4430	CALL fonte_neige_get_vars(pctsrf, &
4431	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
4432
4433
4434	!=============================================================
4435	! Separation entre thermiques et non thermiques dans les sorties
4436	! de fisrtilp
4437	!=============================================================
4438
4439	IF (iflag_thermals>=1) THEN
4440	d_t_lscth=0.
4441	d_t_lscst=0.
4442	d_q_lscth=0.
4443	d_q_lscst=0.
4444	DO k=1,klev
4445	DO i=1,klon
4446	IF (ptconvth(i,k)) THEN
4447	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4448	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4449	ELSE
4450	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4451	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4452	ENDIF
4453	ENDDO
4454	ENDDO
4455
4456	DO i=1,klon
4457	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4458	plul_th(i)=prfl(i,1)+psfl(i,1)
4459	ENDDO
4460	ENDIF
4461
4462	!On effectue les sorties:
4463
4464	#ifdef CPP_Dust
4465	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
4466	pplay, lmax_th, aerosol_couple, &
4467	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
4468	ptconv, read_climoz, clevSTD, &
4469	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
4470	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4471	#else
4472	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
4473	pplay, lmax_th, aerosol_couple, &
4474	ok_ade, ok_aie, ivap, iliq, isol, new_aod, &
4475	ok_sync, ptconv, read_climoz, clevSTD, &
4476	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
4477	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4478	#endif
4479
4480	#ifndef CPP_XIOS
4481	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
4482	#endif
4483
4484	#endif
4485
4486
4487	!====================================================================
4488	! Arret du modele apres hgardfou en cas de detection d'un
4489	! plantage par hgardfou
4490	!====================================================================
4491
4492	IF (abortphy==1) THEN
4493	abort_message ='Plantage hgardfou'
4494	CALL abort_physic (modname,abort_message,1)
4495	ENDIF
4496
4497	! 22.03.04 END
4498	!
4499	!====================================================================
4500	! Si c'est la fin, il faut conserver l'etat de redemarrage
4501	!====================================================================
4502	!
4503
4504	IF (lafin) THEN
4505	itau_phy = itau_phy + itap
4506	CALL phyredem ("restartphy.nc")
4507	! open(97,form="unformatted",file="finbin")
4508	! write(97) u_seri,v_seri,t_seri,q_seri
4509	! close(97)
4510	!$OMP MASTER
4511	IF (read_climoz >= 1) THEN
4512	IF (is_mpi_root) THEN
4513	CALL nf95_close(ncid_climoz)
4514	ENDIF
4515	DEALLOCATE(press_edg_climoz) ! pointer
4516	DEALLOCATE(press_cen_climoz) ! pointer
4517	ENDIF
4518	!$OMP END MASTER
4519	ENDIF
4520
4521	! first=.false.
4522
4523
4524	END SUBROUTINE physiq
4525
4526	END MODULE physiq_mod

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

Download in other formats:

Original Format