Context Navigation

physiq.F90 @ 3871

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

Avoid circular dependency with module

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	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3069	tau_aero(:,:,:,:) = 1.e-15
3070	piz_aero(:,:,:,:) = 1.
3071	cg_aero(:,:,:,:) = 0.
3072	ELSE
3073	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
3074	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3075	piz_aero_sw_rrtm(:,:,:,:) = 1.0
3076	cg_aero_sw_rrtm(:,:,:,:) = 0.0
3077	ENDIF
3078	ENDIF
3079	!
3080	!--STRAT AEROSOL
3081	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
3082	IF (flag_aerosol_strat) THEN
3083	IF (prt_level .GE.10) THEN
3084	PRINT *,'appel a readaerosolstrat', mth_cur
3085	ENDIF
3086	IF (iflag_rrtm.EQ.0) THEN
3087	CALL readaerosolstrato(debut)
3088	ELSE
3089	#ifdef CPP_RRTM
3090	CALL readaerosolstrato_rrtm(debut)
3091	#else
3092
3093	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
3094	call abort_physic(modname,abort_message,1)
3095	#endif
3096	ENDIF
3097	ENDIF
3098	!--fin STRAT AEROSOL
3099
3100
3101	! On prend la somme des fractions nuageuses et des contenus en eau
3102
3103	if (iflag_cld_th>=5) then
3104
3105	do k=1,klev
3106	ptconvth(:,k)=fm_therm(:,k+1)>0.
3107	enddo
3108
3109	if (iflag_coupl==4) then
3110
3111	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3112	! convectives et lsc dans la partie des thermiques
3113	! Le controle par iflag_coupl est peut etre provisoire.
3114	do k=1,klev
3115	do i=1,klon
3116	if (ptconv(i,k).and.ptconvth(i,k)) then
3117	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3118	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3119	else if (ptconv(i,k)) then
3120	cldfra(i,k)=rnebcon(i,k)
3121	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3122	endif
3123	enddo
3124	enddo
3125
3126	else if (iflag_coupl==5) then
3127	do k=1,klev
3128	do i=1,klon
3129	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3130	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3131	enddo
3132	enddo
3133
3134	else
3135
3136	! Si on est sur un point touche par la convection profonde et pas
3137	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
3138	! de la convection profonde.
3139
3140	!IM/FH: 2011/02/23
3141	! definition des points sur lesquels ls thermiques sont actifs
3142
3143	do k=1,klev
3144	do i=1,klon
3145	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
3146	cldfra(i,k)=rnebcon(i,k)
3147	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3148	endif
3149	enddo
3150	enddo
3151
3152	endif
3153
3154	else
3155
3156	! Ancienne version
3157	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3158	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3159	endif
3160
3161	ENDIF
3162
3163	! plulsc(:)=0.
3164	! do k=1,klev,-1
3165	! do i=1,klon
3166	! zzz=prfl(:,k)+psfl(:,k)
3167	! if (.not.ptconvth.zzz.gt.0.)
3168	! enddo prfl, psfl,
3169	! enddo
3170	!
3171	! 2. NUAGES STARTIFORMES
3172	!
3173	IF (ok_stratus) THEN
3174	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3175	DO k = 1, klev
3176	DO i = 1, klon
3177	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3178	cldliq(i,k) = dialiq(i,k)
3179	cldfra(i,k) = diafra(i,k)
3180	ENDIF
3181	ENDDO
3182	ENDDO
3183	ENDIF
3184	!
3185	! Precipitation totale
3186	!
3187	DO i = 1, klon
3188	rain_fall(i) = rain_con(i) + rain_lsc(i)
3189	snow_fall(i) = snow_con(i) + snow_lsc(i)
3190	ENDDO
3191	!IM
3192	IF (ip_ebil_phy.ge.2) THEN
3193	ztit="after diagcld"
3194	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,pdtphys &
3195	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3196	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3197	call diagphy(cell_area,ztit,ip_ebil_phy &
3198	, zero_v, zero_v, zero_v, zero_v, zero_v &
3199	, zero_v, zero_v, zero_v, ztsol &
3200	, d_h_vcol, d_qt, d_ec &
3201	, fs_bound, fq_bound )
3202	END IF
3203	!
3204	! Calculer l'humidite relative pour diagnostique
3205	!
3206	DO k = 1, klev
3207	DO i = 1, klon
3208	zx_t = t_seri(i,k)
3209	IF (thermcep) THEN
3210	if (iflag_ice_thermo.eq.0) then
3211	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3212	else
3213	zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))
3214	endif
3215	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3216	zx_qs = MIN(0.5,zx_qs)
3217	zcor = 1./(1.-retv*zx_qs)
3218	zx_qs = zx_qs*zcor
3219	ELSE
3220	IF (zx_t.LT.t_coup) THEN
3221	zx_qs = qsats(zx_t)/pplay(i,k)
3222	ELSE
3223	zx_qs = qsatl(zx_t)/pplay(i,k)
3224	ENDIF
3225	ENDIF
3226	zx_rh(i,k) = q_seri(i,k)/zx_qs
3227	zqsat(i,k)=zx_qs
3228	ENDDO
3229	ENDDO
3230
3231	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3232	! equivalente a 2m (tpote) pour diagnostique
3233	!
3234	DO i = 1, klon
3235	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3236	IF (thermcep) THEN
3237	IF(zt2m(i).LT.RTT) then
3238	Lheat=RLSTT
3239	ELSE
3240	Lheat=RLVTT
3241	ENDIF
3242	ELSE
3243	IF (zt2m(i).LT.RTT) THEN
3244	Lheat=RLSTT
3245	ELSE
3246	Lheat=RLVTT
3247	ENDIF
3248	ENDIF
3249	tpote(i) = tpot(i)* &
3250	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3251	ENDDO
3252
3253	IF (type_trac == 'inca') THEN
3254	#ifdef INCA
3255	CALL VTe(VTphysiq)
3256	CALL VTb(VTinca)
3257	calday = REAL(days_elapsed + 1) + jH_cur
3258
3259	call chemtime(itap+itau_phy-1, date0, pdtphys, itap)
3260	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
3261	CALL AEROSOL_METEO_CALC( &
3262	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3263	prfl,psfl,pctsrf,cell_area,lat_degrees,lon_degrees,u10m,v10m)
3264	END IF
3265
3266	zxsnow_dummy(:) = 0.0
3267
3268	CALL chemhook_begin (calday, &
3269	days_elapsed+1, &
3270	jH_cur, &
3271	pctsrf(1,1), &
3272	lat_degrees, &
3273	lon_degrees, &
3274	cell_area, &
3275	paprs, &
3276	pplay, &
3277	coefh(1:klon,1:klev,is_ave), &
3278	pphi, &
3279	t_seri, &
3280	u, &
3281	v, &
3282	wo(:, :, 1), &
3283	q_seri, &
3284	zxtsol, &
3285	zxsnow_dummy, &
3286	solsw, &
3287	albsol1, &
3288	rain_fall, &
3289	snow_fall, &
3290	itop_con, &
3291	ibas_con, &
3292	cldfra, &
3293	iim, &
3294	jjm, &
3295	tr_seri, &
3296	ftsol, &
3297	paprs, &
3298	cdragh, &
3299	cdragm, &
3300	pctsrf, &
3301	pdtphys, &
3302	itap)
3303
3304	CALL VTe(VTinca)
3305	CALL VTb(VTphysiq)
3306	#endif
3307	END IF !type_trac = inca
3308	!
3309	! Calculer les parametres optiques des nuages et quelques
3310	! parametres pour diagnostiques:
3311	!
3312
3313	IF (aerosol_couple) THEN
3314	mass_solu_aero(:,:) = ccm(:,:,1)
3315	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3316	END IF
3317
3318	if (ok_newmicro) then
3319	IF (iflag_rrtm.NE.0) THEN
3320	#ifdef CPP_RRTM
3321	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3322	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 pour ok_cdnc'
3323	call abort_physic(modname,abort_message,1)
3324	endif
3325	#else
3326
3327	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
3328	call abort_physic(modname,abort_message,1)
3329	#endif
3330	ENDIF
3331	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3332	paprs, pplay, t_seri, cldliq, cldfra, &
3333	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3334	flwp, fiwp, flwc, fiwc, &
3335	mass_solu_aero, mass_solu_aero_pi, &
3336	cldtaupi, re, fl, ref_liq, ref_ice, &
3337	ref_liq_pi, ref_ice_pi)
3338	else
3339	CALL nuage (paprs, pplay, &
3340	t_seri, cldliq, cldfra, cldtau, cldemi, &
3341	cldh, cldl, cldm, cldt, cldq, &
3342	ok_aie, &
3343	mass_solu_aero, mass_solu_aero_pi, &
3344	bl95_b0, bl95_b1, &
3345	cldtaupi, re, fl)
3346	endif
3347	!
3348	!IM betaCRF
3349	!
3350	cldtaurad = cldtau
3351	cldtaupirad = cldtaupi
3352	cldemirad = cldemi
3353
3354	!
3355	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3356	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3357	!
3358	! global
3359	!
3360	DO k=1, klev
3361	DO i=1, klon
3362	if (pplay(i,k).GE.pfree) THEN
3363	beta(i,k) = beta_pbl
3364	else
3365	beta(i,k) = beta_free
3366	endif
3367	if (mskocean_beta) THEN
3368	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3369	endif
3370	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3371	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3372	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3373	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3374	ENDDO
3375	ENDDO
3376	!
3377	else
3378	!
3379	! regional
3380	!
3381	DO k=1, klev
3382	DO i=1,klon
3383	!
3384	if (lon_degrees(i).ge.lon1_beta.AND.lat_degrees(i).le.lon2_beta.AND. &
3385	lat_degrees(i).le.lat1_beta.AND.lat_degrees(i).ge.lat2_beta) THEN
3386	if (pplay(i,k).GE.pfree) THEN
3387	beta(i,k) = beta_pbl
3388	else
3389	beta(i,k) = beta_free
3390	endif
3391	if (mskocean_beta) THEN
3392	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3393	endif
3394	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3395	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3396	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3397	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3398	endif
3399	!
3400	ENDDO
3401	ENDDO
3402	!
3403	endif
3404	!
3405	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3406	!
3407	IF (MOD(itaprad,radpas).EQ.0) THEN
3408
3409	!albedo SB >>>
3410	if(ok_chlorophyll)then
3411	print*,"-- reading chlorophyll"
3412	call readchlorophyll(debut)
3413	endif
3414	!do i=1,klon
3415	!if(chl_con(i)>1.) print*,i,chl_con(i),pctsrf(i,is_ter)
3416	!enddo
3417	!albedo SB <<<
3418
3419
3420	if (mydebug) then
3421	call writefield_phy('u_seri',u_seri,nbp_lev)
3422	call writefield_phy('v_seri',v_seri,nbp_lev)
3423	call writefield_phy('t_seri',t_seri,nbp_lev)
3424	call writefield_phy('q_seri',q_seri,nbp_lev)
3425	endif
3426
3427	IF (aerosol_couple) THEN
3428	#ifdef INCA
3429	CALL radlwsw_inca &
3430	(kdlon,kflev,dist, rmu0, fract, solaire, &
3431	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3432	wo(:, :, 1), &
3433	cldfrarad, cldemirad, cldtaurad, &
3434	heat,heat0,cool,cool0,radsol,albpla, &
3435	topsw,toplw,solsw,sollw, &
3436	sollwdown, &
3437	topsw0,toplw0,solsw0,sollw0, &
3438	lwdn0, lwdn, lwup0, lwup, &
3439	swdn0, swdn, swup0, swup, &
3440	ok_ade, ok_aie, &
3441	tau_aero, piz_aero, cg_aero, &
3442	topswad_aero, solswad_aero, &
3443	topswad0_aero, solswad0_aero, &
3444	topsw_aero, topsw0_aero, &
3445	solsw_aero, solsw0_aero, &
3446	cldtaupirad, &
3447	topswai_aero, solswai_aero)
3448
3449	#endif
3450	ELSE
3451	!
3452	!IM calcul radiatif pour le cas actuel
3453	!
3454	RCO2 = RCO2_act
3455	RCH4 = RCH4_act
3456	RN2O = RN2O_act
3457	RCFC11 = RCFC11_act
3458	RCFC12 = RCFC12_act
3459	!
3460	IF (prt_level .GE.10) THEN
3461	print *,' ->radlwsw, number 1 '
3462	ENDIF
3463	!
3464	CALL radlwsw &
3465	(dist, rmu0, fract, &
3466	!albedo SB >>>
3467	! paprs, pplay,zxtsol,albsol1, albsol2, &
3468	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3469	!albedo SB <<<
3470	t_seri,q_seri,wo, &
3471	cldfrarad, cldemirad, cldtaurad, &
3472	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3473	flag_aerosol_strat, &
3474	tau_aero, piz_aero, cg_aero, &
3475	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! Rajoute par OB pour RRTM
3476	tau_aero_lw_rrtm, &
3477	cldtaupirad,new_aod, &
3478	zqsat, flwc, fiwc, &
3479	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3480	heat,heat0,cool,cool0,radsol,albpla, &
3481	topsw,toplw,solsw,sollw, &
3482	sollwdown, &
3483	topsw0,toplw0,solsw0,sollw0, &
3484	lwdn0, lwdn, lwup0, lwup, &
3485	swdn0, swdn, swup0, swup, &
3486	topswad_aero, solswad_aero, &
3487	topswai_aero, solswai_aero, &
3488	topswad0_aero, solswad0_aero, &
3489	topsw_aero, topsw0_aero, &
3490	solsw_aero, solsw0_aero, &
3491	topswcf_aero, solswcf_aero, &
3492	!-C. Kleinschmitt for LW diagnostics
3493	toplwad_aero, sollwad_aero,&
3494	toplwai_aero, sollwai_aero, &
3495	toplwad0_aero, sollwad0_aero,&
3496	!-end
3497	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3498	ZSWFT0_i, ZFSDN0, ZFSUP0)
3499
3500	!
3501	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3502	!IM des taux doit etre different du taux actuel
3503	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3504	!
3505	if (ok_4xCO2atm) then
3506	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3507	RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3508	RCFC12_per.NE.RCFC12_act) THEN
3509	!
3510	RCO2 = RCO2_per
3511	RCH4 = RCH4_per
3512	RN2O = RN2O_per
3513	RCFC11 = RCFC11_per
3514	RCFC12 = RCFC12_per
3515	!
3516	IF (prt_level .GE.10) THEN
3517	print *,' ->radlwsw, number 2 '
3518	ENDIF
3519	!
3520	CALL radlwsw &
3521	(dist, rmu0, fract, &
3522	!albedo SB >>>
3523	! paprs, pplay,zxtsol,albsol1, albsol2, &
3524	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3525	!albedo SB <<<
3526	t_seri,q_seri,wo, &
3527	cldfra, cldemi, cldtau, &
3528	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3529	flag_aerosol_strat, &
3530	tau_aero, piz_aero, cg_aero, &
3531	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! Rajoute par OB pour RRTM
3532	tau_aero_lw_rrtm, &
3533	cldtaupi,new_aod, &
3534	zqsat, flwc, fiwc, &
3535	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3536	heatp,heat0p,coolp,cool0p,radsolp,albplap, &
3537	topswp,toplwp,solswp,sollwp, &
3538	sollwdownp, &
3539	topsw0p,toplw0p,solsw0p,sollw0p, &
3540	lwdn0p, lwdnp, lwup0p, lwupp, &
3541	swdn0p, swdnp, swup0p, swupp, &
3542	topswad_aerop, solswad_aerop, &
3543	topswai_aerop, solswai_aerop, &
3544	topswad0_aerop, solswad0_aerop, &
3545	topsw_aerop, topsw0_aerop, &
3546	solsw_aerop, solsw0_aerop, &
3547	topswcf_aerop, solswcf_aerop, &
3548	!-C. Kleinschmitt for LW diagnostics
3549	toplwad_aerop, sollwad_aerop,&
3550	toplwai_aerop, sollwai_aerop, &
3551	toplwad0_aerop, sollwad0_aerop,&
3552	!-end
3553	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3554	ZSWFT0_i, ZFSDN0, ZFSUP0)
3555	endif
3556	endif
3557	!
3558	ENDIF ! aerosol_couple
3559	itaprad = 0
3560	ENDIF ! MOD(itaprad,radpas)
3561	itaprad = itaprad + 1
3562
3563	IF (iflag_radia.eq.0) THEN
3564	IF (prt_level.ge.9) THEN
3565	PRINT *,'--------------------------------------------------'
3566	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3567	PRINT *,'>>>> heat et cool mis a zero '
3568	PRINT *,'--------------------------------------------------'
3569	END IF
3570	heat=0.
3571	cool=0.
3572	sollw=0. ! MPL 01032011
3573	solsw=0.
3574	radsol=0.
3575	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3576	swup0=0.
3577	swdn=0.
3578	swdn0=0.
3579	lwup=0.
3580	lwup0=0.
3581	lwdn=0.
3582	lwdn0=0.
3583	END IF
3584
3585	!
3586	! Ajouter la tendance des rayonnements (tous les pas)
3587	!
3588	d_t_swr(:,:)=heat(:,:)*pdtphys/RDAY
3589	d_t_lwr(:,:)=-cool(:,:)*pdtphys/RDAY
3590	d_t_sw0(:,:)=heat0(:,:)*pdtphys/RDAY
3591	d_t_lw0(:,:)=-cool0(:,:)*pdtphys/RDAY
3592	CALL add_phys_tend(du0,dv0,d_t_swr,dq0,dql0,dqi0,paprs,'SW',abortphy)
3593	CALL add_phys_tend(du0,dv0,d_t_lwr,dq0,dql0,dqi0,paprs,'LW',abortphy)
3594
3595	!
3596	if (mydebug) then
3597	call writefield_phy('u_seri',u_seri,nbp_lev)
3598	call writefield_phy('v_seri',v_seri,nbp_lev)
3599	call writefield_phy('t_seri',t_seri,nbp_lev)
3600	call writefield_phy('q_seri',q_seri,nbp_lev)
3601	endif
3602
3603	!IM
3604	IF (ip_ebil_phy.ge.2) THEN
3605	ztit='after rad'
3606	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,pdtphys &
3607	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3608	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3609	call diagphy(cell_area,ztit,ip_ebil_phy &
3610	, topsw, toplw, solsw, sollw, zero_v &
3611	, zero_v, zero_v, zero_v, ztsol &
3612	, d_h_vcol, d_qt, d_ec &
3613	, fs_bound, fq_bound )
3614	END IF
3615	!
3616	!
3617	! Calculer l'hydrologie de la surface
3618	!
3619	! CALL hydrol(pdtphys,pctsrf,rain_fall, snow_fall, zxevap,
3620	! . agesno, ftsol,fqsurf,fsnow, ruis)
3621	!
3622
3623	!
3624	! Calculer le bilan du sol et la derive de temperature (couplage)
3625	!
3626	DO i = 1, klon
3627	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3628	! a la demande de JLD
3629	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3630	ENDDO
3631	!
3632	!moddeblott(jan95)
3633	! Appeler le programme de parametrisation de l'orographie
3634	! a l'echelle sous-maille:
3635	!
3636	IF (prt_level .GE.10) THEN
3637	print *,' call orography ? ', ok_orodr
3638	ENDIF
3639	!
3640	IF (ok_orodr) THEN
3641	!
3642	! selection des points pour lesquels le shema est actif:
3643	igwd=0
3644	DO i=1,klon
3645	itest(i)=0
3646	! IF ((zstd(i).gt.10.0)) THEN
3647	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3648	itest(i)=1
3649	igwd=igwd+1
3650	idx(igwd)=i
3651	ENDIF
3652	ENDDO
3653	! igwdim=MAX(1,igwd)
3654	!
3655	IF (ok_strato) THEN
3656
3657	CALL drag_noro_strato(klon,klev,pdtphys,paprs,pplay, &
3658	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3659	igwd,idx,itest, &
3660	t_seri, u_seri, v_seri, &
3661	zulow, zvlow, zustrdr, zvstrdr, &
3662	d_t_oro, d_u_oro, d_v_oro)
3663
3664	ELSE
3665	CALL drag_noro(klon,klev,pdtphys,paprs,pplay, &
3666	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3667	igwd,idx,itest, &
3668	t_seri, u_seri, v_seri, &
3669	zulow, zvlow, zustrdr, zvstrdr, &
3670	d_t_oro, d_u_oro, d_v_oro)
3671	ENDIF
3672	!
3673	! ajout des tendances
3674	!-----------------------------------------------------------------------------------------
3675	! ajout des tendances de la trainee de l'orographie
3676	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro',abortphy)
3677	!-----------------------------------------------------------------------------------------
3678	!
3679	ENDIF ! fin de test sur ok_orodr
3680	!
3681	if (mydebug) then
3682	call writefield_phy('u_seri',u_seri,nbp_lev)
3683	call writefield_phy('v_seri',v_seri,nbp_lev)
3684	call writefield_phy('t_seri',t_seri,nbp_lev)
3685	call writefield_phy('q_seri',q_seri,nbp_lev)
3686	endif
3687
3688	IF (ok_orolf) THEN
3689	!
3690	! selection des points pour lesquels le shema est actif:
3691	igwd=0
3692	DO i=1,klon
3693	itest(i)=0
3694	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3695	itest(i)=1
3696	igwd=igwd+1
3697	idx(igwd)=i
3698	ENDIF
3699	ENDDO
3700	! igwdim=MAX(1,igwd)
3701	!
3702	IF (ok_strato) THEN
3703
3704	CALL lift_noro_strato(klon,klev,pdtphys,paprs,pplay, &
3705	lat_degrees,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3706	igwd,idx,itest, &
3707	t_seri, u_seri, v_seri, &
3708	zulow, zvlow, zustrli, zvstrli, &
3709	d_t_lif, d_u_lif, d_v_lif )
3710
3711	ELSE
3712	CALL lift_noro(klon,klev,pdtphys,paprs,pplay, &
3713	lat_degrees,zmea,zstd,zpic, &
3714	itest, &
3715	t_seri, u_seri, v_seri, &
3716	zulow, zvlow, zustrli, zvstrli, &
3717	d_t_lif, d_u_lif, d_v_lif)
3718	ENDIF
3719	!
3720	!-----------------------------------------------------------------------------------------
3721	! ajout des tendances de la portance de l'orographie
3722	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,dqi0,paprs,'lif',abortphy)
3723	!-----------------------------------------------------------------------------------------
3724	!
3725	ENDIF ! fin de test sur ok_orolf
3726	! HINES GWD PARAMETRIZATION
3727
3728	IF (ok_hines) then
3729
3730	CALL hines_gwd(klon,klev,pdtphys,paprs,pplay, &
3731	lat_degrees,t_seri,u_seri,v_seri, &
3732	zustrhi,zvstrhi, &
3733	d_t_hin, d_u_hin, d_v_hin)
3734	!
3735	! ajout des tendances
3736	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,dqi0,paprs,'hin',abortphy)
3737
3738	ENDIF
3739
3740	if (ok_gwd_rando) then
3741	call FLOTT_GWD_rando(pdtphys, pplay, t_seri, u_seri, v_seri, &
3742	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
3743	du_gwd_rando, dv_gwd_rando)
3744	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,dqi0,paprs, &
3745	'flott_gwd_rando',abortphy)
3746	end if
3747
3748	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3749
3750	if (mydebug) then
3751	call writefield_phy('u_seri',u_seri,nbp_lev)
3752	call writefield_phy('v_seri',v_seri,nbp_lev)
3753	call writefield_phy('t_seri',t_seri,nbp_lev)
3754	call writefield_phy('q_seri',q_seri,nbp_lev)
3755	endif
3756
3757	DO i = 1, klon
3758	zustrph(i)=0.
3759	zvstrph(i)=0.
3760	ENDDO
3761	DO k = 1, klev
3762	DO i = 1, klon
3763	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/pdtphys* &
3764	(paprs(i,k)-paprs(i,k+1))/rg
3765	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/pdtphys* &
3766	(paprs(i,k)-paprs(i,k+1))/rg
3767	ENDDO
3768	ENDDO
3769	!
3770	!IM calcul composantes axiales du moment angulaire et couple des montagnes
3771	!
3772	IF (is_sequential .and. ok_orodr) THEN
3773	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
3774	ra,rg,romega, &
3775	lat_degrees,lon_degrees,pphis, &
3776	zustrdr,zustrli,zustrph, &
3777	zvstrdr,zvstrli,zvstrph, &
3778	paprs,u,v, &
3779	aam, torsfc)
3780	ENDIF
3781	!IM cf. FLott END
3782	!IM
3783	IF (ip_ebil_phy.ge.2) THEN
3784	ztit='after orography'
3785	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,pdtphys &
3786	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3787	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3788	call diagphy(cell_area,ztit,ip_ebil_phy &
3789	, zero_v, zero_v, zero_v, zero_v, zero_v &
3790	, zero_v, zero_v, zero_v, ztsol &
3791	, d_h_vcol, d_qt, d_ec &
3792	, fs_bound, fq_bound )
3793	END IF
3794
3795	!DC Calcul de la tendance due au methane
3796	IF(ok_qch4) THEN
3797	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
3798	! ajout de la tendance d'humidite due au methane
3799	CALL add_phys_tend(du0,dv0,dt0,d_q_ch4*pdtphys,dql0,'q_ch4',abortphy)
3800	END IF
3801	!
3802	!
3803	!====================================================================
3804	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3805	!====================================================================
3806	! Abderrahmane 24.08.09
3807
3808	IF (ok_cosp) THEN
3809	! adeclarer
3810	#ifdef CPP_COSP
3811	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/pdtphys)).EQ.0) THEN
3812
3813	IF (prt_level .GE.10) THEN
3814	print*,'freq_cosp',freq_cosp
3815	ENDIF
3816	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3817	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3818	! s ref_liq,ref_ice
3819	call phys_cosp(itap,pdtphys,freq_cosp, &
3820	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
3821	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
3822	klon,klev,lon_degrees,lat_degrees,presnivs,overlap, &
3823	fract,ref_liq,ref_ice, &
3824	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
3825	zu10m,zv10m,pphis, &
3826	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
3827	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
3828	prfl(:,1:klev),psfl(:,1:klev), &
3829	pmflxr(:,1:klev),pmflxs(:,1:klev), &
3830	mr_ozone,cldtau, cldemi)
3831
3832	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3833	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
3834	! M clMISR,
3835	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3836	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
3837
3838	ENDIF
3839
3840	#endif
3841	ENDIF !ok_cosp
3842	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3843	!AA
3844	!AA Installation de l'interface online-offline pour traceurs
3845	!AA
3846	!====================================================================
3847	! Calcul des tendances traceurs
3848	!====================================================================
3849	!
3850
3851	IF (type_trac=='repr') THEN
3852	sh_in(:,:) = q_seri(:,:)
3853	ELSE
3854	sh_in(:,:) = qx(:,:,ivap)
3855	END IF
3856
3857	call phytrac ( &
3858	itap, days_elapsed+1, jH_cur, debut, &
3859	lafin, pdtphys, u, v, t, &
3860	paprs, pplay, pmfu, pmfd, &
3861	pen_u, pde_u, pen_d, pde_d, &
3862	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
3863	u1, v1, ftsol, pctsrf, &
3864	zustar, zu10m, zv10m, &
3865	wstar(:,is_ave), ale_bl, ale_wake, &
3866	lat_degrees, lon_degrees, &
3867	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
3868	presnivs, pphis, pphi, albsol1, &
3869	sh_in, rhcl, cldfra, rneb, &
3870	diafra, cldliq, itop_con, ibas_con, &
3871	pmflxr, pmflxs, prfl, psfl, &
3872	da, phi, mp, upwd, &
3873	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
3874	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
3875	ev, ep, epmlmMm, eplaMm, & !<<RomP
3876	dnwd, aerosol_couple, flxmass_w, &
3877	tau_aero, piz_aero, cg_aero, ccm, &
3878	rfname, &
3879	d_tr_dyn, & !<<RomP
3880	tr_seri)
3881
3882	IF (offline) THEN
3883
3884	IF (prt_level.ge.9) &
3885	print*,'Attention on met a 0 les thermiques pour phystoke'
3886	call phystokenc ( &
3887	nlon,klev,pdtphys,lon_degrees,lat_degrees, &
3888	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
3889	fm_therm,entr_therm, &
3890	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
3891	frac_impa, frac_nucl, &
3892	pphis,cell_area,pdtphys,itap, &
3893	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
3894
3895
3896	ENDIF
3897
3898	!
3899	! Calculer le transport de l'eau et de l'energie (diagnostique)
3900	!
3901	CALL transp (paprs,zxtsol, &
3902	t_seri, q_seri, u_seri, v_seri, zphi, &
3903	ve, vq, ue, uq)
3904	!
3905	!IM global posePB BEG
3906	IF(1.EQ.0) THEN
3907	!
3908	CALL transp_lay (paprs,zxtsol, &
3909	t_seri, q_seri, u_seri, v_seri, zphi, &
3910	ve_lay, vq_lay, ue_lay, uq_lay)
3911	!
3912	ENDIF !(1.EQ.0) THEN
3913	!IM global posePB END
3914	! Accumuler les variables a stocker dans les fichiers histoire:
3915	!
3916
3917	!================================================================
3918	! Conversion of kinetic and potential energy into heat, for
3919	! parameterisation of subgrid-scale motions
3920	!================================================================
3921
3922	d_t_ec(:,:)=0.
3923	forall (k=1: nbp_lev) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
3924	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
3925	u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
3926	zmasse,exner,d_t_ec)
3927	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
3928
3929	!IM
3930	IF (ip_ebil_phy.ge.1) THEN
3931	ztit='after physic'
3932	CALL diagetpq(cell_area,ztit,ip_ebil_phy,1,1,pdtphys &
3933	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3934	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3935	! Comme les tendances de la physique sont ajoute dans la dynamique,
3936	! on devrait avoir que la variation d'entalpie par la dynamique
3937	! est egale a la variation de la physique au pas de temps precedent.
3938	! Donc la somme de ces 2 variations devrait etre nulle.
3939
3940	call diagphy(cell_area,ztit,ip_ebil_phy &
3941	, topsw, toplw, solsw, sollw, sens &
3942	, evap, rain_fall, snow_fall, ztsol &
3943	, d_h_vcol, d_qt, d_ec &
3944	, fs_bound, fq_bound )
3945	!
3946	d_h_vcol_phy=d_h_vcol
3947	!
3948	END IF
3949	!
3950	!=======================================================================
3951	! SORTIES
3952	!=======================================================================
3953	!
3954	!IM initialisation + calculs divers diag AMIP2
3955	!
3956	include "calcul_divers.h"
3957	!
3958	!IM Interpolation sur les niveaux de pression du NMC
3959	! -------------------------------------------------
3960	!
3961	include "calcul_STDlev.h"
3962	!
3963	! slp sea level pressure
3964	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3965	!
3966	!cc prw = eau precipitable
3967	DO i = 1, klon
3968	prw(i) = 0.
3969	DO k = 1, klev
3970	prw(i) = prw(i) + &
3971	q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3972	ENDDO
3973	ENDDO
3974	!
3975	IF (type_trac == 'inca') THEN
3976	#ifdef INCA
3977	CALL VTe(VTphysiq)
3978	CALL VTb(VTinca)
3979
3980	CALL chemhook_end ( &
3981	pdtphys, &
3982	pplay, &
3983	t_seri, &
3984	tr_seri, &
3985	nbtr, &
3986	paprs, &
3987	q_seri, &
3988	cell_area, &
3989	pphi, &
3990	pphis, &
3991	zx_rh)
3992
3993	CALL VTe(VTinca)
3994	CALL VTb(VTphysiq)
3995	#endif
3996	END IF
3997
3998
3999	!
4000	! Convertir les incrementations en tendances
4001	!
4002	IF (prt_level .GE.10) THEN
4003	print *,'Convertir les incrementations en tendances '
4004	ENDIF
4005	!
4006	if (mydebug) then
4007	call writefield_phy('u_seri',u_seri,nbp_lev)
4008	call writefield_phy('v_seri',v_seri,nbp_lev)
4009	call writefield_phy('t_seri',t_seri,nbp_lev)
4010	call writefield_phy('q_seri',q_seri,nbp_lev)
4011	endif
4012
4013	DO k = 1, klev
4014	DO i = 1, klon
4015	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / pdtphys
4016	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / pdtphys
4017	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / pdtphys
4018	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / pdtphys
4019	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / pdtphys
4020	!CR: on ajoute le contenu en glace
4021	if (nqo.eq.3) then
4022	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / pdtphys
4023	endif
4024	ENDDO
4025	ENDDO
4026	!
4027	!CR: nb de traceurs eau: nqo
4028	! IF (nqtot.GE.3) THEN
4029	IF (nqtot.GE.(nqo+1)) THEN
4030	! DO iq = 3, nqtot
4031	DO iq = nqo+1, nqtot
4032	DO k = 1, klev
4033	DO i = 1, klon
4034	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / pdtphys
4035	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / pdtphys
4036	ENDDO
4037	ENDDO
4038	ENDDO
4039	ENDIF
4040	!
4041	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
4042	!IM global posePB include "write_bilKP_ins.h"
4043	!IM global posePB include "write_bilKP_ave.h"
4044	!
4045
4046	! Sauvegarder les valeurs de t et q a la fin de la physique:
4047	!
4048	DO k = 1, klev
4049	DO i = 1, klon
4050	u_ancien(i,k) = u_seri(i,k)
4051	v_ancien(i,k) = v_seri(i,k)
4052	t_ancien(i,k) = t_seri(i,k)
4053	q_ancien(i,k) = q_seri(i,k)
4054	ENDDO
4055	ENDDO
4056
4057	!!! RomP >>>
4058	!CR: nb de traceurs eau: nqo
4059	! IF (nqtot.GE.3) THEN
4060	IF (nqtot.GE.(nqo+1)) THEN
4061	! DO iq = 3, nqtot
4062	DO iq = nqo+1, nqtot
4063	DO k = 1, klev
4064	DO i = 1, klon
4065	! tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
4066	tr_ancien(i,k,iq-nqo) = tr_seri(i,k,iq-nqo)
4067	ENDDO
4068	ENDDO
4069	ENDDO
4070	ENDIF
4071	!!! RomP <<<
4072	!==========================================================================
4073	! Sorties des tendances pour un point particulier
4074	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
4075	! pour le debug
4076	! La valeur de igout est attribuee plus haut dans le programme
4077	!==========================================================================
4078
4079	if (prt_level.ge.1) then
4080	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
4081	write(lunout,*) &
4082	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
4083	write(lunout,*) &
4084	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
4085	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
4086	pctsrf(igout,is_sic)
4087	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
4088	do k=1,klev
4089	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
4090	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
4091	d_t_eva(igout,k)
4092	enddo
4093	write(lunout,*) 'cool,heat'
4094	do k=1,klev
4095	write(lunout,*) cool(igout,k),heat(igout,k)
4096	enddo
4097
4098	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
4099	do k=1,klev
4100	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
4101	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
4102	enddo
4103
4104	write(lunout,*) 'd_ps ',d_ps(igout)
4105	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
4106	do k=1,klev
4107	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
4108	d_qx(igout,k,1),d_qx(igout,k,2)
4109	enddo
4110	endif
4111
4112	!==========================================================================
4113
4114	!============================================================
4115	! Calcul de la temperature potentielle
4116	!============================================================
4117	DO k = 1, klev
4118	DO i = 1, klon
4119	!JYG/IM theta en debut du pas de temps
4120	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
4121	!JYG/IM theta en fin de pas de temps de physique
4122	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
4123	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers MPL 20130625
4124	! fth_fonctions.F90 et parkind1.F90
4125	! sinon thetal=theta
4126	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
4127	! : ql_seri(i,k))
4128	thetal(i,k)=theta(i,k)
4129	ENDDO
4130	ENDDO
4131	!
4132
4133	! 22.03.04 BEG
4134	!=============================================================
4135	! Ecriture des sorties
4136	!=============================================================
4137	#ifdef CPP_IOIPSL
4138
4139	! Recupere des varibles calcule dans differents modules
4140	! pour ecriture dans histxxx.nc
4141
4142	! Get some variables from module fonte_neige_mod
4143	CALL fonte_neige_get_vars(pctsrf, &
4144	zxfqcalving, zxfqfonte, zxffonte)
4145
4146
4147
4148
4149	!=============================================================
4150	! Separation entre thermiques et non thermiques dans les sorties
4151	! de fisrtilp
4152	!=============================================================
4153
4154	if (iflag_thermals>=1) then
4155	d_t_lscth=0.
4156	d_t_lscst=0.
4157	d_q_lscth=0.
4158	d_q_lscst=0.
4159	do k=1,klev
4160	do i=1,klon
4161	if (ptconvth(i,k)) then
4162	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4163	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4164	else
4165	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
4166	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
4167	endif
4168	enddo
4169	enddo
4170
4171	do i=1,klon
4172	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
4173	plul_th(i)=prfl(i,1)+psfl(i,1)
4174	enddo
4175	endif
4176
4177
4178	!On effectue les sorties:
4179
4180	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
4181	pplay, lmax_th, aerosol_couple, &
4182	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
4183	ptconv, read_climoz, clevSTD, &
4184	ptconvth, d_t, qx, d_qx, zmasse, &
4185	flag_aerosol, flag_aerosol_strat, ok_cdnc)
4186
4187
4188
4189	include "write_histday_seri.h"
4190
4191	include "write_paramLMDZ_phy.h"
4192
4193	#endif
4194
4195
4196	!====================================================================
4197	! Arret du modele apres hgardfou en cas de detection d'un
4198	! plantage par hgardfou
4199	!====================================================================
4200
4201	IF (abortphy==1) THEN
4202	abort_message ='Plantage hgardfou'
4203	CALL abort_physic (modname,abort_message,1)
4204	ENDIF
4205
4206
4207	! 22.03.04 END
4208	!
4209	!====================================================================
4210	! Si c'est la fin, il faut conserver l'etat de redemarrage
4211	!====================================================================
4212	!
4213
4214	! -----------------------------------------------------------------
4215	! WSTATS: Saving statistics
4216	! -----------------------------------------------------------------
4217	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
4218	! which can later be used to make the statistic files of the run:
4219	! "stats") only possible in 3D runs !
4220
4221	#ifdef YM_TO_DO_LATER
4222	IF (callstats) THEN
4223
4224	call wstats(klon,o_psol%name,"Surface pressure","Pa" &
4225	,2,paprs(:,1))
4226	call wstats(klon,o_tsol%name,"Surface temperature","K", &
4227	2,zxtsol)
4228	zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
4229	call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
4230	"kg/(s*m2)",2,zx_tmp_fi2d)
4231	zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
4232	call wstats(klon,o_plul%name,"Large-scale Precip", &
4233	"kg/(s*m2)",2,zx_tmp_fi2d)
4234	zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
4235	call wstats(klon,o_pluc%name,"Convective Precip", &
4236	"kg/(s*m2)",2,zx_tmp_fi2d)
4237	call wstats(klon,o_sols%name,"Solar rad. at surf.", &
4238	"W/m2",2,solsw)
4239	call wstats(klon,o_soll%name,"IR rad. at surf.", &
4240	"W/m2",2,sollw)
4241	zx_tmp_fi2d(:) = topsw(:)-toplw(:)
4242	call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
4243	"W/m2",2,zx_tmp_fi2d)
4244
4245
4246
4247	call wstats(klon,o_temp%name,"Air temperature","K", &
4248	3,t_seri)
4249	call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
4250	3,u_seri)
4251	call wstats(klon,o_vitv%name,"Meridional wind", &
4252	"m.s-1",3,v_seri)
4253	call wstats(klon,o_vitw%name,"Vertical wind", &
4254	"m.s-1",3,omega)
4255	call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
4256	3,q_seri)
4257
4258
4259
4260	IF(lafin) THEN
4261	write (,) "Writing stats..."
4262	call mkstats(ierr)
4263	ENDIF
4264
4265	ENDIF !if callstats
4266	#endif
4267
4268	IF (lafin) THEN
4269	itau_phy = itau_phy + itap
4270	CALL phyredem ("restartphy.nc")
4271	! open(97,form="unformatted",file="finbin")
4272	! write(97) u_seri,v_seri,t_seri,q_seri
4273	! close(97)
4274	!$OMP MASTER
4275	if (read_climoz >= 1) then
4276	if (is_mpi_root) then
4277	call nf95_close(ncid_climoz)
4278	end if
4279	deallocate(press_climoz) ! pointer
4280	end if
4281	!$OMP END MASTER
4282	ENDIF
4283
4284	! first=.false.
4285
4286	RETURN
4287	END SUBROUTINE physiq
4288	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
4289	IMPLICIT none
4290	!
4291	! Calculer et imprimer l'eau totale. A utiliser pour verifier
4292	! la conservation de l'eau
4293	!
4294	include "YOMCST.h"
4295	INTEGER klon,klev
4296	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
4297	REAL aire(klon)
4298	REAL qtotal, zx, qcheck
4299	INTEGER i, k
4300	!
4301	zx = 0.0
4302	DO i = 1, klon
4303	zx = zx + aire(i)
4304	ENDDO
4305	qtotal = 0.0
4306	DO k = 1, klev
4307	DO i = 1, klon
4308	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
4309	*(paprs(i,k)-paprs(i,k+1))/RG
4310	ENDDO
4311	ENDDO
4312	!
4313	qcheck = qtotal/zx
4314	!
4315	RETURN
4316	END FUNCTION qcheck

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

Download in other formats:

Original Format