Context Navigation

physiq_mod.F90 @ 4202

Last change on this file since 4202 was 4188, checked in by idelkadi, 2 years ago
Implementation of the Ecrad radiative transfer code in the LMD model (continued) : Integration of aerosols (direct effect)
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: 191.5 KB

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

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

Download in other formats:

Original Format