Context Navigation

physiq.F90 @ 1907

Last change on this file since 1907 was 1907, checked in by lguez, 10 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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.1 KB

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

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

Download in other formats:

Original Format