Context Navigation

physiq_mod.F90 @ 3956

Last change on this file since 3956 was 3956, checked in by jyg, 3 years ago
Bug fixes concerning various variables ill-initialized, ill-used, ill-printed, or ill-placed. + cv_gen moved from phys_local_var_mod.F90 to phys_state_var_mod.F90; ==> changes in physiq_mod.F90 and phys_output_write.F90 + awake_dens added in phys_state_var_mod.F90 + cv_gen and awake_dens now initialized in phyetat0.F90 and written in phyredem.F90 + cv_gen, awake_dens, and solswfdiff now initialized in old_lmdz1d.F90 and scm.F90 + useless variables suppressed in pbl_surface_mod.F90.
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: 186.5 KB

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

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

Download in other formats:

Original Format