Context Navigation

physiq.F90 @ 3881

Last change on this file since 3881 was 3881, checked in by ymipsl, 9 years ago

fix uninitailized array when not using aerosol in LMDZ5

File size: 138.0 KB

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

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

Download in other formats:

Original Format