Context Navigation

physiq_mod.F90 @ 5342

Last change on this file since 5342 was 3000, checked in by jyg, 7 years ago

(i) Implementation of the conditionning of the Alp
provided by the wakes: when iflag_alp_wk_cond > 0,
Alp_wk is divided by the probability that there is
a gust front within the grid cell.

New subroutine alpale_wk.F90 called by physiq.

(ii) Some changes concerning the initialization of
the variable wake_dens.

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: 162.8 KB

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

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

Download in other formats:

Original Format