Context Navigation

physiq_mod.F90 @ 3071

Last change on this file since 3071 was 3071, checked in by jyg, 7 years ago

Correction in ener_conserv.F90: bils_tke is
set to zero if iflag_pbl <= 1.

Cosmetic correction in physiq_mod.F90: pbl_tke is
initialized to zero if iflag_pbl <= 1 and
klon_glo=1.

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: 164.4 KB

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

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

Download in other formats:

Original Format