Context Navigation

physiq.F90 @ 1921

Last change on this file since 1921 was 1921, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1909:1920 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 119.7 KB

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

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

Download in other formats:

Original Format