Context Navigation

physiq_mod.F90 @ 3412

Last change on this file since 3412 was 3412, checked in by oboucher, 6 years ago
Externalizing a controlling parameter for activating aerosol radiative feedback Need to pass the info all the way to recmwf_aero. Default is True as before.
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.5 KB

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

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

Download in other formats:

Original Format