Context Navigation

physiq_mod.F90 @ 3489

Last change on this file since 3489 was 3489, checked in by musat, 5 years ago
Ajout bornage a 2m
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 179.5 KB

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

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

Download in other formats:

Original Format