Context Navigation

physiq_mod.F90 @ 2994

Last change on this file since 2994 was 2994, checked in by oboucher, 7 years ago
Addition to previous commit. Call of stratomask to diagnose tropopause level according to WMO definition Previous commit aimed at diagnosing tropospheric and stratospheric O3 columns in DU.
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: 162.0 KB

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

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

Download in other formats:

Original Format