Context Navigation

physiq_mod.F90 @ 3106

Last change on this file since 3106 was 3106, checked in by oboucher, 6 years ago
Adding clear-sky clean (no aerosol) LW radiative fluxes at TOA and SFC as diagnostic
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: 165.0 KB

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

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

Download in other formats:

Original Format