Context Navigation

physiq_mod.F90 @ 3460

Last change on this file since 3460 was 3460, checked in by Laurent Fairhead, 5 years ago
Need initialisation of RRTM before starting to write outputs
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: 177.1 KB

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

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

Download in other formats:

Original Format