Context Navigation

physiq.F90 @ 3825

Last change on this file since 3825 was 3825, checked in by ymipsl, 10 years ago

Reorganize geometry and grid modules. Prepare physics for unstructutured grid support. Simplify initialization of physics from dynamic.
Compiled only with dynd3dmem, but not tested for moment.

File size: 136.7 KB

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

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

Download in other formats:

Original Format