Context Navigation

physiq.F90 @ 2309

Last change on this file since 2309 was 2309, checked in by jyg, 9 years ago
Modification of physiq.F90: when wakes are present, an adiabatic adjustment is applied to the wake region before calling the Emanuel convection scheme. The switch 'ok_adjwk', read from .def files, controls this adjustment.
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: 140.6 KB

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

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

Download in other formats:

Original Format