Context Navigation

physiq_mod.F90 @ 3908

Last change on this file since 3908 was 3908, checked in by idelkadi, 3 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

Inclusion of the ecrad directory containing the sources of the ECRAD code
- interface routine : radiation_scheme.F90
Adaptation of compilation scripts :
- compilation under CPP key CPP_ECRAD
- compilation with option "-rad ecard" or "-ecard true"
- The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
Runing LMDZ simulations with ecrad, you need :
- logical key iflag_rrtm = 2 in physiq.def
- namelist_ecrad (DefLists?)
- the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Id

File size: 186.0 KB

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

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

Download in other formats:

Original Format