Context Navigation

physiq_mod.F90 @ 3296

Last change on this file since 3296 was 3274, checked in by oboucher, 7 years ago
Implementing the MACspV2 aerosol plume climatology which can be called by setting flag_aerosol=7 and aerosols1980.nc pointing to aerosols.nat.nc. Also requires the MACv2.0-SP_v1.nc file.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 170.6 KB

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

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

Download in other formats:

Original Format