Context Navigation

physiq_mod.F90 @ 3482

Last change on this file since 3482 was 3482, checked in by acozic, 6 years ago
call to aerosol_meteo_calc is now done for all inca chemistries configurations
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 161.8 KB

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

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

Download in other formats:

Original Format