Context Navigation

physiq_mod.F90 @ 2695

Last change on this file since 2695 was 2692, checked in by oboucher, 8 years ago
Adding a call to stratosphere_mask in the case of StratAer?. A lot of cosmetic changes on if/endif do/enddo and calls
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: 157.2 KB

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

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

Download in other formats:

Original Format