Context Navigation

physiq.F90 @ 2340

Last change on this file since 2340 was 2340, checked in by Ehouarn Millour, 9 years ago
Bug fix: Adding the tendencies from wakes to dry adjustement tendencies should only be done if wakes are computed. EM + JYG
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.8 KB

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

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

Download in other formats:

Original Format