Context Navigation

physiq_mod.F90 @ 3086

Last change on this file since 3086 was 3086, checked in by dcugnet, 6 years ago

Modification of the ozone interpolation routine to ensure more realistic
profiles, in particular when the tropopause adjustment is activated, as
well as an improved accuracy on the chemical tropopause of streched ozone
field:

stretching function has no longer a trapezoidal shape, but is a hat

function, which top is at the target (ldmz) tropopause, with lower /
upper bounds respectively lower / higher than the lowest / highest
tropopause (lmdz and file). This helps to avoid stretching the profile
too far from the strictly required zone.

revert to a classical linear vertical interpolation on log-pressure

because the stairs-shape effect of the conservative interpolation can not
be avoided, even at second order method. Note that the second order
conservative method is still used for Cariolle linearized ozone
coefficients interpolation: the regr_pr_time routine now switches from
one method to the other depending on flag "Ploc" (pressures must be given
at layers interfaces - Ploc=='I' - for conservative method and layers
centers - Ploc=='C' - for the regular one).

force the stretching function to be discretized on at least three

levels to avoid the effect of the stretching to be reduced to a single
point shifting of the ozone profile when discretized.

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

File size: 164.8 KB

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

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

Download in other formats:

Original Format