Context Navigation

physiq_mod.F90 @ 4120

Last change on this file since 4120 was 4120, checked in by dcugnet, 2 years ago
New water names: H2Ov, H2Ol, H2Oi, H2Or -> H2O_g, H2O_l, H2O_s, H2O_r. Corrections for the lOldCode=.FALSE., not activated yet.
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: 190.0 KB

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

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

Download in other formats:

Original Format