Context Navigation

physiq_mod.F90 @ 3029

Last change on this file since 3029 was 3029, checked in by Laurent Fairhead, 7 years ago
Change of location for xios missing value definition as the previous location did not work for the coupled model
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: 163.0 KB

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

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

Download in other formats:

Original Format