Context Navigation

physiq_mod.F90 @ 2824

Last change on this file since 2824 was 2824, checked in by Laurent Fairhead, 7 years ago
Introduction of convoccur variable and modification of the wbeff variable in the convection so that wbeff/convoccur is the mean value of wbeff when convection occurs Jyg
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.2 KB

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

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

Download in other formats:

Original Format