Context Navigation

physiq_mod.F90 @ 5198

Last change on this file since 5198 was 2906, checked in by acozic, 7 years ago
Add an argument to initialisation of INCA
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: 161.0 KB

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

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

Download in other formats:

Original Format