Context Navigation

physiq.F90 @ 2234

Last change on this file since 2234 was 2231, checked in by oboucher, 9 years ago
Putting minimum values for SW and LW aerosol optical depth values Putting default value for aerosol single scattering albedo if no aerosol Doing this consistently in all cases covered by flag_aerosol and flag_aerosol_strat
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.0 KB

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

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

Download in other formats:

Original Format