Context Navigation

physiq.F90 @ 1917

Last change on this file since 1917 was 1917, checked in by acozic, 11 years ago
Change coefh dimensions in several call of subroutine (to fit with dimension declared in these routines)
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.4 KB

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

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

Download in other formats:

Original Format