Context Navigation

physiq_mod.F90 @ 4677

Last change on this file since 4677 was 4677, checked in by idelkadi, 8 months ago

Implementation in the LMDZ code of the double call of the ECRAD radiative transfer code to estimate the 3D radiative effect of clouds.

This double call of Ecrad is controlled by the ok_3Deffect logic key.
If this key is enabled, 2 files of parameter configuration "namelists" for ECRAD are required at runtime: namelist_ecrad and namelist_ecrad_s2.
If this key is deactivated, the configuration and initialization part (reading namelist and netcdf files) is performed only once during simulation (1st call to ECRAD). Otherwise, configuration and initialization are performed each time Ecrad is called.

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

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

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

Download in other formats:

Original Format