Context Navigation

physiq_mod.F90 @ 3301

Last change on this file since 3301 was 3301, checked in by Laurent Fairhead, 6 years ago
Disabling vertical interpolation at each physics timestep as it is not needed for CMIP6 workflow
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 169.4 KB

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

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

Download in other formats:

Original Format