Context Navigation

physiq_mod.F90 @ 2630

Last change on this file since 2630 was 2630, checked in by fhourdin, 8 years ago
Importation du modèle d'aérosols de Boucher, Escribano et al.
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: 162.1 KB

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

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

Download in other formats:

Original Format