Context Navigation

physiq.F90 @ 5223

Last change on this file since 5223 was 2303, checked in by jescribano, 9 years ago
Bugs corrections, control vector is now fine mode+coarse mode and seasalt coarse+fine, change in emission scheme parameters, more outputs at 10h30 and 13h30 LT. (Pending correct optical and sedimentation parameters)
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: 141.1 KB

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

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

Download in other formats:

Original Format