Context Navigation

physiq_mod.F90 @ 3629

Last change on this file since 3629 was 3629, checked in by acozic, 4 years ago

Add new grid, new axis and new variables for cmip protocole and dr2xml

field_group id="coord_hyb"
grid_ref="klevp1_bnds"
grid_ref="klev_bnds"
domain id="greordered"
axis id="axis_lat"
axis id="bnds"
axis id="klevp1"
axis id="klev"

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: 161.4 KB

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

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

Download in other formats:

Original Format