Context Navigation

physiq_mod.F90 @ 3110

Last change on this file since 3110 was 3110, checked in by Laurent Fairhead, 6 years ago
Added sza, solar zenithal angle, to outputs for CMIP6 dataresquest
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: 165.0 KB

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

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

Download in other formats:

Original Format