Context Navigation

physiq.F90 @ 2105

Last change on this file since 2105 was 2105, checked in by fhourdin, 10 years ago

Modification des options pour inclure un aspect grande Ã©chelle
dans le calcul de Alp et pour tenir compte des prÃ©cpitatios
sratiformes pour alimenter les poches froides.
ContrÃ´lÃ©s respectivement par les flags alp_ofet=X, avec X<0.
(la valeur par dÃ©faut est 0) et iflag_wake=3.

Pour le calcul de Alp :
! Ajout d'une composante 3 * A * w w'2 a w'3 avec w=www : w max sous pbase
! ou A est la fraction couverte par les ascendances w'
! on utilise le fait que A * w'3 = ALP
! et donc A * w'2 ~ ALP / sqrt(ALE) (on ajoute 0.1 pour les
! singularites)

Modified options to account for large scale velocity in ALP
and for feeding of wakes by large scale rainfall.
Controled respectively by alp_offset=X, with X<0. (default being 0) and iflag_wake=3.

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: 127.3 KB

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

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

Download in other formats:

Original Format