Context Navigation

physiq.F90 @ 2363

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

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

Download in other formats:

Original Format