Context Navigation

physiq_mod.F90 @ 3776

Last change on this file since 3776 was 3776, checked in by asima, 4 years ago
Preparing physiq_mod.def for LMDZ-SPLA update. Also, section "USE" re-arranged in alphabetical order.
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 `Author Date Id Revision`
File size: 186.5 KB

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

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

Download in other formats:

Original Format