Context Navigation

physiq.F90 @ 2328

Last change on this file since 2328 was 2328, checked in by Laurent Fairhead, 9 years ago
Special option for Dice 1D case with a new value of iflag_radia=2. It allows to force LWUP computation via ground temperature MPL
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: 140.2 KB

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

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

Download in other formats:

Original Format