Context Navigation

physiq_mod.F90 @ 4056

Last change on this file since 4056 was 4056, checked in by dcugnet, 2 years ago

Most of the changes are intended to help to eventually remove the constraints about the tracers assumptions, in particular water tracers.

Remove index tables itr_indice and niadv, replaced by tracers(:)%isAdvected and tracers(:)%isH2OFamily. Most of the loops are now from 1 to nqtot:
- DO iq=nqo+1,nqtot loops are replaced with: DO iq=1,nqtot

IF(tracers(iq)%isH2Ofamily) CYCLE

DO it=1,nbtr; iq=niadv(it+nqo)

and DO it=1,nqtottr; iq=itr_indice(it) loops are replaced with:

it = 0
DO iq = 1, nqtot

IF(.NOT.tracers(iq)%isAdvected .OR. tracers(iq)%isH2Ofamily) CYCLE
it = it+1

Move some StratAer? related code from infotrac to infotrac_phy
Remove "nqperes" variable:

DO iq=1,nqpere loops are replaced with:
DO iq=1,nqtot

IF(tracers(iq)%parent/='air') CYCLE

Cosmetic changes (justification, SELECT CASE instead of multiple IF...) mostly in advtrac* routines.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Id

File size: 187.7 KB

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

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

Download in other formats:

Original Format