Context Navigation

physiq_mod.F90 @ 4412

Last change on this file since 4412 was 4412, checked in by evignon, 16 months ago
travail de l'atelier nuages du 30/01/23: on renomme cldliq (ou radliq) en radocond car le nom est tres trompeur + on ajoute des commentaires dans lscp_mod
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 `Id`
File size: 192.2 KB

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

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

Download in other formats:

Original Format