Context Navigation

physiq.F90 @ 2194

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

Use of add_phys_tend to add the tendencies from radiation.
In order to control temperature range after radiation.

Utilisation de add_phys_tend pour ajouter les tendances
du rayonnement, de façon à contrôler les températures après
cet ajout.

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

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

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

Download in other formats:

Original Format