Context Navigation

physiq_mod.F90 @ 3777

Last change on this file since 3777 was 3777, checked in by lmdz-user, 4 years ago
Moved the opening and reading of the physics restart file before the initialisation of the physics output as the output procedure needs to know the timestep at the beginning of the run to initialise the calendar correctly. That timestep is read in the physics restart file. This should close tickets #109 and #114
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: 187.6 KB

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

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

Download in other formats:

Original Format