Context Navigation

physiq_mod.F90 @ 3351

Last change on this file since 3351 was 3351, checked in by acozic, 6 years ago
correct a bug on a call for Inca
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.2 KB

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

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

Download in other formats:

Original Format