Context Navigation

physiq_mod.F90 @ 4648

Last change on this file since 4648 was 4647, checked in by idelkadi, 14 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: 196.1 KB

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

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

Download in other formats:

Original Format