Context Navigation

physiq_mod.F90 @ 3082

Last change on this file since 3082 was 3082, checked in by oboucher, 6 years ago
Adding new diagnostics for clean (no aerosol) clear sky SW radiative fluxes
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: 164.8 KB

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

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

Download in other formats:

Original Format