Context Navigation

physiq_mod.F90 @ 2831

Last change on this file since 2831 was 2831, checked in by acozic, 7 years ago
add an argument date0 to inca initialisation
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: 156.2 KB

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

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

Download in other formats:

Original Format