Context Navigation

physiq_mod.F90 @ 3597

Last change on this file since 3597 was 3597, checked in by Laurent Fairhead, 5 years ago
Temporary fix for tr_ancien initialisation see ticket 104 LF
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 183.0 KB

Line
1	!
2	! $Id: physiq_mod.F90 3597 2019-11-04 15:56:08Z 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	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
341	PARAMETER (ok_cvl=.TRUE.)
342	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
343	PARAMETER (ok_gust=.FALSE.)
344	INTEGER iflag_radia ! active ou non le rayonnement (MPL)
345	SAVE iflag_radia
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	LOGICAL, SAVE :: new_aod
1068	!$OMP THREADPRIVATE(new_aod)
1069	!
1070	!--STRAT AEROSOL
1071	INTEGER, SAVE :: flag_aerosol_strat
1072	!$OMP THREADPRIVATE(flag_aerosol_strat)
1073	!
1074	!--INTERACTIVE AEROSOL FEEDBACK ON RADIATION
1075	LOGICAL, SAVE :: flag_aer_feedback
1076	!$OMP THREADPRIVATE(flag_aer_feedback)
1077
1078	!c-fin STRAT AEROSOL
1079	!
1080	! Declaration des constantes et des fonctions thermodynamiques
1081	!
1082	LOGICAL,SAVE :: first=.TRUE.
1083	!$OMP THREADPRIVATE(first)
1084
1085	! VARIABLES RELATED TO OZONE CLIMATOLOGIES ; all are OpenMP shared
1086	! Note that pressure vectors are in Pa and in stricly ascending order
1087	INTEGER,SAVE :: read_climoz ! Read ozone climatology
1088	! (let it keep the default OpenMP shared attribute)
1089	! Allowed values are 0, 1 and 2
1090	! 0: do not read an ozone climatology
1091	! 1: read a single ozone climatology that will be used day and night
1092	! 2: read two ozone climatologies, the average day and night
1093	! climatology and the daylight climatology
1094	INTEGER,SAVE :: ncid_climoz ! NetCDF file identifier
1095	REAL, POINTER, SAVE :: press_cen_climoz(:) ! Pressure levels
1096	REAL, POINTER, SAVE :: press_edg_climoz(:) ! Edges of pressure intervals
1097	REAL, POINTER, SAVE :: time_climoz(:) ! Time vector
1098	CHARACTER(LEN=13), PARAMETER :: vars_climoz(2) &
1099	= ["tro3 ","tro3_daylight"]
1100	! vars_climoz(1:read_climoz): variables names in climoz file.
1101	! vars_climoz(1:read_climoz-2) if read_climoz>2 (temporary)
1102	REAL :: ro3i ! 0<=ro3i<=360 ; required time index in NetCDF file for
1103	! the ozone fields, old method.
1104
1105	include "YOMCST.h"
1106	include "YOETHF.h"
1107	include "FCTTRE.h"
1108	!IM 100106 BEG : pouvoir sortir les ctes de la physique
1109	include "conema3.h"
1110	include "fisrtilp.h"
1111	include "nuage.h"
1112	include "compbl.h"
1113	!IM 100106 END : pouvoir sortir les ctes de la physique
1114	!
1115	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1116	! Declarations pour Simulateur COSP
1117	!============================================================
1118	real :: mr_ozone(klon,klev), phicosp(klon,klev)
1119
1120	!IM stations CFMIP
1121	INTEGER, SAVE :: nCFMIP
1122	!$OMP THREADPRIVATE(nCFMIP)
1123	INTEGER, PARAMETER :: npCFMIP=120
1124	INTEGER, ALLOCATABLE, SAVE :: tabCFMIP(:)
1125	REAL, ALLOCATABLE, SAVE :: lonCFMIP(:), latCFMIP(:)
1126	!$OMP THREADPRIVATE(tabCFMIP, lonCFMIP, latCFMIP)
1127	INTEGER, ALLOCATABLE, SAVE :: tabijGCM(:)
1128	REAL, ALLOCATABLE, SAVE :: lonGCM(:), latGCM(:)
1129	!$OMP THREADPRIVATE(tabijGCM, lonGCM, latGCM)
1130	INTEGER, ALLOCATABLE, SAVE :: iGCM(:), jGCM(:)
1131	!$OMP THREADPRIVATE(iGCM, jGCM)
1132	logical, dimension(nfiles) :: phys_out_filestations
1133	logical, parameter :: lNMC=.FALSE.
1134
1135	!IM betaCRF
1136	REAL, SAVE :: pfree, beta_pbl, beta_free
1137	!$OMP THREADPRIVATE(pfree, beta_pbl, beta_free)
1138	REAL, SAVE :: lon1_beta, lon2_beta, lat1_beta, lat2_beta
1139	!$OMP THREADPRIVATE(lon1_beta, lon2_beta, lat1_beta, lat2_beta)
1140	LOGICAL, SAVE :: mskocean_beta
1141	!$OMP THREADPRIVATE(mskocean_beta)
1142	REAL, dimension(klon, klev) :: beta ! facteur sur cldtaurad et
1143	! cldemirad pour evaluer les
1144	! retros liees aux CRF
1145	REAL, dimension(klon, klev) :: cldtaurad ! epaisseur optique
1146	! pour radlwsw pour
1147	! tester "CRF off"
1148	REAL, dimension(klon, klev) :: cldtaupirad ! epaisseur optique
1149	! pour radlwsw pour
1150	! tester "CRF off"
1151	REAL, dimension(klon, klev) :: cldemirad ! emissivite pour
1152	! radlwsw pour tester
1153	! "CRF off"
1154	REAL, dimension(klon, klev) :: cldfrarad ! fraction nuageuse
1155
1156	#ifdef INCA
1157	! set de variables utilisees pour l'initialisation des valeurs provenant de INCA
1158	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_tauinca
1159	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_pizinca
1160	REAL, DIMENSION(klon,klev,naero_grp,nbands) :: init_cginca
1161	REAL, DIMENSION(klon,klev,nbands) :: init_ccminca
1162	#endif
1163	REAL, DIMENSION(klon,nbtr) :: init_source
1164
1165	!lwoff=y : offset LW CRE for radiation code and other schemes
1166	REAL, SAVE :: betalwoff
1167	!OMP THREADPRIVATE(betalwoff)
1168	!
1169	INTEGER :: nbtr_tmp ! Number of tracer inside concvl
1170	REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
1171	REAL, dimension(klon,klev) :: ch_in ! Condensed humidity entering in phytrac (eau liquide)
1172	integer iostat
1173
1174	REAL zzz
1175	!albedo SB >>>
1176	REAL,DIMENSION(6), SAVE :: SFRWL
1177	!$OMP THREADPRIVATE(SFRWL)
1178	!albedo SB <<<
1179
1180	!--OB variables for mass fixer (hard coded for now)
1181	LOGICAL, PARAMETER :: mass_fixer=.FALSE.
1182	REAL qql1(klon),qql2(klon),corrqql
1183
1184	REAL pi
1185
1186	pi = 4. * ATAN(1.)
1187
1188	! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
1189	jjmp1=nbp_lat
1190
1191	!======================================================================
1192	! Gestion calendrier : mise a jour du module phys_cal_mod
1193	!
1194	pdtphys=pdtphys_
1195	CALL update_time(pdtphys)
1196	phys_tstep=NINT(pdtphys)
1197	#ifdef CPP_XIOS
1198	IF (.NOT. debut .AND. is_omp_master) CALL xios_update_calendar(itap+1)
1199	#endif
1200
1201	!======================================================================
1202	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1203	! Utilise notamment en 1D mais peut etre active egalement en 3D
1204	! en imposant la valeur de igout.
1205	!======================================================================d
1206	IF (prt_level.ge.1) THEN
1207	igout=klon/2+1/klon
1208	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1209	write(lunout,*) 'igout, lat, lon ',igout, latitude_deg(igout), &
1210	longitude_deg(igout)
1211	write(lunout,*) &
1212	'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1213	write(lunout,*) &
1214	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1215
1216	write(lunout,*) 'paprs, play, phi, u, v, t'
1217	DO k=1,klev
1218	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
1219	u(igout,k),v(igout,k),t(igout,k)
1220	ENDDO
1221	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1222	DO k=1,klev
1223	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1224	ENDDO
1225	ENDIF
1226
1227	! Quick check on pressure levels:
1228	CALL assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
1229	"physiq_mod paprs bad order")
1230
1231	IF (first) THEN
1232	CALL init_etat0_limit_unstruct
1233	IF (.NOT. create_etat0_limit) CALL init_limit_read(days_elapsed)
1234	!CR:nvelles variables convection/poches froides
1235
1236	WRITE(lunout,*) '================================================='
1237	WRITE(lunout,*) 'Allocation des variables locales et sauvegardees'
1238	WRITE(lunout,*) '================================================='
1239	CALL phys_local_var_init
1240	!
1241	! appel a la lecture du run.def physique
1242	CALL conf_phys(ok_journe, ok_mensuel, &
1243	ok_instan, ok_hf, &
1244	ok_LES, &
1245	callstats, &
1246	solarlong0,seuil_inversion, &
1247	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
1248	iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
1249	ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, aerosol_couple, &
1250	chemistry_couple, &
1251	flag_aerosol, flag_aerosol_strat, flag_aer_feedback, new_aod, &
1252	flag_bc_internal_mixture, bl95_b0, bl95_b1, &
1253	! nv flags pour la convection et les
1254	! poches froides
1255	read_climoz, &
1256	alp_offset)
1257	CALL phys_state_var_init(read_climoz)
1258	CALL phys_output_var_init
1259	IF (read_climoz>=1 .AND. create_etat0_limit .AND. grid_type==unstructured) &
1260	CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
1261
1262	#ifdef CPP_StratAer
1263	CALL strataer_init
1264	#endif
1265
1266	print*, '================================================='
1267	!
1268	!CR: check sur le nb de traceurs de l eau
1269	IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
1270	WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
1271	'(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
1272	abort_message='see above'
1273	CALL abort_physic(modname,abort_message,1)
1274	ENDIF
1275
1276	Ncvpaseq1 = 0
1277	dnwd0=0.0
1278	ftd=0.0
1279	fqd=0.0
1280	cin=0.
1281	!ym Attention pbase pas initialise dans concvl !!!!
1282	pbase=0
1283	!IM 180608
1284
1285	itau_con=0
1286	first=.FALSE.
1287
1288	ENDIF ! first
1289
1290	!ym => necessaire pour iflag_con != 2
1291	pmfd(:,:) = 0.
1292	pen_u(:,:) = 0.
1293	pen_d(:,:) = 0.
1294	pde_d(:,:) = 0.
1295	pde_u(:,:) = 0.
1296	aam=0.
1297	d_t_adjwk(:,:)=0
1298	d_q_adjwk(:,:)=0
1299
1300	alp_bl_conv(:)=0.
1301
1302	torsfc=0.
1303	forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1304
1305	modname = 'physiq'
1306
1307	IF (debut) THEN
1308	CALL suphel ! initialiser constantes et parametres phys.
1309	! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1310	tau_gl=5.
1311	CALL getin_p('tau_gl', tau_gl)
1312	! tau_gl : constante de rappel de la temperature a la surface de la glace - en
1313	! secondes
1314	tau_gl=86400.*tau_gl
1315	WRITE(lunout,*) 'debut physiq_mod tau_gl=',tau_gl
1316
1317	iflag_bug_t2m_ipslcm61 = 1
1318	CALL getin_p('iflag_bug_t2m_ipslcm61', iflag_bug_t2m_ipslcm61)
1319	iflag_bug_t2m_stab_ipslcm61 = -1
1320	CALL getin_p('iflag_bug_t2m_stab_ipslcm61', iflag_bug_t2m_stab_ipslcm61)
1321
1322	CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
1323	CALL getin_p('random_notrig_max',random_notrig_max)
1324	CALL getin_p('ok_adjwk',ok_adjwk)
1325	IF (ok_adjwk) iflag_adjwk=2 ! for compatibility with older versions
1326	! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
1327	! 1 => convective adjustment but state variables are unchanged
1328	! 2 => convective adjustment and state variables are changed
1329	CALL getin_p('iflag_adjwk',iflag_adjwk)
1330	CALL getin_p('dtcon_multistep_max',dtcon_multistep_max)
1331	CALL getin_p('dqcon_multistep_max',dqcon_multistep_max)
1332	CALL getin_p('oliqmax',oliqmax)
1333	CALL getin_p('oicemax',oicemax)
1334	CALL getin_p('ratqsp0',ratqsp0)
1335	CALL getin_p('ratqsdp',ratqsdp)
1336	iflag_wake_tend = 0
1337	CALL getin_p('iflag_wake_tend',iflag_wake_tend)
1338	ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
1339	! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
1340	CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
1341	CALL getin_p('ok_bug_cv_trac',ok_bug_cv_trac)
1342	CALL getin_p('ok_bug_split_th',ok_bug_split_th)
1343	fl_ebil = 0 ! by default, conservation diagnostics are desactivated
1344	CALL getin_p('fl_ebil',fl_ebil)
1345	fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
1346	CALL getin_p('fl_cor_ebil',fl_cor_ebil)
1347	iflag_phytrac = 1 ! by default we do want to call phytrac
1348	CALL getin_p('iflag_phytrac',iflag_phytrac)
1349	nvm_lmdz = 13
1350	CALL getin_p('NVM',nvm_lmdz)
1351
1352	WRITE(lunout,*) 'iflag_alp_wk_cond=', iflag_alp_wk_cond
1353	WRITE(lunout,*) 'random_ntrig_max=', random_notrig_max
1354	WRITE(lunout,*) 'ok_adjwk=', ok_adjwk
1355	WRITE(lunout,*) 'iflag_adjwk=', iflag_adjwk
1356	WRITE(lunout,*) 'qtcon_multistep_max=',dtcon_multistep_max
1357	WRITE(lunout,*) 'qdcon_multistep_max=',dqcon_multistep_max
1358	WRITE(lunout,*) 'ratqsp0=', ratqsp0
1359	WRITE(lunout,*) 'ratqsdp=', ratqsdp
1360	WRITE(lunout,*) 'iflag_wake_tend=', iflag_wake_tend
1361	WRITE(lunout,*) 'ok_bad_ecmwf_thermo=',ok_bad_ecmwf_thermo
1362	WRITE(lunout,*) 'ok_bug_cv_trac=', ok_bug_cv_trac
1363	WRITE(lunout,*) 'ok_bug_split_th=', ok_bug_split_th
1364	WRITE(lunout,*) 'fl_ebil=', fl_ebil
1365	WRITE(lunout,*) 'fl_cor_ebil=', fl_cor_ebil
1366	WRITE(lunout,*) 'iflag_phytrac=', iflag_phytrac
1367	WRITE(lunout,*) 'NVM=', nvm_lmdz
1368
1369	!--PC: defining fields to be exchanged between LMDz, ORCHIDEE and NEMO
1370	WRITE(lunout,*) 'Call to infocfields from physiq'
1371	CALL infocfields_init
1372
1373	ENDIF
1374
1375	IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
1376
1377	!======================================================================
1378	! Gestion calendrier : mise a jour du module phys_cal_mod
1379	!
1380	! CALL phys_cal_update(jD_cur,jH_cur)
1381
1382	!
1383	! Si c'est le debut, il faut initialiser plusieurs choses
1384	! ********
1385	!
1386	IF (debut) THEN
1387	!rv CRinitialisation de wght_th et lalim_conv pour la
1388	!definition de la couche alimentation de la convection a partir
1389	!des caracteristiques du thermique
1390	wght_th(:,:)=1.
1391	lalim_conv(:)=1
1392	!RC
1393	ustar(:,:)=0.
1394	! u10m(:,:)=0.
1395	! v10m(:,:)=0.
1396	rain_con(:)=0.
1397	snow_con(:)=0.
1398	topswai(:)=0.
1399	topswad(:)=0.
1400	solswai(:)=0.
1401	solswad(:)=0.
1402
1403	wmax_th(:)=0.
1404	tau_overturning_th(:)=0.
1405
1406	IF (type_trac == 'inca') THEN
1407	! jg : initialisation jusqu'au ces variables sont dans restart
1408	ccm(:,:,:) = 0.
1409	tau_aero(:,:,:,:) = 0.
1410	piz_aero(:,:,:,:) = 0.
1411	cg_aero(:,:,:,:) = 0.
1412
1413	config_inca='none' ! default
1414	CALL getin_p('config_inca',config_inca)
1415
1416	ELSE
1417	config_inca='none' ! default
1418	ENDIF
1419
1420	tau_aero(:,:,:,:) = 1.e-15
1421	piz_aero(:,:,:,:) = 1.
1422	cg_aero(:,:,:,:) = 0.
1423
1424	IF (aerosol_couple .AND. (config_inca /= "aero" &
1425	.AND. config_inca /= "aeNP ")) THEN
1426	abort_message &
1427	= 'if aerosol_couple is activated, config_inca need to be ' &
1428	// 'aero or aeNP'
1429	CALL abort_physic (modname,abort_message,1)
1430	ENDIF
1431
1432	rnebcon0(:,:) = 0.0
1433	clwcon0(:,:) = 0.0
1434	rnebcon(:,:) = 0.0
1435	clwcon(:,:) = 0.0
1436
1437	!
1438	print*,'iflag_coupl,iflag_clos,iflag_wake', &
1439	iflag_coupl,iflag_clos,iflag_wake
1440	print*,'iflag_cycle_diurne', iflag_cycle_diurne
1441	!
1442	IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
1443	abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
1444	CALL abort_physic (modname,abort_message,1)
1445	ENDIF
1446	!
1447	!
1448	! Initialiser les compteurs:
1449	!
1450	itap = 0
1451	itaprad = 0
1452	itapcv = 0
1453	itapwk = 0
1454
1455	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1456	!! Un petit travail \`a faire ici.
1457	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1458
1459	IF (iflag_pbl>1) THEN
1460	PRINT*, "Using method MELLOR&YAMADA"
1461	ENDIF
1462
1463	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1464	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
1465	! phylmd plutot que dyn3d
1466	! Attention : la version precedente n'etait pas tres propre.
1467	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1468	! pour obtenir le meme resultat.
1469	!jyg for fh<
1470	WRITE(lunout,*) 'Pas de temps phys_tstep pdtphys ',phys_tstep,pdtphys
1471	IF (abs(phys_tstep-pdtphys)>1.e-10) THEN
1472	abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
1473	CALL abort_physic(modname,abort_message,1)
1474	ENDIF
1475	!>jyg
1476	IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
1477	radpas = NINT( 86400./phys_tstep)/nbapp_rad
1478	ELSE
1479	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1480	'multiple de nbapp_rad'
1481	WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
1482	'mais 1+1<>2'
1483	abort_message='nbre de pas de temps physique n est pas multiple ' &
1484	// 'de nbapp_rad'
1485	CALL abort_physic(modname,abort_message,1)
1486	ENDIF
1487	IF (nbapp_cv .EQ. 0) nbapp_cv=86400./phys_tstep
1488	IF (nbapp_wk .EQ. 0) nbapp_wk=86400./phys_tstep
1489	print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
1490	IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
1491	cvpas_0 = NINT( 86400./phys_tstep)/nbapp_cv
1492	cvpas = cvpas_0
1493	print *,'physiq, cvpas ',cvpas
1494	ELSE
1495	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1496	'multiple de nbapp_cv'
1497	WRITE(lunout,*) 'changer nbapp_cv ou alors commenter ce test ', &
1498	'mais 1+1<>2'
1499	abort_message='nbre de pas de temps physique n est pas multiple ' &
1500	// 'de nbapp_cv'
1501	CALL abort_physic(modname,abort_message,1)
1502	ENDIF
1503	IF (MOD(NINT(86400./phys_tstep),nbapp_wk).EQ.0) THEN
1504	wkpas = NINT( 86400./phys_tstep)/nbapp_wk
1505	! print *,'physiq, wkpas ',wkpas
1506	ELSE
1507	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1508	'multiple de nbapp_wk'
1509	WRITE(lunout,*) 'changer nbapp_wk ou alors commenter ce test ', &
1510	'mais 1+1<>2'
1511	abort_message='nbre de pas de temps physique n est pas multiple ' &
1512	// 'de nbapp_wk'
1513	CALL abort_physic(modname,abort_message,1)
1514	ENDIF
1515	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1516	CALL init_iophy_new(latitude_deg,longitude_deg)
1517
1518	!===================================================================
1519	!IM stations CFMIP
1520	nCFMIP=npCFMIP
1521	OPEN(98,file='npCFMIP_param.data',status='old', &
1522	form='formatted',iostat=iostat)
1523	IF (iostat == 0) THEN
1524	READ(98,*,end=998) nCFMIP
1525	998 CONTINUE
1526	CLOSE(98)
1527	CONTINUE
1528	IF(nCFMIP.GT.npCFMIP) THEN
1529	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1530	CALL abort_physic("physiq", "", 1)
1531	ELSE
1532	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1533	ENDIF
1534
1535	!
1536	ALLOCATE(tabCFMIP(nCFMIP))
1537	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1538	ALLOCATE(tabijGCM(nCFMIP))
1539	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1540	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1541	!
1542	! lecture des nCFMIP stations CFMIP, de leur numero
1543	! et des coordonnees geographiques lonCFMIP, latCFMIP
1544	!
1545	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
1546	lonCFMIP, latCFMIP)
1547	!
1548	! identification des
1549	! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
1550	! grille de LMDZ
1551	! 2) indices points tabijGCM de la grille physique 1d sur
1552	! klon points
1553	! 3) indices iGCM, jGCM de la grille physique 2d
1554	!
1555	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1556	tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1557	!
1558	ELSE
1559	ALLOCATE(tabijGCM(0))
1560	ALLOCATE(lonGCM(0), latGCM(0))
1561	ALLOCATE(iGCM(0), jGCM(0))
1562	ENDIF
1563
1564	#ifdef CPP_IOIPSL
1565
1566	!$OMP MASTER
1567	! FH : if ok_sync=.true. , the time axis is written at each time step
1568	! in the output files. Only at the end in the opposite case
1569	ok_sync_omp=.FALSE.
1570	CALL getin('ok_sync',ok_sync_omp)
1571	CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
1572	iGCM,jGCM,lonGCM,latGCM, &
1573	jjmp1,nlevSTD,clevSTD,rlevSTD, phys_tstep,ok_veget, &
1574	type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1575	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1576	read_climoz, phys_out_filestations, &
1577	new_aod, aerosol_couple, &
1578	flag_aerosol_strat, pdtphys, paprs, pphis, &
1579	pplay, lmax_th, ptconv, ptconvth, ivap, &
1580	d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
1581	!$OMP END MASTER
1582	!$OMP BARRIER
1583	ok_sync=ok_sync_omp
1584
1585	freq_outNMC(1) = ecrit_files(7)
1586	freq_outNMC(2) = ecrit_files(8)
1587	freq_outNMC(3) = ecrit_files(9)
1588	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1589	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1590	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1591
1592	#ifndef CPP_XIOS
1593	CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
1594	#endif
1595
1596	#endif
1597	ecrit_reg = ecrit_reg * un_jour
1598	ecrit_tra = ecrit_tra * un_jour
1599
1600	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
1601	date0 = jD_ref
1602	WRITE(,) 'physiq date0 : ',date0
1603	!
1604
1605	! CALL create_climoz(read_climoz)
1606	IF (.NOT. create_etat0_limit) CALL init_aero_fromfile(flag_aerosol) !! initialise aero from file for XIOS interpolation (unstructured_grid)
1607	IF (.NOT. create_etat0_limit) CALL init_readaerosolstrato(flag_aerosol_strat) !! initialise aero strato from file for XIOS interpolation (unstructured_grid)
1608
1609	#ifdef CPP_COSP
1610	IF (ok_cosp) THEN
1611	DO k = 1, klev
1612	DO i = 1, klon
1613	phicosp(i,k) = pphi(i,k) + pphis(i)
1614	ENDDO
1615	ENDDO
1616	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
1617	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1618	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1619	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1620	JrNt,ref_liq,ref_ice, &
1621	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1622	zu10m,zv10m,pphis, &
1623	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1624	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1625	prfl(:,1:klev),psfl(:,1:klev), &
1626	pmflxr(:,1:klev),pmflxs(:,1:klev), &
1627	mr_ozone,cldtau, cldemi)
1628	ENDIF
1629	#endif
1630
1631	#ifdef CPP_COSP2
1632	IF (ok_cosp) THEN
1633	DO k = 1, klev
1634	DO i = 1, klon
1635	phicosp(i,k) = pphi(i,k) + pphis(i)
1636	ENDDO
1637	ENDDO
1638	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
1639	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1640	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1641	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1642	JrNt,ref_liq,ref_ice, &
1643	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1644	zu10m,zv10m,pphis, &
1645	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1646	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1647	prfl(:,1:klev),psfl(:,1:klev), &
1648	pmflxr(:,1:klev),pmflxs(:,1:klev), &
1649	mr_ozone,cldtau, cldemi)
1650	ENDIF
1651	#endif
1652
1653	#ifdef CPP_COSPV2
1654	IF (ok_cosp) THEN
1655	DO k = 1, klev
1656	DO i = 1, klon
1657	phicosp(i,k) = pphi(i,k) + pphis(i)
1658	ENDDO
1659	ENDDO
1660	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
1661	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
1662	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
1663	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
1664	JrNt,ref_liq,ref_ice, &
1665	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
1666	zu10m,zv10m,pphis, &
1667	phicosp,paprs(:,1:klev),pplay,zxtsol,t_seri, &
1668	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
1669	prfl(:,1:klev),psfl(:,1:klev), &
1670	pmflxr(:,1:klev),pmflxs(:,1:klev), &
1671	mr_ozone,cldtau, cldemi)
1672	ENDIF
1673	#endif
1674
1675	!
1676	!
1677	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1678	! Nouvelle initialisation pour le rayonnement RRTM
1679	!
1680	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1681
1682	CALL iniradia(klon,klev,paprs(1,1:klev+1))
1683	! Initialisation des champs dans phytrac qui sont utilisés par phys_output_write
1684	IF (iflag_phytrac == 1 ) THEN
1685	CALL phytrac_init()
1686	ENDIF
1687
1688	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
1689	pplay, lmax_th, aerosol_couple, &
1690	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod, ok_sync,&
1691	ptconv, read_climoz, clevSTD, &
1692	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
1693	flag_aerosol, flag_aerosol_strat, ok_cdnc)
1694
1695	#ifdef CPP_XIOS
1696	IF (is_omp_master) CALL xios_update_calendar(1)
1697	#endif
1698	IF(read_climoz>=1 .AND. create_etat0_limit) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
1699	CALL create_etat0_limit_unstruct
1700	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1701
1702	!jyg<
1703	IF (iflag_pbl<=1) THEN
1704	! No TKE for Standard Physics
1705	pbl_tke(:,:,:)=0.
1706
1707	ELSE IF (klon_glo==1) THEN
1708	pbl_tke(:,:,is_ave) = 0.
1709	DO nsrf=1,nbsrf
1710	DO k = 1,klev+1
1711	pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
1712	+pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
1713	ENDDO
1714	ENDDO
1715	ELSE
1716	pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ??
1717	!>jyg
1718	ENDIF
1719	!IM begin
1720	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1721	,ratqs(1,1)
1722	!IM end
1723
1724
1725	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1726	!
1727	! on remet le calendrier a zero
1728	!
1729	IF (raz_date .eq. 1) THEN
1730	itau_phy = 0
1731	ENDIF
1732
1733	! IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
1734	! WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
1735	! pdtphys
1736	! abort_message='Pas physique n est pas correct '
1737	! ! call abort_physic(modname,abort_message,1)
1738	! phys_tstep=pdtphys
1739	! ENDIF
1740	IF (nlon .NE. klon) THEN
1741	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1742	klon
1743	abort_message='nlon et klon ne sont pas coherents'
1744	CALL abort_physic(modname,abort_message,1)
1745	ENDIF
1746	IF (nlev .NE. klev) THEN
1747	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1748	klev
1749	abort_message='nlev et klev ne sont pas coherents'
1750	CALL abort_physic(modname,abort_message,1)
1751	ENDIF
1752	!
1753	IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
1754	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1755	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1756	abort_message='Nbre d appels au rayonnement insuffisant'
1757	CALL abort_physic(modname,abort_message,1)
1758	ENDIF
1759	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1760	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1761	ok_cvl
1762	!
1763	!KE43
1764	! Initialisation pour la convection de K.E. (sb):
1765	IF (iflag_con.GE.3) THEN
1766
1767	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1768	WRITE(lunout,*) &
1769	"On va utiliser le melange convectif des traceurs qui"
1770	WRITE(lunout,*)"est calcule dans convect4.3"
1771	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1772
1773	DO i = 1, klon
1774	ema_cbmf(i) = 0.
1775	ema_pcb(i) = 0.
1776	ema_pct(i) = 0.
1777	! ema_workcbmf(i) = 0.
1778	ENDDO
1779	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1780	DO i = 1, klon
1781	ibas_con(i) = 1
1782	itop_con(i) = 1
1783	ENDDO
1784	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1785	!================================================================
1786	!CR:04.12.07: initialisations poches froides
1787	! Controle de ALE et ALP pour la fermeture convective (jyg)
1788	IF (iflag_wake>=1) THEN
1789	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1790	,alp_bl_prescr, ale_bl_prescr)
1791	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1792	! print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
1793	!
1794	! Initialize tendencies of wake state variables (for some flag values
1795	! they are not computed).
1796	d_deltat_wk(:,:) = 0.
1797	d_deltaq_wk(:,:) = 0.
1798	d_deltat_wk_gw(:,:) = 0.
1799	d_deltaq_wk_gw(:,:) = 0.
1800	d_deltat_vdf(:,:) = 0.
1801	d_deltaq_vdf(:,:) = 0.
1802	d_deltat_the(:,:) = 0.
1803	d_deltaq_the(:,:) = 0.
1804	d_deltat_ajs_cv(:,:) = 0.
1805	d_deltaq_ajs_cv(:,:) = 0.
1806	d_s_wk(:) = 0.
1807	d_dens_wk(:) = 0.
1808	ENDIF
1809
1810	! do i = 1,klon
1811	! Ale_bl(i)=0.
1812	! Alp_bl(i)=0.
1813	! enddo
1814
1815	!ELSE
1816	! ALLOCATE(tabijGCM(0))
1817	! ALLOCATE(lonGCM(0), latGCM(0))
1818	! ALLOCATE(iGCM(0), jGCM(0))
1819	ENDIF
1820
1821	DO i=1,klon
1822	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1823	ENDDO
1824
1825	!34EK
1826	IF (ok_orodr) THEN
1827
1828	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1829	! FH sans doute a enlever de finitivement ou, si on le
1830	! garde, l'activer justement quand ok_orodr = false.
1831	! ce rugoro est utilise par la couche limite et fait double emploi
1832	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1833	! DO i=1,klon
1834	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1835	! ENDDO
1836	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1837	IF (ok_strato) THEN
1838	CALL SUGWD_strato(klon,klev,paprs,pplay)
1839	ELSE
1840	CALL SUGWD(klon,klev,paprs,pplay)
1841	ENDIF
1842
1843	DO i=1,klon
1844	zuthe(i)=0.
1845	zvthe(i)=0.
1846	IF (zstd(i).gt.10.) THEN
1847	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1848	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1849	ENDIF
1850	ENDDO
1851	ENDIF
1852	!
1853	!
1854	lmt_pas = NINT(86400./phys_tstep * 1.0) ! tous les jours
1855	WRITE(lunout,*)'La frequence de lecture surface est de ', &
1856	lmt_pas
1857	!
1858	capemaxcels = 't_max(X)'
1859	t2mincels = 't_min(X)'
1860	t2maxcels = 't_max(X)'
1861	tinst = 'inst(X)'
1862	tave = 'ave(X)'
1863	!IM cf. AM 081204 BEG
1864	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1865	!IM cf. AM 081204 END
1866	!
1867	!=============================================================
1868	! Initialisation des sorties
1869	!=============================================================
1870
1871	#ifdef CPP_XIOS
1872	! Get "missing_val" value from XML files (from temperature variable)
1873	!$OMP MASTER
1874	CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1875	!$OMP END MASTER
1876	!$OMP BARRIER
1877	missing_val=missing_val_omp
1878	#endif
1879
1880	#ifdef CPP_XIOS
1881	! Need to put this initialisation after phyetat0 as in the coupled model the XIOS context is only
1882	! initialised at that moment
1883	! Get "missing_val" value from XML files (from temperature variable)
1884	!$OMP MASTER
1885	CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1886	!$OMP END MASTER
1887	!$OMP BARRIER
1888	missing_val=missing_val_omp
1889	#endif
1890
1891
1892	CALL printflag( tabcntr0,radpas,ok_journe, &
1893	ok_instan, ok_region )
1894	!
1895	!
1896	!
1897	! Prescrire l'ozone dans l'atmosphere
1898	!
1899	!
1900	!c DO i = 1, klon
1901	!c DO k = 1, klev
1902	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1903	!c ENDDO
1904	!c ENDDO
1905	!
1906	IF (type_trac == 'inca') THEN
1907	#ifdef INCA
1908	CALL VTe(VTphysiq)
1909	CALL VTb(VTinca)
1910	calday = REAL(days_elapsed) + jH_cur
1911	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1912
1913	CALL chemini( &
1914	rg, &
1915	ra, &
1916	cell_area, &
1917	latitude_deg, &
1918	longitude_deg, &
1919	presnivs, &
1920	calday, &
1921	klon, &
1922	nqtot, &
1923	nqo, &
1924	pdtphys, &
1925	annee_ref, &
1926	year_cur, &
1927	day_ref, &
1928	day_ini, &
1929	start_time, &
1930	itau_phy, &
1931	date0, &
1932	io_lon, &
1933	io_lat, &
1934	chemistry_couple, &
1935	init_source, &
1936	init_tauinca, &
1937	init_pizinca, &
1938	init_cginca, &
1939	init_ccminca)
1940
1941
1942	! initialisation des variables depuis le restart de inca
1943	ccm(:,:,:) = init_ccminca
1944	tau_aero(:,:,:,:) = init_tauinca
1945	piz_aero(:,:,:,:) = init_pizinca
1946	cg_aero(:,:,:,:) = init_cginca
1947	!
1948
1949
1950	CALL VTe(VTinca)
1951	CALL VTb(VTphysiq)
1952	#endif
1953	ENDIF
1954
1955	!$omp single
1956	IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, &
1957	press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
1958	!$omp end single
1959	!
1960	!IM betaCRF
1961	pfree=70000. !Pa
1962	beta_pbl=1.
1963	beta_free=1.
1964	lon1_beta=-180.
1965	lon2_beta=+180.
1966	lat1_beta=90.
1967	lat2_beta=-90.
1968	mskocean_beta=.FALSE.
1969
1970	!albedo SB >>>
1971	SELECT CASE(nsw)
1972	CASE(2)
1973	SFRWL(1)=0.45538747
1974	SFRWL(2)=0.54461211
1975	CASE(4)
1976	SFRWL(1)=0.45538747
1977	SFRWL(2)=0.32870591
1978	SFRWL(3)=0.18568763
1979	SFRWL(4)=3.02191470E-02
1980	CASE(6)
1981	SFRWL(1)=1.28432794E-03
1982	SFRWL(2)=0.12304168
1983	SFRWL(3)=0.33106142
1984	SFRWL(4)=0.32870591
1985	SFRWL(5)=0.18568763
1986	SFRWL(6)=3.02191470E-02
1987	END SELECT
1988
1989
1990	!albedo SB <<<
1991
1992	OPEN(99,file='beta_crf.data',status='old', &
1993	form='formatted',err=9999)
1994	READ(99,*,end=9998) pfree
1995	READ(99,*,end=9998) beta_pbl
1996	READ(99,*,end=9998) beta_free
1997	READ(99,*,end=9998) lon1_beta
1998	READ(99,*,end=9998) lon2_beta
1999	READ(99,*,end=9998) lat1_beta
2000	READ(99,*,end=9998) lat2_beta
2001	READ(99,*,end=9998) mskocean_beta
2002	9998 Continue
2003	CLOSE(99)
2004	9999 Continue
2005	WRITE(,)'pfree=',pfree
2006	WRITE(,)'beta_pbl=',beta_pbl
2007	WRITE(,)'beta_free=',beta_free
2008	WRITE(,)'lon1_beta=',lon1_beta
2009	WRITE(,)'lon2_beta=',lon2_beta
2010	WRITE(,)'lat1_beta=',lat1_beta
2011	WRITE(,)'lat2_beta=',lat2_beta
2012	WRITE(,)'mskocean_beta=',mskocean_beta
2013
2014	!lwoff=y : offset LW CRE for radiation code and other schemes
2015	!lwoff=y : betalwoff=1.
2016	betalwoff=0.
2017	IF (ok_lwoff) THEN
2018	betalwoff=1.
2019	ENDIF
2020	WRITE(,)'ok_lwoff=',ok_lwoff
2021	!
2022	!lwoff=y to begin only sollw and sollwdown are set up to CS values
2023	sollw = sollw + betalwoff * (sollw0 - sollw)
2024	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
2025	sollwdown(:))
2026	!
2027	! Temporary solutions adressing ticket #104 and the non initialisation of tr_ancien
2028	! LF
2029	WRITE(lunout,*)' WARNING: tr_ancien initialised to tr_seri'
2030	DO iq = nqo+1, nqtot
2031	tr_ancien(:,:,iq-nqo)=tr_seri(:,:,iq-nqo)
2032	ENDDO
2033
2034
2035	ENDIF
2036	!
2037	! ************** Fin de IF ( debut ) *************
2038	!
2039	!
2040	! Incrementer le compteur de la physique
2041	!
2042	itap = itap + 1
2043	IF (is_master .OR. prt_level > 9) THEN
2044	IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
2045	WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
2046	WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
2047	100 FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
2048	ENDIF
2049	ENDIF
2050	!
2051	!
2052	! Update fraction of the sub-surfaces (pctsrf) and
2053	! initialize, where a new fraction has appeared, all variables depending
2054	! on the surface fraction.
2055	!
2056	CALL change_srf_frac(itap, phys_tstep, days_elapsed+1, &
2057	pctsrf, fevap, z0m, z0h, agesno, &
2058	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
2059
2060	! Update time and other variables in Reprobus
2061	IF (type_trac == 'repr') THEN
2062	#ifdef REPROBUS
2063	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
2064	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
2065	CALL Rtime(debut)
2066	#endif
2067	ENDIF
2068
2069	! Tendances bidons pour les processus qui n'affectent pas certaines
2070	! variables.
2071	du0(:,:)=0.
2072	dv0(:,:)=0.
2073	dt0 = 0.
2074	dq0(:,:)=0.
2075	dql0(:,:)=0.
2076	dqi0(:,:)=0.
2077	dsig0(:) = 0.
2078	ddens0(:) = 0.
2079	wkoccur1(:)=1
2080	!
2081	! Mettre a zero des variables de sortie (pour securite)
2082	!
2083	DO i = 1, klon
2084	d_ps(i) = 0.0
2085	ENDDO
2086	DO k = 1, klev
2087	DO i = 1, klon
2088	d_t(i,k) = 0.0
2089	d_u(i,k) = 0.0
2090	d_v(i,k) = 0.0
2091	ENDDO
2092	ENDDO
2093	DO iq = 1, nqtot
2094	DO k = 1, klev
2095	DO i = 1, klon
2096	d_qx(i,k,iq) = 0.0
2097	ENDDO
2098	ENDDO
2099	ENDDO
2100	beta_prec_fisrt(:,:)=0.
2101	beta_prec(:,:)=0.
2102	!
2103	! Output variables from the convective scheme should not be set to 0
2104	! since convection is not always called at every time step.
2105	IF (ok_bug_cv_trac) THEN
2106	da(:,:)=0.
2107	mp(:,:)=0.
2108	phi(:,:,:)=0.
2109	! RomP >>>
2110	phi2(:,:,:)=0.
2111	epmlmMm(:,:,:)=0.
2112	eplaMm(:,:)=0.
2113	d1a(:,:)=0.
2114	dam(:,:)=0.
2115	pmflxr(:,:)=0.
2116	pmflxs(:,:)=0.
2117	! RomP <<<
2118	ENDIF
2119
2120	!
2121	! Ne pas affecter les valeurs entrees de u, v, h, et q
2122	!
2123	DO k = 1, klev
2124	DO i = 1, klon
2125	t_seri(i,k) = t(i,k)
2126	u_seri(i,k) = u(i,k)
2127	v_seri(i,k) = v(i,k)
2128	q_seri(i,k) = qx(i,k,ivap)
2129	ql_seri(i,k) = qx(i,k,iliq)
2130	!CR: ATTENTION, on rajoute la variable glace
2131	IF (nqo.eq.2) THEN
2132	qs_seri(i,k) = 0.
2133	ELSE IF (nqo.eq.3) THEN
2134	qs_seri(i,k) = qx(i,k,isol)
2135	ENDIF
2136	ENDDO
2137	ENDDO
2138	!
2139	!--OB mass fixer
2140	IF (mass_fixer) THEN
2141	!--store initial water burden
2142	qql1(:)=0.0
2143	DO k = 1, klev
2144	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
2145	ENDDO
2146	ENDIF
2147	!--fin mass fixer
2148
2149	tke0(:,:)=pbl_tke(:,:,is_ave)
2150	!CR:Nombre de traceurs de l'eau: nqo
2151	! IF (nqtot.GE.3) THEN
2152	IF (nqtot.GE.(nqo+1)) THEN
2153	! DO iq = 3, nqtot
2154	DO iq = nqo+1, nqtot
2155	DO k = 1, klev
2156	DO i = 1, klon
2157	! tr_seri(i,k,iq-2) = qx(i,k,iq)
2158	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
2159	ENDDO
2160	ENDDO
2161	ENDDO
2162	ELSE
2163	DO k = 1, klev
2164	DO i = 1, klon
2165	tr_seri(i,k,1) = 0.0
2166	ENDDO
2167	ENDDO
2168	ENDIF
2169	!
2170	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	wo(:, :, 1), &
3729	q_seri, &
3730	zxtsol, &
3731	zxsnow_dummy, &
3732	solsw, &
3733	albsol1, &
3734	rain_fall, &
3735	snow_fall, &
3736	itop_con, &
3737	ibas_con, &
3738	cldfra, &
3739	nbp_lon, &
3740	nbp_lat-1, &
3741	tr_seri, &
3742	ftsol, &
3743	paprs, &
3744	cdragh, &
3745	cdragm, &
3746	pctsrf, &
3747	pdtphys, &
3748	itap)
3749
3750	CALL VTe(VTinca)
3751	CALL VTb(VTphysiq)
3752	#endif
3753	ENDIF !type_trac = inca
3754
3755
3756	!
3757	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3758	!
3759	IF (MOD(itaprad,radpas).EQ.0) THEN
3760
3761	!
3762	!jq - introduce the aerosol direct and first indirect radiative forcings
3763	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3764	IF (flag_aerosol .GT. 0) THEN
3765	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3766	IF (.NOT. aerosol_couple) THEN
3767	!
3768	CALL readaerosol_optic( &
3769	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3770	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3771	mass_solu_aero, mass_solu_aero_pi, &
3772	tau_aero, piz_aero, cg_aero, &
3773	tausum_aero, tau3d_aero)
3774	ENDIF
3775	ELSE ! RRTM radiation
3776	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3777	abort_message='config_inca=aero et rrtm=1 impossible'
3778	CALL abort_physic(modname,abort_message,1)
3779	ELSE
3780	!
3781	#ifdef CPP_RRTM
3782	IF (NSW.EQ.6) THEN
3783	!--new aerosol properties SW and LW
3784	!
3785	#ifdef CPP_Dust
3786	!--SPL aerosol model
3787	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
3788	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3789	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3790	tausum_aero, tau3d_aero)
3791	#else
3792	!--climatologies or INCA aerosols
3793	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, ok_volcan, &
3794	new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
3795	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3796	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3797	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3798	tausum_aero, drytausum_aero, tau3d_aero)
3799	#endif
3800
3801	IF (flag_aerosol .EQ. 7) THEN
3802	CALL MACv2SP(pphis,pplay,paprs,longitude_deg,latitude_deg, &
3803	tau_aero_sw_rrtm,piz_aero_sw_rrtm,cg_aero_sw_rrtm)
3804	ENDIF
3805
3806	!
3807	ELSE IF (NSW.EQ.2) THEN
3808	!--for now we use the old aerosol properties
3809	!
3810	CALL readaerosol_optic( &
3811	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3812	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3813	mass_solu_aero, mass_solu_aero_pi, &
3814	tau_aero, piz_aero, cg_aero, &
3815	tausum_aero, tau3d_aero)
3816	!
3817	!--natural aerosols
3818	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
3819	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
3820	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
3821	!--all aerosols
3822	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
3823	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
3824	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
3825	!
3826	!--no LW optics
3827	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3828	!
3829	ELSE
3830	abort_message='Only NSW=2 or 6 are possible with ' &
3831	// 'aerosols and iflag_rrtm=1'
3832	CALL abort_physic(modname,abort_message,1)
3833	ENDIF
3834	#else
3835	abort_message='You should compile with -rrtm if running ' &
3836	// 'with iflag_rrtm=1'
3837	CALL abort_physic(modname,abort_message,1)
3838	#endif
3839	!
3840	ENDIF
3841	ENDIF
3842	ELSE !--flag_aerosol = 0
3843	tausum_aero(:,:,:) = 0.
3844	drytausum_aero(:,:) = 0.
3845	mass_solu_aero(:,:) = 0.
3846	mass_solu_aero_pi(:,:) = 0.
3847	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3848	tau_aero(:,:,:,:) = 1.e-15
3849	piz_aero(:,:,:,:) = 1.
3850	cg_aero(:,:,:,:) = 0.
3851	ELSE
3852	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
3853	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3854	piz_aero_sw_rrtm(:,:,:,:) = 1.0
3855	cg_aero_sw_rrtm(:,:,:,:) = 0.0
3856	ENDIF
3857	ENDIF
3858	!
3859	!--WMO criterion to determine tropopause
3860	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3861	!
3862	!--STRAT AEROSOL
3863	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
3864	IF (flag_aerosol_strat.GT.0) THEN
3865	IF (prt_level .GE.10) THEN
3866	PRINT *,'appel a readaerosolstrat', mth_cur
3867	ENDIF
3868	IF (iflag_rrtm.EQ.0) THEN
3869	IF (flag_aerosol_strat.EQ.1) THEN
3870	CALL readaerosolstrato(debut)
3871	ELSE
3872	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
3873	CALL abort_physic(modname,abort_message,1)
3874	ENDIF
3875	ELSE
3876	#ifdef CPP_RRTM
3877	#ifndef CPP_StratAer
3878	!--prescribed strat aerosols
3879	!--only in the case of non-interactive strat aerosols
3880	IF (flag_aerosol_strat.EQ.1) THEN
3881	CALL readaerosolstrato1_rrtm(debut)
3882	ELSEIF (flag_aerosol_strat.EQ.2) THEN
3883	CALL readaerosolstrato2_rrtm(debut, ok_volcan)
3884	ELSE
3885	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
3886	CALL abort_physic(modname,abort_message,1)
3887	ENDIF
3888	#endif
3889	#else
3890	abort_message='You should compile with -rrtm if running ' &
3891	// 'with iflag_rrtm=1'
3892	CALL abort_physic(modname,abort_message,1)
3893	#endif
3894	ENDIF
3895	ELSE
3896	tausum_aero(:,:,id_STRAT_phy) = 0.
3897	ENDIF
3898	!
3899	#ifdef CPP_RRTM
3900	#ifdef CPP_StratAer
3901	!--compute stratospheric mask
3902	CALL stratosphere_mask(missing_val, pphis, t_seri, pplay, latitude_deg)
3903	!--interactive strat aerosols
3904	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
3905	#endif
3906	#endif
3907	!--fin STRAT AEROSOL
3908	!
3909
3910	! Calculer les parametres optiques des nuages et quelques
3911	! parametres pour diagnostiques:
3912	!
3913	IF (aerosol_couple.AND.config_inca=='aero') THEN
3914	mass_solu_aero(:,:) = ccm(:,:,1)
3915	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3916	ENDIF
3917
3918	IF (ok_newmicro) then
3919	IF (iflag_rrtm.NE.0) THEN
3920	#ifdef CPP_RRTM
3921	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3922	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
3923	// 'pour ok_cdnc'
3924	CALL abort_physic(modname,abort_message,1)
3925	ENDIF
3926	#else
3927
3928	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
3929	CALL abort_physic(modname,abort_message,1)
3930	#endif
3931	ENDIF
3932	CALL newmicro (flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, &
3933	paprs, pplay, t_seri, cldliq, cldfra, &
3934	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3935	flwp, fiwp, flwc, fiwc, &
3936	mass_solu_aero, mass_solu_aero_pi, &
3937	cldtaupi, re, fl, ref_liq, ref_ice, &
3938	ref_liq_pi, ref_ice_pi)
3939	ELSE
3940	CALL nuage (paprs, pplay, &
3941	t_seri, cldliq, cldfra, cldtau, cldemi, &
3942	cldh, cldl, cldm, cldt, cldq, &
3943	ok_aie, &
3944	mass_solu_aero, mass_solu_aero_pi, &
3945	bl95_b0, bl95_b1, &
3946	cldtaupi, re, fl)
3947	ENDIF
3948	!
3949	!IM betaCRF
3950	!
3951	cldtaurad = cldtau
3952	cldtaupirad = cldtaupi
3953	cldemirad = cldemi
3954	cldfrarad = cldfra
3955
3956	!
3957	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3958	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3959	!
3960	! global
3961	!
3962	!IM 251017 begin
3963	! print*,'physiq betaCRF global zdtime=',zdtime
3964	!IM 251017 end
3965	DO k=1, klev
3966	DO i=1, klon
3967	IF (pplay(i,k).GE.pfree) THEN
3968	beta(i,k) = beta_pbl
3969	ELSE
3970	beta(i,k) = beta_free
3971	ENDIF
3972	IF (mskocean_beta) THEN
3973	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3974	ENDIF
3975	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3976	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3977	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3978	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3979	ENDDO
3980	ENDDO
3981	!
3982	ELSE
3983	!
3984	! regional
3985	!
3986	DO k=1, klev
3987	DO i=1,klon
3988	!
3989	IF (longitude_deg(i).ge.lon1_beta.AND. &
3990	longitude_deg(i).le.lon2_beta.AND. &
3991	latitude_deg(i).le.lat1_beta.AND. &
3992	latitude_deg(i).ge.lat2_beta) THEN
3993	IF (pplay(i,k).GE.pfree) THEN
3994	beta(i,k) = beta_pbl
3995	ELSE
3996	beta(i,k) = beta_free
3997	ENDIF
3998	IF (mskocean_beta) THEN
3999	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
4000	ENDIF
4001	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
4002	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
4003	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
4004	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
4005	ENDIF
4006	!
4007	ENDDO
4008	ENDDO
4009	!
4010	ENDIF
4011
4012	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
4013	IF (ok_chlorophyll) THEN
4014	print*,"-- reading chlorophyll"
4015	CALL readchlorophyll(debut)
4016	ENDIF
4017
4018	!--if ok_suntime_rrtm we use ancillay data for RSUN
4019	!--previous values are therefore overwritten
4020	!--this is needed for CMIP6 runs
4021	!--and only possible for new radiation scheme
4022	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
4023	#ifdef CPP_RRTM
4024	CALL read_rsun_rrtm(debut)
4025	#endif
4026	ENDIF
4027
4028	IF (mydebug) THEN
4029	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4030	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4031	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4032	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4033	ENDIF
4034
4035	!
4036	!sonia : If Iflag_radia >=2, pertubation of some variables
4037	!input to radiation (DICE)
4038	!
4039	IF (iflag_radia .ge. 2) THEN
4040	zsav_tsol (:) = zxtsol(:)
4041	CALL perturb_radlwsw(zxtsol,iflag_radia)
4042	ENDIF
4043
4044	IF (aerosol_couple.AND.config_inca=='aero') THEN
4045	#ifdef INCA
4046	CALL radlwsw_inca &
4047	(chemistry_couple, kdlon,kflev,dist, rmu0, fract, solaire, &
4048	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
4049	size(wo,3), wo, &
4050	cldfrarad, cldemirad, cldtaurad, &
4051	heat,heat0,cool,cool0,albpla, &
4052	topsw,toplw,solsw,sollw, &
4053	sollwdown, &
4054	topsw0,toplw0,solsw0,sollw0, &
4055	lwdn0, lwdn, lwup0, lwup, &
4056	swdn0, swdn, swup0, swup, &
4057	ok_ade, ok_aie, &
4058	tau_aero, piz_aero, cg_aero, &
4059	topswad_aero, solswad_aero, &
4060	topswad0_aero, solswad0_aero, &
4061	topsw_aero, topsw0_aero, &
4062	solsw_aero, solsw0_aero, &
4063	cldtaupirad, &
4064	topswai_aero, solswai_aero)
4065	#endif
4066	ELSE
4067	!
4068	!IM calcul radiatif pour le cas actuel
4069	!
4070	RCO2 = RCO2_act
4071	RCH4 = RCH4_act
4072	RN2O = RN2O_act
4073	RCFC11 = RCFC11_act
4074	RCFC12 = RCFC12_act
4075	!
4076	!--interactive CO2 in ppm from carbon cycle
4077	IF (carbon_cycle_rad.AND..NOT.debut) THEN
4078	RCO2=RCO2_glo
4079	ENDIF
4080	!
4081	IF (prt_level .GE.10) THEN
4082	print *,' ->radlwsw, number 1 '
4083	ENDIF
4084	!
4085	CALL radlwsw &
4086	(dist, rmu0, fract, &
4087	!albedo SB >>>
4088	! paprs, pplay,zxtsol,albsol1, albsol2, &
4089	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4090	!albedo SB <<<
4091	t_seri,q_seri,wo, &
4092	cldfrarad, cldemirad, cldtaurad, &
4093	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, &
4094	flag_aerosol, &
4095	flag_aerosol_strat, flag_aer_feedback, &
4096	tau_aero, piz_aero, cg_aero, &
4097	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4098	! Rajoute par OB pour RRTM
4099	tau_aero_lw_rrtm, &
4100	cldtaupirad,new_aod, &
4101	! zqsat, flwcrad, fiwcrad, &
4102	zqsat, flwc, fiwc, &
4103	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4104	heat,heat0,cool,cool0,albpla, &
4105	heat_volc,cool_volc, &
4106	topsw,toplw,solsw,sollw, &
4107	sollwdown, &
4108	topsw0,toplw0,solsw0,sollw0, &
4109	lwdnc0, lwdn0, lwdn, lwupc0, lwup0, lwup, &
4110	swdnc0, swdn0, swdn, swupc0, swup0, swup, &
4111	topswad_aero, solswad_aero, &
4112	topswai_aero, solswai_aero, &
4113	topswad0_aero, solswad0_aero, &
4114	topsw_aero, topsw0_aero, &
4115	solsw_aero, solsw0_aero, &
4116	topswcf_aero, solswcf_aero, &
4117	!-C. Kleinschmitt for LW diagnostics
4118	toplwad_aero, sollwad_aero,&
4119	toplwai_aero, sollwai_aero, &
4120	toplwad0_aero, sollwad0_aero,&
4121	!-end
4122	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4123	ZSWFT0_i, ZFSDN0, ZFSUP0)
4124
4125	!lwoff=y, betalwoff=1. : offset LW CRE for radiation code and other
4126	!schemes
4127	toplw = toplw + betalwoff * (toplw0 - toplw)
4128	sollw = sollw + betalwoff * (sollw0 - sollw)
4129	lwdn = lwdn + betalwoff * (lwdn0 - lwdn)
4130	lwup = lwup + betalwoff * (lwup0 - lwup)
4131	sollwdown(:)= sollwdown(:) + betalwoff (-1.ZFLDN0(:,1) - &
4132	sollwdown(:))
4133	cool = cool + betalwoff * (cool0 - cool)
4134
4135	#ifndef CPP_XIOS
4136	!--OB 30/05/2016 modified 21/10/2016
4137	!--here we return swaero_diag and dryaod_diag to FALSE
4138	!--and histdef will switch it back to TRUE if necessary
4139	!--this is necessary to get the right swaero at first step
4140	!--but only in the case of no XIOS as XIOS is covered elsewhere
4141	IF (debut) swaerofree_diag = .FALSE.
4142	IF (debut) swaero_diag = .FALSE.
4143	IF (debut) dryaod_diag = .FALSE.
4144	!--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
4145	!--as for swaero_diag, see above
4146	IF (debut) ok_4xCO2atm = .FALSE.
4147
4148	!
4149	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
4150	!IM des taux doit etre different du taux actuel
4151	!IM Par defaut on a les taux perturbes egaux aux taux actuels
4152	!
4153	IF (RCO2_per.NE.RCO2_act.OR. &
4154	RCH4_per.NE.RCH4_act.OR. &
4155	RN2O_per.NE.RN2O_act.OR. &
4156	RCFC11_per.NE.RCFC11_act.OR. &
4157	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
4158	#endif
4159	!
4160	IF (ok_4xCO2atm) THEN
4161	!
4162	RCO2 = RCO2_per
4163	RCH4 = RCH4_per
4164	RN2O = RN2O_per
4165	RCFC11 = RCFC11_per
4166	RCFC12 = RCFC12_per
4167	!
4168	IF (prt_level .GE.10) THEN
4169	print *,' ->radlwsw, number 2 '
4170	ENDIF
4171	!
4172	CALL radlwsw &
4173	(dist, rmu0, fract, &
4174	!albedo SB >>>
4175	! paprs, pplay,zxtsol,albsol1, albsol2, &
4176	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
4177	!albedo SB <<<
4178	t_seri,q_seri,wo, &
4179	cldfrarad, cldemirad, cldtaurad, &
4180	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, ok_volcan, &
4181	flag_aerosol, &
4182	flag_aerosol_strat, flag_aer_feedback, &
4183	tau_aero, piz_aero, cg_aero, &
4184	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
4185	! Rajoute par OB pour RRTM
4186	tau_aero_lw_rrtm, &
4187	cldtaupi,new_aod, &
4188	! zqsat, flwcrad, fiwcrad, &
4189	zqsat, flwc, fiwc, &
4190	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
4191	heatp,heat0p,coolp,cool0p,albplap, &
4192	heat_volc,cool_volc, &
4193	topswp,toplwp,solswp,sollwp, &
4194	sollwdownp, &
4195	topsw0p,toplw0p,solsw0p,sollw0p, &
4196	lwdnc0p, lwdn0p, lwdnp, lwupc0p, lwup0p, lwupp, &
4197	swdnc0p, swdn0p, swdnp, swupc0p, swup0p, swupp, &
4198	topswad_aerop, solswad_aerop, &
4199	topswai_aerop, solswai_aerop, &
4200	topswad0_aerop, solswad0_aerop, &
4201	topsw_aerop, topsw0_aerop, &
4202	solsw_aerop, solsw0_aerop, &
4203	topswcf_aerop, solswcf_aerop, &
4204	!-C. Kleinschmitt for LW diagnostics
4205	toplwad_aerop, sollwad_aerop,&
4206	toplwai_aerop, sollwai_aerop, &
4207	toplwad0_aerop, sollwad0_aerop,&
4208	!-end
4209	ZLWFT0_i, ZFLDN0, ZFLUP0, &
4210	ZSWFT0_i, ZFSDN0, ZFSUP0)
4211	endif !ok_4xCO2atm
4212	ENDIF ! aerosol_couple
4213	itaprad = 0
4214	!
4215	! If Iflag_radia >=2, reset pertubed variables
4216	!
4217	IF (iflag_radia .ge. 2) THEN
4218	zxtsol(:) = zsav_tsol (:)
4219	ENDIF
4220	ENDIF ! MOD(itaprad,radpas)
4221	itaprad = itaprad + 1
4222
4223	IF (iflag_radia.eq.0) THEN
4224	IF (prt_level.ge.9) THEN
4225	PRINT *,'--------------------------------------------------'
4226	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
4227	PRINT *,'>>>> heat et cool mis a zero '
4228	PRINT *,'--------------------------------------------------'
4229	ENDIF
4230	heat=0.
4231	cool=0.
4232	sollw=0. ! MPL 01032011
4233	solsw=0.
4234	radsol=0.
4235	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
4236	swup0=0.
4237	lwup=0.
4238	lwup0=0.
4239	lwdn=0.
4240	lwdn0=0.
4241	ENDIF
4242
4243	!
4244	! Calculer radsol a l'exterieur de radlwsw
4245	! pour prendre en compte le cycle diurne
4246	! recode par Olivier Boucher en sept 2015
4247	!
4248	radsol=solsw*swradcorr+sollw
4249
4250	IF (ok_4xCO2atm) THEN
4251	radsolp=solswp*swradcorr+sollwp
4252	ENDIF
4253
4254	!
4255	! Ajouter la tendance des rayonnements (tous les pas)
4256	! avec une correction pour le cycle diurne dans le SW
4257	!
4258
4259	DO k=1, klev
4260	d_t_swr(:,k)=swradcorr(:)heat(:,k)phys_tstep/RDAY
4261	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)phys_tstep/RDAY
4262	d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
4263	d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
4264	ENDDO
4265
4266	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
4267	CALL prt_enerbil('SW',itap)
4268	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
4269	CALL prt_enerbil('LW',itap)
4270
4271	!
4272	IF (mydebug) THEN
4273	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4274	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4275	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4276	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4277	ENDIF
4278
4279	! Calculer l'hydrologie de la surface
4280	!
4281	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
4282	! . agesno, ftsol,fqsurf,fsnow, ruis)
4283	!
4284
4285	!
4286	! Calculer le bilan du sol et la derive de temperature (couplage)
4287	!
4288	DO i = 1, klon
4289	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
4290	! a la demande de JLD
4291	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
4292	ENDDO
4293	!
4294	!moddeblott(jan95)
4295	! Appeler le programme de parametrisation de l'orographie
4296	! a l'echelle sous-maille:
4297	!
4298	IF (prt_level .GE.10) THEN
4299	print *,' call orography ? ', ok_orodr
4300	ENDIF
4301	!
4302	IF (ok_orodr) THEN
4303	!
4304	! selection des points pour lesquels le shema est actif:
4305	igwd=0
4306	DO i=1,klon
4307	itest(i)=0
4308	! IF ((zstd(i).gt.10.0)) THEN
4309	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4310	itest(i)=1
4311	igwd=igwd+1
4312	idx(igwd)=i
4313	ENDIF
4314	ENDDO
4315	! igwdim=MAX(1,igwd)
4316	!
4317	IF (ok_strato) THEN
4318
4319	CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
4320	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4321	igwd,idx,itest, &
4322	t_seri, u_seri, v_seri, &
4323	zulow, zvlow, zustrdr, zvstrdr, &
4324	d_t_oro, d_u_oro, d_v_oro)
4325
4326	ELSE
4327	CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
4328	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4329	igwd,idx,itest, &
4330	t_seri, u_seri, v_seri, &
4331	zulow, zvlow, zustrdr, zvstrdr, &
4332	d_t_oro, d_u_oro, d_v_oro)
4333	ENDIF
4334	!
4335	! ajout des tendances
4336	!-----------------------------------------------------------------------
4337	! ajout des tendances de la trainee de l'orographie
4338	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
4339	abortphy,flag_inhib_tend,itap,0)
4340	CALL prt_enerbil('oro',itap)
4341	!----------------------------------------------------------------------
4342	!
4343	ENDIF ! fin de test sur ok_orodr
4344	!
4345	IF (mydebug) THEN
4346	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4347	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4348	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4349	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4350	ENDIF
4351
4352	IF (ok_orolf) THEN
4353	!
4354	! selection des points pour lesquels le shema est actif:
4355	igwd=0
4356	DO i=1,klon
4357	itest(i)=0
4358	IF ((zpic(i)-zmea(i)).GT.100.) THEN
4359	itest(i)=1
4360	igwd=igwd+1
4361	idx(igwd)=i
4362	ENDIF
4363	ENDDO
4364	! igwdim=MAX(1,igwd)
4365	!
4366	IF (ok_strato) THEN
4367
4368	CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
4369	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
4370	igwd,idx,itest, &
4371	t_seri, u_seri, v_seri, &
4372	zulow, zvlow, zustrli, zvstrli, &
4373	d_t_lif, d_u_lif, d_v_lif )
4374
4375	ELSE
4376	CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
4377	latitude_deg,zmea,zstd,zpic, &
4378	itest, &
4379	t_seri, u_seri, v_seri, &
4380	zulow, zvlow, zustrli, zvstrli, &
4381	d_t_lif, d_u_lif, d_v_lif)
4382	ENDIF
4383
4384	! ajout des tendances de la portance de l'orographie
4385	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
4386	'lif', abortphy,flag_inhib_tend,itap,0)
4387	CALL prt_enerbil('lif',itap)
4388	ENDIF ! fin de test sur ok_orolf
4389
4390	IF (ok_hines) then
4391	! HINES GWD PARAMETRIZATION
4392	east_gwstress=0.
4393	west_gwstress=0.
4394	du_gwd_hines=0.
4395	dv_gwd_hines=0.
4396	CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
4397	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4398	du_gwd_hines, dv_gwd_hines)
4399	zustr_gwd_hines=0.
4400	zvstr_gwd_hines=0.
4401	DO k = 1, klev
4402	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
4403	* (paprs(:, k)-paprs(:, k+1))/rg
4404	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
4405	* (paprs(:, k)-paprs(:, k+1))/rg
4406	ENDDO
4407
4408	d_t_hin(:, :)=0.
4409	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4410	dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4411	CALL prt_enerbil('hin',itap)
4412	ENDIF
4413
4414	IF (.not. ok_hines .and. ok_gwd_rando) then
4415	! ym missing init for east_gwstress & west_gwstress -> added in phys_local_var_mod
4416	CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
4417	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4418	dv_gwd_front, east_gwstress, west_gwstress)
4419	zustr_gwd_front=0.
4420	zvstr_gwd_front=0.
4421	DO k = 1, klev
4422	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
4423	* (paprs(:, k)-paprs(:, k+1))/rg
4424	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
4425	* (paprs(:, k)-paprs(:, k+1))/rg
4426	ENDDO
4427
4428	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
4429	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4430	CALL prt_enerbil('front_gwd_rando',itap)
4431	ENDIF
4432
4433	IF (ok_gwd_rando) THEN
4434	CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
4435	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4436	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4437	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
4438	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4439	CALL prt_enerbil('flott_gwd_rando',itap)
4440	zustr_gwd_rando=0.
4441	zvstr_gwd_rando=0.
4442	DO k = 1, klev
4443	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
4444	* (paprs(:, k)-paprs(:, k+1))/rg
4445	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
4446	* (paprs(:, k)-paprs(:, k+1))/rg
4447	ENDDO
4448	ENDIF
4449
4450	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4451
4452	IF (mydebug) THEN
4453	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4454	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4455	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4456	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4457	ENDIF
4458
4459	DO i = 1, klon
4460	zustrph(i)=0.
4461	zvstrph(i)=0.
4462	ENDDO
4463	DO k = 1, klev
4464	DO i = 1, klon
4465	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
4466	(paprs(i,k)-paprs(i,k+1))/rg
4467	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
4468	(paprs(i,k)-paprs(i,k+1))/rg
4469	ENDDO
4470	ENDDO
4471	!
4472	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4473	!
4474	IF (is_sequential .and. ok_orodr) THEN
4475	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4476	ra,rg,romega, &
4477	latitude_deg,longitude_deg,pphis, &
4478	zustrdr,zustrli,zustrph, &
4479	zvstrdr,zvstrli,zvstrph, &
4480	paprs,u,v, &
4481	aam, torsfc)
4482	ENDIF
4483	!IM cf. FLott END
4484	!DC Calcul de la tendance due au methane
4485	IF (ok_qch4) THEN
4486	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4487	! ajout de la tendance d'humidite due au methane
4488	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
4489	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
4490	'q_ch4', abortphy,flag_inhib_tend,itap,0)
4491	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
4492	ENDIF
4493	!
4494	!
4495
4496	!===============================================================
4497	! Additional tendency of TKE due to orography
4498	!===============================================================
4499	!
4500	! Inititialization
4501	!------------------
4502
4503	addtkeoro=0
4504	CALL getin_p('addtkeoro',addtkeoro)
4505
4506	IF (prt_level.ge.5) &
4507	print*,'addtkeoro', addtkeoro
4508
4509	alphatkeoro=1.
4510	CALL getin_p('alphatkeoro',alphatkeoro)
4511	alphatkeoro=min(max(0.,alphatkeoro),1.)
4512
4513	smallscales_tkeoro=.FALSE.
4514	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4515
4516
4517	dtadd(:,:)=0.
4518	duadd(:,:)=0.
4519	dvadd(:,:)=0.
4520
4521	! Choices for addtkeoro:
4522	! ** 0 no TKE tendency from orography
4523	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4524	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4525	!
4526
4527	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4528	! -------------------------------------------
4529
4530
4531	! selection des points pour lesquels le schema est actif:
4532
4533
4534	IF (addtkeoro .EQ. 1 ) THEN
4535
4536	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4537	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4538
4539	ELSE IF (addtkeoro .EQ. 2) THEN
4540
4541	IF (smallscales_tkeoro) THEN
4542	igwd=0
4543	DO i=1,klon
4544	itest(i)=0
4545	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4546	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4547	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4548	IF (zstd(i).GT.1.0) THEN
4549	itest(i)=1
4550	igwd=igwd+1
4551	idx(igwd)=i
4552	ENDIF
4553	ENDDO
4554
4555	ELSE
4556
4557	igwd=0
4558	DO i=1,klon
4559	itest(i)=0
4560	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4561	itest(i)=1
4562	igwd=igwd+1
4563	idx(igwd)=i
4564	ENDIF
4565	ENDDO
4566
4567	ENDIF
4568
4569	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
4570	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4571	igwd,idx,itest, &
4572	t_seri, u_seri, v_seri, &
4573	zulow, zvlow, zustrdr, zvstrdr, &
4574	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4575
4576	zustrdr(:)=0.
4577	zvstrdr(:)=0.
4578	zulow(:)=0.
4579	zvlow(:)=0.
4580
4581	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4582	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4583	ENDIF
4584
4585
4586	! TKE update from subgrid temperature and wind tendencies
4587	!----------------------------------------------------------
4588	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4589
4590
4591	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pctsrf,pbl_tke)
4592
4593
4594	ENDIF
4595	! -----
4596	!===============================================================
4597
4598
4599	!====================================================================
4600	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4601	!====================================================================
4602	! Abderrahmane 24.08.09
4603
4604	IF (ok_cosp) THEN
4605	! adeclarer
4606	#ifdef CPP_COSP
4607	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4608
4609	IF (prt_level .GE.10) THEN
4610	print*,'freq_cosp',freq_cosp
4611	ENDIF
4612	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4613	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4614	! s ref_liq,ref_ice
4615	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
4616	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4617	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4618	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4619	JrNt,ref_liq,ref_ice, &
4620	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4621	zu10m,zv10m,pphis, &
4622	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4623	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4624	prfl(:,1:klev),psfl(:,1:klev), &
4625	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4626	mr_ozone,cldtau, cldemi)
4627
4628	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4629	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4630	! M clMISR,
4631	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4632	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4633
4634	ENDIF
4635	#endif
4636
4637	#ifdef CPP_COSP2
4638	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4639
4640	IF (prt_level .GE.10) THEN
4641	print*,'freq_cosp',freq_cosp
4642	ENDIF
4643	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4644	print*,'Dans physiq.F avant appel '
4645	! s ref_liq,ref_ice
4646	CALL phys_cosp2(itap,phys_tstep,freq_cosp, &
4647	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4648	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4649	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4650	JrNt,ref_liq,ref_ice, &
4651	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4652	zu10m,zv10m,pphis, &
4653	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4654	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4655	prfl(:,1:klev),psfl(:,1:klev), &
4656	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4657	mr_ozone,cldtau, cldemi)
4658	ENDIF
4659	#endif
4660
4661	#ifdef CPP_COSPV2
4662	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4663
4664	IF (prt_level .GE.10) THEN
4665	print*,'freq_cosp',freq_cosp
4666	ENDIF
4667	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4668	print*,'Dans physiq.F avant appel '
4669	! s ref_liq,ref_ice
4670	CALL lmdz_cosp_interface(itap,phys_tstep,freq_cosp, &
4671	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4672	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4673	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4674	JrNt,ref_liq,ref_ice, &
4675	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4676	zu10m,zv10m,pphis, &
4677	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4678	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4679	prfl(:,1:klev),psfl(:,1:klev), &
4680	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4681	mr_ozone,cldtau, cldemi)
4682	ENDIF
4683	#endif
4684
4685	ENDIF !ok_cosp
4686
4687
4688	! Marine
4689
4690	IF (ok_airs) then
4691
4692	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
4693	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4694	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4695	& map_prop_hc,map_prop_hist,&
4696	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4697	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4698	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4699	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4700	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4701	& map_ntot,map_hc,map_hist,&
4702	& map_Cb,map_ThCi,map_Anv,&
4703	& alt_tropo )
4704	ENDIF
4705
4706	ENDIF ! ok_airs
4707
4708
4709	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4710	!AA
4711	!AA Installation de l'interface online-offline pour traceurs
4712	!AA
4713	!====================================================================
4714	! Calcul des tendances traceurs
4715	!====================================================================
4716	!
4717
4718	IF (type_trac=='repr') THEN
4719	sh_in(:,:) = q_seri(:,:)
4720	ELSE
4721	sh_in(:,:) = qx(:,:,ivap)
4722	ch_in(:,:) = qx(:,:,iliq)
4723	ENDIF
4724
4725	IF (iflag_phytrac == 1 ) THEN
4726
4727	#ifdef CPP_Dust
4728	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
4729	pdtphys,ftsol, & ! I
4730	t,q_seri,paprs,pplay,RHcl, & ! I
4731	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
4732	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
4733	u_seri, v_seri, latitude_deg, longitude_deg, &
4734	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
4735	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
4736	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
4737	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
4738	fm_therm, entr_therm, rneb, & ! I
4739	beta_prec_fisrt,beta_prec, & !I
4740	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
4741	d_tr_dyn,tr_seri)
4742
4743	#else
4744
4745	CALL phytrac ( &
4746	itap, days_elapsed+1, jH_cur, debut, &
4747	lafin, phys_tstep, u, v, t, &
4748	paprs, pplay, pmfu, pmfd, &
4749	pen_u, pde_u, pen_d, pde_d, &
4750	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
4751	u1, v1, ftsol, pctsrf, &
4752	zustar, zu10m, zv10m, &
4753	wstar(:,is_ave), ale_bl, ale_wake, &
4754	latitude_deg, longitude_deg, &
4755	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4756	presnivs, pphis, pphi, albsol1, &
4757	sh_in, ch_in, rhcl, cldfra, rneb, &
4758	diafra, cldliq, itop_con, ibas_con, &
4759	pmflxr, pmflxs, prfl, psfl, &
4760	da, phi, mp, upwd, &
4761	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
4762	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
4763	ev, ep, epmlmMm, eplaMm, & !<<RomP
4764	dnwd, aerosol_couple, flxmass_w, &
4765	tau_aero, piz_aero, cg_aero, ccm, &
4766	rfname, &
4767	d_tr_dyn, & !<<RomP
4768	tr_seri, init_source)
4769	#endif
4770	ENDIF ! (iflag_phytrac=1)
4771
4772	IF (offline) THEN
4773
4774	IF (prt_level.ge.9) &
4775	print*,'Attention on met a 0 les thermiques pour phystoke'
4776	CALL phystokenc ( &
4777	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4778	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4779	fm_therm,entr_therm, &
4780	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4781	frac_impa, frac_nucl, &
4782	pphis,cell_area,phys_tstep,itap, &
4783	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4784
4785
4786	ENDIF
4787
4788	!
4789	! Calculer le transport de l'eau et de l'energie (diagnostique)
4790	!
4791	CALL transp (paprs,zxtsol, &
4792	t_seri, q_seri, ql_seri, qs_seri, u_seri, v_seri, zphi, &
4793	ve, vq, ue, uq, vwat, uwat)
4794	!
4795	!IM global posePB BEG
4796	IF(1.EQ.0) THEN
4797	!
4798	CALL transp_lay (paprs,zxtsol, &
4799	t_seri, q_seri, u_seri, v_seri, zphi, &
4800	ve_lay, vq_lay, ue_lay, uq_lay)
4801	!
4802	ENDIF !(1.EQ.0) THEN
4803	!IM global posePB END
4804	! Accumuler les variables a stocker dans les fichiers histoire:
4805	!
4806
4807	!================================================================
4808	! Conversion of kinetic and potential energy into heat, for
4809	! parameterisation of subgrid-scale motions
4810	!================================================================
4811
4812	d_t_ec(:,:)=0.
4813	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4814	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
4815	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4816	zmasse,exner,d_t_ec)
4817	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4818
4819	!=======================================================================
4820	! SORTIES
4821	!=======================================================================
4822	!
4823	!IM initialisation + calculs divers diag AMIP2
4824	!
4825	include "calcul_divers.h"
4826	!
4827	!IM Interpolation sur les niveaux de pression du NMC
4828	! -------------------------------------------------
4829	!
4830	include "calcul_STDlev.h"
4831	!
4832	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
4833	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
4834	!
4835	!cc prw = eau precipitable
4836	! prlw = colonne eau liquide
4837	! prlw = colonne eau solide
4838	prw(:) = 0.
4839	prlw(:) = 0.
4840	prsw(:) = 0.
4841	DO k = 1, klev
4842	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
4843	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
4844	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
4845	ENDDO
4846	!
4847	IF (type_trac == 'inca') THEN
4848	#ifdef INCA
4849	CALL VTe(VTphysiq)
4850	CALL VTb(VTinca)
4851
4852	CALL chemhook_end ( &
4853	phys_tstep, &
4854	pplay, &
4855	t_seri, &
4856	tr_seri, &
4857	nbtr, &
4858	paprs, &
4859	q_seri, &
4860	cell_area, &
4861	pphi, &
4862	pphis, &
4863	zx_rh, &
4864	aps, bps)
4865
4866	CALL VTe(VTinca)
4867	CALL VTb(VTphysiq)
4868	#endif
4869	ENDIF
4870
4871
4872	!
4873	! Convertir les incrementations en tendances
4874	!
4875	IF (prt_level .GE.10) THEN
4876	print *,'Convertir les incrementations en tendances '
4877	ENDIF
4878	!
4879	IF (mydebug) THEN
4880	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4881	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4882	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4883	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4884	ENDIF
4885
4886	DO k = 1, klev
4887	DO i = 1, klon
4888	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
4889	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
4890	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
4891	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
4892	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
4893	!CR: on ajoute le contenu en glace
4894	IF (nqo.eq.3) THEN
4895	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
4896	ENDIF
4897	ENDDO
4898	ENDDO
4899	!
4900	!CR: nb de traceurs eau: nqo
4901	! IF (nqtot.GE.3) THEN
4902	IF (nqtot.GE.(nqo+1)) THEN
4903	! DO iq = 3, nqtot
4904	DO iq = nqo+1, nqtot
4905	DO k = 1, klev
4906	DO i = 1, klon
4907	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / phys_tstep
4908	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / phys_tstep
4909	ENDDO
4910	ENDDO
4911	ENDDO
4912	ENDIF
4913	!
4914	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4915	!IM global posePB include "write_bilKP_ins.h"
4916	!IM global posePB include "write_bilKP_ave.h"
4917	!
4918
4919	!--OB mass fixer
4920	!--profile is corrected to force mass conservation of water
4921	IF (mass_fixer) THEN
4922	qql2(:)=0.0
4923	DO k = 1, klev
4924	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
4925	ENDDO
4926	DO i = 1, klon
4927	!--compute ratio of what q+ql should be with conservation to what it is
4928	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
4929	DO k = 1, klev
4930	q_seri(i,k) =q_seri(i,k)*corrqql
4931	ql_seri(i,k)=ql_seri(i,k)*corrqql
4932	ENDDO
4933	ENDDO
4934	ENDIF
4935	!--fin mass fixer
4936
4937	! Sauvegarder les valeurs de t et q a la fin de la physique:
4938	!
4939	u_ancien(:,:) = u_seri(:,:)
4940	v_ancien(:,:) = v_seri(:,:)
4941	t_ancien(:,:) = t_seri(:,:)
4942	q_ancien(:,:) = q_seri(:,:)
4943	ql_ancien(:,:) = ql_seri(:,:)
4944	qs_ancien(:,:) = qs_seri(:,:)
4945	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
4946	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
4947	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
4948	! !! RomP >>>
4949	!CR: nb de traceurs eau: nqo
4950	IF (nqtot.GT.nqo) THEN
4951	DO iq = nqo+1, nqtot
4952	tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
4953	ENDDO
4954	ENDIF
4955	! !! RomP <<<
4956	!==========================================================================
4957	! Sorties des tendances pour un point particulier
4958	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4959	! pour le debug
4960	! La valeur de igout est attribuee plus haut dans le programme
4961	!==========================================================================
4962
4963	IF (prt_level.ge.1) THEN
4964	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4965	write(lunout,*) &
4966	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4967	write(lunout,*) &
4968	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4969	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4970	pctsrf(igout,is_sic)
4971	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4972	DO k=1,klev
4973	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4974	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4975	d_t_eva(igout,k)
4976	ENDDO
4977	write(lunout,*) 'cool,heat'
4978	DO k=1,klev
4979	write(lunout,*) cool(igout,k),heat(igout,k)
4980	ENDDO
4981
4982	!jyg< (En attendant de statuer sur le sort de d_t_oli)
4983	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4984	!jyg! do k=1,klev
4985	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4986	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4987	!jyg! enddo
4988	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4989	DO k=1,klev
4990	write(lunout,*) d_t_vdf(igout,k), &
4991	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4992	ENDDO
4993	!>jyg
4994
4995	write(lunout,*) 'd_ps ',d_ps(igout)
4996	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4997	DO k=1,klev
4998	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4999	d_qx(igout,k,1),d_qx(igout,k,2)
5000	ENDDO
5001	ENDIF
5002
5003	!============================================================
5004	! Calcul de la temperature potentielle
5005	!============================================================
5006	DO k = 1, klev
5007	DO i = 1, klon
5008	!JYG/IM theta en debut du pas de temps
5009	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
5010	!JYG/IM theta en fin de pas de temps de physique
5011	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
5012	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
5013	! MPL 20130625
5014	! fth_fonctions.F90 et parkind1.F90
5015	! sinon thetal=theta
5016	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
5017	! : ql_seri(i,k))
5018	thetal(i,k)=theta(i,k)
5019	ENDDO
5020	ENDDO
5021	!
5022
5023	! 22.03.04 BEG
5024	!=============================================================
5025	! Ecriture des sorties
5026	!=============================================================
5027	#ifdef CPP_IOIPSL
5028
5029	! Recupere des varibles calcule dans differents modules
5030	! pour ecriture dans histxxx.nc
5031
5032	! Get some variables from module fonte_neige_mod
5033	CALL fonte_neige_get_vars(pctsrf, &
5034	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
5035
5036
5037	!=============================================================
5038	! Separation entre thermiques et non thermiques dans les sorties
5039	! de fisrtilp
5040	!=============================================================
5041
5042	IF (iflag_thermals>=1) THEN
5043	d_t_lscth=0.
5044	d_t_lscst=0.
5045	d_q_lscth=0.
5046	d_q_lscst=0.
5047	DO k=1,klev
5048	DO i=1,klon
5049	IF (ptconvth(i,k)) THEN
5050	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5051	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5052	ELSE
5053	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
5054	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
5055	ENDIF
5056	ENDDO
5057	ENDDO
5058
5059	DO i=1,klon
5060	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
5061	plul_th(i)=prfl(i,1)+psfl(i,1)
5062	ENDDO
5063	ENDIF
5064
5065	!On effectue les sorties:
5066
5067	#ifdef CPP_Dust
5068	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
5069	pplay, lmax_th, aerosol_couple, &
5070	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
5071	ptconv, read_climoz, clevSTD, &
5072	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
5073	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5074	#else
5075	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
5076	pplay, lmax_th, aerosol_couple, &
5077	ok_ade, ok_aie, ok_volcan, ivap, iliq, isol, new_aod, &
5078	ok_sync, ptconv, read_climoz, clevSTD, &
5079	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
5080	flag_aerosol, flag_aerosol_strat, ok_cdnc)
5081	#endif
5082
5083	#ifndef CPP_XIOS
5084	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
5085	#endif
5086
5087	#endif
5088
5089	! On remet des variables a .false. apres un premier appel
5090	IF (debut) THEN
5091	#ifdef CPP_XIOS
5092	swaero_diag=.FALSE.
5093	swaerofree_diag=.FALSE.
5094	dryaod_diag=.FALSE.
5095	ok_4xCO2atm= .FALSE.
5096	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
5097
5098	IF (is_master) THEN
5099	!--setting up swaero_diag to TRUE in XIOS case
5100	IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
5101	xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
5102	xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR. &
5103	(iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
5104	xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0")))) &
5105	!!!--for now these fields are not in the XML files so they are omitted
5106	!!! xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
5107	swaero_diag=.TRUE.
5108
5109	!--setting up swaerofree_diag to TRUE in XIOS case
5110	IF (xios_field_is_active("SWdnSFCcleanclr").OR.xios_field_is_active("SWupSFCcleanclr").OR. &
5111	xios_field_is_active("SWupTOAcleanclr").OR.xios_field_is_active("rsucsaf").OR. &
5112	xios_field_is_active("rsdcsaf") .OR. xios_field_is_active("LWdnSFCcleanclr").OR. &
5113	xios_field_is_active("LWupTOAcleanclr")) &
5114	swaerofree_diag=.TRUE.
5115
5116	!--setting up dryaod_diag to TRUE in XIOS case
5117	DO naero = 1, naero_tot-1
5118	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
5119	ENDDO
5120	!
5121	!--setting up ok_4xCO2atm to TRUE in XIOS case
5122	IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
5123	xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
5124	xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
5125	xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
5126	xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
5127	xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
5128	ok_4xCO2atm=.TRUE.
5129	ENDIF
5130	!$OMP BARRIER
5131	CALL bcast(swaero_diag)
5132	CALL bcast(swaerofree_diag)
5133	CALL bcast(dryaod_diag)
5134	CALL bcast(ok_4xCO2atm)
5135	! write (lunout,*)'ok_4xCO2atm= ',swaero_diag, swaerofree_diag, dryaod_diag, ok_4xCO2atm
5136	#endif
5137	ENDIF
5138
5139	!====================================================================
5140	! Arret du modele apres hgardfou en cas de detection d'un
5141	! plantage par hgardfou
5142	!====================================================================
5143
5144	IF (abortphy==1) THEN
5145	abort_message ='Plantage hgardfou'
5146	CALL abort_physic (modname,abort_message,1)
5147	ENDIF
5148
5149	! 22.03.04 END
5150	!
5151	!====================================================================
5152	! Si c'est la fin, il faut conserver l'etat de redemarrage
5153	!====================================================================
5154	!
5155
5156	IF (lafin) THEN
5157	itau_phy = itau_phy + itap
5158	CALL phyredem ("restartphy.nc")
5159	! open(97,form="unformatted",file="finbin")
5160	! write(97) u_seri,v_seri,t_seri,q_seri
5161	! close(97)
5162
5163	IF (is_omp_master) THEN
5164
5165	IF (read_climoz >= 1) THEN
5166	IF (is_mpi_root) CALL nf95_close(ncid_climoz)
5167	DEALLOCATE(press_edg_climoz) ! pointer
5168	DEALLOCATE(press_cen_climoz) ! pointer
5169	ENDIF
5170
5171	ENDIF
5172	#ifdef CPP_XIOS
5173	IF (is_omp_master) CALL xios_context_finalize
5174	#endif
5175	WRITE(lunout,*) ' physiq fin, nombre de steps ou cvpas = 1 : ', Ncvpaseq1
5176	ENDIF
5177
5178	! first=.false.
5179
5180	END SUBROUTINE physiq
5181
5182	END MODULE physiq_mod

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

Download in other formats:

Original Format