Context Navigation

physiq.F90 @ 2318

Last change on this file since 2318 was 2318, checked in by jyg, 9 years ago
Bug fixing concerning the adiabatic adjustment of the wake region (svn2309).
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.8 KB

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

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

Download in other formats:

Original Format