Context Navigation

physiq_mod.F90 @ 2675

Last change on this file since 2675 was 2675, checked in by oboucher, 8 years ago
Dealing with initialisation of swaero_diag in the case of XIOS outputs.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 157.2 KB

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

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

Download in other formats:

Original Format