Context Navigation

physiq.F90 @ 2240

Last change on this file since 2240 was 2240, checked in by fhourdin, 9 years ago

Revisite de la formule des flux de surface
(en priorité sur l'océan) en tenant compte des bourrasques de
vent et de la différence entre les hauteurs de rugosités pour
la quantité de mouvement, l'enthalpie et éventuellement l'humidité.

Etape 1 :
Introduction d'un calcul de gustiness dans la physique
gustiness(:)=f_gust_bl * ale_bl + f_gust_wk * ame_wk
Cette variable est passée ensuite jusqu'au fin fond de la couche limite.
L'étape 1 est prête à commettre, ne nécessite pas de nouvelles
variables dans les startphy et assure la convergence numérique.

Introduction of gustiness in the surface flux computation.
Gustiness is computed from as
gustiness(:)=f_gust_bl * ale_bl + f_gust_wk * ame_wk
and pass through pbl_surface down to the routines that compute
surface fluxes.

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: 137.0 KB

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

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

Download in other formats:

Original Format