Context Navigation

physiq_mod.F90 @ 5214

Last change on this file since 5214 was 2787, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2727:2785 into testing branch
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: 153.5 KB

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

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

Download in other formats:

Original Format