Context Navigation

physiq_mod.F90 @ 3826

Last change on this file since 3826 was 3826, checked in by musat, 3 years ago
Nouveaux calculs 2m/10m activables par iflag_new_t2mq2m=1 Valeur par defaut iflag_new_t2mq2m=0 (<=> anciens calculs CMIP6)
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: 171.9 KB

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

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

Download in other formats:

Original Format