Context Navigation

physiq.F90 @ 2359

Last change on this file since 2359 was 2359, checked in by oboucher, 10 years ago
Recoded diurnal cycle in physiq.F90 iflag_diurnal_cycle is the new flag =0 for no diurnal cycle =1 for old diurnal cycle =2 for new diurnal cycle with updates in between two calls to radiation according to how rmu0 has changed over the previous timestep of physiq
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 144.6 KB

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

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

Download in other formats:

Original Format