Context Navigation

physiq_mod.F90 @ 3513

Last change on this file since 3513 was 3512, checked in by idelkadi, 6 years ago
Implementation de la version cospv2 dans LMDZ (suite)
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: 182.5 KB

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

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

Download in other formats:

Original Format