Context Navigation

physiq.F90 @ 2221

Last change on this file since 2221 was 2221, checked in by Ehouarn Millour, 10 years ago
Some cleanup: remove (unused) clesph0 from dynamics. EM
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.5 KB

Line
1	! $Id: physiq.F90 2221 2015-03-09 06:38:03Z emillour $
2	!#define IO_DEBUG
3
4	SUBROUTINE physiq (nlon,nlev, &
5	debut,lafin,jD_cur, jH_cur,pdtphys, &
6	paprs,pplay,pphi,pphis,presnivs, &
7	u,v,t,qx, &
8	flxmass_w, &
9	d_u, d_v, d_t, d_qx, d_ps &
10	, dudyn)
11
12	USE ioipsl, only: histbeg, histvert, histdef, histend, histsync, &
13	histwrite, ju2ymds, ymds2ju, getin
14	USE comgeomphy
15	USE phys_cal_mod, only: year_len, mth_len, days_elapsed, jh_1jan, year_cur, &
16	mth_cur, phys_cal_update
17	USE write_field_phy
18	USE dimphy
19	USE infotrac
20	USE mod_grid_phy_lmdz
21	USE mod_phys_lmdz_para
22	USE iophy
23	USE misc_mod, mydebug=>debug
24	USE vampir
25	USE pbl_surface_mod, ONLY : pbl_surface
26	USE change_srf_frac_mod
27	USE surface_data, ONLY : type_ocean, ok_veget, ok_snow
28	USE phys_local_var_mod ! Variables internes non sauvegardees de la physique
29	USE phys_state_var_mod ! Variables sauvegardees de la physique
30	USE phys_output_var_mod ! Variables pour les ecritures des sorties
31	USE phys_output_write_mod
32	USE fonte_neige_mod, ONLY : fonte_neige_get_vars
33	USE phys_output_mod
34	USE phys_output_ctrlout_mod
35	USE iophy
36	use open_climoz_m, only: open_climoz ! ozone climatology from a file
37	use regr_pr_av_m, only: regr_pr_av
38	use netcdf95, only: nf95_close
39	!IM for NMC files
40	! use netcdf, only: nf90_fill_real
41	use netcdf
42	use mod_phys_lmdz_mpi_data, only: is_mpi_root
43	USE aero_mod
44	use ozonecm_m, only: ozonecm ! ozone of J.-F. Royer
45	use conf_phys_m, only: conf_phys
46	use radlwsw_m, only: radlwsw
47	use phyaqua_mod, only: zenang_an
48	USE control_mod
49	#ifdef REPROBUS
50	USE CHEM_REP, ONLY : Init_chem_rep_xjour
51	#endif
52	USE indice_sol_mod
53	USE phytrac_mod, ONLY : phytrac
54
55	#ifdef CPP_RRTM
56	USE YOERAD , ONLY : NRADLP
57	#endif
58
59	!IM stations CFMIP
60	USE CFMIP_point_locations
61	use FLOTT_GWD_rando_m, only: FLOTT_GWD_rando
62
63	IMPLICIT none
64	!>======================================================================
65	!!
66	!! Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
67	!!
68	!! Objet: Moniteur general de la physique du modele
69	!!AA Modifications quant aux traceurs :
70	!!AA - uniformisation des parametrisations ds phytrac
71	!!AA - stockage des moyennes des champs necessaires
72	!!AA en mode traceur off-line
73	!!======================================================================
74	!! CLEFS CPP POUR LES IO
75	!! =====================
76	#define histNMC
77	!!======================================================================
78	!! modif ( P. Le Van , 12/10/98 )
79	!!
80	!! Arguments:
81	!!
82	!! nlon----input-I-nombre de points horizontaux
83	!! nlev----input-I-nombre de couches verticales, doit etre egale a klev
84	!! debut---input-L-variable logique indiquant le premier passage
85	!! lafin---input-L-variable logique indiquant le dernier passage
86	!! jD_cur -R-jour courant a l'appel de la physique (jour julien)
87	!! jH_cur -R-heure courante a l'appel de la physique (jour julien)
88	!! pdtphys-input-R-pas d'integration pour la physique (seconde)
89	!! paprs---input-R-pression pour chaque inter-couche (en Pa)
90	!! pplay---input-R-pression pour le mileu de chaque couche (en Pa)
91	!! pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
92	!! pphis---input-R-geopotentiel du sol
93	!! presnivs-input_R_pressions approximat. des milieux couches ( en PA)
94	!! u-------input-R-vitesse dans la direction X (de O a E) en m/s
95	!! v-------input-R-vitesse Y (de S a N) en m/s
96	!! t-------input-R-temperature (K)
97	!! qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
98	!! d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
99	!! d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
100	!! flxmass_w -input-R- flux de masse verticale
101	!! d_u-----output-R-tendance physique de "u" (m/s/s)
102	!! d_v-----output-R-tendance physique de "v" (m/s/s)
103	!! d_t-----output-R-tendance physique de "t" (K/s)
104	!! d_qx----output-R-tendance physique de "qx" (kg/kg/s)
105	!! d_ps----output-R-tendance physique de la pression au sol
106	!!======================================================================
107	include "dimensions.h"
108	integer jjmp1
109	parameter (jjmp1=jjm+1-1/jjm)
110	integer iip1
111	parameter (iip1=iim+1)
112
113	include "regdim.h"
114	include "dimsoil.h"
115	include "clesphys.h"
116	include "temps.h"
117	include "iniprint.h"
118	include "thermcell.h"
119	!======================================================================
120	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
121	PARAMETER (ok_cvl=.TRUE.)
122	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
123	PARAMETER (ok_gust=.FALSE.)
124	integer iflag_radia ! active ou non le rayonnement (MPL)
125	save iflag_radia
126	!$OMP THREADPRIVATE(iflag_radia)
127	!======================================================================
128	LOGICAL check ! Verifier la conservation du modele en eau
129	PARAMETER (check=.FALSE.)
130	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
131	PARAMETER (ok_stratus=.FALSE.)
132	!======================================================================
133	REAL amn, amx
134	INTEGER igout
135	!======================================================================
136	! Clef controlant l'activation du cycle diurne:
137	!cc LOGICAL cycle_diurne
138	!cc PARAMETER (cycle_diurne=.FALSE.)
139	!======================================================================
140	! Modele thermique du sol, a activer pour le cycle diurne:
141	!cc LOGICAL soil_model
142	!cc PARAMETER (soil_model=.FALSE.)
143	!======================================================================
144	! Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
145	! le calcul du rayonnement est celle apres la precipitation des nuages.
146	! Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
147	! la condensation et la precipitation. Cette cle augmente les impacts
148	! radiatifs des nuages.
149	!cc LOGICAL new_oliq
150	!cc PARAMETER (new_oliq=.FALSE.)
151	!======================================================================
152	! Clefs controlant deux parametrisations de l'orographie:
153	!c LOGICAL ok_orodr
154	!cc PARAMETER (ok_orodr=.FALSE.)
155	!cc LOGICAL ok_orolf
156	!cc PARAMETER (ok_orolf=.FALSE.)
157	!======================================================================
158	LOGICAL ok_journe ! sortir le fichier journalier
159	save ok_journe
160	!$OMP THREADPRIVATE(ok_journe)
161	!
162	LOGICAL ok_mensuel ! sortir le fichier mensuel
163	save ok_mensuel
164	!$OMP THREADPRIVATE(ok_mensuel)
165	!
166	LOGICAL ok_instan ! sortir le fichier instantane
167	save ok_instan
168	!$OMP THREADPRIVATE(ok_instan)
169	!
170	LOGICAL ok_LES ! sortir le fichier LES
171	save ok_LES
172	!$OMP THREADPRIVATE(ok_LES)
173	!
174	LOGICAL callstats ! sortir le fichier stats
175	save callstats
176	!$OMP THREADPRIVATE(callstats)
177	!
178	LOGICAL ok_region ! sortir le fichier regional
179	PARAMETER (ok_region=.FALSE.)
180	!======================================================================
181	real seuil_inversion
182	save seuil_inversion
183	!$OMP THREADPRIVATE(seuil_inversion)
184	integer iflag_ratqs
185	save iflag_ratqs
186	!$OMP THREADPRIVATE(iflag_ratqs)
187	real facteur
188
189	REAL wmax_th(klon)
190	REAL tau_overturning_th(klon)
191
192	integer lmax_th(klon)
193	integer limbas(klon)
194	real ratqscth(klon,klev)
195	real ratqsdiff(klon,klev)
196	real zqsatth(klon,klev)
197
198	!======================================================================
199	!
200	INTEGER ivap ! indice de traceurs pour vapeur d'eau
201	PARAMETER (ivap=1)
202	INTEGER iliq ! indice de traceurs pour eau liquide
203	PARAMETER (iliq=2)
204	!CR: on ajoute la phase glace
205	INTEGER isol ! indice de traceurs pour eau glace
206	PARAMETER (isol=3)
207	!
208	!
209	! Variables argument:
210	!
211	INTEGER nlon
212	INTEGER nlev
213	REAL, intent(in):: jD_cur, jH_cur
214
215	REAL pdtphys
216	LOGICAL debut, lafin
217	REAL paprs(klon,klev+1)
218	REAL pplay(klon,klev)
219	REAL pphi(klon,klev)
220	REAL pphis(klon)
221	REAL presnivs(klev)
222	REAL znivsig(klev)
223	real pir
224
225	REAL u(klon,klev)
226	REAL v(klon,klev)
227	REAL t(klon,klev),thetal(klon,klev)
228	! thetal: ligne suivante a decommenter si vous avez les fichiers MPL 20130625
229	! fth_fonctions.F90 et parkind1.F90
230	! sinon thetal=theta
231	! REAL fth_thetae,fth_thetav,fth_thetal
232	REAL qx(klon,klev,nqtot)
233	REAL flxmass_w(klon,klev)
234	REAL d_u(klon,klev)
235	REAL d_v(klon,klev)
236	REAL d_t(klon,klev)
237	REAL d_qx(klon,klev,nqtot)
238	REAL d_ps(klon)
239	! Variables pour le transport convectif
240	real da(klon,klev),phi(klon,klev,klev),mp(klon,klev)
241	real wght_cvfd(klon,klev)
242	! Variables pour le lessivage convectif
243	! RomP >>>
244	real phi2(klon,klev,klev)
245	real d1a(klon,klev),dam(klon,klev)
246	real ev(klon,klev),ep(klon,klev)
247	real clw(klon,klev),elij(klon,klev,klev)
248	real epmlmMm(klon,klev,klev),eplaMm(klon,klev)
249	! RomP <<<
250	!IM definition dynamique o_trac dans phys_output_open
251	! type(ctrl_out) :: o_trac(nqtot)
252
253	! variables a une pression donnee
254	!
255	include "declare_STDlev.h"
256	!
257	!
258	include "radopt.h"
259	!
260	!
261
262
263	INTEGER debug
264	INTEGER n
265	!ym INTEGER npoints
266	!ym PARAMETER(npoints=klon)
267	!
268	INTEGER nregISCtot
269	PARAMETER(nregISCtot=1)
270	!
271	! imin_debut, nbpti, jmin_debut, nbptj : parametres pour sorties sur 1 region rectangulaire
272	! y compris pour 1 point
273	! imin_debut : indice minimum de i; nbpti : nombre de points en direction i (longitude)
274	! jmin_debut : indice minimum de j; nbptj : nombre de points en direction j (latitude)
275	INTEGER imin_debut, nbpti
276	INTEGER jmin_debut, nbptj
277	!IM: region='3d' <==> sorties en global
278	CHARACTER*3 region
279	PARAMETER(region='3d')
280	logical ok_hf
281	!
282	save ok_hf
283	!$OMP THREADPRIVATE(ok_hf)
284
285	INTEGER,PARAMETER :: longcles=20
286	REAL,SAVE :: clesphy0(longcles)
287	!$OMP THREADPRIVATE(clesphy0)
288	!
289	! Variables propres a la physique
290	INTEGER itap
291	SAVE itap ! compteur pour la physique
292	!$OMP THREADPRIVATE(itap)
293	!
294	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_cldth
639	save iflag_cldth
640	!$OMP THREADPRIVATE(iflag_cldth)
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_cldth,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_cldth.GT.-1) THEN
1017	abort_message = 'Tiedtke needs iflag_cldth=-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_cldth=', iflag_cldth
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	day_ini, &
1311	start_time, &
1312	itau_phy, &
1313	io_lon, &
1314	io_lat)
1315
1316	CALL VTe(VTinca)
1317	CALL VTb(VTphysiq)
1318	#endif
1319	END IF
1320	!
1321	!
1322	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1323	! Nouvelle initialisation pour le rayonnement RRTM
1324	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1325
1326	call iniradia(klon,klev,paprs(1,1:klev+1))
1327
1328	!$omp single
1329	if (read_climoz >= 1) then
1330	call open_climoz(ncid_climoz, press_climoz)
1331	END IF
1332	!$omp end single
1333	!
1334	!IM betaCRF
1335	pfree=70000. !Pa
1336	beta_pbl=1.
1337	beta_free=1.
1338	lon1_beta=-180.
1339	lon2_beta=+180.
1340	lat1_beta=90.
1341	lat2_beta=-90.
1342	mskocean_beta=.FALSE.
1343
1344	OPEN(99,file='beta_crf.data',status='old', &
1345	form='formatted',err=9999)
1346	READ(99,*,end=9998) pfree
1347	READ(99,*,end=9998) beta_pbl
1348	READ(99,*,end=9998) beta_free
1349	READ(99,*,end=9998) lon1_beta
1350	READ(99,*,end=9998) lon2_beta
1351	READ(99,*,end=9998) lat1_beta
1352	READ(99,*,end=9998) lat2_beta
1353	READ(99,*,end=9998) mskocean_beta
1354	9998 Continue
1355	CLOSE(99)
1356	9999 Continue
1357	WRITE(,)'pfree=',pfree
1358	WRITE(,)'beta_pbl=',beta_pbl
1359	WRITE(,)'beta_free=',beta_free
1360	WRITE(,)'lon1_beta=',lon1_beta
1361	WRITE(,)'lon2_beta=',lon2_beta
1362	WRITE(,)'lat1_beta=',lat1_beta
1363	WRITE(,)'lat2_beta=',lat2_beta
1364	WRITE(,)'mskocean_beta=',mskocean_beta
1365	ENDIF
1366	!
1367	! ************** Fin de IF ( debut ) *************
1368	!
1369	!
1370	! Incrementer le compteur de la physique
1371	!
1372	itap = itap + 1
1373	!
1374	!
1375	! Update fraction of the sub-surfaces (pctsrf) and
1376	! initialize, where a new fraction has appeared, all variables depending
1377	! on the surface fraction.
1378	!
1379	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
1380	pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
1381
1382
1383	! Update time and other variables in Reprobus
1384	IF (type_trac == 'repr') THEN
1385	#ifdef REPROBUS
1386	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1387	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1388	CALL Rtime(debut)
1389	#endif
1390	END IF
1391
1392
1393	! Tendances bidons pour les processus qui n'affectent pas certaines
1394	! variables.
1395	du0(:,:)=0.
1396	dv0(:,:)=0.
1397	dt0 = 0.
1398	dq0(:,:)=0.
1399	dql0(:,:)=0.
1400	dqi0(:,:)=0.
1401	!
1402	! Mettre a zero des variables de sortie (pour securite)
1403	!
1404	DO i = 1, klon
1405	d_ps(i) = 0.0
1406	ENDDO
1407	DO k = 1, klev
1408	DO i = 1, klon
1409	d_t(i,k) = 0.0
1410	d_u(i,k) = 0.0
1411	d_v(i,k) = 0.0
1412	ENDDO
1413	ENDDO
1414	DO iq = 1, nqtot
1415	DO k = 1, klev
1416	DO i = 1, klon
1417	d_qx(i,k,iq) = 0.0
1418	ENDDO
1419	ENDDO
1420	ENDDO
1421	da(:,:)=0.
1422	mp(:,:)=0.
1423	phi(:,:,:)=0.
1424	! RomP >>>
1425	phi2(:,:,:)=0.
1426	beta_prec_fisrt(:,:)=0.
1427	beta_prec(:,:)=0.
1428	epmlmMm(:,:,:)=0.
1429	eplaMm(:,:)=0.
1430	d1a(:,:)=0.
1431	dam(:,:)=0.
1432	pmflxr=0.
1433	pmflxs=0.
1434	! RomP <<<
1435
1436	!
1437	! Ne pas affecter les valeurs entrees de u, v, h, et q
1438	!
1439	DO k = 1, klev
1440	DO i = 1, klon
1441	t_seri(i,k) = t(i,k)
1442	u_seri(i,k) = u(i,k)
1443	v_seri(i,k) = v(i,k)
1444	q_seri(i,k) = qx(i,k,ivap)
1445	ql_seri(i,k) = qx(i,k,iliq)
1446	!CR: ATTENTION, on rajoute la variable glace
1447	if (nqo.eq.2) then
1448	qs_seri(i,k) = 0.
1449	else if (nqo.eq.3) then
1450	qs_seri(i,k) = qx(i,k,isol)
1451	endif
1452	ENDDO
1453	ENDDO
1454	tke0(:,:)=pbl_tke(:,:,is_ave)
1455	!CR:Nombre de traceurs de l'eau: nqo
1456	! IF (nqtot.GE.3) THEN
1457	IF (nqtot.GE.(nqo+1)) THEN
1458	! DO iq = 3, nqtot
1459	DO iq = nqo+1, nqtot
1460	DO k = 1, klev
1461	DO i = 1, klon
1462	! tr_seri(i,k,iq-2) = qx(i,k,iq)
1463	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
1464	ENDDO
1465	ENDDO
1466	ENDDO
1467	ELSE
1468	DO k = 1, klev
1469	DO i = 1, klon
1470	tr_seri(i,k,1) = 0.0
1471	ENDDO
1472	ENDDO
1473	ENDIF
1474	!
1475	DO i = 1, klon
1476	ztsol(i) = 0.
1477	ENDDO
1478	DO nsrf = 1, nbsrf
1479	DO i = 1, klon
1480	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1481	ENDDO
1482	ENDDO
1483	!IM
1484	IF (ip_ebil_phy.ge.1) THEN
1485	ztit='after dynamic'
1486	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
1487	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1488	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1489	! Comme les tendances de la physique sont ajoute dans la dynamique,
1490	! on devrait avoir que la variation d'entalpie par la dynamique
1491	! est egale a la variation de la physique au pas de temps precedent.
1492	! Donc la somme de ces 2 variations devrait etre nulle.
1493	call diagphy(airephy,ztit,ip_ebil_phy &
1494	, zero_v, zero_v, zero_v, zero_v, zero_v &
1495	, zero_v, zero_v, zero_v, ztsol &
1496	, d_h_vcol+d_h_vcol_phy, d_qt, 0. &
1497	, fs_bound, fq_bound )
1498	END IF
1499
1500	! Diagnostiquer la tendance dynamique
1501	!
1502	IF (ancien_ok) THEN
1503	DO k = 1, klev
1504	DO i = 1, klon
1505	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1506	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1507	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1508	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1509	ENDDO
1510	ENDDO
1511	!!! RomP >>> td dyn traceur
1512	IF (nqtot.GE.3) THEN
1513	DO iq = 3, nqtot
1514	DO k = 1, klev
1515	DO i = 1, klon
1516	d_tr_dyn(i,k,iq-2)= &
1517	(tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
1518	! iiq=niadv(iq)
1519	! print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
1520	ENDDO
1521	ENDDO
1522	ENDDO
1523	ENDIF
1524	!!! RomP <<<
1525	ELSE
1526	DO k = 1, klev
1527	DO i = 1, klon
1528	d_u_dyn(i,k) = 0.0
1529	d_v_dyn(i,k) = 0.0
1530	d_t_dyn(i,k) = 0.0
1531	d_q_dyn(i,k) = 0.0
1532	ENDDO
1533	ENDDO
1534	!!! RomP >>> td dyn traceur
1535	IF (nqtot.GE.3) THEN
1536	DO iq = 3, nqtot
1537	DO k = 1, klev
1538	DO i = 1, klon
1539	d_tr_dyn(i,k,iq-2)= 0.0
1540	ENDDO
1541	ENDDO
1542	ENDDO
1543	ENDIF
1544	!!! RomP <<<
1545	ancien_ok = .TRUE.
1546	ENDIF
1547	!
1548	! Ajouter le geopotentiel du sol:
1549	!
1550	DO k = 1, klev
1551	DO i = 1, klon
1552	zphi(i,k) = pphi(i,k) + pphis(i)
1553	ENDDO
1554	ENDDO
1555	!
1556	! Verifier les temperatures
1557	!
1558	!IM BEG
1559	IF (check) THEN
1560	amn=MIN(ftsol(1,is_ter),1000.)
1561	amx=MAX(ftsol(1,is_ter),-1000.)
1562	DO i=2, klon
1563	amn=MIN(ftsol(i,is_ter),amn)
1564	amx=MAX(ftsol(i,is_ter),amx)
1565	ENDDO
1566	!
1567	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1568	ENDIF !(check) THEN
1569	!IM END
1570	!
1571	CALL hgardfou(t_seri,ftsol,'debutphy')
1572	!
1573	!IM BEG
1574	IF (check) THEN
1575	amn=MIN(ftsol(1,is_ter),1000.)
1576	amx=MAX(ftsol(1,is_ter),-1000.)
1577	DO i=2, klon
1578	amn=MIN(ftsol(i,is_ter),amn)
1579	amx=MAX(ftsol(i,is_ter),amx)
1580	ENDDO
1581	!
1582	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1583	ENDIF !(check) THEN
1584	!IM END
1585	!
1586	! Mettre en action les conditions aux limites (albedo, sst, etc.).
1587	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
1588	!
1589	if (read_climoz >= 1) then
1590	! Ozone from a file
1591	! Update required ozone index:
1592	ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
1593	if (ro3i == 361) ro3i = 360
1594	! (This should never occur, except perhaps because of roundup
1595	! error. See documentation.)
1596	if (ro3i /= co3i) then
1597	! Update ozone field:
1598	if (read_climoz == 1) then
1599	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
1600	press_in_edg=press_climoz, paprs=paprs, v3=wo)
1601	else
1602	! read_climoz == 2
1603	call regr_pr_av(ncid_climoz, (/"tro3 ", "tro3_daylight"/), &
1604	julien=ro3i, press_in_edg=press_climoz, paprs=paprs, v3=wo)
1605	end if
1606	! Convert from mole fraction of ozone to column density of ozone in a
1607	! cell, in kDU:
1608	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
1609	* zmasse / dobson_u / 1e3
1610	! (By regridding ozone values for LMDZ only once every 360th of
1611	! year, we have already neglected the variation of pressure in one
1612	! 360th of year. So do not recompute "wo" at each time step even if
1613	! "zmasse" changes a little.)
1614	co3i = ro3i
1615	end if
1616	ELSEIF (MOD(itap-1,lmt_pas) == 0) THEN
1617	! Once per day, update ozone from Royer:
1618
1619	IF (solarlong0<-999.) then
1620	! Generic case with evolvoing season
1621	zzz=real(days_elapsed+1)
1622	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
1623	! Particular case with annual mean insolation
1624	zzz=real(90) ! could be revisited
1625	IF (read_climoz/=-1) THEN
1626	abort_message ='read_climoz=-1 is recommended when solarlong0=1000.'
1627	CALL abort_gcm (modname,abort_message,1)
1628	ENDIF
1629	ELSE
1630	! Case where the season is imposed with solarlong0
1631	zzz=real(90) ! could be revisited
1632	ENDIF
1633	wo(:,:,1)=ozonecm(rlat, paprs,read_climoz,rjour=zzz)
1634	ENDIF
1635	!
1636	! Re-evaporer l'eau liquide nuageuse
1637	!
1638	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1639	DO i = 1, klon
1640	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1641	!jyg<
1642	! Attention : Arnaud a propose des formules completement differentes
1643	! A verifier !!!
1644	zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1645	IF (iflag_ice_thermo .EQ. 0) THEN
1646	zlsdcp=zlvdcp
1647	ENDIF
1648	!>jyg
1649
1650	if (iflag_ice_thermo.eq.0) then
1651	!pas necessaire a priori
1652
1653	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1654	zb = MAX(0.0,ql_seri(i,k))
1655	za = - MAX(0.0,ql_seri(i,k)) &
1656	* (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1657	t_seri(i,k) = t_seri(i,k) + za
1658	q_seri(i,k) = q_seri(i,k) + zb
1659	ql_seri(i,k) = 0.0
1660	d_t_eva(i,k) = za
1661	d_q_eva(i,k) = zb
1662
1663	else
1664
1665	!CR: on rÃ©-Ã©vapore eau liquide et glace
1666
1667	! zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1668	! zb = MAX(0.0,ql_seri(i,k))
1669	! za = - MAX(0.0,ql_seri(i,k)) &
1670	! * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1671	zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
1672	za = - MAX(0.0,ql_seri(i,k))*zlvdcp &
1673	- MAX(0.0,qs_seri(i,k))*zlsdcp
1674	t_seri(i,k) = t_seri(i,k) + za
1675	q_seri(i,k) = q_seri(i,k) + zb
1676	ql_seri(i,k) = 0.0
1677	!on Ã©vapore la glace
1678	qs_seri(i,k) = 0.0
1679	d_t_eva(i,k) = za
1680	d_q_eva(i,k) = zb
1681	endif
1682
1683	ENDDO
1684	ENDDO
1685	!IM
1686	IF (ip_ebil_phy.ge.2) THEN
1687	ztit='after reevap'
1688	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
1689	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1690	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1691	call diagphy(airephy,ztit,ip_ebil_phy &
1692	, zero_v, zero_v, zero_v, zero_v, zero_v &
1693	, zero_v, zero_v, zero_v, ztsol &
1694	, d_h_vcol, d_qt, d_ec &
1695	, fs_bound, fq_bound )
1696	!
1697	END IF
1698
1699	!
1700	!=========================================================================
1701	! Calculs de l'orbite.
1702	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1703	! doit donc etre plac\'e avant radlwsw et pbl_surface
1704
1705	!!! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1706	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1707	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1708	!
1709	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1710	! solarlong0
1711	if (solarlong0<-999.) then
1712	if (new_orbit) then
1713	! calcul selon la routine utilisee pour les planetes
1714	call solarlong(day_since_equinox, zlongi, dist)
1715	else
1716	! calcul selon la routine utilisee pour l'AR4
1717	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1718	endif
1719	else
1720	zlongi=solarlong0 ! longitude solaire vraie
1721	dist=1. ! distance au soleil / moyenne
1722	endif
1723	if(prt_level.ge.1) &
1724	write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1725
1726
1727	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1728	! Calcul de l'ensoleillement :
1729	! ============================
1730	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
1731	! l'annee a partir d'une formule analytique.
1732	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
1733	! non nul aux poles.
1734	IF (abs(solarlong0-1000.)<1.e-4) then
1735	call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
1736	ELSE
1737	! Avec ou sans cycle diurne
1738	IF (cycle_diurne) THEN
1739	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1740	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1741	ELSE
1742	CALL angle(zlongi, rlat, fract, rmu0)
1743	ENDIF
1744	ENDIF
1745
1746	! AI Janv 2014
1747	do i = 1, klon
1748	if (fract(i).le.0.) then
1749	JrNt(i)=0.
1750	else
1751	JrNt(i)=1.
1752	endif
1753	enddo
1754
1755	if (mydebug) then
1756	call writefield_phy('u_seri',u_seri,llm)
1757	call writefield_phy('v_seri',v_seri,llm)
1758	call writefield_phy('t_seri',t_seri,llm)
1759	call writefield_phy('q_seri',q_seri,llm)
1760	endif
1761
1762	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1763	! Appel au pbl_surface : Planetary Boudary Layer et Surface
1764	! Cela implique tous les interactions des sous-surfaces et la partie diffusion
1765	! turbulent du couche limit.
1766	!
1767	! Certains varibales de sorties de pbl_surface sont utiliser que pour
1768	! ecriture des fihiers hist_XXXX.nc, ces sont :
1769	! qsol, zq2m, s_pblh, s_lcl,
1770	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1771	! s_therm, s_trmb1, s_trmb2, s_trmb3,
1772	! zxrugs, zu10m, zv10m, fder,
1773	! zxqsurf, rh2m, zxfluxu, zxfluxv,
1774	! frugs, agesno, fsollw, fsolsw,
1775	! d_ts, fevap, fluxlat, t2m,
1776	! wfbils, wfbilo, fluxt, fluxu, fluxv,
1777	!
1778	! Certains ne sont pas utiliser du tout :
1779	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1780	!
1781
1782	! Calcul de l'humidite de saturation au niveau du sol
1783
1784
1785
1786	if (iflag_pbl/=0) then
1787
1788	!jyg+nrlmd<
1789	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
1790	print *,'debut du splitting de la PBL'
1791	ENDIF
1792	!!!
1793	!=================================================================
1794	! PROVISOIRE : DECOUPLAGE PBL/WAKE
1795	! --------------------------------
1796	!
1797	!! wake_deltat_sav(:,:)=wake_deltat(:,:)
1798	!! wake_deltaq_sav(:,:)=wake_deltaq(:,:)
1799	!! wake_deltat(:,:)=0.
1800	!! wake_deltaq(:,:)=0.
1801	!=================================================================
1802	!>jyg+nrlmd
1803	!
1804	CALL pbl_surface( &
1805	dtime, date0, itap, days_elapsed+1, &
1806	debut, lafin, &
1807	rlon, rlat, rugoro, rmu0, &
1808	zsig, sollwdown, pphi, cldt, &
1809	rain_fall, snow_fall, solsw, sollw, &
1810	t_seri, q_seri, u_seri, v_seri, &
1811	!nrlmd+jyg<
1812	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
1813	!>nrlmd+jyg
1814	pplay, paprs, pctsrf, &
1815	ftsol,falb1,falb2,ustar,u10m,v10m,wstar, &
1816	cdragh, cdragm, u1, v1, &
1817	albsol1, albsol2, sens, evap, &
1818	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
1819	zxtsol, zxfluxlat, zt2m, qsat2m, &
1820	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
1821	!nrlmd<
1822	!jyg<
1823	d_t_vdf_w, d_q_vdf_w, &
1824	d_t_vdf_x, d_q_vdf_x, &
1825	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
1826	!>jyg
1827	delta_tsurf,wake_dens, &
1828	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
1829	kh,kh_x,kh_w, &
1830	!>nrlmd
1831	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
1832	slab_wfbils, &
1833	qsol, zq2m, s_pblh, s_lcl, &
1834	!jyg<
1835	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
1836	!>jyg
1837	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
1838	s_therm, s_trmb1, s_trmb2, s_trmb3, &
1839	zxrugs, zustar, zu10m, zv10m, fder, &
1840	zxqsurf, rh2m, zxfluxu, zxfluxv, &
1841	frugs, agesno, fsollw, fsolsw, &
1842	d_ts, fevap, fluxlat, t2m, &
1843	wfbils, wfbilo, fluxt, fluxu, fluxv, &
1844	dsens, devap, zxsnow, &
1845	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
1846	!nrlmd+jyg<
1847	wake_delta_pbl_TKE &
1848	!>nrlmd+jyg
1849	)
1850	!
1851	!=================================================================
1852	! PROVISOIRE : DECOUPLAGE PBL/WAKE
1853	! --------------------------------
1854	!
1855	!! wake_deltat(:,:)=wake_deltat_sav(:,:)
1856	!! wake_deltaq(:,:)=wake_deltaq_sav(:,:)
1857	!=================================================================
1858	!
1859	! Add turbulent diffusion tendency to the wake difference variables
1860	IF (mod(iflag_pbl_split,2) .NE. 0) THEN
1861	wake_deltat(:,:) = wake_deltat(:,:) + (d_t_vdf_w(:,:)-d_t_vdf_x(:,:))
1862	wake_deltaq(:,:) = wake_deltaq(:,:) + (d_q_vdf_w(:,:)-d_q_vdf_x(:,:))
1863	ENDIF
1864
1865
1866	!---------------------------------------------------------------------
1867	! ajout des tendances de la diffusion turbulente
1868	IF (klon_glo==1) THEN
1869	CALL add_pbl_tend &
1870	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
1871	ELSE
1872	CALL add_phys_tend &
1873	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
1874	ENDIF
1875	!--------------------------------------------------------------------
1876
1877	if (mydebug) then
1878	call writefield_phy('u_seri',u_seri,llm)
1879	call writefield_phy('v_seri',v_seri,llm)
1880	call writefield_phy('t_seri',t_seri,llm)
1881	call writefield_phy('q_seri',q_seri,llm)
1882	endif
1883
1884	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
1885	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
1886
1887
1888	IF (ip_ebil_phy.ge.2) THEN
1889	ztit='after surface_main'
1890	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
1891	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1892	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1893	call diagphy(airephy,ztit,ip_ebil_phy &
1894	, zero_v, zero_v, zero_v, zero_v, sens &
1895	, evap , zero_v, zero_v, ztsol &
1896	, d_h_vcol, d_qt, d_ec &
1897	, fs_bound, fq_bound )
1898	END IF
1899
1900	ENDIF
1901	! =================================================================== c
1902	! Calcul de Qsat
1903
1904	DO k = 1, klev
1905	DO i = 1, klon
1906	zx_t = t_seri(i,k)
1907	IF (thermcep) THEN
1908	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1909	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1910	zx_qs = MIN(0.5,zx_qs)
1911	zcor = 1./(1.-retv*zx_qs)
1912	zx_qs = zx_qs*zcor
1913	ELSE
1914	IF (zx_t.LT.t_coup) THEN
1915	zx_qs = qsats(zx_t)/pplay(i,k)
1916	ELSE
1917	zx_qs = qsatl(zx_t)/pplay(i,k)
1918	ENDIF
1919	ENDIF
1920	zqsat(i,k)=zx_qs
1921	ENDDO
1922	ENDDO
1923
1924	if (prt_level.ge.1) then
1925	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1926	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1927	endif
1928	!
1929	! Appeler la convection (au choix)
1930	!
1931	DO k = 1, klev
1932	DO i = 1, klon
1933	conv_q(i,k) = d_q_dyn(i,k) &
1934	+ d_q_vdf(i,k)/dtime
1935	conv_t(i,k) = d_t_dyn(i,k) &
1936	+ d_t_vdf(i,k)/dtime
1937	ENDDO
1938	ENDDO
1939	IF (check) THEN
1940	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1941	WRITE(lunout,*) "avantcon=", za
1942	ENDIF
1943	zx_ajustq = .FALSE.
1944	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1945	IF (zx_ajustq) THEN
1946	DO i = 1, klon
1947	z_avant(i) = 0.0
1948	ENDDO
1949	DO k = 1, klev
1950	DO i = 1, klon
1951	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
1952	*(paprs(i,k)-paprs(i,k+1))/RG
1953	ENDDO
1954	ENDDO
1955	ENDIF
1956
1957	! Calcule de vitesse verticale a partir de flux de masse verticale
1958	DO k = 1, klev
1959	DO i = 1, klon
1960	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1961	END DO
1962	END DO
1963	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
1964	omega(igout, :)
1965
1966	IF (iflag_con.EQ.1) THEN
1967	abort_message ='reactiver le call conlmd dans physiq.F'
1968	CALL abort_gcm (modname,abort_message,1)
1969	! CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1970	! . d_t_con, d_q_con,
1971	! . rain_con, snow_con, ibas_con, itop_con)
1972	ELSE IF (iflag_con.EQ.2) THEN
1973	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
1974	conv_t, conv_q, -evap, omega, &
1975	d_t_con, d_q_con, rain_con, snow_con, &
1976	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
1977	kcbot, kctop, kdtop, pmflxr, pmflxs)
1978	d_u_con = 0.
1979	d_v_con = 0.
1980
1981	WHERE (rain_con < 0.) rain_con = 0.
1982	WHERE (snow_con < 0.) snow_con = 0.
1983	DO i = 1, klon
1984	ibas_con(i) = klev+1 - kcbot(i)
1985	itop_con(i) = klev+1 - kctop(i)
1986	ENDDO
1987	ELSE IF (iflag_con.GE.3) THEN
1988	! nb of tracers for the KE convection:
1989	! MAF la partie traceurs est faite dans phytrac
1990	! on met ntra=1 pour limiter les appels mais on peut
1991	! supprimer les calculs / ftra.
1992	ntra = 1
1993
1994	!=========================================================================
1995	!ajout pour la parametrisation des poches froides: calcul de
1996	!t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1997	do k=1,klev
1998	do i=1,klon
1999	if (iflag_wake>=1) then
2000	t_wake(i,k) = t_seri(i,k) &
2001	+(1-wake_s(i))*wake_deltat(i,k)
2002	q_wake(i,k) = q_seri(i,k) &
2003	+(1-wake_s(i))*wake_deltaq(i,k)
2004	t_undi(i,k) = t_seri(i,k) &
2005	-wake_s(i)*wake_deltat(i,k)
2006	q_undi(i,k) = q_seri(i,k) &
2007	-wake_s(i)*wake_deltaq(i,k)
2008	else
2009	t_wake(i,k) = t_seri(i,k)
2010	q_wake(i,k) = q_seri(i,k)
2011	t_undi(i,k) = t_seri(i,k)
2012	q_undi(i,k) = q_seri(i,k)
2013	endif
2014	enddo
2015	enddo
2016
2017	! Calcul de l'energie disponible ALE (J/kg) et de la puissance
2018	! disponible ALP (W/m2) pour le soulevement des particules dans
2019	! le modele convectif
2020	!
2021	do i = 1,klon
2022	ALE(i) = 0.
2023	ALP(i) = 0.
2024	enddo
2025	!
2026	!calcul de ale_wake et alp_wake
2027	if (iflag_wake>=1) then
2028	if (itap .le. it_wape_prescr) then
2029	do i = 1,klon
2030	ale_wake(i) = wape_prescr
2031	alp_wake(i) = fip_prescr
2032	enddo
2033	else
2034	do i = 1,klon
2035	!jyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2036	!cc ale_wake(i) = 0.5wake_cstar(i)*2
2037	ale_wake(i) = wake_pe(i)
2038	alp_wake(i) = wake_fip(i)
2039	enddo
2040	endif
2041	else
2042	do i = 1,klon
2043	ale_wake(i) = 0.
2044	alp_wake(i) = 0.
2045	enddo
2046	endif
2047	!combinaison avec ale et alp de couche limite: constantes si pas
2048	!de couplage, valeurs calculees dans le thermique sinon
2049	if (iflag_coupl.eq.0) then
2050	if (debut.and.prt_level.gt.9) &
2051	WRITE(lunout,*)'ALE et ALP imposes'
2052	do i = 1,klon
2053	!on ne couple que ale
2054	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2055	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2056	!on ne couple que alp
2057	! ALP(i) = alp_wake(i) + Alp_bl(i)
2058	ALP(i) = alp_wake(i) + alp_bl_prescr
2059	enddo
2060	else
2061	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2062	! do i = 1,klon
2063	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2064	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
2065	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2066	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2067	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2068	! enddo
2069
2070	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2071	! Modif FH 2010/04/27. Sans doute temporaire.
2072	! Deux options pour le alp_offset : constant si >?? 0 ou
2073	! proportionnel ??a w si <0
2074	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2075	! Estimation d'une vitesse verticale effective pour ALP
2076	www(1:klon)=0.
2077	do k=2,klev-1
2078	do i=1,klon
2079	www(i)=max(www(i),-omega(i,k)RDt_seri(i,k)/(RG*paprs(i,k)) &
2080	& zw2(i,k)zw2(i,k))
2081	! if (paprs(i,k)>pbase(i)) then
2082	! calcul approche de la vitesse verticale en m/s
2083	! www(i)=max(www(i),-omega(i,k)RDtemp(i,k)/(RG*paprs(i,k))
2084	! endif
2085	! Le 0.1 est en gros H / ps = 1e5 / 1e4
2086	enddo
2087	enddo
2088	do i=1,klon
2089	if (www(i)>0. .and. ale_bl(i)>0. ) www(i)=www(i)/ale_bl(i)
2090	enddo
2091
2092
2093	do i = 1,klon
2094	ALE(i) = max(ale_wake(i),Ale_bl(i))
2095	!cc nrlmd le 10/04/2012----------Stochastic triggering--------------
2096	if (iflag_trig_bl.ge.1) then
2097	ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
2098	endif
2099	!cc fin nrlmd le 10/04/2012
2100	if (alp_offset>=0.) then
2101	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2102	else
2103
2104	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2105	! _ _
2106	! Ajout d'une composante 3 * A * w w'2 a w'3 avec w=www : w max sous pbase
2107	! ou A est la fraction couverte par les ascendances w'
2108	! on utilise le fait que A * w'3 = ALP
2109	! et donc A * w'2 ~ ALP / sqrt(ALE) (on ajoute 0.1 pour les
2110	! singularites)
2111	ALP(i)=alp_wake(i)(1.+3.www(i)/( sqrt(ale_wake(i))+0.1) ) &
2112	& +alp_bl(i) (1.+3.www(i)/( sqrt(ale_bl(i)) +0.1) )
2113	! ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2114	! if (alp(i)<0.) then
2115	! print*,'ALP ',alp(i),alp_wake(i) &
2116	! ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2117	! endif
2118	endif
2119	enddo
2120	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2121
2122	endif
2123	do i=1,klon
2124	if (alp(i)>alp_max) then
2125	IF(prt_level>9)WRITE(lunout,*) &
2126	'WARNING SUPER ALP (seuil=',alp_max, &
2127	'): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2128	alp(i)=alp_max
2129	endif
2130	if (ale(i)>ale_max) then
2131	IF(prt_level>9)WRITE(lunout,*) &
2132	'WARNING SUPER ALE (seuil=',ale_max, &
2133	'): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2134	ale(i)=ale_max
2135	endif
2136	enddo
2137
2138	!fin calcul ale et alp
2139	!=======================================================================
2140
2141
2142	! sb, oct02:
2143	! Schema de convection modularise et vectorise:
2144	! (driver commun aux versions 3 et 4)
2145	!
2146	IF (ok_cvl) THEN ! new driver for convectL
2147
2148	IF (type_trac == 'repr') THEN
2149	nbtr_tmp=ntra
2150	ELSE
2151	nbtr_tmp=nbtr
2152	END IF
2153	!jyg iflag_con est dans clesphys
2154	!c CALL concvl (iflag_con,iflag_clos,
2155	CALL concvl (iflag_clos, &
2156	dtime,paprs,pplay,t_undi,q_undi, &
2157	t_wake,q_wake,wake_s, &
2158	u_seri,v_seri,tr_seri,nbtr_tmp, &
2159	ALE,ALP, &
2160	sig1,w01, &
2161	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2162	rain_con, snow_con, ibas_con, itop_con, sigd, &
2163	ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
2164	Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2165	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2166	! RomP >>>
2167	!! . pmflxr,pmflxs,da,phi,mp,
2168	!! . ftd,fqd,lalim_conv,wght_th)
2169	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2170	ftd,fqd,lalim_conv,wght_th, &
2171	ev, ep,epmlmMm,eplaMm, &
2172	wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2173	tau_cld_cv,coefw_cld_cv)
2174	! RomP <<<
2175
2176	!IM begin
2177	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2178	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2179	!IM end
2180	!IM cf. FH
2181	clwcon0=qcondc
2182	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2183
2184	do i = 1, klon
2185	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2186	enddo
2187
2188	ELSE ! ok_cvl
2189
2190	! MAF conema3 ne contient pas les traceurs
2191	CALL conema3 (dtime, &
2192	paprs,pplay,t_seri,q_seri, &
2193	u_seri,v_seri,tr_seri,ntra, &
2194	sig1,w01, &
2195	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2196	rain_con, snow_con, ibas_con, itop_con, &
2197	upwd,dnwd,dnwd0,bas,top, &
2198	Ma,cape,tvp,rflag, &
2199	pbase &
2200	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2201	,clwcon0)
2202
2203	ENDIF ! ok_cvl
2204
2205	!
2206	! Correction precip
2207	rain_con = rain_con * cvl_corr
2208	snow_con = snow_con * cvl_corr
2209	!
2210
2211	IF (.NOT. ok_gust) THEN
2212	do i = 1, klon
2213	wd(i)=0.0
2214	enddo
2215	ENDIF
2216
2217	! =================================================================== c
2218	! Calcul des proprietes des nuages convectifs
2219	!
2220
2221	! calcul des proprietes des nuages convectifs
2222	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2223	IF (iflag_cld_cv <= 1) THEN
2224	call clouds_gno &
2225	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2226	ELSE
2227	call clouds_bigauss &
2228	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
2229	ENDIF
2230
2231
2232	! =================================================================== c
2233
2234	DO i = 1, klon
2235	itop_con(i) = min(max(itop_con(i),1),klev)
2236	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2237	ENDDO
2238
2239	DO i = 1, klon
2240	ema_pcb(i) = paprs(i,ibas_con(i))
2241	ENDDO
2242	DO i = 1, klon
2243	! L'idicage de itop_con peut cacher un pb potentiel
2244	! FH sous la dictee de JYG, CR
2245	ema_pct(i) = paprs(i,itop_con(i)+1)
2246
2247	if (itop_con(i).gt.klev-3) then
2248	if(prt_level >= 9) then
2249	write(lunout,*)'La convection monte trop haut '
2250	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2251	endif
2252	endif
2253	ENDDO
2254	ELSE IF (iflag_con.eq.0) THEN
2255	write(lunout,*) 'On n appelle pas la convection'
2256	clwcon0=0.
2257	rnebcon0=0.
2258	d_t_con=0.
2259	d_q_con=0.
2260	d_u_con=0.
2261	d_v_con=0.
2262	rain_con=0.
2263	snow_con=0.
2264	bas=1
2265	top=1
2266	ELSE
2267	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2268	call abort_gcm("physiq", "", 1)
2269	ENDIF
2270
2271	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2272	! . d_u_con, d_v_con)
2273
2274	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
2275	'convection')
2276	!----------------------------------------------------------------------------
2277
2278	if (mydebug) then
2279	call writefield_phy('u_seri',u_seri,llm)
2280	call writefield_phy('v_seri',v_seri,llm)
2281	call writefield_phy('t_seri',t_seri,llm)
2282	call writefield_phy('q_seri',q_seri,llm)
2283	endif
2284
2285	!IM
2286	IF (ip_ebil_phy.ge.2) THEN
2287	ztit='after convect'
2288	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2289	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2290	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2291	call diagphy(airephy,ztit,ip_ebil_phy &
2292	, zero_v, zero_v, zero_v, zero_v, zero_v &
2293	, zero_v, rain_con, snow_con, ztsol &
2294	, d_h_vcol, d_qt, d_ec &
2295	, fs_bound, fq_bound )
2296	END IF
2297	!
2298	IF (check) THEN
2299	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2300	WRITE(lunout,*)"aprescon=", za
2301	zx_t = 0.0
2302	za = 0.0
2303	DO i = 1, klon
2304	za = za + airephy(i)/REAL(klon)
2305	zx_t = zx_t + (rain_con(i)+ &
2306	snow_con(i))*airephy(i)/REAL(klon)
2307	ENDDO
2308	zx_t = zx_t/za*dtime
2309	WRITE(lunout,*)"Precip=", zx_t
2310	ENDIF
2311	IF (zx_ajustq) THEN
2312	DO i = 1, klon
2313	z_apres(i) = 0.0
2314	ENDDO
2315	DO k = 1, klev
2316	DO i = 1, klon
2317	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2318	*(paprs(i,k)-paprs(i,k+1))/RG
2319	ENDDO
2320	ENDDO
2321	DO i = 1, klon
2322	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2323	/z_apres(i)
2324	ENDDO
2325	DO k = 1, klev
2326	DO i = 1, klon
2327	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2328	z_factor(i).LT.(1.0-1.0E-08)) THEN
2329	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2330	ENDIF
2331	ENDDO
2332	ENDDO
2333	ENDIF
2334	zx_ajustq=.FALSE.
2335
2336	!
2337	!=============================================================================
2338	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
2339	!pour la couche limite diffuse pour l instant
2340	!
2341	!
2342	!!! nrlmd le 22/03/2011---Si on met les poches hors des thermiques il faut rajouter cette
2343	!------------------------- tendance calculée hors des poches froides
2344	!
2345	if (iflag_wake>=1) then
2346	DO k=1,klev
2347	DO i=1,klon
2348	dt_dwn(i,k) = ftd(i,k)
2349	dq_dwn(i,k) = fqd(i,k)
2350	M_dwn(i,k) = dnwd0(i,k)
2351	M_up(i,k) = upwd(i,k)
2352	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2353	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2354	ENDDO
2355	ENDDO
2356	!nrlmd+jyg<
2357	DO k=1,klev
2358	DO i=1,klon
2359	wdt_PBL(i,k) = 0.
2360	wdq_PBL(i,k) = 0.
2361	udt_PBL(i,k) = 0.
2362	udq_PBL(i,k) = 0.
2363	ENDDO
2364	ENDDO
2365	!
2366	IF (mod(iflag_pbl_split,2) .EQ. 1) THEN
2367	DO k=1,klev
2368	DO i=1,klon
2369	wdt_PBL(i,k) = wdt_PBL(i,k) + d_t_vdf_w(i,k)/dtime
2370	wdq_PBL(i,k) = wdq_PBL(i,k) + d_q_vdf_w(i,k)/dtime
2371	udt_PBL(i,k) = udt_PBL(i,k) + d_t_vdf_x(i,k)/dtime
2372	udq_PBL(i,k) = udq_PBL(i,k) + d_q_vdf_x(i,k)/dtime
2373	!! dt_dwn(i,k) = dt_dwn(i,k) + d_t_vdf_w(i,k)/dtime
2374	!! dq_dwn(i,k) = dq_dwn(i,k) + d_q_vdf_w(i,k)/dtime
2375	!! dt_a (i,k) = dt_a(i,k) + d_t_vdf_x(i,k)/dtime
2376	!! dq_a (i,k) = dq_a(i,k) + d_q_vdf_x(i,k)/dtime
2377	ENDDO
2378	ENDDO
2379	ENDIF
2380	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2381	DO k=1,klev
2382	DO i=1,klon
2383	!! dt_dwn(i,k) = dt_dwn(i,k) + 0.
2384	!! dq_dwn(i,k) = dq_dwn(i,k) + 0.
2385	!! dt_a(i,k) = dt_a(i,k) + d_t_ajs(i,k)/dtime
2386	!! dq_a(i,k) = dq_a(i,k) + d_q_ajs(i,k)/dtime
2387	udt_PBL(i,k) = udt_PBL(i,k) + d_t_ajs(i,k)/dtime
2388	udq_PBL(i,k) = udq_PBL(i,k) + d_q_ajs(i,k)/dtime
2389	ENDDO
2390	ENDDO
2391	ENDIF
2392	!>nrlmd+jyg
2393
2394	IF (iflag_wake==2) THEN
2395	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2396	DO k = 1,klev
2397	dt_dwn(:,k)= dt_dwn(:,k)+ &
2398	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2399	dq_dwn(:,k)= dq_dwn(:,k)+ &
2400	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2401	ENDDO
2402	ELSEIF (iflag_wake==3) THEN
2403	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2404	DO k = 1,klev
2405	DO i=1,klon
2406	IF (rneb(i,k)==0.) THEN
2407	! On ne tient compte des tendances qu'en dehors des nuages (c'est ï¿½\| dire
2408	! a priri dans une region ou l'eau se reevapore).
2409	dt_dwn(i,k)= dt_dwn(i,k)+ &
2410	ok_wk_lsp(i)*d_t_lsc(i,k)/dtime
2411	dq_dwn(i,k)= dq_dwn(i,k)+ &
2412	ok_wk_lsp(i)*d_q_lsc(i,k)/dtime
2413	ENDIF
2414	ENDDO
2415	ENDDO
2416	ENDIF
2417
2418	!
2419	!calcul caracteristiques de la poche froide
2420	call calWAKE (paprs,pplay,dtime &
2421	,t_seri,q_seri,omega &
2422	,dt_dwn,dq_dwn,M_dwn,M_up &
2423	,dt_a,dq_a,sigd &
2424	,wdt_PBL,wdq_PBL &
2425	,udt_PBL,udq_PBL &
2426	,wake_deltat,wake_deltaq,wake_dth &
2427	,wake_h,wake_s,wake_dens &
2428	,wake_pe,wake_fip,wake_gfl &
2429	,dt_wake,dq_wake &
2430	,wake_k, t_undi,q_undi &
2431	,wake_omgbdth,wake_dp_omgb &
2432	,wake_dtKE,wake_dqKE &
2433	,wake_dtPBL,wake_dqPBL &
2434	,wake_omg,wake_dp_deltomg &
2435	,wake_spread,wake_Cstar,wake_d_deltat_gw &
2436	,wake_ddeltat,wake_ddeltaq)
2437	!
2438	!-------------------------------------------------------------------------
2439	! ajout des tendances des poches froides
2440	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2441	! coherent avec les autres d_t_...
2442	d_t_wake(:,:)=dt_wake(:,:)*dtime
2443	d_q_wake(:,:)=dq_wake(:,:)*dtime
2444	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake')
2445	!------------------------------------------------------------------------
2446
2447	endif ! (iflag_wake>=1)
2448	!
2449	!===================================================================
2450	!JYG
2451	IF (ip_ebil_phy.ge.2) THEN
2452	ztit='after wake'
2453	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2454	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2455	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2456	call diagphy(airephy,ztit,ip_ebil_phy &
2457	, zero_v, zero_v, zero_v, zero_v, zero_v &
2458	, zero_v, zero_v, zero_v, ztsol &
2459	, d_h_vcol, d_qt, d_ec &
2460	, fs_bound, fq_bound )
2461	END IF
2462
2463	! print*,'apres callwake iflag_cldth=', iflag_cldth
2464	!
2465	!===================================================================
2466	! Convection seche (thermiques ou ajustement)
2467	!===================================================================
2468	!
2469	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
2470	,seuil_inversion,weak_inversion,dthmin)
2471
2472
2473
2474	d_t_ajsb(:,:)=0.
2475	d_q_ajsb(:,:)=0.
2476	d_t_ajs(:,:)=0.
2477	d_u_ajs(:,:)=0.
2478	d_v_ajs(:,:)=0.
2479	d_q_ajs(:,:)=0.
2480	clwcon0th(:,:)=0.
2481	!
2482	! fm_therm(:,:)=0.
2483	! entr_therm(:,:)=0.
2484	! detr_therm(:,:)=0.
2485	!
2486	IF(prt_level>9)WRITE(lunout,*) &
2487	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
2488	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2489	if(iflag_thermals<0) then
2490	! Rien
2491	! ====
2492	IF(prt_level>9)WRITE(lunout,*)'pas de convection seche'
2493
2494
2495	else
2496
2497	! Thermiques
2498	! ==========
2499	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
2500	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2501
2502
2503	!cc nrlmd le 10/04/2012
2504	DO k=1,klev+1
2505	DO i=1,klon
2506	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2507	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2508	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2509	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2510	ENDDO
2511	ENDDO
2512	!cc fin nrlmd le 10/04/2012
2513
2514	if (iflag_thermals>=1) then
2515	!jyg<
2516	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2517	! Appel des thermiques avec les profils exterieurs aux poches
2518	DO k=1,klev
2519	DO i=1,klon
2520	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2521	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2522	ENDDO
2523	ENDDO
2524	ELSE
2525	! Appel des thermiques avec les profils moyens
2526	DO k=1,klev
2527	DO i=1,klon
2528	t_therm(i,k) = t_seri(i,k)
2529	q_therm(i,k) = q_seri(i,k)
2530	ENDDO
2531	ENDDO
2532	ENDIF
2533	!>jyg
2534	call calltherm(pdtphys &
2535	,pplay,paprs,pphi,weak_inversion &
2536	!! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
2537	,u_seri,v_seri,t_therm,q_therm,zqsat,debut & !jyg
2538	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
2539	,fm_therm,entr_therm,detr_therm &
2540	,zqasc,clwcon0th,lmax_th,ratqscth &
2541	,ratqsdiff,zqsatth &
2542	!on rajoute ale et alp, et les caracteristiques de la couche alim
2543	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
2544	,ztv,zpspsk,ztla,zthl &
2545	!cc nrlmd le 10/04/2012
2546	,pbl_tke_input,pctsrf,omega,airephy &
2547	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
2548	,n2,s2,ale_bl_stat &
2549	,therm_tke_max,env_tke_max &
2550	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
2551	,alp_bl_conv,alp_bl_stat &
2552	!cc fin nrlmd le 10/04/2012
2553	,zqla,ztva )
2554	!
2555	!jyg<
2556	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2557	! Si les thermiques ne sont presents que hors des poches, la tendance moyenne
2558	! associée doit etre multipliee par la fraction surfacique qu'ils couvrent.
2559	DO k=1,klev
2560	DO i=1,klon
2561	!
2562	wake_deltat(i,k) = wake_deltat(i,k) - d_t_ajs(i,k)
2563	wake_deltaq(i,k) = wake_deltaq(i,k) - d_q_ajs(i,k)
2564	t_seri(i,k) = t_therm(i,k) + wake_s(i)*wake_deltat(i,k)
2565	q_seri(i,k) = q_therm(i,k) + wake_s(i)*wake_deltaq(i,k)
2566	!
2567	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
2568	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
2569	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
2570	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
2571	!
2572	ENDDO
2573	ENDDO
2574	ELSE
2575	DO k=1,klev
2576	DO i=1,klon
2577	t_seri(i,k) = t_therm(i,k)
2578	q_seri(i,k) = q_therm(i,k)
2579	ENDDO
2580	ENDDO
2581	ENDIF
2582	!>jyg
2583
2584	!cc nrlmd le 10/04/2012
2585	!-----------Stochastic triggering-----------
2586	if (iflag_trig_bl.ge.1) then
2587	!
2588	IF (prt_level .GE. 10) THEN
2589	print *,'cin, ale_bl_stat, alp_bl_stat ', &
2590	cin, ale_bl_stat, alp_bl_stat
2591	ENDIF
2592
2593
2594	!----Initialisations
2595	do i=1,klon
2596	proba_notrig(i)=1.
2597	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
2598	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
2599	tau_trig(i)=tau_trig_shallow
2600	else
2601	tau_trig(i)=tau_trig_deep
2602	endif
2603	enddo
2604	!
2605	IF (prt_level .GE. 10) THEN
2606	print *,'random_notrig, tau_trig ', &
2607	random_notrig, tau_trig
2608	print *,'s_trig,s2,n2 ', &
2609	s_trig,s2,n2
2610	ENDIF
2611
2612	!Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
2613	IF (iflag_trig_bl.eq.1) then
2614
2615	!----Tirage al\'eatoire et calcul de ale_bl_trig
2616	do i=1,klon
2617	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
2618	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2619	(n2(i)*dtime/tau_trig(i))
2620	! print *, 'proba_notrig(i) ',proba_notrig(i)
2621	if (random_notrig(i) .ge. proba_notrig(i)) then
2622	ale_bl_trig(i)=ale_bl_stat(i)
2623	else
2624	ale_bl_trig(i)=0.
2625	endif
2626	else
2627	proba_notrig(i)=1.
2628	random_notrig(i)=0.
2629	ale_bl_trig(i)=0.
2630	endif
2631	enddo
2632
2633	ELSE IF (iflag_trig_bl.ge.2) then
2634
2635	do i=1,klon
2636	if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) ) then
2637	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2638	(n2(i)*dtime/tau_trig(i))
2639	! print *, 'proba_notrig(i) ',proba_notrig(i)
2640	if (random_notrig(i) .ge. proba_notrig(i)) then
2641	ale_bl_trig(i)=Ale_bl(i)
2642	else
2643	ale_bl_trig(i)=0.
2644	endif
2645	else
2646	proba_notrig(i)=1.
2647	random_notrig(i)=0.
2648	ale_bl_trig(i)=0.
2649	endif
2650	enddo
2651
2652	ENDIF
2653
2654	!
2655	IF (prt_level .GE. 10) THEN
2656	print *,'proba_notrig, ale_bl_trig ', &
2657	proba_notrig, ale_bl_trig
2658	ENDIF
2659
2660	endif !(iflag_trig_bl)
2661
2662	!-----------Statistical closure-----------
2663	if (iflag_clos_bl.eq.1) then
2664
2665	do i=1,klon
2666	!CR: alp probabiliste
2667	if (ale_bl_trig(i).gt.0.) then
2668	alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
2669	endif
2670	enddo
2671
2672	else if (iflag_clos_bl.eq.2) then
2673
2674	!CR: alp calculee dans thermcell_main
2675	do i=1,klon
2676	alp_bl(i)=alp_bl_stat(i)
2677	enddo
2678
2679	else
2680
2681	alp_bl_stat(:)=0.
2682
2683	endif !(iflag_clos_bl)
2684
2685	IF (prt_level .GE. 10) THEN
2686	print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
2687	ENDIF
2688
2689	!cc fin nrlmd le 10/04/2012
2690
2691	! ----------------------------------------------------------------------
2692	! Transport de la TKE par les panaches thermiques.
2693	! FH : 2010/02/01
2694	! if (iflag_pbl.eq.10) then
2695	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2696	! s rg,paprs,pbl_tke)
2697	! endif
2698	! ----------------------------------------------------------------------
2699	!IM/FH: 2011/02/23
2700	! Couplage Thermiques/Emanuel seulement si T<0
2701	if (iflag_coupl==2) then
2702	IF (prt_level .GE. 10) THEN
2703	print*,'Couplage Thermiques/Emanuel seulement si T<0'
2704	ENDIF
2705	do i=1,klon
2706	if (t_seri(i,lmax_th(i))>273.) then
2707	Ale_bl(i)=0.
2708	endif
2709	enddo
2710	endif
2711
2712	do i=1,klon
2713	! zmax_th(i)=pphi(i,lmax_th(i))/rg
2714	!CR:04/05/12:correction calcul zmax
2715	zmax_th(i)=zmax0(i)
2716	enddo
2717
2718	endif
2719
2720
2721	! Ajustement sec
2722	! ==============
2723
2724	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
2725	! a partir du sommet des thermiques.
2726	! Dans le cas contraire, on demarre au niveau 1.
2727
2728	if (iflag_thermals>=13.or.iflag_thermals<=0) then
2729
2730	if(iflag_thermals.eq.0) then
2731	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2732	limbas(:)=1
2733	else
2734	limbas(:)=lmax_th(:)
2735	endif
2736
2737	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2738	! pour des test de convergence numerique.
2739	! Le nouveau ajsec est a priori mieux, meme pour le cas
2740	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2741	! non nulles numeriquement pour des mailles non concernees.
2742
2743	if (iflag_thermals==0) then
2744	! Calling adjustment alone (but not the thermal plume model)
2745	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
2746	, d_t_ajsb, d_q_ajsb)
2747	else if (iflag_thermals>0) then
2748	! Calling adjustment above the top of thermal plumes
2749	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
2750	, d_t_ajsb, d_q_ajsb)
2751	endif
2752
2753	!-----------------------------------------------------------------------
2754	! ajout des tendances de l'ajustement sec ou des thermiques
2755	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs,'ajsb')
2756	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2757	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2758
2759	!---------------------------------------------------------------------
2760
2761	endif
2762
2763	endif
2764	!
2765	!===================================================================
2766	!IM
2767	IF (ip_ebil_phy.ge.2) THEN
2768	ztit='after dry_adjust'
2769	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2770	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2771	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2772	call diagphy(airephy,ztit,ip_ebil_phy &
2773	, zero_v, zero_v, zero_v, zero_v, zero_v &
2774	, zero_v, zero_v, zero_v, ztsol &
2775	, d_h_vcol, d_qt, d_ec &
2776	, fs_bound, fq_bound )
2777	END IF
2778
2779
2780	!-------------------------------------------------------------------------
2781	! Computation of ratqs, the width (normalized) of the subrid scale
2782	! water distribution
2783	CALL calcratqs(klon,klev,prt_level,lunout, &
2784	iflag_ratqs,iflag_con,iflag_cldth,pdtphys, &
2785	ratqsbas,ratqshaut,tau_ratqs,fact_cldcon, &
2786	ptconv,ptconvth,clwcon0th, rnebcon0th, &
2787	paprs,pplay,q_seri,zqsat,fm_therm, &
2788	ratqs,ratqsc)
2789
2790
2791	!
2792	! Appeler le processus de condensation a grande echelle
2793	! et le processus de precipitation
2794	!-------------------------------------------------------------------------
2795	IF (prt_level .GE.10) THEN
2796	print *,'itap, ->fisrtilp ',itap
2797	ENDIF
2798	!
2799	CALL fisrtilp(dtime,paprs,pplay, &
2800	t_seri, q_seri,ptconv,ratqs, &
2801	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
2802	rain_lsc, snow_lsc, &
2803	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
2804	frac_impa, frac_nucl, beta_prec_fisrt, &
2805	prfl, psfl, rhcl, &
2806	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldth, &
2807	iflag_ice_thermo)
2808	!
2809	WHERE (rain_lsc < 0) rain_lsc = 0.
2810	WHERE (snow_lsc < 0) snow_lsc = 0.
2811
2812	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs,'lsc')
2813	!---------------------------------------------------------------------------
2814	DO k = 1, klev
2815	DO i = 1, klon
2816	cldfra(i,k) = rneb(i,k)
2817	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
2818	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2819	ENDDO
2820	ENDDO
2821	IF (check) THEN
2822	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2823	WRITE(lunout,*)"apresilp=", za
2824	zx_t = 0.0
2825	za = 0.0
2826	DO i = 1, klon
2827	za = za + airephy(i)/REAL(klon)
2828	zx_t = zx_t + (rain_lsc(i) &
2829	+ snow_lsc(i))*airephy(i)/REAL(klon)
2830	ENDDO
2831	zx_t = zx_t/za*dtime
2832	WRITE(lunout,*)"Precip=", zx_t
2833	ENDIF
2834	!IM
2835	IF (ip_ebil_phy.ge.2) THEN
2836	ztit='after fisrt'
2837	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2838	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2839	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2840	call diagphy(airephy,ztit,ip_ebil_phy &
2841	, zero_v, zero_v, zero_v, zero_v, zero_v &
2842	, zero_v, rain_lsc, snow_lsc, ztsol &
2843	, d_h_vcol, d_qt, d_ec &
2844	, fs_bound, fq_bound )
2845	END IF
2846
2847	if (mydebug) then
2848	call writefield_phy('u_seri',u_seri,llm)
2849	call writefield_phy('v_seri',v_seri,llm)
2850	call writefield_phy('t_seri',t_seri,llm)
2851	call writefield_phy('q_seri',q_seri,llm)
2852	endif
2853
2854	!
2855	!-------------------------------------------------------------------
2856	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2857	!-------------------------------------------------------------------
2858
2859	! 1. NUAGES CONVECTIFS
2860	!
2861	!IM cf FH
2862	! IF (iflag_cldth.eq.-1) THEN ! seulement pour Tiedtke
2863	IF (iflag_cldth.le.-1) THEN ! seulement pour Tiedtke
2864	snow_tiedtke=0.
2865	! print*,'avant calcul de la pseudo precip '
2866	! print*,'iflag_cldth',iflag_cldth
2867	if (iflag_cldth.eq.-1) then
2868	rain_tiedtke=rain_con
2869	else
2870	! print*,'calcul de la pseudo precip '
2871	rain_tiedtke=0.
2872	! print*,'calcul de la pseudo precip 0'
2873	do k=1,klev
2874	do i=1,klon
2875	if (d_q_con(i,k).lt.0.) then
2876	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
2877	*(paprs(i,k)-paprs(i,k+1))/rg
2878	endif
2879	enddo
2880	enddo
2881	endif
2882	!
2883	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2884	!
2885
2886	! Nuages diagnostiques pour Tiedtke
2887	CALL diagcld1(paprs,pplay, &
2888	!IM cf FH . rain_con,snow_con,ibas_con,itop_con,
2889	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
2890	diafra,dialiq)
2891	DO k = 1, klev
2892	DO i = 1, klon
2893	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2894	cldliq(i,k) = dialiq(i,k)
2895	cldfra(i,k) = diafra(i,k)
2896	ENDIF
2897	ENDDO
2898	ENDDO
2899
2900	ELSE IF (iflag_cldth.ge.3) THEN
2901	! On prend pour les nuages convectifs le max du calcul de la
2902	! convection et du calcul du pas de temps precedent diminue d'un facteur
2903	! facttemps
2904	facteur = pdtphys *facttemps
2905	do k=1,klev
2906	do i=1,klon
2907	rnebcon(i,k)=rnebcon(i,k)*facteur
2908	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k)) &
2909	then
2910	rnebcon(i,k)=rnebcon0(i,k)
2911	clwcon(i,k)=clwcon0(i,k)
2912	endif
2913	enddo
2914	enddo
2915
2916	!
2917	!jq - introduce the aerosol direct and first indirect radiative forcings
2918	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2919	IF (flag_aerosol .gt. 0) THEN
2920	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
2921	IF (.NOT. aerosol_couple) THEN
2922	!
2923	CALL readaerosol_optic( &
2924	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
2925	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
2926	mass_solu_aero, mass_solu_aero_pi, &
2927	tau_aero, piz_aero, cg_aero, &
2928	tausum_aero, tau3d_aero)
2929	ENDIF
2930	ELSE ! RRTM radiation
2931	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
2932	abort_message='config_inca=aero et rrtm=1 impossible'
2933	call abort_gcm(modname,abort_message,1)
2934	ELSE
2935	!
2936	#ifdef CPP_RRTM
2937	CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
2938	new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
2939	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
2940	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
2941	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
2942	tausum_aero, tau3d_aero)
2943
2944	CALL aeropt_lw_rrtm
2945	#else
2946
2947	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
2948	call abort_gcm(modname,abort_message,1)
2949	#endif
2950	!
2951	ENDIF
2952	ENDIF
2953	ELSE
2954	tausum_aero(:,:,:) = 0.
2955	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
2956	tau_aero(:,:,:,:) = 0.
2957	piz_aero(:,:,:,:) = 0.
2958	cg_aero(:,:,:,:) = 0.
2959	ELSE
2960	tau_aero_sw_rrtm(:,:,:,:)=0.0
2961	piz_aero_sw_rrtm(:,:,:,:)=0.0
2962	cg_aero_sw_rrtm(:,:,:,:)=0.0
2963	ENDIF
2964	ENDIF
2965	!
2966	!--STRAT AEROSOL
2967	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
2968	IF (flag_aerosol_strat) THEN
2969	IF (prt_level .GE.10) THEN
2970	PRINT *,'appel a readaerosolstrat', mth_cur
2971	ENDIF
2972	IF (iflag_rrtm.EQ.0) THEN
2973	CALL readaerosolstrato(debut)
2974	ELSE
2975	#ifdef CPP_RRTM
2976	CALL readaerosolstrato_rrtm(debut)
2977	#else
2978
2979	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
2980	call abort_gcm(modname,abort_message,1)
2981	#endif
2982	ENDIF
2983	ENDIF
2984	!--fin STRAT AEROSOL
2985
2986
2987	! On prend la somme des fractions nuageuses et des contenus en eau
2988
2989	if (iflag_cldth>=5) then
2990
2991	do k=1,klev
2992	ptconvth(:,k)=fm_therm(:,k+1)>0.
2993	enddo
2994
2995	if (iflag_coupl==4) then
2996
2997	! Dans le cas iflag_coupl==4, on prend la somme des convertures
2998	! convectives et lsc dans la partie des thermiques
2999	! Le controle par iflag_coupl est peut etre provisoire.
3000	do k=1,klev
3001	do i=1,klon
3002	if (ptconv(i,k).and.ptconvth(i,k)) then
3003	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3004	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3005	else if (ptconv(i,k)) then
3006	cldfra(i,k)=rnebcon(i,k)
3007	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3008	endif
3009	enddo
3010	enddo
3011
3012	else if (iflag_coupl==5) then
3013	do k=1,klev
3014	do i=1,klon
3015	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3016	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3017	enddo
3018	enddo
3019
3020	else
3021
3022	! Si on est sur un point touche par la convection profonde et pas
3023	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
3024	! de la convection profonde.
3025
3026	!IM/FH: 2011/02/23
3027	! definition des points sur lesquels ls thermiques sont actifs
3028
3029	do k=1,klev
3030	do i=1,klon
3031	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
3032	cldfra(i,k)=rnebcon(i,k)
3033	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3034	endif
3035	enddo
3036	enddo
3037
3038	endif
3039
3040	else
3041
3042	! Ancienne version
3043	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3044	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3045	endif
3046
3047	ENDIF
3048
3049	! plulsc(:)=0.
3050	! do k=1,klev,-1
3051	! do i=1,klon
3052	! zzz=prfl(:,k)+psfl(:,k)
3053	! if (.not.ptconvth.zzz.gt.0.)
3054	! enddo prfl, psfl,
3055	! enddo
3056	!
3057	! 2. NUAGES STARTIFORMES
3058	!
3059	IF (ok_stratus) THEN
3060	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3061	DO k = 1, klev
3062	DO i = 1, klon
3063	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3064	cldliq(i,k) = dialiq(i,k)
3065	cldfra(i,k) = diafra(i,k)
3066	ENDIF
3067	ENDDO
3068	ENDDO
3069	ENDIF
3070	!
3071	! Precipitation totale
3072	!
3073	DO i = 1, klon
3074	rain_fall(i) = rain_con(i) + rain_lsc(i)
3075	snow_fall(i) = snow_con(i) + snow_lsc(i)
3076	ENDDO
3077	!IM
3078	IF (ip_ebil_phy.ge.2) THEN
3079	ztit="after diagcld"
3080	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3081	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3082	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3083	call diagphy(airephy,ztit,ip_ebil_phy &
3084	, zero_v, zero_v, zero_v, zero_v, zero_v &
3085	, zero_v, zero_v, zero_v, ztsol &
3086	, d_h_vcol, d_qt, d_ec &
3087	, fs_bound, fq_bound )
3088	END IF
3089	!
3090	! Calculer l'humidite relative pour diagnostique
3091	!
3092	DO k = 1, klev
3093	DO i = 1, klon
3094	zx_t = t_seri(i,k)
3095	IF (thermcep) THEN
3096	if (iflag_ice_thermo.eq.0) then
3097	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3098	else
3099	zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))
3100	endif
3101	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3102	zx_qs = MIN(0.5,zx_qs)
3103	zcor = 1./(1.-retv*zx_qs)
3104	zx_qs = zx_qs*zcor
3105	ELSE
3106	IF (zx_t.LT.t_coup) THEN
3107	zx_qs = qsats(zx_t)/pplay(i,k)
3108	ELSE
3109	zx_qs = qsatl(zx_t)/pplay(i,k)
3110	ENDIF
3111	ENDIF
3112	zx_rh(i,k) = q_seri(i,k)/zx_qs
3113	zqsat(i,k)=zx_qs
3114	ENDDO
3115	ENDDO
3116
3117	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3118	! equivalente a 2m (tpote) pour diagnostique
3119	!
3120	DO i = 1, klon
3121	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3122	IF (thermcep) THEN
3123	IF(zt2m(i).LT.RTT) then
3124	Lheat=RLSTT
3125	ELSE
3126	Lheat=RLVTT
3127	ENDIF
3128	ELSE
3129	IF (zt2m(i).LT.RTT) THEN
3130	Lheat=RLSTT
3131	ELSE
3132	Lheat=RLVTT
3133	ENDIF
3134	ENDIF
3135	tpote(i) = tpot(i)* &
3136	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3137	ENDDO
3138
3139	IF (type_trac == 'inca') THEN
3140	#ifdef INCA
3141	CALL VTe(VTphysiq)
3142	CALL VTb(VTinca)
3143	calday = REAL(days_elapsed + 1) + jH_cur
3144
3145	call chemtime(itap+itau_phy-1, date0, dtime, itap)
3146	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
3147	CALL AEROSOL_METEO_CALC( &
3148	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3149	prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
3150	END IF
3151
3152	zxsnow_dummy(:) = 0.0
3153
3154	CALL chemhook_begin (calday, &
3155	days_elapsed+1, &
3156	jH_cur, &
3157	pctsrf(1,1), &
3158	rlat, &
3159	rlon, &
3160	airephy, &
3161	paprs, &
3162	pplay, &
3163	coefh(1:klon,1:klev,is_ave), &
3164	pphi, &
3165	t_seri, &
3166	u, &
3167	v, &
3168	wo(:, :, 1), &
3169	q_seri, &
3170	zxtsol, &
3171	zxsnow_dummy, &
3172	solsw, &
3173	albsol1, &
3174	rain_fall, &
3175	snow_fall, &
3176	itop_con, &
3177	ibas_con, &
3178	cldfra, &
3179	iim, &
3180	jjm, &
3181	tr_seri, &
3182	ftsol, &
3183	paprs, &
3184	cdragh, &
3185	cdragm, &
3186	pctsrf, &
3187	pdtphys, &
3188	itap)
3189
3190	CALL VTe(VTinca)
3191	CALL VTb(VTphysiq)
3192	#endif
3193	END IF !type_trac = inca
3194	!
3195	! Calculer les parametres optiques des nuages et quelques
3196	! parametres pour diagnostiques:
3197	!
3198
3199	IF (aerosol_couple) THEN
3200	mass_solu_aero(:,:) = ccm(:,:,1)
3201	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3202	END IF
3203
3204	if (ok_newmicro) then
3205	IF (iflag_rrtm.NE.0) THEN
3206	#ifdef CPP_RRTM
3207	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3208	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 pour ok_cdnc'
3209	call abort_gcm(modname,abort_message,1)
3210	endif
3211	#else
3212
3213	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
3214	call abort_gcm(modname,abort_message,1)
3215	#endif
3216	ENDIF
3217	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3218	paprs, pplay, t_seri, cldliq, cldfra, &
3219	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3220	flwp, fiwp, flwc, fiwc, &
3221	mass_solu_aero, mass_solu_aero_pi, &
3222	cldtaupi, re, fl, ref_liq, ref_ice, &
3223	ref_liq_pi, ref_ice_pi)
3224	else
3225	CALL nuage (paprs, pplay, &
3226	t_seri, cldliq, cldfra, cldtau, cldemi, &
3227	cldh, cldl, cldm, cldt, cldq, &
3228	ok_aie, &
3229	mass_solu_aero, mass_solu_aero_pi, &
3230	bl95_b0, bl95_b1, &
3231	cldtaupi, re, fl)
3232	endif
3233	!
3234	!IM betaCRF
3235	!
3236	cldtaurad = cldtau
3237	cldtaupirad = cldtaupi
3238	cldemirad = cldemi
3239
3240	!
3241	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3242	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3243	!
3244	! global
3245	!
3246	DO k=1, klev
3247	DO i=1, klon
3248	if (pplay(i,k).GE.pfree) THEN
3249	beta(i,k) = beta_pbl
3250	else
3251	beta(i,k) = beta_free
3252	endif
3253	if (mskocean_beta) THEN
3254	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3255	endif
3256	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3257	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3258	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3259	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3260	ENDDO
3261	ENDDO
3262	!
3263	else
3264	!
3265	! regional
3266	!
3267	DO k=1, klev
3268	DO i=1,klon
3269	!
3270	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
3271	rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
3272	if (pplay(i,k).GE.pfree) THEN
3273	beta(i,k) = beta_pbl
3274	else
3275	beta(i,k) = beta_free
3276	endif
3277	if (mskocean_beta) THEN
3278	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3279	endif
3280	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3281	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3282	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3283	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3284	endif
3285	!
3286	ENDDO
3287	ENDDO
3288	!
3289	endif
3290	!
3291	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3292	!
3293	IF (MOD(itaprad,radpas).EQ.0) THEN
3294
3295	DO i = 1, klon
3296	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce) &
3297	+ falb1(i,is_lic) * pctsrf(i,is_lic) &
3298	+ falb1(i,is_ter) * pctsrf(i,is_ter) &
3299	+ falb1(i,is_sic) * pctsrf(i,is_sic)
3300	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce) &
3301	+ falb2(i,is_lic) * pctsrf(i,is_lic) &
3302	+ falb2(i,is_ter) * pctsrf(i,is_ter) &
3303	+ falb2(i,is_sic) * pctsrf(i,is_sic)
3304	ENDDO
3305
3306	if (mydebug) then
3307	call writefield_phy('u_seri',u_seri,llm)
3308	call writefield_phy('v_seri',v_seri,llm)
3309	call writefield_phy('t_seri',t_seri,llm)
3310	call writefield_phy('q_seri',q_seri,llm)
3311	endif
3312
3313	IF (aerosol_couple) THEN
3314	#ifdef INCA
3315	CALL radlwsw_inca &
3316	(kdlon,kflev,dist, rmu0, fract, solaire, &
3317	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3318	wo(:, :, 1), &
3319	cldfrarad, cldemirad, cldtaurad, &
3320	heat,heat0,cool,cool0,radsol,albpla, &
3321	topsw,toplw,solsw,sollw, &
3322	sollwdown, &
3323	topsw0,toplw0,solsw0,sollw0, &
3324	lwdn0, lwdn, lwup0, lwup, &
3325	swdn0, swdn, swup0, swup, &
3326	ok_ade, ok_aie, &
3327	tau_aero, piz_aero, cg_aero, &
3328	topswad_aero, solswad_aero, &
3329	topswad0_aero, solswad0_aero, &
3330	topsw_aero, topsw0_aero, &
3331	solsw_aero, solsw0_aero, &
3332	cldtaupirad, &
3333	topswai_aero, solswai_aero)
3334
3335	#endif
3336	ELSE
3337	!
3338	!IM calcul radiatif pour le cas actuel
3339	!
3340	RCO2 = RCO2_act
3341	RCH4 = RCH4_act
3342	RN2O = RN2O_act
3343	RCFC11 = RCFC11_act
3344	RCFC12 = RCFC12_act
3345	!
3346	IF (prt_level .GE.10) THEN
3347	print *,' ->radlwsw, number 1 '
3348	ENDIF
3349	!
3350	CALL radlwsw &
3351	(dist, rmu0, fract, &
3352	paprs, pplay,zxtsol,albsol1, albsol2, &
3353	t_seri,q_seri,wo, &
3354	cldfrarad, cldemirad, cldtaurad, &
3355	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3356	flag_aerosol_strat, &
3357	tau_aero, piz_aero, cg_aero, &
3358	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! Rajoute par OB pour RRTM
3359	tau_aero_lw_rrtm, &
3360	cldtaupirad,new_aod, &
3361	zqsat, flwc, fiwc, &
3362	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3363	heat,heat0,cool,cool0,radsol,albpla, &
3364	topsw,toplw,solsw,sollw, &
3365	sollwdown, &
3366	topsw0,toplw0,solsw0,sollw0, &
3367	lwdn0, lwdn, lwup0, lwup, &
3368	swdn0, swdn, swup0, swup, &
3369	topswad_aero, solswad_aero, &
3370	topswai_aero, solswai_aero, &
3371	topswad0_aero, solswad0_aero, &
3372	topsw_aero, topsw0_aero, &
3373	solsw_aero, solsw0_aero, &
3374	topswcf_aero, solswcf_aero, &
3375	!-C. Kleinschmitt for LW diagnostics
3376	toplwad_aero, sollwad_aero,&
3377	toplwai_aero, sollwai_aero, &
3378	toplwad0_aero, sollwad0_aero,&
3379	!-end
3380	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3381	ZSWFT0_i, ZFSDN0, ZFSUP0)
3382
3383	!
3384	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3385	!IM des taux doit etre different du taux actuel
3386	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3387	!
3388	if (ok_4xCO2atm) then
3389	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3390	RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3391	RCFC12_per.NE.RCFC12_act) THEN
3392	!
3393	RCO2 = RCO2_per
3394	RCH4 = RCH4_per
3395	RN2O = RN2O_per
3396	RCFC11 = RCFC11_per
3397	RCFC12 = RCFC12_per
3398	!
3399	IF (prt_level .GE.10) THEN
3400	print *,' ->radlwsw, number 2 '
3401	ENDIF
3402	!
3403	CALL radlwsw &
3404	(dist, rmu0, fract, &
3405	paprs, pplay,zxtsol,albsol1, albsol2, &
3406	t_seri,q_seri,wo, &
3407	cldfra, cldemi, cldtau, &
3408	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3409	flag_aerosol_strat, &
3410	tau_aero, piz_aero, cg_aero, &
3411	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! Rajoute par OB pour RRTM
3412	tau_aero_lw_rrtm, &
3413	cldtaupi,new_aod, &
3414	zqsat, flwc, fiwc, &
3415	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3416	heatp,heat0p,coolp,cool0p,radsolp,albplap, &
3417	topswp,toplwp,solswp,sollwp, &
3418	sollwdownp, &
3419	topsw0p,toplw0p,solsw0p,sollw0p, &
3420	lwdn0p, lwdnp, lwup0p, lwupp, &
3421	swdn0p, swdnp, swup0p, swupp, &
3422	topswad_aerop, solswad_aerop, &
3423	topswai_aerop, solswai_aerop, &
3424	topswad0_aerop, solswad0_aerop, &
3425	topsw_aerop, topsw0_aerop, &
3426	solsw_aerop, solsw0_aerop, &
3427	topswcf_aerop, solswcf_aerop, &
3428	!-C. Kleinschmitt for LW diagnostics
3429	toplwad_aerop, sollwad_aerop,&
3430	toplwai_aerop, sollwai_aerop, &
3431	toplwad0_aerop, sollwad0_aerop,&
3432	!-end
3433	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3434	ZSWFT0_i, ZFSDN0, ZFSUP0)
3435	endif
3436	endif
3437	!
3438	ENDIF ! aerosol_couple
3439	itaprad = 0
3440	ENDIF ! MOD(itaprad,radpas)
3441	itaprad = itaprad + 1
3442
3443	IF (iflag_radia.eq.0) THEN
3444	IF (prt_level.ge.9) THEN
3445	PRINT *,'--------------------------------------------------'
3446	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3447	PRINT *,'>>>> heat et cool mis a zero '
3448	PRINT *,'--------------------------------------------------'
3449	END IF
3450	heat=0.
3451	cool=0.
3452	sollw=0. ! MPL 01032011
3453	solsw=0.
3454	radsol=0.
3455	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3456	swup0=0.
3457	swdn=0.
3458	swdn0=0.
3459	lwup=0.
3460	lwup0=0.
3461	lwdn=0.
3462	lwdn0=0.
3463	END IF
3464
3465	!
3466	! Ajouter la tendance des rayonnements (tous les pas)
3467	!
3468	d_t_swr(:,:)=heat(:,:)*dtime/RDAY
3469	d_t_lwr(:,:)=-cool(:,:)*dtime/RDAY
3470	d_t_sw0(:,:)=heat0(:,:)*dtime/RDAY
3471	d_t_lw0(:,:)=-cool0(:,:)*dtime/RDAY
3472	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW')
3473	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW')
3474
3475	!
3476	if (mydebug) then
3477	call writefield_phy('u_seri',u_seri,llm)
3478	call writefield_phy('v_seri',v_seri,llm)
3479	call writefield_phy('t_seri',t_seri,llm)
3480	call writefield_phy('q_seri',q_seri,llm)
3481	endif
3482
3483	!IM
3484	IF (ip_ebil_phy.ge.2) THEN
3485	ztit='after rad'
3486	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3487	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3488	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3489	call diagphy(airephy,ztit,ip_ebil_phy &
3490	, topsw, toplw, solsw, sollw, zero_v &
3491	, zero_v, zero_v, zero_v, ztsol &
3492	, d_h_vcol, d_qt, d_ec &
3493	, fs_bound, fq_bound )
3494	END IF
3495	!
3496	!
3497	! Calculer l'hydrologie de la surface
3498	!
3499	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3500	! . agesno, ftsol,fqsurf,fsnow, ruis)
3501	!
3502
3503	!
3504	! Calculer le bilan du sol et la derive de temperature (couplage)
3505	!
3506	DO i = 1, klon
3507	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3508	! a la demande de JLD
3509	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3510	ENDDO
3511	!
3512	!moddeblott(jan95)
3513	! Appeler le programme de parametrisation de l'orographie
3514	! a l'echelle sous-maille:
3515	!
3516	IF (prt_level .GE.10) THEN
3517	print *,' call orography ? ', ok_orodr
3518	ENDIF
3519	!
3520	IF (ok_orodr) THEN
3521	!
3522	! selection des points pour lesquels le shema est actif:
3523	igwd=0
3524	DO i=1,klon
3525	itest(i)=0
3526	! IF ((zstd(i).gt.10.0)) THEN
3527	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3528	itest(i)=1
3529	igwd=igwd+1
3530	idx(igwd)=i
3531	ENDIF
3532	ENDDO
3533	! igwdim=MAX(1,igwd)
3534	!
3535	IF (ok_strato) THEN
3536
3537	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
3538	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3539	igwd,idx,itest, &
3540	t_seri, u_seri, v_seri, &
3541	zulow, zvlow, zustrdr, zvstrdr, &
3542	d_t_oro, d_u_oro, d_v_oro)
3543
3544	ELSE
3545	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
3546	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3547	igwd,idx,itest, &
3548	t_seri, u_seri, v_seri, &
3549	zulow, zvlow, zustrdr, zvstrdr, &
3550	d_t_oro, d_u_oro, d_v_oro)
3551	ENDIF
3552	!
3553	! ajout des tendances
3554	!-----------------------------------------------------------------------------------------
3555	! ajout des tendances de la trainee de l'orographie
3556	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro')
3557	!-----------------------------------------------------------------------------------------
3558	!
3559	ENDIF ! fin de test sur ok_orodr
3560	!
3561	if (mydebug) then
3562	call writefield_phy('u_seri',u_seri,llm)
3563	call writefield_phy('v_seri',v_seri,llm)
3564	call writefield_phy('t_seri',t_seri,llm)
3565	call writefield_phy('q_seri',q_seri,llm)
3566	endif
3567
3568	IF (ok_orolf) THEN
3569	!
3570	! selection des points pour lesquels le shema est actif:
3571	igwd=0
3572	DO i=1,klon
3573	itest(i)=0
3574	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3575	itest(i)=1
3576	igwd=igwd+1
3577	idx(igwd)=i
3578	ENDIF
3579	ENDDO
3580	! igwdim=MAX(1,igwd)
3581	!
3582	IF (ok_strato) THEN
3583
3584	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
3585	rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3586	igwd,idx,itest, &
3587	t_seri, u_seri, v_seri, &
3588	zulow, zvlow, zustrli, zvstrli, &
3589	d_t_lif, d_u_lif, d_v_lif )
3590
3591	ELSE
3592	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
3593	rlat,zmea,zstd,zpic, &
3594	itest, &
3595	t_seri, u_seri, v_seri, &
3596	zulow, zvlow, zustrli, zvstrli, &
3597	d_t_lif, d_u_lif, d_v_lif)
3598	ENDIF
3599	!
3600	!-----------------------------------------------------------------------------------------
3601	! ajout des tendances de la portance de l'orographie
3602	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,dqi0,paprs,'lif')
3603	!-----------------------------------------------------------------------------------------
3604	!
3605	ENDIF ! fin de test sur ok_orolf
3606	! HINES GWD PARAMETRIZATION
3607
3608	IF (ok_hines) then
3609
3610	CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
3611	rlat,t_seri,u_seri,v_seri, &
3612	zustrhi,zvstrhi, &
3613	d_t_hin, d_u_hin, d_v_hin)
3614	!
3615	! ajout des tendances
3616	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,dqi0,paprs,'hin')
3617
3618	ENDIF
3619
3620	if (ok_gwd_rando) then
3621	call FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
3622	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
3623	du_gwd_rando, dv_gwd_rando)
3624	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,dqi0,paprs, &
3625	'flott_gwd_rando')
3626	end if
3627
3628	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3629
3630	if (mydebug) then
3631	call writefield_phy('u_seri',u_seri,llm)
3632	call writefield_phy('v_seri',v_seri,llm)
3633	call writefield_phy('t_seri',t_seri,llm)
3634	call writefield_phy('q_seri',q_seri,llm)
3635	endif
3636
3637	DO i = 1, klon
3638	zustrph(i)=0.
3639	zvstrph(i)=0.
3640	ENDDO
3641	DO k = 1, klev
3642	DO i = 1, klon
3643	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
3644	(paprs(i,k)-paprs(i,k+1))/rg
3645	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
3646	(paprs(i,k)-paprs(i,k+1))/rg
3647	ENDDO
3648	ENDDO
3649	!
3650	!IM calcul composantes axiales du moment angulaire et couple des montagnes
3651	!
3652	IF (is_sequential .and. ok_orodr) THEN
3653	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
3654	ra,rg,romega, &
3655	rlat,rlon,pphis, &
3656	zustrdr,zustrli,zustrph, &
3657	zvstrdr,zvstrli,zvstrph, &
3658	paprs,u,v, &
3659	aam, torsfc)
3660	ENDIF
3661	!IM cf. FLott END
3662	!IM
3663	IF (ip_ebil_phy.ge.2) THEN
3664	ztit='after orography'
3665	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3666	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3667	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3668	call diagphy(airephy,ztit,ip_ebil_phy &
3669	, zero_v, zero_v, zero_v, zero_v, zero_v &
3670	, zero_v, zero_v, zero_v, ztsol &
3671	, d_h_vcol, d_qt, d_ec &
3672	, fs_bound, fq_bound )
3673	END IF
3674
3675	!DC Calcul de la tendance due au methane
3676	IF(ok_qch4) THEN
3677	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
3678	! ajout de la tendance d'humidite due au methane
3679	CALL add_phys_tend(du0,dv0,dt0,d_q_ch4*dtime,dql0,'q_ch4')
3680	END IF
3681	!
3682	!
3683	!====================================================================
3684	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3685	!====================================================================
3686	! Abderrahmane 24.08.09
3687
3688	IF (ok_cosp) THEN
3689	! adeclarer
3690	#ifdef CPP_COSP
3691	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3692
3693	IF (prt_level .GE.10) THEN
3694	print*,'freq_cosp',freq_cosp
3695	ENDIF
3696	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3697	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3698	! s ref_liq,ref_ice
3699	call phys_cosp(itap,dtime,freq_cosp, &
3700	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
3701	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
3702	klon,klev,rlon,rlat,presnivs,overlap, &
3703	fract,ref_liq,ref_ice, &
3704	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
3705	zu10m,zv10m,pphis, &
3706	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
3707	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
3708	prfl(:,1:klev),psfl(:,1:klev), &
3709	pmflxr(:,1:klev),pmflxs(:,1:klev), &
3710	mr_ozone,cldtau, cldemi)
3711
3712	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3713	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
3714	! M clMISR,
3715	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3716	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
3717
3718	ENDIF
3719
3720	#endif
3721	ENDIF !ok_cosp
3722	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3723	!AA
3724	!AA Installation de l'interface online-offline pour traceurs
3725	!AA
3726	!====================================================================
3727	! Calcul des tendances traceurs
3728	!====================================================================
3729	!
3730
3731	IF (type_trac=='repr') THEN
3732	sh_in(:,:) = q_seri(:,:)
3733	ELSE
3734	sh_in(:,:) = qx(:,:,ivap)
3735	END IF
3736
3737	call phytrac ( &
3738	itap, days_elapsed+1, jH_cur, debut, &
3739	lafin, dtime, u, v, t, &
3740	paprs, pplay, pmfu, pmfd, &
3741	pen_u, pde_u, pen_d, pde_d, &
3742	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
3743	u1, v1, ftsol, pctsrf, &
3744	zustar, zu10m, zv10m, &
3745	wstar(:,is_ave), ale_bl, ale_wake, &
3746	rlat, rlon, &
3747	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
3748	presnivs, pphis, pphi, albsol1, &
3749	sh_in, rhcl, cldfra, rneb, &
3750	diafra, cldliq, itop_con, ibas_con, &
3751	pmflxr, pmflxs, prfl, psfl, &
3752	da, phi, mp, upwd, &
3753	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
3754	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
3755	ev, ep, epmlmMm, eplaMm, & !<<RomP
3756	dnwd, aerosol_couple, flxmass_w, &
3757	tau_aero, piz_aero, cg_aero, ccm, &
3758	rfname, &
3759	d_tr_dyn, & !<<RomP
3760	tr_seri)
3761
3762	IF (offline) THEN
3763
3764	IF (prt_level.ge.9) &
3765	print*,'Attention on met a 0 les thermiques pour phystoke'
3766	call phystokenc ( &
3767	nlon,klev,pdtphys,rlon,rlat, &
3768	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
3769	fm_therm,entr_therm, &
3770	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
3771	frac_impa, frac_nucl, &
3772	pphis,airephy,dtime,itap, &
3773	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
3774
3775
3776	ENDIF
3777
3778	!
3779	! Calculer le transport de l'eau et de l'energie (diagnostique)
3780	!
3781	CALL transp (paprs,zxtsol, &
3782	t_seri, q_seri, u_seri, v_seri, zphi, &
3783	ve, vq, ue, uq)
3784	!
3785	!IM global posePB BEG
3786	IF(1.EQ.0) THEN
3787	!
3788	CALL transp_lay (paprs,zxtsol, &
3789	t_seri, q_seri, u_seri, v_seri, zphi, &
3790	ve_lay, vq_lay, ue_lay, uq_lay)
3791	!
3792	ENDIF !(1.EQ.0) THEN
3793	!IM global posePB END
3794	! Accumuler les variables a stocker dans les fichiers histoire:
3795	!
3796
3797	!================================================================
3798	! Conversion of kinetic and potential energy into heat, for
3799	! parameterisation of subgrid-scale motions
3800	!================================================================
3801
3802	d_t_ec(:,:)=0.
3803	forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
3804	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
3805	u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
3806	zmasse,exner,d_t_ec)
3807	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
3808
3809	!IM
3810	IF (ip_ebil_phy.ge.1) THEN
3811	ztit='after physic'
3812	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
3813	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3814	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3815	! Comme les tendances de la physique sont ajoute dans la dynamique,
3816	! on devrait avoir que la variation d'entalpie par la dynamique
3817	! est egale a la variation de la physique au pas de temps precedent.
3818	! Donc la somme de ces 2 variations devrait etre nulle.
3819
3820	call diagphy(airephy,ztit,ip_ebil_phy &
3821	, topsw, toplw, solsw, sollw, sens &
3822	, evap, rain_fall, snow_fall, ztsol &
3823	, d_h_vcol, d_qt, d_ec &
3824	, fs_bound, fq_bound )
3825	!
3826	d_h_vcol_phy=d_h_vcol
3827	!
3828	END IF
3829	!
3830	!=======================================================================
3831	! SORTIES
3832	!=======================================================================
3833	!
3834	!IM initialisation + calculs divers diag AMIP2
3835	!
3836	include "calcul_divers.h"
3837	!
3838	!IM Interpolation sur les niveaux de pression du NMC
3839	! -------------------------------------------------
3840	!
3841	include "calcul_STDlev.h"
3842	!
3843	! slp sea level pressure
3844	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3845	!
3846	!cc prw = eau precipitable
3847	DO i = 1, klon
3848	prw(i) = 0.
3849	DO k = 1, klev
3850	prw(i) = prw(i) + &
3851	q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3852	ENDDO
3853	ENDDO
3854	!
3855	IF (type_trac == 'inca') THEN
3856	#ifdef INCA
3857	CALL VTe(VTphysiq)
3858	CALL VTb(VTinca)
3859
3860	CALL chemhook_end ( &
3861	dtime, &
3862	pplay, &
3863	t_seri, &
3864	tr_seri, &
3865	nbtr, &
3866	paprs, &
3867	q_seri, &
3868	airephy, &
3869	pphi, &
3870	pphis, &
3871	zx_rh)
3872
3873	CALL VTe(VTinca)
3874	CALL VTb(VTphysiq)
3875	#endif
3876	END IF
3877
3878
3879	!
3880	! Convertir les incrementations en tendances
3881	!
3882	IF (prt_level .GE.10) THEN
3883	print *,'Convertir les incrementations en tendances '
3884	ENDIF
3885	!
3886	if (mydebug) then
3887	call writefield_phy('u_seri',u_seri,llm)
3888	call writefield_phy('v_seri',v_seri,llm)
3889	call writefield_phy('t_seri',t_seri,llm)
3890	call writefield_phy('q_seri',q_seri,llm)
3891	endif
3892
3893	DO k = 1, klev
3894	DO i = 1, klon
3895	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3896	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3897	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3898	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3899	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3900	!CR: on ajoute le contenu en glace
3901	if (nqo.eq.3) then
3902	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
3903	endif
3904	ENDDO
3905	ENDDO
3906	!
3907	!CR: nb de traceurs eau: nqo
3908	! IF (nqtot.GE.3) THEN
3909	IF (nqtot.GE.(nqo+1)) THEN
3910	! DO iq = 3, nqtot
3911	DO iq = nqo+1, nqtot
3912	DO k = 1, klev
3913	DO i = 1, klon
3914	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3915	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
3916	ENDDO
3917	ENDDO
3918	ENDDO
3919	ENDIF
3920	!
3921	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3922	!IM global posePB include "write_bilKP_ins.h"
3923	!IM global posePB include "write_bilKP_ave.h"
3924	!
3925
3926	! Sauvegarder les valeurs de t et q a la fin de la physique:
3927	!
3928	DO k = 1, klev
3929	DO i = 1, klon
3930	u_ancien(i,k) = u_seri(i,k)
3931	v_ancien(i,k) = v_seri(i,k)
3932	t_ancien(i,k) = t_seri(i,k)
3933	q_ancien(i,k) = q_seri(i,k)
3934	ENDDO
3935	ENDDO
3936
3937	!!! RomP >>>
3938	!CR: nb de traceurs eau: nqo
3939	! IF (nqtot.GE.3) THEN
3940	IF (nqtot.GE.(nqo+1)) THEN
3941	! DO iq = 3, nqtot
3942	DO iq = nqo+1, nqtot
3943	DO k = 1, klev
3944	DO i = 1, klon
3945	! tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
3946	tr_ancien(i,k,iq-nqo) = tr_seri(i,k,iq-nqo)
3947	ENDDO
3948	ENDDO
3949	ENDDO
3950	ENDIF
3951	!!! RomP <<<
3952	!==========================================================================
3953	! Sorties des tendances pour un point particulier
3954	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3955	! pour le debug
3956	! La valeur de igout est attribuee plus haut dans le programme
3957	!==========================================================================
3958
3959	if (prt_level.ge.1) then
3960	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3961	write(lunout,*) &
3962	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3963	write(lunout,*) &
3964	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
3965	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
3966	pctsrf(igout,is_sic)
3967	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3968	do k=1,klev
3969	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
3970	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
3971	d_t_eva(igout,k)
3972	enddo
3973	write(lunout,*) 'cool,heat'
3974	do k=1,klev
3975	write(lunout,*) cool(igout,k),heat(igout,k)
3976	enddo
3977
3978	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3979	do k=1,klev
3980	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
3981	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3982	enddo
3983
3984	write(lunout,*) 'd_ps ',d_ps(igout)
3985	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3986	do k=1,klev
3987	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
3988	d_qx(igout,k,1),d_qx(igout,k,2)
3989	enddo
3990	endif
3991
3992	!==========================================================================
3993
3994	!============================================================
3995	! Calcul de la temperature potentielle
3996	!============================================================
3997	DO k = 1, klev
3998	DO i = 1, klon
3999	!JYG/IM theta en debut du pas de temps
4000	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4001	!JYG/IM theta en fin de pas de temps de physique
4002	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4003	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers MPL 20130625
4004	! fth_fonctions.F90 et parkind1.F90
4005	! sinon thetal=theta
4006	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4007	! : ql_seri(i,k))
4008	thetal(i,k)=theta(i,k)
4009	ENDDO
4010	ENDDO
4011	!
4012
4013	! 22.03.04 BEG
4014	!=============================================================
4015	! Ecriture des sorties
4016	!=============================================================
4017	#ifdef CPP_IOIPSL
4018
4019	! Recupere des varibles calcule dans differents modules
4020	! pour ecriture dans histxxx.nc
4021
4022	! Get some variables from module fonte_neige_mod
4023	CALL fonte_neige_get_vars(pctsrf, &
4024	zxfqcalving, zxfqfonte, zxffonte)
4025
4026
4027
4028
4029	!=============================================================
4030	! Separation entre thermiques et non thermiques dans les sorties
4031	! de fisrtilp
4032	!=============================================================
4033
4034	if (iflag_thermals>=1) then
4035	d_t_lscth=0.
4036	d_t_lscst=0.
4037	d_q_lscth=0.
4038	d_q_lscst=0.
4039	do k=1,klev
4040	do i=1,klon
4041	if (ptconvth(i,k)) then
4042	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4043	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4044	else
4045	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4046	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4047	endif
4048	enddo
4049	enddo
4050
4051	do i=1,klon
4052	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4053	plul_th(i)=prfl(i,1)+psfl(i,1)
4054	enddo
4055	endif
4056
4057
4058	!On effectue les sorties:
4059
4060	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
4061	pplay, lmax_th, aerosol_couple, &
4062	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
4063	ptconv, read_climoz, clevSTD, &
4064	ptconvth, d_t, qx, d_qx, zmasse, &
4065	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4066
4067
4068
4069
4070	include "write_histday_seri.h"
4071
4072	include "write_paramLMDZ_phy.h"
4073
4074	#endif
4075
4076	! 22.03.04 END
4077	!
4078	!====================================================================
4079	! Si c'est la fin, il faut conserver l'etat de redemarrage
4080	!====================================================================
4081	!
4082
4083	! -----------------------------------------------------------------
4084	! WSTATS: Saving statistics
4085	! -----------------------------------------------------------------
4086	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
4087	! which can later be used to make the statistic files of the run:
4088	! "stats") only possible in 3D runs !
4089
4090
4091	IF (callstats) THEN
4092
4093	call wstats(klon,o_psol%name,"Surface pressure","Pa" &
4094	,2,paprs(:,1))
4095	call wstats(klon,o_tsol%name,"Surface temperature","K", &
4096	2,zxtsol)
4097	zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
4098	call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
4099	"kg/(s*m2)",2,zx_tmp_fi2d)
4100	zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
4101	call wstats(klon,o_plul%name,"Large-scale Precip", &
4102	"kg/(s*m2)",2,zx_tmp_fi2d)
4103	zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
4104	call wstats(klon,o_pluc%name,"Convective Precip", &
4105	"kg/(s*m2)",2,zx_tmp_fi2d)
4106	call wstats(klon,o_sols%name,"Solar rad. at surf.", &
4107	"W/m2",2,solsw)
4108	call wstats(klon,o_soll%name,"IR rad. at surf.", &
4109	"W/m2",2,sollw)
4110	zx_tmp_fi2d(:) = topsw(:)-toplw(:)
4111	call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
4112	"W/m2",2,zx_tmp_fi2d)
4113
4114
4115
4116	call wstats(klon,o_temp%name,"Air temperature","K", &
4117	3,t_seri)
4118	call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
4119	3,u_seri)
4120	call wstats(klon,o_vitv%name,"Meridional wind", &
4121	"m.s-1",3,v_seri)
4122	call wstats(klon,o_vitw%name,"Vertical wind", &
4123	"m.s-1",3,omega)
4124	call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
4125	3,q_seri)
4126
4127
4128
4129	IF(lafin) THEN
4130	write (,) "Writing stats..."
4131	call mkstats(ierr)
4132	ENDIF
4133
4134	ENDIF !if callstats
4135
4136	IF (lafin) THEN
4137	itau_phy = itau_phy + itap
4138	CALL phyredem ("restartphy.nc")
4139	! open(97,form="unformatted",file="finbin")
4140	! write(97) u_seri,v_seri,t_seri,q_seri
4141	! close(97)
4142	!$OMP MASTER
4143	if (read_climoz >= 1) then
4144	if (is_mpi_root) then
4145	call nf95_close(ncid_climoz)
4146	end if
4147	deallocate(press_climoz) ! pointer
4148	end if
4149	!$OMP END MASTER
4150	ENDIF
4151
4152	! first=.false.
4153
4154	RETURN
4155	END SUBROUTINE physiq
4156	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
4157	IMPLICIT none
4158	!
4159	! Calculer et imprimer l'eau totale. A utiliser pour verifier
4160	! la conservation de l'eau
4161	!
4162	include "YOMCST.h"
4163	INTEGER klon,klev
4164	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
4165	REAL aire(klon)
4166	REAL qtotal, zx, qcheck
4167	INTEGER i, k
4168	!
4169	zx = 0.0
4170	DO i = 1, klon
4171	zx = zx + aire(i)
4172	ENDDO
4173	qtotal = 0.0
4174	DO k = 1, klev
4175	DO i = 1, klon
4176	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
4177	*(paprs(i,k)-paprs(i,k+1))/RG
4178	ENDDO
4179	ENDDO
4180	!
4181	qcheck = qtotal/zx
4182	!
4183	RETURN
4184	END FUNCTION qcheck
4185	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
4186	IMPLICIT none
4187	!
4188	! Tranformer une variable de la grille physique a
4189	! la grille d'ecriture
4190	!
4191	INTEGER nfield,nlon,iim,jjmp1, jjm
4192	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
4193	!
4194	INTEGER i, n, ig
4195	!
4196	jjm = jjmp1 - 1
4197	DO n = 1, nfield
4198	DO i=1,iim
4199	ecrit(i,n) = fi(1,n)
4200	ecrit(i+jjm*iim,n) = fi(nlon,n)
4201	ENDDO
4202	DO ig = 1, nlon - 2
4203	ecrit(iim+ig,n) = fi(1+ig,n)
4204	ENDDO
4205	ENDDO
4206	RETURN
4207	END SUBROUTINE gr_fi_ecrit

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

Download in other formats:

Original Format