Context Navigation

physiq_mod.F90 @ 2878

Last change on this file since 2878 was 2877, checked in by musat, 7 years ago
Bug correction linked to smaller call frequency of convection with respect to the physic'
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: 157.0 KB

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

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

Download in other formats:

Original Format