Context Navigation

physiq_mod.F90 @ 3630

Last change on this file since 3630 was 3630, checked in by Laurent Fairhead, 4 years ago
Parameter new_aod is not needed anymore as it is assumed to be true all the time. This means that we cannot replay AR4 simulations with new LMDZ sources (we probably couldn't anyway) LF, OB
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: 183.0 KB

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

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

Download in other formats:

Original Format