Context Navigation

physiq_mod.F90 @ 2973

Last change on this file since 2973 was 2973, checked in by Laurent Fairhead, 7 years ago
Bypassing the call to phytrac when not needed. Controlled by iflag_phytrac read in physiq.def and which is 0 by default (thus no call to phytrac by default) LF
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: 161.1 KB

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

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

Download in other formats:

Original Format