Context Navigation

physiq.F90 @ 1948

Last change on this file since 1948 was 1948, checked in by idelkadi, 10 years ago

Option pour tourner sans thermiques ni ajustement sec si ifalg_thermals<0 (phylmd/physiq.F90)

Modifications permettant d'imposer des profils initiaux de traceurs

lus dans un tracer.inp.001 (phylmd/1DUTILS.h, 1D_read_forc_cases.h et lmdz1d.F)

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

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

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

Download in other formats:

Original Format