Context Navigation

physiq_mod.F90 @ 2974

Last change on this file since 2974 was 2974, checked in by Laurent Fairhead, 7 years ago
Inclusion of iflag_phytrac in 6.0.12 branch 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.2 KB

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

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

Download in other formats:

Original Format