Context Navigation

physiq.F90 @ 2942

Last change on this file since 2942 was 1877, checked in by Laurent Fairhead, 11 years ago

Ménage sur le code pour éliminer les calculs spécifiques ISCCP hors COSP
Un certain nombre de variables non utilisées dans physiq.F90 ont aussi
été supprimée. Environ 840 lignes supprimées du code physiq.F90

Code cleanup to eliminate specific references to ISCCP outside the COSP
library. Unused variables have been cleaned up from the physiq.F90 routine
as well.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 119.1 KB

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

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

Download in other formats:

Original Format