Context Navigation

physiq_mod.F90 @ 3861

Last change on this file since 3861 was 3861, checked in by acozic, 3 years ago
merge with rev 3860
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: 185.6 KB

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

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

Download in other formats:

Original Format