Context Navigation

physiq.F90 @ 2225

Last change on this file since 2225 was 2224, checked in by crio, 10 years ago
cv3p1_closure: Bug fix in the toothpaste closure computation (dimension of klfc). physiq: check number of H2O tracers if thermodynamical effect of ice is included. Stop simulation if ice water is missing.
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.8 KB

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

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

Download in other formats:

Original Format