Context Navigation

physiq_mod.F90 @ 3999

Last change on this file since 3999 was 3999, checked in by evignon, 3 years ago

commission de la nouvelle routine de condensation
grande echelle simplifiee (lscp, version epuree de fistilp)
et du schema de nuages de phase mixte (en developpement)

La routine lscp n'est active que sous flag
ok_new_lscp=y

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Id

File size: 186.5 KB

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

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

Download in other formats:

Original Format