Context Navigation

physiq_mod.F90 @ 3370

Last change on this file since 3370 was 3370, checked in by idelkadi, 6 years ago
Implementation de COSPv2 dans LMDZ.
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: 172.1 KB

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

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

Download in other formats:

Original Format