Context Navigation

physiq_mod.F90 @ 5445

Last change on this file since 5445 was 3406, checked in by jghattas, 6 years ago
Added all modifications in the model code that were used for the simulations with DYANMICO during the Grand Challeng 2018. Modifications done by Y. Meurdesoif, L. Fairhead and A.K. Traore
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 `Id`
File size: 165.2 KB

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

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

Download in other formats:

Original Format