Context Navigation

physiq_mod.F90 @ 3355

Last change on this file since 3355 was 3355, checked in by Laurent Fairhead, 6 years ago
Commiting variables initialisation as found by Yann LF
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 `Id`
File size: 165.0 KB

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

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

Download in other formats:

Original Format