Context Navigation

physiq_mod.F90 @ 3011

Last change on this file since 3011 was 3011, checked in by Laurent Fairhead, 7 years ago
Reading the number of PFTs from the orchidee.def file JG/LF
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 162.9 KB

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

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

Download in other formats:

Original Format