Context Navigation

physiq_mod.F90 @ 3712

Last change on this file since 3712 was 3706, checked in by adurocher, 4 years ago
Added timers for physiq and display physic profiling
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: 184.7 KB

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

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

Download in other formats:

Original Format