Context Navigation

physiq_mod.F90 @ 2418

Last change on this file since 2418 was 2418, checked in by Ehouarn Millour, 8 years ago

Improving the physics/dynamics interface:

added module callphysiq_mod.F90 in dynphy_lonlat/phy* which contains the routine "call_physiq" which is called by calfis* and calls the physics. This way different "physiq" routine from different physics packages may be called: The calfis* routines now exposes all available fields that might be transmitted to physiq but which is actually send (ie: expected/needed by physiq) is decided in call_physiq.
turned "physiq.F90" into module "physiq_mod.F90" for better control of "physiq" arguments. Extracted embeded "gr_fi_ecrit" as self-standing routine (but note that this routine actually only works in serial mode).

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 145.0 KB

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

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

Download in other formats:

Original Format