Context Navigation

physiq_mod.F90 @ 2984

Last change on this file since 2984 was 2984, checked in by Laurent Fairhead, 7 years ago
We actually do want to call phytrac by default LF
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 161.2 KB

Line
1	!
2	! $Id: physiq_mod.F90 2984 2017-09-08 08:39:25Z fairhead $
3	!
4	!#define IO_DEBUG
5	MODULE physiq_mod
6
7	IMPLICIT NONE
8
9	CONTAINS
10
11	SUBROUTINE physiq (nlon,nlev, &
12	debut,lafin,pdtphys_, &
13	paprs,pplay,pphi,pphis,presnivs, &
14	u,v,rot,t,qx, &
15	flxmass_w, &
16	d_u, d_v, d_t, d_qx, d_ps)
17
18	use assert_m, only: assert
19	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
20	histwrite, ju2ymds, ymds2ju, getin
21	USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
22	USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
23	year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour
24	USE write_field_phy
25	USE dimphy
26	USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac
27	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo
28	USE mod_phys_lmdz_para
29	USE iophy
30	USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
31	USE phystokenc_mod, ONLY: offline, phystokenc
32	USE time_phylmdz_mod, only: raz_date, day_step_phy, update_time,current_time
33	USE vampir
34	USE pbl_surface_mod, ONLY : pbl_surface
35	USE change_srf_frac_mod
36	USE surface_data, ONLY : type_ocean, ok_veget, ok_snow
37	USE tropopause_m, ONLY: dyn_tropopause
38	#ifdef CPP_Dust
39	USE phytracr_spl_mod, ONLY: phytracr_spl
40	#endif
41	USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
42	! [Variables internes non sauvegardees de la physique]
43	! Variables locales pour effectuer les appels en serie
44	t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri, &
45	! Dynamic tendencies (diagnostics)
46	d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn, &
47	d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
48	! Physic tendencies
49	d_t_con,d_q_con,d_u_con,d_v_con, &
50	d_tr, & !! to be removed?? (jyg)
51	d_t_wake,d_q_wake, &
52	d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
53	d_t_ajsb,d_q_ajsb, &
54	d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
55	d_t_ajs_w,d_q_ajs_w, &
56	d_t_ajs_x,d_q_ajs_x, &
57	!
58	d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
59	d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
60	d_t_lscst,d_q_lscst, &
61	d_t_lscth,d_q_lscth, &
62	plul_st,plul_th, &
63	!
64	d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
65	d_t_vdf_w,d_q_vdf_w, &
66	d_t_vdf_x,d_q_vdf_x, &
67	d_ts, &
68	!
69	d_t_oli,d_u_oli,d_v_oli, &
70	d_t_oro,d_u_oro,d_v_oro, &
71	d_t_oro_gw,d_u_oro_gw,d_v_oro_gw, &
72	d_t_lif,d_u_lif,d_v_lif, &
73	d_t_ec, &
74	!
75	du_gwd_hines,dv_gwd_hines,d_t_hin, &
76	dv_gwd_rando,dv_gwd_front, &
77	east_gwstress,west_gwstress, &
78	d_q_ch4, &
79	! Special RRTM
80	ZLWFT0_i,ZSWFT0_i,ZFLDN0, &
81	ZFLUP0,ZFSDN0,ZFSUP0, &
82	!
83	topswad_aero,solswad_aero, &
84	topswai_aero,solswai_aero, &
85	topswad0_aero,solswad0_aero, &
86	!LW additional
87	toplwad_aero,sollwad_aero, &
88	toplwai_aero,sollwai_aero, &
89	toplwad0_aero,sollwad0_aero, &
90	!
91	topsw_aero,solsw_aero, &
92	topsw0_aero,solsw0_aero, &
93	topswcf_aero,solswcf_aero, &
94	tausum_aero,tau3d_aero, &
95	drytausum_aero, &
96	!
97	!variables CFMIP2/CMIP5
98	topswad_aerop, solswad_aerop, &
99	topswai_aerop, solswai_aerop, &
100	topswad0_aerop, solswad0_aerop, &
101	topsw_aerop, topsw0_aerop, &
102	solsw_aerop, solsw0_aerop, &
103	topswcf_aerop, solswcf_aerop, &
104	!LW diagnostics
105	toplwad_aerop, sollwad_aerop, &
106	toplwai_aerop, sollwai_aerop, &
107	toplwad0_aerop, sollwad0_aerop, &
108	!
109	ptstar, pt0, slp, &
110	!
111	bils, &
112	!
113	cldh, cldl,cldm, cldq, cldt, &
114	JrNt, &
115	dthmin, evap, fder, plcl, plfc, &
116	prw, prlw, prsw, &
117	s_lcl, s_pblh, s_pblt, s_therm, &
118	cdragm, cdragh, &
119	zustar, zu10m, zv10m, rh2m, qsat2m, &
120	zq2m, zt2m, weak_inversion, &
121	zt2m_min_mon, zt2m_max_mon, & ! pour calcul_divers.h
122	t2m_min_mon, t2m_max_mon, & ! pour calcul_divers.h
123	!
124	s_pblh_x, s_pblh_w, &
125	s_lcl_x, s_lcl_w, &
126	!
127	slab_wfbils, tpot, tpote, &
128	ue, uq, ve, vq, zxffonte, &
129	zxfqcalving, zxfluxlat, &
130	zxrunofflic, &
131	zxtsol, snow_lsc, zxfqfonte, zxqsurf, &
132	rain_lsc, rain_num, &
133	!
134	sens_x, sens_w, &
135	zxfluxlat_x, zxfluxlat_w, &
136	!
137	dtvdf_x, dtvdf_w, &
138	dqvdf_x, dqvdf_w, &
139	pbl_tke_input, &
140	t_therm, q_therm, u_therm, v_therm, &
141	cdragh_x, cdragh_w, &
142	cdragm_x, cdragm_w, &
143	kh, kh_x, kh_w, &
144	!
145	wake_k, &
146	ale_wake, alp_wake, &
147	wake_h, wake_omg, &
148	! tendencies of delta T and delta q:
149	d_deltat_wk, d_deltaq_wk, & ! due to wakes
150	d_deltat_wk_gw, d_deltaq_wk_gw, & ! due to wake induced gravity waves
151	d_deltat_vdf, d_deltaq_vdf, & ! due to vertical diffusion
152	d_deltat_the, d_deltaq_the, & ! due to thermals
153	d_deltat_ajs_cv, d_deltaq_ajs_cv, & ! due to dry adjustment of (w) before convection
154	! tendencies of wake fractional area and wake number per unit area:
155	d_s_wk, d_dens_wk, & ! due to wakes
156	!!! d_s_vdf, d_dens_vdf, & ! due to vertical diffusion
157	!!! d_s_the, d_dens_the, & ! due to thermals
158	!
159	ptconv, &
160	wbeff, convoccur, zmax_th, &
161	sens, flwp, fiwp, &
162	ale_bl_stat,alp_bl_conv,alp_bl_det, &
163	alp_bl_fluct_m,alp_bl_fluct_tke, &
164	alp_bl_stat, n2, s2, &
165	proba_notrig, random_notrig, &
166	!
167	dnwd, dnwd0, &
168	upwd, omega, &
169	epmax_diag, &
170	ep, &
171	cldemi, &
172	cldfra, cldtau, fiwc, &
173	fl, re, flwc, &
174	ref_liq, ref_ice, theta, &
175	ref_liq_pi, ref_ice_pi, &
176	zphi, zx_rh, &
177	pmfd, pmfu, &
178	!
179	t2m, fluxlat, &
180	fsollw, evap_pot, &
181	fsolsw, wfbils, wfbilo, &
182	wfevap, wfrain, wfsnow, &
183	pmflxr, pmflxs, prfl, &
184	psfl, fraca, Vprecip, &
185	zw2, &
186
187	fluxu, fluxv, &
188	fluxt, &
189
190	uwriteSTD, vwriteSTD, & !pour calcul_STDlev.h
191	wwriteSTD, phiwriteSTD, & !pour calcul_STDlev.h
192	qwriteSTD, twriteSTD, rhwriteSTD, & !pour calcul_STDlev.h
193
194	wdtrainA, wdtrainM, &
195	beta_prec, &
196	rneb, &
197	zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic
198	!
199	USE phys_state_var_mod ! Variables sauvegardees de la physique
200	#ifdef CPP_Dust
201	USE phys_output_write_spl_mod
202	#else
203	USE phys_output_var_mod ! Variables pour les ecritures des sorties
204	#endif
205
206	USE phys_output_write_mod
207	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
208	USE phys_output_mod
209	USE phys_output_ctrlout_mod
210	use open_climoz_m, only: open_climoz ! ozone climatology from a file
211	use regr_pr_time_av_m, only: regr_pr_time_av
212	use netcdf95, only: nf95_close
213	!IM for NMC files
214	! use netcdf, only: nf90_fill_real
215	use netcdf
216	use mod_phys_lmdz_mpi_data, only: is_mpi_root
217	USE aero_mod
218	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
219	use conf_phys_m, only: conf_phys
220	use radlwsw_m, only: radlwsw
221	use phyaqua_mod, only: zenang_an
222	USE time_phylmdz_mod, only: day_step_phy, annee_ref, day_ref, itau_phy, &
223	start_time, pdtphys, day_ini
224	USE tracinca_mod, ONLY: config_inca
225	#ifdef CPP_XIOS
226	USE wxios, ONLY: missing_val, missing_val_omp
227	USE xios, ONLY: xios_get_field_attr, xios_field_is_active
228	#endif
229	#ifdef REPROBUS
230	USE CHEM_REP, ONLY : Init_chem_rep_xjour
231	#endif
232	USE indice_sol_mod
233	USE phytrac_mod, ONLY : phytrac
234
235	#ifdef CPP_RRTM
236	USE YOERAD, ONLY : NRADLP
237	USE YOESW, ONLY : RSUN
238	#endif
239	USE ioipsl_getin_p_mod, ONLY : getin_p
240
241	#ifndef CPP_XIOS
242	USE paramLMDZ_phy_mod
243	#endif
244
245	USE cmp_seri_mod
246	USE add_phys_tend_mod, only : add_pbl_tend, add_phys_tend, diag_phys_tend, prt_enerbil, &
247	& fl_ebil, fl_cor_ebil
248
249	!IM stations CFMIP
250	USE CFMIP_point_locations
251	use FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
252	use ACAMA_GWD_rando_m, only: ACAMA_GWD_rando
253	USE VERTICAL_LAYERS_MOD, ONLY: aps,bps
254
255
256	IMPLICIT none
257	!>======================================================================
258	!!
259	!! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
260	!!
261	!! Objet: Moniteur general de la physique du modele
262	!!AA Modifications quant aux traceurs :
263	!!AA - uniformisation des parametrisations ds phytrac
264	!!AA - stockage des moyennes des champs necessaires
265	!!AA en mode traceur off-line
266	!!======================================================================
267	!! CLEFS CPP POUR LES IO
268	!! =====================
269	#define histNMC
270	!!======================================================================
271	!! modif ( P. Le Van , 12/10/98 )
272	!!
273	!! Arguments:
274	!!
275	!! nlon----input-I-nombre de points horizontaux
276	!! nlev----input-I-nombre de couches verticales, doit etre egale a klev
277	!! debut---input-L-variable logique indiquant le premier passage
278	!! lafin---input-L-variable logique indiquant le dernier passage
279	!! jD_cur -R-jour courant a l'appel de la physique (jour julien)
280	!! jH_cur -R-heure courante a l'appel de la physique (jour julien)
281	!! pdtphys-input-R-pas d'integration pour la physique (seconde)
282	!! paprs---input-R-pression pour chaque inter-couche (en Pa)
283	!! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
284	!! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
285	!! pphis---input-R-geopotentiel du sol
286	!! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
287	!! u-------input-R-vitesse dans la direction X (de O a E) en m/s
288	!! v-------input-R-vitesse Y (de S a N) en m/s
289	!! t-------input-R-temperature (K)
290	!! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
291	!! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
292	!! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
293	!! d_ql_dyn-input-R-tendance dynamique pour "ql" (kg/kg/s)
294	!! d_qs_dyn-input-R-tendance dynamique pour "qs" (kg/kg/s)
295	!! flxmass_w -input-R- flux de masse verticale
296	!! d_u-----output-R-tendance physique de "u" (m/s/s)
297	!! d_v-----output-R-tendance physique de "v" (m/s/s)
298	!! d_t-----output-R-tendance physique de "t" (K/s)
299	!! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
300	!! d_ps----output-R-tendance physique de la pression au sol
301	!!======================================================================
302	integer jjmp1
303	! parameter (jjmp1=jjm+1-1/jjm) ! => (jjmp1=nbp_lat-1/(nbp_lat-1))
304	! integer iip1
305	! parameter (iip1=iim+1)
306
307	include "regdim.h"
308	include "dimsoil.h"
309	include "clesphys.h"
310	include "thermcell.h"
311	!======================================================================
312	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
313	PARAMETER (ok_cvl=.TRUE.)
314	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
315	PARAMETER (ok_gust=.FALSE.)
316	integer iflag_radia ! active ou non le rayonnement (MPL)
317	save iflag_radia
318	!$OMP THREADPRIVATE(iflag_radia)
319	!======================================================================
320	LOGICAL check ! Verifier la conservation du modele en eau
321	PARAMETER (check=.FALSE.)
322	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
323	PARAMETER (ok_stratus=.FALSE.)
324	!======================================================================
325	REAL amn, amx
326	INTEGER igout
327	!======================================================================
328	! Clef controlant l'activation du cycle diurne:
329	! en attente du codage des cles par Fred
330	INTEGER iflag_cycle_diurne
331	PARAMETER (iflag_cycle_diurne=1)
332	!======================================================================
333	! Modele thermique du sol, a activer pour le cycle diurne:
334	!cc LOGICAL soil_model
335	!cc PARAMETER (soil_model=.FALSE.)
336	!======================================================================
337	! Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
338	! le calcul du rayonnement est celle apres la precipitation des nuages.
339	! Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
340	! la condensation et la precipitation. Cette cle augmente les impacts
341	! radiatifs des nuages.
342	!cc LOGICAL new_oliq
343	!cc PARAMETER (new_oliq=.FALSE.)
344	!======================================================================
345	! Clefs controlant deux parametrisations de l'orographie:
346	!c LOGICAL ok_orodr
347	!cc PARAMETER (ok_orodr=.FALSE.)
348	!cc LOGICAL ok_orolf
349	!cc PARAMETER (ok_orolf=.FALSE.)
350	!======================================================================
351	LOGICAL ok_journe ! sortir le fichier journalier
352	save ok_journe
353	!$OMP THREADPRIVATE(ok_journe)
354	!
355	LOGICAL ok_mensuel ! sortir le fichier mensuel
356	save ok_mensuel
357	!$OMP THREADPRIVATE(ok_mensuel)
358	!
359	LOGICAL ok_instan ! sortir le fichier instantane
360	save ok_instan
361	!$OMP THREADPRIVATE(ok_instan)
362	!
363	LOGICAL ok_LES ! sortir le fichier LES
364	save ok_LES
365	!$OMP THREADPRIVATE(ok_LES)
366	!
367	LOGICAL callstats ! sortir le fichier stats
368	save callstats
369	!$OMP THREADPRIVATE(callstats)
370	!
371	LOGICAL ok_region ! sortir le fichier regional
372	PARAMETER (ok_region=.FALSE.)
373	!======================================================================
374	real seuil_inversion
375	save seuil_inversion
376	!$OMP THREADPRIVATE(seuil_inversion)
377	integer iflag_ratqs
378	save iflag_ratqs
379	!$OMP THREADPRIVATE(iflag_ratqs)
380	real facteur
381
382	REAL wmax_th(klon)
383	REAL tau_overturning_th(klon)
384
385	integer lmax_th(klon)
386	integer limbas(klon)
387	real ratqscth(klon,klev)
388	real ratqsdiff(klon,klev)
389	real zqsatth(klon,klev)
390
391	!======================================================================
392	!
393	INTEGER ivap ! indice de traceurs pour vapeur d'eau
394	PARAMETER (ivap=1)
395	INTEGER iliq ! indice de traceurs pour eau liquide
396	PARAMETER (iliq=2)
397	!CR: on ajoute la phase glace
398	INTEGER isol ! indice de traceurs pour eau glace
399	PARAMETER (isol=3)
400	!
401	!
402	! Variables argument:
403	!
404	INTEGER nlon
405	INTEGER nlev
406	REAL,INTENT(IN) :: pdtphys_
407	! NB: pdtphys to be used in physics is in time_phylmdz_mod
408	LOGICAL debut, lafin
409	REAL paprs(klon,klev+1)
410	REAL pplay(klon,klev)
411	REAL pphi(klon,klev)
412	REAL pphis(klon)
413	REAL presnivs(klev)
414	!JLD REAL znivsig(klev)
415	!JLD real pir
416
417	REAL u(klon,klev)
418	REAL v(klon,klev)
419
420	REAL, intent(in):: rot(klon, klev)
421	! relative vorticity, in s-1, needed for frontal waves
422
423	REAL t(klon,klev),thetal(klon,klev)
424	! thetal: ligne suivante a decommenter si vous avez les fichiers
425	! MPL 20130625
426	! fth_fonctions.F90 et parkind1.F90
427	! sinon thetal=theta
428	! REAL fth_thetae,fth_thetav,fth_thetal
429	REAL qx(klon,klev,nqtot)
430	REAL flxmass_w(klon,klev)
431	REAL d_u(klon,klev)
432	REAL d_v(klon,klev)
433	REAL d_t(klon,klev)
434	REAL d_qx(klon,klev,nqtot)
435	REAL d_ps(klon)
436	! variables pour tend_to_tke
437	REAL duadd(klon,klev)
438	REAL dvadd(klon,klev)
439	REAL dtadd(klon,klev)
440
441	! Variables pour le transport convectif
442	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
443	real wght_cvfd(klon,klev)
444	#ifndef CPP_XIOS
445	REAL, SAVE :: missing_val
446	#endif
447	! Variables pour le lessivage convectif
448	! RomP >>>
449	real phi2(klon,klev,klev)
450	real d1a(klon,klev),dam(klon,klev)
451	real ev(klon,klev)
452	real clw(klon,klev),elij(klon,klev,klev)
453	real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
454	! RomP <<<
455	!IM definition dynamique o_trac dans phys_output_open
456	! type(ctrl_out) :: o_trac(nqtot)
457
458	! variables a une pression donnee
459	!
460	include "declare_STDlev.h"
461	!
462	!
463	include "radopt.h"
464	!
465	!
466	INTEGER debug
467	INTEGER n
468	!ym INTEGER npoints
469	!ym PARAMETER(npoints=klon)
470	!
471	INTEGER nregISCtot
472	PARAMETER(nregISCtot=1)
473	!
474	! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties
475	! sur 1 region rectangulaire y compris pour 1 point
476	! imin_debut : indice minimum de i; nbpti : nombre de points en
477	! direction i (longitude)
478	! jmin_debut : indice minimum de j; nbptj : nombre de points en
479	! direction j (latitude)
480	!JLD INTEGER imin_debut, nbpti
481	!JLD INTEGER jmin_debut, nbptj
482	!IM: region='3d' <==> sorties en global
483	CHARACTER*3 region
484	PARAMETER(region='3d')
485	logical ok_hf
486	!
487	save ok_hf
488	!$OMP THREADPRIVATE(ok_hf)
489
490	INTEGER,PARAMETER :: longcles=20
491	REAL,SAVE :: clesphy0(longcles)
492	!$OMP THREADPRIVATE(clesphy0)
493	!
494	! Variables propres a la physique
495	INTEGER itap
496	SAVE itap ! compteur pour la physique
497	!$OMP THREADPRIVATE(itap)
498
499	INTEGER, SAVE :: abortphy=0 ! Reprere si on doit arreter en fin de phys
500	!$OMP THREADPRIVATE(abortphy)
501	!
502	REAL,save :: solarlong0
503	!$OMP THREADPRIVATE(solarlong0)
504
505	!
506	! Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
507	!
508	!IM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
509	REAL zulow(klon),zvlow(klon)
510	!
511	INTEGER igwd,idx(klon),itest(klon)
512	!
513	! REAL,allocatable,save :: run_off_lic_0(:)
514	! !$OMP THREADPRIVATE(run_off_lic_0)
515	!ym SAVE run_off_lic_0
516	!KE43
517	! Variables liees a la convection de K. Emanuel (sb):
518	!
519	REAL bas, top ! cloud base and top levels
520	SAVE bas
521	SAVE top
522	!$OMP THREADPRIVATE(bas, top)
523	!------------------------------------------------------------------
524	! Upmost level reached by deep convection and related variable (jyg)
525	!
526	INTEGER izero
527	INTEGER k_upper_cv
528	!------------------------------------------------------------------
529	!
530	!==========================================================================
531	!CR04.12.07: on ajoute les nouvelles variables du nouveau schema
532	!de convection avec poches froides
533	! Variables li\'ees \`a la poche froide (jyg)
534
535	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
536	!
537	REAL wape_prescr, fip_prescr
538	INTEGER it_wape_prescr
539	SAVE wape_prescr, fip_prescr, it_wape_prescr
540	!$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
541	!
542	! variables supplementaires de concvl
543	REAL Tconv(klon,klev)
544	REAL sij(klon,klev,klev)
545	!! !
546	!! ! variables pour tester la conservation de l'energie dans concvl
547	!! REAL, DIMENSION(klon,klev) :: d_t_con_sat
548	!! REAL, DIMENSION(klon,klev) :: d_q_con_sat
549	!! REAL, DIMENSION(klon,klev) :: dql_sat
550
551	real, save :: alp_bl_prescr=0.
552	real, save :: ale_bl_prescr=0.
553
554	real, save :: wake_s_min_lsp=0.1
555
556	!$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
557	!$OMP THREADPRIVATE(wake_s_min_lsp)
558
559
560	real ok_wk_lsp(klon)
561
562	!RC
563	! Variables li\'ees \`a la poche froide (jyg et rr)
564
565	INTEGER, SAVE :: iflag_wake_tend ! wake: if =0, then wake state variables are
566	! updated within calwake
567	!$OMP THREADPRIVATE(iflag_wake_tend)
568	REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
569	REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
570
571	REAL wake_dth(klon,klev) ! wake : temp pot difference
572
573	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta
574	! transported by LS omega
575	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of
576	! large scale omega
577	REAL wake_dtKE(klon,klev) ! Wake : differential heating
578	! (wake - unpertubed) CONV
579	REAL wake_dqKE(klon,klev) ! Wake : differential moistening
580	! (wake - unpertubed) CONV
581	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
582	REAL wake_spread(klon,klev) ! spreading term in wake_delt
583	!
584	!pourquoi y'a pas de save??
585	!
586	!!! INTEGER, SAVE, DIMENSION(klon) :: wake_k
587	!!! !$OMP THREADPRIVATE(wake_k)
588	!
589	!jyg<
590	!cc REAL wake_pe(klon) ! Wake potential energy - WAPE
591	!>jyg
592
593	REAL wake_gfl(klon) ! Gust Front Length
594	!!! REAL wake_dens(klon) ! moved to phys_state_var_mod
595	!
596	!
597	REAL dt_dwn(klon,klev)
598	REAL dq_dwn(klon,klev)
599	REAL M_dwn(klon,klev)
600	REAL M_up(klon,klev)
601	REAL dt_a(klon,klev)
602	REAL dq_a(klon,klev)
603	REAL d_t_adjwk(klon,klev) !jyg
604	REAL d_q_adjwk(klon,klev) !jyg
605	LOGICAL,SAVE :: ok_adjwk=.FALSE.
606	!$OMP THREADPRIVATE(ok_adjwk)
607	INTEGER,SAVE :: iflag_adjwk=0 !jyg
608	!$OMP THREADPRIVATE(iflag_adjwk) !jyg
609	REAL,SAVE :: oliqmax=999.,oicemax=999.
610	!$OMP THREADPRIVATE(oliqmax,oicemax)
611	REAL, SAVE :: alp_offset
612	!$OMP THREADPRIVATE(alp_offset)
613
614	!
615	!RR:fin declarations poches froides
616	!==========================================================================
617
618	REAL ztv(klon,klev),ztva(klon,klev)
619	REAL zpspsk(klon,klev)
620	REAL ztla(klon,klev),zqla(klon,klev)
621	REAL zthl(klon,klev)
622
623	!cc nrlmd le 10/04/2012
624
625	!--------Stochastic Boundary Layer Triggering: ALE_BL--------
626	!---Propri\'et\'es du thermiques au LCL
627	real zlcl_th(klon) ! Altitude du LCL calcul\'e
628	! continument (pcon dans
629	! thermcell_main.F90)
630	real fraca0(klon) ! Fraction des thermiques au LCL
631	real w0(klon) ! Vitesse des thermiques au LCL
632	real w_conv(klon) ! Vitesse verticale de grande \'echelle au LCL
633	real tke0(klon,klev+1) ! TKE au d\'ebut du pas de temps
634	real therm_tke_max0(klon) ! TKE dans les thermiques au LCL
635	real env_tke_max0(klon) ! TKE dans l'environnement au LCL
636
637	!JLD !---D\'eclenchement stochastique
638	!JLD integer :: tau_trig(klon)
639
640	REAL,SAVE :: random_notrig_max=1.
641	!$OMP THREADPRIVATE(random_notrig_max)
642
643	!--------Statistical Boundary Layer Closure: ALP_BL--------
644	!---Profils de TKE dans et hors du thermique
645	real therm_tke_max(klon,klev) ! Profil de TKE dans les thermiques
646	real env_tke_max(klon,klev) ! Profil de TKE dans l'environnement
647
648	!-------Activer les tendances de TKE due a l'orograp??ie---------
649	INTEGER, SAVE :: addtkeoro
650	!$OMP THREADPRIVATE(addtkeoro)
651	REAL, SAVE :: alphatkeoro
652	!$OMP THREADPRIVATE(alphatkeoro)
653	LOGICAL, SAVE :: smallscales_tkeoro
654	!$OMP THREADPRIVATE(smallscales_tkeoro)
655
656
657
658	!cc fin nrlmd le 10/04/2012
659
660	! Variables locales pour la couche limite (al1):
661	!
662	!Al1 REAL pblh(klon) ! Hauteur de couche limite
663	!Al1 SAVE pblh
664	!34EK
665	!
666	! Variables locales:
667	!
668	!AA
669	!AA Pour phytrac
670	REAL u1(klon) ! vents dans la premiere couche U
671	REAL v1(klon) ! vents dans la premiere couche V
672
673	!@$$ LOGICAL offline ! Controle du stockage ds "physique"
674	!@$$ PARAMETER (offline=.false.)
675	!@$$ INTEGER physid
676	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
677	REAL frac_nucl(klon,klev) ! idem (nucleation)
678	! RomP >>>
679	REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
680	! RomP <<<
681	REAL :: calday
682
683	!IM cf FH pour Tiedtke 080604
684	REAL rain_tiedtke(klon),snow_tiedtke(klon)
685	!
686	!IM 050204 END
687	REAL devap(klon) ! evaporation et sa derivee
688	REAL dsens(klon) ! chaleur sensible et sa derivee
689
690	!
691	! Conditions aux limites
692	!
693	!
694	REAL :: day_since_equinox
695	! Date de l'equinoxe de printemps
696	INTEGER, parameter :: mth_eq=3, day_eq=21
697	REAL :: jD_eq
698
699	LOGICAL, parameter :: new_orbit = .true.
700
701	!
702	INTEGER lmt_pas
703	SAVE lmt_pas ! frequence de mise a jour
704	!$OMP THREADPRIVATE(lmt_pas)
705	real zmasse(klon, nbp_lev),exner(klon, nbp_lev)
706	! (column-density of mass of air in a cell, in kg m-2)
707	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
708
709	!IM sorties
710	REAL un_jour
711	PARAMETER(un_jour=86400.)
712	INTEGER itapm1 !pas de temps de la physique du(es) mois precedents
713	SAVE itapm1 !mis a jour le dernier pas de temps du mois en cours
714	!$OMP THREADPRIVATE(itapm1)
715	!======================================================================
716	!
717	! Declaration des procedures appelees
718	!
719	EXTERNAL angle ! calculer angle zenithal du soleil
720	EXTERNAL alboc ! calculer l'albedo sur ocean
721	EXTERNAL ajsec ! ajustement sec
722	EXTERNAL conlmd ! convection (schema LMD)
723	!KE43
724	EXTERNAL conema3 ! convect4.3
725	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
726	!AA
727	! JBM (3/14) fisrtilp_tr not loaded
728	! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
729	! ! stockage des coefficients necessaires au
730	! ! lessivage OFF-LINE et ON-LINE
731	EXTERNAL hgardfou ! verifier les temperatures
732	EXTERNAL nuage ! calculer les proprietes radiatives
733	!C EXTERNAL o3cm ! initialiser l'ozone
734	EXTERNAL orbite ! calculer l'orbite terrestre
735	EXTERNAL phyetat0 ! lire l'etat initial de la physique
736	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
737	EXTERNAL suphel ! initialiser certaines constantes
738	EXTERNAL transp ! transport total de l'eau et de l'energie
739	!IM
740	EXTERNAL haut2bas !variables de haut en bas
741	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
742	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
743	! EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
744	! EXTERNAL moyglo_aire
745	! moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
746	! par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
747	!
748	!
749	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
750	! Local variables
751	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
752	!
753	REAL rhcl(klon,klev) ! humiditi relative ciel clair
754	REAL dialiq(klon,klev) ! eau liquide nuageuse
755	REAL diafra(klon,klev) ! fraction nuageuse
756	REAL cldliq(klon,klev) ! eau liquide nuageuse
757	!
758	!XXX PB
759	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
760	!
761	REAL zxfluxt(klon, klev)
762	REAL zxfluxq(klon, klev)
763	REAL zxfluxu(klon, klev)
764	REAL zxfluxv(klon, klev)
765
766	! Le rayonnement n'est pas calcule tous les pas, il faut donc
767	! sauvegarder les sorties du rayonnement
768	!ym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
769	!ym SAVE sollwdownclr, toplwdown, toplwdownclr
770	!ym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
771	!
772	INTEGER itaprad
773	SAVE itaprad
774	!$OMP THREADPRIVATE(itaprad)
775	!
776	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
777	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
778	!
779	#ifdef INCA
780	REAL zxsnow_dummy(klon)
781	#endif
782	REAL zsav_tsol(klon)
783	!
784	REAL dist, rmu0(klon), fract(klon)
785	REAL zrmu0(klon), zfract(klon)
786	REAL zdtime, zdtime1, zdtime2, zlongi
787	!
788	REAL qcheck
789	REAL z_avant(klon), z_apres(klon), z_factor(klon)
790	LOGICAL zx_ajustq
791	!
792	REAL za
793	REAL zx_t, zx_qs, zdelta, zcor
794	real zqsat(klon,klev)
795	!
796	INTEGER i, k, iq, j, nsrf, ll, l
797	!
798	REAL t_coup
799	PARAMETER (t_coup=234.0)
800
801	!ym A voir plus tard !!
802	!ym REAL zx_relief(iim,jjmp1)
803	!ym REAL zx_aire(iim,jjmp1)
804	!
805	! Grandeurs de sorties
806	REAL s_capCL(klon)
807	REAL s_oliqCL(klon), s_cteiCL(klon)
808	REAL s_trmb1(klon), s_trmb2(klon)
809	REAL s_trmb3(klon)
810
811	! La convection n'est pas calculee tous les pas, il faut donc
812	! sauvegarder les sorties de la convection
813	!ym SAVE
814	!ym SAVE
815	!ym SAVE
816	!
817	INTEGER itapcv, itapwk
818	SAVE itapcv, itapwk
819	!$OMP THREADPRIVATE(itapcv, itapwk)
820
821	!KE43
822	! Variables locales pour la convection de K. Emanuel (sb):
823
824	REAL tvp(klon,klev) ! virtual temp of lifted parcel
825	CHARACTER*40 capemaxcels !max(CAPE)
826
827	REAL rflag(klon) ! flag fonctionnement de convect
828	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
829
830	! -- convect43:
831	INTEGER ntra ! nb traceurs pour convect4.3
832	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
833	REAL dplcldt(klon), dplcldr(klon)
834	!? . condm_con(klon,klev),conda_con(klon,klev),
835	!? . mr_con(klon,klev),ep_con(klon,klev)
836	!? . ,sadiab(klon,klev),wadiab(klon,klev)
837	! --
838	!34EK
839	!
840	! Variables du changement
841	!
842	! con: convection
843	! lsc: condensation a grande echelle (Large-Scale-Condensation)
844	! ajs: ajustement sec
845	! eva: evaporation de l'eau liquide nuageuse
846	! vdf: couche limite (Vertical DiFfusion)
847	!
848	! tendance nulles
849	REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
850	REAL, dimension(klon) :: dsig0, ddens0
851	INTEGER, dimension(klon) :: wkoccur1
852	! tendance buffer pour appel de add_phys_tend
853	REAL, DIMENSION(klon,klev) :: d_q_ch4_dtime
854	!
855	! Flag pour pouvoir ne pas ajouter les tendances.
856	! Par defaut, les tendances doivente etre ajoutees et
857	! flag_inhib_tend = 0
858	! flag_inhib_tend > 0 : tendances non ajoutees, avec un nombre
859	! croissant de print quand la valeur du flag augmente
860	!!! attention, ce flag doit etre change avec prudence !!!
861	INTEGER :: flag_inhib_tend = 0 ! 0 is the default value
862	!! INTEGER :: flag_inhib_tend = 2
863
864	!
865	!********************************************************
866	! declarations
867
868	!********************************************************
869	!IM 081204 END
870	!
871	REAL pen_u(klon,klev), pen_d(klon,klev)
872	REAL pde_u(klon,klev), pde_d(klon,klev)
873	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
874	!
875	REAL ratqsc(klon,klev)
876	real ratqsbas,ratqshaut,tau_ratqs
877	save ratqsbas,ratqshaut,tau_ratqs
878	!$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
879	REAL, SAVE :: ratqsp0=50000., ratqsdp=20000.
880	!$OMP THREADPRIVATE(ratqsp0, ratqsdp)
881
882	! Parametres lies au nouveau schema de nuages (SB, PDF)
883	real fact_cldcon
884	real facttemps
885	logical ok_newmicro
886	save ok_newmicro
887	!$OMP THREADPRIVATE(ok_newmicro)
888	!real ref_liq_pi(klon,klev), ref_ice_pi(klon,klev)
889	save fact_cldcon,facttemps
890	!$OMP THREADPRIVATE(fact_cldcon,facttemps)
891
892	integer iflag_cld_th
893	save iflag_cld_th
894	!$OMP THREADPRIVATE(iflag_cld_th)
895	!IM logical ptconv(klon,klev) !passe dans phys_local_var_mod
896	!IM cf. AM 081204 BEG
897	logical ptconvth(klon,klev)
898	!IM cf. AM 081204 END
899	!
900	! Variables liees a l'ecriture de la bande histoire physique
901	!
902	!======================================================================
903	!
904
905	!
906	!JLD integer itau_w ! pas de temps ecriture = itap + itau_phy
907	!
908	!
909	! Variables locales pour effectuer les appels en serie
910	!
911	!IM RH a 2m (la surface)
912	REAL Lheat
913
914	INTEGER length
915	PARAMETER ( length = 100 )
916	REAL tabcntr0( length )
917	!
918	!JLD INTEGER ndex2d(nbp_lon*nbp_lat)
919	!IM
920	!
921	!IM AMIP2 BEG
922	!JLD REAL moyglo, mountor
923	!IM 141004 BEG
924	REAL zustrdr(klon), zvstrdr(klon)
925	REAL zustrli(klon), zvstrli(klon)
926	REAL zustrph(klon), zvstrph(klon)
927	REAL aam, torsfc
928	!IM 141004 END
929	!IM 190504 BEG
930	! INTEGER imp1jmp1
931	! PARAMETER(imp1jmp1=(iim+1)*jjmp1)
932	!ym A voir plus tard
933	! REAL zx_tmp((nbp_lon+1)*nbp_lat)
934	! REAL airedyn(nbp_lon+1,nbp_lat)
935	!IM 190504 END
936	!JLD LOGICAL ok_msk
937	!JLD REAL msk(klon)
938	!ym A voir plus tard
939	!ym REAL zm_wo(jjmp1, klev)
940	!IM AMIP2 END
941	!
942	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
943	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
944	!JLD REAL zx_tmp_2d(nbp_lon,nbp_lat)
945	!JLD REAL zx_lon(nbp_lon,nbp_lat)
946	!JLD REAL zx_lat(nbp_lon,nbp_lat)
947	!
948	INTEGER nid_ctesGCM
949	SAVE nid_ctesGCM
950	!$OMP THREADPRIVATE(nid_ctesGCM)
951	!
952	!IM 280405 BEG
953	! INTEGER nid_bilKPins, nid_bilKPave
954	! SAVE nid_bilKPins, nid_bilKPave
955	! !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
956	!
957	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
958	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
959	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
960	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
961	!
962	!JLD REAL zjulian
963	!JLD SAVE zjulian
964	!JLD!$OMP THREADPRIVATE(zjulian)
965
966	!JLD INTEGER nhori, nvert
967	!JLD REAL zsto
968	!JLD REAL zstophy, zout
969
970	character*20 modname
971	character*80 abort_message
972	logical, save :: ok_sync, ok_sync_omp
973	!$OMP THREADPRIVATE(ok_sync)
974	real date0
975
976	! essai writephys
977	integer fid_day, fid_mth, fid_ins
978	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
979	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
980	parameter (prof2d_on = 1, prof3d_on = 2, &
981	prof2d_av = 3, prof3d_av = 4)
982	REAL ztsol(klon)
983	REAL q2m(klon,nbsrf) ! humidite a 2m
984
985	!IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
986	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
987	CHARACTER*40 tinst, tave
988	REAL cldtaupi(klon,klev) ! Cloud optical thickness for
989	! pre-industrial (pi) aerosols
990
991	INTEGER :: naero
992	! Aerosol optical properties
993	CHARACTER*4, DIMENSION(naero_grp) :: rfname
994	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass
995	! concentration
996	! for all soluble
997	! aerosols[ug/m3]
998	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi
999	! - " - (pre-industrial value)
1000
1001	! Parameters
1002	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1003	LOGICAL ok_alw ! Apply aerosol LW effect or not
1004	LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
1005	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1006	SAVE ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1
1007	!$OMP THREADPRIVATE(ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1)
1008	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
1009	! false : lecture des aerosol dans un fichier
1010	!$OMP THREADPRIVATE(aerosol_couple)
1011	INTEGER, SAVE :: flag_aerosol
1012	!$OMP THREADPRIVATE(flag_aerosol)
1013	LOGICAL, SAVE :: flag_bc_internal_mixture
1014	!$OMP THREADPRIVATE(flag_bc_internal_mixture)
1015	LOGICAL, SAVE :: new_aod
1016	!$OMP THREADPRIVATE(new_aod)
1017	!
1018	!--STRAT AEROSOL
1019	INTEGER, SAVE :: flag_aerosol_strat
1020	!$OMP THREADPRIVATE(flag_aerosol_strat)
1021	!c-fin STRAT AEROSOL
1022	!
1023	! Declaration des constantes et des fonctions thermodynamiques
1024	!
1025	LOGICAL,SAVE :: first=.true.
1026	!$OMP THREADPRIVATE(first)
1027
1028	! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
1029	! Note that pressure vectors are in Pa and in stricly ascending order
1030	INTEGER,SAVE :: read_climoz ! Read ozone climatology
1031	! (let it keep the default OpenMP shared attribute)
1032	! Allowed values are 0, 1 and 2
1033	! 0: do not read an ozone climatology
1034	! 1: read a single ozone climatology that will be used day and night
1035	! 2: read two ozone climatologies, the average day and night
1036	! climatology and the daylight climatology
1037	INTEGER,SAVE :: ncid_climoz ! NetCDF file identifier
1038	REAL, POINTER, SAVE :: press_cen_climoz(:) ! Pressure levels
1039	REAL, POINTER, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
1040	REAL, POINTER, SAVE :: time_climoz(:) ! Time vector
1041	CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
1042	= ["tro3 ","tro3_daylight"]
1043	! vars_climoz(1:read_climoz): variables names in climoz file.
1044	! vars_climoz(1:read_climoz-2) if read_climoz>2 (temporary)
1045	REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
1046	! the ozone fields, old method.
1047
1048	include "YOMCST.h"
1049	include "YOETHF.h"
1050	include "FCTTRE.h"
1051	!IM 100106 BEG : pouvoir sortir les ctes de la physique
1052	include "conema3.h"
1053	include "fisrtilp.h"
1054	include "nuage.h"
1055	include "compbl.h"
1056	!IM 100106 END : pouvoir sortir les ctes de la physique
1057	!
1058	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1059	! Declarations pour Simulateur COSP
1060	!============================================================
1061	#ifdef CPP_COSP
1062	real :: mr_ozone(klon,klev)
1063	#endif
1064	!IM stations CFMIP
1065	INTEGER, SAVE :: nCFMIP
1066	!$OMP THREADPRIVATE(nCFMIP)
1067	INTEGER, PARAMETER :: npCFMIP=120
1068	INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
1069	REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
1070	!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
1071	INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
1072	REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
1073	!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
1074	INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
1075	!$OMP THREADPRIVATE(iGCM, jGCM)
1076	logical, dimension(nfiles) :: phys_out_filestations
1077	logical, parameter :: lNMC=.FALSE.
1078
1079	!IM betaCRF
1080	REAL, SAVE :: pfree, beta_pbl, beta_free
1081	!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
1082	REAL, SAVE :: lon1_beta, lon2_beta, lat1_beta, lat2_beta
1083	!$OMP THREADPRIVATE(lon1_beta, lon2_beta, lat1_beta, lat2_beta)
1084	LOGICAL, SAVE :: mskocean_beta
1085	!$OMP THREADPRIVATE(mskocean_beta)
1086	REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et
1087	! cldemirad pour evaluer les
1088	! retros liees aux CRF
1089	REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique
1090	! pour radlwsw pour
1091	! tester "CRF off"
1092	REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique
1093	! pour radlwsw pour
1094	! tester "CRF off"
1095	REAL, dimension(klon, klev) :: cldemirad ! emissivite pour
1096	! radlwsw pour tester
1097	! "CRF off"
1098	REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse
1099
1100	INTEGER :: nbtr_tmp ! Number of tracer inside concvl
1101	REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
1102	REAL, dimension(klon,klev) :: ch_in ! Condensed humidity entering in phytrac (eau liquide)
1103	integer iostat
1104
1105	REAL zzz
1106	!albedo SB >>>
1107	real,dimension(6),save :: SFRWL
1108	!albedo SB <<<
1109
1110	!--OB variables for mass fixer (hard coded for now)
1111	logical, parameter :: mass_fixer=.false.
1112	real qql1(klon),qql2(klon),corrqql
1113
1114	! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
1115	jjmp1=nbp_lat
1116
1117	!======================================================================
1118	! Gestion calendrier : mise a jour du module phys_cal_mod
1119	!
1120	pdtphys=pdtphys_
1121	CALL update_time(pdtphys)
1122
1123	!======================================================================
1124	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1125	! Utilise notamment en 1D mais peut etre active egalement en 3D
1126	! en imposant la valeur de igout.
1127	!======================================================================d
1128	IF (prt_level.ge.1) THEN
1129	igout=klon/2+1/klon
1130	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1131	write(lunout,*) 'igout, lat, lon ',igout, latitude_deg(igout), &
1132	longitude_deg(igout)
1133	write(lunout,*) &
1134	'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1135	write(lunout,*) &
1136	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1137
1138	write(lunout,*) 'paprs, play, phi, u, v, t'
1139	DO k=1,klev
1140	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
1141	u(igout,k),v(igout,k),t(igout,k)
1142	ENDDO
1143	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1144	DO k=1,klev
1145	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1146	ENDDO
1147	ENDIF
1148
1149	! Quick check on pressure levels:
1150	call assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
1151	"physiq_mod paprs bad order")
1152
1153	IF (first) THEN
1154	!CR:nvelles variables convection/poches froides
1155
1156	print*, '================================================='
1157	print*, 'Allocation des variables locales et sauvegardees'
1158	CALL phys_local_var_init
1159	!
1160	pasphys=pdtphys
1161	! appel a la lecture du run.def physique
1162	CALL conf_phys(ok_journe, ok_mensuel, &
1163	ok_instan, ok_hf, &
1164	ok_LES, &
1165	callstats, &
1166	solarlong0,seuil_inversion, &
1167	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
1168	iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
1169	ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple, &
1170	flag_aerosol, flag_aerosol_strat, new_aod, &
1171	flag_bc_internal_mixture, bl95_b0, bl95_b1, &
1172	! nv flags pour la convection et les
1173	! poches froides
1174	read_climoz, &
1175	alp_offset)
1176	CALL phys_state_var_init(read_climoz)
1177	CALL phys_output_var_init
1178	print*, '================================================='
1179	!
1180	!CR: check sur le nb de traceurs de l eau
1181	IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
1182	WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
1183	'(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
1184	STOP
1185	ENDIF
1186
1187	dnwd0=0.0
1188	ftd=0.0
1189	fqd=0.0
1190	cin=0.
1191	!ym Attention pbase pas initialise dans concvl !!!!
1192	pbase=0
1193	!IM 180608
1194
1195	itau_con=0
1196	first=.false.
1197
1198	ENDIF ! first
1199
1200	!ym => necessaire pour iflag_con != 2
1201	pmfd(:,:) = 0.
1202	pen_u(:,:) = 0.
1203	pen_d(:,:) = 0.
1204	pde_d(:,:) = 0.
1205	pde_u(:,:) = 0.
1206	aam=0.
1207	d_t_adjwk(:,:)=0
1208	d_q_adjwk(:,:)=0
1209
1210	alp_bl_conv(:)=0.
1211
1212	torsfc=0.
1213	forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1214
1215	modname = 'physiq'
1216
1217	IF (debut) THEN
1218	CALL suphel ! initialiser constantes et parametres phys.
1219	CALL getin_p('random_notrig_max',random_notrig_max)
1220	CALL getin_p('ok_adjwk',ok_adjwk)
1221	IF (ok_adjwk) iflag_adjwk=2 ! for compatibility with older versions
1222	! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
1223	! 1 => convective adjustment but state variables are unchanged
1224	! 2 => convective adjustment and state variables are changed
1225	CALL getin_p('iflag_adjwk',iflag_adjwk)
1226	CALL getin_p('oliqmax',oliqmax)
1227	CALL getin_p('oicemax',oicemax)
1228	CALL getin_p('ratqsp0',ratqsp0)
1229	CALL getin_p('ratqsdp',ratqsdp)
1230	iflag_wake_tend = 0
1231	CALL getin_p('iflag_wake_tend',iflag_wake_tend)
1232	ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
1233	! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
1234	CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
1235	fl_ebil = 0 ! by default, conservation diagnostics are desactivated
1236	CALL getin_p('fl_ebil',fl_ebil)
1237	fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
1238	CALL getin_p('fl_cor_ebil',fl_cor_ebil)
1239	iflag_phytrac = 1 ! by default we do want to call phytrac
1240	CALL getin_p('iflag_phytrac',iflag_phytrac)
1241	ENDIF
1242
1243	IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
1244
1245
1246	!======================================================================
1247	! Gestion calendrier : mise a jour du module phys_cal_mod
1248	!
1249	! CALL phys_cal_update(jD_cur,jH_cur)
1250
1251	!
1252	! Si c'est le debut, il faut initialiser plusieurs choses
1253	! ********
1254	!
1255	IF (debut) THEN
1256	!rv CRinitialisation de wght_th et lalim_conv pour la
1257	!definition de la couche alimentation de la convection a partir
1258	!des caracteristiques du thermique
1259	wght_th(:,:)=1.
1260	lalim_conv(:)=1
1261	!RC
1262	ustar(:,:)=0.
1263	! u10m(:,:)=0.
1264	! v10m(:,:)=0.
1265	rain_con(:)=0.
1266	snow_con(:)=0.
1267	topswai(:)=0.
1268	topswad(:)=0.
1269	solswai(:)=0.
1270	solswad(:)=0.
1271
1272	wmax_th(:)=0.
1273	tau_overturning_th(:)=0.
1274
1275	IF (type_trac == 'inca') THEN
1276	! jg : initialisation jusqu'au ces variables sont dans restart
1277	ccm(:,:,:) = 0.
1278	tau_aero(:,:,:,:) = 0.
1279	piz_aero(:,:,:,:) = 0.
1280	cg_aero(:,:,:,:) = 0.
1281
1282	config_inca='none' ! default
1283	CALL getin_p('config_inca',config_inca)
1284
1285	ELSE
1286	config_inca='none' ! default
1287	ENDIF
1288
1289	IF (aerosol_couple .AND. (config_inca /= "aero" &
1290	.AND. config_inca /= "aeNP ")) THEN
1291	abort_message &
1292	= 'if aerosol_couple is activated, config_inca need to be ' &
1293	// 'aero or aeNP'
1294	CALL abort_physic (modname,abort_message,1)
1295	ENDIF
1296
1297
1298
1299	rnebcon0(:,:) = 0.0
1300	clwcon0(:,:) = 0.0
1301	rnebcon(:,:) = 0.0
1302	clwcon(:,:) = 0.0
1303
1304	!
1305	print*,'iflag_coupl,iflag_clos,iflag_wake', &
1306	iflag_coupl,iflag_clos,iflag_wake
1307	print*,'iflag_CYCLE_DIURNE', iflag_cycle_diurne
1308	!
1309	IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
1310	abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
1311	CALL abort_physic (modname,abort_message,1)
1312	ENDIF
1313	!
1314	!
1315	! Initialiser les compteurs:
1316	!
1317	itap = 0
1318	itaprad = 0
1319	itapcv = 0
1320	itapwk = 0
1321
1322	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1323	!! Un petit travail \`a faire ici.
1324	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1325
1326	IF (iflag_pbl>1) THEN
1327	PRINT*, "Using method MELLOR&YAMADA"
1328	ENDIF
1329
1330	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1331	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
1332	! phylmd plutot que dyn3d
1333	! Attention : la version precedente n'etait pas tres propre.
1334	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1335	! pour obtenir le meme resultat.
1336	!jyg for fh<
1337	!! dtime=pdtphys
1338	dtime=NINT(pdtphys)
1339	WRITE(lunout,*) 'Pas de temps dtime pdtphys ',dtime,pdtphys
1340	IF (abs(dtime-pdtphys)>1.e-10) THEN
1341	abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
1342	CALL abort_physic(modname,abort_message,1)
1343	ENDIF
1344	!>jyg
1345	IF (MOD(NINT(86400./dtime),nbapp_rad).EQ.0) THEN
1346	radpas = NINT( 86400./dtime)/nbapp_rad
1347	ELSE
1348	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1349	'multiple de nbapp_rad'
1350	WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
1351	'mais 1+1<>2'
1352	abort_message='nbre de pas de temps physique n est pas multiple ' &
1353	// 'de nbapp_rad'
1354	CALL abort_physic(modname,abort_message,1)
1355	ENDIF
1356	IF (nbapp_cv .EQ. 0) nbapp_cv=86400./dtime
1357	IF (nbapp_wk .EQ. 0) nbapp_wk=86400./dtime
1358	print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
1359	IF (MOD(NINT(86400./dtime),nbapp_cv).EQ.0) THEN
1360	cvpas = NINT( 86400./dtime)/nbapp_cv
1361	print *,'physiq, cvpas ',cvpas
1362	ELSE
1363	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1364	'multiple de nbapp_cv'
1365	WRITE(lunout,*) 'changer nbapp_cv ou alors commenter ce test ', &
1366	'mais 1+1<>2'
1367	abort_message='nbre de pas de temps physique n est pas multiple ' &
1368	// 'de nbapp_cv'
1369	call abort_physic(modname,abort_message,1)
1370	ENDIF
1371	IF (MOD(NINT(86400./dtime),nbapp_wk).EQ.0) THEN
1372	wkpas = NINT( 86400./dtime)/nbapp_wk
1373	print *,'physiq, wkpas ',wkpas
1374	ELSE
1375	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1376	'multiple de nbapp_wk'
1377	WRITE(lunout,*) 'changer nbapp_wk ou alors commenter ce test ', &
1378	'mais 1+1<>2'
1379	abort_message='nbre de pas de temps physique n est pas multiple ' &
1380	// 'de nbapp_wk'
1381	call abort_physic(modname,abort_message,1)
1382	ENDIF
1383	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1384
1385	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1386	!jyg<
1387	IF (klon_glo==1) THEN
1388	pbl_tke(:,:,is_ave) = 0.
1389	DO nsrf=1,nbsrf
1390	DO k = 1,klev+1
1391	pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
1392	+pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
1393	ENDDO
1394	ENDDO
1395	!>jyg
1396	ENDIF
1397	!IM begin
1398	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1399	,ratqs(1,1)
1400	!IM end
1401
1402
1403	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1404	!
1405	! on remet le calendrier a zero
1406	!
1407	IF (raz_date .eq. 1) THEN
1408	itau_phy = 0
1409	ENDIF
1410
1411	CALL printflag( tabcntr0,radpas,ok_journe, &
1412	ok_instan, ok_region )
1413	!
1414	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1415	WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
1416	pdtphys
1417	abort_message='Pas physique n est pas correct '
1418	! call abort_physic(modname,abort_message,1)
1419	dtime=pdtphys
1420	ENDIF
1421	IF (nlon .NE. klon) THEN
1422	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1423	klon
1424	abort_message='nlon et klon ne sont pas coherents'
1425	CALL abort_physic(modname,abort_message,1)
1426	ENDIF
1427	IF (nlev .NE. klev) THEN
1428	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1429	klev
1430	abort_message='nlev et klev ne sont pas coherents'
1431	CALL abort_physic(modname,abort_message,1)
1432	ENDIF
1433	!
1434	IF (dtime*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
1435	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1436	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1437	abort_message='Nbre d appels au rayonnement insuffisant'
1438	CALL abort_physic(modname,abort_message,1)
1439	ENDIF
1440	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1441	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1442	ok_cvl
1443	!
1444	!KE43
1445	! Initialisation pour la convection de K.E. (sb):
1446	IF (iflag_con.GE.3) THEN
1447
1448	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1449	WRITE(lunout,*) &
1450	"On va utiliser le melange convectif des traceurs qui"
1451	WRITE(lunout,*)"est calcule dans convect4.3"
1452	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1453
1454	DO i = 1, klon
1455	ema_cbmf(i) = 0.
1456	ema_pcb(i) = 0.
1457	ema_pct(i) = 0.
1458	! ema_workcbmf(i) = 0.
1459	ENDDO
1460	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1461	DO i = 1, klon
1462	ibas_con(i) = 1
1463	itop_con(i) = 1
1464	ENDDO
1465	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1466	!================================================================
1467	!CR:04.12.07: initialisations poches froides
1468	! Controle de ALE et ALP pour la fermeture convective (jyg)
1469	IF (iflag_wake>=1) THEN
1470	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1471	,alp_bl_prescr, ale_bl_prescr)
1472	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1473	! print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
1474	!
1475	! Initialize tendencies of wake state variables (for some flag values
1476	! they are not computed).
1477	d_deltat_wk(:,:) = 0.
1478	d_deltaq_wk(:,:) = 0.
1479	d_deltat_wk_gw(:,:) = 0.
1480	d_deltaq_wk_gw(:,:) = 0.
1481	d_deltat_vdf(:,:) = 0.
1482	d_deltaq_vdf(:,:) = 0.
1483	d_deltat_the(:,:) = 0.
1484	d_deltaq_the(:,:) = 0.
1485	d_deltat_ajs_cv(:,:) = 0.
1486	d_deltaq_ajs_cv(:,:) = 0.
1487	d_s_wk(:) = 0.
1488	d_dens_wk(:) = 0.
1489	ENDIF
1490
1491	! do i = 1,klon
1492	! Ale_bl(i)=0.
1493	! Alp_bl(i)=0.
1494	! enddo
1495
1496	!===================================================================
1497	!IM stations CFMIP
1498	nCFMIP=npCFMIP
1499	OPEN(98,file='npCFMIP_param.data',status='old', &
1500	form='formatted',iostat=iostat)
1501	IF (iostat == 0) THEN
1502	READ(98,*,end=998) nCFMIP
1503	998 CONTINUE
1504	CLOSE(98)
1505	CONTINUE
1506	IF(nCFMIP.GT.npCFMIP) THEN
1507	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1508	CALL abort_physic("physiq", "", 1)
1509	ELSE
1510	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1511	ENDIF
1512
1513	!
1514	ALLOCATE(tabCFMIP(nCFMIP))
1515	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1516	ALLOCATE(tabijGCM(nCFMIP))
1517	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1518	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1519	!
1520	! lecture des nCFMIP stations CFMIP, de leur numero
1521	! et des coordonnees geographiques lonCFMIP, latCFMIP
1522	!
1523	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
1524	lonCFMIP, latCFMIP)
1525	!
1526	! identification des
1527	! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
1528	! grille de LMDZ
1529	! 2) indices points tabijGCM de la grille physique 1d sur
1530	! klon points
1531	! 3) indices iGCM, jGCM de la grille physique 2d
1532	!
1533	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1534	tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1535	!
1536	ELSE
1537	ALLOCATE(tabijGCM(0))
1538	ALLOCATE(lonGCM(0), latGCM(0))
1539	ALLOCATE(iGCM(0), jGCM(0))
1540	ENDIF
1541	ELSE
1542	ALLOCATE(tabijGCM(0))
1543	ALLOCATE(lonGCM(0), latGCM(0))
1544	ALLOCATE(iGCM(0), jGCM(0))
1545	ENDIF
1546
1547	DO i=1,klon
1548	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1549	ENDDO
1550
1551	!34EK
1552	IF (ok_orodr) THEN
1553
1554	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1555	! FH sans doute a enlever de finitivement ou, si on le
1556	! garde, l'activer justement quand ok_orodr = false.
1557	! ce rugoro est utilise par la couche limite et fait double emploi
1558	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1559	! DO i=1,klon
1560	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1561	! ENDDO
1562	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1563	IF (ok_strato) THEN
1564	CALL SUGWD_strato(klon,klev,paprs,pplay)
1565	ELSE
1566	CALL SUGWD(klon,klev,paprs,pplay)
1567	ENDIF
1568
1569	DO i=1,klon
1570	zuthe(i)=0.
1571	zvthe(i)=0.
1572	IF (zstd(i).gt.10.) THEN
1573	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1574	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1575	ENDIF
1576	ENDDO
1577	ENDIF
1578	!
1579	!
1580	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1581	WRITE(lunout,*)'La frequence de lecture surface est de ', &
1582	lmt_pas
1583	!
1584	capemaxcels = 't_max(X)'
1585	t2mincels = 't_min(X)'
1586	t2maxcels = 't_max(X)'
1587	tinst = 'inst(X)'
1588	tave = 'ave(X)'
1589	!IM cf. AM 081204 BEG
1590	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1591	!IM cf. AM 081204 END
1592	!
1593	!=============================================================
1594	! Initialisation des sorties
1595	!=============================================================
1596
1597	#ifdef CPP_XIOS
1598	!--setting up swaero_diag to TRUE in XIOS case
1599	IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
1600	xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
1601	xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR. &
1602	(iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
1603	xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0")))) &
1604	!!!--for now these fields are not in the XML files so they are omitted
1605	!!! xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
1606	swaero_diag=.TRUE.
1607
1608	!--setting up dryaod_diag to TRUE in XIOS case
1609	DO naero = 1, naero_tot-1
1610	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
1611	ENDDO
1612	#endif
1613
1614	#ifdef CPP_IOIPSL
1615
1616	!$OMP MASTER
1617	! FH : if ok_sync=.true. , the time axis is written at each time step
1618	! in the output files. Only at the end in the opposite case
1619	ok_sync_omp=.false.
1620	CALL getin('ok_sync',ok_sync_omp)
1621	CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
1622	iGCM,jGCM,lonGCM,latGCM, &
1623	jjmp1,nlevSTD,clevSTD,rlevSTD, dtime,ok_veget, &
1624	type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1625	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1626	read_climoz, phys_out_filestations, &
1627	new_aod, aerosol_couple, &
1628	flag_aerosol_strat, pdtphys, paprs, pphis, &
1629	pplay, lmax_th, ptconv, ptconvth, ivap, &
1630	d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
1631	!$OMP END MASTER
1632	!$OMP BARRIER
1633	ok_sync=ok_sync_omp
1634
1635	freq_outNMC(1) = ecrit_files(7)
1636	freq_outNMC(2) = ecrit_files(8)
1637	freq_outNMC(3) = ecrit_files(9)
1638	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1639	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1640	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1641
1642	#ifndef CPP_XIOS
1643	CALL ini_paramLMDZ_phy(dtime,nid_ctesGCM)
1644	#endif
1645
1646	#endif
1647	ecrit_reg = ecrit_reg * un_jour
1648	ecrit_tra = ecrit_tra * un_jour
1649
1650	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
1651	date0 = jD_ref
1652	WRITE(,) 'physiq date0 : ',date0
1653	!
1654	!
1655	!
1656	! Prescrire l'ozone dans l'atmosphere
1657	!
1658	!
1659	!c DO i = 1, klon
1660	!c DO k = 1, klev
1661	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1662	!c ENDDO
1663	!c ENDDO
1664	!
1665	IF (type_trac == 'inca') THEN
1666	#ifdef INCA
1667	CALL VTe(VTphysiq)
1668	CALL VTb(VTinca)
1669	calday = REAL(days_elapsed) + jH_cur
1670	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1671
1672	CALL chemini( &
1673	rg, &
1674	ra, &
1675	cell_area, &
1676	latitude_deg, &
1677	longitude_deg, &
1678	presnivs, &
1679	calday, &
1680	klon, &
1681	nqtot, &
1682	nqo, &
1683	pdtphys, &
1684	annee_ref, &
1685	year_cur, &
1686	day_ref, &
1687	day_ini, &
1688	start_time, &
1689	itau_phy, &
1690	date0, &
1691	io_lon, &
1692	io_lat)
1693
1694	CALL VTe(VTinca)
1695	CALL VTb(VTphysiq)
1696	#endif
1697	ENDIF
1698	!
1699	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1700	! Nouvelle initialisation pour le rayonnement RRTM
1701	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1702
1703	CALL iniradia(klon,klev,paprs(1,1:klev+1))
1704
1705	!$omp single
1706	IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, &
1707	press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
1708	!$omp end single
1709	!
1710	!IM betaCRF
1711	pfree=70000. !Pa
1712	beta_pbl=1.
1713	beta_free=1.
1714	lon1_beta=-180.
1715	lon2_beta=+180.
1716	lat1_beta=90.
1717	lat2_beta=-90.
1718	mskocean_beta=.FALSE.
1719
1720	!albedo SB >>>
1721	select case(nsw)
1722	case(2)
1723	SFRWL(1)=0.45538747
1724	SFRWL(2)=0.54461211
1725	case(4)
1726	SFRWL(1)=0.45538747
1727	SFRWL(2)=0.32870591
1728	SFRWL(3)=0.18568763
1729	SFRWL(4)=3.02191470E-02
1730	case(6)
1731	SFRWL(1)=1.28432794E-03
1732	SFRWL(2)=0.12304168
1733	SFRWL(3)=0.33106142
1734	SFRWL(4)=0.32870591
1735	SFRWL(5)=0.18568763
1736	SFRWL(6)=3.02191470E-02
1737	end select
1738
1739
1740	!albedo SB <<<
1741
1742	OPEN(99,file='beta_crf.data',status='old', &
1743	form='formatted',err=9999)
1744	READ(99,*,end=9998) pfree
1745	READ(99,*,end=9998) beta_pbl
1746	READ(99,*,end=9998) beta_free
1747	READ(99,*,end=9998) lon1_beta
1748	READ(99,*,end=9998) lon2_beta
1749	READ(99,*,end=9998) lat1_beta
1750	READ(99,*,end=9998) lat2_beta
1751	READ(99,*,end=9998) mskocean_beta
1752	9998 Continue
1753	CLOSE(99)
1754	9999 Continue
1755	WRITE(,)'pfree=',pfree
1756	WRITE(,)'beta_pbl=',beta_pbl
1757	WRITE(,)'beta_free=',beta_free
1758	WRITE(,)'lon1_beta=',lon1_beta
1759	WRITE(,)'lon2_beta=',lon2_beta
1760	WRITE(,)'lat1_beta=',lat1_beta
1761	WRITE(,)'lat2_beta=',lat2_beta
1762	WRITE(,)'mskocean_beta=',mskocean_beta
1763	ENDIF
1764	!
1765	! ************** Fin de IF ( debut ) *************
1766	!
1767	!
1768	! Incrementer le compteur de la physique
1769	!
1770	itap = itap + 1
1771	IF (is_master .OR. prt_level > 9) THEN
1772	IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
1773	WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
1774	WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
1775	100 FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
1776	ENDIF
1777	ENDIF
1778	!
1779	!
1780	! Update fraction of the sub-surfaces (pctsrf) and
1781	! initialize, where a new fraction has appeared, all variables depending
1782	! on the surface fraction.
1783	!
1784	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
1785	pctsrf, fevap, z0m, z0h, agesno, &
1786	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
1787
1788	! Update time and other variables in Reprobus
1789	IF (type_trac == 'repr') THEN
1790	#ifdef REPROBUS
1791	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1792	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1793	CALL Rtime(debut)
1794	#endif
1795	ENDIF
1796
1797
1798	! Tendances bidons pour les processus qui n'affectent pas certaines
1799	! variables.
1800	du0(:,:)=0.
1801	dv0(:,:)=0.
1802	dt0 = 0.
1803	dq0(:,:)=0.
1804	dql0(:,:)=0.
1805	dqi0(:,:)=0.
1806	dsig0(:) = 0.
1807	ddens0(:) = 0.
1808	wkoccur1(:)=1
1809	!
1810	! Mettre a zero des variables de sortie (pour securite)
1811	!
1812	DO i = 1, klon
1813	d_ps(i) = 0.0
1814	ENDDO
1815	DO k = 1, klev
1816	DO i = 1, klon
1817	d_t(i,k) = 0.0
1818	d_u(i,k) = 0.0
1819	d_v(i,k) = 0.0
1820	ENDDO
1821	ENDDO
1822	DO iq = 1, nqtot
1823	DO k = 1, klev
1824	DO i = 1, klon
1825	d_qx(i,k,iq) = 0.0
1826	ENDDO
1827	ENDDO
1828	ENDDO
1829	da(:,:)=0.
1830	mp(:,:)=0.
1831	phi(:,:,:)=0.
1832	! RomP >>>
1833	phi2(:,:,:)=0.
1834	beta_prec_fisrt(:,:)=0.
1835	beta_prec(:,:)=0.
1836	epmlmMm(:,:,:)=0.
1837	eplaMm(:,:)=0.
1838	d1a(:,:)=0.
1839	dam(:,:)=0.
1840	pmflxr=0.
1841	pmflxs=0.
1842	! RomP <<<
1843
1844	!
1845	! Ne pas affecter les valeurs entrees de u, v, h, et q
1846	!
1847	DO k = 1, klev
1848	DO i = 1, klon
1849	t_seri(i,k) = t(i,k)
1850	u_seri(i,k) = u(i,k)
1851	v_seri(i,k) = v(i,k)
1852	q_seri(i,k) = qx(i,k,ivap)
1853	ql_seri(i,k) = qx(i,k,iliq)
1854	!CR: ATTENTION, on rajoute la variable glace
1855	IF (nqo.eq.2) THEN
1856	qs_seri(i,k) = 0.
1857	ELSE IF (nqo.eq.3) THEN
1858	qs_seri(i,k) = qx(i,k,isol)
1859	ENDIF
1860	ENDDO
1861	ENDDO
1862	!
1863	!--OB mass fixer
1864	IF (mass_fixer) THEN
1865	!--store initial water burden
1866	qql1(:)=0.0
1867	DO k = 1, klev
1868	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
1869	ENDDO
1870	ENDIF
1871	!--fin mass fixer
1872
1873	tke0(:,:)=pbl_tke(:,:,is_ave)
1874	!CR:Nombre de traceurs de l'eau: nqo
1875	! IF (nqtot.GE.3) THEN
1876	IF (nqtot.GE.(nqo+1)) THEN
1877	! DO iq = 3, nqtot
1878	DO iq = nqo+1, nqtot
1879	DO k = 1, klev
1880	DO i = 1, klon
1881	! tr_seri(i,k,iq-2) = qx(i,k,iq)
1882	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
1883	ENDDO
1884	ENDDO
1885	ENDDO
1886	ELSE
1887	DO k = 1, klev
1888	DO i = 1, klon
1889	tr_seri(i,k,1) = 0.0
1890	ENDDO
1891	ENDDO
1892	ENDIF
1893	!
1894	DO i = 1, klon
1895	ztsol(i) = 0.
1896	ENDDO
1897	DO nsrf = 1, nbsrf
1898	DO i = 1, klon
1899	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1900	ENDDO
1901	ENDDO
1902	! Initialize variables used for diagnostic purpose
1903	IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
1904
1905	! Diagnostiquer la tendance dynamique
1906	!
1907	IF (ancien_ok) THEN
1908	!
1909	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/dtime
1910	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/dtime
1911	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/dtime
1912	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/dtime
1913	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/dtime
1914	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/dtime
1915	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
1916	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/dtime
1917	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
1918	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/dtime
1919	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
1920	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/dtime
1921	! !! RomP >>> td dyn traceur
1922	IF (nqtot.GT.nqo) THEN ! jyg
1923	DO iq = nqo+1, nqtot ! jyg
1924	d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/dtime ! jyg
1925	ENDDO
1926	ENDIF
1927	! !! RomP <<<
1928	ELSE
1929	d_u_dyn(:,:) = 0.0
1930	d_v_dyn(:,:) = 0.0
1931	d_t_dyn(:,:) = 0.0
1932	d_q_dyn(:,:) = 0.0
1933	d_ql_dyn(:,:) = 0.0
1934	d_qs_dyn(:,:) = 0.0
1935	d_q_dyn2d(:) = 0.0
1936	d_ql_dyn2d(:) = 0.0
1937	d_qs_dyn2d(:) = 0.0
1938	! !! RomP >>> td dyn traceur
1939	IF (nqtot.GT.nqo) THEN ! jyg
1940	DO iq = nqo+1, nqtot ! jyg
1941	d_tr_dyn(:,:,iq-nqo)= 0.0 ! jyg
1942	ENDDO
1943	ENDIF
1944	! !! RomP <<<
1945	ancien_ok = .TRUE.
1946	ENDIF
1947	!
1948	! Ajouter le geopotentiel du sol:
1949	!
1950	DO k = 1, klev
1951	DO i = 1, klon
1952	zphi(i,k) = pphi(i,k) + pphis(i)
1953	ENDDO
1954	ENDDO
1955	!
1956	! Verifier les temperatures
1957	!
1958	!IM BEG
1959	IF (check) THEN
1960	amn=MIN(ftsol(1,is_ter),1000.)
1961	amx=MAX(ftsol(1,is_ter),-1000.)
1962	DO i=2, klon
1963	amn=MIN(ftsol(i,is_ter),amn)
1964	amx=MAX(ftsol(i,is_ter),amx)
1965	ENDDO
1966	!
1967	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1968	ENDIF !(check) THEN
1969	!IM END
1970	!
1971	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
1972	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
1973
1974	!
1975	!IM BEG
1976	IF (check) THEN
1977	amn=MIN(ftsol(1,is_ter),1000.)
1978	amx=MAX(ftsol(1,is_ter),-1000.)
1979	DO i=2, klon
1980	amn=MIN(ftsol(i,is_ter),amn)
1981	amx=MAX(ftsol(i,is_ter),amx)
1982	ENDDO
1983	!
1984	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1985	ENDIF !(check) THEN
1986	!IM END
1987	!
1988	! Mettre en action les conditions aux limites (albedo, sst, etc.).
1989	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
1990	!
1991	! Update ozone if day change
1992	IF (MOD(itap-1,lmt_pas) == 0) THEN
1993	IF (read_climoz <= 0) THEN
1994	! Once per day, update ozone from Royer:
1995	IF (solarlong0<-999.) then
1996	! Generic case with evolvoing season
1997	zzz=real(days_elapsed+1)
1998	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
1999	! Particular case with annual mean insolation
2000	zzz=real(90) ! could be revisited
2001	IF (read_climoz/=-1) THEN
2002	abort_message ='read_climoz=-1 is recommended when ' &
2003	// 'solarlong0=1000.'
2004	CALL abort_physic (modname,abort_message,1)
2005	ENDIF
2006	ELSE
2007	! Case where the season is imposed with solarlong0
2008	zzz=real(90) ! could be revisited
2009	ENDIF
2010
2011	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
2012	ELSE
2013	!--- ro3i = elapsed days number since current year 1st january, 0h
2014	ro3i=days_elapsed+jh_cur-jh_1jan
2015	!--- scaling for old style files (360 records)
2016	IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
2017	IF(adjust_tropopause) THEN
2018	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2019	ro3i, press_edg_climoz, paprs, wo, time_climoz, &
2020	longitude_deg, latitude_deg, press_cen_climoz, &
2021	dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
2022	ELSE
2023	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2024	ro3i, press_edg_climoz, paprs, wo, time_climoz)
2025	END IF
2026	! Convert from mole fraction of ozone to column density of ozone in a
2027	! cell, in kDU:
2028	FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2029	* zmasse / dobson_u / 1e3
2030	! (By regridding ozone values for LMDZ only once a day, we
2031	! have already neglected the variation of pressure in one
2032	! day. So do not recompute "wo" at each time step even if
2033	! "zmasse" changes a little.)
2034	ENDIF
2035	ENDIF
2036	!
2037	! Re-evaporer l'eau liquide nuageuse
2038	!
2039	CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2040	& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
2041
2042	CALL add_phys_tend &
2043	(du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
2044	'eva',abortphy,flag_inhib_tend,itap,0)
2045	call prt_enerbil('eva',itap)
2046
2047	!=========================================================================
2048	! Calculs de l'orbite.
2049	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2050	! doit donc etre plac\'e avant radlwsw et pbl_surface
2051
2052	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2053	CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2054	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2055	!
2056	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2057	! solarlong0
2058	IF (solarlong0<-999.) THEN
2059	IF (new_orbit) THEN
2060	! calcul selon la routine utilisee pour les planetes
2061	CALL solarlong(day_since_equinox, zlongi, dist)
2062	ELSE
2063	! calcul selon la routine utilisee pour l'AR4
2064	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2065	ENDIF
2066	ELSE
2067	zlongi=solarlong0 ! longitude solaire vraie
2068	dist=1. ! distance au soleil / moyenne
2069	ENDIF
2070
2071	IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2072
2073
2074	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2075	! Calcul de l'ensoleillement :
2076	! ============================
2077	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2078	! l'annee a partir d'une formule analytique.
2079	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2080	! non nul aux poles.
2081	IF (abs(solarlong0-1000.)<1.e-4) THEN
2082	CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
2083	latitude_deg,longitude_deg,rmu0,fract)
2084	swradcorr(:) = 1.0
2085	JrNt(:) = 1.0
2086	zrmu0(:) = rmu0(:)
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