Context Navigation

physiq_mod.F90 @ 3465

Last change on this file since 3465 was 3465, checked in by Laurent Fairhead, 5 years ago
Further modifications for DYNAMICO/LMDZ convergence. These are based on Yann's LMDZ6_V2 sources. Compiles on irene and converges with revision 3459 in a bucket configuration YM/LF
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: 178.0 KB

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

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

Download in other formats:

Original Format