Context Navigation

physiq.F90 @ 2399

Last change on this file since 2399 was 2399, checked in by Ehouarn Millour, 9 years ago
Follow-up from commit 2395: get rid of rlon and rlat, longitude_deg and latitude_deg (from module geometry_mod) should be used instead. Longitudes and latitudes are no longer loaded from startphy.nc but inherited from dynamics (and compatibility with values in startphy.nc is checked). This will change bench results because of roundoffs differences between the two. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 145.4 KB

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

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

Download in other formats:

Original Format