Context Navigation

physiq.F90 @ 2343

Last change on this file since 2343 was 2343, checked in by Ehouarn Millour, 9 years ago
Another step towards a clean separation between physics and dynamics: adapted read_pstoke.F90, read_pstoke0.F90 initphysto.F90 and phystokenc.F90 (now module phystokenc_mod.F90) to not explicitely include/use file/modules from the dynamics. In the process, added module "time_phylmdz_mod.F90" in the physics, which contains the information otherwise found in "temps.h" (which is in the dynamics) and should be used instead. EM
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: 141.9 KB

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

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

Download in other formats:

Original Format