Context Navigation

physiq_mod.F90 @ 2443

Last change on this file since 2443 was 2443, checked in by acozic, 8 years ago
add some validation tests for aerosols parameters (config_inca, aerosol_couple, flag_aerosol)
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: 144.5 KB

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

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

Download in other formats:

Original Format