Context Navigation

physiq.F90 @ 1956

Last change on this file since 1956 was 1956, checked in by lguez, 11 years ago
Correcting bad spacing left over from fixed source form.
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: 120.7 KB

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

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

Download in other formats:

Original Format