Context Navigation

physiq_mod.F90 @ 3820

Last change on this file since 3820 was 3820, checked in by lmdz-users, 3 years ago
to be able to work with dynamico's dynamic, we modify call for the inca initialisation
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 186.8 KB

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

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

Download in other formats:

Original Format