Context Navigation

physiq.F90 @ 2245

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

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

Download in other formats:

Original Format