Context Navigation

physiq_mod.F90 @ 3849

Last change on this file since 3849 was 3849, checked in by acozic, 3 years ago
add a test for calling finalize_inca in physiq_mod only in case of unstructured grid
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: 187.1 KB

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

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

Download in other formats:

Original Format