Context Navigation

physiq_mod.F90 @ 2651

Last change on this file since 2651 was 2651, checked in by musat, 8 years ago
Supression fichier paramLMDZ_py.nc si pas de XIOS.
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: 156.4 KB

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

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

Download in other formats:

Original Format