Context Navigation

physiq_mod.F90 @ 3481

Last change on this file since 3481 was 3481, checked in by acozic, 5 years ago
call to aerosol_meteo_calc is now done for all inca chemistries configurations
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: 171.9 KB

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

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

Download in other formats:

Original Format