Context Navigation

physiq_mod.F90 @ 3239

Last change on this file since 3239 was 3239, checked in by Laurent Fairhead, 6 years ago
Moving some xios_field_is_active tests as they could not give the right results at their previous location (the context was not closed yet). Changing some variables levels in the histmth.xml file to prevent the unecessary multiple calls to the radiation scheme
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: 170.1 KB

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

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

Download in other formats:

Original Format