Context Navigation

physiq.F90 @ 2227

Last change on this file since 2227 was 2227, checked in by Laurent Fairhead, 9 years ago

New ocean albedo.

To activate the new scheme, put iflag_albedo=1 in physiq.def

To activate chlorophyll concentration effect on albedo,
put ok_chlorophyll=y in def file

and download file named chlorophyll.nc
chlorophyll.nc has the same dimension as the model grid with 12 months data,
(i=lon, j=lat, L=1:12) and can be degraded from the original file of dimension
i=1:4320 , j=1:2160 , L=1:12
ada:/workgpfs/rech/psl/rpsl949/clima/chlor_seasonal_clim_seawifs.nc

For 96X96 resolution, chlorophyll.nc file is in
ada:/workgpfs/rech/psl/rpsl949/clima/chlorophyll.nc

Baek

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 2227 2015-03-12 11:07:43Z fairhead $
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(:,:,:,:) = 0.
3029	piz_aero(:,:,:,:) = 0.
3030	cg_aero(:,:,:,:) = 0.
3031	ELSE
3032	tau_aero_sw_rrtm(:,:,:,:)=0.0
3033	piz_aero_sw_rrtm(:,:,:,:)=0.0
3034	cg_aero_sw_rrtm(:,:,:,:)=0.0
3035	ENDIF
3036	ENDIF
3037	!
3038	!--STRAT AEROSOL
3039	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
3040	IF (flag_aerosol_strat) THEN
3041	IF (prt_level .GE.10) THEN
3042	PRINT *,'appel a readaerosolstrat', mth_cur
3043	ENDIF
3044	IF (iflag_rrtm.EQ.0) THEN
3045	CALL readaerosolstrato(debut)
3046	ELSE
3047	#ifdef CPP_RRTM
3048	CALL readaerosolstrato_rrtm(debut)
3049	#else
3050
3051	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
3052	call abort_gcm(modname,abort_message,1)
3053	#endif
3054	ENDIF
3055	ENDIF
3056	!--fin STRAT AEROSOL
3057
3058
3059	! On prend la somme des fractions nuageuses et des contenus en eau
3060
3061	if (iflag_cldth>=5) then
3062
3063	do k=1,klev
3064	ptconvth(:,k)=fm_therm(:,k+1)>0.
3065	enddo
3066
3067	if (iflag_coupl==4) then
3068
3069	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3070	! convectives et lsc dans la partie des thermiques
3071	! Le controle par iflag_coupl est peut etre provisoire.
3072	do k=1,klev
3073	do i=1,klon
3074	if (ptconv(i,k).and.ptconvth(i,k)) then
3075	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3076	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3077	else if (ptconv(i,k)) then
3078	cldfra(i,k)=rnebcon(i,k)
3079	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3080	endif
3081	enddo
3082	enddo
3083
3084	else if (iflag_coupl==5) then
3085	do k=1,klev
3086	do i=1,klon
3087	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3088	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3089	enddo
3090	enddo
3091
3092	else
3093
3094	! Si on est sur un point touche par la convection profonde et pas
3095	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
3096	! de la convection profonde.
3097
3098	!IM/FH: 2011/02/23
3099	! definition des points sur lesquels ls thermiques sont actifs
3100
3101	do k=1,klev
3102	do i=1,klon
3103	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
3104	cldfra(i,k)=rnebcon(i,k)
3105	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3106	endif
3107	enddo
3108	enddo
3109
3110	endif
3111
3112	else
3113
3114	! Ancienne version
3115	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3116	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3117	endif
3118
3119	ENDIF
3120
3121	! plulsc(:)=0.
3122	! do k=1,klev,-1
3123	! do i=1,klon
3124	! zzz=prfl(:,k)+psfl(:,k)
3125	! if (.not.ptconvth.zzz.gt.0.)
3126	! enddo prfl, psfl,
3127	! enddo
3128	!
3129	! 2. NUAGES STARTIFORMES
3130	!
3131	IF (ok_stratus) THEN
3132	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3133	DO k = 1, klev
3134	DO i = 1, klon
3135	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3136	cldliq(i,k) = dialiq(i,k)
3137	cldfra(i,k) = diafra(i,k)
3138	ENDIF
3139	ENDDO
3140	ENDDO
3141	ENDIF
3142	!
3143	! Precipitation totale
3144	!
3145	DO i = 1, klon
3146	rain_fall(i) = rain_con(i) + rain_lsc(i)
3147	snow_fall(i) = snow_con(i) + snow_lsc(i)
3148	ENDDO
3149	!IM
3150	IF (ip_ebil_phy.ge.2) THEN
3151	ztit="after diagcld"
3152	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3153	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3154	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3155	call diagphy(airephy,ztit,ip_ebil_phy &
3156	, zero_v, zero_v, zero_v, zero_v, zero_v &
3157	, zero_v, zero_v, zero_v, ztsol &
3158	, d_h_vcol, d_qt, d_ec &
3159	, fs_bound, fq_bound )
3160	END IF
3161	!
3162	! Calculer l'humidite relative pour diagnostique
3163	!
3164	DO k = 1, klev
3165	DO i = 1, klon
3166	zx_t = t_seri(i,k)
3167	IF (thermcep) THEN
3168	if (iflag_ice_thermo.eq.0) then
3169	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3170	else
3171	zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))
3172	endif
3173	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3174	zx_qs = MIN(0.5,zx_qs)
3175	zcor = 1./(1.-retv*zx_qs)
3176	zx_qs = zx_qs*zcor
3177	ELSE
3178	IF (zx_t.LT.t_coup) THEN
3179	zx_qs = qsats(zx_t)/pplay(i,k)
3180	ELSE
3181	zx_qs = qsatl(zx_t)/pplay(i,k)
3182	ENDIF
3183	ENDIF
3184	zx_rh(i,k) = q_seri(i,k)/zx_qs
3185	zqsat(i,k)=zx_qs
3186	ENDDO
3187	ENDDO
3188
3189	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3190	! equivalente a 2m (tpote) pour diagnostique
3191	!
3192	DO i = 1, klon
3193	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3194	IF (thermcep) THEN
3195	IF(zt2m(i).LT.RTT) then
3196	Lheat=RLSTT
3197	ELSE
3198	Lheat=RLVTT
3199	ENDIF
3200	ELSE
3201	IF (zt2m(i).LT.RTT) THEN
3202	Lheat=RLSTT
3203	ELSE
3204	Lheat=RLVTT
3205	ENDIF
3206	ENDIF
3207	tpote(i) = tpot(i)* &
3208	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3209	ENDDO
3210
3211	IF (type_trac == 'inca') THEN
3212	#ifdef INCA
3213	CALL VTe(VTphysiq)
3214	CALL VTb(VTinca)
3215	calday = REAL(days_elapsed + 1) + jH_cur
3216
3217	call chemtime(itap+itau_phy-1, date0, dtime, itap)
3218	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
3219	CALL AEROSOL_METEO_CALC( &
3220	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3221	prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
3222	END IF
3223
3224	zxsnow_dummy(:) = 0.0
3225
3226	CALL chemhook_begin (calday, &
3227	days_elapsed+1, &
3228	jH_cur, &
3229	pctsrf(1,1), &
3230	rlat, &
3231	rlon, &
3232	airephy, &
3233	paprs, &
3234	pplay, &
3235	coefh(1:klon,1:klev,is_ave), &
3236	pphi, &
3237	t_seri, &
3238	u, &
3239	v, &
3240	wo(:, :, 1), &
3241	q_seri, &
3242	zxtsol, &
3243	zxsnow_dummy, &
3244	solsw, &
3245	albsol1, &
3246	rain_fall, &
3247	snow_fall, &
3248	itop_con, &
3249	ibas_con, &
3250	cldfra, &
3251	iim, &
3252	jjm, &
3253	tr_seri, &
3254	ftsol, &
3255	paprs, &
3256	cdragh, &
3257	cdragm, &
3258	pctsrf, &
3259	pdtphys, &
3260	itap)
3261
3262	CALL VTe(VTinca)
3263	CALL VTb(VTphysiq)
3264	#endif
3265	END IF !type_trac = inca
3266	!
3267	! Calculer les parametres optiques des nuages et quelques
3268	! parametres pour diagnostiques:
3269	!
3270
3271	IF (aerosol_couple) THEN
3272	mass_solu_aero(:,:) = ccm(:,:,1)
3273	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3274	END IF
3275
3276	if (ok_newmicro) then
3277	IF (iflag_rrtm.NE.0) THEN
3278	#ifdef CPP_RRTM
3279	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3280	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 pour ok_cdnc'
3281	call abort_gcm(modname,abort_message,1)
3282	endif
3283	#else
3284
3285	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
3286	call abort_gcm(modname,abort_message,1)
3287	#endif
3288	ENDIF
3289	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3290	paprs, pplay, t_seri, cldliq, cldfra, &
3291	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3292	flwp, fiwp, flwc, fiwc, &
3293	mass_solu_aero, mass_solu_aero_pi, &
3294	cldtaupi, re, fl, ref_liq, ref_ice, &
3295	ref_liq_pi, ref_ice_pi)
3296	else
3297	CALL nuage (paprs, pplay, &
3298	t_seri, cldliq, cldfra, cldtau, cldemi, &
3299	cldh, cldl, cldm, cldt, cldq, &
3300	ok_aie, &
3301	mass_solu_aero, mass_solu_aero_pi, &
3302	bl95_b0, bl95_b1, &
3303	cldtaupi, re, fl)
3304	endif
3305	!
3306	!IM betaCRF
3307	!
3308	cldtaurad = cldtau
3309	cldtaupirad = cldtaupi
3310	cldemirad = cldemi
3311
3312	!
3313	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3314	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3315	!
3316	! global
3317	!
3318	DO k=1, klev
3319	DO i=1, klon
3320	if (pplay(i,k).GE.pfree) THEN
3321	beta(i,k) = beta_pbl
3322	else
3323	beta(i,k) = beta_free
3324	endif
3325	if (mskocean_beta) THEN
3326	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3327	endif
3328	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3329	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3330	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3331	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3332	ENDDO
3333	ENDDO
3334	!
3335	else
3336	!
3337	! regional
3338	!
3339	DO k=1, klev
3340	DO i=1,klon
3341	!
3342	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
3343	rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
3344	if (pplay(i,k).GE.pfree) THEN
3345	beta(i,k) = beta_pbl
3346	else
3347	beta(i,k) = beta_free
3348	endif
3349	if (mskocean_beta) THEN
3350	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3351	endif
3352	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3353	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3354	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3355	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3356	endif
3357	!
3358	ENDDO
3359	ENDDO
3360	!
3361	endif
3362	!
3363	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3364	!
3365	IF (MOD(itaprad,radpas).EQ.0) THEN
3366
3367	!albedo SB >>>
3368	if(ok_chlorophyll)then
3369	print*,"-- reading chlorophyll"
3370	call readchlorophyll(debut)
3371	endif
3372	!do i=1,klon
3373	!if(chl_con(i)>1.) print*,i,chl_con(i),pctsrf(i,is_ter)
3374	!enddo
3375	!albedo SB <<<
3376
3377	!albedo SB >>>
3378	! DO i = 1, klon
3379	! albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce) &
3380	! + falb1(i,is_lic) * pctsrf(i,is_lic) &
3381	! + falb1(i,is_ter) * pctsrf(i,is_ter) &
3382	! + falb1(i,is_sic) * pctsrf(i,is_sic)
3383	! albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce) &
3384	! + falb2(i,is_lic) * pctsrf(i,is_lic) &
3385	! + falb2(i,is_ter) * pctsrf(i,is_ter) &
3386	! + falb2(i,is_sic) * pctsrf(i,is_sic)
3387	! ENDDO
3388	!albedo SB <<<
3389
3390	if (mydebug) then
3391	call writefield_phy('u_seri',u_seri,llm)
3392	call writefield_phy('v_seri',v_seri,llm)
3393	call writefield_phy('t_seri',t_seri,llm)
3394	call writefield_phy('q_seri',q_seri,llm)
3395	endif
3396
3397	IF (aerosol_couple) THEN
3398	#ifdef INCA
3399	CALL radlwsw_inca &
3400	(kdlon,kflev,dist, rmu0, fract, solaire, &
3401	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3402	wo(:, :, 1), &
3403	cldfrarad, cldemirad, cldtaurad, &
3404	heat,heat0,cool,cool0,radsol,albpla, &
3405	topsw,toplw,solsw,sollw, &
3406	sollwdown, &
3407	topsw0,toplw0,solsw0,sollw0, &
3408	lwdn0, lwdn, lwup0, lwup, &
3409	swdn0, swdn, swup0, swup, &
3410	ok_ade, ok_aie, &
3411	tau_aero, piz_aero, cg_aero, &
3412	topswad_aero, solswad_aero, &
3413	topswad0_aero, solswad0_aero, &
3414	topsw_aero, topsw0_aero, &
3415	solsw_aero, solsw0_aero, &
3416	cldtaupirad, &
3417	topswai_aero, solswai_aero)
3418
3419	#endif
3420	ELSE
3421	!
3422	!IM calcul radiatif pour le cas actuel
3423	!
3424	RCO2 = RCO2_act
3425	RCH4 = RCH4_act
3426	RN2O = RN2O_act
3427	RCFC11 = RCFC11_act
3428	RCFC12 = RCFC12_act
3429	!
3430	IF (prt_level .GE.10) THEN
3431	print *,' ->radlwsw, number 1 '
3432	ENDIF
3433	!
3434	CALL radlwsw &
3435	(dist, rmu0, fract, &
3436	!albedo SB >>>
3437	! paprs, pplay,zxtsol,albsol1, albsol2, &
3438	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3439	!albedo SB <<<
3440	t_seri,q_seri,wo, &
3441	cldfrarad, cldemirad, cldtaurad, &
3442	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3443	flag_aerosol_strat, &
3444	tau_aero, piz_aero, cg_aero, &
3445	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! Rajoute par OB pour RRTM
3446	tau_aero_lw_rrtm, &
3447	cldtaupirad,new_aod, &
3448	zqsat, flwc, fiwc, &
3449	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3450	heat,heat0,cool,cool0,radsol,albpla, &
3451	topsw,toplw,solsw,sollw, &
3452	sollwdown, &
3453	topsw0,toplw0,solsw0,sollw0, &
3454	lwdn0, lwdn, lwup0, lwup, &
3455	swdn0, swdn, swup0, swup, &
3456	topswad_aero, solswad_aero, &
3457	topswai_aero, solswai_aero, &
3458	topswad0_aero, solswad0_aero, &
3459	topsw_aero, topsw0_aero, &
3460	solsw_aero, solsw0_aero, &
3461	topswcf_aero, solswcf_aero, &
3462	!-C. Kleinschmitt for LW diagnostics
3463	toplwad_aero, sollwad_aero,&
3464	toplwai_aero, sollwai_aero, &
3465	toplwad0_aero, sollwad0_aero,&
3466	!-end
3467	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3468	ZSWFT0_i, ZFSDN0, ZFSUP0)
3469
3470	!
3471	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3472	!IM des taux doit etre different du taux actuel
3473	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3474	!
3475	if (ok_4xCO2atm) then
3476	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3477	RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3478	RCFC12_per.NE.RCFC12_act) THEN
3479	!
3480	RCO2 = RCO2_per
3481	RCH4 = RCH4_per
3482	RN2O = RN2O_per
3483	RCFC11 = RCFC11_per
3484	RCFC12 = RCFC12_per
3485	!
3486	IF (prt_level .GE.10) THEN
3487	print *,' ->radlwsw, number 2 '
3488	ENDIF
3489	!
3490	CALL radlwsw &
3491	(dist, rmu0, fract, &
3492	!albedo SB >>>
3493	! paprs, pplay,zxtsol,albsol1, albsol2, &
3494	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3495	!albedo SB <<<
3496	t_seri,q_seri,wo, &
3497	cldfra, cldemi, cldtau, &
3498	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3499	flag_aerosol_strat, &
3500	tau_aero, piz_aero, cg_aero, &
3501	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! Rajoute par OB pour RRTM
3502	tau_aero_lw_rrtm, &
3503	cldtaupi,new_aod, &
3504	zqsat, flwc, fiwc, &
3505	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3506	heatp,heat0p,coolp,cool0p,radsolp,albplap, &
3507	topswp,toplwp,solswp,sollwp, &
3508	sollwdownp, &
3509	topsw0p,toplw0p,solsw0p,sollw0p, &
3510	lwdn0p, lwdnp, lwup0p, lwupp, &
3511	swdn0p, swdnp, swup0p, swupp, &
3512	topswad_aerop, solswad_aerop, &
3513	topswai_aerop, solswai_aerop, &
3514	topswad0_aerop, solswad0_aerop, &
3515	topsw_aerop, topsw0_aerop, &
3516	solsw_aerop, solsw0_aerop, &
3517	topswcf_aerop, solswcf_aerop, &
3518	!-C. Kleinschmitt for LW diagnostics
3519	toplwad_aerop, sollwad_aerop,&
3520	toplwai_aerop, sollwai_aerop, &
3521	toplwad0_aerop, sollwad0_aerop,&
3522	!-end
3523	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3524	ZSWFT0_i, ZFSDN0, ZFSUP0)
3525	endif
3526	endif
3527	!
3528	ENDIF ! aerosol_couple
3529	itaprad = 0
3530	ENDIF ! MOD(itaprad,radpas)
3531	itaprad = itaprad + 1
3532
3533	IF (iflag_radia.eq.0) THEN
3534	IF (prt_level.ge.9) THEN
3535	PRINT *,'--------------------------------------------------'
3536	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3537	PRINT *,'>>>> heat et cool mis a zero '
3538	PRINT *,'--------------------------------------------------'
3539	END IF
3540	heat=0.
3541	cool=0.
3542	sollw=0. ! MPL 01032011
3543	solsw=0.
3544	radsol=0.
3545	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3546	swup0=0.
3547	swdn=0.
3548	swdn0=0.
3549	lwup=0.
3550	lwup0=0.
3551	lwdn=0.
3552	lwdn0=0.
3553	END IF
3554
3555	!
3556	! Ajouter la tendance des rayonnements (tous les pas)
3557	!
3558	d_t_swr(:,:)=heat(:,:)*dtime/RDAY
3559	d_t_lwr(:,:)=-cool(:,:)*dtime/RDAY
3560	d_t_sw0(:,:)=heat0(:,:)*dtime/RDAY
3561	d_t_lw0(:,:)=-cool0(:,:)*dtime/RDAY
3562	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW')
3563	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW')
3564
3565	!
3566	if (mydebug) then
3567	call writefield_phy('u_seri',u_seri,llm)
3568	call writefield_phy('v_seri',v_seri,llm)
3569	call writefield_phy('t_seri',t_seri,llm)
3570	call writefield_phy('q_seri',q_seri,llm)
3571	endif
3572
3573	!IM
3574	IF (ip_ebil_phy.ge.2) THEN
3575	ztit='after rad'
3576	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3577	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3578	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3579	call diagphy(airephy,ztit,ip_ebil_phy &
3580	, topsw, toplw, solsw, sollw, zero_v &
3581	, zero_v, zero_v, zero_v, ztsol &
3582	, d_h_vcol, d_qt, d_ec &
3583	, fs_bound, fq_bound )
3584	END IF
3585	!
3586	!
3587	! Calculer l'hydrologie de la surface
3588	!
3589	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3590	! . agesno, ftsol,fqsurf,fsnow, ruis)
3591	!
3592
3593	!
3594	! Calculer le bilan du sol et la derive de temperature (couplage)
3595	!
3596	DO i = 1, klon
3597	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3598	! a la demande de JLD
3599	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3600	ENDDO
3601	!
3602	!moddeblott(jan95)
3603	! Appeler le programme de parametrisation de l'orographie
3604	! a l'echelle sous-maille:
3605	!
3606	IF (prt_level .GE.10) THEN
3607	print *,' call orography ? ', ok_orodr
3608	ENDIF
3609	!
3610	IF (ok_orodr) THEN
3611	!
3612	! selection des points pour lesquels le shema est actif:
3613	igwd=0
3614	DO i=1,klon
3615	itest(i)=0
3616	! IF ((zstd(i).gt.10.0)) THEN
3617	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3618	itest(i)=1
3619	igwd=igwd+1
3620	idx(igwd)=i
3621	ENDIF
3622	ENDDO
3623	! igwdim=MAX(1,igwd)
3624	!
3625	IF (ok_strato) THEN
3626
3627	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
3628	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3629	igwd,idx,itest, &
3630	t_seri, u_seri, v_seri, &
3631	zulow, zvlow, zustrdr, zvstrdr, &
3632	d_t_oro, d_u_oro, d_v_oro)
3633
3634	ELSE
3635	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
3636	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3637	igwd,idx,itest, &
3638	t_seri, u_seri, v_seri, &
3639	zulow, zvlow, zustrdr, zvstrdr, &
3640	d_t_oro, d_u_oro, d_v_oro)
3641	ENDIF
3642	!
3643	! ajout des tendances
3644	!-----------------------------------------------------------------------------------------
3645	! ajout des tendances de la trainee de l'orographie
3646	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro')
3647	!-----------------------------------------------------------------------------------------
3648	!
3649	ENDIF ! fin de test sur ok_orodr
3650	!
3651	if (mydebug) then
3652	call writefield_phy('u_seri',u_seri,llm)
3653	call writefield_phy('v_seri',v_seri,llm)
3654	call writefield_phy('t_seri',t_seri,llm)
3655	call writefield_phy('q_seri',q_seri,llm)
3656	endif
3657
3658	IF (ok_orolf) THEN
3659	!
3660	! selection des points pour lesquels le shema est actif:
3661	igwd=0
3662	DO i=1,klon
3663	itest(i)=0
3664	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3665	itest(i)=1
3666	igwd=igwd+1
3667	idx(igwd)=i
3668	ENDIF
3669	ENDDO
3670	! igwdim=MAX(1,igwd)
3671	!
3672	IF (ok_strato) THEN
3673
3674	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
3675	rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3676	igwd,idx,itest, &
3677	t_seri, u_seri, v_seri, &
3678	zulow, zvlow, zustrli, zvstrli, &
3679	d_t_lif, d_u_lif, d_v_lif )
3680
3681	ELSE
3682	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
3683	rlat,zmea,zstd,zpic, &
3684	itest, &
3685	t_seri, u_seri, v_seri, &
3686	zulow, zvlow, zustrli, zvstrli, &
3687	d_t_lif, d_u_lif, d_v_lif)
3688	ENDIF
3689	!
3690	!-----------------------------------------------------------------------------------------
3691	! ajout des tendances de la portance de l'orographie
3692	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,dqi0,paprs,'lif')
3693	!-----------------------------------------------------------------------------------------
3694	!
3695	ENDIF ! fin de test sur ok_orolf
3696	! HINES GWD PARAMETRIZATION
3697
3698	IF (ok_hines) then
3699
3700	CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
3701	rlat,t_seri,u_seri,v_seri, &
3702	zustrhi,zvstrhi, &
3703	d_t_hin, d_u_hin, d_v_hin)
3704	!
3705	! ajout des tendances
3706	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,dqi0,paprs,'hin')
3707
3708	ENDIF
3709
3710	if (ok_gwd_rando) then
3711	call FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
3712	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
3713	du_gwd_rando, dv_gwd_rando)
3714	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,dqi0,paprs, &
3715	'flott_gwd_rando')
3716	end if
3717
3718	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3719
3720	if (mydebug) then
3721	call writefield_phy('u_seri',u_seri,llm)
3722	call writefield_phy('v_seri',v_seri,llm)
3723	call writefield_phy('t_seri',t_seri,llm)
3724	call writefield_phy('q_seri',q_seri,llm)
3725	endif
3726
3727	DO i = 1, klon
3728	zustrph(i)=0.
3729	zvstrph(i)=0.
3730	ENDDO
3731	DO k = 1, klev
3732	DO i = 1, klon
3733	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
3734	(paprs(i,k)-paprs(i,k+1))/rg
3735	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
3736	(paprs(i,k)-paprs(i,k+1))/rg
3737	ENDDO
3738	ENDDO
3739	!
3740	!IM calcul composantes axiales du moment angulaire et couple des montagnes
3741	!
3742	IF (is_sequential .and. ok_orodr) THEN
3743	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
3744	ra,rg,romega, &
3745	rlat,rlon,pphis, &
3746	zustrdr,zustrli,zustrph, &
3747	zvstrdr,zvstrli,zvstrph, &
3748	paprs,u,v, &
3749	aam, torsfc)
3750	ENDIF
3751	!IM cf. FLott END
3752	!IM
3753	IF (ip_ebil_phy.ge.2) THEN
3754	ztit='after orography'
3755	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3756	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3757	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3758	call diagphy(airephy,ztit,ip_ebil_phy &
3759	, zero_v, zero_v, zero_v, zero_v, zero_v &
3760	, zero_v, zero_v, zero_v, ztsol &
3761	, d_h_vcol, d_qt, d_ec &
3762	, fs_bound, fq_bound )
3763	END IF
3764
3765	!DC Calcul de la tendance due au methane
3766	IF(ok_qch4) THEN
3767	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
3768	! ajout de la tendance d'humidite due au methane
3769	CALL add_phys_tend(du0,dv0,dt0,d_q_ch4*dtime,dql0,'q_ch4')
3770	END IF
3771	!
3772	!
3773	!====================================================================
3774	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3775	!====================================================================
3776	! Abderrahmane 24.08.09
3777
3778	IF (ok_cosp) THEN
3779	! adeclarer
3780	#ifdef CPP_COSP
3781	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3782
3783	IF (prt_level .GE.10) THEN
3784	print*,'freq_cosp',freq_cosp
3785	ENDIF
3786	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3787	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3788	! s ref_liq,ref_ice
3789	call phys_cosp(itap,dtime,freq_cosp, &
3790	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
3791	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
3792	klon,klev,rlon,rlat,presnivs,overlap, &
3793	fract,ref_liq,ref_ice, &
3794	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
3795	zu10m,zv10m,pphis, &
3796	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
3797	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
3798	prfl(:,1:klev),psfl(:,1:klev), &
3799	pmflxr(:,1:klev),pmflxs(:,1:klev), &
3800	mr_ozone,cldtau, cldemi)
3801
3802	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3803	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
3804	! M clMISR,
3805	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3806	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
3807
3808	ENDIF
3809
3810	#endif
3811	ENDIF !ok_cosp
3812	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3813	!AA
3814	!AA Installation de l'interface online-offline pour traceurs
3815	!AA
3816	!====================================================================
3817	! Calcul des tendances traceurs
3818	!====================================================================
3819	!
3820
3821	IF (type_trac=='repr') THEN
3822	sh_in(:,:) = q_seri(:,:)
3823	ELSE
3824	sh_in(:,:) = qx(:,:,ivap)
3825	END IF
3826
3827	call phytrac ( &
3828	itap, days_elapsed+1, jH_cur, debut, &
3829	lafin, dtime, u, v, t, &
3830	paprs, pplay, pmfu, pmfd, &
3831	pen_u, pde_u, pen_d, pde_d, &
3832	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
3833	u1, v1, ftsol, pctsrf, &
3834	zustar, zu10m, zv10m, &
3835	wstar(:,is_ave), ale_bl, ale_wake, &
3836	rlat, rlon, &
3837	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
3838	presnivs, pphis, pphi, albsol1, &
3839	sh_in, rhcl, cldfra, rneb, &
3840	diafra, cldliq, itop_con, ibas_con, &
3841	pmflxr, pmflxs, prfl, psfl, &
3842	da, phi, mp, upwd, &
3843	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
3844	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
3845	ev, ep, epmlmMm, eplaMm, & !<<RomP
3846	dnwd, aerosol_couple, flxmass_w, &
3847	tau_aero, piz_aero, cg_aero, ccm, &
3848	rfname, &
3849	d_tr_dyn, & !<<RomP
3850	tr_seri)
3851
3852	IF (offline) THEN
3853
3854	IF (prt_level.ge.9) &
3855	print*,'Attention on met a 0 les thermiques pour phystoke'
3856	call phystokenc ( &
3857	nlon,klev,pdtphys,rlon,rlat, &
3858	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
3859	fm_therm,entr_therm, &
3860	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
3861	frac_impa, frac_nucl, &
3862	pphis,airephy,dtime,itap, &
3863	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
3864
3865
3866	ENDIF
3867
3868	!
3869	! Calculer le transport de l'eau et de l'energie (diagnostique)
3870	!
3871	CALL transp (paprs,zxtsol, &
3872	t_seri, q_seri, u_seri, v_seri, zphi, &
3873	ve, vq, ue, uq)
3874	!
3875	!IM global posePB BEG
3876	IF(1.EQ.0) THEN
3877	!
3878	CALL transp_lay (paprs,zxtsol, &
3879	t_seri, q_seri, u_seri, v_seri, zphi, &
3880	ve_lay, vq_lay, ue_lay, uq_lay)
3881	!
3882	ENDIF !(1.EQ.0) THEN
3883	!IM global posePB END
3884	! Accumuler les variables a stocker dans les fichiers histoire:
3885	!
3886
3887	!================================================================
3888	! Conversion of kinetic and potential energy into heat, for
3889	! parameterisation of subgrid-scale motions
3890	!================================================================
3891
3892	d_t_ec(:,:)=0.
3893	forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
3894	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
3895	u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
3896	zmasse,exner,d_t_ec)
3897	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
3898
3899	!IM
3900	IF (ip_ebil_phy.ge.1) THEN
3901	ztit='after physic'
3902	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
3903	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3904	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3905	! Comme les tendances de la physique sont ajoute dans la dynamique,
3906	! on devrait avoir que la variation d'entalpie par la dynamique
3907	! est egale a la variation de la physique au pas de temps precedent.
3908	! Donc la somme de ces 2 variations devrait etre nulle.
3909
3910	call diagphy(airephy,ztit,ip_ebil_phy &
3911	, topsw, toplw, solsw, sollw, sens &
3912	, evap, rain_fall, snow_fall, ztsol &
3913	, d_h_vcol, d_qt, d_ec &
3914	, fs_bound, fq_bound )
3915	!
3916	d_h_vcol_phy=d_h_vcol
3917	!
3918	END IF
3919	!
3920	!=======================================================================
3921	! SORTIES
3922	!=======================================================================
3923	!
3924	!IM initialisation + calculs divers diag AMIP2
3925	!
3926	include "calcul_divers.h"
3927	!
3928	!IM Interpolation sur les niveaux de pression du NMC
3929	! -------------------------------------------------
3930	!
3931	include "calcul_STDlev.h"
3932	!
3933	! slp sea level pressure
3934	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3935	!
3936	!cc prw = eau precipitable
3937	DO i = 1, klon
3938	prw(i) = 0.
3939	DO k = 1, klev
3940	prw(i) = prw(i) + &
3941	q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3942	ENDDO
3943	ENDDO
3944	!
3945	IF (type_trac == 'inca') THEN
3946	#ifdef INCA
3947	CALL VTe(VTphysiq)
3948	CALL VTb(VTinca)
3949
3950	CALL chemhook_end ( &
3951	dtime, &
3952	pplay, &
3953	t_seri, &
3954	tr_seri, &
3955	nbtr, &
3956	paprs, &
3957	q_seri, &
3958	airephy, &
3959	pphi, &
3960	pphis, &
3961	zx_rh)
3962
3963	CALL VTe(VTinca)
3964	CALL VTb(VTphysiq)
3965	#endif
3966	END IF
3967
3968
3969	!
3970	! Convertir les incrementations en tendances
3971	!
3972	IF (prt_level .GE.10) THEN
3973	print *,'Convertir les incrementations en tendances '
3974	ENDIF
3975	!
3976	if (mydebug) then
3977	call writefield_phy('u_seri',u_seri,llm)
3978	call writefield_phy('v_seri',v_seri,llm)
3979	call writefield_phy('t_seri',t_seri,llm)
3980	call writefield_phy('q_seri',q_seri,llm)
3981	endif
3982
3983	DO k = 1, klev
3984	DO i = 1, klon
3985	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3986	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3987	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3988	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3989	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3990	!CR: on ajoute le contenu en glace
3991	if (nqo.eq.3) then
3992	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
3993	endif
3994	ENDDO
3995	ENDDO
3996	!
3997	!CR: nb de traceurs eau: nqo
3998	! IF (nqtot.GE.3) THEN
3999	IF (nqtot.GE.(nqo+1)) THEN
4000	! DO iq = 3, nqtot
4001	DO iq = nqo+1, nqtot
4002	DO k = 1, klev
4003	DO i = 1, klon
4004	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
4005	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
4006	ENDDO
4007	ENDDO
4008	ENDDO
4009	ENDIF
4010	!
4011	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4012	!IM global posePB include "write_bilKP_ins.h"
4013	!IM global posePB include "write_bilKP_ave.h"
4014	!
4015
4016	! Sauvegarder les valeurs de t et q a la fin de la physique:
4017	!
4018	DO k = 1, klev
4019	DO i = 1, klon
4020	u_ancien(i,k) = u_seri(i,k)
4021	v_ancien(i,k) = v_seri(i,k)
4022	t_ancien(i,k) = t_seri(i,k)
4023	q_ancien(i,k) = q_seri(i,k)
4024	ENDDO
4025	ENDDO
4026
4027	!!! RomP >>>
4028	!CR: nb de traceurs eau: nqo
4029	! IF (nqtot.GE.3) THEN
4030	IF (nqtot.GE.(nqo+1)) THEN
4031	! DO iq = 3, nqtot
4032	DO iq = nqo+1, nqtot
4033	DO k = 1, klev
4034	DO i = 1, klon
4035	! tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
4036	tr_ancien(i,k,iq-nqo) = tr_seri(i,k,iq-nqo)
4037	ENDDO
4038	ENDDO
4039	ENDDO
4040	ENDIF
4041	!!! RomP <<<
4042	!==========================================================================
4043	! Sorties des tendances pour un point particulier
4044	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4045	! pour le debug
4046	! La valeur de igout est attribuee plus haut dans le programme
4047	!==========================================================================
4048
4049	if (prt_level.ge.1) then
4050	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4051	write(lunout,*) &
4052	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4053	write(lunout,*) &
4054	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4055	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4056	pctsrf(igout,is_sic)
4057	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4058	do k=1,klev
4059	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4060	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4061	d_t_eva(igout,k)
4062	enddo
4063	write(lunout,*) 'cool,heat'
4064	do k=1,klev
4065	write(lunout,*) cool(igout,k),heat(igout,k)
4066	enddo
4067
4068	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4069	do k=1,klev
4070	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4071	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4072	enddo
4073
4074	write(lunout,*) 'd_ps ',d_ps(igout)
4075	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4076	do k=1,klev
4077	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4078	d_qx(igout,k,1),d_qx(igout,k,2)
4079	enddo
4080	endif
4081
4082	!==========================================================================
4083
4084	!============================================================
4085	! Calcul de la temperature potentielle
4086	!============================================================
4087	DO k = 1, klev
4088	DO i = 1, klon
4089	!JYG/IM theta en debut du pas de temps
4090	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4091	!JYG/IM theta en fin de pas de temps de physique
4092	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4093	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers MPL 20130625
4094	! fth_fonctions.F90 et parkind1.F90
4095	! sinon thetal=theta
4096	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4097	! : ql_seri(i,k))
4098	thetal(i,k)=theta(i,k)
4099	ENDDO
4100	ENDDO
4101	!
4102
4103	! 22.03.04 BEG
4104	!=============================================================
4105	! Ecriture des sorties
4106	!=============================================================
4107	#ifdef CPP_IOIPSL
4108
4109	! Recupere des varibles calcule dans differents modules
4110	! pour ecriture dans histxxx.nc
4111
4112	! Get some variables from module fonte_neige_mod
4113	CALL fonte_neige_get_vars(pctsrf, &
4114	zxfqcalving, zxfqfonte, zxffonte)
4115
4116
4117
4118
4119	!=============================================================
4120	! Separation entre thermiques et non thermiques dans les sorties
4121	! de fisrtilp
4122	!=============================================================
4123
4124	if (iflag_thermals>=1) then
4125	d_t_lscth=0.
4126	d_t_lscst=0.
4127	d_q_lscth=0.
4128	d_q_lscst=0.
4129	do k=1,klev
4130	do i=1,klon
4131	if (ptconvth(i,k)) then
4132	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4133	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4134	else
4135	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4136	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4137	endif
4138	enddo
4139	enddo
4140
4141	do i=1,klon
4142	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4143	plul_th(i)=prfl(i,1)+psfl(i,1)
4144	enddo
4145	endif
4146
4147
4148	!On effectue les sorties:
4149
4150	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
4151	pplay, lmax_th, aerosol_couple, &
4152	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
4153	ptconv, read_climoz, clevSTD, &
4154	ptconvth, d_t, qx, d_qx, zmasse, &
4155	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4156
4157
4158
4159
4160	include "write_histday_seri.h"
4161
4162	include "write_paramLMDZ_phy.h"
4163
4164	#endif
4165
4166	! 22.03.04 END
4167	!
4168	!====================================================================
4169	! Si c'est la fin, il faut conserver l'etat de redemarrage
4170	!====================================================================
4171	!
4172
4173	! -----------------------------------------------------------------
4174	! WSTATS: Saving statistics
4175	! -----------------------------------------------------------------
4176	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
4177	! which can later be used to make the statistic files of the run:
4178	! "stats") only possible in 3D runs !
4179
4180
4181	IF (callstats) THEN
4182
4183	call wstats(klon,o_psol%name,"Surface pressure","Pa" &
4184	,2,paprs(:,1))
4185	call wstats(klon,o_tsol%name,"Surface temperature","K", &
4186	2,zxtsol)
4187	zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
4188	call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
4189	"kg/(s*m2)",2,zx_tmp_fi2d)
4190	zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
4191	call wstats(klon,o_plul%name,"Large-scale Precip", &
4192	"kg/(s*m2)",2,zx_tmp_fi2d)
4193	zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
4194	call wstats(klon,o_pluc%name,"Convective Precip", &
4195	"kg/(s*m2)",2,zx_tmp_fi2d)
4196	call wstats(klon,o_sols%name,"Solar rad. at surf.", &
4197	"W/m2",2,solsw)
4198	call wstats(klon,o_soll%name,"IR rad. at surf.", &
4199	"W/m2",2,sollw)
4200	zx_tmp_fi2d(:) = topsw(:)-toplw(:)
4201	call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
4202	"W/m2",2,zx_tmp_fi2d)
4203
4204
4205
4206	call wstats(klon,o_temp%name,"Air temperature","K", &
4207	3,t_seri)
4208	call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
4209	3,u_seri)
4210	call wstats(klon,o_vitv%name,"Meridional wind", &
4211	"m.s-1",3,v_seri)
4212	call wstats(klon,o_vitw%name,"Vertical wind", &
4213	"m.s-1",3,omega)
4214	call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
4215	3,q_seri)
4216
4217
4218
4219	IF(lafin) THEN
4220	write (,) "Writing stats..."
4221	call mkstats(ierr)
4222	ENDIF
4223
4224	ENDIF !if callstats
4225
4226	IF (lafin) THEN
4227	itau_phy = itau_phy + itap
4228	CALL phyredem ("restartphy.nc")
4229	! open(97,form="unformatted",file="finbin")
4230	! write(97) u_seri,v_seri,t_seri,q_seri
4231	! close(97)
4232	!$OMP MASTER
4233	if (read_climoz >= 1) then
4234	if (is_mpi_root) then
4235	call nf95_close(ncid_climoz)
4236	end if
4237	deallocate(press_climoz) ! pointer
4238	end if
4239	!$OMP END MASTER
4240	ENDIF
4241
4242	! first=.false.
4243
4244	RETURN
4245	END SUBROUTINE physiq
4246	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
4247	IMPLICIT none
4248	!
4249	! Calculer et imprimer l'eau totale. A utiliser pour verifier
4250	! la conservation de l'eau
4251	!
4252	include "YOMCST.h"
4253	INTEGER klon,klev
4254	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
4255	REAL aire(klon)
4256	REAL qtotal, zx, qcheck
4257	INTEGER i, k
4258	!
4259	zx = 0.0
4260	DO i = 1, klon
4261	zx = zx + aire(i)
4262	ENDDO
4263	qtotal = 0.0
4264	DO k = 1, klev
4265	DO i = 1, klon
4266	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
4267	*(paprs(i,k)-paprs(i,k+1))/RG
4268	ENDDO
4269	ENDDO
4270	!
4271	qcheck = qtotal/zx
4272	!
4273	RETURN
4274	END FUNCTION qcheck
4275	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
4276	IMPLICIT none
4277	!
4278	! Tranformer une variable de la grille physique a
4279	! la grille d'ecriture
4280	!
4281	INTEGER nfield,nlon,iim,jjmp1, jjm
4282	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
4283	!
4284	INTEGER i, n, ig
4285	!
4286	jjm = jjmp1 - 1
4287	DO n = 1, nfield
4288	DO i=1,iim
4289	ecrit(i,n) = fi(1,n)
4290	ecrit(i+jjm*iim,n) = fi(nlon,n)
4291	ENDDO
4292	DO ig = 1, nlon - 2
4293	ecrit(iim+ig,n) = fi(1+ig,n)
4294	ENDDO
4295	ENDDO
4296	RETURN
4297	END SUBROUTINE gr_fi_ecrit

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

Download in other formats:

Original Format