Context Navigation

physiq_mod.F90 @ 2629

Last change on this file since 2629 was 2629, checked in by musat, 8 years ago
Remove sequential daily TS budgets' files ini_histday_seri.h/write_histday_seri.h
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.6 KB

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

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

Download in other formats:

Original Format