Context Navigation

physiq_mod.F90 @ 2866

Last change on this file since 2866 was 2618, checked in by oboucher, 8 years ago
I have moved the computation of cloud optical properties into the radiation if .. .endif sequence, after the aerosol optical properties. This follows up on rev 2614 and should fix the 1+1=2 problem that was introduced back then.
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: 160.6 KB

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

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

Download in other formats:

Original Format