Context Navigation

physiq_mod.F90 @ 3462

Last change on this file since 3462 was 3462, checked in by fhourdin, 5 years ago
Changing the call to iophy FH
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: 178.1 KB

Line
1	!
2	! $Id: physiq_mod.F90 3462 2019-02-02 00:57:31Z fhourdin $
3	!
4	!#define IO_DEBUG
5	MODULE physiq_mod
6
7	IMPLICIT NONE
8
9	CONTAINS
10
11	SUBROUTINE physiq (nlon,nlev, &
12	debut,lafin,pdtphys_, &
13	paprs,pplay,pphi,pphis,presnivs, &
14	u,v,rot,t,qx, &
15	flxmass_w, &
16	d_u, d_v, d_t, d_qx, d_ps)
17
18	USE assert_m, only: assert
19	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
20	histwrite, ju2ymds, ymds2ju, getin
21	USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
22	USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
23	year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour
24	USE write_field_phy
25	USE dimphy
26	USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac
27	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo, 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, SAVE :: itap ! compteur pour la physique
518	!$OMP THREADPRIVATE(itap)
519
520	INTEGER, SAVE :: abortphy=0 ! Reprere si on doit arreter en fin de phys
521	!$OMP THREADPRIVATE(abortphy)
522	!
523	REAL,SAVE :: solarlong0
524	!$OMP THREADPRIVATE(solarlong0)
525
526	!
527	! Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
528	!
529	!IM 141004 REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
530	REAL zulow(klon),zvlow(klon)
531	!
532	INTEGER igwd,idx(klon),itest(klon)
533	!
534	! REAL,allocatable,save :: run_off_lic_0(:)
535	! !$OMP THREADPRIVATE(run_off_lic_0)
536	!ym SAVE run_off_lic_0
537	!KE43
538	! Variables liees a la convection de K. Emanuel (sb):
539	!
540	REAL, SAVE :: bas, top ! cloud base and top levels
541	!$OMP THREADPRIVATE(bas, top)
542	!------------------------------------------------------------------
543	! Upmost level reached by deep convection and related variable (jyg)
544	!
545	INTEGER izero
546	INTEGER k_upper_cv
547	!------------------------------------------------------------------
548	! Compteur de l'occurence de cvpas=1
549	INTEGER Ncvpaseq1
550	SAVE Ncvpaseq1
551	!$OMP THREADPRIVATE(Ncvpaseq1)
552	!
553	!==========================================================================
554	!CR04.12.07: on ajoute les nouvelles variables du nouveau schema
555	!de convection avec poches froides
556	! Variables li\'ees \`a la poche froide (jyg)
557
558	REAL mip(klon,klev) ! mass flux shed by the adiab ascent at each level
559	!
560	REAL wape_prescr, fip_prescr
561	INTEGER it_wape_prescr
562	SAVE wape_prescr, fip_prescr, it_wape_prescr
563	!$OMP THREADPRIVATE(wape_prescr, fip_prescr, it_wape_prescr)
564	!
565	! variables supplementaires de concvl
566	REAL Tconv(klon,klev)
567	!! variable moved to phys_local_var_mod
568	!! REAL sij(klon,klev,klev)
569	!! !
570	!! ! variables pour tester la conservation de l'energie dans concvl
571	!! REAL, DIMENSION(klon,klev) :: d_t_con_sat
572	!! REAL, DIMENSION(klon,klev) :: d_q_con_sat
573	!! REAL, DIMENSION(klon,klev) :: dql_sat
574
575	REAL, SAVE :: alp_bl_prescr=0.
576	REAL, SAVE :: ale_bl_prescr=0.
577	REAL, SAVE :: wake_s_min_lsp=0.1
578	!$OMP THREADPRIVATE(alp_bl_prescr,ale_bl_prescr)
579	!$OMP THREADPRIVATE(wake_s_min_lsp)
580
581	REAL ok_wk_lsp(klon)
582
583	!RC
584	! Variables li\'ees \`a la poche froide (jyg et rr)
585
586	INTEGER, SAVE :: iflag_wake_tend ! wake: if =0, then wake state variables are
587	! updated within calwake
588	!$OMP THREADPRIVATE(iflag_wake_tend)
589	INTEGER, SAVE :: iflag_alp_wk_cond=0 ! wake: if =0, then Alp_wk is the average lifting
590	! power provided by the wakes; else, Alp_wk is the
591	! lifting power conditionned on the presence of a
592	! gust-front in the grid cell.
593	!$OMP THREADPRIVATE(iflag_alp_wk_cond)
594	REAL t_w(klon,klev),q_w(klon,klev) ! temperature and moisture profiles in the wake region
595	REAL t_x(klon,klev),q_x(klon,klev) ! temperature and moisture profiles in the off-wake region
596
597	REAL wake_dth(klon,klev) ! wake : temp pot difference
598
599	REAL wake_omgbdth(klon,klev) ! Wake : flux of Delta_Theta
600	! transported by LS omega
601	REAL wake_dp_omgb(klon,klev) ! Wake : vertical gradient of
602	! large scale omega
603	REAL wake_dtKE(klon,klev) ! Wake : differential heating
604	! (wake - unpertubed) CONV
605	REAL wake_dqKE(klon,klev) ! Wake : differential moistening
606	! (wake - unpertubed) CONV
607	REAL wake_dp_deltomg(klon,klev) ! Wake : gradient vertical de wake_omg
608	REAL wake_spread(klon,klev) ! spreading term in wake_delt
609	!
610	!pourquoi y'a pas de save??
611	!
612	!!! INTEGER, SAVE, DIMENSION(klon) :: wake_k
613	!!! !$OMP THREADPRIVATE(wake_k)
614	!
615	!jyg<
616	!cc REAL wake_pe(klon) ! Wake potential energy - WAPE
617	!>jyg
618
619	REAL wake_fip_0(klon) ! Average Front Incoming Power (unconditionned)
620	REAL wake_gfl(klon) ! Gust Front Length
621	!!! REAL wake_dens(klon) ! moved to phys_state_var_mod
622	!
623	!
624	REAL dt_dwn(klon,klev)
625	REAL dq_dwn(klon,klev)
626	REAL M_dwn(klon,klev)
627	REAL M_up(klon,klev)
628	REAL dt_a(klon,klev)
629	REAL dq_a(klon,klev)
630	REAL d_t_adjwk(klon,klev) !jyg
631	REAL d_q_adjwk(klon,klev) !jyg
632	LOGICAL,SAVE :: ok_adjwk=.FALSE.
633	!$OMP THREADPRIVATE(ok_adjwk)
634	INTEGER,SAVE :: iflag_adjwk=0 !jyg
635	!$OMP THREADPRIVATE(iflag_adjwk) !jyg
636	REAL,SAVE :: oliqmax=999.,oicemax=999.
637	!$OMP THREADPRIVATE(oliqmax,oicemax)
638	REAL, SAVE :: alp_offset
639	!$OMP THREADPRIVATE(alp_offset)
640	REAL, SAVE :: dtcon_multistep_max=1.e6
641	!$OMP THREADPRIVATE(dtcon_multistep_max)
642	REAL, SAVE :: dqcon_multistep_max=1.e6
643	!$OMP THREADPRIVATE(dqcon_multistep_max)
644
645
646	!
647	!RR:fin declarations poches froides
648	!==========================================================================
649
650	REAL ztv(klon,klev),ztva(klon,klev)
651	REAL zpspsk(klon,klev)
652	REAL ztla(klon,klev),zqla(klon,klev)
653	REAL zthl(klon,klev)
654
655	!cc nrlmd le 10/04/2012
656
657	!--------Stochastic Boundary Layer Triggering: ALE_BL--------
658	!---Propri\'et\'es du thermiques au LCL
659	real zlcl_th(klon) ! Altitude du LCL calcul\'e
660	! continument (pcon dans
661	! thermcell_main.F90)
662	real fraca0(klon) ! Fraction des thermiques au LCL
663	real w0(klon) ! Vitesse des thermiques au LCL
664	real w_conv(klon) ! Vitesse verticale de grande \'echelle au LCL
665	real tke0(klon,klev+1) ! TKE au d\'ebut du pas de temps
666	real therm_tke_max0(klon) ! TKE dans les thermiques au LCL
667	real env_tke_max0(klon) ! TKE dans l'environnement au LCL
668
669	!JLD !---D\'eclenchement stochastique
670	!JLD integer :: tau_trig(klon)
671
672	REAL,SAVE :: random_notrig_max=1.
673	!$OMP THREADPRIVATE(random_notrig_max)
674
675	!--------Statistical Boundary Layer Closure: ALP_BL--------
676	!---Profils de TKE dans et hors du thermique
677	real therm_tke_max(klon,klev) ! Profil de TKE dans les thermiques
678	real env_tke_max(klon,klev) ! Profil de TKE dans l'environnement
679
680	!-------Activer les tendances de TKE due a l'orograp??ie---------
681	INTEGER, SAVE :: addtkeoro
682	!$OMP THREADPRIVATE(addtkeoro)
683	REAL, SAVE :: alphatkeoro
684	!$OMP THREADPRIVATE(alphatkeoro)
685	LOGICAL, SAVE :: smallscales_tkeoro
686	!$OMP THREADPRIVATE(smallscales_tkeoro)
687
688
689
690	!cc fin nrlmd le 10/04/2012
691
692	! Variables locales pour la couche limite (al1):
693	!
694	!Al1 REAL pblh(klon) ! Hauteur de couche limite
695	!Al1 SAVE pblh
696	!34EK
697	!
698	! Variables locales:
699	!
700	!AA
701	!AA Pour phytrac
702	REAL u1(klon) ! vents dans la premiere couche U
703	REAL v1(klon) ! vents dans la premiere couche V
704
705	!@$$ LOGICAL offline ! Controle du stockage ds "physique"
706	!@$$ PARAMETER (offline=.false.)
707	!@$$ INTEGER physid
708	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
709	REAL frac_nucl(klon,klev) ! idem (nucleation)
710	! RomP >>>
711	REAL beta_prec_fisrt(klon,klev) ! taux de conv de l'eau cond (fisrt)
712	! RomP <<<
713	REAL :: calday
714
715	!IM cf FH pour Tiedtke 080604
716	REAL rain_tiedtke(klon),snow_tiedtke(klon)
717	!
718	!IM 050204 END
719	REAL devap(klon) ! evaporation et sa derivee
720	REAL dsens(klon) ! chaleur sensible et sa derivee
721
722	!
723	! Conditions aux limites
724	!
725	!
726	REAL :: day_since_equinox
727	! Date de l'equinoxe de printemps
728	INTEGER, parameter :: mth_eq=3, day_eq=21
729	REAL :: jD_eq
730
731	LOGICAL, parameter :: new_orbit = .TRUE.
732
733	!
734	INTEGER lmt_pas
735	SAVE lmt_pas ! frequence de mise a jour
736	!$OMP THREADPRIVATE(lmt_pas)
737	real zmasse(klon, nbp_lev),exner(klon, nbp_lev)
738	! (column-density of mass of air in a cell, in kg m-2)
739	real, parameter:: dobson_u = 2.1415e-05 ! Dobson unit, in kg m-2
740
741	!IM sorties
742	REAL un_jour
743	PARAMETER(un_jour=86400.)
744	INTEGER itapm1 !pas de temps de la physique du(es) mois precedents
745	SAVE itapm1 !mis a jour le dernier pas de temps du mois en cours
746	!$OMP THREADPRIVATE(itapm1)
747	!======================================================================
748	!
749	! Declaration des procedures appelees
750	!
751	EXTERNAL angle ! calculer angle zenithal du soleil
752	EXTERNAL alboc ! calculer l'albedo sur ocean
753	EXTERNAL ajsec ! ajustement sec
754	EXTERNAL conlmd ! convection (schema LMD)
755	!KE43
756	EXTERNAL conema3 ! convect4.3
757	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
758	!AA
759	! JBM (3/14) fisrtilp_tr not loaded
760	! EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
761	! ! stockage des coefficients necessaires au
762	! ! lessivage OFF-LINE et ON-LINE
763	EXTERNAL hgardfou ! verifier les temperatures
764	EXTERNAL nuage ! calculer les proprietes radiatives
765	!C EXTERNAL o3cm ! initialiser l'ozone
766	EXTERNAL orbite ! calculer l'orbite terrestre
767	EXTERNAL phyetat0 ! lire l'etat initial de la physique
768	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
769	EXTERNAL suphel ! initialiser certaines constantes
770	EXTERNAL transp ! transport total de l'eau et de l'energie
771	!IM
772	EXTERNAL haut2bas !variables de haut en bas
773	EXTERNAL ini_undefSTD !initialise a 0 une variable a 1 niveau de pression
774	EXTERNAL undefSTD !somme les valeurs definies d'1 var a 1 niveau de pression
775	! EXTERNAL moy_undefSTD !moyenne d'1 var a 1 niveau de pression
776	! EXTERNAL moyglo_aire
777	! moyenne globale d'1 var ponderee par l'aire de la maille (moyglo_pondaire)
778	! par la masse/airetot (moyglo_pondaima) et la vraie masse (moyglo_pondmass)
779	!
780	!
781	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
782	! Local variables
783	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
784	!
785	REAL rhcl(klon,klev) ! humiditi relative ciel clair
786	REAL dialiq(klon,klev) ! eau liquide nuageuse
787	REAL diafra(klon,klev) ! fraction nuageuse
788	REAL cldliq(klon,klev) ! eau liquide nuageuse
789	!
790	!XXX PB
791	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
792	!
793	REAL zxfluxt(klon, klev)
794	REAL zxfluxq(klon, klev)
795	REAL zxfluxu(klon, klev)
796	REAL zxfluxv(klon, klev)
797
798	! Le rayonnement n'est pas calcule tous les pas, il faut donc
799	! sauvegarder les sorties du rayonnement
800	!ym SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
801	!ym SAVE sollwdownclr, toplwdown, toplwdownclr
802	!ym SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
803	!
804	INTEGER itaprad
805	SAVE itaprad
806	!$OMP THREADPRIVATE(itaprad)
807	!
808	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
809	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
810	!
811	#ifdef INCA
812	REAL zxsnow_dummy(klon)
813	#endif
814	REAL zsav_tsol(klon)
815	!
816	REAL dist, rmu0(klon), fract(klon)
817	REAL zrmu0(klon), zfract(klon)
818	REAL zdtime, zdtime1, zdtime2, zlongi
819	!
820	REAL qcheck
821	REAL z_avant(klon), z_apres(klon), z_factor(klon)
822	LOGICAL zx_ajustq
823	!
824	REAL za
825	REAL zx_t, zx_qs, zdelta, zcor
826	real zqsat(klon,klev)
827	!
828	INTEGER i, k, iq, j, nsrf, ll, l
829	!
830	REAL t_coup
831	PARAMETER (t_coup=234.0)
832
833	!ym A voir plus tard !!
834	!ym REAL zx_relief(iim,jjmp1)
835	!ym REAL zx_aire(iim,jjmp1)
836	!
837	! Grandeurs de sorties
838	REAL s_capCL(klon)
839	REAL s_oliqCL(klon), s_cteiCL(klon)
840	REAL s_trmb1(klon), s_trmb2(klon)
841	REAL s_trmb3(klon)
842
843	! La convection n'est pas calculee tous les pas, il faut donc
844	! sauvegarder les sorties de la convection
845	!ym SAVE
846	!ym SAVE
847	!ym SAVE
848	!
849	INTEGER itapcv, itapwk
850	SAVE itapcv, itapwk
851	!$OMP THREADPRIVATE(itapcv, itapwk)
852
853	!KE43
854	! Variables locales pour la convection de K. Emanuel (sb):
855
856	REAL tvp(klon,klev) ! virtual temp of lifted parcel
857	CHARACTER*40 capemaxcels !max(CAPE)
858
859	REAL rflag(klon) ! flag fonctionnement de convect
860	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
861
862	! -- convect43:
863	INTEGER ntra ! nb traceurs pour convect4.3
864	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
865	REAL dplcldt(klon), dplcldr(klon)
866	!? . condm_con(klon,klev),conda_con(klon,klev),
867	!? . mr_con(klon,klev),ep_con(klon,klev)
868	!? . ,sadiab(klon,klev),wadiab(klon,klev)
869	! --
870	!34EK
871	!
872	! Variables du changement
873	!
874	! con: convection
875	! lsc: condensation a grande echelle (Large-Scale-Condensation)
876	! ajs: ajustement sec
877	! eva: evaporation de l'eau liquide nuageuse
878	! vdf: couche limite (Vertical DiFfusion)
879	!
880	! tendance nulles
881	REAL, dimension(klon,klev):: du0, dv0, dt0, dq0, dql0, dqi0
882	REAL, dimension(klon) :: dsig0, ddens0
883	INTEGER, dimension(klon) :: wkoccur1
884	! tendance buffer pour appel de add_phys_tend
885	REAL, DIMENSION(klon,klev) :: d_q_ch4_dtime
886	!
887	! Flag pour pouvoir ne pas ajouter les tendances.
888	! Par defaut, les tendances doivente etre ajoutees et
889	! flag_inhib_tend = 0
890	! flag_inhib_tend > 0 : tendances non ajoutees, avec un nombre
891	! croissant de print quand la valeur du flag augmente
892	!!! attention, ce flag doit etre change avec prudence !!!
893	INTEGER :: flag_inhib_tend = 0 ! 0 is the default value
894	!! INTEGER :: flag_inhib_tend = 2
895	!
896	! Logical switch to a bug : reseting to 0 convective variables at the
897	! begining of physiq.
898	LOGICAL, SAVE :: ok_bug_cv_trac = .TRUE.
899	!$OMP THREADPRIVATE(ok_bug_cv_trac)
900	!
901	! Logical switch to a bug : changing wake_deltat when thermals are active
902	! even when there are no wakes.
903	LOGICAL, SAVE :: ok_bug_split_th = .TRUE.
904	!$OMP THREADPRIVATE(ok_bug_split_th)
905
906	!
907	!********************************************************
908	! declarations
909
910	!********************************************************
911	!IM 081204 END
912	!
913	REAL pen_u(klon,klev), pen_d(klon,klev)
914	REAL pde_u(klon,klev), pde_d(klon,klev)
915	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
916	!
917	REAL ratqsbas,ratqshaut,tau_ratqs
918	SAVE ratqsbas,ratqshaut,tau_ratqs
919	!$OMP THREADPRIVATE(ratqsbas,ratqshaut,tau_ratqs)
920	REAL, SAVE :: ratqsp0=50000., ratqsdp=20000.
921	!$OMP THREADPRIVATE(ratqsp0, ratqsdp)
922
923	! Parametres lies au nouveau schema de nuages (SB, PDF)
924	REAL, SAVE :: fact_cldcon
925	REAL, SAVE :: facttemps
926	!$OMP THREADPRIVATE(fact_cldcon,facttemps)
927	LOGICAL, SAVE :: ok_newmicro
928	!$OMP THREADPRIVATE(ok_newmicro)
929
930	INTEGER, SAVE :: iflag_cld_th
931	!$OMP THREADPRIVATE(iflag_cld_th)
932	!IM logical ptconv(klon,klev) !passe dans phys_local_var_mod
933	!IM cf. AM 081204 BEG
934	LOGICAL ptconvth(klon,klev)
935	!IM cf. AM 081204 END
936	!
937	! Variables liees a l'ecriture de la bande histoire physique
938	!
939	!======================================================================
940	!
941	!
942	!JLD integer itau_w ! pas de temps ecriture = itap + itau_phy
943	!
944	!
945	! Variables locales pour effectuer les appels en serie
946	!
947	!IM RH a 2m (la surface)
948	REAL Lheat
949
950	INTEGER length
951	PARAMETER ( length = 100 )
952	REAL tabcntr0( length )
953	!
954	!JLD INTEGER ndex2d(nbp_lon*nbp_lat)
955	!IM
956	!
957	!IM AMIP2 BEG
958	!JLD REAL moyglo, mountor
959	!IM 141004 BEG
960	REAL zustrdr(klon), zvstrdr(klon)
961	REAL zustrli(klon), zvstrli(klon)
962	REAL zustrph(klon), zvstrph(klon)
963	REAL aam, torsfc
964	!IM 141004 END
965	!IM 190504 BEG
966	! INTEGER imp1jmp1
967	! PARAMETER(imp1jmp1=(iim+1)*jjmp1)
968	!ym A voir plus tard
969	! REAL zx_tmp((nbp_lon+1)*nbp_lat)
970	! REAL airedyn(nbp_lon+1,nbp_lat)
971	!IM 190504 END
972	!JLD LOGICAL ok_msk
973	!JLD REAL msk(klon)
974	!ym A voir plus tard
975	!ym REAL zm_wo(jjmp1, klev)
976	!IM AMIP2 END
977	!
978	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
979	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
980	!JLD REAL zx_tmp_2d(nbp_lon,nbp_lat)
981	!JLD REAL zx_lon(nbp_lon,nbp_lat)
982	!JLD REAL zx_lat(nbp_lon,nbp_lat)
983	!
984	INTEGER nid_ctesGCM
985	SAVE nid_ctesGCM
986	!$OMP THREADPRIVATE(nid_ctesGCM)
987	!
988	!IM 280405 BEG
989	! INTEGER nid_bilKPins, nid_bilKPave
990	! SAVE nid_bilKPins, nid_bilKPave
991	! !$OMP THREADPRIVATE(nid_bilKPins, nid_bilKPave)
992	!
993	REAL ve_lay(klon,klev) ! transport meri. de l'energie a chaque niveau vert.
994	REAL vq_lay(klon,klev) ! transport meri. de l'eau a chaque niveau vert.
995	REAL ue_lay(klon,klev) ! transport zonal de l'energie a chaque niveau vert.
996	REAL uq_lay(klon,klev) ! transport zonal de l'eau a chaque niveau vert.
997	!
998	!JLD REAL zjulian
999	!JLD SAVE zjulian
1000	!JLD!$OMP THREADPRIVATE(zjulian)
1001
1002	!JLD INTEGER nhori, nvert
1003	!JLD REAL zsto
1004	!JLD REAL zstophy, zout
1005
1006	CHARACTER*20 modname
1007	CHARACTER*80 abort_message
1008	LOGICAL, SAVE :: ok_sync, ok_sync_omp
1009	!$OMP THREADPRIVATE(ok_sync)
1010	REAL date0
1011
1012	! essai writephys
1013	INTEGER fid_day, fid_mth, fid_ins
1014	PARAMETER (fid_ins = 1, fid_day = 2, fid_mth = 3)
1015	INTEGER prof2d_on, prof3d_on, prof2d_av, prof3d_av
1016	PARAMETER (prof2d_on = 1, prof3d_on = 2, prof2d_av = 3, prof3d_av = 4)
1017	REAL ztsol(klon)
1018	REAL q2m(klon,nbsrf) ! humidite a 2m
1019
1020	!IM: t2m, q2m, ustar, u10m, v10m et t2mincels, t2maxcels
1021	CHARACTER*40 t2mincels, t2maxcels !t2m min., t2m max
1022	CHARACTER*40 tinst, tave
1023	REAL cldtaupi(klon,klev) ! Cloud optical thickness for
1024	! pre-industrial (pi) aerosols
1025
1026	INTEGER :: naero
1027	! Aerosol optical properties
1028	CHARACTER*4, DIMENSION(naero_grp) :: rfname
1029	REAL, DIMENSION(klon,klev) :: mass_solu_aero ! total mass
1030	! concentration
1031	! for all soluble
1032	! aerosols[ug/m3]
1033	REAL, DIMENSION(klon,klev) :: mass_solu_aero_pi
1034	! - " - (pre-industrial value)
1035
1036	! Parameters
1037	LOGICAL ok_ade, ok_aie ! Apply aerosol (in)direct effects or not
1038	LOGICAL ok_alw ! Apply aerosol LW effect or not
1039	LOGICAL ok_cdnc ! ok cloud droplet number concentration (O. Boucher 01-2013)
1040	REAL bl95_b0, bl95_b1 ! Parameter in Boucher and Lohmann (1995)
1041	SAVE ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1
1042	!$OMP THREADPRIVATE(ok_ade, ok_aie, ok_alw, ok_cdnc, bl95_b0, bl95_b1)
1043	LOGICAL, SAVE :: aerosol_couple ! true : calcul des aerosols dans INCA
1044	! false : lecture des aerosol dans un fichier
1045	!$OMP THREADPRIVATE(aerosol_couple)
1046	LOGICAL, SAVE :: chemistry_couple ! true : use INCA chemistry O3
1047	! false : use offline chemistry O3
1048	!$OMP THREADPRIVATE(chemistry_couple)
1049	INTEGER, SAVE :: flag_aerosol
1050	!$OMP THREADPRIVATE(flag_aerosol)
1051	LOGICAL, SAVE :: flag_bc_internal_mixture
1052	!$OMP THREADPRIVATE(flag_bc_internal_mixture)
1053	LOGICAL, SAVE :: new_aod
1054	!$OMP THREADPRIVATE(new_aod)
1055	!
1056	!--STRAT AEROSOL
1057	INTEGER, SAVE :: flag_aerosol_strat
1058	!$OMP THREADPRIVATE(flag_aerosol_strat)
1059	!
1060	!--INTERACTIVE AEROSOL FEEDBACK ON RADIATION
1061	LOGICAL, SAVE :: flag_aer_feedback
1062	!$OMP THREADPRIVATE(flag_aer_feedback)
1063
1064	!c-fin STRAT AEROSOL
1065	!
1066	! Declaration des constantes et des fonctions thermodynamiques
1067	!
1068	LOGICAL,SAVE :: first=.TRUE.
1069	!$OMP THREADPRIVATE(first)
1070
1071	! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
1072	! Note that pressure vectors are in Pa and in stricly ascending order
1073	INTEGER,SAVE :: read_climoz ! Read ozone climatology
1074	! (let it keep the default OpenMP shared attribute)
1075	! Allowed values are 0, 1 and 2
1076	! 0: do not read an ozone climatology
1077	! 1: read a single ozone climatology that will be used day and night
1078	! 2: read two ozone climatologies, the average day and night
1079	! climatology and the daylight climatology
1080	INTEGER,SAVE :: ncid_climoz ! NetCDF file identifier
1081	REAL, POINTER, SAVE :: press_cen_climoz(:) ! Pressure levels
1082	REAL, POINTER, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
1083	REAL, POINTER, SAVE :: time_climoz(:) ! Time vector
1084	CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
1085	= ["tro3 ","tro3_daylight"]
1086	! vars_climoz(1:read_climoz): variables names in climoz file.
1087	! vars_climoz(1:read_climoz-2) if read_climoz>2 (temporary)
1088	REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
1089	! the ozone fields, old method.
1090
1091	include "YOMCST.h"
1092	include "YOETHF.h"
1093	include "FCTTRE.h"
1094	!IM 100106 BEG : pouvoir sortir les ctes de la physique
1095	include "conema3.h"
1096	include "fisrtilp.h"
1097	include "nuage.h"
1098	include "compbl.h"
1099	!IM 100106 END : pouvoir sortir les ctes de la physique
1100	!
1101	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1102	! Declarations pour Simulateur COSP
1103	!============================================================
1104	real :: mr_ozone(klon,klev)
1105
1106	!IM stations CFMIP
1107	INTEGER, SAVE :: nCFMIP
1108	!$OMP THREADPRIVATE(nCFMIP)
1109	INTEGER, PARAMETER :: npCFMIP=120
1110	INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
1111	REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
1112	!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
1113	INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
1114	REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
1115	!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
1116	INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
1117	!$OMP THREADPRIVATE(iGCM, jGCM)
1118	logical, dimension(nfiles) :: phys_out_filestations
1119	logical, parameter :: lNMC=.FALSE.
1120
1121	!IM betaCRF
1122	REAL, SAVE :: pfree, beta_pbl, beta_free
1123	!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
1124	REAL, SAVE :: lon1_beta, lon2_beta, lat1_beta, lat2_beta
1125	!$OMP THREADPRIVATE(lon1_beta, lon2_beta, lat1_beta, lat2_beta)
1126	LOGICAL, SAVE :: mskocean_beta
1127	!$OMP THREADPRIVATE(mskocean_beta)
1128	REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et
1129	! cldemirad pour evaluer les
1130	! retros liees aux CRF
1131	REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique
1132	! pour radlwsw pour
1133	! tester "CRF off"
1134	REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique
1135	! pour radlwsw pour
1136	! tester "CRF off"
1137	REAL, dimension(klon, klev) :: cldemirad ! emissivite pour
1138	! radlwsw pour tester
1139	! "CRF off"
1140	REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse
1141
1142	#ifdef INCA
1143	! set de variables utilisees pour l'initialisation des valeurs provenant de INCA
1144	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_tauinca
1145	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_pizinca
1146	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_cginca
1147	REAL, DIMENSION(klon,klev,nbands) :: init_ccminca
1148	#endif
1149	REAL, DIMENSION(klon,nbtr) :: init_source
1150
1151	!lwoff=y : offset LW CRE for radiation code and other schemes
1152	REAL, SAVE :: betalwoff
1153	!OMP THREADPRIVATE(betalwoff)
1154	!
1155	INTEGER :: nbtr_tmp ! Number of tracer inside concvl
1156	REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
1157	REAL, dimension(klon,klev) :: ch_in ! Condensed humidity entering in phytrac (eau liquide)
1158	integer iostat
1159
1160	REAL zzz
1161	!albedo SB >>>
1162	REAL,DIMENSION(6), SAVE :: SFRWL
1163	!$OMP THREADPRIVATE(SFRWL)
1164	!albedo SB <<<
1165
1166	!--OB variables for mass fixer (hard coded for now)
1167	LOGICAL, PARAMETER :: mass_fixer=.FALSE.
1168	REAL qql1(klon),qql2(klon),corrqql
1169
1170	REAL pi
1171
1172	pi = 4. * ATAN(1.)
1173
1174	! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
1175	jjmp1=nbp_lat
1176
1177	!======================================================================
1178	! Gestion calendrier : mise a jour du module phys_cal_mod
1179	!
1180	pdtphys=pdtphys_
1181	CALL update_time(pdtphys)
1182	phys_tstep=NINT(pdtphys)
1183	#ifdef CPP_XIOS
1184	IF (.NOT. debut .AND. is_omp_master) CALL xios_update_calendar(itap+1)
1185	#endif
1186
1187	!======================================================================
1188	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1189	! Utilise notamment en 1D mais peut etre active egalement en 3D
1190	! en imposant la valeur de igout.
1191	!======================================================================d
1192	IF (prt_level.ge.1) THEN
1193	igout=klon/2+1/klon
1194	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1195	write(lunout,*) 'igout, lat, lon ',igout, latitude_deg(igout), &
1196	longitude_deg(igout)
1197	write(lunout,*) &
1198	'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1199	write(lunout,*) &
1200	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1201
1202	write(lunout,*) 'paprs, play, phi, u, v, t'
1203	DO k=1,klev
1204	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
1205	u(igout,k),v(igout,k),t(igout,k)
1206	ENDDO
1207	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1208	DO k=1,klev
1209	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1210	ENDDO
1211	ENDIF
1212
1213	! Quick check on pressure levels:
1214	CALL assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
1215	"physiq_mod paprs bad order")
1216
1217	IF (first) THEN
1218	CALL init_etat0_limit_unstruct
1219	IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
1220	!CR:nvelles variables convection/poches froides
1221
1222	WRITE(lunout,*) '================================================='
1223	WRITE(lunout,*) 'Allocation des variables locales et sauvegardees'
1224	WRITE(lunout,*) '================================================='
1225	CALL phys_local_var_init
1226	!
1227	! appel a la lecture du run.def physique
1228	CALL conf_phys(ok_journe, ok_mensuel, &
1229	ok_instan, ok_hf, &
1230	ok_LES, &
1231	callstats, &
1232	solarlong0,seuil_inversion, &
1233	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
1234	iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
1235	ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple, &
1236	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, new_aod, &
1237	flag_bc_internal_mixture, bl95_b0, bl95_b1, &
1238	! nv flags pour la convection et les
1239	! poches froides
1240	read_climoz, &
1241	alp_offset)
1242	CALL phys_state_var_init(read_climoz)
1243	CALL phys_output_var_init
1244	IF(read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
1245
1246	print*, '================================================='
1247	!
1248	!CR: check sur le nb de traceurs de l eau
1249	IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
1250	WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
1251	'(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
1252	STOP
1253	ENDIF
1254
1255	Ncvpaseq1 = 0
1256	dnwd0=0.0
1257	ftd=0.0
1258	fqd=0.0
1259	cin=0.
1260	!ym Attention pbase pas initialise dans concvl !!!!
1261	pbase=0
1262	!IM 180608
1263
1264	itau_con=0
1265	first=.FALSE.
1266
1267	ENDIF ! first
1268
1269	!ym => necessaire pour iflag_con != 2
1270	pmfd(:,:) = 0.
1271	pen_u(:,:) = 0.
1272	pen_d(:,:) = 0.
1273	pde_d(:,:) = 0.
1274	pde_u(:,:) = 0.
1275	aam=0.
1276	d_t_adjwk(:,:)=0
1277	d_q_adjwk(:,:)=0
1278
1279	alp_bl_conv(:)=0.
1280
1281	torsfc=0.
1282	forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1283
1284	modname = 'physiq'
1285
1286	IF (debut) THEN
1287	CALL suphel ! initialiser constantes et parametres phys.
1288	! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1289	tau_gl=5.
1290	CALL getin_p('tau_gl', tau_gl)
1291	! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1292	! secondes
1293	tau_gl=86400.*tau_gl
1294	WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
1295
1296	CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
1297	CALL getin_p('random_notrig_max',random_notrig_max)
1298	CALL getin_p('ok_adjwk',ok_adjwk)
1299	IF (ok_adjwk) iflag_adjwk=2 ! for compatibility with older versions
1300	! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
1301	! 1 => convective adjustment but state variables are unchanged
1302	! 2 => convective adjustment and state variables are changed
1303	CALL getin_p('iflag_adjwk',iflag_adjwk)
1304	CALL getin_p('dtcon_multistep_max',dtcon_multistep_max)
1305	CALL getin_p('dqcon_multistep_max',dqcon_multistep_max)
1306	CALL getin_p('oliqmax',oliqmax)
1307	CALL getin_p('oicemax',oicemax)
1308	CALL getin_p('ratqsp0',ratqsp0)
1309	CALL getin_p('ratqsdp',ratqsdp)
1310	iflag_wake_tend = 0
1311	CALL getin_p('iflag_wake_tend',iflag_wake_tend)
1312	ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
1313	! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
1314	CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
1315	CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
1316	CALL getin_p('ok_bug_split_th',ok_bug_split_th)
1317	fl_ebil = 0 ! by default, conservation diagnostics are desactivated
1318	CALL getin_p('fl_ebil',fl_ebil)
1319	fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
1320	CALL getin_p('fl_cor_ebil',fl_cor_ebil)
1321	iflag_phytrac = 1 ! by default we do want to call phytrac
1322	CALL getin_p('iflag_phytrac',iflag_phytrac)
1323	nvm_lmdz = 13
1324	CALL getin_p('NVM',nvm_lmdz)
1325
1326	WRITE(lunout,*) 'iflag_alp_wk_cond=', iflag_alp_wk_cond
1327	WRITE(lunout,*) 'random_ntrig_max=', random_notrig_max
1328	WRITE(lunout,*) 'ok_adjwk=', ok_adjwk
1329	WRITE(lunout,*) 'iflag_adjwk=', iflag_adjwk
1330	WRITE(lunout,*) 'qtcon_multistep_max=',dtcon_multistep_max
1331	WRITE(lunout,*) 'qdcon_multistep_max=',dqcon_multistep_max
1332	WRITE(lunout,*) 'ratqsp0=', ratqsp0
1333	WRITE(lunout,*) 'ratqsdp=', ratqsdp
1334	WRITE(lunout,*) 'iflag_wake_tend=', iflag_wake_tend
1335	WRITE(lunout,*) 'ok_bad_ecmwf_thermo=',ok_bad_ecmwf_thermo
1336	WRITE(lunout,*) 'ok_bug_cv_trac=', ok_bug_cv_trac
1337	WRITE(lunout,*) 'ok_bug_split_th=', ok_bug_split_th
1338	WRITE(lunout,*) 'fl_ebil=', fl_ebil
1339	WRITE(lunout,*) 'fl_cor_ebil=', fl_cor_ebil
1340	WRITE(lunout,*) 'iflag_phytrac=', iflag_phytrac
1341	WRITE(lunout,*) 'NVM=', nvm_lmdz
1342
1343	!--PC: defining fields to be exchanged between LMDz, ORCHIDEE and NEMO
1344	WRITE(lunout,*) 'Call to infocfields from physiq'
1345	CALL infocfields_init
1346
1347	ENDIF
1348
1349	IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
1350
1351	!======================================================================
1352	! Gestion calendrier : mise a jour du module phys_cal_mod
1353	!
1354	! CALL phys_cal_update(jD_cur,jH_cur)
1355
1356	!
1357	! Si c'est le debut, il faut initialiser plusieurs choses
1358	! ********
1359	!
1360	IF (debut) THEN
1361	!rv CRinitialisation de wght_th et lalim_conv pour la
1362	!definition de la couche alimentation de la convection a partir
1363	!des caracteristiques du thermique
1364	wght_th(:,:)=1.
1365	lalim_conv(:)=1
1366	!RC
1367	ustar(:,:)=0.
1368	! u10m(:,:)=0.
1369	! v10m(:,:)=0.
1370	rain_con(:)=0.
1371	snow_con(:)=0.
1372	topswai(:)=0.
1373	topswad(:)=0.
1374	solswai(:)=0.
1375	solswad(:)=0.
1376
1377	wmax_th(:)=0.
1378	tau_overturning_th(:)=0.
1379
1380	IF (type_trac == 'inca') THEN
1381	! jg : initialisation jusqu'au ces variables sont dans restart
1382	ccm(:,:,:) = 0.
1383	tau_aero(:,:,:,:) = 0.
1384	piz_aero(:,:,:,:) = 0.
1385	cg_aero(:,:,:,:) = 0.
1386
1387	config_inca='none' ! default
1388	CALL getin_p('config_inca',config_inca)
1389
1390	ELSE
1391	config_inca='none' ! default
1392	ENDIF
1393
1394	tau_aero(:,:,:,:) = 1.e-15
1395	piz_aero(:,:,:,:) = 1.
1396	cg_aero(:,:,:,:) = 0.
1397
1398	IF (aerosol_couple .AND. (config_inca /= "aero" &
1399	.AND. config_inca /= "aeNP ")) THEN
1400	abort_message &
1401	= 'if aerosol_couple is activated, config_inca need to be ' &
1402	// 'aero or aeNP'
1403	CALL abort_physic (modname,abort_message,1)
1404	ENDIF
1405
1406	rnebcon0(:,:) = 0.0
1407	clwcon0(:,:) = 0.0
1408	rnebcon(:,:) = 0.0
1409	clwcon(:,:) = 0.0
1410
1411	!
1412	print*,'iflag_coupl,iflag_clos,iflag_wake', &
1413	iflag_coupl,iflag_clos,iflag_wake
1414	print*,'iflag_cycle_diurne', iflag_cycle_diurne
1415	!
1416	IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
1417	abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
1418	CALL abort_physic (modname,abort_message,1)
1419	ENDIF
1420	!
1421	!
1422	! Initialiser les compteurs:
1423	!
1424	itap = 0
1425	itaprad = 0
1426	itapcv = 0
1427	itapwk = 0
1428
1429	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1430	!! Un petit travail \`a faire ici.
1431	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1432
1433	IF (iflag_pbl>1) THEN
1434	PRINT*, "Using method MELLOR&YAMADA"
1435	ENDIF
1436
1437	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1438	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
1439	! phylmd plutot que dyn3d
1440	! Attention : la version precedente n'etait pas tres propre.
1441	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1442	! pour obtenir le meme resultat.
1443	!jyg for fh<
1444	WRITE(lunout,*) 'Pas de temps phys_tstep pdtphys ',phys_tstep,pdtphys
1445	IF (abs(phys_tstep-pdtphys)>1.e-10) THEN
1446	abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
1447	CALL abort_physic(modname,abort_message,1)
1448	ENDIF
1449	!>jyg
1450	IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
1451	radpas = NINT( 86400./phys_tstep)/nbapp_rad
1452	ELSE
1453	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1454	'multiple de nbapp_rad'
1455	WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
1456	'mais 1+1<>2'
1457	abort_message='nbre de pas de temps physique n est pas multiple ' &
1458	// 'de nbapp_rad'
1459	CALL abort_physic(modname,abort_message,1)
1460	ENDIF
1461	IF (nbapp_cv .EQ. 0) nbapp_cv=86400./phys_tstep
1462	IF (nbapp_wk .EQ. 0) nbapp_wk=86400./phys_tstep
1463	print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
1464	IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
1465	cvpas_0 = NINT( 86400./phys_tstep)/nbapp_cv
1466	cvpas = cvpas_0
1467	print *,'physiq, cvpas ',cvpas
1468	ELSE
1469	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1470	'multiple de nbapp_cv'
1471	WRITE(lunout,*) 'changer nbapp_cv ou alors commenter ce test ', &
1472	'mais 1+1<>2'
1473	abort_message='nbre de pas de temps physique n est pas multiple ' &
1474	// 'de nbapp_cv'
1475	CALL abort_physic(modname,abort_message,1)
1476	ENDIF
1477	IF (MOD(NINT(86400./phys_tstep),nbapp_wk).EQ.0) THEN
1478	wkpas = NINT( 86400./phys_tstep)/nbapp_wk
1479	! print *,'physiq, wkpas ',wkpas
1480	ELSE
1481	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1482	'multiple de nbapp_wk'
1483	WRITE(lunout,*) 'changer nbapp_wk ou alors commenter ce test ', &
1484	'mais 1+1<>2'
1485	abort_message='nbre de pas de temps physique n est pas multiple ' &
1486	// 'de nbapp_wk'
1487	CALL abort_physic(modname,abort_message,1)
1488	ENDIF
1489	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1490	IF (grid_type/=unstructured) CALL init_iophy_new(latitude_deg,longitude_deg)
1491
1492	!===================================================================
1493	!IM stations CFMIP
1494	nCFMIP=npCFMIP
1495	OPEN(98,file='npCFMIP_param.data',status='old', &
1496	form='formatted',iostat=iostat)
1497	IF (iostat == 0) THEN
1498	READ(98,*,end=998) nCFMIP
1499	998 CONTINUE
1500	CLOSE(98)
1501	CONTINUE
1502	IF(nCFMIP.GT.npCFMIP) THEN
1503	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1504	CALL abort_physic("physiq", "", 1)
1505	ELSE
1506	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1507	ENDIF
1508
1509	!
1510	ALLOCATE(tabCFMIP(nCFMIP))
1511	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1512	ALLOCATE(tabijGCM(nCFMIP))
1513	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1514	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1515	!
1516	! lecture des nCFMIP stations CFMIP, de leur numero
1517	! et des coordonnees geographiques lonCFMIP, latCFMIP
1518	!
1519	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
1520	lonCFMIP, latCFMIP)
1521	!
1522	! identification des
1523	! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
1524	! grille de LMDZ
1525	! 2) indices points tabijGCM de la grille physique 1d sur
1526	! klon points
1527	! 3) indices iGCM, jGCM de la grille physique 2d
1528	!
1529	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1530	tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1531	!
1532	ELSE
1533	ALLOCATE(tabijGCM(0))
1534	ALLOCATE(lonGCM(0), latGCM(0))
1535	ALLOCATE(iGCM(0), jGCM(0))
1536	ENDIF
1537
1538	#ifdef CPP_IOIPSL
1539
1540	!$OMP MASTER
1541	! FH : if ok_sync=.true. , the time axis is written at each time step
1542	! in the output files. Only at the end in the opposite case
1543	ok_sync_omp=.FALSE.
1544	CALL getin('ok_sync',ok_sync_omp)
1545	CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
1546	iGCM,jGCM,lonGCM,latGCM, &
1547	jjmp1,nlevSTD,clevSTD,rlevSTD, phys_tstep,ok_veget, &
1548	type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1549	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1550	read_climoz, phys_out_filestations, &
1551	new_aod, aerosol_couple, &
1552	flag_aerosol_strat, pdtphys, paprs, pphis, &
1553	pplay, lmax_th, ptconv, ptconvth, ivap, &
1554	d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
1555	!$OMP END MASTER
1556	!$OMP BARRIER
1557	ok_sync=ok_sync_omp
1558
1559	freq_outNMC(1) = ecrit_files(7)
1560	freq_outNMC(2) = ecrit_files(8)
1561	freq_outNMC(3) = ecrit_files(9)
1562	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1563	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1564	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1565
1566	#ifndef CPP_XIOS
1567	CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
1568	#endif
1569
1570	#endif
1571	ecrit_reg = ecrit_reg * un_jour
1572	ecrit_tra = ecrit_tra * un_jour
1573
1574	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
1575	date0 = jD_ref
1576	WRITE(,) 'physiq date0 : ',date0
1577	!
1578
1579	! CALL create_climoz(read_climoz)
1580	CALL init_aero_fromfile(flag_aerosol) !! initialise aero from file for XIOS interpolation (unstructured_grid)
1581	CALL init_readaerosolstrato(flag_aerosol_strat) !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
1582
1583	IF(read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
1584	CALL create_etat0_limit_unstruct
1585	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1586	IF (grid_type==unstructured) CALL init_iophy_new(latitude_deg,longitude_deg)
1587
1588	!jyg<
1589	IF (klon_glo==1) THEN
1590	IF (iflag_pbl > 1) THEN
1591	pbl_tke(:,:,is_ave) = 0.
1592	DO nsrf=1,nbsrf
1593	DO k = 1,klev+1
1594	pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
1595	+pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
1596	ENDDO
1597	ENDDO
1598	ELSE ! (iflag_pbl > 1)
1599	pbl_tke(:,:,:) = 0.
1600	ENDIF ! (iflag_pbl > 1)
1601	ELSE
1602	pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
1603	!>jyg
1604	ENDIF
1605	#ifdef CPP_COSP
1606
1607	IF (ok_cosp) THEN
1608	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
1609	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1610	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1611	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1612	JrNt,ref_liq,ref_ice, &
1613	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1614	zu10m,zv10m,pphis, &
1615	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1616	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1617	prfl(:,1:klev),psfl(:,1:klev), &
1618	pmflxr(:,1:klev),pmflxs(:,1:klev), &
1619	mr_ozone,cldtau, cldemi)
1620	ENDIF
1621	#endif
1622
1623	!
1624	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1625	! Nouvelle initialisation pour le rayonnement RRTM
1626	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1627
1628	CALL iniradia(klon,klev,paprs(1,1:klev+1))
1629
1630
1631	#ifdef CPP_XIOS
1632	IF (is_omp_master) CALL xios_update_calendar(1)
1633	#endif
1634
1635	!IM begin
1636	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1637	,ratqs(1,1)
1638	!IM end
1639
1640
1641	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1642	!
1643	! on remet le calendrier a zero
1644	!
1645	IF (raz_date .eq. 1) THEN
1646	itau_phy = 0
1647	ENDIF
1648
1649	! IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
1650	! WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
1651	! pdtphys
1652	! abort_message='Pas physique n est pas correct '
1653	! ! call abort_physic(modname,abort_message,1)
1654	! phys_tstep=pdtphys
1655	! ENDIF
1656	IF (nlon .NE. klon) THEN
1657	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1658	klon
1659	abort_message='nlon et klon ne sont pas coherents'
1660	CALL abort_physic(modname,abort_message,1)
1661	ENDIF
1662	IF (nlev .NE. klev) THEN
1663	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1664	klev
1665	abort_message='nlev et klev ne sont pas coherents'
1666	CALL abort_physic(modname,abort_message,1)
1667	ENDIF
1668	!
1669	IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
1670	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1671	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1672	abort_message='Nbre d appels au rayonnement insuffisant'
1673	CALL abort_physic(modname,abort_message,1)
1674	ENDIF
1675	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1676	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1677	ok_cvl
1678	!
1679	!KE43
1680	! Initialisation pour la convection de K.E. (sb):
1681	IF (iflag_con.GE.3) THEN
1682
1683	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1684	WRITE(lunout,*) &
1685	"On va utiliser le melange convectif des traceurs qui"
1686	WRITE(lunout,*)"est calcule dans convect4.3"
1687	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1688
1689	DO i = 1, klon
1690	ema_cbmf(i) = 0.
1691	ema_pcb(i) = 0.
1692	ema_pct(i) = 0.
1693	! ema_workcbmf(i) = 0.
1694	ENDDO
1695	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1696	DO i = 1, klon
1697	ibas_con(i) = 1
1698	itop_con(i) = 1
1699	ENDDO
1700	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1701	!================================================================
1702	!CR:04.12.07: initialisations poches froides
1703	! Controle de ALE et ALP pour la fermeture convective (jyg)
1704	IF (iflag_wake>=1) THEN
1705	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1706	,alp_bl_prescr, ale_bl_prescr)
1707	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1708	! print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
1709	!
1710	! Initialize tendencies of wake state variables (for some flag values
1711	! they are not computed).
1712	d_deltat_wk(:,:) = 0.
1713	d_deltaq_wk(:,:) = 0.
1714	d_deltat_wk_gw(:,:) = 0.
1715	d_deltaq_wk_gw(:,:) = 0.
1716	d_deltat_vdf(:,:) = 0.
1717	d_deltaq_vdf(:,:) = 0.
1718	d_deltat_the(:,:) = 0.
1719	d_deltaq_the(:,:) = 0.
1720	d_deltat_ajs_cv(:,:) = 0.
1721	d_deltaq_ajs_cv(:,:) = 0.
1722	d_s_wk(:) = 0.
1723	d_dens_wk(:) = 0.
1724	ENDIF
1725
1726	! do i = 1,klon
1727	! Ale_bl(i)=0.
1728	! Alp_bl(i)=0.
1729	! enddo
1730
1731	!ELSE
1732	! ALLOCATE(tabijGCM(0))
1733	! ALLOCATE(lonGCM(0), latGCM(0))
1734	! ALLOCATE(iGCM(0), jGCM(0))
1735	ENDIF
1736
1737	DO i=1,klon
1738	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1739	ENDDO
1740
1741	!34EK
1742	IF (ok_orodr) THEN
1743
1744	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1745	! FH sans doute a enlever de finitivement ou, si on le
1746	! garde, l'activer justement quand ok_orodr = false.
1747	! ce rugoro est utilise par la couche limite et fait double emploi
1748	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1749	! DO i=1,klon
1750	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1751	! ENDDO
1752	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1753	IF (ok_strato) THEN
1754	CALL SUGWD_strato(klon,klev,paprs,pplay)
1755	ELSE
1756	CALL SUGWD(klon,klev,paprs,pplay)
1757	ENDIF
1758
1759	DO i=1,klon
1760	zuthe(i)=0.
1761	zvthe(i)=0.
1762	IF (zstd(i).gt.10.) THEN
1763	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1764	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1765	ENDIF
1766	ENDDO
1767	ENDIF
1768	!
1769	!
1770	lmt_pas = NINT(86400./phys_tstep * 1.0) ! tous les jours
1771	WRITE(lunout,*)'La frequence de lecture surface est de ', &
1772	lmt_pas
1773	!
1774	capemaxcels = 't_max(X)'
1775	t2mincels = 't_min(X)'
1776	t2maxcels = 't_max(X)'
1777	tinst = 'inst(X)'
1778	tave = 'ave(X)'
1779	!IM cf. AM 081204 BEG
1780	!IM cf. AM 081204 END
1781	!
1782	!=============================================================
1783	! Initialisation des sorties
1784	!=============================================================
1785
1786	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
1787	pplay, lmax_th, aerosol_couple, &
1788	ok_ade, ok_aie, ivap, iliq, isol, new_aod, ok_sync,&
1789	ptconv, read_climoz, clevSTD, &
1790	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
1791	flag_aerosol, flag_aerosol_strat, ok_cdnc)
1792
1793	#ifdef CPP_XIOS
1794	! Get "missing_val" value from XML files (from temperature variable)
1795	!$OMP MASTER
1796	CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1797	!$OMP END MASTER
1798	!$OMP BARRIER
1799	missing_val=missing_val_omp
1800	#endif
1801
1802	#ifdef CPP_XIOS
1803	! Need to put this initialisation after phyetat0 as in the coupled model the XIOS context is only
1804	! initialised at that moment
1805	! Get "missing_val" value from XML files (from temperature variable)
1806	!$OMP MASTER
1807	CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1808	!$OMP END MASTER
1809	!$OMP BARRIER
1810	missing_val=missing_val_omp
1811	#endif
1812
1813
1814	CALL printflag( tabcntr0,radpas,ok_journe, &
1815	ok_instan, ok_region )
1816	!
1817	!
1818	!
1819	! Prescrire l'ozone dans l'atmosphere
1820	!
1821	!
1822	!c DO i = 1, klon
1823	!c DO k = 1, klev
1824	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1825	!c ENDDO
1826	!c ENDDO
1827	!
1828	IF (type_trac == 'inca') THEN
1829	#ifdef INCA
1830	CALL VTe(VTphysiq)
1831	CALL VTb(VTinca)
1832	calday = REAL(days_elapsed) + jH_cur
1833	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1834
1835	CALL chemini( &
1836	rg, &
1837	ra, &
1838	cell_area, &
1839	latitude_deg, &
1840	longitude_deg, &
1841	presnivs, &
1842	calday, &
1843	klon, &
1844	nqtot, &
1845	nqo, &
1846	pdtphys, &
1847	annee_ref, &
1848	year_cur, &
1849	day_ref, &
1850	day_ini, &
1851	start_time, &
1852	itau_phy, &
1853	date0, &
1854	io_lon, &
1855	io_lat, &
1856	chemistry_couple, &
1857	init_source, &
1858	init_tauinca, &
1859	init_pizinca, &
1860	init_cginca, &
1861	init_ccminca)
1862
1863
1864	! initialisation des variables depuis le restart de inca
1865	ccm(:,:,:) = init_ccminca
1866	tau_aero(:,:,:,:) = init_tauinca
1867	piz_aero(:,:,:,:) = init_pizinca
1868	cg_aero(:,:,:,:) = init_cginca
1869	!
1870
1871
1872	CALL VTe(VTinca)
1873	CALL VTb(VTphysiq)
1874	#endif
1875	ENDIF
1876	!$omp single
1877	IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, &
1878	press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
1879	!$omp end single
1880	!
1881	!IM betaCRF
1882	pfree=70000. !Pa
1883	beta_pbl=1.
1884	beta_free=1.
1885	lon1_beta=-180.
1886	lon2_beta=+180.
1887	lat1_beta=90.
1888	lat2_beta=-90.
1889	mskocean_beta=.FALSE.
1890
1891	!albedo SB >>>
1892	SELECT CASE(nsw)
1893	CASE(2)
1894	SFRWL(1)=0.45538747
1895	SFRWL(2)=0.54461211
1896	CASE(4)
1897	SFRWL(1)=0.45538747
1898	SFRWL(2)=0.32870591
1899	SFRWL(3)=0.18568763
1900	SFRWL(4)=3.02191470E-02
1901	CASE(6)
1902	SFRWL(1)=1.28432794E-03
1903	SFRWL(2)=0.12304168
1904	SFRWL(3)=0.33106142
1905	SFRWL(4)=0.32870591
1906	SFRWL(5)=0.18568763
1907	SFRWL(6)=3.02191470E-02
1908	END SELECT
1909
1910
1911	!albedo SB <<<
1912
1913	OPEN(99,file='beta_crf.data',status='old', &
1914	form='formatted',err=9999)
1915	READ(99,*,end=9998) pfree
1916	READ(99,*,end=9998) beta_pbl
1917	READ(99,*,end=9998) beta_free
1918	READ(99,*,end=9998) lon1_beta
1919	READ(99,*,end=9998) lon2_beta
1920	READ(99,*,end=9998) lat1_beta
1921	READ(99,*,end=9998) lat2_beta
1922	READ(99,*,end=9998) mskocean_beta
1923	9998 Continue
1924	CLOSE(99)
1925	9999 Continue
1926	WRITE(,)'pfree=',pfree
1927	WRITE(,)'beta_pbl=',beta_pbl
1928	WRITE(,)'beta_free=',beta_free
1929	WRITE(,)'lon1_beta=',lon1_beta
1930	WRITE(,)'lon2_beta=',lon2_beta
1931	WRITE(,)'lat1_beta=',lat1_beta
1932	WRITE(,)'lat2_beta=',lat2_beta
1933	WRITE(,)'mskocean_beta=',mskocean_beta
1934
1935	!lwoff=y : offset LW CRE for radiation code and other schemes
1936	!lwoff=y : betalwoff=1.
1937	betalwoff=0.
1938	IF (ok_lwoff) THEN
1939	betalwoff=1.
1940	ENDIF
1941	WRITE(,)'ok_lwoff=',ok_lwoff
1942	!
1943	!lwoff=y to begin only sollw and sollwdown are set up to CS values
1944	sollw = sollw + betalwoff * (sollw0 - sollw)
1945	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
1946	sollwdown(:))
1947	ENDIF
1948	!
1949	! ************** Fin de IF ( debut ) *************
1950	!
1951	!
1952	! Incrementer le compteur de la physique
1953	!
1954	itap = itap + 1
1955	IF (is_master .OR. prt_level > 9) THEN
1956	IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
1957	WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
1958	WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
1959	100 FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
1960	ENDIF
1961	ENDIF
1962	!
1963	!
1964	! Update fraction of the sub-surfaces (pctsrf) and
1965	! initialize, where a new fraction has appeared, all variables depending
1966	! on the surface fraction.
1967	!
1968	CALL change_srf_frac(itap, phys_tstep, days_elapsed+1, &
1969	pctsrf, fevap, z0m, z0h, agesno, &
1970	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
1971
1972	! Update time and other variables in Reprobus
1973	IF (type_trac == 'repr') THEN
1974	#ifdef REPROBUS
1975	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1976	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1977	CALL Rtime(debut)
1978	#endif
1979	ENDIF
1980
1981	! Tendances bidons pour les processus qui n'affectent pas certaines
1982	! variables.
1983	du0(:,:)=0.
1984	dv0(:,:)=0.
1985	dt0 = 0.
1986	dq0(:,:)=0.
1987	dql0(:,:)=0.
1988	dqi0(:,:)=0.
1989	dsig0(:) = 0.
1990	ddens0(:) = 0.
1991	wkoccur1(:)=1
1992	!
1993	! Mettre a zero des variables de sortie (pour securite)
1994	!
1995	DO i = 1, klon
1996	d_ps(i) = 0.0
1997	ENDDO
1998	DO k = 1, klev
1999	DO i = 1, klon
2000	d_t(i,k) = 0.0
2001	d_u(i,k) = 0.0
2002	d_v(i,k) = 0.0
2003	ENDDO
2004	ENDDO
2005	DO iq = 1, nqtot
2006	DO k = 1, klev
2007	DO i = 1, klon
2008	d_qx(i,k,iq) = 0.0
2009	ENDDO
2010	ENDDO
2011	ENDDO
2012	beta_prec_fisrt(:,:)=0.
2013	beta_prec(:,:)=0.
2014	!
2015	! Output variables from the convective scheme should not be set to 0
2016	! since convection is not always called at every time step.
2017	IF (ok_bug_cv_trac) THEN
2018	da(:,:)=0.
2019	mp(:,:)=0.
2020	phi(:,:,:)=0.
2021	! RomP >>>
2022	phi2(:,:,:)=0.
2023	epmlmMm(:,:,:)=0.
2024	eplaMm(:,:)=0.
2025	d1a(:,:)=0.
2026	dam(:,:)=0.
2027	pmflxr(:,:)=0.
2028	pmflxs(:,:)=0.
2029	! RomP <<<
2030	ENDIF
2031
2032	!
2033	! Ne pas affecter les valeurs entrees de u, v, h, et q
2034	!
2035	DO k = 1, klev
2036	DO i = 1, klon
2037	t_seri(i,k) = t(i,k)
2038	u_seri(i,k) = u(i,k)
2039	v_seri(i,k) = v(i,k)
2040	q_seri(i,k) = qx(i,k,ivap)
2041	ql_seri(i,k) = qx(i,k,iliq)
2042	!CR: ATTENTION, on rajoute la variable glace
2043	IF (nqo.eq.2) THEN
2044	qs_seri(i,k) = 0.
2045	ELSE IF (nqo.eq.3) THEN
2046	qs_seri(i,k) = qx(i,k,isol)
2047	ENDIF
2048	ENDDO
2049	ENDDO
2050	!
2051	!--OB mass fixer
2052	IF (mass_fixer) THEN
2053	!--store initial water burden
2054	qql1(:)=0.0
2055	DO k = 1, klev
2056	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
2057	ENDDO
2058	ENDIF
2059	!--fin mass fixer
2060
2061	tke0(:,:)=pbl_tke(:,:,is_ave)
2062	!CR:Nombre de traceurs de l'eau: nqo
2063	! IF (nqtot.GE.3) THEN
2064	IF (nqtot.GE.(nqo+1)) THEN
2065	! DO iq = 3, nqtot
2066	DO iq = nqo+1, nqtot
2067	DO k = 1, klev
2068	DO i = 1, klon
2069	! tr_seri(i,k,iq-2) = qx(i,k,iq)
2070	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
2071	ENDDO
2072	ENDDO
2073	ENDDO
2074	ELSE
2075	DO k = 1, klev
2076	DO i = 1, klon
2077	tr_seri(i,k,1) = 0.0
2078	ENDDO
2079	ENDDO
2080	ENDIF
2081	!
2082	DO i = 1, klon
2083	ztsol(i) = 0.
2084	ENDDO
2085	DO nsrf = 1, nbsrf
2086	DO i = 1, klon
2087	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
2088	ENDDO
2089	ENDDO
2090	! Initialize variables used for diagnostic purpose
2091	IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
2092
2093	! Diagnostiquer la tendance dynamique
2094	!
2095	IF (ancien_ok) THEN
2096	!
2097	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
2098	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
2099	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
2100	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
2101	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
2102	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
2103	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
2104	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
2105	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
2106	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
2107	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
2108	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
2109	! !! RomP >>> td dyn traceur
2110	IF (nqtot.GT.nqo) THEN ! jyg
2111	DO iq = nqo+1, nqtot ! jyg
2112	d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/phys_tstep ! jyg
2113	ENDDO
2114	ENDIF
2115	! !! RomP <<<
2116	ELSE
2117	d_u_dyn(:,:) = 0.0
2118	d_v_dyn(:,:) = 0.0
2119	d_t_dyn(:,:) = 0.0
2120	d_q_dyn(:,:) = 0.0
2121	d_ql_dyn(:,:) = 0.0
2122	d_qs_dyn(:,:) = 0.0
2123	d_q_dyn2d(:) = 0.0
2124	d_ql_dyn2d(:) = 0.0
2125	d_qs_dyn2d(:) = 0.0
2126	! !! RomP >>> td dyn traceur
2127	IF (nqtot.GT.nqo) THEN ! jyg
2128	DO iq = nqo+1, nqtot ! jyg
2129	d_tr_dyn(:,:,iq-nqo)= 0.0 ! jyg
2130	ENDDO
2131	ENDIF
2132	! !! RomP <<<
2133	ancien_ok = .TRUE.
2134	ENDIF
2135	!
2136	! Ajouter le geopotentiel du sol:
2137	!
2138	DO k = 1, klev
2139	DO i = 1, klon
2140	zphi(i,k) = pphi(i,k) + pphis(i)
2141	ENDDO
2142	ENDDO
2143	!
2144	! Verifier les temperatures
2145	!
2146	!IM BEG
2147	IF (check) THEN
2148	amn=MIN(ftsol(1,is_ter),1000.)
2149	amx=MAX(ftsol(1,is_ter),-1000.)
2150	DO i=2, klon
2151	amn=MIN(ftsol(i,is_ter),amn)
2152	amx=MAX(ftsol(i,is_ter),amx)
2153	ENDDO
2154	!
2155	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
2156	ENDIF !(check) THEN
2157	!IM END
2158	!
2159	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
2160	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
2161
2162	!
2163	!IM BEG
2164	IF (check) THEN
2165	amn=MIN(ftsol(1,is_ter),1000.)
2166	amx=MAX(ftsol(1,is_ter),-1000.)
2167	DO i=2, klon
2168	amn=MIN(ftsol(i,is_ter),amn)
2169	amx=MAX(ftsol(i,is_ter),amx)
2170	ENDDO
2171	!
2172	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2173	ENDIF !(check) THEN
2174	!IM END
2175	!
2176	! Mettre en action les conditions aux limites (albedo, sst, etc.).
2177	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
2178	!
2179	! Update ozone if day change
2180	IF (MOD(itap-1,lmt_pas) == 0) THEN
2181	IF (read_climoz <= 0) THEN
2182	! Once per day, update ozone from Royer:
2183	IF (solarlong0<-999.) then
2184	! Generic case with evolvoing season
2185	zzz=real(days_elapsed+1)
2186	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
2187	! Particular case with annual mean insolation
2188	zzz=real(90) ! could be revisited
2189	IF (read_climoz/=-1) THEN
2190	abort_message ='read_climoz=-1 is recommended when ' &
2191	// 'solarlong0=1000.'
2192	CALL abort_physic (modname,abort_message,1)
2193	ENDIF
2194	ELSE
2195	! Case where the season is imposed with solarlong0
2196	zzz=real(90) ! could be revisited
2197	ENDIF
2198
2199	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
2200	ELSE
2201	!--- ro3i = elapsed days number since current year 1st january, 0h
2202	ro3i=days_elapsed+jh_cur-jh_1jan
2203	!--- scaling for old style files (360 records)
2204	IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
2205	IF(adjust_tropopause) THEN
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 , longitude_deg, latitude_deg, &
2209	dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
2210	ELSE
2211	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2212	ro3i, 'C', press_cen_climoz, pplay, wo, paprs(:,1), &
2213	time_climoz )
2214	ENDIF
2215	! Convert from mole fraction of ozone to column density of ozone in a
2216	! cell, in kDU:
2217	FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2218	* zmasse / dobson_u / 1e3
2219	! (By regridding ozone values for LMDZ only once a day, we
2220	! have already neglected the variation of pressure in one
2221	! day. So do not recompute "wo" at each time step even if
2222	! "zmasse" changes a little.)
2223	ENDIF
2224	ENDIF
2225	!
2226	! Re-evaporer l'eau liquide nuageuse
2227	!
2228	CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2229	& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
2230
2231	CALL add_phys_tend &
2232	(du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
2233	'eva',abortphy,flag_inhib_tend,itap,0)
2234	CALL prt_enerbil('eva',itap)
2235
2236	!=========================================================================
2237	! Calculs de l'orbite.
2238	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2239	! doit donc etre plac\'e avant radlwsw et pbl_surface
2240
2241	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2242	CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2243	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2244	!
2245	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2246	! solarlong0
2247	IF (solarlong0<-999.) THEN
2248	IF (new_orbit) THEN
2249	! calcul selon la routine utilisee pour les planetes
2250	CALL solarlong(day_since_equinox, zlongi, dist)
2251	ELSE
2252	! calcul selon la routine utilisee pour l'AR4
2253	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2254	ENDIF
2255	ELSE
2256	zlongi=solarlong0 ! longitude solaire vraie
2257	dist=1. ! distance au soleil / moyenne
2258	ENDIF
2259
2260	IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2261
2262
2263	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2264	! Calcul de l'ensoleillement :
2265	! ============================
2266	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2267	! l'annee a partir d'une formule analytique.
2268	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2269	! non nul aux poles.
2270	IF (abs(solarlong0-1000.)<1.e-4) THEN
2271	CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
2272	latitude_deg,longitude_deg,rmu0,fract)
2273	swradcorr(:) = 1.0
2274	JrNt(:) = 1.0
2275	zrmu0(:) = rmu0(:)
2276	ELSE
2277	! recode par Olivier Boucher en sept 2015
2278	SELECT CASE (iflag_cycle_diurne)
2279	CASE(0)
2280	! Sans cycle diurne
2281	CALL angle(zlongi, latitude_deg, fract, rmu0)
2282	swradcorr = 1.0
2283	JrNt = 1.0
2284	zrmu0 = rmu0
2285	CASE(1)
2286	! Avec cycle diurne sans application des poids
2287	! bit comparable a l ancienne formulation cycle_diurne=true
2288	! on integre entre gmtime et gmtime+radpas
2289	zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
2290	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2291	latitude_deg,longitude_deg,rmu0,fract)
2292	zrmu0 = rmu0
2293	swradcorr = 1.0
2294	! Calcul du flag jour-nuit
2295	JrNt = 0.0
2296	WHERE (fract.GT.0.0) JrNt = 1.0
2297	CASE(2)
2298	! Avec cycle diurne sans application des poids
2299	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2300	! Comme cette routine est appele a tous les pas de temps de
2301	! la physique meme si le rayonnement n'est pas appele je
2302	! remonte en arriere les radpas-1 pas de temps
2303	! suivant. Petite ruse avec MOD pour prendre en compte le
2304	! premier pas de temps de la physique pendant lequel
2305	! itaprad=0
2306	zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)
2307	zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
2308	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2309	latitude_deg,longitude_deg,rmu0,fract)
2310	!
2311	! Calcul des poids
2312	!
2313	zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
2314	zdtime2=0.0 !--pas de temps de la physique qui se termine
2315	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2316	latitude_deg,longitude_deg,zrmu0,zfract)
2317	swradcorr = 0.0
2318	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2319	swradcorr=zfract/fract*zrmu0/rmu0
2320	! Calcul du flag jour-nuit
2321	JrNt = 0.0
2322	WHERE (zfract.GT.0.0) JrNt = 1.0
2323	END SELECT
2324	ENDIF
2325	sza_o = ACOS (rmu0) *180./pi
2326
2327	IF (mydebug) THEN
2328	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2329	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2330	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2331	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2332	ENDIF
2333
2334	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2335	! Appel au pbl_surface : Planetary Boudary Layer et Surface
2336	! Cela implique tous les interactions des sous-surfaces et la
2337	! partie diffusion turbulent du couche limit.
2338	!
2339	! Certains varibales de sorties de pbl_surface sont utiliser que pour
2340	! ecriture des fihiers hist_XXXX.nc, ces sont :
2341	! qsol, zq2m, s_pblh, s_lcl,
2342	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2343	! s_therm, s_trmb1, s_trmb2, s_trmb3,
2344	! zu10m, zv10m, fder,
2345	! zxqsurf, rh2m, zxfluxu, zxfluxv,
2346	! frugs, agesno, fsollw, fsolsw,
2347	! d_ts, fevap, fluxlat, t2m,
2348	! wfbils, wfbilo, fluxt, fluxu, fluxv,
2349	!
2350	! Certains ne sont pas utiliser du tout :
2351	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2352	!
2353
2354	! Calcul de l'humidite de saturation au niveau du sol
2355
2356
2357
2358	IF (iflag_pbl/=0) THEN
2359
2360	!jyg+nrlmd<
2361	!!jyg IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2362	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
2363	print *,'debut du splitting de la PBL, wake_s = ', wake_s(:)
2364	print *,'debut du splitting de la PBL, wake_deltat = ', wake_deltat(:,1)
2365	print *,'debut du splitting de la PBL, wake_deltaq = ', wake_deltaq(:,1)
2366	ENDIF
2367	! !!
2368	!>jyg+nrlmd
2369	!
2370	!-------gustiness calculation-------!
2371	!ym : Warning gustiness non inialized for iflag_gusts=2 & iflag_gusts=3
2372	gustiness=0 !ym missing init
2373
2374	IF (iflag_gusts==0) THEN
2375	gustiness(1:klon)=0
2376	ELSE IF (iflag_gusts==1) THEN
2377	gustiness(1:klon)=f_gust_blale_bl(1:klon)+f_gust_wkale_wake(1:klon)
2378	ELSE IF (iflag_gusts==2) THEN
2379	gustiness(1:klon)=f_gust_blale_bl_stat(1:klon)+f_gust_wkale_wake(1:klon)
2380	! ELSE IF (iflag_gusts==2) THEN
2381	! do i = 1, klon
2382	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
2383	! *ale_wake(i) !! need to make sigma_wk accessible here
2384	! enddo
2385	! ELSE IF (iflag_gusts==3) THEN
2386	! do i = 1, klon
2387	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
2388	! enddo
2389	ENDIF
2390
2391	CALL pbl_surface( &
2392	phys_tstep, date0, itap, days_elapsed+1, &
2393	debut, lafin, &
2394	longitude_deg, latitude_deg, rugoro, zrmu0, &
2395	zsig, sollwdown, pphi, cldt, &
2396	rain_fall, snow_fall, solsw, sollw, &
2397	gustiness, &
2398	t_seri, q_seri, u_seri, v_seri, &
2399	!nrlmd+jyg<
2400	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2401	!>nrlmd+jyg
2402	pplay, paprs, pctsrf, &
2403	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2404	!albedo SB <<<
2405	cdragh, cdragm, u1, v1, &
2406	!albedo SB >>>
2407	! albsol1, albsol2, sens, evap, &
2408	albsol_dir, albsol_dif, sens, evap, &
2409	!albedo SB <<<
2410	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
2411	zxtsol, zxfluxlat, zt2m, qsat2m, &
2412	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2413	!nrlmd<
2414	!jyg<
2415	d_t_vdf_w, d_q_vdf_w, &
2416	d_t_vdf_x, d_q_vdf_x, &
2417	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2418	!>jyg
2419	delta_tsurf,wake_dens, &
2420	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2421	kh,kh_x,kh_w, &
2422	!>nrlmd
2423	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2424	slab_wfbils, &
2425	qsol, zq2m, s_pblh, s_lcl, &
2426	!jyg<
2427	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2428	!>jyg
2429	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2430	s_therm, s_trmb1, s_trmb2, s_trmb3, &
2431	zustar, zu10m, zv10m, fder, &
2432	zxqsurf, rh2m, zxfluxu, zxfluxv, &
2433	z0m, z0h, agesno, fsollw, fsolsw, &
2434	d_ts, fevap, fluxlat, t2m, &
2435	wfbils, wfbilo, wfevap, wfrain, wfsnow, &
2436	fluxt, fluxu, fluxv, &
2437	dsens, devap, zxsnow, &
2438	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
2439	!nrlmd+jyg<
2440	wake_delta_pbl_TKE, &
2441	!>nrlmd+jyg
2442	treedrg )
2443	!FC
2444	!
2445	! Add turbulent diffusion tendency to the wake difference variables
2446	!!jyg IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2447	IF (mod(iflag_pbl_split,10) .NE. 0) THEN
2448	!jyg<
2449	d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2450	d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2451	CALL add_wake_tend &
2452	(d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, ddens0, wkoccur1, 'vdf', abortphy)
2453	ELSE
2454	d_deltat_vdf(:,:) = 0.
2455	d_deltaq_vdf(:,:) = 0.
2456	!>jyg
2457	ENDIF
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	addtkeoro=0
4394	CALL getin_p('addtkeoro',addtkeoro)
4395
4396	IF (prt_level.ge.5) &
4397	print*,'addtkeoro', addtkeoro
4398
4399	alphatkeoro=1.
4400	CALL getin_p('alphatkeoro',alphatkeoro)
4401	alphatkeoro=min(max(0.,alphatkeoro),1.)
4402
4403	smallscales_tkeoro=.FALSE.
4404	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4405
4406
4407	dtadd(:,:)=0.
4408	duadd(:,:)=0.
4409	dvadd(:,:)=0.
4410
4411	! Choices for addtkeoro:
4412	! ** 0 no TKE tendency from orography
4413	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4414	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4415	!
4416
4417	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4418	! -------------------------------------------
4419
4420
4421	! selection des points pour lesquels le schema est actif:
4422
4423
4424	IF (addtkeoro .EQ. 1 ) THEN
4425
4426	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4427	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4428
4429	ELSE IF (addtkeoro .EQ. 2) THEN
4430
4431	IF (smallscales_tkeoro) THEN
4432	igwd=0
4433	DO i=1,klon
4434	itest(i)=0
4435	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4436	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4437	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4438	IF (zstd(i).GT.1.0) THEN
4439	itest(i)=1
4440	igwd=igwd+1
4441	idx(igwd)=i
4442	ENDIF
4443	ENDDO
4444
4445	ELSE
4446
4447	igwd=0
4448	DO i=1,klon
4449	itest(i)=0
4450	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4451	itest(i)=1
4452	igwd=igwd+1
4453	idx(igwd)=i
4454	ENDIF
4455	ENDDO
4456
4457	ENDIF
4458
4459	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
4460	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4461	igwd,idx,itest, &
4462	t_seri, u_seri, v_seri, &
4463	zulow, zvlow, zustrdr, zvstrdr, &
4464	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4465
4466	zustrdr(:)=0.
4467	zvstrdr(:)=0.
4468	zulow(:)=0.
4469	zvlow(:)=0.
4470
4471	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4472	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4473	ENDIF
4474
4475
4476	! TKE update from subgrid temperature and wind tendencies
4477	!----------------------------------------------------------
4478	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4479
4480
4481	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
4482
4483
4484	ENDIF
4485	! -----
4486	!===============================================================
4487
4488
4489	!====================================================================
4490	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4491	!====================================================================
4492	! Abderrahmane 24.08.09
4493
4494	IF (ok_cosp) THEN
4495	! adeclarer
4496	#ifdef CPP_COSP
4497	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4498
4499	IF (prt_level .GE.10) THEN
4500	print*,'freq_cosp',freq_cosp
4501	ENDIF
4502	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4503	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4504	! s ref_liq,ref_ice
4505	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
4506	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4507	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4508	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4509	JrNt,ref_liq,ref_ice, &
4510	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4511	zu10m,zv10m,pphis, &
4512	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4513	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4514	prfl(:,1:klev),psfl(:,1:klev), &
4515	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4516	mr_ozone,cldtau, cldemi)
4517
4518	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4519	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4520	! M clMISR,
4521	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4522	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4523
4524	ENDIF
4525	#endif
4526
4527	#ifdef CPP_COSP2
4528	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4529
4530	IF (prt_level .GE.10) THEN
4531	print*,'freq_cosp',freq_cosp
4532	ENDIF
4533	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4534	print*,'Dans physiq.F avant appel '
4535	! s ref_liq,ref_ice
4536	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
4537	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4538	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4539	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4540	JrNt,ref_liq,ref_ice, &
4541	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4542	zu10m,zv10m,pphis, &
4543	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4544	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4545	prfl(:,1:klev),psfl(:,1:klev), &
4546	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4547	mr_ozone,cldtau, cldemi)
4548	ENDIF
4549	#endif
4550
4551	ENDIF !ok_cosp
4552
4553
4554	! Marine
4555
4556	IF (ok_airs) then
4557
4558	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
4559	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4560	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4561	& map_prop_hc,map_prop_hist,&
4562	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4563	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4564	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4565	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4566	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4567	& map_ntot,map_hc,map_hist,&
4568	& map_Cb,map_ThCi,map_Anv,&
4569	& alt_tropo )
4570	ENDIF
4571
4572	ENDIF ! ok_airs
4573
4574
4575	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4576	!AA
4577	!AA Installation de l'interface online-offline pour traceurs
4578	!AA
4579	!====================================================================
4580	! Calcul des tendances traceurs
4581	!====================================================================
4582	!
4583
4584	IF (type_trac=='repr') THEN
4585	sh_in(:,:) = q_seri(:,:)
4586	ELSE
4587	sh_in(:,:) = qx(:,:,ivap)
4588	ch_in(:,:) = qx(:,:,iliq)
4589	ENDIF
4590
4591	IF (iflag_phytrac == 1 ) THEN
4592
4593	#ifdef CPP_Dust
4594	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
4595	pdtphys,ftsol, & ! I
4596	t,q_seri,paprs,pplay,RHcl, & ! I
4597	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
4598	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
4599	u_seri, v_seri, latitude_deg, longitude_deg, &
4600	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
4601	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
4602	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
4603	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
4604	fm_therm, entr_therm, rneb, & ! I
4605	beta_prec_fisrt,beta_prec, & !I
4606	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
4607	d_tr_dyn,tr_seri)
4608
4609	#else
4610
4611	CALL phytrac ( &
4612	itap, days_elapsed+1, jH_cur, debut, &
4613	lafin, phys_tstep, u, v, t, &
4614	paprs, pplay, pmfu, pmfd, &
4615	pen_u, pde_u, pen_d, pde_d, &
4616	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
4617	u1, v1, ftsol, pctsrf, &
4618	zustar, zu10m, zv10m, &
4619	wstar(:,is_ave), ale_bl, ale_wake, &
4620	latitude_deg, longitude_deg, &
4621	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4622	presnivs, pphis, pphi, albsol1, &
4623	sh_in, ch_in, rhcl, cldfra, rneb, &
4624	diafra, cldliq, itop_con, ibas_con, &
4625	pmflxr, pmflxs, prfl, psfl, &
4626	da, phi, mp, upwd, &
4627	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
4628	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
4629	ev, ep, epmlmMm, eplaMm, & !<<RomP
4630	dnwd, aerosol_couple, flxmass_w, &
4631	tau_aero, piz_aero, cg_aero, ccm, &
4632	rfname, &
4633	d_tr_dyn, & !<<RomP
4634	tr_seri, init_source)
4635	#endif
4636	ENDIF ! (iflag_phytrac=1)
4637
4638	IF (offline) THEN
4639
4640	IF (prt_level.ge.9) &
4641	print*,'Attention on met a 0 les thermiques pour phystoke'
4642	CALL phystokenc ( &
4643	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4644	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4645	fm_therm,entr_therm, &
4646	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4647	frac_impa, frac_nucl, &
4648	pphis,cell_area,phys_tstep,itap, &
4649	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4650
4651
4652	ENDIF
4653
4654	!
4655	! Calculer le transport de l'eau et de l'energie (diagnostique)
4656	!
4657	CALL transp (paprs,zxtsol, &
4658	t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
4659	ve, vq, ue, uq, vwat, uwat)
4660	!
4661	!IM global posePB BEG
4662	IF(1.EQ.0) THEN
4663	!
4664	CALL transp_lay (paprs,zxtsol, &
4665	t_seri, q_seri, u_seri, v_seri, zphi, &
4666	ve_lay, vq_lay, ue_lay, uq_lay)
4667	!
4668	ENDIF !(1.EQ.0) THEN
4669	!IM global posePB END
4670	! Accumuler les variables a stocker dans les fichiers histoire:
4671	!
4672
4673	!================================================================
4674	! Conversion of kinetic and potential energy into heat, for
4675	! parameterisation of subgrid-scale motions
4676	!================================================================
4677
4678	d_t_ec(:,:)=0.
4679	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4680	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
4681	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4682	zmasse,exner,d_t_ec)
4683	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4684
4685	!=======================================================================
4686	! SORTIES
4687	!=======================================================================
4688	!
4689	!IM initialisation + calculs divers diag AMIP2
4690	!
4691	include "calcul_divers.h"
4692	!
4693	!IM Interpolation sur les niveaux de pression du NMC
4694	! -------------------------------------------------
4695	!
4696	include "calcul_STDlev.h"
4697	!
4698	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
4699	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
4700	!
4701	!cc prw = eau precipitable
4702	! prlw = colonne eau liquide
4703	! prlw = colonne eau solide
4704	prw(:) = 0.
4705	prlw(:) = 0.
4706	prsw(:) = 0.
4707	DO k = 1, klev
4708	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
4709	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
4710	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
4711	ENDDO
4712	!
4713	IF (type_trac == 'inca') THEN
4714	#ifdef INCA
4715	CALL VTe(VTphysiq)
4716	CALL VTb(VTinca)
4717
4718	CALL chemhook_end ( &
4719	phys_tstep, &
4720	pplay, &
4721	t_seri, &
4722	tr_seri, &
4723	nbtr, &
4724	paprs, &
4725	q_seri, &
4726	cell_area, &
4727	pphi, &
4728	pphis, &
4729	zx_rh, &
4730	aps, bps)
4731
4732	CALL VTe(VTinca)
4733	CALL VTb(VTphysiq)
4734	#endif
4735	ENDIF
4736
4737
4738	!
4739	! Convertir les incrementations en tendances
4740	!
4741	IF (prt_level .GE.10) THEN
4742	print *,'Convertir les incrementations en tendances '
4743	ENDIF
4744	!
4745	IF (mydebug) THEN
4746	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4747	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4748	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4749	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4750	ENDIF
4751
4752	DO k = 1, klev
4753	DO i = 1, klon
4754	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
4755	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
4756	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
4757	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
4758	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
4759	!CR: on ajoute le contenu en glace
4760	IF (nqo.eq.3) THEN
4761	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
4762	ENDIF
4763	ENDDO
4764	ENDDO
4765	!
4766	!CR: nb de traceurs eau: nqo
4767	! IF (nqtot.GE.3) THEN
4768	IF (nqtot.GE.(nqo+1)) THEN
4769	! DO iq = 3, nqtot
4770	DO iq = nqo+1, nqtot
4771	DO k = 1, klev
4772	DO i = 1, klon
4773	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / phys_tstep
4774	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / phys_tstep
4775	ENDDO
4776	ENDDO
4777	ENDDO
4778	ENDIF
4779	!
4780	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4781	!IM global posePB include "write_bilKP_ins.h"
4782	!IM global posePB include "write_bilKP_ave.h"
4783	!
4784
4785	!--OB mass fixer
4786	!--profile is corrected to force mass conservation of water
4787	IF (mass_fixer) THEN
4788	qql2(:)=0.0
4789	DO k = 1, klev
4790	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
4791	ENDDO
4792	DO i = 1, klon
4793	!--compute ratio of what q+ql should be with conservation to what it is
4794	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
4795	DO k = 1, klev
4796	q_seri(i,k) =q_seri(i,k)*corrqql
4797	ql_seri(i,k)=ql_seri(i,k)*corrqql
4798	ENDDO
4799	ENDDO
4800	ENDIF
4801	!--fin mass fixer
4802
4803	! Sauvegarder les valeurs de t et q a la fin de la physique:
4804	!
4805	u_ancien(:,:) = u_seri(:,:)
4806	v_ancien(:,:) = v_seri(:,:)
4807	t_ancien(:,:) = t_seri(:,:)
4808	q_ancien(:,:) = q_seri(:,:)
4809	ql_ancien(:,:) = ql_seri(:,:)
4810	qs_ancien(:,:) = qs_seri(:,:)
4811	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
4812	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
4813	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
4814	! !! RomP >>>
4815	!CR: nb de traceurs eau: nqo
4816	IF (nqtot.GT.nqo) THEN
4817	DO iq = nqo+1, nqtot
4818	tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
4819	ENDDO
4820	ENDIF
4821	! !! RomP <<<
4822	!==========================================================================
4823	! Sorties des tendances pour un point particulier
4824	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4825	! pour le debug
4826	! La valeur de igout est attribuee plus haut dans le programme
4827	!==========================================================================
4828
4829	IF (prt_level.ge.1) THEN
4830	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4831	write(lunout,*) &
4832	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4833	write(lunout,*) &
4834	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4835	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4836	pctsrf(igout,is_sic)
4837	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4838	DO k=1,klev
4839	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4840	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4841	d_t_eva(igout,k)
4842	ENDDO
4843	write(lunout,*) 'cool,heat'
4844	DO k=1,klev
4845	write(lunout,*) cool(igout,k),heat(igout,k)
4846	ENDDO
4847
4848	!jyg< (En attendant de statuer sur le sort de d_t_oli)
4849	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4850	!jyg! do k=1,klev
4851	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4852	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4853	!jyg! enddo
4854	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4855	DO k=1,klev
4856	write(lunout,*) d_t_vdf(igout,k), &
4857	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4858	ENDDO
4859	!>jyg
4860
4861	write(lunout,*) 'd_ps ',d_ps(igout)
4862	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4863	DO k=1,klev
4864	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4865	d_qx(igout,k,1),d_qx(igout,k,2)
4866	ENDDO
4867	ENDIF
4868
4869	!============================================================
4870	! Calcul de la temperature potentielle
4871	!============================================================
4872	DO k = 1, klev
4873	DO i = 1, klon
4874	!JYG/IM theta en debut du pas de temps
4875	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4876	!JYG/IM theta en fin de pas de temps de physique
4877	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4878	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
4879	! MPL 20130625
4880	! fth_fonctions.F90 et parkind1.F90
4881	! sinon thetal=theta
4882	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4883	! : ql_seri(i,k))
4884	thetal(i,k)=theta(i,k)
4885	ENDDO
4886	ENDDO
4887	!
4888
4889	! 22.03.04 BEG
4890	!=============================================================
4891	! Ecriture des sorties
4892	!=============================================================
4893	#ifdef CPP_IOIPSL
4894
4895	! Recupere des varibles calcule dans differents modules
4896	! pour ecriture dans histxxx.nc
4897
4898	! Get some variables from module fonte_neige_mod
4899	CALL fonte_neige_get_vars(pctsrf, &
4900	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
4901
4902
4903	!=============================================================
4904	! Separation entre thermiques et non thermiques dans les sorties
4905	! de fisrtilp
4906	!=============================================================
4907
4908	IF (iflag_thermals>=1) THEN
4909	d_t_lscth=0.
4910	d_t_lscst=0.
4911	d_q_lscth=0.
4912	d_q_lscst=0.
4913	DO k=1,klev
4914	DO i=1,klon
4915	IF (ptconvth(i,k)) THEN
4916	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4917	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4918	ELSE
4919	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4920	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4921	ENDIF
4922	ENDDO
4923	ENDDO
4924
4925	DO i=1,klon
4926	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4927	plul_th(i)=prfl(i,1)+psfl(i,1)
4928	ENDDO
4929	ENDIF
4930
4931	!On effectue les sorties:
4932
4933	#ifdef CPP_Dust
4934	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
4935	pplay, lmax_th, aerosol_couple, &
4936	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
4937	ptconv, read_climoz, clevSTD, &
4938	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
4939	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4940	#else
4941	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
4942	pplay, lmax_th, aerosol_couple, &
4943	ok_ade, ok_aie, ivap, iliq, isol, new_aod, &
4944	ok_sync, ptconv, read_climoz, clevSTD, &
4945	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
4946	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4947	#endif
4948
4949	#ifndef CPP_XIOS
4950	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
4951	#endif
4952
4953	#endif
4954
4955	! On remet des variables a .false. apres un premier appel
4956	IF (debut) THEN
4957	#ifdef CPP_XIOS
4958	swaero_diag=.FALSE.
4959	swaerofree_diag=.FALSE.
4960	dryaod_diag=.FALSE.
4961	ok_4xCO2atm= .FALSE.
4962	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
4963
4964	IF (is_master) THEN
4965	!--setting up swaero_diag to TRUE in XIOS case
4966	IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
4967	xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
4968	xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR. &
4969	(iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
4970	xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0")))) &
4971	!!!--for now these fields are not in the XML files so they are omitted
4972	!!! xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
4973	swaero_diag=.TRUE.
4974
4975	!--setting up swaerofree_diag to TRUE in XIOS case
4976	IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
4977	xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR. &
4978	xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
4979	xios_field_is_active("LWupTOAcleanclr")) &
4980	swaerofree_diag=.TRUE.
4981
4982	!--setting up dryaod_diag to TRUE in XIOS case
4983	DO naero = 1, naero_tot-1
4984	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
4985	ENDDO
4986	!
4987	!--setting up ok_4xCO2atm to TRUE in XIOS case
4988	IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
4989	xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
4990	xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
4991	xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
4992	xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
4993	xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
4994	ok_4xCO2atm=.TRUE.
4995	ENDIF
4996	!$OMP BARRIER
4997	CALL bcast(swaero_diag)
4998	CALL bcast(swaerofree_diag)
4999	CALL bcast(dryaod_diag)
5000	CALL bcast(ok_4xCO2atm)
5001	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
5002	#endif
5003	ENDIF
5004
5005	!====================================================================
5006	! Arret du modele apres hgardfou en cas de detection d'un
5007	! plantage par hgardfou
5008	!====================================================================
5009
5010	IF (abortphy==1) THEN
5011	abort_message ='Plantage hgardfou'
5012	CALL abort_physic (modname,abort_message,1)
5013	ENDIF
5014
5015	! 22.03.04 END
5016	!
5017	!====================================================================
5018	! Si c'est la fin, il faut conserver l'etat de redemarrage
5019	!====================================================================
5020	!
5021
5022	IF (lafin) THEN
5023	itau_phy = itau_phy + itap
5024	CALL phyredem ("restartphy.nc")
5025	! open(97,form="unformatted",file="finbin")
5026	! write(97) u_seri,v_seri,t_seri,q_seri
5027	! close(97)
5028
5029	IF (is_omp_master) THEN
5030
5031	IF (read_climoz >= 1) THEN
5032	IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5033	DEALLOCATE(press_edg_climoz) ! pointer
5034	DEALLOCATE(press_cen_climoz) ! pointer
5035	ENDIF
5036
5037	ENDIF
5038	#ifdef CPP_XIOS
5039	IF (is_omp_master) CALL xios_context_finalize
5040	#endif
5041	WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
5042	ENDIF
5043
5044	! first=.false.
5045
5046	END SUBROUTINE physiq
5047
5048	END MODULE physiq_mod

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

Download in other formats:

Original Format