Context Navigation

physiq_mod.F90 @ 2635

Last change on this file since 2635 was 2635, checked in by jyg, 8 years ago
Some cleaning in the wake routines: (i) the wake number per unit area (wdens) is now a state variable (held constant for the time being); (ii) wake state variable changes are computed in subroutine 'physiq' if iflag_wake_tend=1 (it is computed within wake routines if iflag_wake_tend=0, consistent with earlier versions); (iii) the new routine 'add_wake_tend' adds tendencies to wake state variables; (iv) tendencies due to various processes (pbl, wakes, thermals) are named and added separately.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 160.1 KB

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

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

Download in other formats:

Original Format