Context Navigation

physiq_mod.F90 @ 5423

Last change on this file since 5423 was 3899, checked in by lguez, 4 years ago
Old work on cold pool death.
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.7 KB

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

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

Download in other formats:

Original Format