Context Navigation

physiq_mod.F90 @ 2999

Last change on this file since 2999 was 2999, checked in by Ehouarn Millour, 7 years ago
Follow-up to rev 2997; in the XIOS case "missing_val" should also be set before being used... Thus moved its initialization to the beginning of physiq. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 162.0 KB

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

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

Download in other formats:

Original Format