Context Navigation

physiq_mod.F90 @ 2627

Last change on this file since 2627 was 2623, checked in by Ehouarn Millour, 8 years ago
Adding some non-critical initializations in physiq (in practice, mass_solu_aero and mass_solu_aero_pi are in some cases simply not used, and when used, properly initialized elsewhere). EM
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: 160.7 KB

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

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

Download in other formats:

Original Format