Context Navigation

physiq_mod.F90 @ 4098

Last change on this file since 4098 was 4098, checked in by oboucher, 2 years ago
tidying up and correcting a few omissions for nqo > 3
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 `Id`
File size: 190.1 KB

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

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

Download in other formats:

Original Format