Context Navigation

physiq_mod.F90 @ 3340

Last change on this file since 3340 was 3340, checked in by acozic, 6 years ago
merge with trunk 3338 Add a flag chemistry_couple
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: 161.1 KB

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

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

Download in other formats:

Original Format