Context Navigation

physiq_mod.F90 @ 3048

Last change on this file since 3048 was 3048, checked in by musat, 7 years ago
Set up lwoff configuration : add ok_lwoff=y to activate it. Add solbnd CMIP6 output. Enable 3D outputs of SW/LW total sky and clear sky. IM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 164.2 KB

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

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

Download in other formats:

Original Format