Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

physiq_mod.F90

Last change on this file was 4369, checked in by lguez, 18 months ago
Sync latest trunk changes to branch Ocean_skin
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: 192.5 KB

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

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

Download in other formats:

Original Format