Context Navigation

physiq_mod.F90 @ 4009

Last change on this file since 4009 was 4009, checked in by evignon, 3 years ago
! Prise en compte de l'orographie sous maille, des heterogeneites de surface ! sur le ratqs + nouvelle version des ratqs interactifs de Louis ! Le tout est dan un module: calcratqs_multi_mod. ! Pour l'instant, les nouvelles contributions peuvent s'activer ! uniquement de facon separee
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 187.1 KB

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

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

Download in other formats:

Original Format