Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

physiq.F90

Last change on this file was 2381, checked in by acozic, 9 years ago

Make some commit to fit with INCA coupling

merge with rev 2180 on trunk
merge with rev 2185 on trunk
merge with rev 2200 on trunk
change gregorian calendar in wxios

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: 128.9 KB

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

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

Download in other formats:

Original Format