Context Navigation

physiq.F90 @ 2372

Last change on this file since 2372 was 2372, checked in by acozic, 9 years ago
Change call to inca initialisation to fit with new sections dynamique/physic
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.7 KB

Line
1	! $Id: physiq.F90 2372 2015-10-13 12:41:47Z 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)
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
1063	config_inca='none' ! default
1064	CALL getin_p('config_inca',config_inca)
1065
1066	END IF
1067
1068	rnebcon0(:,:) = 0.0
1069	clwcon0(:,:) = 0.0
1070	rnebcon(:,:) = 0.0
1071	clwcon(:,:) = 0.0
1072
1073	!IM
1074	IF (ip_ebil_phy.ge.1) d_h_vcol_phy=0.
1075	!
1076	print*,'iflag_coupl,iflag_clos,iflag_wake', &
1077	iflag_coupl,iflag_clos,iflag_wake
1078	print*,'iflag_CYCLE_DIURNE', iflag_cycle_diurne
1079	!
1080	IF (iflag_con.EQ.2.AND.iflag_cld_th.GT.-1) THEN
1081	abort_message = 'Tiedtke needs iflag_cld_th=-2 or -1'
1082	CALL abort_physic (modname,abort_message,1)
1083	ENDIF
1084	!
1085	!
1086	! Initialiser les compteurs:
1087	!
1088	itap = 0
1089	itaprad = 0
1090
1091	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1092	!! Un petit travail \`a faire ici.
1093	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1094
1095	if (iflag_pbl>1) then
1096	PRINT*, "Using method MELLOR&YAMADA"
1097	endif
1098
1099	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1100	! FH 2008/05/02 changement lie a la lecture de nbapp_rad dans phylmd plutot que
1101	! dyn3d
1102	! Attention : la version precedente n'etait pas tres propre.
1103	! Il se peut qu'il faille prendre une valeur differente de nbapp_rad
1104	! pour obtenir le meme resultat.
1105	dtime=pdtphys
1106	IF (MOD(INT(86400./dtime),nbapp_rad).EQ.0) THEN
1107	radpas = NINT( 86400./dtime/nbapp_rad)
1108	ELSE
1109	WRITE(lunout,*) 'le nombre de pas de temps physique doit etre un multiple de nbapp_rad'
1110	WRITE(lunout,*) 'changer nbapp_rad ou alors commenter ce test mais 1+1<>2'
1111	abort_message='nbre de pas de temps physique n est pas multiple de nbapp_rad'
1112	call abort_physic(modname,abort_message,1)
1113	ENDIF
1114	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1115
1116	CALL phyetat0 ("startphy.nc",clesphy0,tabcntr0)
1117	IF (klon_glo==1) THEN
1118	coefh=0. ; coefm=0. ; pbl_tke=0.
1119	coefh(:,2,:)=1.e-2 ; coefm(:,2,:)=1.e-2 ; pbl_tke(:,2,:)=1.e-2
1120	PRINT*,'FH WARNING : lignes a supprimer'
1121	ENDIF
1122	!IM begin
1123	print*,'physiq: clwcon rnebcon ratqs',clwcon(1,1),rnebcon(1,1) &
1124	,ratqs(1,1)
1125	!IM end
1126
1127
1128
1129	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1130	!
1131	! on remet le calendrier a zero
1132	!
1133	IF (raz_date .eq. 1) THEN
1134	itau_phy = 0
1135	ENDIF
1136
1137	CALL printflag( tabcntr0,radpas,ok_journe, &
1138	ok_instan, ok_region )
1139	!
1140	IF (ABS(dtime-pdtphys).GT.0.001) THEN
1141	WRITE(lunout,*) 'Pas physique n est pas correct',dtime, &
1142	pdtphys
1143	abort_message='Pas physique n est pas correct '
1144	! call abort_physic(modname,abort_message,1)
1145	dtime=pdtphys
1146	ENDIF
1147	IF (nlon .NE. klon) THEN
1148	WRITE(lunout,*)'nlon et klon ne sont pas coherents', nlon, &
1149	klon
1150	abort_message='nlon et klon ne sont pas coherents'
1151	call abort_physic(modname,abort_message,1)
1152	ENDIF
1153	IF (nlev .NE. klev) THEN
1154	WRITE(lunout,*)'nlev et klev ne sont pas coherents', nlev, &
1155	klev
1156	abort_message='nlev et klev ne sont pas coherents'
1157	call abort_physic(modname,abort_message,1)
1158	ENDIF
1159	!
1160	IF (dtime*REAL(radpas).GT.21600..AND.iflag_cycle_diurne.GE.1) THEN
1161	WRITE(lunout,*)'Nbre d appels au rayonnement insuffisant'
1162	WRITE(lunout,*)"Au minimum 4 appels par jour si cycle diurne"
1163	abort_message='Nbre d appels au rayonnement insuffisant'
1164	call abort_physic(modname,abort_message,1)
1165	ENDIF
1166	WRITE(lunout,*)"Clef pour la convection, iflag_con=", iflag_con
1167	WRITE(lunout,*)"Clef pour le driver de la convection, ok_cvl=", &
1168	ok_cvl
1169	!
1170	!KE43
1171	! Initialisation pour la convection de K.E. (sb):
1172	IF (iflag_con.GE.3) THEN
1173
1174	WRITE(lunout,)"** Convection de Kerry Emanuel 4.3 "
1175	WRITE(lunout,*) &
1176	"On va utiliser le melange convectif des traceurs qui"
1177	WRITE(lunout,*)"est calcule dans convect4.3"
1178	WRITE(lunout,*)" !!! penser aux logical flags de phytrac"
1179
1180	DO i = 1, klon
1181	ema_cbmf(i) = 0.
1182	ema_pcb(i) = 0.
1183	ema_pct(i) = 0.
1184	! ema_workcbmf(i) = 0.
1185	ENDDO
1186	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
1187	DO i = 1, klon
1188	ibas_con(i) = 1
1189	itop_con(i) = 1
1190	ENDDO
1191	!IM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
1192	!===============================================================================
1193	!CR:04.12.07: initialisations poches froides
1194	! Controle de ALE et ALP pour la fermeture convective (jyg)
1195	if (iflag_wake>=1) then
1196	CALL ini_wake(0.,0.,it_wape_prescr,wape_prescr,fip_prescr &
1197	,alp_bl_prescr, ale_bl_prescr)
1198	! 11/09/06 rajout initialisation ALE et ALP du wake et PBL(YU)
1199	! print*,'apres ini_wake iflag_cld_th=', iflag_cld_th
1200	endif
1201
1202	! do i = 1,klon
1203	! Ale_bl(i)=0.
1204	! Alp_bl(i)=0.
1205	! enddo
1206
1207	!================================================================================
1208	!IM stations CFMIP
1209	nCFMIP=npCFMIP
1210	OPEN(98,file='npCFMIP_param.data',status='old', &
1211	form='formatted',iostat=iostat)
1212	if (iostat == 0) then
1213	READ(98,*,end=998) nCFMIP
1214	998 CONTINUE
1215	CLOSE(98)
1216	CONTINUE
1217	IF(nCFMIP.GT.npCFMIP) THEN
1218	print*,'nCFMIP > npCFMIP : augmenter npCFMIP et recompiler'
1219	call abort_physic("physiq", "", 1)
1220	else
1221	print*,'physiq npCFMIP=',npCFMIP,'nCFMIP=',nCFMIP
1222	ENDIF
1223
1224	!
1225	ALLOCATE(tabCFMIP(nCFMIP))
1226	ALLOCATE(lonCFMIP(nCFMIP), latCFMIP(nCFMIP))
1227	ALLOCATE(tabijGCM(nCFMIP))
1228	ALLOCATE(lonGCM(nCFMIP), latGCM(nCFMIP))
1229	ALLOCATE(iGCM(nCFMIP), jGCM(nCFMIP))
1230	!
1231	! lecture des nCFMIP stations CFMIP, de leur numero
1232	! et des coordonnees geographiques lonCFMIP, latCFMIP
1233	!
1234	CALL read_CFMIP_point_locations(nCFMIP, tabCFMIP, &
1235	lonCFMIP, latCFMIP)
1236	!
1237	! identification des
1238	! 1) coordonnees lonGCM, latGCM des points CFMIP dans la grille de LMDZ
1239	! 2) indices points tabijGCM de la grille physique 1d sur klon points
1240	! 3) indices iGCM, jGCM de la grille physique 2d
1241	!
1242	CALL LMDZ_CFMIP_point_locations(nCFMIP, lonCFMIP, latCFMIP, &
1243	tabijGCM, lonGCM, latGCM, iGCM, jGCM)
1244	!
1245	else
1246	ALLOCATE(tabijGCM(0))
1247	ALLOCATE(lonGCM(0), latGCM(0))
1248	ALLOCATE(iGCM(0), jGCM(0))
1249	end if
1250	else
1251	ALLOCATE(tabijGCM(0))
1252	ALLOCATE(lonGCM(0), latGCM(0))
1253	ALLOCATE(iGCM(0), jGCM(0))
1254	ENDIF
1255
1256	DO i=1,klon
1257	rugoro(i) = f_rugoro * MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1258	ENDDO
1259
1260	!34EK
1261	IF (ok_orodr) THEN
1262
1263	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1264	! FH sans doute a enlever de finitivement ou, si on le garde, l'activer
1265	! justement quand ok_orodr = false.
1266	! ce rugoro est utilise par la couche limite et fait double emploi
1267	! avec les param\'etrisations sp\'ecifiques de Francois Lott.
1268	! DO i=1,klon
1269	! rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
1270	! ENDDO
1271	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1272	IF (ok_strato) THEN
1273	CALL SUGWD_strato(klon,klev,paprs,pplay)
1274	ELSE
1275	CALL SUGWD(klon,klev,paprs,pplay)
1276	ENDIF
1277
1278	DO i=1,klon
1279	zuthe(i)=0.
1280	zvthe(i)=0.
1281	if(zstd(i).gt.10.)then
1282	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
1283	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
1284	endif
1285	ENDDO
1286	ENDIF
1287	!
1288	!
1289	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
1290	WRITE(lunout,*)'La frequence de lecture surface est de ', &
1291	lmt_pas
1292	!
1293	capemaxcels = 't_max(X)'
1294	t2mincels = 't_min(X)'
1295	t2maxcels = 't_max(X)'
1296	tinst = 'inst(X)'
1297	tave = 'ave(X)'
1298	!IM cf. AM 081204 BEG
1299	write(lunout,*)'AVANT HIST IFLAG_CON=',iflag_con
1300	!IM cf. AM 081204 END
1301	!
1302	!=============================================================
1303	! Initialisation des sorties
1304	!=============================================================
1305
1306	#ifdef CPP_IOIPSL
1307
1308	!$OMP MASTER
1309	! FH : if ok_sync=.true. , the time axis is written at each time step
1310	! in the output files. Only at the end in the opposite case
1311	ok_sync_omp=.false.
1312	CALL getin('ok_sync',ok_sync_omp)
1313	call phys_output_open(rlon,rlat,nCFMIP,tabijGCM, &
1314	iGCM,jGCM,lonGCM,latGCM, &
1315	jjmp1,nlevSTD,clevSTD,rlevSTD, dtime,ok_veget, &
1316	type_ocean,iflag_pbl,iflag_pbl_split,ok_mensuel,ok_journe, &
1317	ok_hf,ok_instan,ok_LES,ok_ade,ok_aie, &
1318	read_climoz, phys_out_filestations, &
1319	new_aod, aerosol_couple, &
1320	flag_aerosol_strat, pdtphys, paprs, pphis, &
1321	pplay, lmax_th, ptconv, ptconvth, ivap, &
1322	d_t, qx, d_qx, zmasse, ok_sync_omp)
1323	!$OMP END MASTER
1324	!$OMP BARRIER
1325	ok_sync=ok_sync_omp
1326
1327	freq_outNMC(1) = ecrit_files(7)
1328	freq_outNMC(2) = ecrit_files(8)
1329	freq_outNMC(3) = ecrit_files(9)
1330	WRITE(lunout,*)'OK freq_outNMC(1)=',freq_outNMC(1)
1331	WRITE(lunout,*)'OK freq_outNMC(2)=',freq_outNMC(2)
1332	WRITE(lunout,*)'OK freq_outNMC(3)=',freq_outNMC(3)
1333
1334	include "ini_histday_seri.h"
1335
1336	include "ini_paramLMDZ_phy.h"
1337
1338	#endif
1339	ecrit_reg = ecrit_reg * un_jour
1340	ecrit_tra = ecrit_tra * un_jour
1341
1342	!XXXPB Positionner date0 pour initialisation de ORCHIDEE
1343	date0 = jD_ref
1344	WRITE(,) 'physiq date0 : ',date0
1345	!
1346	!
1347	!
1348	! Prescrire l'ozone dans l'atmosphere
1349	!
1350	!
1351	!c DO i = 1, klon
1352	!c DO k = 1, klev
1353	!c CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1354	!c ENDDO
1355	!c ENDDO
1356	!
1357	IF (type_trac == 'inca') THEN
1358	#ifdef INCA
1359	CALL VTe(VTphysiq)
1360	CALL VTb(VTinca)
1361	calday = REAL(days_elapsed) + jH_cur
1362	WRITE(lunout,*) 'initial time chemini', days_elapsed, calday
1363
1364	CALL chemini( &
1365	rg, &
1366	ra, &
1367	cell_area, &
1368	rlat, &
1369	rlon, &
1370	presnivs, &
1371	calday, &
1372	klon, &
1373	nqtot, &
1374	pdtphys, &
1375	annee_ref, &
1376	day_ref, &
1377	day_ini, &
1378	start_time, &
1379	itau_phy, &
1380	io_lon, &
1381	io_lat)
1382
1383	CALL VTe(VTinca)
1384	CALL VTb(VTphysiq)
1385	#endif
1386	END IF
1387	!
1388	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1389	! Nouvelle initialisation pour le rayonnement RRTM
1390	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1391
1392	call iniradia(klon,klev,paprs(1,1:klev+1))
1393
1394	!$omp single
1395	if (read_climoz >= 1) then
1396	call open_climoz(ncid_climoz, press_climoz)
1397	END IF
1398	!$omp end single
1399	!
1400	!IM betaCRF
1401	pfree=70000. !Pa
1402	beta_pbl=1.
1403	beta_free=1.
1404	lon1_beta=-180.
1405	lon2_beta=+180.
1406	lat1_beta=90.
1407	lat2_beta=-90.
1408	mskocean_beta=.FALSE.
1409
1410	!albedo SB >>>
1411	select case(nsw)
1412	case(2)
1413	SFRWL(1)=0.45538747
1414	SFRWL(2)=0.54461211
1415	case(4)
1416	SFRWL(1)=0.45538747
1417	SFRWL(2)=0.32870591
1418	SFRWL(3)=0.18568763
1419	SFRWL(4)=3.02191470E-02
1420	case(6)
1421	SFRWL(1)=1.28432794E-03
1422	SFRWL(2)=0.12304168
1423	SFRWL(3)=0.33106142
1424	SFRWL(4)=0.32870591
1425	SFRWL(5)=0.18568763
1426	SFRWL(6)=3.02191470E-02
1427	end select
1428
1429
1430	!albedo SB <<<
1431
1432	OPEN(99,file='beta_crf.data',status='old', &
1433	form='formatted',err=9999)
1434	READ(99,*,end=9998) pfree
1435	READ(99,*,end=9998) beta_pbl
1436	READ(99,*,end=9998) beta_free
1437	READ(99,*,end=9998) lon1_beta
1438	READ(99,*,end=9998) lon2_beta
1439	READ(99,*,end=9998) lat1_beta
1440	READ(99,*,end=9998) lat2_beta
1441	READ(99,*,end=9998) mskocean_beta
1442	9998 Continue
1443	CLOSE(99)
1444	9999 Continue
1445	WRITE(,)'pfree=',pfree
1446	WRITE(,)'beta_pbl=',beta_pbl
1447	WRITE(,)'beta_free=',beta_free
1448	WRITE(,)'lon1_beta=',lon1_beta
1449	WRITE(,)'lon2_beta=',lon2_beta
1450	WRITE(,)'lat1_beta=',lat1_beta
1451	WRITE(,)'lat2_beta=',lat2_beta
1452	WRITE(,)'mskocean_beta=',mskocean_beta
1453	ENDIF
1454	!
1455	! ************** Fin de IF ( debut ) *************
1456	!
1457	!
1458	! Incrementer le compteur de la physique
1459	!
1460	itap = itap + 1
1461	!
1462	!
1463	! Update fraction of the sub-surfaces (pctsrf) and
1464	! initialize, where a new fraction has appeared, all variables depending
1465	! on the surface fraction.
1466	!
1467	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
1468	pctsrf, fevap, z0m, z0h, agesno, &
1469	falb_dir, falb_dif, ftsol, ustar, u10m, v10m, pbl_tke)
1470
1471	! Update time and other variables in Reprobus
1472	IF (type_trac == 'repr') THEN
1473	#ifdef REPROBUS
1474	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1475	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1476	CALL Rtime(debut)
1477	#endif
1478	END IF
1479
1480
1481	! Tendances bidons pour les processus qui n'affectent pas certaines
1482	! variables.
1483	du0(:,:)=0.
1484	dv0(:,:)=0.
1485	dt0 = 0.
1486	dq0(:,:)=0.
1487	dql0(:,:)=0.
1488	dqi0(:,:)=0.
1489	!
1490	! Mettre a zero des variables de sortie (pour securite)
1491	!
1492	DO i = 1, klon
1493	d_ps(i) = 0.0
1494	ENDDO
1495	DO k = 1, klev
1496	DO i = 1, klon
1497	d_t(i,k) = 0.0
1498	d_u(i,k) = 0.0
1499	d_v(i,k) = 0.0
1500	ENDDO
1501	ENDDO
1502	DO iq = 1, nqtot
1503	DO k = 1, klev
1504	DO i = 1, klon
1505	d_qx(i,k,iq) = 0.0
1506	ENDDO
1507	ENDDO
1508	ENDDO
1509	da(:,:)=0.
1510	mp(:,:)=0.
1511	phi(:,:,:)=0.
1512	! RomP >>>
1513	phi2(:,:,:)=0.
1514	beta_prec_fisrt(:,:)=0.
1515	beta_prec(:,:)=0.
1516	epmlmMm(:,:,:)=0.
1517	eplaMm(:,:)=0.
1518	d1a(:,:)=0.
1519	dam(:,:)=0.
1520	pmflxr=0.
1521	pmflxs=0.
1522	! RomP <<<
1523
1524	!
1525	! Ne pas affecter les valeurs entrees de u, v, h, et q
1526	!
1527	DO k = 1, klev
1528	DO i = 1, klon
1529	t_seri(i,k) = t(i,k)
1530	u_seri(i,k) = u(i,k)
1531	v_seri(i,k) = v(i,k)
1532	q_seri(i,k) = qx(i,k,ivap)
1533	ql_seri(i,k) = qx(i,k,iliq)
1534	!CR: ATTENTION, on rajoute la variable glace
1535	if (nqo.eq.2) then
1536	qs_seri(i,k) = 0.
1537	else if (nqo.eq.3) then
1538	qs_seri(i,k) = qx(i,k,isol)
1539	endif
1540	ENDDO
1541	ENDDO
1542	tke0(:,:)=pbl_tke(:,:,is_ave)
1543	!CR:Nombre de traceurs de l'eau: nqo
1544	! IF (nqtot.GE.3) THEN
1545	IF (nqtot.GE.(nqo+1)) THEN
1546	! DO iq = 3, nqtot
1547	DO iq = nqo+1, nqtot
1548	DO k = 1, klev
1549	DO i = 1, klon
1550	! tr_seri(i,k,iq-2) = qx(i,k,iq)
1551	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
1552	ENDDO
1553	ENDDO
1554	ENDDO
1555	ELSE
1556	DO k = 1, klev
1557	DO i = 1, klon
1558	tr_seri(i,k,1) = 0.0
1559	ENDDO
1560	ENDDO
1561	ENDIF
1562	!
1563	DO i = 1, klon
1564	ztsol(i) = 0.
1565	ENDDO
1566	DO nsrf = 1, nbsrf
1567	DO i = 1, klon
1568	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1569	ENDDO
1570	ENDDO
1571	!IM
1572	IF (ip_ebil_phy.ge.1) THEN
1573	ztit='after dynamic'
1574	CALL diagetpq(cell_area,ztit,ip_ebil_phy,1,1,dtime &
1575	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1576	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1577	! Comme les tendances de la physique sont ajoute dans la dynamique,
1578	! on devrait avoir que la variation d'entalpie par la dynamique
1579	! est egale a la variation de la physique au pas de temps precedent.
1580	! Donc la somme de ces 2 variations devrait etre nulle.
1581	call diagphy(cell_area,ztit,ip_ebil_phy &
1582	, zero_v, zero_v, zero_v, zero_v, zero_v &
1583	, zero_v, zero_v, zero_v, ztsol &
1584	, d_h_vcol+d_h_vcol_phy, d_qt, 0. &
1585	, fs_bound, fq_bound )
1586	END IF
1587
1588	! Diagnostiquer la tendance dynamique
1589	!
1590	IF (ancien_ok) THEN
1591	DO k = 1, klev
1592	DO i = 1, klon
1593	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1594	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1595	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1596	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1597	ENDDO
1598	ENDDO
1599	!!! RomP >>> td dyn traceur
1600	!! IF (nqtot.GE.3) THEN ! jyg
1601	!! DO iq = 3, nqtot ! jyg
1602	IF (nqtot.GE.nqo+1) THEN ! jyg
1603	DO iq = nqo+1, nqtot ! jyg
1604	DO k = 1, klev
1605	DO i = 1, klon
1606	!! d_tr_dyn(i,k,iq-2)= & ! jyg
1607	!! (tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime ! jyg
1608	d_tr_dyn(i,k,iq-nqo)= & ! jyg
1609	(tr_seri(i,k,iq-nqo)-tr_ancien(i,k,iq-nqo))/dtime ! jyg
1610	! iiq=niadv(iq)
1611	! print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-nqo),"tra:",iq,tname(iiq)
1612	ENDDO
1613	ENDDO
1614	ENDDO
1615	ENDIF
1616	!!! RomP <<<
1617	ELSE
1618	DO k = 1, klev
1619	DO i = 1, klon
1620	d_u_dyn(i,k) = 0.0
1621	d_v_dyn(i,k) = 0.0
1622	d_t_dyn(i,k) = 0.0
1623	d_q_dyn(i,k) = 0.0
1624	ENDDO
1625	ENDDO
1626	!!! RomP >>> td dyn traceur
1627	!! IF (nqtot.GE.3) THEN ! jyg
1628	!! DO iq = 3, nqtot ! jyg
1629	IF (nqtot.GE.nqo+1) THEN ! jyg
1630	DO iq = nqo+1, nqtot ! jyg
1631	DO k = 1, klev
1632	DO i = 1, klon
1633	!! d_tr_dyn(i,k,iq-2)= 0.0 ! jyg
1634	d_tr_dyn(i,k,iq-nqo)= 0.0 ! jyg
1635	ENDDO
1636	ENDDO
1637	ENDDO
1638	ENDIF
1639	!!! RomP <<<
1640	ancien_ok = .TRUE.
1641	ENDIF
1642	!
1643	! Ajouter le geopotentiel du sol:
1644	!
1645	DO k = 1, klev
1646	DO i = 1, klon
1647	zphi(i,k) = pphi(i,k) + pphis(i)
1648	ENDDO
1649	ENDDO
1650	!
1651	! Verifier les temperatures
1652	!
1653	!IM BEG
1654	IF (check) THEN
1655	amn=MIN(ftsol(1,is_ter),1000.)
1656	amx=MAX(ftsol(1,is_ter),-1000.)
1657	DO i=2, klon
1658	amn=MIN(ftsol(i,is_ter),amn)
1659	amx=MAX(ftsol(i,is_ter),amx)
1660	ENDDO
1661	!
1662	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1663	ENDIF !(check) THEN
1664	!IM END
1665	!
1666	CALL hgardfou(t_seri,ftsol,'debutphy',abortphy)
1667	IF (abortphy==1) Print*,'ERROR ABORT hgardfou debutphy'
1668
1669	!
1670	!IM BEG
1671	IF (check) THEN
1672	amn=MIN(ftsol(1,is_ter),1000.)
1673	amx=MAX(ftsol(1,is_ter),-1000.)
1674	DO i=2, klon
1675	amn=MIN(ftsol(i,is_ter),amn)
1676	amx=MAX(ftsol(i,is_ter),amx)
1677	ENDDO
1678	!
1679	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1680	ENDIF !(check) THEN
1681	!IM END
1682	!
1683	! Mettre en action les conditions aux limites (albedo, sst, etc.).
1684	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
1685	!
1686	if (read_climoz >= 1) then
1687	! Ozone from a file
1688	! Update required ozone index:
1689	ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
1690	if (ro3i == 361) ro3i = 360
1691	! (This should never occur, except perhaps because of roundup
1692	! error. See documentation.)
1693	if (ro3i /= co3i) then
1694	! Update ozone field:
1695	if (read_climoz == 1) then
1696	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
1697	press_in_edg=press_climoz, paprs=paprs, v3=wo)
1698	else
1699	! read_climoz == 2
1700	call regr_pr_av(ncid_climoz, (/"tro3 ", "tro3_daylight"/), &
1701	julien=ro3i, press_in_edg=press_climoz, paprs=paprs, v3=wo)
1702	end if
1703	! Convert from mole fraction of ozone to column density of ozone in a
1704	! cell, in kDU:
1705	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
1706	* zmasse / dobson_u / 1e3
1707	! (By regridding ozone values for LMDZ only once every 360th of
1708	! year, we have already neglected the variation of pressure in one
1709	! 360th of year. So do not recompute "wo" at each time step even if
1710	! "zmasse" changes a little.)
1711	co3i = ro3i
1712	end if
1713	ELSEIF (MOD(itap-1,lmt_pas) == 0) THEN
1714	! Once per day, update ozone from Royer:
1715
1716	IF (solarlong0<-999.) then
1717	! Generic case with evolvoing season
1718	zzz=real(days_elapsed+1)
1719	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
1720	! Particular case with annual mean insolation
1721	zzz=real(90) ! could be revisited
1722	IF (read_climoz/=-1) THEN
1723	abort_message ='read_climoz=-1 is recommended when solarlong0=1000.'
1724	CALL abort_physic (modname,abort_message,1)
1725	ENDIF
1726	ELSE
1727	! Case where the season is imposed with solarlong0
1728	zzz=real(90) ! could be revisited
1729	ENDIF
1730	wo(:,:,1)=ozonecm(rlat, paprs,read_climoz,rjour=zzz)
1731	ENDIF
1732	!
1733	! Re-evaporer l'eau liquide nuageuse
1734	!
1735	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1736	DO i = 1, klon
1737	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1738	!jyg<
1739	! Attention : Arnaud a propose des formules completement differentes
1740	! A verifier !!!
1741	zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1742	IF (iflag_ice_thermo .EQ. 0) THEN
1743	zlsdcp=zlvdcp
1744	ENDIF
1745	!>jyg
1746
1747	if (iflag_ice_thermo.eq.0) then
1748	!pas necessaire a priori
1749
1750	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1751	zb = MAX(0.0,ql_seri(i,k))
1752	za = - MAX(0.0,ql_seri(i,k)) &
1753	* (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1754	t_seri(i,k) = t_seri(i,k) + za
1755	q_seri(i,k) = q_seri(i,k) + zb
1756	ql_seri(i,k) = 0.0
1757	d_t_eva(i,k) = za
1758	d_q_eva(i,k) = zb
1759
1760	else
1761
1762	!CR: on rÃ©-Ã©vapore eau liquide et glace
1763
1764	! zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1765	! zb = MAX(0.0,ql_seri(i,k))
1766	! za = - MAX(0.0,ql_seri(i,k)) &
1767	! * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1768	zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
1769	za = - MAX(0.0,ql_seri(i,k))*zlvdcp &
1770	- MAX(0.0,qs_seri(i,k))*zlsdcp
1771	t_seri(i,k) = t_seri(i,k) + za
1772	q_seri(i,k) = q_seri(i,k) + zb
1773	ql_seri(i,k) = 0.0
1774	!on Ã©vapore la glace
1775	qs_seri(i,k) = 0.0
1776	d_t_eva(i,k) = za
1777	d_q_eva(i,k) = zb
1778	endif
1779
1780	ENDDO
1781	ENDDO
1782	!IM
1783	IF (ip_ebil_phy.ge.2) THEN
1784	ztit='after reevap'
1785	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,1,dtime &
1786	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1787	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1788	call diagphy(cell_area,ztit,ip_ebil_phy &
1789	, zero_v, zero_v, zero_v, zero_v, zero_v &
1790	, zero_v, zero_v, zero_v, ztsol &
1791	, d_h_vcol, d_qt, d_ec &
1792	, fs_bound, fq_bound )
1793	!
1794	END IF
1795
1796	!
1797	!=========================================================================
1798	! Calculs de l'orbite.
1799	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1800	! doit donc etre plac\'e avant radlwsw et pbl_surface
1801
1802	!!! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1803	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1804	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1805	!
1806	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1807	! solarlong0
1808	if (solarlong0<-999.) then
1809	if (new_orbit) then
1810	! calcul selon la routine utilisee pour les planetes
1811	call solarlong(day_since_equinox, zlongi, dist)
1812	else
1813	! calcul selon la routine utilisee pour l'AR4
1814	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1815	endif
1816	else
1817	zlongi=solarlong0 ! longitude solaire vraie
1818	dist=1. ! distance au soleil / moyenne
1819	endif
1820	if(prt_level.ge.1) &
1821	write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1822
1823
1824	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1825	! Calcul de l'ensoleillement :
1826	! ============================
1827	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
1828	! l'annee a partir d'une formule analytique.
1829	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
1830	! non nul aux poles.
1831	IF (abs(solarlong0-1000.)<1.e-4) then
1832	call zenang_an(iflag_cycle_diurne.GE.1,jH_cur,rlat,rlon,rmu0,fract)
1833	JrNt = 1.0
1834	ELSE
1835	! recode par Olivier Boucher en sept 2015
1836	SELECT CASE (iflag_cycle_diurne)
1837	CASE(0)
1838	! Sans cycle diurne
1839	CALL angle(zlongi, rlat, fract, rmu0)
1840	swradcorr = 1.0
1841	JrNt = 1.0
1842	zrmu0 = rmu0
1843	CASE(1)
1844	! Avec cycle diurne sans application des poids
1845	! bit comparable a l ancienne formulation cycle_diurne=true
1846	! on integre entre gmtime et gmtime+radpas
1847	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1848	CALL zenang(zlongi,jH_cur,0.0,zdtime,rlat,rlon,rmu0,fract)
1849	zrmu0 = rmu0
1850	swradcorr = 1.0
1851	! Calcul du flag jour-nuit
1852	JrNt = 0.0
1853	WHERE (fract.GT.0.0) JrNt = 1.0
1854	CASE(2)
1855	! Avec cycle diurne sans application des poids
1856	! On integre entre gmtime-pdtphys et gmtime+pdtphys*(radpas-1)
1857	! Comme cette routine est appele a tous les pas de temps de la physique
1858	! meme si le rayonnement n'est pas appele je remonte en arriere les
1859	! radpas-1 pas de temps suivant. Petite ruse avec MOD pour prendre en
1860	! compte le premier pas de temps de la physique pendant lequel itaprad=0
1861	zdtime1=dtime*REAL(-MOD(itaprad,radpas)-1)
1862	zdtime2=dtime*REAL(radpas-MOD(itaprad,radpas)-1)
1863	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2,rlat,rlon,rmu0,fract)
1864	!
1865	! Calcul des poids
1866	!
1867	zdtime1=-dtime !--on corrige le rayonnement pour representer le
1868	zdtime2=0.0 !--pas de temps de la physique qui se termine
1869	CALL zenang(zlongi,jH_cur,zdtime1,zdtime2,rlat,rlon,zrmu0,zfract)
1870	swradcorr = 0.0
1871	WHERE (rmu0.GE.1.e-10 .OR. fract.GE.1.e-10) swradcorr=zfract/fract*zrmu0/rmu0
1872	! Calcul du flag jour-nuit
1873	JrNt = 0.0
1874	WHERE (zfract.GT.0.0) JrNt = 1.0
1875	END SELECT
1876	ENDIF
1877
1878	if (mydebug) then
1879	call writefield_phy('u_seri',u_seri,nbp_lev)
1880	call writefield_phy('v_seri',v_seri,nbp_lev)
1881	call writefield_phy('t_seri',t_seri,nbp_lev)
1882	call writefield_phy('q_seri',q_seri,nbp_lev)
1883	endif
1884
1885	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1886	! Appel au pbl_surface : Planetary Boudary Layer et Surface
1887	! Cela implique tous les interactions des sous-surfaces et la partie diffusion
1888	! turbulent du couche limit.
1889	!
1890	! Certains varibales de sorties de pbl_surface sont utiliser que pour
1891	! ecriture des fihiers hist_XXXX.nc, ces sont :
1892	! qsol, zq2m, s_pblh, s_lcl,
1893	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1894	! s_therm, s_trmb1, s_trmb2, s_trmb3,
1895	! zu10m, zv10m, fder,
1896	! zxqsurf, rh2m, zxfluxu, zxfluxv,
1897	! frugs, agesno, fsollw, fsolsw,
1898	! d_ts, fevap, fluxlat, t2m,
1899	! wfbils, wfbilo, fluxt, fluxu, fluxv,
1900	!
1901	! Certains ne sont pas utiliser du tout :
1902	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1903	!
1904
1905	! Calcul de l'humidite de saturation au niveau du sol
1906
1907
1908
1909	if (iflag_pbl/=0) then
1910
1911	!jyg+nrlmd<
1912	IF (prt_level .ge. 2 .and. mod(iflag_pbl_split,2) .eq. 1) THEN
1913	print *,'debut du splitting de la PBL'
1914	ENDIF
1915	!!!
1916	!=================================================================
1917	! PROVISOIRE : DECOUPLAGE PBL/WAKE
1918	! --------------------------------
1919	!
1920	!! wake_deltat_sav(:,:)=wake_deltat(:,:)
1921	!! wake_deltaq_sav(:,:)=wake_deltaq(:,:)
1922	!! wake_deltat(:,:)=0.
1923	!! wake_deltaq(:,:)=0.
1924	!=================================================================
1925	!>jyg+nrlmd
1926	!
1927	!-------gustiness calculation-------!
1928	IF (iflag_gusts==0) THEN
1929	gustiness(1:klon)=0
1930	ELSE IF (iflag_gusts==1) THEN
1931	do i = 1, klon
1932	gustiness(i)=f_gust_blale_bl(i)+f_gust_wkale_wake(i)
1933	enddo
1934	! ELSE IF (iflag_gusts==2) THEN
1935	! do i = 1, klon
1936	! gustiness(i)=f_gust_blale_bl(i)+sigma_wk(i)f_gust_wk*ale_wake(i) !! need to make sigma_wk accessible here
1937	! enddo
1938	! ELSE IF (iflag_gusts==3) THEN
1939	! do i = 1, klon
1940	! gustiness(i)=f_gust_blalp_bl(i)+f_gust_wkalp_wake(i)
1941	! enddo
1942	ENDIF
1943
1944
1945
1946	CALL pbl_surface( &
1947	dtime, date0, itap, days_elapsed+1, &
1948	debut, lafin, &
1949	rlon, rlat, rugoro, zrmu0, &
1950	zsig, sollwdown, pphi, cldt, &
1951	rain_fall, snow_fall, solsw, sollw, &
1952	gustiness, &
1953	t_seri, q_seri, u_seri, v_seri, &
1954	!nrlmd+jyg<
1955	wake_deltat, wake_deltaq, wake_cstar, wake_s, &
1956	!>nrlmd+jyg
1957	pplay, paprs, pctsrf, &
1958	ftsol,SFRWL,falb_dir,falb_dif,ustar,u10m,v10m,wstar, &
1959	!albedo SB <<<
1960	cdragh, cdragm, u1, v1, &
1961	!albedo SB >>>
1962	! albsol1, albsol2, sens, evap, &
1963	albsol_dir, albsol_dif, sens, evap, &
1964	!albedo SB <<<
1965	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
1966	zxtsol, zxfluxlat, zt2m, qsat2m, &
1967	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
1968	!nrlmd<
1969	!jyg<
1970	d_t_vdf_w, d_q_vdf_w, &
1971	d_t_vdf_x, d_q_vdf_x, &
1972	sens_x, zxfluxlat_x, sens_w, zxfluxlat_w, &
1973	!>jyg
1974	delta_tsurf,wake_dens, &
1975	cdragh_x,cdragh_w,cdragm_x,cdragm_w, &
1976	kh,kh_x,kh_w, &
1977	!>nrlmd
1978	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
1979	slab_wfbils, &
1980	qsol, zq2m, s_pblh, s_lcl, &
1981	!jyg<
1982	s_pblh_x, s_lcl_x, s_pblh_w, s_lcl_w, &
1983	!>jyg
1984	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
1985	s_therm, s_trmb1, s_trmb2, s_trmb3, &
1986	zustar, zu10m, zv10m, fder, &
1987	zxqsurf, rh2m, zxfluxu, zxfluxv, &
1988	z0m, z0h, agesno, fsollw, fsolsw, &
1989	d_ts, fevap, fluxlat, t2m, &
1990	wfbils, wfbilo, fluxt, fluxu, fluxv, &
1991	dsens, devap, zxsnow, &
1992	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke, &
1993	!nrlmd+jyg<
1994	wake_delta_pbl_TKE &
1995	!>nrlmd+jyg
1996	)
1997	!
1998	!=================================================================
1999	! PROVISOIRE : DECOUPLAGE PBL/WAKE
2000	! --------------------------------
2001	!
2002	!! wake_deltat(:,:)=wake_deltat_sav(:,:)
2003	!! wake_deltaq(:,:)=wake_deltaq_sav(:,:)
2004	!=================================================================
2005	!
2006	! Add turbulent diffusion tendency to the wake difference variables
2007	IF (mod(iflag_pbl_split,2) .NE. 0) THEN
2008	wake_deltat(:,:) = wake_deltat(:,:) + (d_t_vdf_w(:,:)-d_t_vdf_x(:,:))
2009	wake_deltaq(:,:) = wake_deltaq(:,:) + (d_q_vdf_w(:,:)-d_q_vdf_x(:,:))
2010	ENDIF
2011
2012
2013	!---------------------------------------------------------------------
2014	! ajout des tendances de la diffusion turbulente
2015	IF (klon_glo==1) THEN
2016	CALL add_pbl_tend &
2017	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf',abortphy)
2018	ELSE
2019	CALL add_phys_tend &
2020	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf',abortphy)
2021	ENDIF
2022	!--------------------------------------------------------------------
2023
2024	if (mydebug) then
2025	call writefield_phy('u_seri',u_seri,nbp_lev)
2026	call writefield_phy('v_seri',v_seri,nbp_lev)
2027	call writefield_phy('t_seri',t_seri,nbp_lev)
2028	call writefield_phy('q_seri',q_seri,nbp_lev)
2029	endif
2030
2031
2032	!albedo SB >>>
2033	albsol1=0.
2034	albsol2=0.
2035	falb1=0.
2036	falb2=0.
2037	select case(nsw)
2038	case(2)
2039	albsol1=albsol_dir(:,1)
2040	albsol2=albsol_dir(:,2)
2041	falb1=falb_dir(:,1,:)
2042	falb2=falb_dir(:,2,:)
2043	case(4)
2044	albsol1=albsol_dir(:,1)
2045	albsol2=albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)SFRWL(3)+albsol_dir(:,4)*SFRWL(4)
2046	albsol2=albsol2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2047	falb1=falb_dir(:,1,:)
2048	falb2=falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)SFRWL(3)+falb_dir(:,4,:)*SFRWL(4)
2049	falb2=falb2/(SFRWL(2)+SFRWL(3)+SFRWL(4))
2050	case(6)
2051	albsol1=albsol_dir(:,1)SFRWL(1)+albsol_dir(:,2)SFRWL(2)+albsol_dir(:,3)*SFRWL(3)
2052	albsol1=albsol1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2053	albsol2=albsol_dir(:,4)SFRWL(4)+albsol_dir(:,5)SFRWL(5)+albsol_dir(:,6)*SFRWL(6)
2054	albsol2=albsol2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2055	falb1=falb_dir(:,1,:)SFRWL(1)+falb_dir(:,2,:)SFRWL(2)+falb_dir(:,3,:)*SFRWL(3)
2056	falb1=falb1/(SFRWL(1)+SFRWL(2)+SFRWL(3))
2057	falb2=falb_dir(:,4,:)SFRWL(4)+falb_dir(:,5,:)SFRWL(5)+falb_dir(:,6,:)*SFRWL(6)
2058	falb2=falb2/(SFRWL(4)+SFRWL(5)+SFRWL(6))
2059	end select
2060	!albedo SB <<<
2061
2062
2063	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
2064	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
2065
2066
2067	IF (ip_ebil_phy.ge.2) THEN
2068	ztit='after surface_main'
2069	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
2070	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2071	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2072	call diagphy(cell_area,ztit,ip_ebil_phy &
2073	, zero_v, zero_v, zero_v, zero_v, sens &
2074	, evap , zero_v, zero_v, ztsol &
2075	, d_h_vcol, d_qt, d_ec &
2076	, fs_bound, fq_bound )
2077	END IF
2078
2079	ENDIF
2080	! =================================================================== c
2081	! Calcul de Qsat
2082
2083	DO k = 1, klev
2084	DO i = 1, klon
2085	zx_t = t_seri(i,k)
2086	IF (thermcep) THEN
2087	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2088	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2089	zx_qs = MIN(0.5,zx_qs)
2090	zcor = 1./(1.-retv*zx_qs)
2091	zx_qs = zx_qs*zcor
2092	ELSE
2093	IF (zx_t.LT.t_coup) THEN
2094	zx_qs = qsats(zx_t)/pplay(i,k)
2095	ELSE
2096	zx_qs = qsatl(zx_t)/pplay(i,k)
2097	ENDIF
2098	ENDIF
2099	zqsat(i,k)=zx_qs
2100	ENDDO
2101	ENDDO
2102
2103	if (prt_level.ge.1) then
2104	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
2105	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
2106	endif
2107	!
2108	! Appeler la convection (au choix)
2109	!
2110	DO k = 1, klev
2111	DO i = 1, klon
2112	conv_q(i,k) = d_q_dyn(i,k) &
2113	+ d_q_vdf(i,k)/dtime
2114	conv_t(i,k) = d_t_dyn(i,k) &
2115	+ d_t_vdf(i,k)/dtime
2116	ENDDO
2117	ENDDO
2118	IF (check) THEN
2119	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2120	WRITE(lunout,*) "avantcon=", za
2121	ENDIF
2122	zx_ajustq = .FALSE.
2123	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
2124	IF (zx_ajustq) THEN
2125	DO i = 1, klon
2126	z_avant(i) = 0.0
2127	ENDDO
2128	DO k = 1, klev
2129	DO i = 1, klon
2130	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
2131	*(paprs(i,k)-paprs(i,k+1))/RG
2132	ENDDO
2133	ENDDO
2134	ENDIF
2135
2136	! Calcule de vitesse verticale a partir de flux de masse verticale
2137	DO k = 1, klev
2138	DO i = 1, klon
2139	omega(i,k) = RG*flxmass_w(i,k) / cell_area(i)
2140	END DO
2141	END DO
2142	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
2143	omega(igout, :)
2144
2145	IF (iflag_con.EQ.1) THEN
2146	abort_message ='reactiver le call conlmd dans physiq.F'
2147	CALL abort_physic (modname,abort_message,1)
2148	! CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
2149	! . d_t_con, d_q_con,
2150	! . rain_con, snow_con, ibas_con, itop_con)
2151	ELSE IF (iflag_con.EQ.2) THEN
2152	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
2153	conv_t, conv_q, -evap, omega, &
2154	d_t_con, d_q_con, rain_con, snow_con, &
2155	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
2156	kcbot, kctop, kdtop, pmflxr, pmflxs)
2157	d_u_con = 0.
2158	d_v_con = 0.
2159
2160	WHERE (rain_con < 0.) rain_con = 0.
2161	WHERE (snow_con < 0.) snow_con = 0.
2162	DO i = 1, klon
2163	ibas_con(i) = klev+1 - kcbot(i)
2164	itop_con(i) = klev+1 - kctop(i)
2165	ENDDO
2166	ELSE IF (iflag_con.GE.3) THEN
2167	! nb of tracers for the KE convection:
2168	! MAF la partie traceurs est faite dans phytrac
2169	! on met ntra=1 pour limiter les appels mais on peut
2170	! supprimer les calculs / ftra.
2171	ntra = 1
2172
2173	!=========================================================================
2174	!ajout pour la parametrisation des poches froides: calcul de
2175	!t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
2176	do k=1,klev
2177	do i=1,klon
2178	if (iflag_wake>=1) then
2179	t_wake(i,k) = t_seri(i,k) &
2180	+(1-wake_s(i))*wake_deltat(i,k)
2181	q_wake(i,k) = q_seri(i,k) &
2182	+(1-wake_s(i))*wake_deltaq(i,k)
2183	t_undi(i,k) = t_seri(i,k) &
2184	-wake_s(i)*wake_deltat(i,k)
2185	q_undi(i,k) = q_seri(i,k) &
2186	-wake_s(i)*wake_deltaq(i,k)
2187	else
2188	t_wake(i,k) = t_seri(i,k)
2189	q_wake(i,k) = q_seri(i,k)
2190	t_undi(i,k) = t_seri(i,k)
2191	q_undi(i,k) = q_seri(i,k)
2192	endif
2193	enddo
2194	enddo
2195	!
2196	!jyg<
2197	! Perform dry adiabatic adjustment on wake profile
2198	! The corresponding tendencies are added to the convective tendencies
2199	! after the call to the convective scheme.
2200	IF (iflag_wake>=1) then
2201	IF (ok_adjwk) THEN
2202	limbas(:) = 1
2203	CALL ajsec(paprs, pplay, t_wake, q_wake, limbas, &
2204	d_t_adjwk, d_q_adjwk)
2205	ENDIF
2206	!
2207	DO k=1,klev
2208	DO i=1,klon
2209	IF (wake_s(i) .GT. 1.e-3) THEN
2210	t_wake(i,k) = t_wake(i,k) + d_t_adjwk(i,k)
2211	q_wake(i,k) = q_wake(i,k) + d_q_adjwk(i,k)
2212	wake_deltat(i,k) = wake_deltat(i,k) + d_t_adjwk(i,k)
2213	wake_deltaq(i,k) = wake_deltaq(i,k) + d_q_adjwk(i,k)
2214	ENDIF
2215	ENDDO
2216	ENDDO
2217	ENDIF ! (iflag_wake>=1)
2218	!>jyg
2219	!
2220
2221	! Calcul de l'energie disponible ALE (J/kg) et de la puissance
2222	! disponible ALP (W/m2) pour le soulevement des particules dans
2223	! le modele convectif
2224	!
2225	do i = 1,klon
2226	ALE(i) = 0.
2227	ALP(i) = 0.
2228	enddo
2229	!
2230	!calcul de ale_wake et alp_wake
2231	if (iflag_wake>=1) then
2232	if (itap .le. it_wape_prescr) then
2233	do i = 1,klon
2234	ale_wake(i) = wape_prescr
2235	alp_wake(i) = fip_prescr
2236	enddo
2237	else
2238	do i = 1,klon
2239	!jyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
2240	!cc ale_wake(i) = 0.5wake_cstar(i)*2
2241	ale_wake(i) = wake_pe(i)
2242	alp_wake(i) = wake_fip(i)
2243	enddo
2244	endif
2245	else
2246	do i = 1,klon
2247	ale_wake(i) = 0.
2248	alp_wake(i) = 0.
2249	enddo
2250	endif
2251	!combinaison avec ale et alp de couche limite: constantes si pas
2252	!de couplage, valeurs calculees dans le thermique sinon
2253	if (iflag_coupl.eq.0) then
2254	if (debut.and.prt_level.gt.9) &
2255	WRITE(lunout,*)'ALE et ALP imposes'
2256	do i = 1,klon
2257	!on ne couple que ale
2258	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2259	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2260	!on ne couple que alp
2261	! ALP(i) = alp_wake(i) + Alp_bl(i)
2262	ALP(i) = alp_wake(i) + alp_bl_prescr
2263	enddo
2264	else
2265	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2266	! do i = 1,klon
2267	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2268	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
2269	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2270	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2271	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2272	! enddo
2273
2274	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2275	! Modif FH 2010/04/27. Sans doute temporaire.
2276	! Deux options pour le alp_offset : constant si >?? 0 ou
2277	! proportionnel ??a w si <0
2278	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2279	! Estimation d'une vitesse verticale effective pour ALP
2280	if (1==0) THEN
2281	www(1:klon)=0.
2282	do k=2,klev-1
2283	do i=1,klon
2284	www(i)=max(www(i),-omega(i,k)RDt_seri(i,k)/(RG*paprs(i,k)) &
2285	& zw2(i,k)zw2(i,k))
2286	! if (paprs(i,k)>pbase(i)) then
2287	! calcul approche de la vitesse verticale en m/s
2288	! www(i)=max(www(i),-omega(i,k)RDtemp(i,k)/(RG*paprs(i,k))
2289	! endif
2290	! Le 0.1 est en gros H / ps = 1e5 / 1e4
2291	enddo
2292	enddo
2293	do i=1,klon
2294	if (www(i)>0. .and. ale_bl(i)>0. ) www(i)=www(i)/ale_bl(i)
2295	enddo
2296	ENDIF
2297
2298
2299	do i = 1,klon
2300	ALE(i) = max(ale_wake(i),Ale_bl(i))
2301	!cc nrlmd le 10/04/2012----------Stochastic triggering--------------
2302	if (iflag_trig_bl.ge.1) then
2303	ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
2304	endif
2305	!cc fin nrlmd le 10/04/2012
2306	if (alp_offset>=0.) then
2307	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2308	else
2309	abort_message ='Ne pas passer la car www non calcule'
2310	CALL abort_physic (modname,abort_message,1)
2311
2312	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2313	! _ _
2314	! Ajout d'une composante 3 * A * w w'2 a w'3 avec w=www : w max sous pbase
2315	! ou A est la fraction couverte par les ascendances w'
2316	! on utilise le fait que A * w'3 = ALP
2317	! et donc A * w'2 ~ ALP / sqrt(ALE) (on ajoute 0.1 pour les
2318	! singularites)
2319	ALP(i)=alp_wake(i)(1.+3.www(i)/( sqrt(ale_wake(i))+0.1) ) &
2320	& +alp_bl(i) (1.+3.www(i)/( sqrt(ale_bl(i)) +0.1) )
2321	! ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2322	! if (alp(i)<0.) then
2323	! print*,'ALP ',alp(i),alp_wake(i) &
2324	! ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2325	! endif
2326	endif
2327	enddo
2328	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2329
2330	endif
2331	do i=1,klon
2332	if (alp(i)>alp_max) then
2333	IF(prt_level>9)WRITE(lunout,*) &
2334	'WARNING SUPER ALP (seuil=',alp_max, &
2335	'): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2336	alp(i)=alp_max
2337	endif
2338	if (ale(i)>ale_max) then
2339	IF(prt_level>9)WRITE(lunout,*) &
2340	'WARNING SUPER ALE (seuil=',ale_max, &
2341	'): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2342	ale(i)=ale_max
2343	endif
2344	enddo
2345
2346	!fin calcul ale et alp
2347	!=======================================================================
2348
2349
2350	! sb, oct02:
2351	! Schema de convection modularise et vectorise:
2352	! (driver commun aux versions 3 et 4)
2353	!
2354	IF (ok_cvl) THEN ! new driver for convectL
2355	!
2356	!jyg<
2357	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2358	! Calculate the upmost level of deep convection loops: k_upper_cv
2359	! (near 22 km)
2360	izero = klon/2+1/klon
2361	k_upper_cv = klev
2362	DO k = klev,1,-1
2363	IF (pphi(izero,k) > 22.e4) k_upper_cv = k
2364	ENDDO
2365	IF (prt_level .ge. 5) THEN
2366	Print *, 'upmost level of deep convection loops: k_upper_cv = ',k_upper_cv
2367	ENDIF
2368	!
2369	!>jyg
2370	IF (type_trac == 'repr') THEN
2371	nbtr_tmp=ntra
2372	ELSE
2373	nbtr_tmp=nbtr
2374	END IF
2375	!jyg iflag_con est dans clesphys
2376	!c CALL concvl (iflag_con,iflag_clos,
2377	CALL concvl (iflag_clos, &
2378	dtime, paprs, pplay, k_upper_cv, t_undi,q_undi, &
2379	t_wake,q_wake,wake_s, &
2380	u_seri,v_seri,tr_seri,nbtr_tmp, &
2381	ALE,ALP, &
2382	sig1,w01, &
2383	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2384	rain_con, snow_con, ibas_con, itop_con, sigd, &
2385	ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
2386	Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2387	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2388	! RomP >>>
2389	!! . pmflxr,pmflxs,da,phi,mp,
2390	!! . ftd,fqd,lalim_conv,wght_th)
2391	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2392	ftd,fqd,lalim_conv,wght_th, &
2393	ev, ep,epmlmMm,eplaMm, &
2394	wdtrainA,wdtrainM,wght_cvfd,qtc_cv,sigt_cv, &
2395	tau_cld_cv,coefw_cld_cv)
2396	! RomP <<<
2397
2398	!IM begin
2399	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2400	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2401	!IM end
2402	!IM cf. FH
2403	clwcon0=qcondc
2404	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2405
2406	do i = 1, klon
2407	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2408	enddo
2409	!
2410	!jyg<
2411	! Add the tendency due to the dry adjustment of the wake profile
2412	IF (iflag_wake>=1) THEN
2413	DO k=1,klev
2414	DO i=1,klon
2415	ftd(i,k) = ftd(i,k) + wake_s(i)*d_t_adjwk(i,k)/dtime
2416	fqd(i,k) = fqd(i,k) + wake_s(i)*d_q_adjwk(i,k)/dtime
2417	d_t_con(i,k) = d_t_con(i,k) + wake_s(i)*d_t_adjwk(i,k)
2418	d_q_con(i,k) = d_q_con(i,k) + wake_s(i)*d_q_adjwk(i,k)
2419	ENDDO
2420	ENDDO
2421	ENDIF
2422	!>jyg
2423	!
2424	ELSE ! ok_cvl
2425
2426	! MAF conema3 ne contient pas les traceurs
2427	CALL conema3 (dtime, &
2428	paprs,pplay,t_seri,q_seri, &
2429	u_seri,v_seri,tr_seri,ntra, &
2430	sig1,w01, &
2431	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2432	rain_con, snow_con, ibas_con, itop_con, &
2433	upwd,dnwd,dnwd0,bas,top, &
2434	Ma,cape,tvp,rflag, &
2435	pbase &
2436	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2437	,clwcon0)
2438
2439	ENDIF ! ok_cvl
2440
2441	!
2442	! Correction precip
2443	rain_con = rain_con * cvl_corr
2444	snow_con = snow_con * cvl_corr
2445	!
2446
2447	IF (.NOT. ok_gust) THEN
2448	do i = 1, klon
2449	wd(i)=0.0
2450	enddo
2451	ENDIF
2452
2453	! =================================================================== c
2454	! Calcul des proprietes des nuages convectifs
2455	!
2456
2457	! calcul des proprietes des nuages convectifs
2458	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2459	IF (iflag_cld_cv == 0) THEN
2460	call clouds_gno &
2461	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2462	ELSE
2463	call clouds_bigauss &
2464	(klon,klev,q_seri,zqsat,qtc_cv,sigt_cv,ptconv,ratqsc,rnebcon0)
2465	ENDIF
2466
2467
2468	! =================================================================== c
2469
2470	DO i = 1, klon
2471	itop_con(i) = min(max(itop_con(i),1),klev)
2472	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2473	ENDDO
2474
2475	DO i = 1, klon
2476	ema_pcb(i) = paprs(i,ibas_con(i))
2477	ENDDO
2478	DO i = 1, klon
2479	! L'idicage de itop_con peut cacher un pb potentiel
2480	! FH sous la dictee de JYG, CR
2481	ema_pct(i) = paprs(i,itop_con(i)+1)
2482
2483	if (itop_con(i).gt.klev-3) then
2484	if(prt_level >= 9) then
2485	write(lunout,*)'La convection monte trop haut '
2486	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2487	endif
2488	endif
2489	ENDDO
2490	ELSE IF (iflag_con.eq.0) THEN
2491	write(lunout,*) 'On n appelle pas la convection'
2492	clwcon0=0.
2493	rnebcon0=0.
2494	d_t_con=0.
2495	d_q_con=0.
2496	d_u_con=0.
2497	d_v_con=0.
2498	rain_con=0.
2499	snow_con=0.
2500	bas=1
2501	top=1
2502	ELSE
2503	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2504	call abort_physic("physiq", "", 1)
2505	ENDIF
2506
2507	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2508	! . d_u_con, d_v_con)
2509
2510	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
2511	'convection',abortphy)
2512
2513	!----------------------------------------------------------------------------
2514
2515	if (mydebug) then
2516	call writefield_phy('u_seri',u_seri,nbp_lev)
2517	call writefield_phy('v_seri',v_seri,nbp_lev)
2518	call writefield_phy('t_seri',t_seri,nbp_lev)
2519	call writefield_phy('q_seri',q_seri,nbp_lev)
2520	endif
2521
2522	!IM
2523	IF (ip_ebil_phy.ge.2) THEN
2524	ztit='after convect'
2525	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
2526	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2527	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2528	call diagphy(cell_area,ztit,ip_ebil_phy &
2529	, zero_v, zero_v, zero_v, zero_v, zero_v &
2530	, zero_v, rain_con, snow_con, ztsol &
2531	, d_h_vcol, d_qt, d_ec &
2532	, fs_bound, fq_bound )
2533	END IF
2534	!
2535	IF (check) THEN
2536	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
2537	WRITE(lunout,*)"aprescon=", za
2538	zx_t = 0.0
2539	za = 0.0
2540	DO i = 1, klon
2541	za = za + cell_area(i)/REAL(klon)
2542	zx_t = zx_t + (rain_con(i)+ &
2543	snow_con(i))*cell_area(i)/REAL(klon)
2544	ENDDO
2545	zx_t = zx_t/za*dtime
2546	WRITE(lunout,*)"Precip=", zx_t
2547	ENDIF
2548	IF (zx_ajustq) THEN
2549	DO i = 1, klon
2550	z_apres(i) = 0.0
2551	ENDDO
2552	DO k = 1, klev
2553	DO i = 1, klon
2554	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2555	*(paprs(i,k)-paprs(i,k+1))/RG
2556	ENDDO
2557	ENDDO
2558	DO i = 1, klon
2559	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2560	/z_apres(i)
2561	ENDDO
2562	DO k = 1, klev
2563	DO i = 1, klon
2564	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2565	z_factor(i).LT.(1.0-1.0E-08)) THEN
2566	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2567	ENDIF
2568	ENDDO
2569	ENDDO
2570	ENDIF
2571	zx_ajustq=.FALSE.
2572
2573	!
2574	!=============================================================================
2575	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
2576	!pour la couche limite diffuse pour l instant
2577	!
2578	!
2579	!!! nrlmd le 22/03/2011---Si on met les poches hors des thermiques il faut rajouter cette
2580	!------------------------- tendance calculée hors des poches froides
2581	!
2582	if (iflag_wake>=1) then
2583	DO k=1,klev
2584	DO i=1,klon
2585	dt_dwn(i,k) = ftd(i,k)
2586	dq_dwn(i,k) = fqd(i,k)
2587	M_dwn(i,k) = dnwd0(i,k)
2588	M_up(i,k) = upwd(i,k)
2589	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2590	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2591	ENDDO
2592	ENDDO
2593	!nrlmd+jyg<
2594	DO k=1,klev
2595	DO i=1,klon
2596	wdt_PBL(i,k) = 0.
2597	wdq_PBL(i,k) = 0.
2598	udt_PBL(i,k) = 0.
2599	udq_PBL(i,k) = 0.
2600	ENDDO
2601	ENDDO
2602	!
2603	IF (mod(iflag_pbl_split,2) .EQ. 1) THEN
2604	DO k=1,klev
2605	DO i=1,klon
2606	wdt_PBL(i,k) = wdt_PBL(i,k) + d_t_vdf_w(i,k)/dtime
2607	wdq_PBL(i,k) = wdq_PBL(i,k) + d_q_vdf_w(i,k)/dtime
2608	udt_PBL(i,k) = udt_PBL(i,k) + d_t_vdf_x(i,k)/dtime
2609	udq_PBL(i,k) = udq_PBL(i,k) + d_q_vdf_x(i,k)/dtime
2610	!! dt_dwn(i,k) = dt_dwn(i,k) + d_t_vdf_w(i,k)/dtime
2611	!! dq_dwn(i,k) = dq_dwn(i,k) + d_q_vdf_w(i,k)/dtime
2612	!! dt_a (i,k) = dt_a(i,k) + d_t_vdf_x(i,k)/dtime
2613	!! dq_a (i,k) = dq_a(i,k) + d_q_vdf_x(i,k)/dtime
2614	ENDDO
2615	ENDDO
2616	ENDIF
2617	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2618	DO k=1,klev
2619	DO i=1,klon
2620	!! dt_dwn(i,k) = dt_dwn(i,k) + 0.
2621	!! dq_dwn(i,k) = dq_dwn(i,k) + 0.
2622	!! dt_a(i,k) = dt_a(i,k) + d_t_ajs(i,k)/dtime
2623	!! dq_a(i,k) = dq_a(i,k) + d_q_ajs(i,k)/dtime
2624	udt_PBL(i,k) = udt_PBL(i,k) + d_t_ajs(i,k)/dtime
2625	udq_PBL(i,k) = udq_PBL(i,k) + d_q_ajs(i,k)/dtime
2626	ENDDO
2627	ENDDO
2628	ENDIF
2629	!>nrlmd+jyg
2630
2631	IF (iflag_wake==2) THEN
2632	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2633	DO k = 1,klev
2634	dt_dwn(:,k)= dt_dwn(:,k)+ &
2635	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2636	dq_dwn(:,k)= dq_dwn(:,k)+ &
2637	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2638	ENDDO
2639	ELSEIF (iflag_wake==3) THEN
2640	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2641	DO k = 1,klev
2642	DO i=1,klon
2643	IF (rneb(i,k)==0.) THEN
2644	! On ne tient compte des tendances qu'en dehors des nuages (c'est ï¿½\| dire
2645	! a priri dans une region ou l'eau se reevapore).
2646	dt_dwn(i,k)= dt_dwn(i,k)+ &
2647	ok_wk_lsp(i)*d_t_lsc(i,k)/dtime
2648	dq_dwn(i,k)= dq_dwn(i,k)+ &
2649	ok_wk_lsp(i)*d_q_lsc(i,k)/dtime
2650	ENDIF
2651	ENDDO
2652	ENDDO
2653	ENDIF
2654
2655	!
2656	!calcul caracteristiques de la poche froide
2657	call calWAKE (paprs,pplay,dtime &
2658	,t_seri,q_seri,omega &
2659	,dt_dwn,dq_dwn,M_dwn,M_up &
2660	,dt_a,dq_a,sigd &
2661	,wdt_PBL,wdq_PBL &
2662	,udt_PBL,udq_PBL &
2663	,wake_deltat,wake_deltaq,wake_dth &
2664	,wake_h,wake_s,wake_dens &
2665	,wake_pe,wake_fip,wake_gfl &
2666	,dt_wake,dq_wake &
2667	,wake_k, t_undi,q_undi &
2668	,wake_omgbdth,wake_dp_omgb &
2669	,wake_dtKE,wake_dqKE &
2670	,wake_dtPBL,wake_dqPBL &
2671	,wake_omg,wake_dp_deltomg &
2672	,wake_spread,wake_Cstar,wake_d_deltat_gw &
2673	,wake_ddeltat,wake_ddeltaq)
2674	!
2675	!-------------------------------------------------------------------------
2676	! ajout des tendances des poches froides
2677	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2678	! coherent avec les autres d_t_...
2679	d_t_wake(:,:)=dt_wake(:,:)*dtime
2680	d_q_wake(:,:)=dq_wake(:,:)*dtime
2681	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake',abortphy)
2682	!------------------------------------------------------------------------
2683
2684	endif ! (iflag_wake>=1)
2685	!
2686	!===================================================================
2687	!JYG
2688	IF (ip_ebil_phy.ge.2) THEN
2689	ztit='after wake'
2690	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
2691	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2692	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2693	call diagphy(cell_area,ztit,ip_ebil_phy &
2694	, zero_v, zero_v, zero_v, zero_v, zero_v &
2695	, zero_v, zero_v, zero_v, ztsol &
2696	, d_h_vcol, d_qt, d_ec &
2697	, fs_bound, fq_bound )
2698	END IF
2699
2700	! print*,'apres callwake iflag_cld_th=', iflag_cld_th
2701	!
2702	!===================================================================
2703	! Convection seche (thermiques ou ajustement)
2704	!===================================================================
2705	!
2706	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
2707	,seuil_inversion,weak_inversion,dthmin)
2708
2709
2710
2711	d_t_ajsb(:,:)=0.
2712	d_q_ajsb(:,:)=0.
2713	d_t_ajs(:,:)=0.
2714	d_u_ajs(:,:)=0.
2715	d_v_ajs(:,:)=0.
2716	d_q_ajs(:,:)=0.
2717	clwcon0th(:,:)=0.
2718	!
2719	! fm_therm(:,:)=0.
2720	! entr_therm(:,:)=0.
2721	! detr_therm(:,:)=0.
2722	!
2723	IF(prt_level>9)WRITE(lunout,*) &
2724	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
2725	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2726	if(iflag_thermals<0) then
2727	! Rien
2728	! ====
2729	IF(prt_level>9)WRITE(lunout,*)'pas de convection seche'
2730
2731
2732	else
2733
2734	! Thermiques
2735	! ==========
2736	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
2737	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2738
2739
2740	!cc nrlmd le 10/04/2012
2741	DO k=1,klev+1
2742	DO i=1,klon
2743	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2744	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2745	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2746	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2747	ENDDO
2748	ENDDO
2749	!cc fin nrlmd le 10/04/2012
2750
2751	if (iflag_thermals>=1) then
2752	!jyg<
2753	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2754	! Appel des thermiques avec les profils exterieurs aux poches
2755	DO k=1,klev
2756	DO i=1,klon
2757	t_therm(i,k) = t_seri(i,k) - wake_s(i)*wake_deltat(i,k)
2758	q_therm(i,k) = q_seri(i,k) - wake_s(i)*wake_deltaq(i,k)
2759	ENDDO
2760	ENDDO
2761	ELSE
2762	! Appel des thermiques avec les profils moyens
2763	DO k=1,klev
2764	DO i=1,klon
2765	t_therm(i,k) = t_seri(i,k)
2766	q_therm(i,k) = q_seri(i,k)
2767	ENDDO
2768	ENDDO
2769	ENDIF
2770	!>jyg
2771	call calltherm(pdtphys &
2772	,pplay,paprs,pphi,weak_inversion &
2773	!! ,u_seri,v_seri,t_seri,q_seri,zqsat,debut & !jyg
2774	,u_seri,v_seri,t_therm,q_therm,zqsat,debut & !jyg
2775	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
2776	,fm_therm,entr_therm,detr_therm &
2777	,zqasc,clwcon0th,lmax_th,ratqscth &
2778	,ratqsdiff,zqsatth &
2779	!on rajoute ale et alp, et les caracteristiques de la couche alim
2780	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
2781	,ztv,zpspsk,ztla,zthl &
2782	!cc nrlmd le 10/04/2012
2783	,pbl_tke_input,pctsrf,omega,cell_area &
2784	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
2785	,n2,s2,ale_bl_stat &
2786	,therm_tke_max,env_tke_max &
2787	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
2788	,alp_bl_conv,alp_bl_stat &
2789	!cc fin nrlmd le 10/04/2012
2790	,zqla,ztva )
2791	!
2792	!jyg<
2793	IF (mod(iflag_pbl_split/2,2) .EQ. 1) THEN
2794	! Si les thermiques ne sont presents que hors des poches, la tendance moyenne
2795	! associée doit etre multipliee par la fraction surfacique qu'ils couvrent.
2796	DO k=1,klev
2797	DO i=1,klon
2798	!
2799	wake_deltat(i,k) = wake_deltat(i,k) - d_t_ajs(i,k)
2800	wake_deltaq(i,k) = wake_deltaq(i,k) - d_q_ajs(i,k)
2801	t_seri(i,k) = t_therm(i,k) + wake_s(i)*wake_deltat(i,k)
2802	q_seri(i,k) = q_therm(i,k) + wake_s(i)*wake_deltaq(i,k)
2803	!
2804	d_u_ajs(i,k) = d_u_ajs(i,k)*(1.-wake_s(i))
2805	d_v_ajs(i,k) = d_v_ajs(i,k)*(1.-wake_s(i))
2806	d_t_ajs(i,k) = d_t_ajs(i,k)*(1.-wake_s(i))
2807	d_q_ajs(i,k) = d_q_ajs(i,k)*(1.-wake_s(i))
2808	!
2809	ENDDO
2810	ENDDO
2811	ELSE
2812	DO k=1,klev
2813	DO i=1,klon
2814	t_seri(i,k) = t_therm(i,k)
2815	q_seri(i,k) = q_therm(i,k)
2816	ENDDO
2817	ENDDO
2818	ENDIF
2819	!>jyg
2820
2821	!cc nrlmd le 10/04/2012
2822	!-----------Stochastic triggering-----------
2823	if (iflag_trig_bl.ge.1) then
2824	!
2825	IF (prt_level .GE. 10) THEN
2826	print *,'cin, ale_bl_stat, alp_bl_stat ', &
2827	cin, ale_bl_stat, alp_bl_stat
2828	ENDIF
2829
2830
2831	!----Initialisations
2832	do i=1,klon
2833	proba_notrig(i)=1.
2834	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
2835	if ( random_notrig(i) > random_notrig_max ) random_notrig(i)=0.
2836	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
2837	tau_trig(i)=tau_trig_shallow
2838	else
2839	tau_trig(i)=tau_trig_deep
2840	endif
2841	enddo
2842	!
2843	IF (prt_level .GE. 10) THEN
2844	print *,'random_notrig, tau_trig ', &
2845	random_notrig, tau_trig
2846	print *,'s_trig,s2,n2 ', &
2847	s_trig,s2,n2
2848	ENDIF
2849
2850	!Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
2851	IF (iflag_trig_bl.eq.1) then
2852
2853	!----Tirage al\'eatoire et calcul de ale_bl_trig
2854	do i=1,klon
2855	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
2856	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2857	(n2(i)*dtime/tau_trig(i))
2858	! print *, 'proba_notrig(i) ',proba_notrig(i)
2859	if (random_notrig(i) .ge. proba_notrig(i)) then
2860	ale_bl_trig(i)=ale_bl_stat(i)
2861	else
2862	ale_bl_trig(i)=0.
2863	endif
2864	else
2865	proba_notrig(i)=1.
2866	random_notrig(i)=0.
2867	ale_bl_trig(i)=0.
2868	endif
2869	enddo
2870
2871	ELSE IF (iflag_trig_bl.ge.2) then
2872
2873	do i=1,klon
2874	if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) ) then
2875	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2876	(n2(i)*dtime/tau_trig(i))
2877	! print *, 'proba_notrig(i) ',proba_notrig(i)
2878	if (random_notrig(i) .ge. proba_notrig(i)) then
2879	ale_bl_trig(i)=Ale_bl(i)
2880	else
2881	ale_bl_trig(i)=0.
2882	endif
2883	else
2884	proba_notrig(i)=1.
2885	random_notrig(i)=0.
2886	ale_bl_trig(i)=0.
2887	endif
2888	enddo
2889
2890	ENDIF
2891
2892	!
2893	IF (prt_level .GE. 10) THEN
2894	print *,'proba_notrig, ale_bl_trig ', &
2895	proba_notrig, ale_bl_trig
2896	ENDIF
2897
2898	endif !(iflag_trig_bl)
2899
2900	!-----------Statistical closure-----------
2901	if (iflag_clos_bl.eq.1) then
2902
2903	do i=1,klon
2904	!CR: alp probabiliste
2905	if (ale_bl_trig(i).gt.0.) then
2906	alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
2907	endif
2908	enddo
2909
2910	else if (iflag_clos_bl.eq.2) then
2911
2912	!CR: alp calculee dans thermcell_main
2913	do i=1,klon
2914	alp_bl(i)=alp_bl_stat(i)
2915	enddo
2916
2917	else
2918
2919	alp_bl_stat(:)=0.
2920
2921	endif !(iflag_clos_bl)
2922
2923	IF (prt_level .GE. 10) THEN
2924	print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
2925	ENDIF
2926
2927	!cc fin nrlmd le 10/04/2012
2928
2929	! ----------------------------------------------------------------------
2930	! Transport de la TKE par les panaches thermiques.
2931	! FH : 2010/02/01
2932	! if (iflag_pbl.eq.10) then
2933	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2934	! s rg,paprs,pbl_tke)
2935	! endif
2936	! ----------------------------------------------------------------------
2937	!IM/FH: 2011/02/23
2938	! Couplage Thermiques/Emanuel seulement si T<0
2939	if (iflag_coupl==2) then
2940	IF (prt_level .GE. 10) THEN
2941	print*,'Couplage Thermiques/Emanuel seulement si T<0'
2942	ENDIF
2943	do i=1,klon
2944	if (t_seri(i,lmax_th(i))>273.) then
2945	Ale_bl(i)=0.
2946	endif
2947	enddo
2948	endif
2949
2950	do i=1,klon
2951	! zmax_th(i)=pphi(i,lmax_th(i))/rg
2952	!CR:04/05/12:correction calcul zmax
2953	zmax_th(i)=zmax0(i)
2954	enddo
2955
2956	endif
2957
2958
2959	! Ajustement sec
2960	! ==============
2961
2962	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
2963	! a partir du sommet des thermiques.
2964	! Dans le cas contraire, on demarre au niveau 1.
2965
2966	if (iflag_thermals>=13.or.iflag_thermals<=0) then
2967
2968	if(iflag_thermals.eq.0) then
2969	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2970	limbas(:)=1
2971	else
2972	limbas(:)=lmax_th(:)
2973	endif
2974
2975	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2976	! pour des test de convergence numerique.
2977	! Le nouveau ajsec est a priori mieux, meme pour le cas
2978	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2979	! non nulles numeriquement pour des mailles non concernees.
2980
2981	if (iflag_thermals==0) then
2982	! Calling adjustment alone (but not the thermal plume model)
2983	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
2984	, d_t_ajsb, d_q_ajsb)
2985	else if (iflag_thermals>0) then
2986	! Calling adjustment above the top of thermal plumes
2987	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
2988	, d_t_ajsb, d_q_ajsb)
2989	endif
2990
2991	!-----------------------------------------------------------------------
2992	! ajout des tendances de l'ajustement sec ou des thermiques
2993	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs,'ajsb',abortphy)
2994	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2995	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2996
2997	!---------------------------------------------------------------------
2998
2999	endif
3000
3001	endif
3002	!
3003	!===================================================================
3004	!IM
3005	IF (ip_ebil_phy.ge.2) THEN
3006	ztit='after dry_adjust'
3007	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
3008	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3009	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3010	call diagphy(cell_area,ztit,ip_ebil_phy &
3011	, zero_v, zero_v, zero_v, zero_v, zero_v &
3012	, zero_v, zero_v, zero_v, ztsol &
3013	, d_h_vcol, d_qt, d_ec &
3014	, fs_bound, fq_bound )
3015	END IF
3016
3017
3018	!-------------------------------------------------------------------------
3019	! Computation of ratqs, the width (normalized) of the subrid scale
3020	! water distribution
3021	CALL calcratqs(klon,klev,prt_level,lunout, &
3022	iflag_ratqs,iflag_con,iflag_cld_th,pdtphys, &
3023	ratqsbas,ratqshaut,tau_ratqs,fact_cldcon, &
3024	ptconv,ptconvth,clwcon0th, rnebcon0th, &
3025	paprs,pplay,q_seri,zqsat,fm_therm, &
3026	ratqs,ratqsc)
3027
3028
3029	!
3030	! Appeler le processus de condensation a grande echelle
3031	! et le processus de precipitation
3032	!-------------------------------------------------------------------------
3033	IF (prt_level .GE.10) THEN
3034	print *,'itap, ->fisrtilp ',itap
3035	ENDIF
3036	!
3037	CALL fisrtilp(dtime,paprs,pplay, &
3038	t_seri, q_seri,ptconv,ratqs, &
3039	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
3040	rain_lsc, snow_lsc, &
3041	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
3042	frac_impa, frac_nucl, beta_prec_fisrt, &
3043	prfl, psfl, rhcl, &
3044	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, &
3045	iflag_ice_thermo)
3046	!
3047	WHERE (rain_lsc < 0) rain_lsc = 0.
3048	WHERE (snow_lsc < 0) snow_lsc = 0.
3049
3050	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs,'lsc',abortphy)
3051	!---------------------------------------------------------------------------
3052	DO k = 1, klev
3053	DO i = 1, klon
3054	cldfra(i,k) = rneb(i,k)
3055	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
3056	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
3057	ENDDO
3058	ENDDO
3059	IF (check) THEN
3060	za = qcheck(klon,klev,paprs,q_seri,ql_seri,cell_area)
3061	WRITE(lunout,*)"apresilp=", za
3062	zx_t = 0.0
3063	za = 0.0
3064	DO i = 1, klon
3065	za = za + cell_area(i)/REAL(klon)
3066	zx_t = zx_t + (rain_lsc(i) &
3067	+ snow_lsc(i))*cell_area(i)/REAL(klon)
3068	ENDDO
3069	zx_t = zx_t/za*dtime
3070	WRITE(lunout,*)"Precip=", zx_t
3071	ENDIF
3072	!IM
3073	IF (ip_ebil_phy.ge.2) THEN
3074	ztit='after fisrt'
3075	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
3076	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3077	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3078	call diagphy(cell_area,ztit,ip_ebil_phy &
3079	, zero_v, zero_v, zero_v, zero_v, zero_v &
3080	, zero_v, rain_lsc, snow_lsc, ztsol &
3081	, d_h_vcol, d_qt, d_ec &
3082	, fs_bound, fq_bound )
3083	END IF
3084
3085	if (mydebug) then
3086	call writefield_phy('u_seri',u_seri,nbp_lev)
3087	call writefield_phy('v_seri',v_seri,nbp_lev)
3088	call writefield_phy('t_seri',t_seri,nbp_lev)
3089	call writefield_phy('q_seri',q_seri,nbp_lev)
3090	endif
3091
3092	!
3093	!-------------------------------------------------------------------
3094	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
3095	!-------------------------------------------------------------------
3096
3097	! 1. NUAGES CONVECTIFS
3098	!
3099	!IM cf FH
3100	! IF (iflag_cld_th.eq.-1) THEN ! seulement pour Tiedtke
3101	IF (iflag_cld_th.le.-1) THEN ! seulement pour Tiedtke
3102	snow_tiedtke=0.
3103	! print*,'avant calcul de la pseudo precip '
3104	! print*,'iflag_cld_th',iflag_cld_th
3105	if (iflag_cld_th.eq.-1) then
3106	rain_tiedtke=rain_con
3107	else
3108	! print*,'calcul de la pseudo precip '
3109	rain_tiedtke=0.
3110	! print*,'calcul de la pseudo precip 0'
3111	do k=1,klev
3112	do i=1,klon
3113	if (d_q_con(i,k).lt.0.) then
3114	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
3115	*(paprs(i,k)-paprs(i,k+1))/rg
3116	endif
3117	enddo
3118	enddo
3119	endif
3120	!
3121	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
3122	!
3123
3124	! Nuages diagnostiques pour Tiedtke
3125	CALL diagcld1(paprs,pplay, &
3126	!IM cf FH . rain_con,snow_con,ibas_con,itop_con,
3127	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
3128	diafra,dialiq)
3129	DO k = 1, klev
3130	DO i = 1, klon
3131	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3132	cldliq(i,k) = dialiq(i,k)
3133	cldfra(i,k) = diafra(i,k)
3134	ENDIF
3135	ENDDO
3136	ENDDO
3137
3138	ELSE IF (iflag_cld_th.ge.3) THEN
3139	! On prend pour les nuages convectifs le max du calcul de la
3140	! convection et du calcul du pas de temps precedent diminue d'un facteur
3141	! facttemps
3142	facteur = pdtphys *facttemps
3143	do k=1,klev
3144	do i=1,klon
3145	rnebcon(i,k)=rnebcon(i,k)*facteur
3146	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k)) &
3147	then
3148	rnebcon(i,k)=rnebcon0(i,k)
3149	clwcon(i,k)=clwcon0(i,k)
3150	endif
3151	enddo
3152	enddo
3153
3154	!
3155	!jq - introduce the aerosol direct and first indirect radiative forcings
3156	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
3157	IF (flag_aerosol .gt. 0) THEN
3158	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3159	IF (.NOT. aerosol_couple) THEN
3160	!
3161	CALL readaerosol_optic( &
3162	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3163	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3164	mass_solu_aero, mass_solu_aero_pi, &
3165	tau_aero, piz_aero, cg_aero, &
3166	tausum_aero, tau3d_aero)
3167	ENDIF
3168	ELSE ! RRTM radiation
3169	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
3170	abort_message='config_inca=aero et rrtm=1 impossible'
3171	call abort_physic(modname,abort_message,1)
3172	ELSE
3173	!
3174	#ifdef CPP_RRTM
3175	IF (NSW.EQ.6) THEN
3176	!--new aerosol properties
3177	!
3178	CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
3179	new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3180	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3181	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
3182	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
3183	tausum_aero, tau3d_aero)
3184
3185	ELSE IF (NSW.EQ.2) THEN
3186	!--for now we use the old aerosol properties
3187	!
3188	CALL readaerosol_optic( &
3189	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
3190	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
3191	mass_solu_aero, mass_solu_aero_pi, &
3192	tau_aero, piz_aero, cg_aero, &
3193	tausum_aero, tau3d_aero)
3194	!
3195	!--natural aerosols
3196	tau_aero_sw_rrtm(:,:,1,:)=tau_aero(:,:,3,:)
3197	piz_aero_sw_rrtm(:,:,1,:)=piz_aero(:,:,3,:)
3198	cg_aero_sw_rrtm (:,:,1,:)=cg_aero (:,:,3,:)
3199	!--all aerosols
3200	tau_aero_sw_rrtm(:,:,2,:)=tau_aero(:,:,2,:)
3201	piz_aero_sw_rrtm(:,:,2,:)=piz_aero(:,:,2,:)
3202	cg_aero_sw_rrtm (:,:,2,:)=cg_aero (:,:,2,:)
3203	ELSE
3204	abort_message='Only NSW=2 or 6 are possible with aerosols and iflag_rrtm=1'
3205	call abort_physic(modname,abort_message,1)
3206	ENDIF
3207
3208	CALL aeropt_lw_rrtm
3209	!
3210	#else
3211	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
3212	call abort_physic(modname,abort_message,1)
3213	#endif
3214	!
3215	ENDIF
3216	ENDIF
3217	ELSE
3218	tausum_aero(:,:,:) = 0.
3219	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
3220	tau_aero(:,:,:,:) = 1.e-15
3221	piz_aero(:,:,:,:) = 1.
3222	cg_aero(:,:,:,:) = 0.
3223	ELSE
3224	tau_aero_sw_rrtm(:,:,:,:) = 1.e-15
3225	tau_aero_lw_rrtm(:,:,:,:) = 1.e-15
3226	piz_aero_sw_rrtm(:,:,:,:) = 1.0
3227	cg_aero_sw_rrtm(:,:,:,:) = 0.0
3228	ENDIF
3229	ENDIF
3230	!
3231	!--STRAT AEROSOL
3232	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
3233	IF (flag_aerosol_strat) THEN
3234	IF (prt_level .GE.10) THEN
3235	PRINT *,'appel a readaerosolstrat', mth_cur
3236	ENDIF
3237	IF (iflag_rrtm.EQ.0) THEN
3238	CALL readaerosolstrato(debut)
3239	ELSE
3240	#ifdef CPP_RRTM
3241	CALL readaerosolstrato_rrtm(debut)
3242	#else
3243
3244	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
3245	call abort_physic(modname,abort_message,1)
3246	#endif
3247	ENDIF
3248	ENDIF
3249	!--fin STRAT AEROSOL
3250
3251
3252	! On prend la somme des fractions nuageuses et des contenus en eau
3253
3254	if (iflag_cld_th>=5) then
3255
3256	do k=1,klev
3257	ptconvth(:,k)=fm_therm(:,k+1)>0.
3258	enddo
3259
3260	if (iflag_coupl==4) then
3261
3262	! Dans le cas iflag_coupl==4, on prend la somme des convertures
3263	! convectives et lsc dans la partie des thermiques
3264	! Le controle par iflag_coupl est peut etre provisoire.
3265	do k=1,klev
3266	do i=1,klon
3267	if (ptconv(i,k).and.ptconvth(i,k)) then
3268	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3269	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3270	else if (ptconv(i,k)) then
3271	cldfra(i,k)=rnebcon(i,k)
3272	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3273	endif
3274	enddo
3275	enddo
3276
3277	else if (iflag_coupl==5) then
3278	do k=1,klev
3279	do i=1,klon
3280	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
3281	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
3282	enddo
3283	enddo
3284
3285	else
3286
3287	! Si on est sur un point touche par la convection profonde et pas
3288	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
3289	! de la convection profonde.
3290
3291	!IM/FH: 2011/02/23
3292	! definition des points sur lesquels ls thermiques sont actifs
3293
3294	do k=1,klev
3295	do i=1,klon
3296	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
3297	cldfra(i,k)=rnebcon(i,k)
3298	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
3299	endif
3300	enddo
3301	enddo
3302
3303	endif
3304
3305	else
3306
3307	! Ancienne version
3308	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
3309	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
3310	endif
3311
3312	ENDIF
3313
3314	! plulsc(:)=0.
3315	! do k=1,klev,-1
3316	! do i=1,klon
3317	! zzz=prfl(:,k)+psfl(:,k)
3318	! if (.not.ptconvth.zzz.gt.0.)
3319	! enddo prfl, psfl,
3320	! enddo
3321	!
3322	! 2. NUAGES STARTIFORMES
3323	!
3324	IF (ok_stratus) THEN
3325	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
3326	DO k = 1, klev
3327	DO i = 1, klon
3328	IF (diafra(i,k).GT.cldfra(i,k)) THEN
3329	cldliq(i,k) = dialiq(i,k)
3330	cldfra(i,k) = diafra(i,k)
3331	ENDIF
3332	ENDDO
3333	ENDDO
3334	ENDIF
3335	!
3336	! Precipitation totale
3337	!
3338	DO i = 1, klon
3339	rain_fall(i) = rain_con(i) + rain_lsc(i)
3340	snow_fall(i) = snow_con(i) + snow_lsc(i)
3341	ENDDO
3342	!IM
3343	IF (ip_ebil_phy.ge.2) THEN
3344	ztit="after diagcld"
3345	CALL diagetpq(cell_area,ztit,ip_ebil_phy,2,2,dtime &
3346	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3347	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3348	call diagphy(cell_area,ztit,ip_ebil_phy &
3349	, zero_v, zero_v, zero_v, zero_v, zero_v &
3350	, zero_v, zero_v, zero_v, ztsol &
3351	, d_h_vcol, d_qt, d_ec &
3352	, fs_bound, fq_bound )
3353	END IF
3354	!
3355	! Calculer l'humidite relative pour diagnostique
3356	!
3357	DO k = 1, klev
3358	DO i = 1, klon
3359	zx_t = t_seri(i,k)
3360	IF (thermcep) THEN
3361	if (iflag_ice_thermo.eq.0) then
3362	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
3363	else
3364	zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))
3365	endif
3366	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
3367	zx_qs = MIN(0.5,zx_qs)
3368	zcor = 1./(1.-retv*zx_qs)
3369	zx_qs = zx_qs*zcor
3370	ELSE
3371	IF (zx_t.LT.t_coup) THEN
3372	zx_qs = qsats(zx_t)/pplay(i,k)
3373	ELSE
3374	zx_qs = qsatl(zx_t)/pplay(i,k)
3375	ENDIF
3376	ENDIF
3377	zx_rh(i,k) = q_seri(i,k)/zx_qs
3378	zqsat(i,k)=zx_qs
3379	ENDDO
3380	ENDDO
3381
3382	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
3383	! equivalente a 2m (tpote) pour diagnostique
3384	!
3385	DO i = 1, klon
3386	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
3387	IF (thermcep) THEN
3388	IF(zt2m(i).LT.RTT) then
3389	Lheat=RLSTT
3390	ELSE
3391	Lheat=RLVTT
3392	ENDIF
3393	ELSE
3394	IF (zt2m(i).LT.RTT) THEN
3395	Lheat=RLSTT
3396	ELSE
3397	Lheat=RLVTT
3398	ENDIF
3399	ENDIF
3400	tpote(i) = tpot(i)* &
3401	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
3402	ENDDO
3403
3404	IF (type_trac == 'inca') THEN
3405	#ifdef INCA
3406	CALL VTe(VTphysiq)
3407	CALL VTb(VTinca)
3408	calday = REAL(days_elapsed + 1) + jH_cur
3409
3410	call chemtime(itap+itau_phy-1, date0, dtime, itap)
3411	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
3412	CALL AEROSOL_METEO_CALC( &
3413	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
3414	prfl,psfl,pctsrf,cell_area,rlat,rlon,u10m,v10m)
3415	END IF
3416
3417	zxsnow_dummy(:) = 0.0
3418
3419	CALL chemhook_begin (calday, &
3420	days_elapsed+1, &
3421	jH_cur, &
3422	pctsrf(1,1), &
3423	rlat, &
3424	rlon, &
3425	cell_area, &
3426	paprs, &
3427	pplay, &
3428	coefh(1:klon,1:klev,is_ave), &
3429	pphi, &
3430	t_seri, &
3431	u, &
3432	v, &
3433	wo(:, :, 1), &
3434	q_seri, &
3435	zxtsol, &
3436	zxsnow_dummy, &
3437	solsw, &
3438	albsol1, &
3439	rain_fall, &
3440	snow_fall, &
3441	itop_con, &
3442	ibas_con, &
3443	cldfra, &
3444	nbp_lon, &
3445	nbp_lat-1, &
3446	tr_seri, &
3447	ftsol, &
3448	paprs, &
3449	cdragh, &
3450	cdragm, &
3451	pctsrf, &
3452	pdtphys, &
3453	itap)
3454
3455	CALL VTe(VTinca)
3456	CALL VTb(VTphysiq)
3457	#endif
3458	END IF !type_trac = inca
3459	!
3460	! Calculer les parametres optiques des nuages et quelques
3461	! parametres pour diagnostiques:
3462	!
3463
3464	IF (aerosol_couple) THEN
3465	mass_solu_aero(:,:) = ccm(:,:,1)
3466	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3467	END IF
3468
3469	if (ok_newmicro) then
3470	IF (iflag_rrtm.NE.0) THEN
3471	#ifdef CPP_RRTM
3472	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3473	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 pour ok_cdnc'
3474	call abort_physic(modname,abort_message,1)
3475	endif
3476	#else
3477
3478	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
3479	call abort_physic(modname,abort_message,1)
3480	#endif
3481	ENDIF
3482	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3483	paprs, pplay, t_seri, cldliq, cldfra, &
3484	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3485	flwp, fiwp, flwc, fiwc, &
3486	mass_solu_aero, mass_solu_aero_pi, &
3487	cldtaupi, re, fl, ref_liq, ref_ice, &
3488	ref_liq_pi, ref_ice_pi)
3489	else
3490	CALL nuage (paprs, pplay, &
3491	t_seri, cldliq, cldfra, cldtau, cldemi, &
3492	cldh, cldl, cldm, cldt, cldq, &
3493	ok_aie, &
3494	mass_solu_aero, mass_solu_aero_pi, &
3495	bl95_b0, bl95_b1, &
3496	cldtaupi, re, fl)
3497	endif
3498	!
3499	!IM betaCRF
3500	!
3501	cldtaurad = cldtau
3502	cldtaupirad = cldtaupi
3503	cldemirad = cldemi
3504
3505	!
3506	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3507	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3508	!
3509	! global
3510	!
3511	DO k=1, klev
3512	DO i=1, klon
3513	if (pplay(i,k).GE.pfree) THEN
3514	beta(i,k) = beta_pbl
3515	else
3516	beta(i,k) = beta_free
3517	endif
3518	if (mskocean_beta) THEN
3519	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3520	endif
3521	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3522	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3523	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3524	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3525	ENDDO
3526	ENDDO
3527	!
3528	else
3529	!
3530	! regional
3531	!
3532	DO k=1, klev
3533	DO i=1,klon
3534	!
3535	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
3536	rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
3537	if (pplay(i,k).GE.pfree) THEN
3538	beta(i,k) = beta_pbl
3539	else
3540	beta(i,k) = beta_free
3541	endif
3542	if (mskocean_beta) THEN
3543	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3544	endif
3545	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3546	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3547	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3548	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3549	endif
3550	!
3551	ENDDO
3552	ENDDO
3553	!
3554	endif
3555	!
3556	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3557	!
3558	IF (MOD(itaprad,radpas).EQ.0) THEN
3559
3560	!albedo SB >>>
3561	if(ok_chlorophyll)then
3562	print*,"-- reading chlorophyll"
3563	call readchlorophyll(debut)
3564	endif
3565	!do i=1,klon
3566	!if(chl_con(i)>1.) print*,i,chl_con(i),pctsrf(i,is_ter)
3567	!enddo
3568	!albedo SB <<<
3569
3570
3571	if (mydebug) then
3572	call writefield_phy('u_seri',u_seri,nbp_lev)
3573	call writefield_phy('v_seri',v_seri,nbp_lev)
3574	call writefield_phy('t_seri',t_seri,nbp_lev)
3575	call writefield_phy('q_seri',q_seri,nbp_lev)
3576	endif
3577
3578	!
3579	!sonia : If Iflag_radia >=2, pertubation of some variables input to radiation
3580	!(DICE)
3581	!
3582	IF (iflag_radia .ge. 2) THEN
3583	zsav_tsol (:) = zxtsol(:)
3584	call perturb_radlwsw(zxtsol,iflag_radia)
3585	ENDIF
3586
3587	IF (aerosol_couple) THEN
3588	#ifdef INCA
3589	CALL radlwsw_inca &
3590	(kdlon,kflev,dist, rmu0, fract, solaire, &
3591	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3592	wo(:, :, 1), &
3593	cldfrarad, cldemirad, cldtaurad, &
3594	heat,heat0,cool,cool0,albpla, &
3595	topsw,toplw,solsw,sollw, &
3596	sollwdown, &
3597	topsw0,toplw0,solsw0,sollw0, &
3598	lwdn0, lwdn, lwup0, lwup, &
3599	swdn0, swdn, swup0, swup, &
3600	ok_ade, ok_aie, &
3601	tau_aero, piz_aero, cg_aero, &
3602	topswad_aero, solswad_aero, &
3603	topswad0_aero, solswad0_aero, &
3604	topsw_aero, topsw0_aero, &
3605	solsw_aero, solsw0_aero, &
3606	cldtaupirad, &
3607	topswai_aero, solswai_aero)
3608
3609	#endif
3610	ELSE
3611	!
3612	!IM calcul radiatif pour le cas actuel
3613	!
3614	RCO2 = RCO2_act
3615	RCH4 = RCH4_act
3616	RN2O = RN2O_act
3617	RCFC11 = RCFC11_act
3618	RCFC12 = RCFC12_act
3619	!
3620	IF (prt_level .GE.10) THEN
3621	print *,' ->radlwsw, number 1 '
3622	ENDIF
3623	!
3624	CALL radlwsw &
3625	(dist, rmu0, fract, &
3626	!albedo SB >>>
3627	! paprs, pplay,zxtsol,albsol1, albsol2, &
3628	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3629	!albedo SB <<<
3630	t_seri,q_seri,wo, &
3631	cldfrarad, cldemirad, cldtaurad, &
3632	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3633	flag_aerosol_strat, &
3634	tau_aero, piz_aero, cg_aero, &
3635	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! Rajoute par OB pour RRTM
3636	tau_aero_lw_rrtm, &
3637	cldtaupirad,new_aod, &
3638	zqsat, flwc, fiwc, &
3639	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3640	heat,heat0,cool,cool0,albpla, &
3641	topsw,toplw,solsw,sollw, &
3642	sollwdown, &
3643	topsw0,toplw0,solsw0,sollw0, &
3644	lwdn0, lwdn, lwup0, lwup, &
3645	swdn0, swdn, swup0, swup, &
3646	topswad_aero, solswad_aero, &
3647	topswai_aero, solswai_aero, &
3648	topswad0_aero, solswad0_aero, &
3649	topsw_aero, topsw0_aero, &
3650	solsw_aero, solsw0_aero, &
3651	topswcf_aero, solswcf_aero, &
3652	!-C. Kleinschmitt for LW diagnostics
3653	toplwad_aero, sollwad_aero,&
3654	toplwai_aero, sollwai_aero, &
3655	toplwad0_aero, sollwad0_aero,&
3656	!-end
3657	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3658	ZSWFT0_i, ZFSDN0, ZFSUP0)
3659
3660	!
3661	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3662	!IM des taux doit etre different du taux actuel
3663	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3664	!
3665	if (ok_4xCO2atm) then
3666	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3667	RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3668	RCFC12_per.NE.RCFC12_act) THEN
3669	!
3670	RCO2 = RCO2_per
3671	RCH4 = RCH4_per
3672	RN2O = RN2O_per
3673	RCFC11 = RCFC11_per
3674	RCFC12 = RCFC12_per
3675	!
3676	IF (prt_level .GE.10) THEN
3677	print *,' ->radlwsw, number 2 '
3678	ENDIF
3679	!
3680	CALL radlwsw &
3681	(dist, rmu0, fract, &
3682	!albedo SB >>>
3683	! paprs, pplay,zxtsol,albsol1, albsol2, &
3684	paprs, pplay,zxtsol,SFRWL,albsol_dir, albsol_dif, &
3685	!albedo SB <<<
3686	t_seri,q_seri,wo, &
3687	cldfra, cldemi, cldtau, &
3688	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3689	flag_aerosol_strat, &
3690	tau_aero, piz_aero, cg_aero, &
3691	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! Rajoute par OB pour RRTM
3692	tau_aero_lw_rrtm, &
3693	cldtaupi,new_aod, &
3694	zqsat, flwc, fiwc, &
3695	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3696	heatp,heat0p,coolp,cool0p,albplap, &
3697	topswp,toplwp,solswp,sollwp, &
3698	sollwdownp, &
3699	topsw0p,toplw0p,solsw0p,sollw0p, &
3700	lwdn0p, lwdnp, lwup0p, lwupp, &
3701	swdn0p, swdnp, swup0p, swupp, &
3702	topswad_aerop, solswad_aerop, &
3703	topswai_aerop, solswai_aerop, &
3704	topswad0_aerop, solswad0_aerop, &
3705	topsw_aerop, topsw0_aerop, &
3706	solsw_aerop, solsw0_aerop, &
3707	topswcf_aerop, solswcf_aerop, &
3708	!-C. Kleinschmitt for LW diagnostics
3709	toplwad_aerop, sollwad_aerop,&
3710	toplwai_aerop, sollwai_aerop, &
3711	toplwad0_aerop, sollwad0_aerop,&
3712	!-end
3713	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3714	ZSWFT0_i, ZFSDN0, ZFSUP0)
3715	endif
3716	endif
3717	!
3718	ENDIF ! aerosol_couple
3719	itaprad = 0
3720	!
3721	! If Iflag_radia >=2, reset pertubed variables
3722	!
3723	IF (iflag_radia .ge. 2) THEN
3724	zxtsol(:) = zsav_tsol (:)
3725	ENDIF
3726	ENDIF ! MOD(itaprad,radpas)
3727	itaprad = itaprad + 1
3728
3729	IF (iflag_radia.eq.0) THEN
3730	IF (prt_level.ge.9) THEN
3731	PRINT *,'--------------------------------------------------'
3732	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3733	PRINT *,'>>>> heat et cool mis a zero '
3734	PRINT *,'--------------------------------------------------'
3735	END IF
3736	heat=0.
3737	cool=0.
3738	sollw=0. ! MPL 01032011
3739	solsw=0.
3740	radsol=0.
3741	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3742	swup0=0.
3743	lwup=0.
3744	lwup0=0.
3745	lwdn=0.
3746	lwdn0=0.
3747	END IF
3748
3749	!
3750	! Calculer radsol a l'exterieur de radlwsw
3751	! pour prendre en compte le cycle diurne
3752	! recode par Olivier Boucher en sept 2015
3753	!
3754	radsol=solsw*swradcorr+sollw
3755	if (ok_4xCO2atm) then
3756	radsolp=solswp*swradcorr+sollwp
3757	endif
3758
3759	!
3760	! Ajouter la tendance des rayonnements (tous les pas)
3761	! avec une correction pour le cycle diurne dans le SW
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