Context Navigation

physiq_mod.F90 @ 3480

Last change on this file since 3480 was 3480, checked in by Laurent Fairhead, 5 years ago
Integration of revisions 3425 through 3427 from IPSLCM6.0.15 branch into the trunk 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.3 KB

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

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

Download in other formats:

Original Format