Context Navigation

physiq_mod.F90 @ 3784

Last change on this file since 3784 was 3784, checked in by evignon, 4 years ago
correction v3783 pour convergence et compilation ancienne physique, Etienne aide par Ehouarn
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.7 KB

Line
1	!
2	! $Id: physiq_mod.F90 3784 2020-11-09 14:09:22Z evignon $
3	!
4	!#define IO_DEBUG
5	MODULE physiq_mod
6
7	IMPLICIT NONE
8
9	CONTAINS
10
11	SUBROUTINE physiq (nlon,nlev, &
12	debut,lafin,pdtphys_, &
13	paprs,pplay,pphi,pphis,presnivs, &
14	u,v,rot,t,qx, &
15	flxmass_w, &
16	d_u, d_v, d_t, d_qx, d_ps)
17
18	! For clarity, the "USE" section is now arranged in alphabetical order,
19	! with a separate section for CPP keys
20	! PLEASE try to follow this rule
21
22	USE ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
23	USE aero_mod
24	USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
25	& fl_ebil, fl_cor_ebil
26	USE assert_m, only: assert
27	USE change_srf_frac_mod
28	USE conf_phys_m, only: conf_phys
29	USE carbon_cycle_mod, ONLY : infocfields_init, RCO2_glo, carbon_cycle_rad
30	USE CFMIP_point_locations ! IM stations CFMIP
31	USE cmp_seri_mod
32	USE dimphy
33	USE etat0_limit_unstruct_mod
34	USE FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
35	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
36	USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
37	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
38	histwrite, ju2ymds, ymds2ju, getin
39	USE ioipsl_getin_p_mod, ONLY : getin_p
40	USE indice_sol_mod
41	USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac
42	USE iophy
43	USE limit_read_mod, ONLY : init_limit_read
44	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, grid1dTo2d_glo, grid_type, unstructured
45	USE mod_phys_lmdz_mpi_data, only: is_mpi_root
46	USE mod_phys_lmdz_para
47	USE netcdf95, only: nf95_close
48	USE netcdf, only: nf90_fill_real ! IM for NMC files
49	USE open_climoz_m, only: open_climoz ! ozone climatology from a file
50	USE ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
51	USE pbl_surface_mod, ONLY : pbl_surface
52	USE phyaqua_mod, only: zenang_an
53	USE phystokenc_mod, ONLY: offline, phystokenc
54	USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
55	year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour
56	!! USE phys_local_var_mod, ONLY : a long list of variables
57	!! ==> see below, after "CPP Keys" section
58	USE phys_state_var_mod ! Variables sauvegardees de la physique
59	USE phys_output_mod
60	USE phys_output_ctrlout_mod
61	USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
62	USE readaerosol_mod, ONLY : init_aero_fromfile
63	USE readaerosolstrato_m, ONLY : init_readaerosolstrato
64	USE radlwsw_m, only: radlwsw
65	USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
66	USE regr_pr_time_av_m, only: regr_pr_time_av
67	USE surface_data, ONLY : type_ocean, ok_veget, ok_snow
68	USE time_phylmdz_mod, only: annee_ref, current_time, day_ini, day_ref, &
69	day_step_phy, itau_phy, pdtphys, raz_date, start_time, update_time
70	USE tracinca_mod, ONLY: config_inca
71	USE tropopause_m, ONLY: dyn_tropopause
72	USE vampir
73	USE VERTICAL_LAYERS_MOD, ONLY: aps,bps, ap, bp
74	USE write_field_phy
75
76	!USE cmp_seri_mod
77	! USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
78	! & fl_ebil, fl_cor_ebil
79
80	!!!!!!!!!!!!!!!!!! "USE" section for CPP keys !!!!!!!!!!!!!!!!!!!!!!!!
81	!
82	!
83	#ifdef CPP_Dust
84	USE phytracr_spl_mod, ONLY: phytracr_spl, phytracr_spl_out_init
85	USE phys_output_write_spl_mod
86	#else
87	USE phytrac_mod, ONLY : phytrac_init, phytrac
88	USE phys_output_write_mod
89	#endif
90
91
92	#ifdef REPROBUS
93	USE CHEM_REP, ONLY : Init_chem_rep_xjour, &
94	d_q_rep,d_ql_rep,d_qi_rep,ptrop,ttrop, &
95	ztrop, gravit,itroprep, Z1,Z2,fac,B
96	#endif
97
98
99	#ifdef CPP_RRTM
100	USE YOERAD, ONLY : NRADLP
101	USE YOESW, ONLY : RSUN
102	#endif
103
104
105	#ifdef CPP_StratAer
106	USE strataer_mod, ONLY: strataer_init
107	#endif
108
109
110	#ifdef CPP_XIOS
111	USE xios, ONLY: xios_update_calendar, xios_context_finalize, &
112	xios_get_field_attr, xios_field_is_active
113	USE wxios, ONLY: missing_val, missing_val_omp
114	#endif
115	#ifndef CPP_XIOS
116	USE paramLMDZ_phy_mod
117	#endif
118	!
119	!
120	!!!!!!!!!!!!!!!!!! END "USE" for CPP keys !!!!!!!!!!!!!!!!!!!!!!
121
122	USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
123	! [Variables internes non sauvegardees de la physique]
124	! Variables locales pour effectuer les appels en serie
125	t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri, &
126	! Dynamic tendencies (diagnostics)
127	d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn, &
128	d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
129	! Physic tendencies
130	d_t_con,d_q_con,d_u_con,d_v_con, &
131	d_tr, & !! to be removed?? (jyg)
132	d_t_wake,d_q_wake, &
133	d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
134	d_t_ajsb,d_q_ajsb, &
135	d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
136	d_t_ajs_w,d_q_ajs_w, &
137	d_t_ajs_x,d_q_ajs_x, &
138	!
139	d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
140	d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
141	d_t_lscst,d_q_lscst, &
142	d_t_lscth,d_q_lscth, &
143	plul_st,plul_th, &
144	!
145	d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
146	d_ts, &
147	!
148	d_t_oli,d_u_oli,d_v_oli, &
149	d_t_oro,d_u_oro,d_v_oro, &
150	d_t_oro_gw,d_u_oro_gw,d_v_oro_gw, &
151	d_t_lif,d_u_lif,d_v_lif, &
152	d_t_ec, &
153	!
154	du_gwd_hines,dv_gwd_hines,d_t_hin, &
155	dv_gwd_rando,dv_gwd_front, &
156	east_gwstress,west_gwstress, &
157	d_q_ch4, &
158	! Special RRTM
159	ZLWFT0_i,ZSWFT0_i,ZFLDN0, &
160	ZFLUP0,ZFSDN0,ZFSUP0, &
161	!
162	topswad_aero,solswad_aero, &
163	topswai_aero,solswai_aero, &
164	topswad0_aero,solswad0_aero, &
165	!LW additional
166	toplwad_aero,sollwad_aero, &
167	toplwai_aero,sollwai_aero, &
168	toplwad0_aero,sollwad0_aero, &
169	!
170	topsw_aero,solsw_aero, &
171	topsw0_aero,solsw0_aero, &
172	topswcf_aero,solswcf_aero, &
173	tausum_aero,tau3d_aero, &
174	drytausum_aero, &
175	!
176	!variables CFMIP2/CMIP5
177	topswad_aerop, solswad_aerop, &
178	topswai_aerop, solswai_aerop, &
179	topswad0_aerop, solswad0_aerop, &
180	topsw_aerop, topsw0_aerop, &
181	solsw_aerop, solsw0_aerop, &
182	topswcf_aerop, solswcf_aerop, &
183	!LW diagnostics
184	toplwad_aerop, sollwad_aerop, &
185	toplwai_aerop, sollwai_aerop, &
186	toplwad0_aerop, sollwad0_aerop, &
187	!
188	ptstar, pt0, slp, &
189	!
190	bils, &
191	!
192	cldh, cldl,cldm, cldq, cldt, &
193	JrNt, &
194	dthmin, evap, fder, plcl, plfc, &
195	prw, prlw, prsw, &
196	s_lcl, s_pblh, s_pblt, s_therm, &
197	cdragm, cdragh, &
198	zustar, zu10m, zv10m, rh2m, qsat2m, &
199	zq2m, zt2m, 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 chemini( &
1948	rg, &
1949	ra, &
1950	cell_area, &
1951	latitude_deg, &
1952	longitude_deg, &
1953	presnivs, &
1954	calday, &
1955	klon, &
1956	nqtot, &
1957	nqo, &
1958	pdtphys, &
1959	annee_ref, &
1960	year_cur, &
1961	day_ref, &
1962	day_ini, &
1963	start_time, &
1964	itau_phy, &
1965	date0, &
1966	io_lon, &
1967	io_lat, &
1968	chemistry_couple, &
1969	init_source, &
1970	init_tauinca, &
1971	init_pizinca, &
1972	init_cginca, &
1973	init_ccminca)
1974
1975
1976	! initialisation des variables depuis le restart de inca
1977	ccm(:,:,:) = init_ccminca
1978	tau_aero(:,:,:,:) = init_tauinca
1979	piz_aero(:,:,:,:) = init_pizinca
1980	cg_aero(:,:,:,:) = init_cginca
1981	!
1982
1983
1984	CALL VTe(VTinca)
1985	CALL VTb(VTphysiq)
1986	#endif
1987	ENDIF
1988	IF (type_trac == 'repr') THEN
1989	#ifdef REPROBUS
1990	CALL chemini_rep( &
1991	presnivs, &
1992	pdtphys, &
1993	annee_ref, &
1994	day_ref, &
1995	day_ini, &
1996	start_time, &
1997	itau_phy, &
1998	io_lon, &
1999	io_lat)
2000	#endif
2001	ENDIF
2002
2003	!$omp single
2004	IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, &
2005	press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
2006	!$omp end single
2007	!
2008	!IM betaCRF
2009	pfree=70000. !Pa
2010	beta_pbl=1.
2011	beta_free=1.
2012	lon1_beta=-180.
2013	lon2_beta=+180.
2014	lat1_beta=90.
2015	lat2_beta=-90.
2016	mskocean_beta=.FALSE.
2017
2018	!albedo SB >>>
2019	SELECT CASE(nsw)
2020	CASE(2)
2021	SFRWL(1)=0.45538747
2022	SFRWL(2)=0.54461211
2023	CASE(4)
2024	SFRWL(1)=0.45538747
2025	SFRWL(2)=0.32870591
2026	SFRWL(3)=0.18568763
2027	SFRWL(4)=3.02191470E-02
2028	CASE(6)
2029	SFRWL(1)=1.28432794E-03
2030	SFRWL(2)=0.12304168
2031	SFRWL(3)=0.33106142
2032	SFRWL(4)=0.32870591
2033	SFRWL(5)=0.18568763
2034	SFRWL(6)=3.02191470E-02
2035	END SELECT
2036
2037
2038	!albedo SB <<<
2039
2040	OPEN(99,file='beta_crf.data',status='old', &
2041	form='formatted',err=9999)
2042	READ(99,*,end=9998) pfree
2043	READ(99,*,end=9998) beta_pbl
2044	READ(99,*,end=9998) beta_free
2045	READ(99,*,end=9998) lon1_beta
2046	READ(99,*,end=9998) lon2_beta
2047	READ(99,*,end=9998) lat1_beta
2048	READ(99,*,end=9998) lat2_beta
2049	READ(99,*,end=9998) mskocean_beta
2050	9998 Continue
2051	CLOSE(99)
2052	9999 Continue
2053	WRITE(,)'pfree=',pfree
2054	WRITE(,)'beta_pbl=',beta_pbl
2055	WRITE(,)'beta_free=',beta_free
2056	WRITE(,)'lon1_beta=',lon1_beta
2057	WRITE(,)'lon2_beta=',lon2_beta
2058	WRITE(,)'lat1_beta=',lat1_beta
2059	WRITE(,)'lat2_beta=',lat2_beta
2060	WRITE(,)'mskocean_beta=',mskocean_beta
2061
2062	!lwoff=y : offset LW CRE for radiation code and other schemes
2063	!lwoff=y : betalwoff=1.
2064	betalwoff=0.
2065	IF (ok_lwoff) THEN
2066	betalwoff=1.
2067	ENDIF
2068	WRITE(,)'ok_lwoff=',ok_lwoff
2069	!
2070	!lwoff=y to begin only sollw and sollwdown are set up to CS values
2071	sollw = sollw + betalwoff * (sollw0 - sollw)
2072	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
2073	sollwdown(:))
2074
2075
2076	ENDIF
2077	!
2078	! ************** Fin de IF ( debut ) *************
2079	!
2080	!
2081	! Incrementer le compteur de la physique
2082	!
2083	itap = itap + 1
2084	IF (is_master .OR. prt_level > 9) THEN
2085	IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
2086	WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
2087	WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
2088	100 FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
2089	ENDIF
2090	ENDIF
2091	!
2092	!
2093	! Update fraction of the sub-surfaces (pctsrf) and
2094	! initialize, where a new fraction has appeared, all variables depending
2095	! on the surface fraction.
2096	!
2097	CALL change_srf_frac(itap, phys_tstep, days_elapsed+1, &
2098	pctsrf, fevap, z0m, z0h, agesno, &
2099	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
2100
2101	! Update time and other variables in Reprobus
2102	IF (type_trac == 'repr') THEN
2103	#ifdef REPROBUS
2104	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
2105	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
2106	CALL Rtime(debut)
2107	#endif
2108	ENDIF
2109
2110	! Tendances bidons pour les processus qui n'affectent pas certaines
2111	! variables.
2112	du0(:,:)=0.
2113	dv0(:,:)=0.
2114	dt0 = 0.
2115	dq0(:,:)=0.
2116	dql0(:,:)=0.
2117	dqi0(:,:)=0.
2118	dsig0(:) = 0.
2119	ddens0(:) = 0.
2120	wkoccur1(:)=1
2121	!
2122	! Mettre a zero des variables de sortie (pour securite)
2123	!
2124	DO i = 1, klon
2125	d_ps(i) = 0.0
2126	ENDDO
2127	DO k = 1, klev
2128	DO i = 1, klon
2129	d_t(i,k) = 0.0
2130	d_u(i,k) = 0.0
2131	d_v(i,k) = 0.0
2132	ENDDO
2133	ENDDO
2134	DO iq = 1, nqtot
2135	DO k = 1, klev
2136	DO i = 1, klon
2137	d_qx(i,k,iq) = 0.0
2138	ENDDO
2139	ENDDO
2140	ENDDO
2141	beta_prec_fisrt(:,:)=0.
2142	beta_prec(:,:)=0.
2143	!
2144	! Output variables from the convective scheme should not be set to 0
2145	! since convection is not always called at every time step.
2146	IF (ok_bug_cv_trac) THEN
2147	da(:,:)=0.
2148	mp(:,:)=0.
2149	phi(:,:,:)=0.
2150	! RomP >>>
2151	phi2(:,:,:)=0.
2152	epmlmMm(:,:,:)=0.
2153	eplaMm(:,:)=0.
2154	d1a(:,:)=0.
2155	dam(:,:)=0.
2156	pmflxr(:,:)=0.
2157	pmflxs(:,:)=0.
2158	! RomP <<<
2159	ENDIF
2160
2161	!
2162	! Ne pas affecter les valeurs entrees de u, v, h, et q
2163	!
2164	DO k = 1, klev
2165	DO i = 1, klon
2166	t_seri(i,k) = t(i,k)
2167	u_seri(i,k) = u(i,k)
2168	v_seri(i,k) = v(i,k)
2169	q_seri(i,k) = qx(i,k,ivap)
2170	ql_seri(i,k) = qx(i,k,iliq)
2171	!CR: ATTENTION, on rajoute la variable glace
2172	IF (nqo.eq.2) THEN
2173	qs_seri(i,k) = 0.
2174	ELSE IF (nqo.eq.3) THEN
2175	qs_seri(i,k) = qx(i,k,isol)
2176	ENDIF
2177	ENDDO
2178	ENDDO
2179	!
2180	!--OB mass fixer
2181	IF (mass_fixer) THEN
2182	!--store initial water burden
2183	qql1(:)=0.0
2184	DO k = 1, klev
2185	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
2186	ENDDO
2187	ENDIF
2188	!--fin mass fixer
2189
2190	tke0(:,:)=pbl_tke(:,:,is_ave)
2191	!CR:Nombre de traceurs de l'eau: nqo
2192	! IF (nqtot.GE.3) THEN
2193	IF (nqtot.GE.(nqo+1)) THEN
2194	! DO iq = 3, nqtot
2195	DO iq = nqo+1, nqtot
2196	DO k = 1, klev
2197	DO i = 1, klon
2198	! tr_seri(i,k,iq-2) = qx(i,k,iq)
2199	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
2200	ENDDO
2201	ENDDO
2202	ENDDO
2203	ELSE
2204	DO k = 1, klev
2205	DO i = 1, klon
2206	tr_seri(i,k,1) = 0.0
2207	ENDDO
2208	ENDDO
2209	ENDIF
2210	!
2211	! Temporary solutions adressing ticket #104 and the non initialisation of tr_ancien
2212	! LF
2213	IF (debut) THEN
2214	WRITE(lunout,*)' WARNING: tr_ancien initialised to tr_seri'
2215	DO iq = nqo+1, nqtot
2216	tr_ancien(:,:,iq-nqo)=tr_seri(:,:,iq-nqo)
2217	ENDDO
2218	ENDIF
2219	!
2220	DO i = 1, klon
2221	ztsol(i) = 0.
2222	ENDDO
2223	DO nsrf = 1, nbsrf
2224	DO i = 1, klon
2225	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2226	ENDDO
2227	ENDDO
2228	! Initialize variables used for diagnostic purpose
2229	IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
2230
2231	! Diagnostiquer la tendance dynamique
2232	!
2233	IF (ancien_ok) THEN
2234	!
2235	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
2236	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
2237	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
2238	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
2239	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
2240	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
2241	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
2242	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
2243	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
2244	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
2245	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
2246	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
2247	! !! RomP >>> td dyn traceur
2248	IF (nqtot.GT.nqo) THEN ! jyg
2249	DO iq = nqo+1, nqtot ! jyg
2250	d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/phys_tstep ! jyg
2251	ENDDO
2252	ENDIF
2253	! !! RomP <<<
2254	ELSE
2255	d_u_dyn(:,:) = 0.0
2256	d_v_dyn(:,:) = 0.0
2257	d_t_dyn(:,:) = 0.0
2258	d_q_dyn(:,:) = 0.0
2259	d_ql_dyn(:,:) = 0.0
2260	d_qs_dyn(:,:) = 0.0
2261	d_q_dyn2d(:) = 0.0
2262	d_ql_dyn2d(:) = 0.0
2263	d_qs_dyn2d(:) = 0.0
2264	! !! RomP >>> td dyn traceur
2265	IF (nqtot.GT.nqo) THEN ! jyg
2266	DO iq = nqo+1, nqtot ! jyg
2267	d_tr_dyn(:,:,iq-nqo)= 0.0 ! jyg
2268	ENDDO
2269	ENDIF
2270	! !! RomP <<<
2271	ancien_ok = .TRUE.
2272	ENDIF
2273	!
2274	! Ajouter le geopotentiel du sol:
2275	!
2276	DO k = 1, klev
2277	DO i = 1, klon
2278	zphi(i,k) = pphi(i,k) + pphis(i)
2279	ENDDO
2280	ENDDO
2281	!
2282	! Verifier les temperatures
2283	!
2284	!IM BEG
2285	IF (check) THEN
2286	amn=MIN(ftsol(1,is_ter),1000.)
2287	amx=MAX(ftsol(1,is_ter),-1000.)
2288	DO i=2, klon
2289	amn=MIN(ftsol(i,is_ter),amn)
2290	amx=MAX(ftsol(i,is_ter),amx)
2291	ENDDO
2292	!
2293	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2294	ENDIF !(check) THEN
2295	!IM END
2296	!
2297	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
2298	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
2299
2300	!
2301	!IM BEG
2302	IF (check) THEN
2303	amn=MIN(ftsol(1,is_ter),1000.)
2304	amx=MAX(ftsol(1,is_ter),-1000.)
2305	DO i=2, klon
2306	amn=MIN(ftsol(i,is_ter),amn)
2307	amx=MAX(ftsol(i,is_ter),amx)
2308	ENDDO
2309	!
2310	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2311	ENDIF !(check) THEN
2312	!IM END
2313	!
2314	! Mettre en action les conditions aux limites (albedo, sst, etc.).
2315	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
2316	!
2317	! Update ozone if day change
2318	IF (MOD(itap-1,lmt_pas) == 0) THEN
2319	IF (read_climoz <= 0) THEN
2320	! Once per day, update ozone from Royer:
2321	IF (solarlong0<-999.) then
2322	! Generic case with evolvoing season
2323	zzz=real(days_elapsed+1)
2324	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
2325	! Particular case with annual mean insolation
2326	zzz=real(90) ! could be revisited
2327	IF (read_climoz/=-1) THEN
2328	abort_message ='read_climoz=-1 is recommended when ' &
2329	// 'solarlong0=1000.'
2330	CALL abort_physic (modname,abort_message,1)
2331	ENDIF
2332	ELSE
2333	! Case where the season is imposed with solarlong0
2334	zzz=real(90) ! could be revisited
2335	ENDIF
2336
2337	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
2338	#ifdef REPROBUS
2339	ptrop=dyn_tropopause(t_seri, ztsol, paprs, pplay, rot)/100.
2340	DO i = 1, klon
2341	Z1=t_seri(i,itroprep(i)+1)
2342	Z2=t_seri(i,itroprep(i))
2343	fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2344	B=Z2-fac*alog(pplay(i,itroprep(i)))
2345	ttrop(i)= fac*alog(ptrop(i))+B
2346	!
2347	Z1= 1.e-3 * ( pphi(i,itroprep(i)+1)+pphis(i) ) / gravit
2348	Z2= 1.e-3 * ( pphi(i,itroprep(i)) +pphis(i) ) / gravit
2349	fac=(Z1-Z2)/alog(pplay(i,itroprep(i)+1)/pplay(i,itroprep(i)))
2350	B=Z2-fac*alog(pplay(i,itroprep(i)))
2351	ztrop(i)=fac*alog(ptrop(i))+B
2352	ENDDO
2353	#endif
2354	ELSE
2355	!--- ro3i = elapsed days number since current year 1st january, 0h
2356	ro3i=days_elapsed+jh_cur-jh_1jan
2357	!--- scaling for old style files (360 records)
2358	IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
2359	IF(adjust_tropopause) THEN
2360	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2361	ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), &
2362	time_climoz , longitude_deg, latitude_deg, &
2363	dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
2364	ELSE
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 )
2368	ENDIF
2369	! Convert from mole fraction of ozone to column density of ozone in a
2370	! cell, in kDU:
2371	FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2372	* zmasse / dobson_u / 1e3
2373	! (By regridding ozone values for LMDZ only once a day, we
2374	! have already neglected the variation of pressure in one
2375	! day. So do not recompute "wo" at each time step even if
2376	! "zmasse" changes a little.)
2377	ENDIF
2378	ENDIF
2379	!
2380	! Re-evaporer l'eau liquide nuageuse
2381	!
2382	CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2383	& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
2384
2385	CALL add_phys_tend &
2386	(du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
2387	'eva',abortphy,flag_inhib_tend,itap,0)
2388	CALL prt_enerbil('eva',itap)
2389
2390	!=========================================================================
2391	! Calculs de l'orbite.
2392	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2393	! doit donc etre plac\'e avant radlwsw et pbl_surface
2394
2395	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2396	CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2397	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2398	!
2399	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2400	! solarlong0
2401	IF (solarlong0<-999.) THEN
2402	IF (new_orbit) THEN
2403	! calcul selon la routine utilisee pour les planetes
2404	CALL solarlong(day_since_equinox, zlongi, dist)
2405	ELSE
2406	! calcul selon la routine utilisee pour l'AR4
2407	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2408	ENDIF
2409	ELSE
2410	zlongi=solarlong0 ! longitude solaire vraie
2411	dist=1. ! distance au soleil / moyenne
2412	ENDIF
2413
2414	IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2415
2416
2417	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2418	! Calcul de l'ensoleillement :
2419	! ============================
2420	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2421	! l'annee a partir d'une formule analytique.
2422	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2423	! non nul aux poles.
2424	IF (abs(solarlong0-1000.)<1.e-4) THEN
2425	CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
2426	latitude_deg,longitude_deg,rmu0,fract)
2427	swradcorr(:) = 1.0
2428	JrNt(:) = 1.0
2429	zrmu0(:) = rmu0(:)
2430	ELSE
2431	! recode par Olivier Boucher en sept 2015
2432	SELECT CASE (iflag_cycle_diurne)
2433	CASE(0)
2434	! Sans cycle diurne
2435	CALL angle(zlongi, latitude_deg, fract, rmu0)
2436	swradcorr = 1.0
2437	JrNt = 1.0
2438	zrmu0 = rmu0
2439	CASE(1)
2440	! Avec cycle diurne sans application des poids
2441	! bit comparable a l ancienne formulation cycle_diurne=true
2442	! on integre entre gmtime et gmtime+radpas
2443	zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
2444	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2445	latitude_deg,longitude_deg,rmu0,fract)
2446	zrmu0 = rmu0
2447	swradcorr = 1.0
2448	! Calcul du flag jour-nuit
2449	JrNt = 0.0
2450	WHERE (fract.GT.0.0) JrNt = 1.0
2451	CASE(2)
2452	! Avec cycle diurne sans application des poids
2453	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2454	! Comme cette routine est appele a tous les pas de temps de
2455	! la physique meme si le rayonnement n'est pas appele je
2456	! remonte en arriere les radpas-1 pas de temps
2457	! suivant. Petite ruse avec MOD pour prendre en compte le
2458	! premier pas de temps de la physique pendant lequel
2459	! itaprad=0
2460	zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)
2461	zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
2462	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2463	latitude_deg,longitude_deg,rmu0,fract)
2464	!
2465	! Calcul des poids
2466	!
2467	zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
2468	zdtime2=0.0 !--pas de temps de la physique qui se termine
2469	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2470	latitude_deg,longitude_deg,zrmu0,zfract)
2471	swradcorr = 0.0
2472	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2473	swradcorr=zfract/fract*zrmu0/rmu0
2474	! Calcul du flag jour-nuit
2475	JrNt = 0.0
2476	WHERE (zfract.GT.0.0) JrNt = 1.0
2477	END SELECT
2478	ENDIF
2479	sza_o = ACOS (rmu0) *180./pi
2480
2481	IF (mydebug) THEN
2482	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2483	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2484	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2485	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2486	ENDIF
2487
2488	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2489	! Appel au pbl_surface : Planetary Boudary Layer et Surface
2490	! Cela implique tous les interactions des sous-surfaces et la
2491	! partie diffusion turbulent du couche limit.
2492	!
2493	! Certains varibales de sorties de pbl_surface sont utiliser que pour
2494	! ecriture des fihiers hist_XXXX.nc, ces sont :
2495	! qsol, zq2m, s_pblh, s_lcl,
2496	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2497	! s_therm, s_trmb1, s_trmb2, s_trmb3,
2498	! zu10m, zv10m, fder,
2499	! zxqsurf, rh2m, zxfluxu, zxfluxv,
2500	! frugs, agesno, fsollw, fsolsw,
2501	! d_ts, fevap, fluxlat, t2m,
2502	! wfbils, wfbilo, fluxt, fluxu, fluxv,
2503	!
2504	! Certains ne sont pas utiliser du tout :
2505	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2506	!
2507
2508	! Calcul de l'humidite de saturation au niveau du sol
2509
2510
2511
2512	IF (iflag_pbl/=0) THEN
2513
2514	!jyg+nrlmd<
2515	!!jyg IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2516	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
2517	print *,'debut du splitting de la PBL, wake_s = ', wake_s(:)
2518	print *,'debut du splitting de la PBL, wake_deltat = ', wake_deltat(:,1)
2519	print *,'debut du splitting de la PBL, wake_deltaq = ', wake_deltaq(:,1)
2520	ENDIF
2521	! !!
2522	!>jyg+nrlmd
2523	!
2524	!-------gustiness calculation-------!
2525	!ym : Warning gustiness non inialized for iflag_gusts=2 & iflag_gusts=3
2526	gustiness=0 !ym missing init
2527
2528	IF (iflag_gusts==0) THEN
2529	gustiness(1:klon)=0
2530	ELSE IF (iflag_gusts==1) THEN
2531	gustiness(1:klon)=f_gust_blale_bl(1:klon)+f_gust_wkale_wake(1:klon)
2532	ELSE IF (iflag_gusts==2) THEN
2533	gustiness(1:klon)=f_gust_blale_bl_stat(1:klon)+f_gust_wkale_wake(1:klon)
2534	! ELSE IF (iflag_gusts==2) THEN
2535	! do i = 1, klon
2536	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
2537	! *ale_wake(i) !! need to make sigma_wk accessible here
2538	! enddo
2539	! ELSE IF (iflag_gusts==3) THEN
2540	! do i = 1, klon
2541	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
2542	! enddo
2543	ENDIF
2544
2545	CALL pbl_surface( &
2546	phys_tstep, date0, itap, days_elapsed+1, &
2547	debut, lafin, &
2548	longitude_deg, latitude_deg, rugoro, zrmu0, &
2549	zsig, sollwdown, pphi, cldt, &
2550	rain_fall, snow_fall, solsw, solswfdiff, sollw, &
2551	gustiness, &
2552	t_seri, q_seri, u_seri, v_seri, &
2553	!nrlmd+jyg<
2554	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2555	!>nrlmd+jyg
2556	pplay, paprs, pctsrf, &
2557	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2558	!albedo SB <<<
2559	cdragh, cdragm, u1, v1, &
2560	!albedo SB >>>
2561	! albsol1, albsol2, sens, evap, &
2562	albsol_dir, albsol_dif, sens, evap, &
2563	!albedo SB <<<
2564	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
2565	zxtsol, zxfluxlat, zt2m, qsat2m, &
2566	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2567	!nrlmd<
2568	!jyg<
2569	d_t_vdf_w, d_q_vdf_w, &
2570	d_t_vdf_x, d_q_vdf_x, &
2571	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2572	!>jyg
2573	delta_tsurf,wake_dens, &
2574	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2575	kh,kh_x,kh_w, &
2576	!>nrlmd
2577	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2578	slab_wfbils, &
2579	qsol, zq2m, s_pblh, s_lcl, &
2580	!jyg<
2581	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2582	!>jyg
2583	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2584	s_therm, s_trmb1, s_trmb2, s_trmb3, &
2585	zustar, zu10m, zv10m, fder, &
2586	zxqsurf, rh2m, zxfluxu, zxfluxv, &
2587	z0m, z0h, agesno, fsollw, fsolsw, &
2588	d_ts, fevap, fluxlat, t2m, &
2589	wfbils, wfbilo, wfevap, wfrain, wfsnow, &
2590	fluxt, fluxu, fluxv, &
2591	dsens, devap, zxsnow, &
2592	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
2593	!nrlmd+jyg<
2594	wake_delta_pbl_TKE, &
2595	!>nrlmd+jyg
2596	treedrg )
2597	!FC
2598	!
2599	! Add turbulent diffusion tendency to the wake difference variables
2600	!!jyg IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2601	IF (mod(iflag_pbl_split,10) .NE. 0) THEN
2602	!jyg<
2603	d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2604	d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2605	CALL add_wake_tend &
2606	(d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
2607	ELSE
2608	d_deltat_vdf(:,:) = 0.
2609	d_deltaq_vdf(:,:) = 0.
2610	!>jyg
2611	ENDIF
2612
2613	!add limitation for t,q at and wind at 10m
2614	if ( iflag_bug_t2m_ipslcm61 == 0 ) THEN
2615	CALL borne_var_surf( klon,klev,nbsrf, &
2616	iflag_bug_t2m_stab_ipslcm61, &
2617	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1), &
2618	ftsol,zxqsurf,pctsrf,paprs, &
2619	t2m, q2m, u10m, v10m, &
2620	zt2m_cor, zq2m_cor, zu10m_cor, zv10m_cor, &
2621	zrh2m_cor, zqsat2m_cor)
2622	ELSE
2623	zt2m_cor(:)=zt2m(:)
2624	zq2m_cor(:)=zq2m(:)
2625	zu10m_cor(:)=zu10m(:)
2626	zv10m_cor(:)=zv10m(:)
2627	zqsat2m_cor=999.999
2628	ENDIF
2629
2630	!---------------------------------------------------------------------
2631	! ajout des tendances de la diffusion turbulente
2632	IF (klon_glo==1) THEN
2633	CALL add_pbl_tend &
2634	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2635	'vdf',abortphy,flag_inhib_tend,itap)
2636	ELSE
2637	CALL add_phys_tend &
2638	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2639	'vdf',abortphy,flag_inhib_tend,itap,0)
2640	ENDIF
2641	CALL prt_enerbil('vdf',itap)
2642	!--------------------------------------------------------------------
2643
2644	IF (mydebug) THEN
2645	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2646	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2647	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2648	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2649	ENDIF
2650
2651	!albedo SB >>>
2652	albsol1=0.
2653	albsol2=0.
2654	falb1=0.
2655	falb2=0.
2656	SELECT CASE(nsw)
2657	CASE(2)
2658	albsol1=albsol_dir(:,1)
2659	albsol2=albsol_dir(:,2)
2660	falb1=falb_dir(:,1,:)
2661	falb2=falb_dir(:,2,:)
2662	CASE(4)
2663	albsol1=albsol_dir(:,1)
2664	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3) &
2665	+albsol_dir(:,4)*SFRWL(4)
2666	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2667	falb1=falb_dir(:,1,:)
2668	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3) &
2669	+falb_dir(:,4,:)*SFRWL(4)
2670	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2671	CASE(6)
2672	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2) &
2673	+albsol_dir(:,3)*SFRWL(3)
2674	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2675	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5) &
2676	+albsol_dir(:,6)*SFRWL(6)
2677	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2678	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2) &
2679	+falb_dir(:,3,:)*SFRWL(3)
2680	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2681	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5) &
2682	+falb_dir(:,6,:)*SFRWL(6)
2683	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2684	END SELECt
2685	!albedo SB <<<
2686
2687
2688	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2689	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2690
2691	ENDIF
2692	! =================================================================== c
2693	! Calcul de Qsat
2694
2695	DO k = 1, klev
2696	DO i = 1, klon
2697	zx_t = t_seri(i,k)
2698	IF (thermcep) THEN
2699	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2700	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2701	zx_qs = MIN(0.5,zx_qs)
2702	zcor = 1./(1.-retv*zx_qs)
2703	zx_qs = zx_qs*zcor
2704	ELSE
2705	!! IF (zx_t.LT.t_coup) THEN !jyg
2706	IF (zx_t.LT.rtt) THEN !jyg
2707	zx_qs = qsats(zx_t)/pplay(i,k)
2708	ELSE
2709	zx_qs = qsatl(zx_t)/pplay(i,k)
2710	ENDIF
2711	ENDIF
2712	zqsat(i,k)=zx_qs
2713	ENDDO
2714	ENDDO
2715
2716	IF (prt_level.ge.1) THEN
2717	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2718	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2719	ENDIF
2720	!
2721	! Appeler la convection (au choix)
2722	!
2723	DO k = 1, klev
2724	DO i = 1, klon
2725	conv_q(i,k) = d_q_dyn(i,k) &
2726	+ d_q_vdf(i,k)/phys_tstep
2727	conv_t(i,k) = d_t_dyn(i,k) &
2728	+ d_t_vdf(i,k)/phys_tstep
2729	ENDDO
2730	ENDDO
2731	IF (check) THEN
2732	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2733	WRITE(lunout,*) "avantcon=", za
2734	ENDIF
2735	zx_ajustq = .FALSE.
2736	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2737	IF (zx_ajustq) THEN
2738	DO i = 1, klon
2739	z_avant(i) = 0.0
2740	ENDDO
2741	DO k = 1, klev
2742	DO i = 1, klon
2743	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2744	*(paprs(i,k)-paprs(i,k+1))/RG
2745	ENDDO
2746	ENDDO
2747	ENDIF
2748
2749	! Calcule de vitesse verticale a partir de flux de masse verticale
2750	DO k = 1, klev
2751	DO i = 1, klon
2752	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
2753	ENDDO
2754	ENDDO
2755
2756	IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2757	omega(igout, :)
2758	!
2759	! Appel de la convection tous les "cvpas"
2760	!
2761	!!jyg IF (MOD(itapcv,cvpas).EQ.0) THEN
2762	!! print *,' physiq : itapcv, cvpas, itap-1, cvpas_0 ', &
2763	!! itapcv, cvpas, itap-1, cvpas_0
2764	IF (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itap-1,cvpas_0).EQ.0) THEN
2765
2766	!
2767	! Mettre a zero des variables de sortie (pour securite)
2768	!
2769	pmflxr(:,:) = 0.
2770	pmflxs(:,:) = 0.
2771	wdtrainA(:,:) = 0.
2772	wdtrainS(:,:) = 0.
2773	wdtrainM(:,:) = 0.
2774	upwd(:,:) = 0.
2775	dnwd(:,:) = 0.
2776	ep(:,:) = 0.
2777	da(:,:)=0.
2778	mp(:,:)=0.
2779	wght_cvfd(:,:)=0.
2780	phi(:,:,:)=0.
2781	phi2(:,:,:)=0.
2782	epmlmMm(:,:,:)=0.
2783	eplaMm(:,:)=0.
2784	d1a(:,:)=0.
2785	dam(:,:)=0.
2786	elij(:,:,:)=0.
2787	ev(:,:)=0.
2788	qtaa(:,:)=0.
2789	clw(:,:)=0.
2790	sij(:,:,:)=0.
2791	!
2792	IF (iflag_con.EQ.1) THEN
2793	abort_message ='reactiver le call conlmd dans physiq.F'
2794	CALL abort_physic (modname,abort_message,1)
2795	! CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
2796	! . d_t_con, d_q_con,
2797	! . rain_con, snow_con, ibas_con, itop_con)
2798	ELSE IF (iflag_con.EQ.2) THEN
2799	CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
2800	conv_t, conv_q, -evap, omega, &
2801	d_t_con, d_q_con, rain_con, snow_con, &
2802	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2803	kcbot, kctop, kdtop, pmflxr, pmflxs)
2804	d_u_con = 0.
2805	d_v_con = 0.
2806
2807	WHERE (rain_con < 0.) rain_con = 0.
2808	WHERE (snow_con < 0.) snow_con = 0.
2809	DO i = 1, klon
2810	ibas_con(i) = klev+1 - kcbot(i)
2811	itop_con(i) = klev+1 - kctop(i)
2812	ENDDO
2813	ELSE IF (iflag_con.GE.3) THEN
2814	! nb of tracers for the KE convection:
2815	! MAF la partie traceurs est faite dans phytrac
2816	! on met ntra=1 pour limiter les appels mais on peut
2817	! supprimer les calculs / ftra.
2818	ntra = 1
2819
2820	!=======================================================================
2821	!ajout pour la parametrisation des poches froides: calcul de
2822	!t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
2823	IF (iflag_wake>=1) THEN
2824	DO k=1,klev
2825	DO i=1,klon
2826	t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
2827	q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
2828	t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2829	q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2830	ENDDO
2831	ENDDO
2832	ELSE
2833	t_w(:,:) = t_seri(:,:)
2834	q_w(:,:) = q_seri(:,:)
2835	t_x(:,:) = t_seri(:,:)
2836	q_x(:,:) = q_seri(:,:)
2837	ENDIF
2838	!
2839	!jyg<
2840	! Perform dry adiabatic adjustment on wake profile
2841	! The corresponding tendencies are added to the convective tendencies
2842	! after the call to the convective scheme.
2843	IF (iflag_wake>=1) then
2844	IF (iflag_adjwk >= 1) THEN
2845	limbas(:) = 1
2846	CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
2847	d_t_adjwk, d_q_adjwk)
2848	!
2849	DO k=1,klev
2850	DO i=1,klon
2851	IF (wake_s(i) .GT. 1.e-3) THEN
2852	t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
2853	q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
2854	d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
2855	d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
2856	ELSE
2857	d_deltat_ajs_cv(i,k) = 0.
2858	d_deltaq_ajs_cv(i,k) = 0.
2859	ENDIF
2860	ENDDO
2861	ENDDO
2862	IF (iflag_adjwk == 2) THEN
2863	CALL add_wake_tend &
2864	(d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, ddens0, wkoccur1, 'ajs_cv', abortphy)
2865	ENDIF ! (iflag_adjwk == 2)
2866	ENDIF ! (iflag_adjwk >= 1)
2867	ENDIF ! (iflag_wake>=1)
2868	!>jyg
2869	!
2870
2871	!! print *,'physiq. q_w(1,k), q_x(1,k) ', &
2872	!! (k, q_w(1,k), q_x(1,k),k=1,25)
2873
2874	!jyg<
2875	CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri, &
2876	alp_offset, it_wape_prescr, wape_prescr, fip_prescr, &
2877	ale_bl_prescr, alp_bl_prescr, &
2878	wake_pe, wake_fip, &
2879	Ale_bl, Ale_bl_trig, Alp_bl, &
2880	Ale, Alp , Ale_wake, Alp_wake)
2881	!>jyg
2882	!
2883	! sb, oct02:
2884	! Schema de convection modularise et vectorise:
2885	! (driver commun aux versions 3 et 4)
2886	!
2887	IF (ok_cvl) THEN ! new driver for convectL
2888	!
2889	!jyg<
2890	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2891	! Calculate the upmost level of deep convection loops: k_upper_cv
2892	! (near 22 km)
2893	k_upper_cv = klev
2894	!izero = klon/2+1/klon
2895	!DO k = klev,1,-1
2896	! IF (pphi(izero,k) > 22.e4) k_upper_cv = k
2897	!ENDDO
2898	! FH : nouveau calcul base sur un profil global sans quoi
2899	! le modele etait sensible au decoupage de domaines
2900	DO k = klev,1,-1
2901	IF (-7*log(presnivs(k)/presnivs(1)) > 25.) k_upper_cv = k
2902	ENDDO
2903	IF (prt_level .ge. 5) THEN
2904	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
2905	k_upper_cv
2906	ENDIF
2907	!
2908	!>jyg
2909	IF (type_trac == 'repr') THEN
2910	nbtr_tmp=ntra
2911	ELSE
2912	nbtr_tmp=nbtr
2913	ENDIF
2914	!jyg iflag_con est dans clesphys
2915	!c CALL concvl (iflag_con,iflag_clos,
2916	CALL concvl (iflag_clos, &
2917	phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
2918	t_w,q_w,wake_s, &
2919	u_seri,v_seri,tr_seri,nbtr_tmp, &
2920	ALE,ALP, &
2921	sig1,w01, &
2922	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2923	rain_con, snow_con, ibas_con, itop_con, sigd, &
2924	ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
2925	Ma,mipsh,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2926	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2927	! RomP >>>
2928	!! . pmflxr,pmflxs,da,phi,mp,
2929	!! . ftd,fqd,lalim_conv,wght_th)
2930	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,qtaa,clw,elij, &
2931	ftd,fqd,lalim_conv,wght_th, &
2932	ev, ep,epmlmMm,eplaMm, &
2933	wdtrainA, wdtrainS, wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2934	tau_cld_cv,coefw_cld_cv,epmax_diag)
2935
2936	! RomP <<<
2937
2938	!IM begin
2939	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2940	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2941	!IM end
2942	!IM cf. FH
2943	clwcon0=qcondc
2944	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2945	!
2946	!jyg<
2947	! If convective tendencies are too large, then call convection
2948	! every time step
2949	cvpas = cvpas_0
2950	DO k=1,k_upper_cv
2951	DO i=1,klon
2952	IF (d_t_con(i,k) > 6.721 .AND. d_t_con(i,k) < 6.722 .AND.&
2953	d_q_con(i,k) > -.0002171 .AND. d_q_con(i,k) < -.0002170) THEN
2954	dtcon_multistep_max = 3.
2955	dqcon_multistep_max = 0.02
2956	ENDIF
2957	ENDDO
2958	ENDDO
2959	!
2960	DO k=1,k_upper_cv
2961	DO i=1,klon
2962	!! IF (abs(d_t_con(i,k)) > 0.24 .OR. &
2963	!! abs(d_q_con(i,k)) > 2.e-2) THEN
2964	IF (abs(d_t_con(i,k)) > dtcon_multistep_max .OR. &
2965	abs(d_q_con(i,k)) > dqcon_multistep_max) THEN
2966	cvpas = 1
2967	!! print *,'physiq1, i,k,d_t_con(i,k),d_q_con(i,k) ', &
2968	!! i,k,d_t_con(i,k),d_q_con(i,k)
2969	ENDIF
2970	ENDDO
2971	ENDDO
2972	!!! Ligne a ne surtout pas remettre sans avoir murement reflechi (jyg)
2973	!!! call bcast(cvpas)
2974	!!! ------------------------------------------------------------
2975	!>jyg
2976	!
2977	DO i = 1, klon
2978	IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+cvpas
2979	ENDDO
2980	!
2981	!jyg<
2982	! Add the tendency due to the dry adjustment of the wake profile
2983	IF (iflag_wake>=1) THEN
2984	IF (iflag_adjwk == 2) THEN
2985	DO k=1,klev
2986	DO i=1,klon
2987	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
2988	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
2989	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
2990	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
2991	ENDDO
2992	ENDDO
2993	ENDIF ! (iflag_adjwk = 2)
2994	ENDIF ! (iflag_wake>=1)
2995	!>jyg
2996	!
2997	ELSE ! ok_cvl
2998
2999	! MAF conema3 ne contient pas les traceurs
3000	CALL conema3 (phys_tstep, &
3001	paprs,pplay,t_seri,q_seri, &
3002	u_seri,v_seri,tr_seri,ntra, &
3003	sig1,w01, &
3004	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
3005	rain_con, snow_con, ibas_con, itop_con, &
3006	upwd,dnwd,dnwd0,bas,top, &
3007	Ma,cape,tvp,rflag, &
3008	pbase &
3009	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
3010	,clwcon0)
3011
3012	ENDIF ! ok_cvl
3013
3014	!
3015	! Correction precip
3016	rain_con = rain_con * cvl_corr
3017	snow_con = snow_con * cvl_corr
3018	!
3019
3020	IF (.NOT. ok_gust) THEN
3021	do i = 1, klon
3022	wd(i)=0.0
3023	enddo
3024	ENDIF
3025
3026	! =================================================================== c
3027	! Calcul des proprietes des nuages convectifs
3028	!
3029
3030	! calcul des proprietes des nuages convectifs
3031	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
3032	IF (iflag_cld_cv == 0) THEN
3033	CALL clouds_gno &
3034	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
3035	ELSE
3036	CALL clouds_bigauss &
3037	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
3038	ENDIF
3039
3040
3041	! =================================================================== c
3042
3043	DO i = 1, klon
3044	itop_con(i) = min(max(itop_con(i),1),klev)
3045	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
3046	ENDDO
3047
3048	DO i = 1, klon
3049	ema_pcb(i) = paprs(i,ibas_con(i))
3050	ENDDO
3051	DO i = 1, klon
3052	! L'idicage de itop_con peut cacher un pb potentiel
3053	! FH sous la dictee de JYG, CR
3054	ema_pct(i) = paprs(i,itop_con(i)+1)
3055
3056	IF (itop_con(i).gt.klev-3) THEN
3057	IF (prt_level >= 9) THEN
3058	write(lunout,*)'La convection monte trop haut '
3059	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
3060	ENDIF
3061	ENDIF
3062	ENDDO
3063	ELSE IF (iflag_con.eq.0) THEN
3064	write(lunout,*) 'On n appelle pas la convection'
3065	clwcon0=0.
3066	rnebcon0=0.
3067	d_t_con=0.
3068	d_q_con=0.
3069	d_u_con=0.
3070	d_v_con=0.
3071	rain_con=0.
3072	snow_con=0.
3073	bas=1
3074	top=1
3075	ELSE
3076	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
3077	CALL abort_physic("physiq", "", 1)
3078	ENDIF
3079
3080	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
3081	! . d_u_con, d_v_con)
3082
3083	!jyg Reinitialize proba_notrig and itapcv when convection has been called
3084	proba_notrig(:) = 1.
3085	itapcv = 0
3086	ENDIF ! (MOD(itapcv,cvpas).EQ.0 .OR. MOD(itapcv,cvpas_0).EQ.0)
3087	!
3088	itapcv = itapcv+1
3089	!
3090	! Compter les steps ou cvpas=1
3091	IF (cvpas == 1) THEN
3092	Ncvpaseq1 = Ncvpaseq1+1
3093	ENDIF
3094	IF (mod(itap,1000) == 0) THEN
3095	print *,' physiq, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
3096	ENDIF
3097
3098	!!!jyg Appel diagnostique a add_phys_tend pour tester la conservation de
3099	!!! l'energie dans les courants satures.
3100	!! d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
3101	!! d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
3102	!! dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
3103	!! CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat, &
3104	!! dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
3105	!! itap, 1)
3106	!! call prt_enerbil('convection_sat',itap)
3107	!!
3108	!!
3109	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
3110	'convection',abortphy,flag_inhib_tend,itap,0)
3111	CALL prt_enerbil('convection',itap)
3112
3113	!-------------------------------------------------------------------------
3114
3115	IF (mydebug) THEN
3116	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3117	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3118	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3119	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3120	ENDIF
3121
3122	IF (check) THEN
3123	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3124	WRITE(lunout,*)"aprescon=", za
3125	zx_t = 0.0
3126	za = 0.0
3127	DO i = 1, klon
3128	za = za + cell_area(i)/REAL(klon)
3129	zx_t = zx_t + (rain_con(i)+ &
3130	snow_con(i))*cell_area(i)/REAL(klon)
3131	ENDDO
3132	zx_t = zx_t/za*phys_tstep
3133	WRITE(lunout,*)"Precip=", zx_t
3134	ENDIF
3135	IF (zx_ajustq) THEN
3136	DO i = 1, klon
3137	z_apres(i) = 0.0
3138	ENDDO
3139	DO k = 1, klev
3140	DO i = 1, klon
3141	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
3142	*(paprs(i,k)-paprs(i,k+1))/RG
3143	ENDDO
3144	ENDDO
3145	DO i = 1, klon
3146	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*phys_tstep) &
3147	/z_apres(i)
3148	ENDDO
3149	DO k = 1, klev
3150	DO i = 1, klon
3151	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
3152	z_factor(i).LT.(1.0-1.0E-08)) THEN
3153	q_seri(i,k) = q_seri(i,k) * z_factor(i)
3154	ENDIF
3155	ENDDO
3156	ENDDO
3157	ENDIF
3158	zx_ajustq=.FALSE.
3159
3160	!
3161	!==========================================================================
3162	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
3163	!pour la couche limite diffuse pour l instant
3164	!
3165	!
3166	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
3167	! il faut rajouter cette tendance calcul\'ee hors des poches
3168	! froides
3169	!
3170	IF (iflag_wake>=1) THEN
3171	!
3172	!
3173	! Call wakes every "wkpas" step
3174	!
3175	IF (MOD(itapwk,wkpas).EQ.0) THEN
3176	!
3177	DO k=1,klev
3178	DO i=1,klon
3179	dt_dwn(i,k) = ftd(i,k)
3180	dq_dwn(i,k) = fqd(i,k)
3181	M_dwn(i,k) = dnwd0(i,k)
3182	M_up(i,k) = upwd(i,k)
3183	dt_a(i,k) = d_t_con(i,k)/phys_tstep - ftd(i,k)
3184	dq_a(i,k) = d_q_con(i,k)/phys_tstep - fqd(i,k)
3185	ENDDO
3186	ENDDO
3187
3188	IF (iflag_wake==2) THEN
3189	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3190	DO k = 1,klev
3191	dt_dwn(:,k)= dt_dwn(:,k)+ &
3192	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
3193	dq_dwn(:,k)= dq_dwn(:,k)+ &
3194	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
3195	ENDDO
3196	ELSEIF (iflag_wake==3) THEN
3197	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
3198	DO k = 1,klev
3199	DO i=1,klon
3200	IF (rneb(i,k)==0.) THEN
3201	! On ne tient compte des tendances qu'en dehors des
3202	! nuages (c'est-\`a-dire a priri dans une region ou
3203	! l'eau se reevapore).
3204	dt_dwn(i,k)= dt_dwn(i,k)+ &
3205	ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
3206	dq_dwn(i,k)= dq_dwn(i,k)+ &
3207	ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
3208	ENDIF
3209	ENDDO
3210	ENDDO
3211	ENDIF
3212
3213	!
3214	!calcul caracteristiques de la poche froide
3215	CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
3216	t_seri, q_seri, omega, &
3217	dt_dwn, dq_dwn, M_dwn, M_up, &
3218	dt_a, dq_a, cv_gen, &
3219	sigd, cin, &
3220	wake_deltat, wake_deltaq, wake_s, awake_dens, wake_dens, &
3221	wake_dth, wake_h, &
3222	!! wake_pe, wake_fip, wake_gfl, &
3223	wake_pe, wake_fip_0, wake_gfl, & !! jyg
3224	d_t_wake, d_q_wake, &
3225	wake_k, t_x, q_x, &
3226	wake_omgbdth, wake_dp_omgb, &
3227	wake_dtKE, wake_dqKE, &
3228	wake_omg, wake_dp_deltomg, &
3229	wake_spread, wake_Cstar, d_deltat_wk_gw, &
3230	d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk)
3231	!
3232	!jyg Reinitialize itapwk when wakes have been called
3233	itapwk = 0
3234	ENDIF ! (MOD(itapwk,wkpas).EQ.0)
3235	!
3236	itapwk = itapwk+1
3237	!
3238	!-----------------------------------------------------------------------
3239	! ajout des tendances des poches froides
3240	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
3241	abortphy,flag_inhib_tend,itap,0)
3242	CALL prt_enerbil('wake',itap)
3243	!------------------------------------------------------------------------
3244
3245	! Increment Wake state variables
3246	IF (iflag_wake_tend .GT. 0.) THEN
3247
3248	CALL add_wake_tend &
3249	(d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_a_wk, d_dens_wk, wake_k, &
3250	'wake', abortphy)
3251	CALL prt_enerbil('wake',itap)
3252	ENDIF ! (iflag_wake_tend .GT. 0.)
3253	!
3254	IF (prt_level .GE. 10) THEN
3255	print *,' physiq, after calwake, wake_s: ',wake_s(:)
3256	print *,' physiq, after calwake, wake_deltat: ',wake_deltat(:,1)
3257	print *,' physiq, after calwake, wake_deltaq: ',wake_deltaq(:,1)
3258	ENDIF
3259
3260	IF (iflag_alp_wk_cond .GT. 0.) THEN
3261
3262	CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
3263	wake_fip)
3264	ELSE
3265	wake_fip(:) = wake_fip_0(:)
3266	ENDIF ! (iflag_alp_wk_cond .GT. 0.)
3267
3268	ENDIF ! (iflag_wake>=1)
3269	!
3270	!===================================================================
3271	! Convection seche (thermiques ou ajustement)
3272	!===================================================================
3273	!
3274	CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
3275	,seuil_inversion,weak_inversion,dthmin)
3276
3277
3278
3279	d_t_ajsb(:,:)=0.
3280	d_q_ajsb(:,:)=0.
3281	d_t_ajs(:,:)=0.
3282	d_u_ajs(:,:)=0.
3283	d_v_ajs(:,:)=0.
3284	d_q_ajs(:,:)=0.
3285	clwcon0th(:,:)=0.
3286	!
3287	! fm_therm(:,:)=0.
3288	! entr_therm(:,:)=0.
3289	! detr_therm(:,:)=0.
3290	!
3291	IF (prt_level>9) WRITE(lunout,*) &
3292	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
3293	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3294	IF (iflag_thermals<0) THEN
3295	! Rien
3296	! ====
3297	IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
3298
3299
3300	ELSE
3301
3302	! Thermiques
3303	! ==========
3304	IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
3305	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
3306
3307
3308	!cc nrlmd le 10/04/2012
3309	DO k=1,klev+1
3310	DO i=1,klon
3311	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
3312	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
3313	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
3314	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
3315	ENDDO
3316	ENDDO
3317	!cc fin nrlmd le 10/04/2012
3318
3319	IF (iflag_thermals>=1) THEN
3320	!jyg<
3321	!! IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3322	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3323	! Appel des thermiques avec les profils exterieurs aux poches
3324	DO k=1,klev
3325	DO i=1,klon
3326	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
3327	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
3328	u_therm(i,k) = u_seri(i,k)
3329	v_therm(i,k) = v_seri(i,k)
3330	ENDDO
3331	ENDDO
3332	ELSE
3333	! Appel des thermiques avec les profils moyens
3334	DO k=1,klev
3335	DO i=1,klon
3336	t_therm(i,k) = t_seri(i,k)
3337	q_therm(i,k) = q_seri(i,k)
3338	u_therm(i,k) = u_seri(i,k)
3339	v_therm(i,k) = v_seri(i,k)
3340	ENDDO
3341	ENDDO
3342	ENDIF
3343	!>jyg
3344	CALL calltherm(pdtphys &
3345	,pplay,paprs,pphi,weak_inversion &
3346	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
3347	,u_therm,v_therm,t_therm,q_therm,zqsat,debut & !jyg
3348	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
3349	,fm_therm,entr_therm,detr_therm &
3350	,zqasc,clwcon0th,lmax_th,ratqscth &
3351	,ratqsdiff,zqsatth &
3352	!on rajoute ale et alp, et les
3353	!caracteristiques de la couche alim
3354	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
3355	,ztv,zpspsk,ztla,zthl &
3356	!cc nrlmd le 10/04/2012
3357	,pbl_tke_input,pctsrf,omega,cell_area &
3358	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
3359	,n2,s2,ale_bl_stat &
3360	,therm_tke_max,env_tke_max &
3361	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
3362	,alp_bl_conv,alp_bl_stat &
3363	!cc fin nrlmd le 10/04/2012
3364	,zqla,ztva )
3365	!
3366	!jyg<
3367	!!jyg IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
3368	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
3369	! Si les thermiques ne sont presents que hors des
3370	! poches, la tendance moyenne associ\'ee doit etre
3371	! multipliee par la fraction surfacique qu'ils couvrent.
3372	DO k=1,klev
3373	DO i=1,klon
3374	!
3375	d_deltat_the(i,k) = - d_t_ajs(i,k)
3376	d_deltaq_the(i,k) = - d_q_ajs(i,k)
3377	!
3378	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
3379	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
3380	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
3381	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
3382	!
3383	ENDDO
3384	ENDDO
3385	!
3386	IF (ok_bug_split_th) THEN
3387	CALL add_wake_tend &
3388	(d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wkoccur1, 'the', abortphy)
3389	ELSE
3390	CALL add_wake_tend &
3391	(d_deltat_the, d_deltaq_the, dsig0, ddens0, ddens0, wake_k, 'the', abortphy)
3392	ENDIF
3393	CALL prt_enerbil('the',itap)
3394	!
3395	ENDIF ! (mod(iflag_pbl_split/10,10) .GE. 1)
3396	!
3397	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
3398	dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
3399	CALL prt_enerbil('thermals',itap)
3400	!
3401	!
3402	CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area, &
3403	cin, s2, n2, &
3404	ale_bl_trig, ale_bl_stat, ale_bl, &
3405	alp_bl, alp_bl_stat, &
3406	proba_notrig, random_notrig, cv_gen)
3407	!>jyg
3408
3409	! ------------------------------------------------------------------
3410	! Transport de la TKE par les panaches thermiques.
3411	! FH : 2010/02/01
3412	! if (iflag_pbl.eq.10) then
3413	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3414	! s rg,paprs,pbl_tke)
3415	! endif
3416	! -------------------------------------------------------------------
3417
3418	DO i=1,klon
3419	! zmax_th(i)=pphi(i,lmax_th(i))/rg
3420	!CR:04/05/12:correction calcul zmax
3421	zmax_th(i)=zmax0(i)
3422	ENDDO
3423
3424	ENDIF
3425
3426	! Ajustement sec
3427	! ==============
3428
3429	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3430	! a partir du sommet des thermiques.
3431	! Dans le cas contraire, on demarre au niveau 1.
3432
3433	IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
3434
3435	IF (iflag_thermals.eq.0) THEN
3436	IF (prt_level>9) WRITE(lunout,*)'ajsec'
3437	limbas(:)=1
3438	ELSE
3439	limbas(:)=lmax_th(:)
3440	ENDIF
3441
3442	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3443	! pour des test de convergence numerique.
3444	! Le nouveau ajsec est a priori mieux, meme pour le cas
3445	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3446	! non nulles numeriquement pour des mailles non concernees.
3447
3448	IF (iflag_thermals==0) THEN
3449	! Calling adjustment alone (but not the thermal plume model)
3450	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3451	, d_t_ajsb, d_q_ajsb)
3452	ELSE IF (iflag_thermals>0) THEN
3453	! Calling adjustment above the top of thermal plumes
3454	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3455	, d_t_ajsb, d_q_ajsb)
3456	ENDIF
3457
3458	!--------------------------------------------------------------------
3459	! ajout des tendances de l'ajustement sec ou des thermiques
3460	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
3461	'ajsb',abortphy,flag_inhib_tend,itap,0)
3462	CALL prt_enerbil('ajsb',itap)
3463	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3464	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3465
3466	!---------------------------------------------------------------------
3467
3468	ENDIF
3469
3470	ENDIF
3471	!
3472	!===================================================================
3473	! Computation of ratqs, the width (normalized) of the subrid scale
3474	! water distribution
3475	CALL calcratqs(klon,klev,prt_level,lunout, &
3476	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3477	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
3478	tau_ratqs,fact_cldcon, &
3479	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3480	paprs,pplay,q_seri,zqsat,fm_therm, &
3481	ratqs,ratqsc)
3482
3483
3484	!
3485	! Appeler le processus de condensation a grande echelle
3486	! et le processus de precipitation
3487	!-------------------------------------------------------------------------
3488	IF (prt_level .GE.10) THEN
3489	print *,'itap, ->fisrtilp ',itap
3490	ENDIF
3491	!
3492	CALL fisrtilp(phys_tstep,paprs,pplay, &
3493	t_seri, q_seri,ptconv,ratqs, &
3494	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3495	rain_lsc, snow_lsc, &
3496	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3497	frac_impa, frac_nucl, beta_prec_fisrt, &
3498	prfl, psfl, rhcl, &
3499	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3500	iflag_ice_thermo)
3501	!
3502	WHERE (rain_lsc < 0) rain_lsc = 0.
3503	WHERE (snow_lsc < 0) snow_lsc = 0.
3504
3505	!+JLD
3506	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3507	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
3508	! & ,rain_lsc,snow_lsc
3509	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3510	!-JLD
3511	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
3512	'lsc',abortphy,flag_inhib_tend,itap,0)
3513	CALL prt_enerbil('lsc',itap)
3514	rain_num(:)=0.
3515	DO k = 1, klev
3516	DO i = 1, klon
3517	IF (ql_seri(i,k)>oliqmax) THEN
3518	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3519	ql_seri(i,k)=oliqmax
3520	ENDIF
3521	ENDDO
3522	ENDDO
3523	IF (nqo==3) THEN
3524	DO k = 1, klev
3525	DO i = 1, klon
3526	IF (qs_seri(i,k)>oicemax) THEN
3527	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3528	qs_seri(i,k)=oicemax
3529	ENDIF
3530	ENDDO
3531	ENDDO
3532	ENDIF
3533
3534	!---------------------------------------------------------------------------
3535	DO k = 1, klev
3536	DO i = 1, klon
3537	cldfra(i,k) = rneb(i,k)
3538	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3539	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3540	ENDDO
3541	ENDDO
3542	IF (check) THEN
3543	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3544	WRITE(lunout,*)"apresilp=", za
3545	zx_t = 0.0
3546	za = 0.0
3547	DO i = 1, klon
3548	za = za + cell_area(i)/REAL(klon)
3549	zx_t = zx_t + (rain_lsc(i) &
3550	+ snow_lsc(i))*cell_area(i)/REAL(klon)
3551	ENDDO
3552	zx_t = zx_t/za*phys_tstep
3553	WRITE(lunout,*)"Precip=", zx_t
3554	ENDIF
3555
3556	IF (mydebug) THEN
3557	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3558	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3559	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3560	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3561	ENDIF
3562
3563	!
3564	!-------------------------------------------------------------------
3565	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3566	!-------------------------------------------------------------------
3567
3568	! 1. NUAGES CONVECTIFS
3569	!
3570	!IM cf FH
3571	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3572	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3573	snow_tiedtke=0.
3574	! print*,'avant calcul de la pseudo precip '
3575	! print*,'iflag_cld_th',iflag_cld_th
3576	IF (iflag_cld_th.eq.-1) THEN
3577	rain_tiedtke=rain_con
3578	ELSE
3579	! print*,'calcul de la pseudo precip '
3580	rain_tiedtke=0.
3581	! print*,'calcul de la pseudo precip 0'
3582	DO k=1,klev
3583	DO i=1,klon
3584	IF (d_q_con(i,k).lt.0.) THEN
3585	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3586	*(paprs(i,k)-paprs(i,k+1))/rg
3587	ENDIF
3588	ENDDO
3589	ENDDO
3590	ENDIF
3591	!
3592	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3593	!
3594
3595	! Nuages diagnostiques pour Tiedtke
3596	CALL diagcld1(paprs,pplay, &
3597	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3598	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3599	diafra,dialiq)
3600	DO k = 1, klev
3601	DO i = 1, klon
3602	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3603	cldliq(i,k) = dialiq(i,k)
3604	cldfra(i,k) = diafra(i,k)
3605	ENDIF
3606	ENDDO
3607	ENDDO
3608
3609	ELSE IF (iflag_cld_th.ge.3) THEN
3610	! On prend pour les nuages convectifs le max du calcul de la
3611	! convection et du calcul du pas de temps precedent diminue d'un facteur
3612	! facttemps
3613	facteur = pdtphys *facttemps
3614	DO k=1,klev
3615	DO i=1,klon
3616	rnebcon(i,k)=rnebcon(i,k)*facteur
3617	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
3618	rnebcon(i,k)=rnebcon0(i,k)
3619	clwcon(i,k)=clwcon0(i,k)
3620	ENDIF
3621	ENDDO
3622	ENDDO
3623
3624	! On prend la somme des fractions nuageuses et des contenus en eau
3625
3626	IF (iflag_cld_th>=5) THEN
3627
3628	DO k=1,klev
3629	ptconvth(:,k)=fm_therm(:,k+1)>0.
3630	ENDDO
3631
3632	IF (iflag_coupl==4) THEN
3633
3634	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3635	! convectives et lsc dans la partie des thermiques
3636	! Le controle par iflag_coupl est peut etre provisoire.
3637	DO k=1,klev
3638	DO i=1,klon
3639	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
3640	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3641	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3642	ELSE IF (ptconv(i,k)) THEN
3643	cldfra(i,k)=rnebcon(i,k)
3644	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3645	ENDIF
3646	ENDDO
3647	ENDDO
3648
3649	ELSE IF (iflag_coupl==5) THEN
3650	DO k=1,klev
3651	DO i=1,klon
3652	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3653	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3654	ENDDO
3655	ENDDO
3656
3657	ELSE
3658
3659	! Si on est sur un point touche par la convection
3660	! profonde et pas par les thermiques, on prend la
3661	! couverture nuageuse et l'eau nuageuse de la convection
3662	! profonde.
3663
3664	!IM/FH: 2011/02/23
3665	! definition des points sur lesquels ls thermiques sont actifs
3666
3667	DO k=1,klev
3668	DO i=1,klon
3669	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
3670	cldfra(i,k)=rnebcon(i,k)
3671	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3672	ENDIF
3673	ENDDO
3674	ENDDO
3675
3676	ENDIF
3677
3678	ELSE
3679
3680	! Ancienne version
3681	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3682	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3683	ENDIF
3684
3685	ENDIF
3686
3687	! plulsc(:)=0.
3688	! do k=1,klev,-1
3689	! do i=1,klon
3690	! zzz=prfl(:,k)+psfl(:,k)
3691	! if (.not.ptconvth.zzz.gt.0.)
3692	! enddo prfl, psfl,
3693	! enddo
3694	!
3695	! 2. NUAGES STARTIFORMES
3696	!
3697	IF (ok_stratus) THEN
3698	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3699	DO k = 1, klev
3700	DO i = 1, klon
3701	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3702	cldliq(i,k) = dialiq(i,k)
3703	cldfra(i,k) = diafra(i,k)
3704	ENDIF
3705	ENDDO
3706	ENDDO
3707	ENDIF
3708	!
3709	! Precipitation totale
3710	!
3711	DO i = 1, klon
3712	rain_fall(i) = rain_con(i) + rain_lsc(i)
3713	snow_fall(i) = snow_con(i) + snow_lsc(i)
3714	ENDDO
3715	!
3716	! Calculer l'humidite relative pour diagnostique
3717	!
3718	DO k = 1, klev
3719	DO i = 1, klon
3720	zx_t = t_seri(i,k)
3721	IF (thermcep) THEN
3722	!! if (iflag_ice_thermo.eq.0) then !jyg
3723	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3724	!! else !jyg
3725	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
3726	!! endif !jyg
3727	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3728	zx_qs = MIN(0.5,zx_qs)
3729	zcor = 1./(1.-retv*zx_qs)
3730	zx_qs = zx_qs*zcor
3731	ELSE
3732	!! IF (zx_t.LT.t_coup) THEN !jyg
3733	IF (zx_t.LT.rtt) THEN !jyg
3734	zx_qs = qsats(zx_t)/pplay(i,k)
3735	ELSE
3736	zx_qs = qsatl(zx_t)/pplay(i,k)
3737	ENDIF
3738	ENDIF
3739	zx_rh(i,k) = q_seri(i,k)/zx_qs
3740	IF (iflag_ice_thermo .GT. 0) THEN
3741	zx_rhl(i,k) = q_seri(i,k)/(qsatl(zx_t)/pplay(i,k))
3742	zx_rhi(i,k) = q_seri(i,k)/(qsats(zx_t)/pplay(i,k))
3743	ENDIF
3744	zqsat(i,k)=zx_qs
3745	ENDDO
3746	ENDDO
3747
3748	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3749	! equivalente a 2m (tpote) pour diagnostique
3750	!
3751	DO i = 1, klon
3752	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3753	IF (thermcep) THEN
3754	IF(zt2m(i).LT.RTT) then
3755	Lheat=RLSTT
3756	ELSE
3757	Lheat=RLVTT
3758	ENDIF
3759	ELSE
3760	IF (zt2m(i).LT.RTT) THEN
3761	Lheat=RLSTT
3762	ELSE
3763	Lheat=RLVTT
3764	ENDIF
3765	ENDIF
3766	tpote(i) = tpot(i)* &
3767	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3768	ENDDO
3769
3770	IF (type_trac == 'inca') THEN
3771	#ifdef INCA
3772	CALL VTe(VTphysiq)
3773	CALL VTb(VTinca)
3774	calday = REAL(days_elapsed + 1) + jH_cur
3775
3776	CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
3777	CALL AEROSOL_METEO_CALC( &
3778	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3779	prfl,psfl,pctsrf,cell_area, &
3780	latitude_deg,longitude_deg,u10m,v10m)
3781
3782	zxsnow_dummy(:) = 0.0
3783
3784	CALL chemhook_begin (calday, &
3785	days_elapsed+1, &
3786	jH_cur, &
3787	pctsrf(1,1), &
3788	latitude_deg, &
3789	longitude_deg, &
3790	cell_area, &
3791	paprs, &
3792	pplay, &
3793	coefh(1:klon,1:klev,is_ave), &
3794	pphi, &
3795	t_seri, &
3796	u, &
3797	v, &
3798	rot, &
3799	wo(:, :, 1), &
3800	q_seri, &
3801	zxtsol, &
3802	zt2m, &
3803	zxsnow_dummy, &
3804	solsw, &
3805	albsol1, &
3806	rain_fall, &
3807	snow_fall, &
3808	itop_con, &
3809	ibas_con, &
3810	cldfra, &
3811	nbp_lon, &
3812	nbp_lat-1, &
3813	tr_seri, &
3814	ftsol, &
3815	paprs, &
3816	cdragh, &
3817	cdragm, &
3818	pctsrf, &
3819	pdtphys, &
3820	itap)
3821
3822	CALL VTe(VTinca)
3823	CALL VTb(VTphysiq)
3824	#endif
3825	ENDIF !type_trac = inca
3826	IF (type_trac == 'repr') THEN
3827	#ifdef REPROBUS
3828	!CALL chemtime_rep(itap+itau_phy-1, date0, dtime, itap)
3829	CALL chemtime_rep(itap+itau_phy-1, date0, phys_tstep, itap)
3830	#endif
3831	ENDIF
3832
3833	!
3834	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3835	!
3836	IF (MOD(itaprad,radpas).EQ.0) THEN
3837
3838	!
3839	!jq - introduce the aerosol direct and first indirect radiative forcings
3840	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3841	IF (flag_aerosol .GT. 0) THEN
3842	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3843	IF (.NOT. aerosol_couple) THEN
3844	!
3845	CALL readaerosol_optic( &
3846	debut, flag_aerosol, itap, jD_cur-jD_ref, &
3847	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3848	mass_solu_aero, mass_solu_aero_pi, &
3849	tau_aero, piz_aero, cg_aero, &
3850	tausum_aero, tau3d_aero)
3851	ENDIF
3852	ELSE ! RRTM radiation
3853	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3854	abort_message='config_inca=aero et rrtm=1 impossible'
3855	CALL abort_physic(modname,abort_message,1)
3856	ELSE
3857	!
3858	#ifdef CPP_RRTM
3859	IF (NSW.EQ.6) THEN
3860	!--new aerosol properties SW and LW
3861	!
3862	#ifdef CPP_Dust
3863	!--SPL aerosol model
3864	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
3865	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3866	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3867	tausum_aero, tau3d_aero)
3868	#else
3869	!--climatologies or INCA aerosols
3870	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
3871	flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
3872	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3873	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3874	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3875	tausum_aero, drytausum_aero, tau3d_aero)
3876	#endif
3877
3878	IF (flag_aerosol .EQ. 7) THEN
3879	CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg, &
3880	tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
3881	ENDIF
3882
3883	!
3884	ELSE IF (NSW.EQ.2) THEN
3885	!--for now we use the old aerosol properties
3886	!
3887	CALL readaerosol_optic( &
3888	debut, flag_aerosol, itap, jD_cur-jD_ref, &
3889	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3890	mass_solu_aero, mass_solu_aero_pi, &
3891	tau_aero, piz_aero, cg_aero, &
3892	tausum_aero, tau3d_aero)
3893	!
3894	!--natural aerosols
3895	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
3896	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
3897	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
3898	!--all aerosols
3899	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
3900	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
3901	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
3902	!
3903	!--no LW optics
3904	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3905	!
3906	ELSE
3907	abort_message='Only NSW=2 or 6 are possible with ' &
3908	// 'aerosols and iflag_rrtm=1'
3909	CALL abort_physic(modname,abort_message,1)
3910	ENDIF
3911	#else
3912	abort_message='You should compile with -rrtm if running ' &
3913	// 'with iflag_rrtm=1'
3914	CALL abort_physic(modname,abort_message,1)
3915	#endif
3916	!
3917	ENDIF
3918	ENDIF
3919	ELSE !--flag_aerosol = 0
3920	tausum_aero(:,:,:) = 0.
3921	drytausum_aero(:,:) = 0.
3922	mass_solu_aero(:,:) = 0.
3923	mass_solu_aero_pi(:,:) = 0.
3924	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3925	tau_aero(:,:,:,:) = 1.e-15
3926	piz_aero(:,:,:,:) = 1.
3927	cg_aero(:,:,:,:) = 0.
3928	ELSE
3929	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
3930	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3931	piz_aero_sw_rrtm(:,:,:,:) = 1.0
3932	cg_aero_sw_rrtm(:,:,:,:) = 0.0
3933	ENDIF
3934	ENDIF
3935	!
3936	!--WMO criterion to determine tropopause
3937	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3938	!
3939	!--STRAT AEROSOL
3940	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
3941	IF (flag_aerosol_strat.GT.0) THEN
3942	IF (prt_level .GE.10) THEN
3943	PRINT *,'appel a readaerosolstrat', mth_cur
3944	ENDIF
3945	IF (iflag_rrtm.EQ.0) THEN
3946	IF (flag_aerosol_strat.EQ.1) THEN
3947	CALL readaerosolstrato(debut)
3948	ELSE
3949	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
3950	CALL abort_physic(modname,abort_message,1)
3951	ENDIF
3952	ELSE
3953	#ifdef CPP_RRTM
3954	#ifndef CPP_StratAer
3955	!--prescribed strat aerosols
3956	!--only in the case of non-interactive strat aerosols
3957	IF (flag_aerosol_strat.EQ.1) THEN
3958	CALL readaerosolstrato1_rrtm(debut)
3959	ELSEIF (flag_aerosol_strat.EQ.2) THEN
3960	CALL readaerosolstrato2_rrtm(debut, ok_volcan)
3961	ELSE
3962	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
3963	CALL abort_physic(modname,abort_message,1)
3964	ENDIF
3965	#endif
3966	#else
3967	abort_message='You should compile with -rrtm if running ' &
3968	// 'with iflag_rrtm=1'
3969	CALL abort_physic(modname,abort_message,1)
3970	#endif
3971	ENDIF
3972	ELSE
3973	tausum_aero(:,:,id_STRAT_phy) = 0.
3974	ENDIF
3975	!
3976	#ifdef CPP_RRTM
3977	#ifdef CPP_StratAer
3978	!--compute stratospheric mask
3979	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3980	!--interactive strat aerosols
3981	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
3982	#endif
3983	#endif
3984	!--fin STRAT AEROSOL
3985	!
3986
3987	! Calculer les parametres optiques des nuages et quelques
3988	! parametres pour diagnostiques:
3989	!
3990	IF (aerosol_couple.AND.config_inca=='aero') THEN
3991	mass_solu_aero(:,:) = ccm(:,:,1)
3992	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3993	ENDIF
3994
3995	IF (ok_newmicro) then
3996	IF (iflag_rrtm.NE.0) THEN
3997	#ifdef CPP_RRTM
3998	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3999	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
4000	// 'pour ok_cdnc'
4001	CALL abort_physic(modname,abort_message,1)
4002	ENDIF
4003	#else
4004
4005	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
4006	CALL abort_physic(modname,abort_message,1)
4007	#endif
4008	ENDIF
4009	CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
4010	paprs, pplay, t_seri, cldliq, cldfra, &
4011	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
4012	flwp, fiwp, flwc, fiwc, &
4013	mass_solu_aero, mass_solu_aero_pi, &
4014	cldtaupi, re, fl, ref_liq, ref_ice, &
4015	ref_liq_pi, ref_ice_pi)
4016	ELSE
4017	CALL nuage (paprs, pplay, &
4018	t_seri, cldliq, cldfra, cldtau, cldemi, &
4019	cldh, cldl, cldm, cldt, cldq, &
4020	ok_aie, &
4021	mass_solu_aero, mass_solu_aero_pi, &
4022	bl95_b0, bl95_b1, &
4023	cldtaupi, re, fl)
4024	ENDIF
4025	!
4026	!IM betaCRF
4027	!
4028	cldtaurad = cldtau
4029	cldtaupirad = cldtaupi
4030	cldemirad = cldemi
4031	cldfrarad = cldfra
4032
4033	!
4034	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
4035	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
4036	!
4037	! global
4038	!
4039	!IM 251017 begin
4040	! print*,'physiq betaCRF global zdtime=',zdtime
4041	!IM 251017 end
4042	DO k=1, klev
4043	DO i=1, klon
4044	IF (pplay(i,k).GE.pfree) THEN
4045	beta(i,k) = beta_pbl
4046	ELSE
4047	beta(i,k) = beta_free
4048	ENDIF
4049	IF (mskocean_beta) THEN
4050	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4051	ENDIF
4052	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4053	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4054	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4055	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4056	ENDDO
4057	ENDDO
4058	!
4059	ELSE
4060	!
4061	! regional
4062	!
4063	DO k=1, klev
4064	DO i=1,klon
4065	!
4066	IF (longitude_deg(i).ge.lon1_beta.AND. &
4067	longitude_deg(i).le.lon2_beta.AND. &
4068	latitude_deg(i).le.lat1_beta.AND. &
4069	latitude_deg(i).ge.lat2_beta) THEN
4070	IF (pplay(i,k).GE.pfree) THEN
4071	beta(i,k) = beta_pbl
4072	ELSE
4073	beta(i,k) = beta_free
4074	ENDIF
4075	IF (mskocean_beta) THEN
4076	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4077	ENDIF
4078	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4079	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4080	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4081	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4082	ENDIF
4083	!
4084	ENDDO
4085	ENDDO
4086	!
4087	ENDIF
4088
4089	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
4090	IF (ok_chlorophyll) THEN
4091	print*,"-- reading chlorophyll"
4092	CALL readchlorophyll(debut)
4093	ENDIF
4094
4095	!--if ok_suntime_rrtm we use ancillay data for RSUN
4096	!--previous values are therefore overwritten
4097	!--this is needed for CMIP6 runs
4098	!--and only possible for new radiation scheme
4099	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
4100	#ifdef CPP_RRTM
4101	CALL read_rsun_rrtm(debut)
4102	#endif
4103	ENDIF
4104
4105	IF (mydebug) THEN
4106	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4107	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4108	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4109	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4110	ENDIF
4111
4112	!
4113	!sonia : If Iflag_radia >=2, pertubation of some variables
4114	!input to radiation (DICE)
4115	!
4116	IF (iflag_radia .ge. 2) THEN
4117	zsav_tsol (:) = zxtsol(:)
4118	CALL perturb_radlwsw(zxtsol,iflag_radia)
4119	ENDIF
4120
4121	IF (aerosol_couple.AND.config_inca=='aero') THEN
4122	#ifdef INCA
4123	CALL radlwsw_inca &
4124	(chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
4125	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
4126	size(wo,3), wo, &
4127	cldfrarad, cldemirad, cldtaurad, &
4128	heat,heat0,cool,cool0,albpla, &
4129	topsw,toplw,solsw,sollw, &
4130	sollwdown, &
4131	topsw0,toplw0,solsw0,sollw0, &
4132	lwdn0, lwdn, lwup0, lwup, &
4133	swdn0, swdn, swup0, swup, &
4134	ok_ade, ok_aie, &
4135	tau_aero, piz_aero, cg_aero, &
4136	topswad_aero, solswad_aero, &
4137	topswad0_aero, solswad0_aero, &
4138	topsw_aero, topsw0_aero, &
4139	solsw_aero, solsw0_aero, &
4140	cldtaupirad, &
4141	topswai_aero, solswai_aero)
4142	#endif
4143	ELSE
4144	!
4145	!IM calcul radiatif pour le cas actuel
4146	!
4147	RCO2 = RCO2_act
4148	RCH4 = RCH4_act
4149	RN2O = RN2O_act
4150	RCFC11 = RCFC11_act
4151	RCFC12 = RCFC12_act
4152	!
4153	!--interactive CO2 in ppm from carbon cycle
4154	IF (carbon_cycle_rad.AND..NOT.debut) THEN
4155	RCO2=RCO2_glo
4156	ENDIF
4157	!
4158	IF (prt_level .GE.10) THEN
4159	print *,' ->radlwsw, number 1 '
4160	ENDIF
4161	!
4162	CALL radlwsw &
4163	(dist, rmu0, fract, &
4164	!albedo SB >>>
4165	! paprs, pplay,zxtsol,albsol1, albsol2, &
4166	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4167	!albedo SB <<<
4168	t_seri,q_seri,wo, &
4169	cldfrarad, cldemirad, cldtaurad, &
4170	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, &
4171	flag_aerosol, &
4172	flag_aerosol_strat, flag_aer_feedback, &
4173	tau_aero, piz_aero, cg_aero, &
4174	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4175	! Rajoute par OB pour RRTM
4176	tau_aero_lw_rrtm, &
4177	cldtaupirad, &
4178	! zqsat, flwcrad, fiwcrad, &
4179	zqsat, flwc, fiwc, &
4180	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4181	heat,heat0,cool,cool0,albpla, &
4182	heat_volc,cool_volc, &
4183	topsw,toplw,solsw,solswfdiff,sollw, &
4184	sollwdown, &
4185	topsw0,toplw0,solsw0,sollw0, &
4186	lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, &
4187	swdnc0, swdn0, swdn, swupc0, swup0, swup, &
4188	topswad_aero, solswad_aero, &
4189	topswai_aero, solswai_aero, &
4190	topswad0_aero, solswad0_aero, &
4191	topsw_aero, topsw0_aero, &
4192	solsw_aero, solsw0_aero, &
4193	topswcf_aero, solswcf_aero, &
4194	!-C. Kleinschmitt for LW diagnostics
4195	toplwad_aero, sollwad_aero,&
4196	toplwai_aero, sollwai_aero, &
4197	toplwad0_aero, sollwad0_aero,&
4198	!-end
4199	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4200	ZSWFT0_i, ZFSDN0, ZFSUP0)
4201
4202	!lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
4203	!schemes
4204	toplw = toplw + betalwoff * (toplw0 - toplw)
4205	sollw = sollw + betalwoff * (sollw0 - sollw)
4206	lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
4207	lwup = lwup + betalwoff * (lwup0 - lwup)
4208	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
4209	sollwdown(:))
4210	cool = cool + betalwoff * (cool0 - cool)
4211
4212	#ifndef CPP_XIOS
4213	!--OB 30/05/2016 modified 21/10/2016
4214	!--here we return swaero_diag and dryaod_diag to FALSE
4215	!--and histdef will switch it back to TRUE if necessary
4216	!--this is necessary to get the right swaero at first step
4217	!--but only in the case of no XIOS as XIOS is covered elsewhere
4218	IF (debut) swaerofree_diag = .FALSE.
4219	IF (debut) swaero_diag = .FALSE.
4220	IF (debut) dryaod_diag = .FALSE.
4221	!--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
4222	!--as for swaero_diag, see above
4223	IF (debut) ok_4xCO2atm = .FALSE.
4224
4225	!
4226	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
4227	!IM des taux doit etre different du taux actuel
4228	!IM Par defaut on a les taux perturbes egaux aux taux actuels
4229	!
4230	IF (RCO2_per.NE.RCO2_act.OR. &
4231	RCH4_per.NE.RCH4_act.OR. &
4232	RN2O_per.NE.RN2O_act.OR. &
4233	RCFC11_per.NE.RCFC11_act.OR. &
4234	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
4235	#endif
4236	!
4237	IF (ok_4xCO2atm) THEN
4238	!
4239	RCO2 = RCO2_per
4240	RCH4 = RCH4_per
4241	RN2O = RN2O_per
4242	RCFC11 = RCFC11_per
4243	RCFC12 = RCFC12_per
4244	!
4245	IF (prt_level .GE.10) THEN
4246	print *,' ->radlwsw, number 2 '
4247	ENDIF
4248	!
4249	CALL radlwsw &
4250	(dist, rmu0, fract, &
4251	!albedo SB >>>
4252	! paprs, pplay,zxtsol,albsol1, albsol2, &
4253	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4254	!albedo SB <<<
4255	t_seri,q_seri,wo, &
4256	cldfrarad, cldemirad, cldtaurad, &
4257	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, &
4258	flag_aerosol, &
4259	flag_aerosol_strat, flag_aer_feedback, &
4260	tau_aero, piz_aero, cg_aero, &
4261	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4262	! Rajoute par OB pour RRTM
4263	tau_aero_lw_rrtm, &
4264	cldtaupi, &
4265	! zqsat, flwcrad, fiwcrad, &
4266	zqsat, flwc, fiwc, &
4267	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4268	heatp,heat0p,coolp,cool0p,albplap, &
4269	heat_volc,cool_volc, &
4270	topswp,toplwp,solswp,solswfdiffp,sollwp, &
4271	sollwdownp, &
4272	topsw0p,toplw0p,solsw0p,sollw0p, &
4273	lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp, &
4274	swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
4275	topswad_aerop, solswad_aerop, &
4276	topswai_aerop, solswai_aerop, &
4277	topswad0_aerop, solswad0_aerop, &
4278	topsw_aerop, topsw0_aerop, &
4279	solsw_aerop, solsw0_aerop, &
4280	topswcf_aerop, solswcf_aerop, &
4281	!-C. Kleinschmitt for LW diagnostics
4282	toplwad_aerop, sollwad_aerop,&
4283	toplwai_aerop, sollwai_aerop, &
4284	toplwad0_aerop, sollwad0_aerop,&
4285	!-end
4286	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4287	ZSWFT0_i, ZFSDN0, ZFSUP0)
4288	endif !ok_4xCO2atm
4289	ENDIF ! aerosol_couple
4290	itaprad = 0
4291	!
4292	! If Iflag_radia >=2, reset pertubed variables
4293	!
4294	IF (iflag_radia .ge. 2) THEN
4295	zxtsol(:) = zsav_tsol (:)
4296	ENDIF
4297	ENDIF ! MOD(itaprad,radpas)
4298	itaprad = itaprad + 1
4299
4300	IF (iflag_radia.eq.0) THEN
4301	IF (prt_level.ge.9) THEN
4302	PRINT *,'--------------------------------------------------'
4303	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4304	PRINT *,'>>>> heat et cool mis a zero '
4305	PRINT *,'--------------------------------------------------'
4306	ENDIF
4307	heat=0.
4308	cool=0.
4309	sollw=0. ! MPL 01032011
4310	solsw=0.
4311	radsol=0.
4312	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4313	swup0=0.
4314	lwup=0.
4315	lwup0=0.
4316	lwdn=0.
4317	lwdn0=0.
4318	ENDIF
4319
4320	!
4321	! Calculer radsol a l'exterieur de radlwsw
4322	! pour prendre en compte le cycle diurne
4323	! recode par Olivier Boucher en sept 2015
4324	!
4325	radsol=solsw*swradcorr+sollw
4326
4327	IF (ok_4xCO2atm) THEN
4328	radsolp=solswp*swradcorr+sollwp
4329	ENDIF
4330
4331	!
4332	! Ajouter la tendance des rayonnements (tous les pas)
4333	! avec une correction pour le cycle diurne dans le SW
4334	!
4335
4336	DO k=1, klev
4337	d_t_swr(:,k)=swradcorr(:)heat(:,k)phys_tstep/RDAY
4338	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)phys_tstep/RDAY
4339	d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
4340	d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
4341	ENDDO
4342
4343	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4344	CALL prt_enerbil('SW',itap)
4345	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4346	CALL prt_enerbil('LW',itap)
4347
4348	!
4349	IF (mydebug) THEN
4350	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4351	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4352	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4353	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4354	ENDIF
4355
4356	! Calculer l'hydrologie de la surface
4357	!
4358	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4359	! . agesno, ftsol,fqsurf,fsnow, ruis)
4360	!
4361
4362	!
4363	! Calculer le bilan du sol et la derive de temperature (couplage)
4364	!
4365	DO i = 1, klon
4366	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4367	! a la demande de JLD
4368	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4369	ENDDO
4370	!
4371	!moddeblott(jan95)
4372	! Appeler le programme de parametrisation de l'orographie
4373	! a l'echelle sous-maille:
4374	!
4375	IF (prt_level .GE.10) THEN
4376	print *,' call orography ? ', ok_orodr
4377	ENDIF
4378	!
4379	IF (ok_orodr) THEN
4380	!
4381	! selection des points pour lesquels le shema est actif:
4382	igwd=0
4383	DO i=1,klon
4384	itest(i)=0
4385	! IF ((zstd(i).gt.10.0)) THEN
4386	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4387	itest(i)=1
4388	igwd=igwd+1
4389	idx(igwd)=i
4390	ENDIF
4391	ENDDO
4392	! igwdim=MAX(1,igwd)
4393	!
4394	IF (ok_strato) THEN
4395
4396	CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
4397	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4398	igwd,idx,itest, &
4399	t_seri, u_seri, v_seri, &
4400	zulow, zvlow, zustrdr, zvstrdr, &
4401	d_t_oro, d_u_oro, d_v_oro)
4402
4403	ELSE
4404	CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
4405	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4406	igwd,idx,itest, &
4407	t_seri, u_seri, v_seri, &
4408	zulow, zvlow, zustrdr, zvstrdr, &
4409	d_t_oro, d_u_oro, d_v_oro)
4410	ENDIF
4411	!
4412	! ajout des tendances
4413	!-----------------------------------------------------------------------
4414	! ajout des tendances de la trainee de l'orographie
4415	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
4416	abortphy,flag_inhib_tend,itap,0)
4417	CALL prt_enerbil('oro',itap)
4418	!----------------------------------------------------------------------
4419	!
4420	ENDIF ! fin de test sur ok_orodr
4421	!
4422	IF (mydebug) THEN
4423	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4424	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4425	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4426	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4427	ENDIF
4428
4429	IF (ok_orolf) THEN
4430	!
4431	! selection des points pour lesquels le shema est actif:
4432	igwd=0
4433	DO i=1,klon
4434	itest(i)=0
4435	IF ((zpic(i)-zmea(i)).GT.100.) THEN
4436	itest(i)=1
4437	igwd=igwd+1
4438	idx(igwd)=i
4439	ENDIF
4440	ENDDO
4441	! igwdim=MAX(1,igwd)
4442	!
4443	IF (ok_strato) THEN
4444
4445	CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
4446	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4447	igwd,idx,itest, &
4448	t_seri, u_seri, v_seri, &
4449	zulow, zvlow, zustrli, zvstrli, &
4450	d_t_lif, d_u_lif, d_v_lif )
4451
4452	ELSE
4453	CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
4454	latitude_deg,zmea,zstd,zpic, &
4455	itest, &
4456	t_seri, u_seri, v_seri, &
4457	zulow, zvlow, zustrli, zvstrli, &
4458	d_t_lif, d_u_lif, d_v_lif)
4459	ENDIF
4460
4461	! ajout des tendances de la portance de l'orographie
4462	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
4463	'lif', abortphy,flag_inhib_tend,itap,0)
4464	CALL prt_enerbil('lif',itap)
4465	ENDIF ! fin de test sur ok_orolf
4466
4467	IF (ok_hines) then
4468	! HINES GWD PARAMETRIZATION
4469	east_gwstress=0.
4470	west_gwstress=0.
4471	du_gwd_hines=0.
4472	dv_gwd_hines=0.
4473	CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
4474	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4475	du_gwd_hines, dv_gwd_hines)
4476	zustr_gwd_hines=0.
4477	zvstr_gwd_hines=0.
4478	DO k = 1, klev
4479	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
4480	* (paprs(:, k)-paprs(:, k+1))/rg
4481	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
4482	* (paprs(:, k)-paprs(:, k+1))/rg
4483	ENDDO
4484
4485	d_t_hin(:, :)=0.
4486	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4487	dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4488	CALL prt_enerbil('hin',itap)
4489	ENDIF
4490
4491	IF (.not. ok_hines .and. ok_gwd_rando) then
4492	! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
4493	CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
4494	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4495	dv_gwd_front, east_gwstress, west_gwstress)
4496	zustr_gwd_front=0.
4497	zvstr_gwd_front=0.
4498	DO k = 1, klev
4499	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
4500	* (paprs(:, k)-paprs(:, k+1))/rg
4501	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
4502	* (paprs(:, k)-paprs(:, k+1))/rg
4503	ENDDO
4504
4505	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
4506	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4507	CALL prt_enerbil('front_gwd_rando',itap)
4508	ENDIF
4509
4510	IF (ok_gwd_rando) THEN
4511	CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
4512	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4513	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4514	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
4515	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4516	CALL prt_enerbil('flott_gwd_rando',itap)
4517	zustr_gwd_rando=0.
4518	zvstr_gwd_rando=0.
4519	DO k = 1, klev
4520	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
4521	* (paprs(:, k)-paprs(:, k+1))/rg
4522	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
4523	* (paprs(:, k)-paprs(:, k+1))/rg
4524	ENDDO
4525	ENDIF
4526
4527	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4528
4529	IF (mydebug) THEN
4530	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4531	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4532	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4533	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4534	ENDIF
4535
4536	DO i = 1, klon
4537	zustrph(i)=0.
4538	zvstrph(i)=0.
4539	ENDDO
4540	DO k = 1, klev
4541	DO i = 1, klon
4542	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
4543	(paprs(i,k)-paprs(i,k+1))/rg
4544	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
4545	(paprs(i,k)-paprs(i,k+1))/rg
4546	ENDDO
4547	ENDDO
4548	!
4549	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4550	!
4551	IF (is_sequential .and. ok_orodr) THEN
4552	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4553	ra,rg,romega, &
4554	latitude_deg,longitude_deg,pphis, &
4555	zustrdr,zustrli,zustrph, &
4556	zvstrdr,zvstrli,zvstrph, &
4557	paprs,u,v, &
4558	aam, torsfc)
4559	ENDIF
4560	!IM cf. FLott END
4561	!DC Calcul de la tendance due au methane
4562	IF (ok_qch4) THEN
4563	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4564	! ajout de la tendance d'humidite due au methane
4565	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
4566	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
4567	'q_ch4', abortphy,flag_inhib_tend,itap,0)
4568	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
4569	ENDIF
4570	!
4571	!
4572
4573	!===============================================================
4574	! Additional tendency of TKE due to orography
4575	!===============================================================
4576	!
4577	! Inititialization
4578	!------------------
4579
4580	addtkeoro=0
4581	CALL getin_p('addtkeoro',addtkeoro)
4582
4583	IF (prt_level.ge.5) &
4584	print*,'addtkeoro', addtkeoro
4585
4586	alphatkeoro=1.
4587	CALL getin_p('alphatkeoro',alphatkeoro)
4588	alphatkeoro=min(max(0.,alphatkeoro),1.)
4589
4590	smallscales_tkeoro=.FALSE.
4591	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4592
4593
4594	dtadd(:,:)=0.
4595	duadd(:,:)=0.
4596	dvadd(:,:)=0.
4597
4598	! Choices for addtkeoro:
4599	! ** 0 no TKE tendency from orography
4600	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4601	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4602	!
4603
4604	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4605	! -------------------------------------------
4606
4607
4608	! selection des points pour lesquels le schema est actif:
4609
4610
4611	IF (addtkeoro .EQ. 1 ) THEN
4612
4613	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4614	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4615
4616	ELSE IF (addtkeoro .EQ. 2) THEN
4617
4618	IF (smallscales_tkeoro) THEN
4619	igwd=0
4620	DO i=1,klon
4621	itest(i)=0
4622	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4623	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4624	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4625	IF (zstd(i).GT.1.0) THEN
4626	itest(i)=1
4627	igwd=igwd+1
4628	idx(igwd)=i
4629	ENDIF
4630	ENDDO
4631
4632	ELSE
4633
4634	igwd=0
4635	DO i=1,klon
4636	itest(i)=0
4637	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4638	itest(i)=1
4639	igwd=igwd+1
4640	idx(igwd)=i
4641	ENDIF
4642	ENDDO
4643
4644	ENDIF
4645
4646	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
4647	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4648	igwd,idx,itest, &
4649	t_seri, u_seri, v_seri, &
4650	zulow, zvlow, zustrdr, zvstrdr, &
4651	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4652
4653	zustrdr(:)=0.
4654	zvstrdr(:)=0.
4655	zulow(:)=0.
4656	zvlow(:)=0.
4657
4658	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4659	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4660	ENDIF
4661
4662
4663	! TKE update from subgrid temperature and wind tendencies
4664	!----------------------------------------------------------
4665	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4666
4667
4668	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
4669
4670
4671	ENDIF
4672	! -----
4673	!===============================================================
4674
4675
4676	!====================================================================
4677	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4678	!====================================================================
4679	! Abderrahmane 24.08.09
4680
4681	IF (ok_cosp) THEN
4682	! adeclarer
4683	#ifdef CPP_COSP
4684	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4685
4686	IF (prt_level .GE.10) THEN
4687	print*,'freq_cosp',freq_cosp
4688	ENDIF
4689	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4690	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4691	! s ref_liq,ref_ice
4692	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
4693	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4694	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4695	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4696	JrNt,ref_liq,ref_ice, &
4697	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4698	zu10m,zv10m,pphis, &
4699	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4700	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4701	prfl(:,1:klev),psfl(:,1:klev), &
4702	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4703	mr_ozone,cldtau, cldemi)
4704
4705	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4706	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4707	! M clMISR,
4708	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4709	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4710
4711	ENDIF
4712	#endif
4713
4714	#ifdef CPP_COSP2
4715	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4716
4717	IF (prt_level .GE.10) THEN
4718	print*,'freq_cosp',freq_cosp
4719	ENDIF
4720	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4721	print*,'Dans physiq.F avant appel '
4722	! s ref_liq,ref_ice
4723	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
4724	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4725	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4726	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4727	JrNt,ref_liq,ref_ice, &
4728	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4729	zu10m,zv10m,pphis, &
4730	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4731	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4732	prfl(:,1:klev),psfl(:,1:klev), &
4733	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4734	mr_ozone,cldtau, cldemi)
4735	ENDIF
4736	#endif
4737
4738	#ifdef CPP_COSPV2
4739	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4740	! IF (MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4741
4742	IF (prt_level .GE.10) THEN
4743	print*,'freq_cosp',freq_cosp
4744	ENDIF
4745	DO k = 1, klev
4746	DO i = 1, klon
4747	phicosp(i,k) = pphi(i,k) + pphis(i)
4748	ENDDO
4749	ENDDO
4750	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4751	print*,'Dans physiq.F avant appel '
4752	! s ref_liq,ref_ice
4753	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
4754	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4755	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4756	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4757	JrNt,ref_liq,ref_ice, &
4758	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4759	zu10m,zv10m,pphis, &
4760	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4761	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4762	prfl(:,1:klev),psfl(:,1:klev), &
4763	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4764	mr_ozone,cldtau, cldemi)
4765	ENDIF
4766	#endif
4767
4768	ENDIF !ok_cosp
4769
4770
4771	! Marine
4772
4773	IF (ok_airs) then
4774
4775	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
4776	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4777	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4778	& map_prop_hc,map_prop_hist,&
4779	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4780	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4781	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4782	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4783	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4784	& map_ntot,map_hc,map_hist,&
4785	& map_Cb,map_ThCi,map_Anv,&
4786	& alt_tropo )
4787	ENDIF
4788
4789	ENDIF ! ok_airs
4790
4791
4792	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4793	!AA
4794	!AA Installation de l'interface online-offline pour traceurs
4795	!AA
4796	!====================================================================
4797	! Calcul des tendances traceurs
4798	!====================================================================
4799	!
4800
4801	IF (type_trac=='repr') THEN
4802	!MM pas d'impact, car on recupere q_seri,tr_seri,t_seri via phys_local_var_mod
4803	!MM dans Reprobus
4804	sh_in(:,:) = q_seri(:,:)
4805	#ifdef REPROBUS
4806	d_q_rep(:,:) = 0.
4807	d_ql_rep(:,:) = 0.
4808	d_qi_rep(:,:) = 0.
4809	#endif
4810	ELSE
4811	sh_in(:,:) = qx(:,:,ivap)
4812	ch_in(:,:) = qx(:,:,iliq)
4813	ENDIF
4814
4815	#ifdef CPP_Dust
4816	! Avec SPLA, iflag_phytrac est forcé =1
4817	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
4818	pdtphys,ftsol, & ! I
4819	t,q_seri,paprs,pplay,RHcl, & ! I
4820	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
4821	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
4822	u_seri, v_seri, latitude_deg, longitude_deg, &
4823	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
4824	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
4825	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
4826	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
4827	fm_therm, entr_therm, rneb, & ! I
4828	beta_prec_fisrt,beta_prec, & !I
4829	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
4830	d_tr_dyn,tr_seri)
4831
4832	#else
4833	IF (iflag_phytrac == 1 ) THEN
4834	CALL phytrac ( &
4835	itap, days_elapsed+1, jH_cur, debut, &
4836	lafin, phys_tstep, u, v, t, &
4837	paprs, pplay, pmfu, pmfd, &
4838	pen_u, pde_u, pen_d, pde_d, &
4839	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
4840	u1, v1, ftsol, pctsrf, &
4841	zustar, zu10m, zv10m, &
4842	wstar(:,is_ave), ale_bl, ale_wake, &
4843	latitude_deg, longitude_deg, &
4844	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4845	presnivs, pphis, pphi, albsol1, &
4846	sh_in, ch_in, rhcl, cldfra, rneb, &
4847	diafra, cldliq, itop_con, ibas_con, &
4848	pmflxr, pmflxs, prfl, psfl, &
4849	da, phi, mp, upwd, &
4850	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
4851	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
4852	ev, ep, epmlmMm, eplaMm, & !<<RomP
4853	dnwd, aerosol_couple, flxmass_w, &
4854	tau_aero, piz_aero, cg_aero, ccm, &
4855	rfname, &
4856	d_tr_dyn, & !<<RomP
4857	tr_seri, init_source)
4858	#ifdef REPROBUS
4859
4860
4861	print*,'avt add phys rep',abortphy
4862
4863	CALL add_phys_tend &
4864	(du0,dv0,dt0,d_q_rep,d_ql_rep,d_qi_rep,paprs,&
4865	'rep',abortphy,flag_inhib_tend,itap,0)
4866	IF (abortphy==1) Print*,'ERROR ABORT REP'
4867
4868	print*,'apr add phys rep',abortphy
4869
4870	#endif
4871	ENDIF ! (iflag_phytrac=1)
4872
4873	#endif
4874	!ENDIF ! (iflag_phytrac=1)
4875
4876	IF (offline) THEN
4877
4878	IF (prt_level.ge.9) &
4879	print*,'Attention on met a 0 les thermiques pour phystoke'
4880	CALL phystokenc ( &
4881	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4882	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4883	fm_therm,entr_therm, &
4884	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4885	frac_impa, frac_nucl, &
4886	pphis,cell_area,phys_tstep,itap, &
4887	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4888
4889
4890	ENDIF
4891
4892	!
4893	! Calculer le transport de l'eau et de l'energie (diagnostique)
4894	!
4895	CALL transp (paprs,zxtsol, &
4896	t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
4897	ve, vq, ue, uq, vwat, uwat)
4898	!
4899	!IM global posePB BEG
4900	IF(1.EQ.0) THEN
4901	!
4902	CALL transp_lay (paprs,zxtsol, &
4903	t_seri, q_seri, u_seri, v_seri, zphi, &
4904	ve_lay, vq_lay, ue_lay, uq_lay)
4905	!
4906	ENDIF !(1.EQ.0) THEN
4907	!IM global posePB END
4908	! Accumuler les variables a stocker dans les fichiers histoire:
4909	!
4910
4911	!================================================================
4912	! Conversion of kinetic and potential energy into heat, for
4913	! parameterisation of subgrid-scale motions
4914	!================================================================
4915
4916	d_t_ec(:,:)=0.
4917	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4918	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
4919	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4920	zmasse,exner,d_t_ec)
4921	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4922
4923	!=======================================================================
4924	! SORTIES
4925	!=======================================================================
4926	!
4927	!IM initialisation + calculs divers diag AMIP2
4928	!
4929	include "calcul_divers.h"
4930	!
4931	!IM Interpolation sur les niveaux de pression du NMC
4932	! -------------------------------------------------
4933	!
4934	include "calcul_STDlev.h"
4935	!
4936	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
4937	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
4938	!
4939	!cc prw = eau precipitable
4940	! prlw = colonne eau liquide
4941	! prlw = colonne eau solide
4942	prw(:) = 0.
4943	prlw(:) = 0.
4944	prsw(:) = 0.
4945	DO k = 1, klev
4946	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
4947	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
4948	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
4949	ENDDO
4950	!
4951	IF (type_trac == 'inca') THEN
4952	#ifdef INCA
4953	CALL VTe(VTphysiq)
4954	CALL VTb(VTinca)
4955
4956	CALL chemhook_end ( &
4957	phys_tstep, &
4958	pplay, &
4959	t_seri, &
4960	tr_seri, &
4961	nbtr, &
4962	paprs, &
4963	q_seri, &
4964	cell_area, &
4965	pphi, &
4966	pphis, &
4967	zx_rh, &
4968	aps, bps, ap, bp)
4969
4970	CALL VTe(VTinca)
4971	CALL VTb(VTphysiq)
4972	#endif
4973	ENDIF
4974
4975
4976	!
4977	! Convertir les incrementations en tendances
4978	!
4979	IF (prt_level .GE.10) THEN
4980	print *,'Convertir les incrementations en tendances '
4981	ENDIF
4982	!
4983	IF (mydebug) THEN
4984	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4985	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4986	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4987	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4988	ENDIF
4989
4990	DO k = 1, klev
4991	DO i = 1, klon
4992	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
4993	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
4994	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
4995	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
4996	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
4997	!CR: on ajoute le contenu en glace
4998	IF (nqo.eq.3) THEN
4999	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
5000	ENDIF
5001	ENDDO
5002	ENDDO
5003	!
5004	!CR: nb de traceurs eau: nqo
5005	! IF (nqtot.GE.3) THEN
5006	IF (nqtot.GE.(nqo+1)) THEN
5007	! DO iq = 3, nqtot
5008	DO iq = nqo+1, nqtot
5009	DO k = 1, klev
5010	DO i = 1, klon
5011	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / phys_tstep
5012	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / phys_tstep
5013	ENDDO
5014	ENDDO
5015	ENDDO
5016	ENDIF
5017	!
5018	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
5019	!IM global posePB include "write_bilKP_ins.h"
5020	!IM global posePB include "write_bilKP_ave.h"
5021	!
5022
5023	!--OB mass fixer
5024	!--profile is corrected to force mass conservation of water
5025	IF (mass_fixer) THEN
5026	qql2(:)=0.0
5027	DO k = 1, klev
5028	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
5029	ENDDO
5030	DO i = 1, klon
5031	!--compute ratio of what q+ql should be with conservation to what it is
5032	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
5033	DO k = 1, klev
5034	q_seri(i,k) =q_seri(i,k)*corrqql
5035	ql_seri(i,k)=ql_seri(i,k)*corrqql
5036	ENDDO
5037	ENDDO
5038	ENDIF
5039	!--fin mass fixer
5040
5041	! Sauvegarder les valeurs de t et q a la fin de la physique:
5042	!
5043	u_ancien(:,:) = u_seri(:,:)
5044	v_ancien(:,:) = v_seri(:,:)
5045	t_ancien(:,:) = t_seri(:,:)
5046	q_ancien(:,:) = q_seri(:,:)
5047	ql_ancien(:,:) = ql_seri(:,:)
5048	qs_ancien(:,:) = qs_seri(:,:)
5049	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
5050	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
5051	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
5052	! !! RomP >>>
5053	!CR: nb de traceurs eau: nqo
5054	IF (nqtot.GT.nqo) THEN
5055	DO iq = nqo+1, nqtot
5056	tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
5057	ENDDO
5058	ENDIF
5059	! !! RomP <<<
5060	!==========================================================================
5061	! Sorties des tendances pour un point particulier
5062	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
5063	! pour le debug
5064	! La valeur de igout est attribuee plus haut dans le programme
5065	!==========================================================================
5066
5067	IF (prt_level.ge.1) THEN
5068	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
5069	write(lunout,*) &
5070	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
5071	write(lunout,*) &
5072	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
5073	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
5074	pctsrf(igout,is_sic)
5075	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
5076	DO k=1,klev
5077	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
5078	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
5079	d_t_eva(igout,k)
5080	ENDDO
5081	write(lunout,*) 'cool,heat'
5082	DO k=1,klev
5083	write(lunout,*) cool(igout,k),heat(igout,k)
5084	ENDDO
5085
5086	!jyg< (En attendant de statuer sur le sort de d_t_oli)
5087	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5088	!jyg! do k=1,klev
5089	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
5090	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5091	!jyg! enddo
5092	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
5093	DO k=1,klev
5094	write(lunout,*) d_t_vdf(igout,k), &
5095	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
5096	ENDDO
5097	!>jyg
5098
5099	write(lunout,*) 'd_ps ',d_ps(igout)
5100	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
5101	DO k=1,klev
5102	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
5103	d_qx(igout,k,1),d_qx(igout,k,2)
5104	ENDDO
5105	ENDIF
5106
5107	!============================================================
5108	! Calcul de la temperature potentielle
5109	!============================================================
5110	DO k = 1, klev
5111	DO i = 1, klon
5112	!JYG/IM theta en debut du pas de temps
5113	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
5114	!JYG/IM theta en fin de pas de temps de physique
5115	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
5116	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
5117	! MPL 20130625
5118	! fth_fonctions.F90 et parkind1.F90
5119	! sinon thetal=theta
5120	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
5121	! : ql_seri(i,k))
5122	thetal(i,k)=theta(i,k)
5123	ENDDO
5124	ENDDO
5125	!
5126
5127	! 22.03.04 BEG
5128	!=============================================================
5129	! Ecriture des sorties
5130	!=============================================================
5131	#ifdef CPP_IOIPSL
5132
5133	! Recupere des varibles calcule dans differents modules
5134	! pour ecriture dans histxxx.nc
5135
5136	! Get some variables from module fonte_neige_mod
5137	CALL fonte_neige_get_vars(pctsrf, &
5138	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
5139
5140
5141	!=============================================================
5142	! Separation entre thermiques et non thermiques dans les sorties
5143	! de fisrtilp
5144	!=============================================================
5145
5146	IF (iflag_thermals>=1) THEN
5147	d_t_lscth=0.
5148	d_t_lscst=0.
5149	d_q_lscth=0.
5150	d_q_lscst=0.
5151	DO k=1,klev
5152	DO i=1,klon
5153	IF (ptconvth(i,k)) THEN
5154	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5155	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5156	ELSE
5157	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5158	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5159	ENDIF
5160	ENDDO
5161	ENDDO
5162
5163	DO i=1,klon
5164	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
5165	plul_th(i)=prfl(i,1)+psfl(i,1)
5166	ENDDO
5167	ENDIF
5168
5169	!On effectue les sorties:
5170
5171	#ifdef CPP_Dust
5172	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
5173	pplay, lmax_th, aerosol_couple, &
5174	ok_ade, ok_aie, ivap, ok_sync, &
5175	ptconv, read_climoz, clevSTD, &
5176	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
5177	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5178	#else
5179	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
5180	pplay, lmax_th, aerosol_couple, &
5181	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, &
5182	ok_sync, ptconv, read_climoz, clevSTD, &
5183	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
5184	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5185	#endif
5186
5187	#ifndef CPP_XIOS
5188	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
5189	#endif
5190
5191	#endif
5192
5193	! Pour XIOS : On remet des variables a .false. apres un premier appel
5194	IF (debut) THEN
5195	#ifdef CPP_XIOS
5196	swaero_diag=.FALSE.
5197	swaerofree_diag=.FALSE.
5198	dryaod_diag=.FALSE.
5199	ok_4xCO2atm= .FALSE.
5200	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
5201
5202	IF (is_master) THEN
5203	!--setting up swaero_diag to TRUE in XIOS case
5204	IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
5205	xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
5206	xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR. &
5207	(iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
5208	xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0")))) &
5209	!!!--for now these fields are not in the XML files so they are omitted
5210	!!! xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
5211	swaero_diag=.TRUE.
5212
5213	!--setting up swaerofree_diag to TRUE in XIOS case
5214	IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
5215	xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR. &
5216	xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
5217	xios_field_is_active("LWupTOAcleanclr")) &
5218	swaerofree_diag=.TRUE.
5219
5220	!--setting up dryaod_diag to TRUE in XIOS case
5221	DO naero = 1, naero_tot-1
5222	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
5223	ENDDO
5224	!
5225	!--setting up ok_4xCO2atm to TRUE in XIOS case
5226	IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
5227	xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
5228	xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
5229	xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
5230	xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
5231	xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
5232	ok_4xCO2atm=.TRUE.
5233	ENDIF
5234	!$OMP BARRIER
5235	CALL bcast(swaero_diag)
5236	CALL bcast(swaerofree_diag)
5237	CALL bcast(dryaod_diag)
5238	CALL bcast(ok_4xCO2atm)
5239	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
5240	#endif
5241	ENDIF
5242
5243	!====================================================================
5244	! Arret du modele apres hgardfou en cas de detection d'un
5245	! plantage par hgardfou
5246	!====================================================================
5247
5248	IF (abortphy==1) THEN
5249	abort_message ='Plantage hgardfou'
5250	CALL abort_physic (modname,abort_message,1)
5251	ENDIF
5252
5253	! 22.03.04 END
5254	!
5255	!====================================================================
5256	! Si c'est la fin, il faut conserver l'etat de redemarrage
5257	!====================================================================
5258	!
5259
5260	IF (lafin) THEN
5261	itau_phy = itau_phy + itap
5262	CALL phyredem ("restartphy.nc")
5263	! open(97,form="unformatted",file="finbin")
5264	! write(97) u_seri,v_seri,t_seri,q_seri
5265	! close(97)
5266
5267	IF (is_omp_master) THEN
5268
5269	IF (read_climoz >= 1) THEN
5270	IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5271	DEALLOCATE(press_edg_climoz) ! pointer
5272	DEALLOCATE(press_cen_climoz) ! pointer
5273	ENDIF
5274
5275	ENDIF
5276	#ifdef CPP_XIOS
5277	IF (is_omp_master) CALL xios_context_finalize
5278	#endif
5279	WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
5280	ENDIF
5281
5282	! first=.false.
5283
5284	END SUBROUTINE physiq
5285
5286	END MODULE physiq_mod

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

Download in other formats:

Original Format