Context Navigation

physiq_mod.F90 @ 3324

Last change on this file since 3324 was 3324, checked in by musat, 6 years ago
tau_gl declare dans clesphys.h et lu dans physiq_mod. Par defaut tau_gl=86400.*5. pour coherence avec versions anciennes.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 170.0 KB

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

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

Download in other formats:

Original Format