Context Navigation

physiq_mod.F90 @ 3322

Last change on this file since 3322 was 3322, checked in by Laurent Fairhead, 6 years ago
Continuing convergence between DYNAMICO and LMDZ new physics. The xios stuff
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 163.4 KB

Line
1	!
2	! $Id: physiq_mod.F90 3322 2018-05-04 13:33:28Z fairhead $
3	!
4	!#define IO_DEBUG
5	MODULE physiq_mod
6
7	IMPLICIT NONE
8
9	CONTAINS
10
11	SUBROUTINE physiq (nlon,nlev, &
12	debut,lafin,pdtphys_, &
13	paprs,pplay,pphi,pphis,presnivs, &
14	u,v,rot,t,qx, &
15	flxmass_w, &
16	d_u, d_v, d_t, d_qx, d_ps)
17
18	use assert_m, only: assert
19	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
20	histwrite, ju2ymds, ymds2ju, getin
21	USE geometry_mod, ONLY: cell_area, latitude_deg, longitude_deg
22	USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, &
23	year_cur, mth_cur,jD_cur, jH_cur, jD_ref, day_cur, hour
24	USE write_field_phy
25	USE dimphy
26	USE infotrac_phy, ONLY: nqtot, nbtr, nqo, type_trac
27	USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, nbp_lev, klon_glo
28	USE mod_phys_lmdz_para
29	USE iophy
30	USE print_control_mod, ONLY: mydebug=>debug , lunout, prt_level
31	USE phystokenc_mod, ONLY: offline, phystokenc
32	USE time_phylmdz_mod, only: raz_date, day_step_phy, update_time,current_time
33	USE vampir
34	USE pbl_surface_mod, ONLY : pbl_surface
35	USE change_srf_frac_mod
36	USE surface_data, ONLY : type_ocean, ok_veget, ok_snow
37	USE tropopause_m, ONLY: dyn_tropopause
38	#ifdef CPP_Dust
39	USE phytracr_spl_mod, ONLY: phytracr_spl
40	#endif
41	USE phys_local_var_mod, ONLY: phys_local_var_init, phys_local_var_end, &
42	! [Variables internes non sauvegardees de la physique]
43	! Variables locales pour effectuer les appels en serie
44	t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,tr_seri, &
45	! Dynamic tendencies (diagnostics)
46	d_t_dyn,d_q_dyn,d_ql_dyn,d_qs_dyn,d_u_dyn,d_v_dyn,d_tr_dyn, &
47	d_q_dyn2d,d_ql_dyn2d,d_qs_dyn2d, &
48	! Physic tendencies
49	d_t_con,d_q_con,d_u_con,d_v_con, &
50	d_tr, & !! to be removed?? (jyg)
51	d_t_wake,d_q_wake, &
52	d_t_lwr,d_t_lw0,d_t_swr,d_t_sw0, &
53	d_t_ajsb,d_q_ajsb, &
54	d_t_ajs,d_q_ajs,d_u_ajs,d_v_ajs, &
55	d_t_ajs_w,d_q_ajs_w, &
56	d_t_ajs_x,d_q_ajs_x, &
57	!
58	d_t_eva,d_q_eva,d_ql_eva,d_qi_eva, &
59	d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc, &
60	d_t_lscst,d_q_lscst, &
61	d_t_lscth,d_q_lscth, &
62	plul_st,plul_th, &
63	!
64	d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_t_diss, &
65	d_t_vdf_w,d_q_vdf_w, &
66	d_t_vdf_x,d_q_vdf_x, &
67	d_ts, &
68	!
69	d_t_oli,d_u_oli,d_v_oli, &
70	d_t_oro,d_u_oro,d_v_oro, &
71	d_t_oro_gw,d_u_oro_gw,d_v_oro_gw, &
72	d_t_lif,d_u_lif,d_v_lif, &
73	d_t_ec, &
74	!
75	du_gwd_hines,dv_gwd_hines,d_t_hin, &
76	dv_gwd_rando,dv_gwd_front, &
77	east_gwstress,west_gwstress, &
78	d_q_ch4, &
79	! Special RRTM
80	ZLWFT0_i,ZSWFT0_i,ZFLDN0, &
81	ZFLUP0,ZFSDN0,ZFSUP0, &
82	!
83	topswad_aero,solswad_aero, &
84	topswai_aero,solswai_aero, &
85	topswad0_aero,solswad0_aero, &
86	!LW additional
87	toplwad_aero,sollwad_aero, &
88	toplwai_aero,sollwai_aero, &
89	toplwad0_aero,sollwad0_aero, &
90	!
91	topsw_aero,solsw_aero, &
92	topsw0_aero,solsw0_aero, &
93	topswcf_aero,solswcf_aero, &
94	tausum_aero,tau3d_aero, &
95	drytausum_aero, &
96	!
97	!variables CFMIP2/CMIP5
98	topswad_aerop, solswad_aerop, &
99	topswai_aerop, solswai_aerop, &
100	topswad0_aerop, solswad0_aerop, &
101	topsw_aerop, topsw0_aerop, &
102	solsw_aerop, solsw0_aerop, &
103	topswcf_aerop, solswcf_aerop, &
104	!LW diagnostics
105	toplwad_aerop, sollwad_aerop, &
106	toplwai_aerop, sollwai_aerop, &
107	toplwad0_aerop, sollwad0_aerop, &
108	!
109	ptstar, pt0, slp, &
110	!
111	bils, &
112	!
113	cldh, cldl,cldm, cldq, cldt, &
114	JrNt, &
115	dthmin, evap, fder, plcl, plfc, &
116	prw, prlw, prsw, &
117	s_lcl, s_pblh, s_pblt, s_therm, &
118	cdragm, cdragh, &
119	zustar, zu10m, zv10m, rh2m, qsat2m, &
120	zq2m, zt2m, weak_inversion, &
121	zt2m_min_mon, zt2m_max_mon, & ! pour calcul_divers.h
122	t2m_min_mon, t2m_max_mon, & ! pour calcul_divers.h
123	!
124	s_pblh_x, s_pblh_w, &
125	s_lcl_x, s_lcl_w, &
126	!
127	slab_wfbils, tpot, tpote, &
128	ue, uq, ve, vq, zxffonte, &
129	zxfqcalving, zxfluxlat, &
130	zxrunofflic, &
131	zxtsol, snow_lsc, zxfqfonte, zxqsurf, &
132	rain_lsc, rain_num, &
133	!
134	sens_x, sens_w, &
135	zxfluxlat_x, zxfluxlat_w, &
136	!
137	dtvdf_x, dtvdf_w, &
138	dqvdf_x, dqvdf_w, &
139	pbl_tke_input, &
140	t_therm, q_therm, u_therm, v_therm, &
141	cdragh_x, cdragh_w, &
142	cdragm_x, cdragm_w, &
143	kh, kh_x, kh_w, &
144	!
145	wake_k, &
146	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	INTEGER :: nbtr_tmp ! Number of tracer inside concvl
1108	REAL, dimension(klon,klev) :: sh_in ! Specific humidity entering in phytrac
1109	REAL, dimension(klon,klev) :: ch_in ! Condensed humidity entering in phytrac (eau liquide)
1110	integer iostat
1111
1112	REAL zzz
1113	!albedo SB >>>
1114	real,dimension(6),save :: SFRWL
1115	!$OMP THREADPRIVATE(SFRWL)
1116	!albedo SB <<<
1117
1118	!--OB variables for mass fixer (hard coded for now)
1119	logical, parameter :: mass_fixer=.false.
1120	real qql1(klon),qql2(klon),corrqql
1121
1122	! Ehouarn: set value of jjmp1 since it is no longer a "fixed parameter"
1123	jjmp1=nbp_lat
1124
1125	!======================================================================
1126	! Gestion calendrier : mise a jour du module phys_cal_mod
1127	!
1128	pdtphys=pdtphys_
1129	CALL update_time(pdtphys)
1130	phys_tstep=NINT(pdtphys)
1131
1132	!======================================================================
1133	! Ecriture eventuelle d'un profil verticale en entree de la physique.
1134	! Utilise notamment en 1D mais peut etre active egalement en 3D
1135	! en imposant la valeur de igout.
1136	!======================================================================d
1137	IF (prt_level.ge.1) THEN
1138	igout=klon/2+1/klon
1139	write(lunout,*) 'DEBUT DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
1140	write(lunout,*) 'igout, lat, lon ',igout, latitude_deg(igout), &
1141	longitude_deg(igout)
1142	write(lunout,*) &
1143	'nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys'
1144	write(lunout,*) &
1145	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur,pdtphys
1146
1147	write(lunout,*) 'paprs, play, phi, u, v, t'
1148	DO k=1,klev
1149	write(lunout,*) paprs(igout,k),pplay(igout,k),pphi(igout,k), &
1150	u(igout,k),v(igout,k),t(igout,k)
1151	ENDDO
1152	write(lunout,*) 'ovap (g/kg), oliq (g/kg)'
1153	DO k=1,klev
1154	write(lunout,) qx(igout,k,1)1000,qx(igout,k,2)*1000.
1155	ENDDO
1156	ENDIF
1157
1158	! Quick check on pressure levels:
1159	call assert(paprs(:, nbp_lev + 1) < paprs(:, nbp_lev), &
1160	"physiq_mod paprs bad order")
1161
1162	IF (first) THEN
1163	!CR:nvelles variables convection/poches froides
1164
1165	print*, '================================================='
1166	print*, 'Allocation des variables locales et sauvegardees'
1167	CALL phys_local_var_init
1168	!
1169	! appel a la lecture du run.def physique
1170	CALL conf_phys(ok_journe, ok_mensuel, &
1171	ok_instan, ok_hf, &
1172	ok_LES, &
1173	callstats, &
1174	solarlong0,seuil_inversion, &
1175	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
1176	iflag_cld_th,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
1177	ok_ade, ok_aie, ok_alw, ok_cdnc, aerosol_couple, &
1178	flag_aerosol, flag_aerosol_strat, new_aod, &
1179	flag_bc_internal_mixture, bl95_b0, bl95_b1, &
1180	! nv flags pour la convection et les
1181	! poches froides
1182	read_climoz, &
1183	alp_offset)
1184	CALL phys_state_var_init(read_climoz)
1185	CALL phys_output_var_init
1186	print*, '================================================='
1187	!
1188	!CR: check sur le nb de traceurs de l eau
1189	IF ((iflag_ice_thermo.gt.0).and.(nqo==2)) THEN
1190	WRITE (lunout, *) ' iflag_ice_thermo==1 requires 3 H2O tracers ', &
1191	'(H2Ov, H2Ol, H2Oi) but nqo=', nqo, '. Might as well stop here.'
1192	STOP
1193	ENDIF
1194
1195	dnwd0=0.0
1196	ftd=0.0
1197	fqd=0.0
1198	cin=0.
1199	!ym Attention pbase pas initialise dans concvl !!!!
1200	pbase=0
1201	!IM 180608
1202
1203
1204	itau_con=0
1205	first=.false.
1206
1207	ENDIF ! first
1208
1209	!ym => necessaire pour iflag_con != 2
1210	pmfd(:,:) = 0.
1211	pen_u(:,:) = 0.
1212	pen_d(:,:) = 0.
1213	pde_d(:,:) = 0.
1214	pde_u(:,:) = 0.
1215	aam=0.
1216	d_t_adjwk(:,:)=0
1217	d_q_adjwk(:,:)=0
1218
1219	alp_bl_conv(:)=0.
1220
1221	torsfc=0.
1222	forall (k=1: nbp_lev) zmasse(:, k) = (paprs(:, k)-paprs(:, k+1)) / rg
1223
1224	modname = 'physiq'
1225
1226	IF (debut) THEN
1227	CALL suphel ! initialiser constantes et parametres phys.
1228	CALL getin_p('iflag_alp_wk_cond', iflag_alp_wk_cond)
1229	CALL getin_p('random_notrig_max',random_notrig_max)
1230	CALL getin_p('ok_adjwk',ok_adjwk)
1231	IF (ok_adjwk) iflag_adjwk=2 ! for compatibility with older versions
1232	! iflag_adjwk: ! 0 = Default: no convective adjustment of w-region
1233	! 1 => convective adjustment but state variables are unchanged
1234	! 2 => convective adjustment and state variables are changed
1235	CALL getin_p('iflag_adjwk',iflag_adjwk)
1236	CALL getin_p('oliqmax',oliqmax)
1237	CALL getin_p('oicemax',oicemax)
1238	CALL getin_p('ratqsp0',ratqsp0)
1239	CALL getin_p('ratqsdp',ratqsdp)
1240	iflag_wake_tend = 0
1241	CALL getin_p('iflag_wake_tend',iflag_wake_tend)
1242	ok_bad_ecmwf_thermo=.TRUE. ! By default thermodynamical constants are set
1243	! in rrtm/suphec.F90 (and rvtmp2 is set to 0).
1244	CALL getin_p('ok_bad_ecmwf_thermo',ok_bad_ecmwf_thermo)
1245	fl_ebil = 0 ! by default, conservation diagnostics are desactivated
1246	CALL getin_p('fl_ebil',fl_ebil)
1247	fl_cor_ebil = 0 ! by default, no correction to ensure energy conservation
1248	CALL getin_p('fl_cor_ebil',fl_cor_ebil)
1249	iflag_phytrac = 1 ! by default we do want to call phytrac
1250	CALL getin_p('iflag_phytrac',iflag_phytrac)
1251	nvm_lmdz = 13
1252	CALL getin_p('NVM',nvm_lmdz)
1253	ENDIF
1254
1255	IF (prt_level.ge.1) print *,'CONVERGENCE PHYSIQUE THERM 1 '
1256
1257
1258	!======================================================================
1259	! Gestion calendrier : mise a jour du module phys_cal_mod
1260	!
1261	! CALL phys_cal_update(jD_cur,jH_cur)
1262
1263	!
1264	! Si c'est le debut, il faut initialiser plusieurs choses
1265	! ********
1266	!
1267	IF (debut) THEN
1268	!rv CRinitialisation de wght_th et lalim_conv pour la
1269	!definition de la couche alimentation de la convection a partir
1270	!des caracteristiques du thermique
1271	wght_th(:,:)=1.
1272	lalim_conv(:)=1
1273	!RC
1274	ustar(:,:)=0.
1275	! u10m(:,:)=0.
1276	! v10m(:,:)=0.
1277	rain_con(:)=0.
1278	snow_con(:)=0.
1279	topswai(:)=0.
1280	topswad(:)=0.
1281	solswai(:)=0.
1282	solswad(:)=0.
1283
1284	wmax_th(:)=0.
1285	tau_overturning_th(:)=0.
1286
1287	IF (type_trac == 'inca') THEN
1288	! jg : initialisation jusqu'au ces variables sont dans restart
1289	ccm(:,:,:) = 0.
1290	tau_aero(:,:,:,:) = 0.
1291	piz_aero(:,:,:,:) = 0.
1292	cg_aero(:,:,:,:) = 0.
1293
1294	config_inca='none' ! default
1295	CALL getin_p('config_inca',config_inca)
1296
1297	ELSE
1298	config_inca='none' ! default
1299	ENDIF
1300
1301	IF (aerosol_couple .AND. (config_inca /= "aero" &
1302	.AND. config_inca /= "aeNP ")) THEN
1303	abort_message &
1304	= 'if aerosol_couple is activated, config_inca need to be ' &
1305	// 'aero or aeNP'
1306	CALL abort_physic (modname,abort_message,1)
1307	ENDIF
1308
1309
1310
1311	rnebcon0(:,:) = 0.0
1312	clwcon0(:,:) = 0.0
1313	rnebcon(:,:) = 0.0
1314	clwcon(:,:) = 0.0
1315
1316	!
1317	print*,'iflag_coupl,iflag_clos,iflag_wake', &
1318	iflag_coupl,iflag_clos,iflag_wake
1319	print*,'iflag_CYCLE_DIURNE', iflag_cycle_diurne
1320	!
1321	IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
1322	abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
1323	CALL abort_physic (modname,abort_message,1)
1324	ENDIF
1325	!
1326	!
1327	! Initialiser les compteurs:
1328	!
1329	itap = 0
1330	itaprad = 0
1331	itapcv = 0
1332	itapwk = 0
1333
1334	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1335	!! Un petit travail \`a faire ici.
1336	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1337
1338	IF (iflag_pbl>1) THEN
1339	PRINT*, "Using method MELLOR&YAMADA"
1340	ENDIF
1341
1342	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1343	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans
1344	! phylmd plutot que dyn3d
1345	! Attention : la version precedente n'etait pas tres propre.
1346	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1347	! pour obtenir le meme resultat.
1348	!jyg for fh<
1349	WRITE(lunout,*) 'Pas de temps phys_tstep pdtphys ',phys_tstep,pdtphys
1350	IF (abs(phys_tstep-pdtphys)>1.e-10) THEN
1351	abort_message='pas de temps doit etre entier en seconde pour orchidee et XIOS'
1352	CALL abort_physic(modname,abort_message,1)
1353	ENDIF
1354	!>jyg
1355	IF (MOD(NINT(86400./phys_tstep),nbapp_rad).EQ.0) THEN
1356	radpas = NINT( 86400./phys_tstep)/nbapp_rad
1357	ELSE
1358	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1359	'multiple de nbapp_rad'
1360	WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test ', &
1361	'mais 1+1<>2'
1362	abort_message='nbre de pas de temps physique n est pas multiple ' &
1363	// 'de nbapp_rad'
1364	CALL abort_physic(modname,abort_message,1)
1365	ENDIF
1366	IF (nbapp_cv .EQ. 0) nbapp_cv=86400./phys_tstep
1367	IF (nbapp_wk .EQ. 0) nbapp_wk=86400./phys_tstep
1368	print *,'physiq, nbapp_cv, nbapp_wk ',nbapp_cv,nbapp_wk
1369	IF (MOD(NINT(86400./phys_tstep),nbapp_cv).EQ.0) THEN
1370	cvpas = NINT( 86400./phys_tstep)/nbapp_cv
1371	print *,'physiq, cvpas ',cvpas
1372	ELSE
1373	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1374	'multiple de nbapp_cv'
1375	WRITE(lunout,*) 'changer nbapp_cv ou alors commenter ce test ', &
1376	'mais 1+1<>2'
1377	abort_message='nbre de pas de temps physique n est pas multiple ' &
1378	// 'de nbapp_cv'
1379	call abort_physic(modname,abort_message,1)
1380	ENDIF
1381	IF (MOD(NINT(86400./phys_tstep),nbapp_wk).EQ.0) THEN
1382	wkpas = NINT( 86400./phys_tstep)/nbapp_wk
1383	print *,'physiq, wkpas ',wkpas
1384	ELSE
1385	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un ', &
1386	'multiple de nbapp_wk'
1387	WRITE(lunout,*) 'changer nbapp_wk ou alors commenter ce test ', &
1388	'mais 1+1<>2'
1389	abort_message='nbre de pas de temps physique n est pas multiple ' &
1390	// 'de nbapp_wk'
1391	call abort_physic(modname,abort_message,1)
1392	ENDIF
1393	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1394	CALL init_iophy_new(latitude_deg,longitude_deg)
1395
1396
1397	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1398	!jyg<
1399	IF (klon_glo==1) THEN
1400	pbl_tke(:,:,is_ave) = 0.
1401	DO nsrf=1,nbsrf
1402	DO k = 1,klev+1
1403	pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) &
1404	+pctsrf(:,nsrf)*pbl_tke(:,k,nsrf)
1405	ENDDO
1406	ENDDO
1407	!>jyg
1408	ENDIF
1409	!IM begin
1410	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1411	,ratqs(1,1)
1412	!IM end
1413
1414
1415	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1416	!
1417	! on remet le calendrier a zero
1418	!
1419	IF (raz_date .eq. 1) THEN
1420	itau_phy = 0
1421	ENDIF
1422
1423	CALL printflag( tabcntr0,radpas,ok_journe, &
1424	ok_instan, ok_region )
1425	!
1426	! IF (ABS(phys_tstep-pdtphys).GT.0.001) THEN
1427	! WRITE(lunout,*) 'Pas physique n est pas correct',phys_tstep, &
1428	! pdtphys
1429	! abort_message='Pas physique n est pas correct '
1430	! ! call abort_physic(modname,abort_message,1)
1431	! phys_tstep=pdtphys
1432	! ENDIF
1433	IF (nlon .NE. klon) THEN
1434	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1435	klon
1436	abort_message='nlon et klon ne sont pas coherents'
1437	CALL abort_physic(modname,abort_message,1)
1438	ENDIF
1439	IF (nlev .NE. klev) THEN
1440	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1441	klev
1442	abort_message='nlev et klev ne sont pas coherents'
1443	CALL abort_physic(modname,abort_message,1)
1444	ENDIF
1445	!
1446	IF (phys_tstep*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
1447	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1448	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1449	abort_message='Nbre d appels au rayonnement insuffisant'
1450	CALL abort_physic(modname,abort_message,1)
1451	ENDIF
1452	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1453	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1454	ok_cvl
1455	!
1456	!KE43
1457	! Initialisation pour la convection de K.E. (sb):
1458	IF (iflag_con.GE.3) THEN
1459
1460	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1461	WRITE(lunout,*) &
1462	"On va utiliser le melange convectif des traceurs qui"
1463	WRITE(lunout,*)"est calcule dans convect4.3"
1464	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1465
1466	DO i = 1, klon
1467	ema_cbmf(i) = 0.
1468	ema_pcb(i) = 0.
1469	ema_pct(i) = 0.
1470	! ema_workcbmf(i) = 0.
1471	ENDDO
1472	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1473	DO i = 1, klon
1474	ibas_con(i) = 1
1475	itop_con(i) = 1
1476	ENDDO
1477	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1478	!================================================================
1479	!CR:04.12.07: initialisations poches froides
1480	! Controle de ALE et ALP pour la fermeture convective (jyg)
1481	IF (iflag_wake>=1) THEN
1482	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1483	,alp_bl_prescr, ale_bl_prescr)
1484	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1485	! print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
1486	!
1487	! Initialize tendencies of wake state variables (for some flag values
1488	! they are not computed).
1489	d_deltat_wk(:,:) = 0.
1490	d_deltaq_wk(:,:) = 0.
1491	d_deltat_wk_gw(:,:) = 0.
1492	d_deltaq_wk_gw(:,:) = 0.
1493	d_deltat_vdf(:,:) = 0.
1494	d_deltaq_vdf(:,:) = 0.
1495	d_deltat_the(:,:) = 0.
1496	d_deltaq_the(:,:) = 0.
1497	d_deltat_ajs_cv(:,:) = 0.
1498	d_deltaq_ajs_cv(:,:) = 0.
1499	d_s_wk(:) = 0.
1500	d_dens_wk(:) = 0.
1501	ENDIF
1502
1503	! do i = 1,klon
1504	! Ale_bl(i)=0.
1505	! Alp_bl(i)=0.
1506	! enddo
1507
1508	!===================================================================
1509	!IM stations CFMIP
1510	nCFMIP=npCFMIP
1511	OPEN(98,file='npCFMIP_param.data',status='old', &
1512	form='formatted',iostat=iostat)
1513	IF (iostat == 0) THEN
1514	READ(98,*,end=998) nCFMIP
1515	998 CONTINUE
1516	CLOSE(98)
1517	CONTINUE
1518	IF(nCFMIP.GT.npCFMIP) THEN
1519	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1520	CALL abort_physic("physiq", "", 1)
1521	ELSE
1522	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1523	ENDIF
1524
1525	!
1526	ALLOCATE(tabCFMIP(nCFMIP))
1527	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1528	ALLOCATE(tabijGCM(nCFMIP))
1529	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1530	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1531	!
1532	! lecture des nCFMIP stations CFMIP, de leur numero
1533	! et des coordonnees geographiques lonCFMIP, latCFMIP
1534	!
1535	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
1536	lonCFMIP, latCFMIP)
1537	!
1538	! identification des
1539	! 1) coordonnees lonGCM, latGCM des points CFMIP dans la
1540	! grille de LMDZ
1541	! 2) indices points tabijGCM de la grille physique 1d sur
1542	! klon points
1543	! 3) indices iGCM, jGCM de la grille physique 2d
1544	!
1545	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1546	tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1547	!
1548	ELSE
1549	ALLOCATE(tabijGCM(0))
1550	ALLOCATE(lonGCM(0), latGCM(0))
1551	ALLOCATE(iGCM(0), jGCM(0))
1552	ENDIF
1553	ELSE
1554	ALLOCATE(tabijGCM(0))
1555	ALLOCATE(lonGCM(0), latGCM(0))
1556	ALLOCATE(iGCM(0), jGCM(0))
1557	ENDIF
1558
1559	DO i=1,klon
1560	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1561	ENDDO
1562
1563	!34EK
1564	IF (ok_orodr) THEN
1565
1566	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1567	! FH sans doute a enlever de finitivement ou, si on le
1568	! garde, l'activer justement quand ok_orodr = false.
1569	! ce rugoro est utilise par la couche limite et fait double emploi
1570	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1571	! DO i=1,klon
1572	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1573	! ENDDO
1574	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1575	IF (ok_strato) THEN
1576	CALL SUGWD_strato(klon,klev,paprs,pplay)
1577	ELSE
1578	CALL SUGWD(klon,klev,paprs,pplay)
1579	ENDIF
1580
1581	DO i=1,klon
1582	zuthe(i)=0.
1583	zvthe(i)=0.
1584	IF (zstd(i).gt.10.) THEN
1585	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1586	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1587	ENDIF
1588	ENDDO
1589	ENDIF
1590	!
1591	!
1592	lmt_pas = NINT(86400./phys_tstep * 1.0) ! tous les jours
1593	WRITE(lunout,*)'La frequence de lecture surface est de ', &
1594	lmt_pas
1595	!
1596	capemaxcels = 't_max(X)'
1597	t2mincels = 't_min(X)'
1598	t2maxcels = 't_max(X)'
1599	tinst = 'inst(X)'
1600	tave = 'ave(X)'
1601	!IM cf. AM 081204 BEG
1602	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1603	!IM cf. AM 081204 END
1604	!
1605	!=============================================================
1606	! Initialisation des sorties
1607	!=============================================================
1608
1609	#ifdef CPP_XIOS
1610	! Get "missing_val" value from XML files (from temperature variable)
1611	!$OMP MASTER
1612	CALL xios_get_field_attr("temp",default_value=missing_val_omp)
1613	!$OMP END MASTER
1614	!$OMP BARRIER
1615	missing_val=missing_val_omp
1616	#endif
1617
1618	#ifdef CPP_XIOS
1619	!--setting up swaero_diag to TRUE in XIOS case
1620	IF (xios_field_is_active("topswad").OR.xios_field_is_active("topswad0").OR. &
1621	xios_field_is_active("solswad").OR.xios_field_is_active("solswad0").OR. &
1622	xios_field_is_active("topswai").OR.xios_field_is_active("solswai").OR. &
1623	(iflag_rrtm==1.AND.(xios_field_is_active("toplwad").OR.xios_field_is_active("toplwad0").OR. &
1624	xios_field_is_active("sollwad").OR.xios_field_is_active("sollwad0")))) &
1625	!!!--for now these fields are not in the XML files so they are omitted
1626	!!! xios_field_is_active("toplwai").OR.xios_field_is_active("sollwai") !))) &
1627	swaero_diag=.TRUE.
1628
1629	!--setting up dryaod_diag to TRUE in XIOS case
1630	DO naero = 1, naero_tot-1
1631	IF (xios_field_is_active("dryod550_"//name_aero_tau(naero))) dryaod_diag=.TRUE.
1632	ENDDO
1633	!
1634	!--setting up ok_4xCO2atm to TRUE in XIOS case
1635	IF (xios_field_is_active("rsut4co2").OR.xios_field_is_active("rlut4co2").OR. &
1636	xios_field_is_active("rsutcs4co2").OR.xios_field_is_active("rlutcs4co2").OR. &
1637	xios_field_is_active("rsu4co2").OR.xios_field_is_active("rsucs4co2").OR. &
1638	xios_field_is_active("rsd4co2").OR.xios_field_is_active("rsdcs4co2").OR. &
1639	xios_field_is_active("rlu4co2").OR.xios_field_is_active("rlucs4co2").OR. &
1640	xios_field_is_active("rld4co2").OR.xios_field_is_active("rldcs4co2")) &
1641	ok_4xCO2atm=.TRUE.
1642	#endif
1643
1644	#ifdef CPP_IOIPSL
1645
1646	!$OMP MASTER
1647	! FH : if ok_sync=.true. , the time axis is written at each time step
1648	! in the output files. Only at the end in the opposite case
1649	ok_sync_omp=.false.
1650	CALL getin('ok_sync',ok_sync_omp)
1651	CALL phys_output_open(longitude_deg,latitude_deg,nCFMIP,tabijGCM, &
1652	iGCM,jGCM,lonGCM,latGCM, &
1653	jjmp1,nlevSTD,clevSTD,rlevSTD, phys_tstep,ok_veget, &
1654	type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1655	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1656	read_climoz, phys_out_filestations, &
1657	new_aod, aerosol_couple, &
1658	flag_aerosol_strat, pdtphys, paprs, pphis, &
1659	pplay, lmax_th, ptconv, ptconvth, ivap, &
1660	d_u, d_t, qx, d_qx, zmasse, ok_sync_omp)
1661	!$OMP END MASTER
1662	!$OMP BARRIER
1663	ok_sync=ok_sync_omp
1664
1665	freq_outNMC(1) = ecrit_files(7)
1666	freq_outNMC(2) = ecrit_files(8)
1667	freq_outNMC(3) = ecrit_files(9)
1668	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1669	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1670	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1671
1672	#ifndef CPP_XIOS
1673	CALL ini_paramLMDZ_phy(phys_tstep,nid_ctesGCM)
1674	#endif
1675
1676	#endif
1677	ecrit_reg = ecrit_reg * un_jour
1678	ecrit_tra = ecrit_tra * un_jour
1679
1680	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
1681	date0 = jD_ref
1682	WRITE(,) 'physiq date0 : ',date0
1683	!
1684	!
1685	!
1686	! Prescrire l'ozone dans l'atmosphere
1687	!
1688	!
1689	!c DO i = 1, klon
1690	!c DO k = 1, klev
1691	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1692	!c ENDDO
1693	!c ENDDO
1694	!
1695	IF (type_trac == 'inca') THEN
1696	#ifdef INCA
1697	CALL VTe(VTphysiq)
1698	CALL VTb(VTinca)
1699	calday = REAL(days_elapsed) + jH_cur
1700	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1701
1702	CALL chemini( &
1703	rg, &
1704	ra, &
1705	cell_area, &
1706	latitude_deg, &
1707	longitude_deg, &
1708	presnivs, &
1709	calday, &
1710	klon, &
1711	nqtot, &
1712	nqo, &
1713	pdtphys, &
1714	annee_ref, &
1715	year_cur, &
1716	day_ref, &
1717	day_ini, &
1718	start_time, &
1719	itau_phy, &
1720	date0, &
1721	io_lon, &
1722	io_lat)
1723
1724	CALL VTe(VTinca)
1725	CALL VTb(VTphysiq)
1726	#endif
1727	ENDIF
1728	!
1729	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1730	! Nouvelle initialisation pour le rayonnement RRTM
1731	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1732
1733	CALL iniradia(klon,klev,paprs(1,1:klev+1))
1734
1735	!$omp single
1736	IF (read_climoz >= 1) CALL open_climoz(ncid_climoz, press_cen_climoz, &
1737	press_edg_climoz, time_climoz, ok_daily_climoz, adjust_tropopause)
1738	!$omp end single
1739	!
1740	!IM betaCRF
1741	pfree=70000. !Pa
1742	beta_pbl=1.
1743	beta_free=1.
1744	lon1_beta=-180.
1745	lon2_beta=+180.
1746	lat1_beta=90.
1747	lat2_beta=-90.
1748	mskocean_beta=.FALSE.
1749
1750	!albedo SB >>>
1751	select case(nsw)
1752	case(2)
1753	SFRWL(1)=0.45538747
1754	SFRWL(2)=0.54461211
1755	case(4)
1756	SFRWL(1)=0.45538747
1757	SFRWL(2)=0.32870591
1758	SFRWL(3)=0.18568763
1759	SFRWL(4)=3.02191470E-02
1760	case(6)
1761	SFRWL(1)=1.28432794E-03
1762	SFRWL(2)=0.12304168
1763	SFRWL(3)=0.33106142
1764	SFRWL(4)=0.32870591
1765	SFRWL(5)=0.18568763
1766	SFRWL(6)=3.02191470E-02
1767	end select
1768
1769
1770	!albedo SB <<<
1771
1772	OPEN(99,file='beta_crf.data',status='old', &
1773	form='formatted',err=9999)
1774	READ(99,*,end=9998) pfree
1775	READ(99,*,end=9998) beta_pbl
1776	READ(99,*,end=9998) beta_free
1777	READ(99,*,end=9998) lon1_beta
1778	READ(99,*,end=9998) lon2_beta
1779	READ(99,*,end=9998) lat1_beta
1780	READ(99,*,end=9998) lat2_beta
1781	READ(99,*,end=9998) mskocean_beta
1782	9998 Continue
1783	CLOSE(99)
1784	9999 Continue
1785	WRITE(,)'pfree=',pfree
1786	WRITE(,)'beta_pbl=',beta_pbl
1787	WRITE(,)'beta_free=',beta_free
1788	WRITE(,)'lon1_beta=',lon1_beta
1789	WRITE(,)'lon2_beta=',lon2_beta
1790	WRITE(,)'lat1_beta=',lat1_beta
1791	WRITE(,)'lat2_beta=',lat2_beta
1792	WRITE(,)'mskocean_beta=',mskocean_beta
1793	ENDIF
1794	!
1795	! ************** Fin de IF ( debut ) *************
1796	!
1797	!
1798	! Incrementer le compteur de la physique
1799	!
1800	itap = itap + 1
1801	IF (is_master .OR. prt_level > 9) THEN
1802	IF (prt_level > 5 .or. MOD(itap,5) == 0) THEN
1803	WRITE(LUNOUT,*)'Entering physics elapsed seconds since start ', current_time
1804	WRITE(LUNOUT,100)year_cur,mth_cur,day_cur,hour/3600.
1805	100 FORMAT('Date = ',i4.4,' / ',i2.2, ' / ',i2.2,' : ',f20.17)
1806	ENDIF
1807	ENDIF
1808	!
1809	!
1810	! Update fraction of the sub-surfaces (pctsrf) and
1811	! initialize, where a new fraction has appeared, all variables depending
1812	! on the surface fraction.
1813	!
1814	CALL change_srf_frac(itap, phys_tstep, days_elapsed+1, &
1815	pctsrf, fevap, z0m, z0h, agesno, &
1816	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
1817
1818	! Update time and other variables in Reprobus
1819	IF (type_trac == 'repr') THEN
1820	#ifdef REPROBUS
1821	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1822	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1823	CALL Rtime(debut)
1824	#endif
1825	ENDIF
1826
1827
1828	! Tendances bidons pour les processus qui n'affectent pas certaines
1829	! variables.
1830	du0(:,:)=0.
1831	dv0(:,:)=0.
1832	dt0 = 0.
1833	dq0(:,:)=0.
1834	dql0(:,:)=0.
1835	dqi0(:,:)=0.
1836	dsig0(:) = 0.
1837	ddens0(:) = 0.
1838	wkoccur1(:)=1
1839	!
1840	! Mettre a zero des variables de sortie (pour securite)
1841	!
1842	DO i = 1, klon
1843	d_ps(i) = 0.0
1844	ENDDO
1845	DO k = 1, klev
1846	DO i = 1, klon
1847	d_t(i,k) = 0.0
1848	d_u(i,k) = 0.0
1849	d_v(i,k) = 0.0
1850	ENDDO
1851	ENDDO
1852	DO iq = 1, nqtot
1853	DO k = 1, klev
1854	DO i = 1, klon
1855	d_qx(i,k,iq) = 0.0
1856	ENDDO
1857	ENDDO
1858	ENDDO
1859	da(:,:)=0.
1860	mp(:,:)=0.
1861	phi(:,:,:)=0.
1862	! RomP >>>
1863	phi2(:,:,:)=0.
1864	beta_prec_fisrt(:,:)=0.
1865	beta_prec(:,:)=0.
1866	epmlmMm(:,:,:)=0.
1867	eplaMm(:,:)=0.
1868	d1a(:,:)=0.
1869	dam(:,:)=0.
1870	pmflxr=0.
1871	pmflxs=0.
1872	! RomP <<<
1873
1874	!
1875	! Ne pas affecter les valeurs entrees de u, v, h, et q
1876	!
1877	DO k = 1, klev
1878	DO i = 1, klon
1879	t_seri(i,k) = t(i,k)
1880	u_seri(i,k) = u(i,k)
1881	v_seri(i,k) = v(i,k)
1882	q_seri(i,k) = qx(i,k,ivap)
1883	ql_seri(i,k) = qx(i,k,iliq)
1884	!CR: ATTENTION, on rajoute la variable glace
1885	IF (nqo.eq.2) THEN
1886	qs_seri(i,k) = 0.
1887	ELSE IF (nqo.eq.3) THEN
1888	qs_seri(i,k) = qx(i,k,isol)
1889	ENDIF
1890	ENDDO
1891	ENDDO
1892	!
1893	!--OB mass fixer
1894	IF (mass_fixer) THEN
1895	!--store initial water burden
1896	qql1(:)=0.0
1897	DO k = 1, klev
1898	qql1(:)=qql1(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
1899	ENDDO
1900	ENDIF
1901	!--fin mass fixer
1902
1903	tke0(:,:)=pbl_tke(:,:,is_ave)
1904	!CR:Nombre de traceurs de l'eau: nqo
1905	! IF (nqtot.GE.3) THEN
1906	IF (nqtot.GE.(nqo+1)) THEN
1907	! DO iq = 3, nqtot
1908	DO iq = nqo+1, nqtot
1909	DO k = 1, klev
1910	DO i = 1, klon
1911	! tr_seri(i,k,iq-2) = qx(i,k,iq)
1912	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
1913	ENDDO
1914	ENDDO
1915	ENDDO
1916	ELSE
1917	DO k = 1, klev
1918	DO i = 1, klon
1919	tr_seri(i,k,1) = 0.0
1920	ENDDO
1921	ENDDO
1922	ENDIF
1923	!
1924	DO i = 1, klon
1925	ztsol(i) = 0.
1926	ENDDO
1927	DO nsrf = 1, nbsrf
1928	DO i = 1, klon
1929	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1930	ENDDO
1931	ENDDO
1932	! Initialize variables used for diagnostic purpose
1933	IF (flag_inhib_tend .ne. 0) CALL init_cmp_seri
1934
1935	! Diagnostiquer la tendance dynamique
1936	!
1937	IF (ancien_ok) THEN
1938	!
1939	d_u_dyn(:,:) = (u_seri(:,:)-u_ancien(:,:))/phys_tstep
1940	d_v_dyn(:,:) = (v_seri(:,:)-v_ancien(:,:))/phys_tstep
1941	d_t_dyn(:,:) = (t_seri(:,:)-t_ancien(:,:))/phys_tstep
1942	d_q_dyn(:,:) = (q_seri(:,:)-q_ancien(:,:))/phys_tstep
1943	d_ql_dyn(:,:) = (ql_seri(:,:)-ql_ancien(:,:))/phys_tstep
1944	d_qs_dyn(:,:) = (qs_seri(:,:)-qs_ancien(:,:))/phys_tstep
1945	CALL water_int(klon,klev,q_seri,zmasse,zx_tmp_fi2d)
1946	d_q_dyn2d(:)=(zx_tmp_fi2d(:)-prw_ancien(:))/phys_tstep
1947	CALL water_int(klon,klev,ql_seri,zmasse,zx_tmp_fi2d)
1948	d_ql_dyn2d(:)=(zx_tmp_fi2d(:)-prlw_ancien(:))/phys_tstep
1949	CALL water_int(klon,klev,qs_seri,zmasse,zx_tmp_fi2d)
1950	d_qs_dyn2d(:)=(zx_tmp_fi2d(:)-prsw_ancien(:))/phys_tstep
1951	! !! RomP >>> td dyn traceur
1952	IF (nqtot.GT.nqo) THEN ! jyg
1953	DO iq = nqo+1, nqtot ! jyg
1954	d_tr_dyn(:,:,iq-nqo)=(tr_seri(:,:,iq-nqo)-tr_ancien(:,:,iq-nqo))/phys_tstep ! jyg
1955	ENDDO
1956	ENDIF
1957	! !! RomP <<<
1958	ELSE
1959	d_u_dyn(:,:) = 0.0
1960	d_v_dyn(:,:) = 0.0
1961	d_t_dyn(:,:) = 0.0
1962	d_q_dyn(:,:) = 0.0
1963	d_ql_dyn(:,:) = 0.0
1964	d_qs_dyn(:,:) = 0.0
1965	d_q_dyn2d(:) = 0.0
1966	d_ql_dyn2d(:) = 0.0
1967	d_qs_dyn2d(:) = 0.0
1968	! !! RomP >>> td dyn traceur
1969	IF (nqtot.GT.nqo) THEN ! jyg
1970	DO iq = nqo+1, nqtot ! jyg
1971	d_tr_dyn(:,:,iq-nqo)= 0.0 ! jyg
1972	ENDDO
1973	ENDIF
1974	! !! RomP <<<
1975	ancien_ok = .TRUE.
1976	ENDIF
1977	!
1978	! Ajouter le geopotentiel du sol:
1979	!
1980	DO k = 1, klev
1981	DO i = 1, klon
1982	zphi(i,k) = pphi(i,k) + pphis(i)
1983	ENDDO
1984	ENDDO
1985	!
1986	! Verifier les temperatures
1987	!
1988	!IM BEG
1989	IF (check) THEN
1990	amn=MIN(ftsol(1,is_ter),1000.)
1991	amx=MAX(ftsol(1,is_ter),-1000.)
1992	DO i=2, klon
1993	amn=MIN(ftsol(i,is_ter),amn)
1994	amx=MAX(ftsol(i,is_ter),amx)
1995	ENDDO
1996	!
1997	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1998	ENDIF !(check) THEN
1999	!IM END
2000	!
2001	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
2002	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
2003
2004	!
2005	!IM BEG
2006	IF (check) THEN
2007	amn=MIN(ftsol(1,is_ter),1000.)
2008	amx=MAX(ftsol(1,is_ter),-1000.)
2009	DO i=2, klon
2010	amn=MIN(ftsol(i,is_ter),amn)
2011	amx=MAX(ftsol(i,is_ter),amx)
2012	ENDDO
2013	!
2014	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
2015	ENDIF !(check) THEN
2016	!IM END
2017	!
2018	! Mettre en action les conditions aux limites (albedo, sst, etc.).
2019	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
2020	!
2021	! Update ozone if day change
2022	IF (MOD(itap-1,lmt_pas) == 0) THEN
2023	IF (read_climoz <= 0) THEN
2024	! Once per day, update ozone from Royer:
2025	IF (solarlong0<-999.) then
2026	! Generic case with evolvoing season
2027	zzz=real(days_elapsed+1)
2028	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
2029	! Particular case with annual mean insolation
2030	zzz=real(90) ! could be revisited
2031	IF (read_climoz/=-1) THEN
2032	abort_message ='read_climoz=-1 is recommended when ' &
2033	// 'solarlong0=1000.'
2034	CALL abort_physic (modname,abort_message,1)
2035	ENDIF
2036	ELSE
2037	! Case where the season is imposed with solarlong0
2038	zzz=real(90) ! could be revisited
2039	ENDIF
2040
2041	wo(:,:,1)=ozonecm(latitude_deg, paprs,read_climoz,rjour=zzz)
2042	ELSE
2043	!--- ro3i = elapsed days number since current year 1st january, 0h
2044	ro3i=days_elapsed+jh_cur-jh_1jan
2045	!--- scaling for old style files (360 records)
2046	IF(SIZE(time_climoz)==360.AND..NOT.ok_daily_climoz) ro3i=ro3i*360./year_len
2047	IF(adjust_tropopause) THEN
2048	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2049	ro3i, press_edg_climoz, paprs, wo, time_climoz, &
2050	longitude_deg, latitude_deg, press_cen_climoz, &
2051	dyn_tropopause(t_seri, ztsol, paprs, pplay, rot))
2052	ELSE
2053	CALL regr_pr_time_av(ncid_climoz, vars_climoz(1:read_climoz), &
2054	ro3i, press_edg_climoz, paprs, wo, time_climoz)
2055	END IF
2056	! Convert from mole fraction of ozone to column density of ozone in a
2057	! cell, in kDU:
2058	FORALL (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
2059	* zmasse / dobson_u / 1e3
2060	! (By regridding ozone values for LMDZ only once a day, we
2061	! have already neglected the variation of pressure in one
2062	! day. So do not recompute "wo" at each time step even if
2063	! "zmasse" changes a little.)
2064	ENDIF
2065	ENDIF
2066	!
2067	! Re-evaporer l'eau liquide nuageuse
2068	!
2069	CALL reevap (klon,klev,iflag_ice_thermo,t_seri,q_seri,ql_seri,qs_seri, &
2070	& d_t_eva,d_q_eva,d_ql_eva,d_qi_eva)
2071
2072	CALL add_phys_tend &
2073	(du0,dv0,d_t_eva,d_q_eva,d_ql_eva,d_qi_eva,paprs,&
2074	'eva',abortphy,flag_inhib_tend,itap,0)
2075	call prt_enerbil('eva',itap)
2076
2077	!=========================================================================
2078	! Calculs de l'orbite.
2079	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
2080	! doit donc etre plac\'e avant radlwsw et pbl_surface
2081
2082	! !! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2083	CALL ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
2084	day_since_equinox = (jD_cur + jH_cur) - jD_eq
2085	!
2086	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
2087	! solarlong0
2088	IF (solarlong0<-999.) THEN
2089	IF (new_orbit) THEN
2090	! calcul selon la routine utilisee pour les planetes
2091	CALL solarlong(day_since_equinox, zlongi, dist)
2092	ELSE
2093	! calcul selon la routine utilisee pour l'AR4
2094	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
2095	ENDIF
2096	ELSE
2097	zlongi=solarlong0 ! longitude solaire vraie
2098	dist=1. ! distance au soleil / moyenne
2099	ENDIF
2100
2101	IF (prt_level.ge.1) write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
2102
2103
2104	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2105	! Calcul de l'ensoleillement :
2106	! ============================
2107	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
2108	! l'annee a partir d'une formule analytique.
2109	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
2110	! non nul aux poles.
2111	IF (abs(solarlong0-1000.)<1.e-4) THEN
2112	CALL zenang_an(iflag_cycle_diurne.GE.1,jH_cur, &
2113	latitude_deg,longitude_deg,rmu0,fract)
2114	swradcorr(:) = 1.0
2115	JrNt(:) = 1.0
2116	zrmu0(:) = rmu0(:)
2117	ELSE
2118	! recode par Olivier Boucher en sept 2015
2119	SELECT CASE (iflag_cycle_diurne)
2120	CASE(0)
2121	! Sans cycle diurne
2122	CALL angle(zlongi, latitude_deg, fract, rmu0)
2123	swradcorr = 1.0
2124	JrNt = 1.0
2125	zrmu0 = rmu0
2126	CASE(1)
2127	! Avec cycle diurne sans application des poids
2128	! bit comparable a l ancienne formulation cycle_diurne=true
2129	! on integre entre gmtime et gmtime+radpas
2130	zdtime=phys_tstep*REAL(radpas) ! pas de temps du rayonnement (s)
2131	CALL zenang(zlongi,jH_cur,0.0,zdtime, &
2132	latitude_deg,longitude_deg,rmu0,fract)
2133	zrmu0 = rmu0
2134	swradcorr = 1.0
2135	! Calcul du flag jour-nuit
2136	JrNt = 0.0
2137	WHERE (fract.GT.0.0) JrNt = 1.0
2138	CASE(2)
2139	! Avec cycle diurne sans application des poids
2140	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
2141	! Comme cette routine est appele a tous les pas de temps de
2142	! la physique meme si le rayonnement n'est pas appele je
2143	! remonte en arriere les radpas-1 pas de temps
2144	! suivant. Petite ruse avec MOD pour prendre en compte le
2145	! premier pas de temps de la physique pendant lequel
2146	! itaprad=0
2147	zdtime1=phys_tstep*REAL(-MOD(itaprad,radpas)-1)
2148	zdtime2=phys_tstep*REAL(radpas-MOD(itaprad,radpas)-1)
2149	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2150	latitude_deg,longitude_deg,rmu0,fract)
2151	!
2152	! Calcul des poids
2153	!
2154	zdtime1=-phys_tstep !--on corrige le rayonnement pour representer le
2155	zdtime2=0.0 !--pas de temps de la physique qui se termine
2156	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2, &
2157	latitude_deg,longitude_deg,zrmu0,zfract)
2158	swradcorr = 0.0
2159	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) &
2160	swradcorr=zfract/fract*zrmu0/rmu0
2161	! Calcul du flag jour-nuit
2162	JrNt = 0.0
2163	WHERE (zfract.GT.0.0) JrNt = 1.0
2164	END SELECT
2165	ENDIF
2166
2167	IF (mydebug) THEN
2168	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2169	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2170	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2171	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2172	ENDIF
2173
2174	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2175	! Appel au pbl_surface : Planetary Boudary Layer et Surface
2176	! Cela implique tous les interactions des sous-surfaces et la
2177	! partie diffusion turbulent du couche limit.
2178	!
2179	! Certains varibales de sorties de pbl_surface sont utiliser que pour
2180	! ecriture des fihiers hist_XXXX.nc, ces sont :
2181	! qsol, zq2m, s_pblh, s_lcl,
2182	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
2183	! s_therm, s_trmb1, s_trmb2, s_trmb3,
2184	! zu10m, zv10m, fder,
2185	! zxqsurf, rh2m, zxfluxu, zxfluxv,
2186	! frugs, agesno, fsollw, fsolsw,
2187	! d_ts, fevap, fluxlat, t2m,
2188	! wfbils, wfbilo, fluxt, fluxu, fluxv,
2189	!
2190	! Certains ne sont pas utiliser du tout :
2191	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
2192	!
2193
2194	! Calcul de l'humidite de saturation au niveau du sol
2195
2196
2197
2198	IF (iflag_pbl/=0) THEN
2199
2200	!jyg+nrlmd<
2201	!!jyg IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
2202	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,10) .ge. 1) THEN
2203	print *,'debut du splitting de la PBL'
2204	ENDIF
2205	! !!
2206	!>jyg+nrlmd
2207	!
2208	!-------gustiness calculation-------!
2209	IF (iflag_gusts==0) THEN
2210	gustiness(1:klon)=0
2211	ELSE IF (iflag_gusts==1) THEN
2212	do i = 1, klon
2213	gustiness(i)=f_gust_blale_bl(i)+f_gust_wkale_wake(i)
2214	enddo
2215	! ELSE IF (iflag_gusts==2) THEN
2216	! do i = 1, klon
2217	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk&
2218	! *ale_wake(i) !! need to make sigma_wk accessible here
2219	! enddo
2220	! ELSE IF (iflag_gusts==3) THEN
2221	! do i = 1, klon
2222	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
2223	! enddo
2224	ENDIF
2225
2226
2227
2228	CALL pbl_surface( &
2229	phys_tstep, date0, itap, days_elapsed+1, &
2230	debut, lafin, &
2231	longitude_deg, latitude_deg, rugoro, zrmu0, &
2232	zsig, sollwdown, pphi, cldt, &
2233	rain_fall, snow_fall, solsw, sollw, &
2234	gustiness, &
2235	t_seri, q_seri, u_seri, v_seri, &
2236	!nrlmd+jyg<
2237	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
2238	!>nrlmd+jyg
2239	pplay, paprs, pctsrf, &
2240	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
2241	!albedo SB <<<
2242	cdragh, cdragm, u1, v1, &
2243	!albedo SB >>>
2244	! albsol1, albsol2, sens, evap, &
2245	albsol_dir, albsol_dif, sens, evap, &
2246	!albedo SB <<<
2247	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
2248	zxtsol, zxfluxlat, zt2m, qsat2m, &
2249	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
2250	!nrlmd<
2251	!jyg<
2252	d_t_vdf_w, d_q_vdf_w, &
2253	d_t_vdf_x, d_q_vdf_x, &
2254	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
2255	!>jyg
2256	delta_tsurf,wake_dens, &
2257	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
2258	kh,kh_x,kh_w, &
2259	!>nrlmd
2260	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
2261	slab_wfbils, &
2262	qsol, zq2m, s_pblh, s_lcl, &
2263	!jyg<
2264	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
2265	!>jyg
2266	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
2267	s_therm, s_trmb1, s_trmb2, s_trmb3, &
2268	zustar, zu10m, zv10m, fder, &
2269	zxqsurf, rh2m, zxfluxu, zxfluxv, &
2270	z0m, z0h, agesno, fsollw, fsolsw, &
2271	d_ts, fevap, fluxlat, t2m, &
2272	wfbils, wfbilo, wfevap, wfrain, wfsnow, &
2273	fluxt, fluxu, fluxv, &
2274	dsens, devap, zxsnow, &
2275	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
2276	!nrlmd+jyg<
2277	wake_delta_pbl_TKE, &
2278	!>nrlmd+jyg
2279	treedrg )
2280	!FC
2281	!
2282	! Add turbulent diffusion tendency to the wake difference variables
2283	!!jyg IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2284	IF (mod(iflag_pbl_split,10) .NE. 0) THEN
2285	!jyg<
2286	d_deltat_vdf(:,:) = d_t_vdf_w(:,:)-d_t_vdf_x(:,:)
2287	d_deltaq_vdf(:,:) = d_q_vdf_w(:,:)-d_q_vdf_x(:,:)
2288	CALL add_wake_tend &
2289	(d_deltat_vdf, d_deltaq_vdf, dsig0, ddens0, wkoccur1, 'vdf', abortphy)
2290	ELSE
2291	d_deltat_vdf(:,:) = 0.
2292	d_deltaq_vdf(:,:) = 0.
2293	!>jyg
2294	ENDIF
2295
2296
2297
2298
2299
2300	!---------------------------------------------------------------------
2301	! ajout des tendances de la diffusion turbulente
2302	IF (klon_glo==1) THEN
2303	CALL add_pbl_tend &
2304	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2305	'vdf',abortphy,flag_inhib_tend,itap)
2306	ELSE
2307	CALL add_phys_tend &
2308	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,&
2309	'vdf',abortphy,flag_inhib_tend,itap,0)
2310	ENDIF
2311	call prt_enerbil('vdf',itap)
2312	!--------------------------------------------------------------------
2313
2314	IF (mydebug) THEN
2315	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2316	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2317	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2318	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2319	ENDIF
2320
2321	!albedo SB >>>
2322	albsol1=0.
2323	albsol2=0.
2324	falb1=0.
2325	falb2=0.
2326	SELECT CASE(nsw)
2327	CASE(2)
2328	albsol1=albsol_dir(:,1)
2329	albsol2=albsol_dir(:,2)
2330	falb1=falb_dir(:,1,:)
2331	falb2=falb_dir(:,2,:)
2332	CASE(4)
2333	albsol1=albsol_dir(:,1)
2334	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3) &
2335	+albsol_dir(:,4)*SFRWL(4)
2336	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2337	falb1=falb_dir(:,1,:)
2338	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3) &
2339	+falb_dir(:,4,:)*SFRWL(4)
2340	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2341	CASE(6)
2342	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2) &
2343	+albsol_dir(:,3)*SFRWL(3)
2344	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2345	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5) &
2346	+albsol_dir(:,6)*SFRWL(6)
2347	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2348	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2) &
2349	+falb_dir(:,3,:)*SFRWL(3)
2350	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2351	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5) &
2352	+falb_dir(:,6,:)*SFRWL(6)
2353	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2354	END SELECt
2355	!albedo SB <<<
2356
2357
2358	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2359	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2360
2361	ENDIF
2362	! =================================================================== c
2363	! Calcul de Qsat
2364
2365	DO k = 1, klev
2366	DO i = 1, klon
2367	zx_t = t_seri(i,k)
2368	IF (thermcep) THEN
2369	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2370	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2371	zx_qs = MIN(0.5,zx_qs)
2372	zcor = 1./(1.-retv*zx_qs)
2373	zx_qs = zx_qs*zcor
2374	ELSE
2375	!! IF (zx_t.LT.t_coup) THEN !jyg
2376	IF (zx_t.LT.rtt) THEN !jyg
2377	zx_qs = qsats(zx_t)/pplay(i,k)
2378	ELSE
2379	zx_qs = qsatl(zx_t)/pplay(i,k)
2380	ENDIF
2381	ENDIF
2382	zqsat(i,k)=zx_qs
2383	ENDDO
2384	ENDDO
2385
2386	IF (prt_level.ge.1) THEN
2387	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2388	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2389	ENDIF
2390	!
2391	! Appeler la convection (au choix)
2392	!
2393	DO k = 1, klev
2394	DO i = 1, klon
2395	conv_q(i,k) = d_q_dyn(i,k) &
2396	+ d_q_vdf(i,k)/phys_tstep
2397	conv_t(i,k) = d_t_dyn(i,k) &
2398	+ d_t_vdf(i,k)/phys_tstep
2399	ENDDO
2400	ENDDO
2401	IF (check) THEN
2402	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2403	WRITE(lunout,*) "avantcon=", za
2404	ENDIF
2405	zx_ajustq = .FALSE.
2406	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2407	IF (zx_ajustq) THEN
2408	DO i = 1, klon
2409	z_avant(i) = 0.0
2410	ENDDO
2411	DO k = 1, klev
2412	DO i = 1, klon
2413	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2414	*(paprs(i,k)-paprs(i,k+1))/RG
2415	ENDDO
2416	ENDDO
2417	ENDIF
2418
2419	! Calcule de vitesse verticale a partir de flux de masse verticale
2420	DO k = 1, klev
2421	DO i = 1, klon
2422	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
2423	ENDDO
2424	ENDDO
2425
2426	IF (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2427	omega(igout, :)
2428	!
2429	! Appel de la convection tous les "cvpas"
2430	!
2431	IF (MOD(itapcv,cvpas).EQ.0) THEN
2432
2433	IF (iflag_con.EQ.1) THEN
2434	abort_message ='reactiver le call conlmd dans physiq.F'
2435	CALL abort_physic (modname,abort_message,1)
2436	! CALL conlmd (phys_tstep, paprs, pplay, t_seri, q_seri, conv_q,
2437	! . d_t_con, d_q_con,
2438	! . rain_con, snow_con, ibas_con, itop_con)
2439	ELSE IF (iflag_con.EQ.2) THEN
2440	CALL conflx(phys_tstep, paprs, pplay, t_seri, q_seri, &
2441	conv_t, conv_q, -evap, omega, &
2442	d_t_con, d_q_con, rain_con, snow_con, &
2443	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2444	kcbot, kctop, kdtop, pmflxr, pmflxs)
2445	d_u_con = 0.
2446	d_v_con = 0.
2447
2448	WHERE (rain_con < 0.) rain_con = 0.
2449	WHERE (snow_con < 0.) snow_con = 0.
2450	DO i = 1, klon
2451	ibas_con(i) = klev+1 - kcbot(i)
2452	itop_con(i) = klev+1 - kctop(i)
2453	ENDDO
2454	ELSE IF (iflag_con.GE.3) THEN
2455	! nb of tracers for the KE convection:
2456	! MAF la partie traceurs est faite dans phytrac
2457	! on met ntra=1 pour limiter les appels mais on peut
2458	! supprimer les calculs / ftra.
2459	ntra = 1
2460
2461	!=======================================================================
2462	!ajout pour la parametrisation des poches froides: calcul de
2463	!t_w et t_x: si pas de poches froides, t_w=t_x=t_seri
2464	IF (iflag_wake>=1) THEN
2465	DO k=1,klev
2466	DO i=1,klon
2467	t_w(i,k) = t_seri(i,k) + (1-wake_s(i))*wake_deltat(i,k)
2468	q_w(i,k) = q_seri(i,k) + (1-wake_s(i))*wake_deltaq(i,k)
2469	t_x(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2470	q_x(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2471	ENDDO
2472	ENDDO
2473	ELSE
2474	t_w(:,:) = t_seri(:,:)
2475	q_w(:,:) = q_seri(:,:)
2476	t_x(:,:) = t_seri(:,:)
2477	q_x(:,:) = q_seri(:,:)
2478	ENDIF
2479	!
2480	!jyg<
2481	! Perform dry adiabatic adjustment on wake profile
2482	! The corresponding tendencies are added to the convective tendencies
2483	! after the call to the convective scheme.
2484	IF (iflag_wake>=1) then
2485	IF (iflag_adjwk >= 1) THEN
2486	limbas(:) = 1
2487	CALL ajsec(paprs, pplay, t_w, q_w, limbas, &
2488	d_t_adjwk, d_q_adjwk)
2489	!
2490	DO k=1,klev
2491	DO i=1,klon
2492	IF (wake_s(i) .GT. 1.e-3) THEN
2493	t_w(i,k) = t_w(i,k) + d_t_adjwk(i,k)
2494	q_w(i,k) = q_w(i,k) + d_q_adjwk(i,k)
2495	d_deltat_ajs_cv(i,k) = d_t_adjwk(i,k)
2496	d_deltaq_ajs_cv(i,k) = d_q_adjwk(i,k)
2497	ELSE
2498	d_deltat_ajs_cv(i,k) = 0.
2499	d_deltaq_ajs_cv(i,k) = 0.
2500	ENDIF
2501	ENDDO
2502	ENDDO
2503	IF (iflag_adjwk == 2) THEN
2504	CALL add_wake_tend &
2505	(d_deltat_ajs_cv, d_deltaq_ajs_cv, dsig0, ddens0, wkoccur1, 'ajs_cv', abortphy)
2506	ENDIF ! (iflag_adjwk == 2)
2507	ENDIF ! (iflag_adjwk >= 1)
2508	ENDIF ! (iflag_wake>=1)
2509	!>jyg
2510	!
2511
2512	!! print *,'physiq. q_w(1,k), q_x(1,k) ', &
2513	!! (k, q_w(1,k), q_x(1,k),k=1,25)
2514
2515	!jyg<
2516	CALL alpale( debut, itap, phys_tstep, paprs, omega, t_seri, &
2517	alp_offset, it_wape_prescr, wape_prescr, fip_prescr, &
2518	ale_bl_prescr, alp_bl_prescr, &
2519	wake_pe, wake_fip, &
2520	Ale_bl, Ale_bl_trig, Alp_bl, &
2521	Ale, Alp , Ale_wake, Alp_wake)
2522	!>jyg
2523	!
2524	! sb, oct02:
2525	! Schema de convection modularise et vectorise:
2526	! (driver commun aux versions 3 et 4)
2527	!
2528	IF (ok_cvl) THEN ! new driver for convectL
2529	!
2530	!jyg<
2531	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2532	! Calculate the upmost level of deep convection loops: k_upper_cv
2533	! (near 22 km)
2534	izero = klon/2+1/klon
2535	k_upper_cv = klev
2536	DO k = klev,1,-1
2537	IF (pphi(izero,k) > 22.e4) k_upper_cv = k
2538	ENDDO
2539	IF (prt_level .ge. 5) THEN
2540	Print *, 'upmost level of deep convection loops: k_upper_cv = ', &
2541	k_upper_cv
2542	ENDIF
2543	!
2544	!>jyg
2545	IF (type_trac == 'repr') THEN
2546	nbtr_tmp=ntra
2547	ELSE
2548	nbtr_tmp=nbtr
2549	ENDIF
2550	!jyg iflag_con est dans clesphys
2551	!c CALL concvl (iflag_con,iflag_clos,
2552	CALL concvl (iflag_clos, &
2553	phys_tstep, paprs, pplay, k_upper_cv, t_x,q_x, &
2554	t_w,q_w,wake_s, &
2555	u_seri,v_seri,tr_seri,nbtr_tmp, &
2556	ALE,ALP, &
2557	sig1,w01, &
2558	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2559	rain_con, snow_con, ibas_con, itop_con, sigd, &
2560	ema_cbmf,plcl,plfc,wbeff,convoccur,upwd,dnwd,dnwd0, &
2561	Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2562	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2563	! RomP >>>
2564	!! . pmflxr,pmflxs,da,phi,mp,
2565	!! . ftd,fqd,lalim_conv,wght_th)
2566	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2567	ftd,fqd,lalim_conv,wght_th, &
2568	ev, ep,epmlmMm,eplaMm, &
2569	wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2570	tau_cld_cv,coefw_cld_cv,epmax_diag)
2571
2572	! RomP <<<
2573
2574	!IM begin
2575	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2576	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2577	!IM end
2578	!IM cf. FH
2579	clwcon0=qcondc
2580	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2581
2582	DO i = 1, klon
2583	IF (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2584	ENDDO
2585	!
2586	!jyg<
2587	! Add the tendency due to the dry adjustment of the wake profile
2588	IF (iflag_wake>=1) THEN
2589	IF (iflag_adjwk == 2) THEN
2590	DO k=1,klev
2591	DO i=1,klon
2592	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/phys_tstep
2593	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/phys_tstep
2594	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
2595	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
2596	ENDDO
2597	ENDDO
2598	ENDIF ! (iflag_adjwk = 2)
2599	ENDIF ! (iflag_wake>=1)
2600	!>jyg
2601	!
2602	ELSE ! ok_cvl
2603
2604	! MAF conema3 ne contient pas les traceurs
2605	CALL conema3 (phys_tstep, &
2606	paprs,pplay,t_seri,q_seri, &
2607	u_seri,v_seri,tr_seri,ntra, &
2608	sig1,w01, &
2609	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2610	rain_con, snow_con, ibas_con, itop_con, &
2611	upwd,dnwd,dnwd0,bas,top, &
2612	Ma,cape,tvp,rflag, &
2613	pbase &
2614	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2615	,clwcon0)
2616
2617	ENDIF ! ok_cvl
2618
2619	!
2620	! Correction precip
2621	rain_con = rain_con * cvl_corr
2622	snow_con = snow_con * cvl_corr
2623	!
2624
2625	IF (.NOT. ok_gust) THEN
2626	do i = 1, klon
2627	wd(i)=0.0
2628	enddo
2629	ENDIF
2630
2631	! =================================================================== c
2632	! Calcul des proprietes des nuages convectifs
2633	!
2634
2635	! calcul des proprietes des nuages convectifs
2636	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2637	IF (iflag_cld_cv == 0) THEN
2638	CALL clouds_gno &
2639	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2640	ELSE
2641	CALL clouds_bigauss &
2642	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
2643	ENDIF
2644
2645
2646	! =================================================================== c
2647
2648	DO i = 1, klon
2649	itop_con(i) = min(max(itop_con(i),1),klev)
2650	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2651	ENDDO
2652
2653	DO i = 1, klon
2654	ema_pcb(i) = paprs(i,ibas_con(i))
2655	ENDDO
2656	DO i = 1, klon
2657	! L'idicage de itop_con peut cacher un pb potentiel
2658	! FH sous la dictee de JYG, CR
2659	ema_pct(i) = paprs(i,itop_con(i)+1)
2660
2661	IF (itop_con(i).gt.klev-3) THEN
2662	IF (prt_level >= 9) THEN
2663	write(lunout,*)'La convection monte trop haut '
2664	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2665	ENDIF
2666	ENDIF
2667	ENDDO
2668	ELSE IF (iflag_con.eq.0) THEN
2669	write(lunout,*) 'On n appelle pas la convection'
2670	clwcon0=0.
2671	rnebcon0=0.
2672	d_t_con=0.
2673	d_q_con=0.
2674	d_u_con=0.
2675	d_v_con=0.
2676	rain_con=0.
2677	snow_con=0.
2678	bas=1
2679	top=1
2680	ELSE
2681	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2682	CALL abort_physic("physiq", "", 1)
2683	ENDIF
2684
2685	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2686	! . d_u_con, d_v_con)
2687
2688	!jyg Reinitialize proba_notrig and itapcv when convection has been called
2689	proba_notrig(:) = 1.
2690	itapcv = 0
2691	ENDIF ! (MOD(itapcv,cvpas).EQ.0)
2692	!
2693	itapcv = itapcv+1
2694
2695	!!!jyg Appel diagnostique a add_phys_tend pour tester la conservation de
2696	!!! l'energie dans les courants satures.
2697	!! d_t_con_sat(:,:) = d_t_con(:,:) - ftd(:,:)*dtime
2698	!! d_q_con_sat(:,:) = d_q_con(:,:) - fqd(:,:)*dtime
2699	!! dql_sat(:,:) = (wdtrainA(:,:)+wdtrainM(:,:))*dtime/zmasse(:,:)
2700	!! CALL add_phys_tend(d_u_con, d_v_con, d_t_con_sat, d_q_con_sat, dql_sat, &
2701	!! dqi0, paprs, 'convection_sat', abortphy, flag_inhib_tend,&
2702	!! itap, 1)
2703	!! call prt_enerbil('convection_sat',itap)
2704	!!
2705	!!
2706	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
2707	'convection',abortphy,flag_inhib_tend,itap,0)
2708	call prt_enerbil('convection',itap)
2709
2710	!-------------------------------------------------------------------------
2711
2712	IF (mydebug) THEN
2713	CALL writefield_phy('u_seri',u_seri,nbp_lev)
2714	CALL writefield_phy('v_seri',v_seri,nbp_lev)
2715	CALL writefield_phy('t_seri',t_seri,nbp_lev)
2716	CALL writefield_phy('q_seri',q_seri,nbp_lev)
2717	ENDIF
2718
2719	IF (check) THEN
2720	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2721	WRITE(lunout,*)"aprescon=", za
2722	zx_t = 0.0
2723	za = 0.0
2724	DO i = 1, klon
2725	za = za + cell_area(i)/REAL(klon)
2726	zx_t = zx_t + (rain_con(i)+ &
2727	snow_con(i))*cell_area(i)/REAL(klon)
2728	ENDDO
2729	zx_t = zx_t/za*phys_tstep
2730	WRITE(lunout,*)"Precip=", zx_t
2731	ENDIF
2732	IF (zx_ajustq) THEN
2733	DO i = 1, klon
2734	z_apres(i) = 0.0
2735	ENDDO
2736	DO k = 1, klev
2737	DO i = 1, klon
2738	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2739	*(paprs(i,k)-paprs(i,k+1))/RG
2740	ENDDO
2741	ENDDO
2742	DO i = 1, klon
2743	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*phys_tstep) &
2744	/z_apres(i)
2745	ENDDO
2746	DO k = 1, klev
2747	DO i = 1, klon
2748	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2749	z_factor(i).LT.(1.0-1.0E-08)) THEN
2750	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2751	ENDIF
2752	ENDDO
2753	ENDDO
2754	ENDIF
2755	zx_ajustq=.FALSE.
2756
2757	!
2758	!==========================================================================
2759	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
2760	!pour la couche limite diffuse pour l instant
2761	!
2762	!
2763	! nrlmd le 22/03/2011---Si on met les poches hors des thermiques
2764	! il faut rajouter cette tendance calcul\'ee hors des poches
2765	! froides
2766	!
2767	IF (iflag_wake>=1) THEN
2768	!
2769	!
2770	! Call wakes every "wkpas" step
2771	!
2772	IF (MOD(itapwk,wkpas).EQ.0) THEN
2773	!
2774	DO k=1,klev
2775	DO i=1,klon
2776	dt_dwn(i,k) = ftd(i,k)
2777	dq_dwn(i,k) = fqd(i,k)
2778	M_dwn(i,k) = dnwd0(i,k)
2779	M_up(i,k) = upwd(i,k)
2780	dt_a(i,k) = d_t_con(i,k)/phys_tstep - ftd(i,k)
2781	dq_a(i,k) = d_q_con(i,k)/phys_tstep - fqd(i,k)
2782	ENDDO
2783	ENDDO
2784
2785	IF (iflag_wake==2) THEN
2786	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2787	DO k = 1,klev
2788	dt_dwn(:,k)= dt_dwn(:,k)+ &
2789	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/phys_tstep
2790	dq_dwn(:,k)= dq_dwn(:,k)+ &
2791	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/phys_tstep
2792	ENDDO
2793	ELSEIF (iflag_wake==3) THEN
2794	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2795	DO k = 1,klev
2796	DO i=1,klon
2797	IF (rneb(i,k)==0.) THEN
2798	! On ne tient compte des tendances qu'en dehors des
2799	! nuages (c'est-\`a-dire a priri dans une region ou
2800	! l'eau se reevapore).
2801	dt_dwn(i,k)= dt_dwn(i,k)+ &
2802	ok_wk_lsp(i)*d_t_lsc(i,k)/phys_tstep
2803	dq_dwn(i,k)= dq_dwn(i,k)+ &
2804	ok_wk_lsp(i)*d_q_lsc(i,k)/phys_tstep
2805	ENDIF
2806	ENDDO
2807	ENDDO
2808	ENDIF
2809
2810	!
2811	!calcul caracteristiques de la poche froide
2812	CALL calWAKE (iflag_wake_tend, paprs, pplay, phys_tstep, &
2813	t_seri, q_seri, omega, &
2814	dt_dwn, dq_dwn, M_dwn, M_up, &
2815	dt_a, dq_a, &
2816	sigd, &
2817	wake_deltat, wake_deltaq, wake_s, wake_dens, &
2818	wake_dth, wake_h, &
2819	!! wake_pe, wake_fip, wake_gfl, &
2820	wake_pe, wake_fip_0, wake_gfl, & !! jyg
2821	d_t_wake, d_q_wake, &
2822	wake_k, t_x, q_x, &
2823	wake_omgbdth, wake_dp_omgb, &
2824	wake_dtKE, wake_dqKE, &
2825	wake_omg, wake_dp_deltomg, &
2826	wake_spread, wake_Cstar, d_deltat_wk_gw, &
2827	d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk)
2828	!
2829	!jyg Reinitialize itapwk when wakes have been called
2830	itapwk = 0
2831	ENDIF ! (MOD(itapwk,wkpas).EQ.0)
2832	!
2833	itapwk = itapwk+1
2834	!
2835	!-----------------------------------------------------------------------
2836	! ajout des tendances des poches froides
2837	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake', &
2838	abortphy,flag_inhib_tend,itap,0)
2839	call prt_enerbil('wake',itap)
2840	!------------------------------------------------------------------------
2841
2842	! Increment Wake state variables
2843	IF (iflag_wake_tend .GT. 0.) THEN
2844
2845	CALL add_wake_tend &
2846	(d_deltat_wk, d_deltaq_wk, d_s_wk, d_dens_wk, wake_k, &
2847	'wake', abortphy)
2848	call prt_enerbil('wake',itap)
2849	ENDIF ! (iflag_wake_tend .GT. 0.)
2850
2851	IF (iflag_alp_wk_cond .GT. 0.) THEN
2852
2853	CALL alpale_wk(phys_tstep, cell_area, wake_k, wake_s, wake_dens, wake_fip_0, &
2854	wake_fip)
2855	ELSE
2856	wake_fip(:) = wake_fip_0(:)
2857	ENDIF ! (iflag_alp_wk_cond .GT. 0.)
2858
2859	ENDIF ! (iflag_wake>=1)
2860	!
2861	!===================================================================
2862	! Convection seche (thermiques ou ajustement)
2863	!===================================================================
2864	!
2865	CALL stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
2866	,seuil_inversion,weak_inversion,dthmin)
2867
2868
2869
2870	d_t_ajsb(:,:)=0.
2871	d_q_ajsb(:,:)=0.
2872	d_t_ajs(:,:)=0.
2873	d_u_ajs(:,:)=0.
2874	d_v_ajs(:,:)=0.
2875	d_q_ajs(:,:)=0.
2876	clwcon0th(:,:)=0.
2877	!
2878	! fm_therm(:,:)=0.
2879	! entr_therm(:,:)=0.
2880	! detr_therm(:,:)=0.
2881	!
2882	IF (prt_level>9) WRITE(lunout,*) &
2883	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
2884	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2885	IF (iflag_thermals<0) THEN
2886	! Rien
2887	! ====
2888	IF (prt_level>9) WRITE(lunout,*)'pas de convection seche'
2889
2890
2891	ELSE
2892
2893	! Thermiques
2894	! ==========
2895	IF (prt_level>9) WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
2896	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2897
2898
2899	!cc nrlmd le 10/04/2012
2900	DO k=1,klev+1
2901	DO i=1,klon
2902	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2903	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2904	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2905	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2906	ENDDO
2907	ENDDO
2908	!cc fin nrlmd le 10/04/2012
2909
2910	IF (iflag_thermals>=1) THEN
2911	!jyg<
2912	!! IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2913	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
2914	! Appel des thermiques avec les profils exterieurs aux poches
2915	DO k=1,klev
2916	DO i=1,klon
2917	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2918	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2919	u_therm(i,k) = u_seri(i,k)
2920	v_therm(i,k) = v_seri(i,k)
2921	ENDDO
2922	ENDDO
2923	ELSE
2924	! Appel des thermiques avec les profils moyens
2925	DO k=1,klev
2926	DO i=1,klon
2927	t_therm(i,k) = t_seri(i,k)
2928	q_therm(i,k) = q_seri(i,k)
2929	u_therm(i,k) = u_seri(i,k)
2930	v_therm(i,k) = v_seri(i,k)
2931	ENDDO
2932	ENDDO
2933	ENDIF
2934	!>jyg
2935	CALL calltherm(pdtphys &
2936	,pplay,paprs,pphi,weak_inversion &
2937	! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
2938	,u_therm,v_therm,t_therm,q_therm,zqsat,debut & !jyg
2939	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
2940	,fm_therm,entr_therm,detr_therm &
2941	,zqasc,clwcon0th,lmax_th,ratqscth &
2942	,ratqsdiff,zqsatth &
2943	!on rajoute ale et alp, et les
2944	!caracteristiques de la couche alim
2945	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
2946	,ztv,zpspsk,ztla,zthl &
2947	!cc nrlmd le 10/04/2012
2948	,pbl_tke_input,pctsrf,omega,cell_area &
2949	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
2950	,n2,s2,ale_bl_stat &
2951	,therm_tke_max,env_tke_max &
2952	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
2953	,alp_bl_conv,alp_bl_stat &
2954	!cc fin nrlmd le 10/04/2012
2955	,zqla,ztva )
2956	!
2957	!jyg<
2958	!!jyg IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2959	IF (mod(iflag_pbl_split/10,10) .GE. 1) THEN
2960	! Si les thermiques ne sont presents que hors des
2961	! poches, la tendance moyenne associ\'ee doit etre
2962	! multipliee par la fraction surfacique qu'ils couvrent.
2963	DO k=1,klev
2964	DO i=1,klon
2965	!
2966	d_deltat_the(i,k) = - d_t_ajs(i,k)
2967	d_deltaq_the(i,k) = - d_q_ajs(i,k)
2968	!
2969	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
2970	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
2971	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
2972	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
2973	!
2974	ENDDO
2975	ENDDO
2976	!
2977	CALL add_wake_tend &
2978	(d_deltat_the, d_deltaq_the, dsig0, ddens0, wkoccur1, 'the', abortphy)
2979	call prt_enerbil('the',itap)
2980	!
2981	ENDIF ! (mod(iflag_pbl_split/10,10) .GE. 1)
2982	!
2983	CALL add_phys_tend(d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs, &
2984	dql0,dqi0,paprs,'thermals', abortphy,flag_inhib_tend,itap,0)
2985	call prt_enerbil('thermals',itap)
2986	!
2987	!
2988	CALL alpale_th( phys_tstep, lmax_th, t_seri, cell_area, &
2989	cin, s2, n2, &
2990	ale_bl_trig, ale_bl_stat, ale_bl, &
2991	alp_bl, alp_bl_stat, &
2992	proba_notrig, random_notrig)
2993	!>jyg
2994
2995	! ------------------------------------------------------------------
2996	! Transport de la TKE par les panaches thermiques.
2997	! FH : 2010/02/01
2998	! if (iflag_pbl.eq.10) then
2999	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
3000	! s rg,paprs,pbl_tke)
3001	! endif
3002	! -------------------------------------------------------------------
3003
3004	DO i=1,klon
3005	! zmax_th(i)=pphi(i,lmax_th(i))/rg
3006	!CR:04/05/12:correction calcul zmax
3007	zmax_th(i)=zmax0(i)
3008	ENDDO
3009
3010	ENDIF
3011
3012	! Ajustement sec
3013	! ==============
3014
3015	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
3016	! a partir du sommet des thermiques.
3017	! Dans le cas contraire, on demarre au niveau 1.
3018
3019	IF (iflag_thermals>=13.or.iflag_thermals<=0) THEN
3020
3021	IF (iflag_thermals.eq.0) THEN
3022	IF (prt_level>9) WRITE(lunout,*)'ajsec'
3023	limbas(:)=1
3024	ELSE
3025	limbas(:)=lmax_th(:)
3026	ENDIF
3027
3028	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
3029	! pour des test de convergence numerique.
3030	! Le nouveau ajsec est a priori mieux, meme pour le cas
3031	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
3032	! non nulles numeriquement pour des mailles non concernees.
3033
3034	IF (iflag_thermals==0) THEN
3035	! Calling adjustment alone (but not the thermal plume model)
3036	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
3037	, d_t_ajsb, d_q_ajsb)
3038	ELSE IF (iflag_thermals>0) THEN
3039	! Calling adjustment above the top of thermal plumes
3040	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
3041	, d_t_ajsb, d_q_ajsb)
3042	ENDIF
3043
3044	!--------------------------------------------------------------------
3045	! ajout des tendances de l'ajustement sec ou des thermiques
3046	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs, &
3047	'ajsb',abortphy,flag_inhib_tend,itap,0)
3048	call prt_enerbil('ajsb',itap)
3049	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
3050	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
3051
3052	!---------------------------------------------------------------------
3053
3054	ENDIF
3055
3056	ENDIF
3057	!
3058	!===================================================================
3059	! Computation of ratqs, the width (normalized) of the subrid scale
3060	! water distribution
3061	CALL calcratqs(klon,klev,prt_level,lunout, &
3062	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3063	ratqsbas,ratqshaut,ratqsp0, ratqsdp, &
3064	tau_ratqs,fact_cldcon, &
3065	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3066	paprs,pplay,q_seri,zqsat,fm_therm, &
3067	ratqs,ratqsc)
3068
3069
3070	!
3071	! Appeler le processus de condensation a grande echelle
3072	! et le processus de precipitation
3073	!-------------------------------------------------------------------------
3074	IF (prt_level .GE.10) THEN
3075	print *,'itap, ->fisrtilp ',itap
3076	ENDIF
3077	!
3078	CALL fisrtilp(phys_tstep,paprs,pplay, &
3079	t_seri, q_seri,ptconv,ratqs, &
3080	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3081	rain_lsc, snow_lsc, &
3082	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3083	frac_impa, frac_nucl, beta_prec_fisrt, &
3084	prfl, psfl, rhcl, &
3085	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3086	iflag_ice_thermo)
3087	!
3088	WHERE (rain_lsc < 0) rain_lsc = 0.
3089	WHERE (snow_lsc < 0) snow_lsc = 0.
3090
3091	!+JLD
3092	! write(*,9000) 'phys lsc',"enerbil: bil_q, bil_e,",rain_lsc+snow_lsc &
3093	! & ,((rcw-rcpd)rain_lsc + (rcs-rcpd)snow_lsc)t_seri(1,1)-rlvttrain_lsc+rlstt*snow_lsc &
3094	! & ,rain_lsc,snow_lsc
3095	! write(*,9000) "rcpv","rcw",rcpv,rcw,rcs,t_seri(1,1)
3096	!-JLD
3097	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs, &
3098	'lsc',abortphy,flag_inhib_tend,itap,0)
3099	call prt_enerbil('lsc',itap)
3100	rain_num(:)=0.
3101	DO k = 1, klev
3102	DO i = 1, klon
3103	IF (ql_seri(i,k)>oliqmax) THEN
3104	rain_num(i)=rain_num(i)+(ql_seri(i,k)-oliqmax)*zmasse(i,k)/pdtphys
3105	ql_seri(i,k)=oliqmax
3106	ENDIF
3107	ENDDO
3108	ENDDO
3109	IF (nqo==3) THEN
3110	DO k = 1, klev
3111	DO i = 1, klon
3112	IF (qs_seri(i,k)>oicemax) THEN
3113	rain_num(i)=rain_num(i)+(qs_seri(i,k)-oicemax)*zmasse(i,k)/pdtphys
3114	qs_seri(i,k)=oicemax
3115	ENDIF
3116	ENDDO
3117	ENDDO
3118	ENDIF
3119
3120	!---------------------------------------------------------------------------
3121	DO k = 1, klev
3122	DO i = 1, klon
3123	cldfra(i,k) = rneb(i,k)
3124	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3125	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3126	ENDDO
3127	ENDDO
3128	IF (check) THEN
3129	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3130	WRITE(lunout,*)"apresilp=", za
3131	zx_t = 0.0
3132	za = 0.0
3133	DO i = 1, klon
3134	za = za + cell_area(i)/REAL(klon)
3135	zx_t = zx_t + (rain_lsc(i) &
3136	+ snow_lsc(i))*cell_area(i)/REAL(klon)
3137	ENDDO
3138	zx_t = zx_t/za*phys_tstep
3139	WRITE(lunout,*)"Precip=", zx_t
3140	ENDIF
3141
3142	IF (mydebug) THEN
3143	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3144	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3145	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3146	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3147	ENDIF
3148
3149	!
3150	!-------------------------------------------------------------------
3151	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3152	!-------------------------------------------------------------------
3153
3154	! 1. NUAGES CONVECTIFS
3155	!
3156	!IM cf FH
3157	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3158	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3159	snow_tiedtke=0.
3160	! print*,'avant calcul de la pseudo precip '
3161	! print*,'iflag_cld_th',iflag_cld_th
3162	IF (iflag_cld_th.eq.-1) THEN
3163	rain_tiedtke=rain_con
3164	ELSE
3165	! print*,'calcul de la pseudo precip '
3166	rain_tiedtke=0.
3167	! print*,'calcul de la pseudo precip 0'
3168	DO k=1,klev
3169	DO i=1,klon
3170	IF (d_q_con(i,k).lt.0.) THEN
3171	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3172	*(paprs(i,k)-paprs(i,k+1))/rg
3173	ENDIF
3174	ENDDO
3175	ENDDO
3176	ENDIF
3177	!
3178	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3179	!
3180
3181	! Nuages diagnostiques pour Tiedtke
3182	CALL diagcld1(paprs,pplay, &
3183	!IM cf FH. rain_con,snow_con,ibas_con,itop_con,
3184	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3185	diafra,dialiq)
3186	DO k = 1, klev
3187	DO i = 1, klon
3188	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3189	cldliq(i,k) = dialiq(i,k)
3190	cldfra(i,k) = diafra(i,k)
3191	ENDIF
3192	ENDDO
3193	ENDDO
3194
3195	ELSE IF (iflag_cld_th.ge.3) THEN
3196	! On prend pour les nuages convectifs le max du calcul de la
3197	! convection et du calcul du pas de temps precedent diminue d'un facteur
3198	! facttemps
3199	facteur = pdtphys *facttemps
3200	DO k=1,klev
3201	DO i=1,klon
3202	rnebcon(i,k)=rnebcon(i,k)*facteur
3203	IF (rnebcon0(i,k)clwcon0(i,k).GT.rnebcon(i,k)clwcon(i,k)) THEN
3204	rnebcon(i,k)=rnebcon0(i,k)
3205	clwcon(i,k)=clwcon0(i,k)
3206	ENDIF
3207	ENDDO
3208	ENDDO
3209
3210	! On prend la somme des fractions nuageuses et des contenus en eau
3211
3212	IF (iflag_cld_th>=5) THEN
3213
3214	DO k=1,klev
3215	ptconvth(:,k)=fm_therm(:,k+1)>0.
3216	ENDDO
3217
3218	IF (iflag_coupl==4) THEN
3219
3220	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3221	! convectives et lsc dans la partie des thermiques
3222	! Le controle par iflag_coupl est peut etre provisoire.
3223	DO k=1,klev
3224	DO i=1,klon
3225	IF (ptconv(i,k).AND.ptconvth(i,k)) THEN
3226	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3227	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3228	ELSE IF (ptconv(i,k)) THEN
3229	cldfra(i,k)=rnebcon(i,k)
3230	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3231	ENDIF
3232	ENDDO
3233	ENDDO
3234
3235	ELSE IF (iflag_coupl==5) THEN
3236	DO k=1,klev
3237	DO i=1,klon
3238	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3239	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3240	ENDDO
3241	ENDDO
3242
3243	ELSE
3244
3245	! Si on est sur un point touche par la convection
3246	! profonde et pas par les thermiques, on prend la
3247	! couverture nuageuse et l'eau nuageuse de la convection
3248	! profonde.
3249
3250	!IM/FH: 2011/02/23
3251	! definition des points sur lesquels ls thermiques sont actifs
3252
3253	DO k=1,klev
3254	DO i=1,klon
3255	IF (ptconv(i,k).AND. .NOT.ptconvth(i,k)) THEN
3256	cldfra(i,k)=rnebcon(i,k)
3257	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3258	ENDIF
3259	ENDDO
3260	ENDDO
3261
3262	ENDIF
3263
3264	ELSE
3265
3266	! Ancienne version
3267	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3268	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3269	ENDIF
3270
3271	ENDIF
3272
3273	! plulsc(:)=0.
3274	! do k=1,klev,-1
3275	! do i=1,klon
3276	! zzz=prfl(:,k)+psfl(:,k)
3277	! if (.not.ptconvth.zzz.gt.0.)
3278	! enddo prfl, psfl,
3279	! enddo
3280	!
3281	! 2. NUAGES STARTIFORMES
3282	!
3283	IF (ok_stratus) THEN
3284	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3285	DO k = 1, klev
3286	DO i = 1, klon
3287	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3288	cldliq(i,k) = dialiq(i,k)
3289	cldfra(i,k) = diafra(i,k)
3290	ENDIF
3291	ENDDO
3292	ENDDO
3293	ENDIF
3294	!
3295	! Precipitation totale
3296	!
3297	DO i = 1, klon
3298	rain_fall(i) = rain_con(i) + rain_lsc(i)
3299	snow_fall(i) = snow_con(i) + snow_lsc(i)
3300	ENDDO
3301	!
3302	! Calculer l'humidite relative pour diagnostique
3303	!
3304	DO k = 1, klev
3305	DO i = 1, klon
3306	zx_t = t_seri(i,k)
3307	IF (thermcep) THEN
3308	!! if (iflag_ice_thermo.eq.0) then !jyg
3309	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3310	!! else !jyg
3311	!! zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t)) !jyg
3312	!! endif !jyg
3313	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3314	zx_qs = MIN(0.5,zx_qs)
3315	zcor = 1./(1.-retv*zx_qs)
3316	zx_qs = zx_qs*zcor
3317	ELSE
3318	!! IF (zx_t.LT.t_coup) THEN !jyg
3319	IF (zx_t.LT.rtt) THEN !jyg
3320	zx_qs = qsats(zx_t)/pplay(i,k)
3321	ELSE
3322	zx_qs = qsatl(zx_t)/pplay(i,k)
3323	ENDIF
3324	ENDIF
3325	zx_rh(i,k) = q_seri(i,k)/zx_qs
3326	zqsat(i,k)=zx_qs
3327	ENDDO
3328	ENDDO
3329
3330	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3331	! equivalente a 2m (tpote) pour diagnostique
3332	!
3333	DO i = 1, klon
3334	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3335	IF (thermcep) THEN
3336	IF(zt2m(i).LT.RTT) then
3337	Lheat=RLSTT
3338	ELSE
3339	Lheat=RLVTT
3340	ENDIF
3341	ELSE
3342	IF (zt2m(i).LT.RTT) THEN
3343	Lheat=RLSTT
3344	ELSE
3345	Lheat=RLVTT
3346	ENDIF
3347	ENDIF
3348	tpote(i) = tpot(i)* &
3349	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3350	ENDDO
3351
3352	IF (type_trac == 'inca') THEN
3353	#ifdef INCA
3354	CALL VTe(VTphysiq)
3355	CALL VTb(VTinca)
3356	calday = REAL(days_elapsed + 1) + jH_cur
3357
3358	CALL chemtime(itap+itau_phy-1, date0, phys_tstep, itap)
3359	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
3360	CALL AEROSOL_METEO_CALC( &
3361	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3362	prfl,psfl,pctsrf,cell_area, &
3363	latitude_deg,longitude_deg,u10m,v10m)
3364	ENDIF
3365
3366	zxsnow_dummy(:) = 0.0
3367
3368	CALL chemhook_begin (calday, &
3369	days_elapsed+1, &
3370	jH_cur, &
3371	pctsrf(1,1), &
3372	latitude_deg, &
3373	longitude_deg, &
3374	cell_area, &
3375	paprs, &
3376	pplay, &
3377	coefh(1:klon,1:klev,is_ave), &
3378	pphi, &
3379	t_seri, &
3380	u, &
3381	v, &
3382	wo(:, :, 1), &
3383	q_seri, &
3384	zxtsol, &
3385	zxsnow_dummy, &
3386	solsw, &
3387	albsol1, &
3388	rain_fall, &
3389	snow_fall, &
3390	itop_con, &
3391	ibas_con, &
3392	cldfra, &
3393	nbp_lon, &
3394	nbp_lat-1, &
3395	tr_seri, &
3396	ftsol, &
3397	paprs, &
3398	cdragh, &
3399	cdragm, &
3400	pctsrf, &
3401	pdtphys, &
3402	itap)
3403
3404	CALL VTe(VTinca)
3405	CALL VTb(VTphysiq)
3406	#endif
3407	ENDIF !type_trac = inca
3408
3409
3410	!
3411	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3412	!
3413	IF (MOD(itaprad,radpas).EQ.0) THEN
3414
3415	!
3416	!jq - introduce the aerosol direct and first indirect radiative forcings
3417	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3418	IF (flag_aerosol .GT. 0) THEN
3419	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3420	IF (.NOT. aerosol_couple) THEN
3421	!
3422	CALL readaerosol_optic( &
3423	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3424	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3425	mass_solu_aero, mass_solu_aero_pi, &
3426	tau_aero, piz_aero, cg_aero, &
3427	tausum_aero, tau3d_aero)
3428	ENDIF
3429	ELSE ! RRTM radiation
3430	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3431	abort_message='config_inca=aero et rrtm=1 impossible'
3432	CALL abort_physic(modname,abort_message,1)
3433	ELSE
3434	!
3435	#ifdef CPP_RRTM
3436	IF (NSW.EQ.6) THEN
3437	!--new aerosol properties SW and LW
3438	!
3439	#ifdef CPP_Dust
3440	!--SPL aerosol model
3441	CALL splaerosol_optic_rrtm( ok_alw, pplay, paprs, t_seri, rhcl, &
3442	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3443	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3444	tausum_aero, tau3d_aero)
3445	#else
3446	!--climatologies or INCA aerosols
3447	CALL readaerosol_optic_rrtm( debut, aerosol_couple, ok_alw, &
3448	new_aod, flag_aerosol, flag_bc_internal_mixture, itap, jD_cur-jD_ref, &
3449	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3450	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3451	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3452	tausum_aero, drytausum_aero, tau3d_aero)
3453	#endif
3454	!
3455	ELSE IF (NSW.EQ.2) THEN
3456	!--for now we use the old aerosol properties
3457	!
3458	CALL readaerosol_optic( &
3459	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3460	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3461	mass_solu_aero, mass_solu_aero_pi, &
3462	tau_aero, piz_aero, cg_aero, &
3463	tausum_aero, tau3d_aero)
3464	!
3465	!--natural aerosols
3466	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
3467	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
3468	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
3469	!--all aerosols
3470	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
3471	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
3472	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
3473	!
3474	!--no LW optics
3475	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3476	!
3477	ELSE
3478	abort_message='Only NSW=2 or 6 are possible with ' &
3479	// 'aerosols and iflag_rrtm=1'
3480	CALL abort_physic(modname,abort_message,1)
3481	ENDIF
3482	#else
3483	abort_message='You should compile with -rrtm if running ' &
3484	// 'with iflag_rrtm=1'
3485	CALL abort_physic(modname,abort_message,1)
3486	#endif
3487	!
3488	ENDIF
3489	ENDIF
3490	ELSE !--flag_aerosol = 0
3491	tausum_aero(:,:,:) = 0.
3492	drytausum_aero(:,:) = 0.
3493	mass_solu_aero(:,:) = 0.
3494	mass_solu_aero_pi(:,:) = 0.
3495	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3496	tau_aero(:,:,:,:) = 1.e-15
3497	piz_aero(:,:,:,:) = 1.
3498	cg_aero(:,:,:,:) = 0.
3499	ELSE
3500	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
3501	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3502	piz_aero_sw_rrtm(:,:,:,:) = 1.0
3503	cg_aero_sw_rrtm(:,:,:,:) = 0.0
3504	ENDIF
3505	ENDIF
3506	!
3507	!--WMO criterion to determine tropopause
3508	CALL stratosphere_mask(missing_val, t_seri, pplay, latitude_deg)
3509	!
3510	!--STRAT AEROSOL
3511	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
3512	IF (flag_aerosol_strat.GT.0) THEN
3513	IF (prt_level .GE.10) THEN
3514	PRINT *,'appel a readaerosolstrat', mth_cur
3515	ENDIF
3516	IF (iflag_rrtm.EQ.0) THEN
3517	IF (flag_aerosol_strat.EQ.1) THEN
3518	CALL readaerosolstrato(debut)
3519	ELSE
3520	abort_message='flag_aerosol_strat must equal 1 for rrtm=0'
3521	CALL abort_physic(modname,abort_message,1)
3522	ENDIF
3523	ELSE
3524	#ifdef CPP_RRTM
3525	#ifndef CPP_StratAer
3526	!--prescribed strat aerosols
3527	!--only in the case of non-interactive strat aerosols
3528	IF (flag_aerosol_strat.EQ.1) THEN
3529	CALL readaerosolstrato1_rrtm(debut)
3530	ELSEIF (flag_aerosol_strat.EQ.2) THEN
3531	CALL readaerosolstrato2_rrtm(debut)
3532	ELSE
3533	abort_message='flag_aerosol_strat must equal 1 or 2 for rrtm=1'
3534	CALL abort_physic(modname,abort_message,1)
3535	ENDIF
3536	#endif
3537	#else
3538	abort_message='You should compile with -rrtm if running ' &
3539	// 'with iflag_rrtm=1'
3540	CALL abort_physic(modname,abort_message,1)
3541	#endif
3542	ENDIF
3543	ENDIF
3544	!
3545	#ifdef CPP_RRTM
3546	#ifdef CPP_StratAer
3547	!--compute stratospheric mask
3548	CALL stratosphere_mask(missing_val, t_seri, pplay, latitude_deg)
3549	!--interactive strat aerosols
3550	CALL calcaerosolstrato_rrtm(pplay,t_seri,paprs,debut)
3551	#endif
3552	#endif
3553	!--fin STRAT AEROSOL
3554	!
3555
3556	! Calculer les parametres optiques des nuages et quelques
3557	! parametres pour diagnostiques:
3558	!
3559	IF (aerosol_couple.AND.config_inca=='aero') THEN
3560	mass_solu_aero(:,:) = ccm(:,:,1)
3561	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3562	ENDIF
3563
3564	IF (ok_newmicro) then
3565	IF (iflag_rrtm.NE.0) THEN
3566	#ifdef CPP_RRTM
3567	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3568	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 ' &
3569	// 'pour ok_cdnc'
3570	CALL abort_physic(modname,abort_message,1)
3571	ENDIF
3572	#else
3573
3574	abort_message='You should compile with -rrtm if running with '//'iflag_rrtm=1'
3575	CALL abort_physic(modname,abort_message,1)
3576	#endif
3577	ENDIF
3578	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3579	paprs, pplay, t_seri, cldliq, cldfra, &
3580	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3581	flwp, fiwp, flwc, fiwc, &
3582	mass_solu_aero, mass_solu_aero_pi, &
3583	cldtaupi, re, fl, ref_liq, ref_ice, &
3584	ref_liq_pi, ref_ice_pi)
3585	ELSE
3586	CALL nuage (paprs, pplay, &
3587	t_seri, cldliq, cldfra, cldtau, cldemi, &
3588	cldh, cldl, cldm, cldt, cldq, &
3589	ok_aie, &
3590	mass_solu_aero, mass_solu_aero_pi, &
3591	bl95_b0, bl95_b1, &
3592	cldtaupi, re, fl)
3593	ENDIF
3594	!
3595	!IM betaCRF
3596	!
3597	cldtaurad = cldtau
3598	cldtaupirad = cldtaupi
3599	cldemirad = cldemi
3600	cldfrarad = cldfra
3601
3602	!
3603	IF (lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3604	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3605	!
3606	! global
3607	!
3608	DO k=1, klev
3609	DO i=1, klon
3610	IF (pplay(i,k).GE.pfree) THEN
3611	beta(i,k) = beta_pbl
3612	ELSE
3613	beta(i,k) = beta_free
3614	ENDIF
3615	IF (mskocean_beta) THEN
3616	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3617	ENDIF
3618	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3619	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3620	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3621	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3622	ENDDO
3623	ENDDO
3624	!
3625	ELSE
3626	!
3627	! regional
3628	!
3629	DO k=1, klev
3630	DO i=1,klon
3631	!
3632	IF (longitude_deg(i).ge.lon1_beta.AND. &
3633	longitude_deg(i).le.lon2_beta.AND. &
3634	latitude_deg(i).le.lat1_beta.AND. &
3635	latitude_deg(i).ge.lat2_beta) THEN
3636	IF (pplay(i,k).GE.pfree) THEN
3637	beta(i,k) = beta_pbl
3638	ELSE
3639	beta(i,k) = beta_free
3640	ENDIF
3641	IF (mskocean_beta) THEN
3642	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3643	ENDIF
3644	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3645	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3646	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3647	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3648	ENDIF
3649	!
3650	ENDDO
3651	ENDDO
3652	!
3653	ENDIF
3654
3655	!lecture de la chlorophylle pour le nouvel albedo de Sunghye Baek
3656	IF (ok_chlorophyll) THEN
3657	print*,"-- reading chlorophyll"
3658	CALL readchlorophyll(debut)
3659	ENDIF
3660
3661	!--if ok_suntime_rrtm we use ancillay data for RSUN
3662	!--previous values are therefore overwritten
3663	!--this is needed for CMIP6 runs
3664	!--and only possible for new radiation scheme
3665	IF (iflag_rrtm.EQ.1.AND.ok_suntime_rrtm) THEN
3666	#ifdef CPP_RRTM
3667	CALL read_rsun_rrtm(debut)
3668	#endif
3669	ENDIF
3670
3671	IF (mydebug) THEN
3672	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3673	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3674	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3675	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3676	ENDIF
3677
3678	!
3679	!sonia : If Iflag_radia >=2, pertubation of some variables
3680	!input to radiation (DICE)
3681	!
3682	IF (iflag_radia .ge. 2) THEN
3683	zsav_tsol (:) = zxtsol(:)
3684	CALL perturb_radlwsw(zxtsol,iflag_radia)
3685	ENDIF
3686
3687	IF (aerosol_couple.AND.config_inca=='aero') THEN
3688	#ifdef INCA
3689	CALL radlwsw_inca &
3690	(kdlon,kflev,dist, rmu0, fract, solaire, &
3691	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3692	size(wo,3), wo, &
3693	cldfrarad, cldemirad, cldtaurad, &
3694	heat,heat0,cool,cool0,albpla, &
3695	topsw,toplw,solsw,sollw, &
3696	sollwdown, &
3697	topsw0,toplw0,solsw0,sollw0, &
3698	lwdn0, lwdn, lwup0, lwup, &
3699	swdn0, swdn, swup0, swup, &
3700	ok_ade, ok_aie, &
3701	tau_aero, piz_aero, cg_aero, &
3702	topswad_aero, solswad_aero, &
3703	topswad0_aero, solswad0_aero, &
3704	topsw_aero, topsw0_aero, &
3705	solsw_aero, solsw0_aero, &
3706	cldtaupirad, &
3707	topswai_aero, solswai_aero)
3708	#endif
3709	ELSE
3710	!
3711	!IM calcul radiatif pour le cas actuel
3712	!
3713	RCO2 = RCO2_act
3714	RCH4 = RCH4_act
3715	RN2O = RN2O_act
3716	RCFC11 = RCFC11_act
3717	RCFC12 = RCFC12_act
3718	!
3719	IF (prt_level .GE.10) THEN
3720	print *,' ->radlwsw, number 1 '
3721	ENDIF
3722	!
3723	CALL radlwsw &
3724	(dist, rmu0, fract, &
3725	!albedo SB >>>
3726	! paprs, pplay,zxtsol,albsol1, albsol2, &
3727	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3728	!albedo SB <<<
3729	t_seri,q_seri,wo, &
3730	cldfrarad, cldemirad, cldtaurad, &
3731	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
3732	flag_aerosol_strat, &
3733	tau_aero, piz_aero, cg_aero, &
3734	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3735	! Rajoute par OB pour RRTM
3736	tau_aero_lw_rrtm, &
3737	cldtaupirad,new_aod, &
3738	zqsat, flwc, fiwc, &
3739	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3740	heat,heat0,cool,cool0,albpla, &
3741	topsw,toplw,solsw,sollw, &
3742	sollwdown, &
3743	topsw0,toplw0,solsw0,sollw0, &
3744	lwdn0, lwdn, lwup0, lwup, &
3745	swdn0, swdn, swup0, swup, &
3746	topswad_aero, solswad_aero, &
3747	topswai_aero, solswai_aero, &
3748	topswad0_aero, solswad0_aero, &
3749	topsw_aero, topsw0_aero, &
3750	solsw_aero, solsw0_aero, &
3751	topswcf_aero, solswcf_aero, &
3752	!-C. Kleinschmitt for LW diagnostics
3753	toplwad_aero, sollwad_aero,&
3754	toplwai_aero, sollwai_aero, &
3755	toplwad0_aero, sollwad0_aero,&
3756	!-end
3757	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3758	ZSWFT0_i, ZFSDN0, ZFSUP0)
3759
3760	#ifndef CPP_XIOS
3761	!--OB 30/05/2016 modified 21/10/2016
3762	!--here we return swaero_diag and dryaod_diag to FALSE
3763	!--and histdef will switch it back to TRUE if necessary
3764	!--this is necessary to get the right swaero at first step
3765	!--but only in the case of no XIOS as XIOS is covered elsewhere
3766	IF (debut) swaero_diag = .FALSE.
3767	IF (debut) dryaod_diag = .FALSE.
3768	!--IM 15/09/2017 here we return ok_4xCO2atm to FALSE
3769	!--as for swaero_diag, see above
3770	IF (debut) ok_4xCO2atm = .FALSE.
3771	#endif
3772	!
3773	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3774	!IM des taux doit etre different du taux actuel
3775	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3776	!
3777	IF (RCO2_per.NE.RCO2_act.OR. &
3778	RCH4_per.NE.RCH4_act.OR. &
3779	RN2O_per.NE.RN2O_act.OR. &
3780	RCFC11_per.NE.RCFC11_act.OR. &
3781	RCFC12_per.NE.RCFC12_act) ok_4xCO2atm =.TRUE.
3782	!
3783	IF (ok_4xCO2atm) THEN
3784	!
3785	RCO2 = RCO2_per
3786	RCH4 = RCH4_per
3787	RN2O = RN2O_per
3788	RCFC11 = RCFC11_per
3789	RCFC12 = RCFC12_per
3790	!
3791	IF (prt_level .GE.10) THEN
3792	print *,' ->radlwsw, number 2 '
3793	ENDIF
3794	!
3795	CALL radlwsw &
3796	(dist, rmu0, fract, &
3797	!albedo SB >>>
3798	! paprs, pplay,zxtsol,albsol1, albsol2, &
3799	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3800	!albedo SB <<<
3801	t_seri,q_seri,wo, &
3802	cldfrarad, cldemirad, cldtaurad, &
3803	ok_ade.OR.flag_aerosol_strat.GT.0, ok_aie, flag_aerosol, &
3804	flag_aerosol_strat, &
3805	tau_aero, piz_aero, cg_aero, &
3806	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3807	! Rajoute par OB pour RRTM
3808	tau_aero_lw_rrtm, &
3809	cldtaupi,new_aod, &
3810	zqsat, flwc, fiwc, &
3811	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3812	heatp,heat0p,coolp,cool0p,albplap, &
3813	topswp,toplwp,solswp,sollwp, &
3814	sollwdownp, &
3815	topsw0p,toplw0p,solsw0p,sollw0p, &
3816	lwdn0p, lwdnp, lwup0p, lwupp, &
3817	swdn0p, swdnp, swup0p, swupp, &
3818	topswad_aerop, solswad_aerop, &
3819	topswai_aerop, solswai_aerop, &
3820	topswad0_aerop, solswad0_aerop, &
3821	topsw_aerop, topsw0_aerop, &
3822	solsw_aerop, solsw0_aerop, &
3823	topswcf_aerop, solswcf_aerop, &
3824	!-C. Kleinschmitt for LW diagnostics
3825	toplwad_aerop, sollwad_aerop,&
3826	toplwai_aerop, sollwai_aerop, &
3827	toplwad0_aerop, sollwad0_aerop,&
3828	!-end
3829	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3830	ZSWFT0_i, ZFSDN0, ZFSUP0)
3831	endif
3832	!
3833	ENDIF ! aerosol_couple
3834	itaprad = 0
3835	!
3836	! If Iflag_radia >=2, reset pertubed variables
3837	!
3838	IF (iflag_radia .ge. 2) THEN
3839	zxtsol(:) = zsav_tsol (:)
3840	ENDIF
3841	ENDIF ! MOD(itaprad,radpas)
3842	itaprad = itaprad + 1
3843
3844	IF (iflag_radia.eq.0) THEN
3845	IF (prt_level.ge.9) THEN
3846	PRINT *,'--------------------------------------------------'
3847	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3848	PRINT *,'>>>> heat et cool mis a zero '
3849	PRINT *,'--------------------------------------------------'
3850	ENDIF
3851	heat=0.
3852	cool=0.
3853	sollw=0. ! MPL 01032011
3854	solsw=0.
3855	radsol=0.
3856	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3857	swup0=0.
3858	lwup=0.
3859	lwup0=0.
3860	lwdn=0.
3861	lwdn0=0.
3862	ENDIF
3863
3864	!
3865	! Calculer radsol a l'exterieur de radlwsw
3866	! pour prendre en compte le cycle diurne
3867	! recode par Olivier Boucher en sept 2015
3868	!
3869	radsol=solsw*swradcorr+sollw
3870
3871	IF (ok_4xCO2atm) THEN
3872	radsolp=solswp*swradcorr+sollwp
3873	ENDIF
3874
3875	!
3876	! Ajouter la tendance des rayonnements (tous les pas)
3877	! avec une correction pour le cycle diurne dans le SW
3878	!
3879
3880	DO k=1, klev
3881	d_t_swr(:,k)=swradcorr(:)heat(:,k)phys_tstep/RDAY
3882	d_t_sw0(:,k)=swradcorr(:)heat0(:,k)phys_tstep/RDAY
3883	d_t_lwr(:,k)=-cool(:,k)*phys_tstep/RDAY
3884	d_t_lw0(:,k)=-cool0(:,k)*phys_tstep/RDAY
3885	ENDDO
3886
3887	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy,flag_inhib_tend,itap,0)
3888	call prt_enerbil('SW',itap)
3889	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy,flag_inhib_tend,itap,0)
3890	call prt_enerbil('LW',itap)
3891
3892	!
3893	IF (mydebug) THEN
3894	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3895	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3896	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3897	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3898	ENDIF
3899
3900	! Calculer l'hydrologie de la surface
3901	!
3902	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3903	! . agesno, ftsol,fqsurf,fsnow, ruis)
3904	!
3905
3906	!
3907	! Calculer le bilan du sol et la derive de temperature (couplage)
3908	!
3909	DO i = 1, klon
3910	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3911	! a la demande de JLD
3912	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3913	ENDDO
3914	!
3915	!moddeblott(jan95)
3916	! Appeler le programme de parametrisation de l'orographie
3917	! a l'echelle sous-maille:
3918	!
3919	IF (prt_level .GE.10) THEN
3920	print *,' call orography ? ', ok_orodr
3921	ENDIF
3922	!
3923	IF (ok_orodr) THEN
3924	!
3925	! selection des points pour lesquels le shema est actif:
3926	igwd=0
3927	DO i=1,klon
3928	itest(i)=0
3929	! IF ((zstd(i).gt.10.0)) THEN
3930	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3931	itest(i)=1
3932	igwd=igwd+1
3933	idx(igwd)=i
3934	ENDIF
3935	ENDDO
3936	! igwdim=MAX(1,igwd)
3937	!
3938	IF (ok_strato) THEN
3939
3940	CALL drag_noro_strato(0,klon,klev,phys_tstep,paprs,pplay, &
3941	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3942	igwd,idx,itest, &
3943	t_seri, u_seri, v_seri, &
3944	zulow, zvlow, zustrdr, zvstrdr, &
3945	d_t_oro, d_u_oro, d_v_oro)
3946
3947	ELSE
3948	CALL drag_noro(klon,klev,phys_tstep,paprs,pplay, &
3949	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3950	igwd,idx,itest, &
3951	t_seri, u_seri, v_seri, &
3952	zulow, zvlow, zustrdr, zvstrdr, &
3953	d_t_oro, d_u_oro, d_v_oro)
3954	ENDIF
3955	!
3956	! ajout des tendances
3957	!-----------------------------------------------------------------------
3958	! ajout des tendances de la trainee de l'orographie
3959	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro', &
3960	abortphy,flag_inhib_tend,itap,0)
3961	call prt_enerbil('oro',itap)
3962	!----------------------------------------------------------------------
3963	!
3964	ENDIF ! fin de test sur ok_orodr
3965	!
3966	IF (mydebug) THEN
3967	CALL writefield_phy('u_seri',u_seri,nbp_lev)
3968	CALL writefield_phy('v_seri',v_seri,nbp_lev)
3969	CALL writefield_phy('t_seri',t_seri,nbp_lev)
3970	CALL writefield_phy('q_seri',q_seri,nbp_lev)
3971	ENDIF
3972
3973	IF (ok_orolf) THEN
3974	!
3975	! selection des points pour lesquels le shema est actif:
3976	igwd=0
3977	DO i=1,klon
3978	itest(i)=0
3979	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3980	itest(i)=1
3981	igwd=igwd+1
3982	idx(igwd)=i
3983	ENDIF
3984	ENDDO
3985	! igwdim=MAX(1,igwd)
3986	!
3987	IF (ok_strato) THEN
3988
3989	CALL lift_noro_strato(klon,klev,phys_tstep,paprs,pplay, &
3990	latitude_deg,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3991	igwd,idx,itest, &
3992	t_seri, u_seri, v_seri, &
3993	zulow, zvlow, zustrli, zvstrli, &
3994	d_t_lif, d_u_lif, d_v_lif )
3995
3996	ELSE
3997	CALL lift_noro(klon,klev,phys_tstep,paprs,pplay, &
3998	latitude_deg,zmea,zstd,zpic, &
3999	itest, &
4000	t_seri, u_seri, v_seri, &
4001	zulow, zvlow, zustrli, zvstrli, &
4002	d_t_lif, d_u_lif, d_v_lif)
4003	ENDIF
4004
4005	! ajout des tendances de la portance de l'orographie
4006	CALL add_phys_tend(d_u_lif, d_v_lif, d_t_lif, dq0, dql0, dqi0, paprs, &
4007	'lif', abortphy,flag_inhib_tend,itap,0)
4008	call prt_enerbil('lif',itap)
4009	ENDIF ! fin de test sur ok_orolf
4010
4011	IF (ok_hines) then
4012	! HINES GWD PARAMETRIZATION
4013	east_gwstress=0.
4014	west_gwstress=0.
4015	du_gwd_hines=0.
4016	dv_gwd_hines=0.
4017	CALL hines_gwd(klon, klev, phys_tstep, paprs, pplay, latitude_deg, t_seri, &
4018	u_seri, v_seri, zustr_gwd_hines, zvstr_gwd_hines, d_t_hin, &
4019	du_gwd_hines, dv_gwd_hines)
4020	zustr_gwd_hines=0.
4021	zvstr_gwd_hines=0.
4022	DO k = 1, klev
4023	zustr_gwd_hines(:)=zustr_gwd_hines(:)+ du_gwd_hines(:, k)/phys_tstep &
4024	* (paprs(:, k)-paprs(:, k+1))/rg
4025	zvstr_gwd_hines(:)=zvstr_gwd_hines(:)+ dv_gwd_hines(:, k)/phys_tstep &
4026	* (paprs(:, k)-paprs(:, k+1))/rg
4027	ENDDO
4028
4029	d_t_hin(:, :)=0.
4030	CALL add_phys_tend(du_gwd_hines, dv_gwd_hines, d_t_hin, dq0, dql0, &
4031	dqi0, paprs, 'hin', abortphy,flag_inhib_tend,itap,0)
4032	call prt_enerbil('hin',itap)
4033	ENDIF
4034
4035	IF (.not. ok_hines .and. ok_gwd_rando) then
4036	CALL acama_GWD_rando(PHYS_TSTEP, pplay, latitude_deg, t_seri, u_seri, &
4037	v_seri, rot, zustr_gwd_front, zvstr_gwd_front, du_gwd_front, &
4038	dv_gwd_front, east_gwstress, west_gwstress)
4039	zustr_gwd_front=0.
4040	zvstr_gwd_front=0.
4041	DO k = 1, klev
4042	zustr_gwd_front(:)=zustr_gwd_front(:)+ du_gwd_front(:, k)/phys_tstep &
4043	* (paprs(:, k)-paprs(:, k+1))/rg
4044	zvstr_gwd_front(:)=zvstr_gwd_front(:)+ dv_gwd_front(:, k)/phys_tstep &
4045	* (paprs(:, k)-paprs(:, k+1))/rg
4046	ENDDO
4047
4048	CALL add_phys_tend(du_gwd_front, dv_gwd_front, dt0, dq0, dql0, dqi0, &
4049	paprs, 'front_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4050	call prt_enerbil('front_gwd_rando',itap)
4051	ENDIF
4052
4053	IF (ok_gwd_rando) THEN
4054	CALL FLOTT_GWD_rando(PHYS_TSTEP, pplay, t_seri, u_seri, v_seri, &
4055	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
4056	du_gwd_rando, dv_gwd_rando, east_gwstress, west_gwstress)
4057	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0, dqi0, &
4058	paprs, 'flott_gwd_rando', abortphy,flag_inhib_tend,itap,0)
4059	call prt_enerbil('flott_gwd_rando',itap)
4060	zustr_gwd_rando=0.
4061	zvstr_gwd_rando=0.
4062	DO k = 1, klev
4063	zustr_gwd_rando(:)=zustr_gwd_rando(:)+ du_gwd_rando(:, k)/phys_tstep &
4064	* (paprs(:, k)-paprs(:, k+1))/rg
4065	zvstr_gwd_rando(:)=zvstr_gwd_rando(:)+ dv_gwd_rando(:, k)/phys_tstep &
4066	* (paprs(:, k)-paprs(:, k+1))/rg
4067	ENDDO
4068	ENDIF
4069
4070	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
4071
4072	IF (mydebug) THEN
4073	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4074	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4075	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4076	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4077	ENDIF
4078
4079	DO i = 1, klon
4080	zustrph(i)=0.
4081	zvstrph(i)=0.
4082	ENDDO
4083	DO k = 1, klev
4084	DO i = 1, klon
4085	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/phys_tstep* &
4086	(paprs(i,k)-paprs(i,k+1))/rg
4087	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/phys_tstep* &
4088	(paprs(i,k)-paprs(i,k+1))/rg
4089	ENDDO
4090	ENDDO
4091	!
4092	!IM calcul composantes axiales du moment angulaire et couple des montagnes
4093	!
4094	IF (is_sequential .and. ok_orodr) THEN
4095	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
4096	ra,rg,romega, &
4097	latitude_deg,longitude_deg,pphis, &
4098	zustrdr,zustrli,zustrph, &
4099	zvstrdr,zvstrli,zvstrph, &
4100	paprs,u,v, &
4101	aam, torsfc)
4102	ENDIF
4103	!IM cf. FLott END
4104	!DC Calcul de la tendance due au methane
4105	IF(ok_qch4) THEN
4106	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
4107	! ajout de la tendance d'humidite due au methane
4108	d_q_ch4_dtime(:,:) = d_q_ch4(:,:)*phys_tstep
4109	CALL add_phys_tend(du0, dv0, dt0, d_q_ch4_dtime, dql0, dqi0, paprs, &
4110	'q_ch4', abortphy,flag_inhib_tend,itap,0)
4111	d_q_ch4(:,:) = d_q_ch4_dtime(:,:)/phys_tstep
4112	ENDIF
4113	!
4114	!
4115
4116	!===============================================================
4117	! Additional tendency of TKE due to orography
4118	!===============================================================
4119	!
4120	! Inititialization
4121	!------------------
4122
4123
4124
4125	addtkeoro=0
4126	CALL getin_p('addtkeoro',addtkeoro)
4127
4128	IF (prt_level.ge.5) &
4129	print*,'addtkeoro', addtkeoro
4130
4131	alphatkeoro=1.
4132	CALL getin_p('alphatkeoro',alphatkeoro)
4133	alphatkeoro=min(max(0.,alphatkeoro),1.)
4134
4135	smallscales_tkeoro=.false.
4136	CALL getin_p('smallscales_tkeoro',smallscales_tkeoro)
4137
4138
4139	dtadd(:,:)=0.
4140	duadd(:,:)=0.
4141	dvadd(:,:)=0.
4142
4143
4144
4145	! Choices for addtkeoro:
4146	! ** 0 no TKE tendency from orography
4147	! ** 1 we include a fraction alphatkeoro of the whole tendency duoro
4148	! ** 2 we include a fraction alphatkeoro of the gravity wave part of duoro
4149	!
4150
4151	IF (addtkeoro .GT. 0 .AND. ok_orodr ) THEN
4152	! -------------------------------------------
4153
4154
4155	! selection des points pour lesquels le schema est actif:
4156
4157
4158
4159	IF (addtkeoro .EQ. 1 ) THEN
4160
4161	duadd(:,:)=alphatkeoro*d_u_oro(:,:)
4162	dvadd(:,:)=alphatkeoro*d_v_oro(:,:)
4163
4164	ELSE IF (addtkeoro .EQ. 2) THEN
4165
4166
4167
4168	IF (smallscales_tkeoro) THEN
4169	igwd=0
4170	DO i=1,klon
4171	itest(i)=0
4172	! Etienne: ici je prends en compte plus de relief que la routine drag_noro_strato
4173	! car on peut s'attendre a ce que les petites echelles produisent aussi de la TKE
4174	! Mais attention, cela ne va pas dans le sens de la conservation de l'energie!
4175	IF (zstd(i).GT.1.0) THEN
4176	itest(i)=1
4177	igwd=igwd+1
4178	idx(igwd)=i
4179	ENDIF
4180	ENDDO
4181
4182	ELSE
4183
4184	igwd=0
4185	DO i=1,klon
4186	itest(i)=0
4187	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
4188	itest(i)=1
4189	igwd=igwd+1
4190	idx(igwd)=i
4191	ENDIF
4192	ENDDO
4193
4194	END IF
4195
4196
4197
4198
4199	CALL drag_noro_strato(addtkeoro,klon,klev,phys_tstep,paprs,pplay, &
4200	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
4201	igwd,idx,itest, &
4202	t_seri, u_seri, v_seri, &
4203	zulow, zvlow, zustrdr, zvstrdr, &
4204	d_t_oro_gw, d_u_oro_gw, d_v_oro_gw)
4205
4206	zustrdr(:)=0.
4207	zvstrdr(:)=0.
4208	zulow(:)=0.
4209	zvlow(:)=0.
4210
4211	duadd(:,:)=alphatkeoro*d_u_oro_gw(:,:)
4212	dvadd(:,:)=alphatkeoro*d_v_oro_gw(:,:)
4213	END IF
4214
4215
4216
4217	! TKE update from subgrid temperature and wind tendencies
4218	!----------------------------------------------------------
4219	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4220
4221
4222	CALL tend_to_tke(pdtphys,paprs,exner,t_seri,u_seri,v_seri,dtadd,duadd,dvadd,pbl_tke)
4223
4224
4225
4226	ENDIF
4227	! -----
4228	!===============================================================
4229
4230
4231
4232	!====================================================================
4233	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
4234	!====================================================================
4235	! Abderrahmane 24.08.09
4236
4237	IF (ok_cosp) THEN
4238	! adeclarer
4239	#ifdef CPP_COSP
4240	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/phys_tstep)).EQ.0) THEN
4241
4242	IF (prt_level .GE.10) THEN
4243	print*,'freq_cosp',freq_cosp
4244	ENDIF
4245	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
4246	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
4247	! s ref_liq,ref_ice
4248	CALL phys_cosp(itap,phys_tstep,freq_cosp, &
4249	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
4250	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, missing_val, &
4251	klon,klev,longitude_deg,latitude_deg,presnivs,overlap, &
4252	JrNt,ref_liq,ref_ice, &
4253	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
4254	zu10m,zv10m,pphis, &
4255	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
4256	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
4257	prfl(:,1:klev),psfl(:,1:klev), &
4258	pmflxr(:,1:klev),pmflxs(:,1:klev), &
4259	mr_ozone,cldtau, cldemi)
4260
4261	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
4262	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
4263	! M clMISR,
4264	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
4265	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
4266
4267	ENDIF
4268
4269	#endif
4270	ENDIF !ok_cosp
4271
4272
4273	! Marine
4274
4275	IF (ok_airs) then
4276
4277	IF (itap.eq.1.or.MOD(itap,NINT(freq_airs/phys_tstep)).EQ.0) THEN
4278	write(,) 'je vais appeler simu_airs, ok_airs, freq_airs=', ok_airs, freq_airs
4279	CALL simu_airs(itap,rneb, t_seri, cldemi, fiwc, ref_ice, pphi, pplay, paprs,&
4280	& map_prop_hc,map_prop_hist,&
4281	& map_emis_hc,map_iwp_hc,map_deltaz_hc,map_pcld_hc,map_tcld_hc,&
4282	& map_emis_Cb,map_pcld_Cb,map_tcld_Cb,&
4283	& map_emis_ThCi,map_pcld_ThCi,map_tcld_ThCi,&
4284	& map_emis_Anv,map_pcld_Anv,map_tcld_Anv,&
4285	& map_emis_hist,map_iwp_hist,map_deltaz_hist,map_rad_hist,&
4286	& map_ntot,map_hc,map_hist,&
4287	& map_Cb,map_ThCi,map_Anv,&
4288	& alt_tropo )
4289	ENDIF
4290
4291	ENDIF ! ok_airs
4292
4293
4294	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4295	!AA
4296	!AA Installation de l'interface online-offline pour traceurs
4297	!AA
4298	!====================================================================
4299	! Calcul des tendances traceurs
4300	!====================================================================
4301	!
4302
4303	IF (type_trac=='repr') THEN
4304	sh_in(:,:) = q_seri(:,:)
4305	ELSE
4306	sh_in(:,:) = qx(:,:,ivap)
4307	ch_in(:,:) = qx(:,:,iliq)
4308	ENDIF
4309
4310	IF (iflag_phytrac == 1 ) THEN
4311
4312	#ifdef CPP_Dust
4313	CALL phytracr_spl ( debut,lafin , jD_cur,jH_cur,iflag_con, & ! I
4314	pdtphys,ftsol, & ! I
4315	t,q_seri,paprs,pplay,RHcl, & ! I
4316	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, & ! I
4317	coefh(1:klon,1:klev,is_ave), cdragh, cdragm, u1, v1, & ! I
4318	u_seri, v_seri, latitude_deg, longitude_deg, &
4319	pphis,pctsrf,pmflxr,pmflxs,prfl,psfl, & ! I
4320	da,phi,phi2,d1a,dam,mp,ep,sigd,sij,clw,elij, & ! I
4321	epmlmMm,eplaMm,upwd,dnwd,itop_con,ibas_con, & ! I
4322	ev,wdtrainA, wdtrainM,wght_cvfd, & ! I
4323	fm_therm, entr_therm, rneb, & ! I
4324	beta_prec_fisrt,beta_prec, & !I
4325	zu10m,zv10m,wstar,ale_bl,ale_wake, & ! I
4326	d_tr_dyn,tr_seri)
4327
4328	#else
4329
4330	CALL phytrac ( &
4331	itap, days_elapsed+1, jH_cur, debut, &
4332	lafin, phys_tstep, u, v, t, &
4333	paprs, pplay, pmfu, pmfd, &
4334	pen_u, pde_u, pen_d, pde_d, &
4335	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
4336	u1, v1, ftsol, pctsrf, &
4337	zustar, zu10m, zv10m, &
4338	wstar(:,is_ave), ale_bl, ale_wake, &
4339	latitude_deg, longitude_deg, &
4340	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
4341	presnivs, pphis, pphi, albsol1, &
4342	sh_in, ch_in, rhcl, cldfra, rneb, &
4343	diafra, cldliq, itop_con, ibas_con, &
4344	pmflxr, pmflxs, prfl, psfl, &
4345	da, phi, mp, upwd, &
4346	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
4347	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
4348	ev, ep, epmlmMm, eplaMm, & !<<RomP
4349	dnwd, aerosol_couple, flxmass_w, &
4350	tau_aero, piz_aero, cg_aero, ccm, &
4351	rfname, &
4352	d_tr_dyn, & !<<RomP
4353	tr_seri)
4354	#endif
4355	ENDIF ! (iflag_phytrac=1)
4356
4357	IF (offline) THEN
4358
4359	IF (prt_level.ge.9) &
4360	print*,'Attention on met a 0 les thermiques pour phystoke'
4361	CALL phystokenc ( &
4362	nlon,klev,pdtphys,longitude_deg,latitude_deg, &
4363	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
4364	fm_therm,entr_therm, &
4365	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
4366	frac_impa, frac_nucl, &
4367	pphis,cell_area,phys_tstep,itap, &
4368	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
4369
4370
4371	ENDIF
4372
4373	!
4374	! Calculer le transport de l'eau et de l'energie (diagnostique)
4375	!
4376	CALL transp (paprs,zxtsol, &
4377	t_seri, q_seri, u_seri, v_seri, zphi, &
4378	ve, vq, ue, uq)
4379	!
4380	!IM global posePB BEG
4381	IF(1.EQ.0) THEN
4382	!
4383	CALL transp_lay (paprs,zxtsol, &
4384	t_seri, q_seri, u_seri, v_seri, zphi, &
4385	ve_lay, vq_lay, ue_lay, uq_lay)
4386	!
4387	ENDIF !(1.EQ.0) THEN
4388	!IM global posePB END
4389	! Accumuler les variables a stocker dans les fichiers histoire:
4390	!
4391
4392	!================================================================
4393	! Conversion of kinetic and potential energy into heat, for
4394	! parameterisation of subgrid-scale motions
4395	!================================================================
4396
4397	d_t_ec(:,:)=0.
4398	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
4399	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap),qx(:,:,iliq),qx(:,:,isol), &
4400	u_seri,v_seri,t_seri,q_seri,ql_seri,qs_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
4401	zmasse,exner,d_t_ec)
4402	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
4403
4404	!=======================================================================
4405	! SORTIES
4406	!=======================================================================
4407	!
4408	!IM initialisation + calculs divers diag AMIP2
4409	!
4410	include "calcul_divers.h"
4411	!
4412	!IM Interpolation sur les niveaux de pression du NMC
4413	! -------------------------------------------------
4414	!
4415	include "calcul_STDlev.h"
4416	!
4417	! slp sea level pressure derived from Arpege-IFS : CALL ctstar + CALL pppmer
4418	CALL diag_slp(klon,t_seri,paprs,pplay,pphis,ptstar,pt0,slp)
4419	!
4420	!cc prw = eau precipitable
4421	! prlw = colonne eau liquide
4422	! prlw = colonne eau solide
4423	prw(:) = 0.
4424	prlw(:) = 0.
4425	prsw(:) = 0.
4426	DO k = 1, klev
4427	prw(:) = prw(:) + q_seri(:,k)*zmasse(:,k)
4428	prlw(:) = prlw(:) + ql_seri(:,k)*zmasse(:,k)
4429	prsw(:) = prsw(:) + qs_seri(:,k)*zmasse(:,k)
4430	ENDDO
4431	!
4432	IF (type_trac == 'inca') THEN
4433	#ifdef INCA
4434	CALL VTe(VTphysiq)
4435	CALL VTb(VTinca)
4436
4437	CALL chemhook_end ( &
4438	phys_tstep, &
4439	pplay, &
4440	t_seri, &
4441	tr_seri, &
4442	nbtr, &
4443	paprs, &
4444	q_seri, &
4445	cell_area, &
4446	pphi, &
4447	pphis, &
4448	zx_rh, &
4449	aps, bps)
4450
4451	CALL VTe(VTinca)
4452	CALL VTb(VTphysiq)
4453	#endif
4454	ENDIF
4455
4456
4457	!
4458	! Convertir les incrementations en tendances
4459	!
4460	IF (prt_level .GE.10) THEN
4461	print *,'Convertir les incrementations en tendances '
4462	ENDIF
4463	!
4464	IF (mydebug) THEN
4465	CALL writefield_phy('u_seri',u_seri,nbp_lev)
4466	CALL writefield_phy('v_seri',v_seri,nbp_lev)
4467	CALL writefield_phy('t_seri',t_seri,nbp_lev)
4468	CALL writefield_phy('q_seri',q_seri,nbp_lev)
4469	ENDIF
4470
4471	DO k = 1, klev
4472	DO i = 1, klon
4473	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / phys_tstep
4474	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / phys_tstep
4475	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / phys_tstep
4476	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / phys_tstep
4477	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / phys_tstep
4478	!CR: on ajoute le contenu en glace
4479	IF (nqo.eq.3) THEN
4480	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / phys_tstep
4481	ENDIF
4482	ENDDO
4483	ENDDO
4484	!
4485	!CR: nb de traceurs eau: nqo
4486	! IF (nqtot.GE.3) THEN
4487	IF (nqtot.GE.(nqo+1)) THEN
4488	! DO iq = 3, nqtot
4489	DO iq = nqo+1, nqtot
4490	DO k = 1, klev
4491	DO i = 1, klon
4492	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / phys_tstep
4493	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / phys_tstep
4494	ENDDO
4495	ENDDO
4496	ENDDO
4497	ENDIF
4498	!
4499	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4500	!IM global posePB include "write_bilKP_ins.h"
4501	!IM global posePB include "write_bilKP_ave.h"
4502	!
4503
4504	!--OB mass fixer
4505	!--profile is corrected to force mass conservation of water
4506	IF (mass_fixer) THEN
4507	qql2(:)=0.0
4508	DO k = 1, klev
4509	qql2(:)=qql2(:)+(q_seri(:,k)+ql_seri(:,k)+qs_seri(:,k))*zmasse(:,k)
4510	ENDDO
4511	DO i = 1, klon
4512	!--compute ratio of what q+ql should be with conservation to what it is
4513	corrqql=(qql1(i)+(evap(i)-rain_fall(i)-snow_fall(i))*pdtphys)/qql2(i)
4514	DO k = 1, klev
4515	q_seri(i,k) =q_seri(i,k)*corrqql
4516	ql_seri(i,k)=ql_seri(i,k)*corrqql
4517	ENDDO
4518	ENDDO
4519	ENDIF
4520	!--fin mass fixer
4521
4522	! Sauvegarder les valeurs de t et q a la fin de la physique:
4523	!
4524	u_ancien(:,:) = u_seri(:,:)
4525	v_ancien(:,:) = v_seri(:,:)
4526	t_ancien(:,:) = t_seri(:,:)
4527	q_ancien(:,:) = q_seri(:,:)
4528	ql_ancien(:,:) = ql_seri(:,:)
4529	qs_ancien(:,:) = qs_seri(:,:)
4530	CALL water_int(klon,klev,q_ancien,zmasse,prw_ancien)
4531	CALL water_int(klon,klev,ql_ancien,zmasse,prlw_ancien)
4532	CALL water_int(klon,klev,qs_ancien,zmasse,prsw_ancien)
4533	! !! RomP >>>
4534	!CR: nb de traceurs eau: nqo
4535	IF (nqtot.GT.nqo) THEN
4536	DO iq = nqo+1, nqtot
4537	tr_ancien(:,:,iq-nqo) = tr_seri(:,:,iq-nqo)
4538	ENDDO
4539	ENDIF
4540	! !! RomP <<<
4541	!==========================================================================
4542	! Sorties des tendances pour un point particulier
4543	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4544	! pour le debug
4545	! La valeur de igout est attribuee plus haut dans le programme
4546	!==========================================================================
4547
4548	IF (prt_level.ge.1) THEN
4549	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4550	write(lunout,*) &
4551	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4552	write(lunout,*) &
4553	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4554	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4555	pctsrf(igout,is_sic)
4556	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4557	DO k=1,klev
4558	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4559	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4560	d_t_eva(igout,k)
4561	ENDDO
4562	write(lunout,*) 'cool,heat'
4563	DO k=1,klev
4564	write(lunout,*) cool(igout,k),heat(igout,k)
4565	ENDDO
4566
4567	!jyg< (En attendant de statuer sur le sort de d_t_oli)
4568	!jyg! write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4569	!jyg! do k=1,klev
4570	!jyg! write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4571	!jyg! d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4572	!jyg! enddo
4573	write(lunout,*) 'd_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4574	DO k=1,klev
4575	write(lunout,*) d_t_vdf(igout,k), &
4576	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4577	ENDDO
4578	!>jyg
4579
4580	write(lunout,*) 'd_ps ',d_ps(igout)
4581	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4582	DO k=1,klev
4583	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4584	d_qx(igout,k,1),d_qx(igout,k,2)
4585	ENDDO
4586	ENDIF
4587
4588	!============================================================
4589	! Calcul de la temperature potentielle
4590	!============================================================
4591	DO k = 1, klev
4592	DO i = 1, klon
4593	!JYG/IM theta en debut du pas de temps
4594	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4595	!JYG/IM theta en fin de pas de temps de physique
4596	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4597	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers
4598	! MPL 20130625
4599	! fth_fonctions.F90 et parkind1.F90
4600	! sinon thetal=theta
4601	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4602	! : ql_seri(i,k))
4603	thetal(i,k)=theta(i,k)
4604	ENDDO
4605	ENDDO
4606	!
4607
4608	! 22.03.04 BEG
4609	!=============================================================
4610	! Ecriture des sorties
4611	!=============================================================
4612	#ifdef CPP_IOIPSL
4613
4614	! Recupere des varibles calcule dans differents modules
4615	! pour ecriture dans histxxx.nc
4616
4617	! Get some variables from module fonte_neige_mod
4618	CALL fonte_neige_get_vars(pctsrf, &
4619	zxfqcalving, zxfqfonte, zxffonte, zxrunofflic)
4620
4621
4622	!=============================================================
4623	! Separation entre thermiques et non thermiques dans les sorties
4624	! de fisrtilp
4625	!=============================================================
4626
4627	IF (iflag_thermals>=1) THEN
4628	d_t_lscth=0.
4629	d_t_lscst=0.
4630	d_q_lscth=0.
4631	d_q_lscst=0.
4632	DO k=1,klev
4633	DO i=1,klon
4634	IF (ptconvth(i,k)) THEN
4635	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4636	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4637	ELSE
4638	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4639	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4640	ENDIF
4641	ENDDO
4642	ENDDO
4643
4644	DO i=1,klon
4645	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4646	plul_th(i)=prfl(i,1)+psfl(i,1)
4647	ENDDO
4648	ENDIF
4649
4650	!On effectue les sorties:
4651
4652	#ifdef CPP_Dust
4653	CALL phys_output_write_spl(itap, pdtphys, paprs, pphis, &
4654	pplay, lmax_th, aerosol_couple, &
4655	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
4656	ptconv, read_climoz, clevSTD, &
4657	ptconvth, d_t, qx, d_qx, d_tr_dyn, zmasse, &
4658	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4659	#else
4660	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
4661	pplay, lmax_th, aerosol_couple, &
4662	ok_ade, ok_aie, ivap, iliq, isol, new_aod, &
4663	ok_sync, ptconv, read_climoz, clevSTD, &
4664	ptconvth, d_u, d_t, qx, d_qx, zmasse, &
4665	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4666	#endif
4667
4668	#ifndef CPP_XIOS
4669	CALL write_paramLMDZ_phy(itap,nid_ctesGCM,ok_sync)
4670	#endif
4671
4672	#endif
4673
4674
4675	!====================================================================
4676	! Arret du modele apres hgardfou en cas de detection d'un
4677	! plantage par hgardfou
4678	!====================================================================
4679
4680	IF (abortphy==1) THEN
4681	abort_message ='Plantage hgardfou'
4682	CALL abort_physic (modname,abort_message,1)
4683	ENDIF
4684
4685	! 22.03.04 END
4686	!
4687	!====================================================================
4688	! Si c'est la fin, il faut conserver l'etat de redemarrage
4689	!====================================================================
4690	!
4691
4692	IF (lafin) THEN
4693	itau_phy = itau_phy + itap
4694	CALL phyredem ("restartphy.nc")
4695	! open(97,form="unformatted",file="finbin")
4696	! write(97) u_seri,v_seri,t_seri,q_seri
4697	! close(97)
4698	!$OMP MASTER
4699	IF (read_climoz >= 1) THEN
4700	IF (is_mpi_root) THEN
4701	CALL nf95_close(ncid_climoz)
4702	ENDIF
4703	DEALLOCATE(press_edg_climoz) ! pointer
4704	DEALLOCATE(press_cen_climoz) ! pointer
4705	ENDIF
4706	!$OMP END MASTER
4707	ENDIF
4708
4709	! first=.false.
4710
4711
4712	END SUBROUTINE physiq
4713
4714	END MODULE physiq_mod

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

Download in other formats:

Original Format