Context Navigation

physiq_mod.F90 @ 2475

Last change on this file since 2475 was 2469, checked in by lguez, 9 years ago
Indented physis_mod.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 155.9 KB

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

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

Download in other formats:

Original Format