Context Navigation

physiq_mod.F90 @ 3335

Last change on this file since 3335 was 3327, checked in by musat, 7 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: 171.0 KB

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

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

Download in other formats:

Original Format