Context Navigation

physiq_mod.F90 @ 4516

Last change on this file since 4516 was 4516, checked in by jyg, 13 months ago

ok_bug_ajs_cv = Logical switch to a bug : modifying directly wake_deltat
by adding the (w) dry adjustment tendency to wake_deltat.

IF ok_bug_ajs_cv is false, then tendencies are not added to wake_deltat

(and are taken into account through only the heating and moistening
profiles).

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: 193.6 KB

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

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

Download in other formats:

Original Format