Context Navigation

physiq_mod.F90 @ 3112

Last change on this file since 3112 was 3112, checked in by musat, 6 years ago
Add LMDZ outputs: icc3dcon, icc3dstra Add CMIP6 outputs : clic, clis, rsdsdiff, rsdscsdiff Correct parasolRefl_sea Add "HomeCMIP6" variable HomeCMIP6_parasolCRefl IM
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: 165.2 KB

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

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

Download in other formats:

Original Format