Context Navigation

physiq_mod.F90 @ 3409

Last change on this file since 3409 was 3387, checked in by oboucher, 6 years ago
adding infocfields_init routine and its call from physiq. This routines reads which fields need to be transferred between model components in ESM configuration. No impact whatsoever in LMDZ mode.
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: 172.3 KB

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

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

Download in other formats:

Original Format