Context Navigation

physiq_mod.F90 @ 3150

Last change on this file since 3150 was 3150, checked in by jyg, 6 years ago
Introducing upper bounds on the convective tendencies above which convection scheme is called every time step.
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: 167.6 KB

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

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

Download in other formats:

Original Format