Context Navigation

physiq_mod.F90 @ 3339

Last change on this file since 3339 was 3339, checked in by acozic, 6 years ago
merge with rev 3338 of the trunk (add a flag chemistry_couple to use Ozone offline or Ozone calcul by Inca)
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 170.2 KB

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

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

Download in other formats:

Original Format