Context Navigation

physiq_mod.F90 @ 2720

Last change on this file since 2720 was 2720, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2664:2719 into testing branch
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: 150.9 KB

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

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

Download in other formats:

Original Format