Context Navigation

physiq.F90 @ 2300

Last change on this file since 2300 was 2283, checked in by Laurent Fairhead, 9 years ago

Backport of trunk revisions 2271, 2279, 2280, 2282 into LMDZ6_rc0 branch:

modifications for NMC/XIOS

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

Line
1	! $Id: physiq.F90 2283 2015-05-18 15:25:29Z fhourdin $
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	itau_phy)
1310
1311	CALL VTe(VTinca)
1312	CALL VTb(VTphysiq)
1313	#endif
1314	END IF
1315	!
1316	!
1317	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1318	! Nouvelle initialisation pour le rayonnement RRTM
1319	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1320
1321	call iniradia(klon,klev,paprs(1,1:klev+1))
1322
1323	!$omp single
1324	if (read_climoz >= 1) then
1325	call open_climoz(ncid_climoz, press_climoz)
1326	END IF
1327	!$omp end single
1328	!
1329	!IM betaCRF
1330	pfree=70000. !Pa
1331	beta_pbl=1.
1332	beta_free=1.
1333	lon1_beta=-180.
1334	lon2_beta=+180.
1335	lat1_beta=90.
1336	lat2_beta=-90.
1337	mskocean_beta=.FALSE.
1338
1339	OPEN(99,file='beta_crf.data',status='old', &
1340	form='formatted',err=9999)
1341	READ(99,*,end=9998) pfree
1342	READ(99,*,end=9998) beta_pbl
1343	READ(99,*,end=9998) beta_free
1344	READ(99,*,end=9998) lon1_beta
1345	READ(99,*,end=9998) lon2_beta
1346	READ(99,*,end=9998) lat1_beta
1347	READ(99,*,end=9998) lat2_beta
1348	READ(99,*,end=9998) mskocean_beta
1349	9998 Continue
1350	CLOSE(99)
1351	9999 Continue
1352	WRITE(,)'pfree=',pfree
1353	WRITE(,)'beta_pbl=',beta_pbl
1354	WRITE(,)'beta_free=',beta_free
1355	WRITE(,)'lon1_beta=',lon1_beta
1356	WRITE(,)'lon2_beta=',lon2_beta
1357	WRITE(,)'lat1_beta=',lat1_beta
1358	WRITE(,)'lat2_beta=',lat2_beta
1359	WRITE(,)'mskocean_beta=',mskocean_beta
1360	ENDIF
1361	!
1362	! ************** Fin de IF ( debut ) *************
1363	!
1364	!
1365	! Incrementer le compteur de la physique
1366	!
1367	itap = itap + 1
1368	!
1369	!
1370	! Update fraction of the sub-surfaces (pctsrf) and
1371	! initialize, where a new fraction has appeared, all variables depending
1372	! on the surface fraction.
1373	!
1374	CALL change_srf_frac(itap, dtime, days_elapsed+1, &
1375	pctsrf, falb1, falb2, ftsol, ustar, u10m, v10m, pbl_tke)
1376
1377
1378	! Update time and other variables in Reprobus
1379	IF (type_trac == 'repr') THEN
1380	#ifdef REPROBUS
1381	CALL Init_chem_rep_xjour(jD_cur-jD_ref+day_ref)
1382	print*,'xjour equivalent rjourvrai',jD_cur-jD_ref+day_ref
1383	CALL Rtime(debut)
1384	#endif
1385	END IF
1386
1387
1388	! Tendances bidons pour les processus qui n'affectent pas certaines
1389	! variables.
1390	du0(:,:)=0.
1391	dv0(:,:)=0.
1392	dt0 = 0.
1393	dq0(:,:)=0.
1394	dql0(:,:)=0.
1395	dqi0(:,:)=0.
1396	!
1397	! Mettre a zero des variables de sortie (pour securite)
1398	!
1399	DO i = 1, klon
1400	d_ps(i) = 0.0
1401	ENDDO
1402	DO k = 1, klev
1403	DO i = 1, klon
1404	d_t(i,k) = 0.0
1405	d_u(i,k) = 0.0
1406	d_v(i,k) = 0.0
1407	ENDDO
1408	ENDDO
1409	DO iq = 1, nqtot
1410	DO k = 1, klev
1411	DO i = 1, klon
1412	d_qx(i,k,iq) = 0.0
1413	ENDDO
1414	ENDDO
1415	ENDDO
1416	da(:,:)=0.
1417	mp(:,:)=0.
1418	phi(:,:,:)=0.
1419	! RomP >>>
1420	phi2(:,:,:)=0.
1421	beta_prec_fisrt(:,:)=0.
1422	beta_prec(:,:)=0.
1423	epmlmMm(:,:,:)=0.
1424	eplaMm(:,:)=0.
1425	d1a(:,:)=0.
1426	dam(:,:)=0.
1427	pmflxr=0.
1428	pmflxs=0.
1429	! RomP <<<
1430
1431	!
1432	! Ne pas affecter les valeurs entrees de u, v, h, et q
1433	!
1434	DO k = 1, klev
1435	DO i = 1, klon
1436	t_seri(i,k) = t(i,k)
1437	u_seri(i,k) = u(i,k)
1438	v_seri(i,k) = v(i,k)
1439	q_seri(i,k) = qx(i,k,ivap)
1440	ql_seri(i,k) = qx(i,k,iliq)
1441	!CR: ATTENTION, on rajoute la variable glace
1442	if (nqo.eq.2) then
1443	qs_seri(i,k) = 0.
1444	else if (nqo.eq.3) then
1445	qs_seri(i,k) = qx(i,k,isol)
1446	endif
1447	ENDDO
1448	ENDDO
1449	tke0(:,:)=pbl_tke(:,:,is_ave)
1450	!CR:Nombre de traceurs de l'eau: nqo
1451	! IF (nqtot.GE.3) THEN
1452	IF (nqtot.GE.(nqo+1)) THEN
1453	! DO iq = 3, nqtot
1454	DO iq = nqo+1, nqtot
1455	DO k = 1, klev
1456	DO i = 1, klon
1457	! tr_seri(i,k,iq-2) = qx(i,k,iq)
1458	tr_seri(i,k,iq-nqo) = qx(i,k,iq)
1459	ENDDO
1460	ENDDO
1461	ENDDO
1462	ELSE
1463	DO k = 1, klev
1464	DO i = 1, klon
1465	tr_seri(i,k,1) = 0.0
1466	ENDDO
1467	ENDDO
1468	ENDIF
1469	!
1470	DO i = 1, klon
1471	ztsol(i) = 0.
1472	ENDDO
1473	DO nsrf = 1, nbsrf
1474	DO i = 1, klon
1475	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1476	ENDDO
1477	ENDDO
1478	!IM
1479	IF (ip_ebil_phy.ge.1) THEN
1480	ztit='after dynamic'
1481	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
1482	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1483	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1484	! Comme les tendances de la physique sont ajoute dans la dynamique,
1485	! on devrait avoir que la variation d'entalpie par la dynamique
1486	! est egale a la variation de la physique au pas de temps precedent.
1487	! Donc la somme de ces 2 variations devrait etre nulle.
1488	call diagphy(airephy,ztit,ip_ebil_phy &
1489	, zero_v, zero_v, zero_v, zero_v, zero_v &
1490	, zero_v, zero_v, zero_v, ztsol &
1491	, d_h_vcol+d_h_vcol_phy, d_qt, 0. &
1492	, fs_bound, fq_bound )
1493	END IF
1494
1495	! Diagnostiquer la tendance dynamique
1496	!
1497	IF (ancien_ok) THEN
1498	DO k = 1, klev
1499	DO i = 1, klon
1500	d_u_dyn(i,k) = (u_seri(i,k)-u_ancien(i,k))/dtime
1501	d_v_dyn(i,k) = (v_seri(i,k)-v_ancien(i,k))/dtime
1502	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1503	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1504	ENDDO
1505	ENDDO
1506	!!! RomP >>> td dyn traceur
1507	IF (nqtot.GE.3) THEN
1508	DO iq = 3, nqtot
1509	DO k = 1, klev
1510	DO i = 1, klon
1511	d_tr_dyn(i,k,iq-2)= &
1512	(tr_seri(i,k,iq-2)-tr_ancien(i,k,iq-2))/dtime
1513	! iiq=niadv(iq)
1514	! print*,i,k," d_tr_dyn",d_tr_dyn(i,k,iq-2),"tra:",iq,tname(iiq)
1515	ENDDO
1516	ENDDO
1517	ENDDO
1518	ENDIF
1519	!!! RomP <<<
1520	ELSE
1521	DO k = 1, klev
1522	DO i = 1, klon
1523	d_u_dyn(i,k) = 0.0
1524	d_v_dyn(i,k) = 0.0
1525	d_t_dyn(i,k) = 0.0
1526	d_q_dyn(i,k) = 0.0
1527	ENDDO
1528	ENDDO
1529	!!! RomP >>> td dyn traceur
1530	IF (nqtot.GE.3) THEN
1531	DO iq = 3, nqtot
1532	DO k = 1, klev
1533	DO i = 1, klon
1534	d_tr_dyn(i,k,iq-2)= 0.0
1535	ENDDO
1536	ENDDO
1537	ENDDO
1538	ENDIF
1539	!!! RomP <<<
1540	ancien_ok = .TRUE.
1541	ENDIF
1542	!
1543	! Ajouter le geopotentiel du sol:
1544	!
1545	DO k = 1, klev
1546	DO i = 1, klon
1547	zphi(i,k) = pphi(i,k) + pphis(i)
1548	ENDDO
1549	ENDDO
1550	!
1551	! Verifier les temperatures
1552	!
1553	!IM BEG
1554	IF (check) THEN
1555	amn=MIN(ftsol(1,is_ter),1000.)
1556	amx=MAX(ftsol(1,is_ter),-1000.)
1557	DO i=2, klon
1558	amn=MIN(ftsol(i,is_ter),amn)
1559	amx=MAX(ftsol(i,is_ter),amx)
1560	ENDDO
1561	!
1562	PRINT*,' debut avant hgardfou min max ftsol',itap,amn,amx
1563	ENDIF !(check) THEN
1564	!IM END
1565	!
1566	CALL hgardfou(t_seri,ftsol,'debutphy')
1567	!
1568	!IM BEG
1569	IF (check) THEN
1570	amn=MIN(ftsol(1,is_ter),1000.)
1571	amx=MAX(ftsol(1,is_ter),-1000.)
1572	DO i=2, klon
1573	amn=MIN(ftsol(i,is_ter),amn)
1574	amx=MAX(ftsol(i,is_ter),amx)
1575	ENDDO
1576	!
1577	PRINT*,' debut apres hgardfou min max ftsol',itap,amn,amx
1578	ENDIF !(check) THEN
1579	!IM END
1580	!
1581	! Mettre en action les conditions aux limites (albedo, sst, etc.).
1582	! Prescrire l'ozone et calculer l'albedo sur l'ocean.
1583	!
1584	if (read_climoz >= 1) then
1585	! Ozone from a file
1586	! Update required ozone index:
1587	ro3i = int((days_elapsed + jh_cur - jh_1jan) / year_len * 360.) + 1
1588	if (ro3i == 361) ro3i = 360
1589	! (This should never occur, except perhaps because of roundup
1590	! error. See documentation.)
1591	if (ro3i /= co3i) then
1592	! Update ozone field:
1593	if (read_climoz == 1) then
1594	call regr_pr_av(ncid_climoz, (/"tro3"/), julien=ro3i, &
1595	press_in_edg=press_climoz, paprs=paprs, v3=wo)
1596	else
1597	! read_climoz == 2
1598	call regr_pr_av(ncid_climoz, (/"tro3 ", "tro3_daylight"/), &
1599	julien=ro3i, press_in_edg=press_climoz, paprs=paprs, v3=wo)
1600	end if
1601	! Convert from mole fraction of ozone to column density of ozone in a
1602	! cell, in kDU:
1603	forall (l = 1: read_climoz) wo(:, :, l) = wo(:, :, l) * rmo3 / rmd &
1604	* zmasse / dobson_u / 1e3
1605	! (By regridding ozone values for LMDZ only once every 360th of
1606	! year, we have already neglected the variation of pressure in one
1607	! 360th of year. So do not recompute "wo" at each time step even if
1608	! "zmasse" changes a little.)
1609	co3i = ro3i
1610	end if
1611	ELSEIF (MOD(itap-1,lmt_pas) == 0) THEN
1612	! Once per day, update ozone from Royer:
1613
1614	IF (solarlong0<-999.) then
1615	! Generic case with evolvoing season
1616	zzz=real(days_elapsed+1)
1617	ELSE IF (abs(solarlong0-1000.)<1.e-4) then
1618	! Particular case with annual mean insolation
1619	zzz=real(90) ! could be revisited
1620	IF (read_climoz/=-1) THEN
1621	abort_message ='read_climoz=-1 is recommended when solarlong0=1000.'
1622	CALL abort_gcm (modname,abort_message,1)
1623	ENDIF
1624	ELSE
1625	! Case where the season is imposed with solarlong0
1626	zzz=real(90) ! could be revisited
1627	ENDIF
1628	wo(:,:,1)=ozonecm(rlat, paprs,read_climoz,rjour=zzz)
1629	ENDIF
1630	!
1631	! Re-evaporer l'eau liquide nuageuse
1632	!
1633	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1634	DO i = 1, klon
1635	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1636	!jyg<
1637	! Attention : Arnaud a propose des formules completement differentes
1638	! A verifier !!!
1639	zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1640	IF (iflag_ice_thermo .EQ. 0) THEN
1641	zlsdcp=zlvdcp
1642	ENDIF
1643	!>jyg
1644
1645	if (iflag_ice_thermo.eq.0) then
1646	!pas necessaire a priori
1647
1648	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1649	zb = MAX(0.0,ql_seri(i,k))
1650	za = - MAX(0.0,ql_seri(i,k)) &
1651	* (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1652	t_seri(i,k) = t_seri(i,k) + za
1653	q_seri(i,k) = q_seri(i,k) + zb
1654	ql_seri(i,k) = 0.0
1655	d_t_eva(i,k) = za
1656	d_q_eva(i,k) = zb
1657
1658	else
1659
1660	!CR: on ré-évapore eau liquide et glace
1661
1662	! zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1663	! zb = MAX(0.0,ql_seri(i,k))
1664	! za = - MAX(0.0,ql_seri(i,k)) &
1665	! * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1666	zb = MAX(0.0,ql_seri(i,k)+qs_seri(i,k))
1667	za = - MAX(0.0,ql_seri(i,k))*zlvdcp &
1668	- MAX(0.0,qs_seri(i,k))*zlsdcp
1669	t_seri(i,k) = t_seri(i,k) + za
1670	q_seri(i,k) = q_seri(i,k) + zb
1671	ql_seri(i,k) = 0.0
1672	!on évapore la glace
1673	qs_seri(i,k) = 0.0
1674	d_t_eva(i,k) = za
1675	d_q_eva(i,k) = zb
1676	endif
1677
1678	ENDDO
1679	ENDDO
1680	!IM
1681	IF (ip_ebil_phy.ge.2) THEN
1682	ztit='after reevap'
1683	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,1,dtime &
1684	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1685	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1686	call diagphy(airephy,ztit,ip_ebil_phy &
1687	, zero_v, zero_v, zero_v, zero_v, zero_v &
1688	, zero_v, zero_v, zero_v, ztsol &
1689	, d_h_vcol, d_qt, d_ec &
1690	, fs_bound, fq_bound )
1691	!
1692	END IF
1693
1694	!
1695	!=========================================================================
1696	! Calculs de l'orbite.
1697	! Necessaires pour le rayonnement et la surface (calcul de l'albedo).
1698	! doit donc etre plac\'e avant radlwsw et pbl_surface
1699
1700	!!! jyg 17 Sep 2010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1701	call ymds2ju(year_cur, mth_eq, day_eq,0., jD_eq)
1702	day_since_equinox = (jD_cur + jH_cur) - jD_eq
1703	!
1704	! choix entre calcul de la longitude solaire vraie ou valeur fixee a
1705	! solarlong0
1706	if (solarlong0<-999.) then
1707	if (new_orbit) then
1708	! calcul selon la routine utilisee pour les planetes
1709	call solarlong(day_since_equinox, zlongi, dist)
1710	else
1711	! calcul selon la routine utilisee pour l'AR4
1712	CALL orbite(REAL(days_elapsed+1),zlongi,dist)
1713	endif
1714	else
1715	zlongi=solarlong0 ! longitude solaire vraie
1716	dist=1. ! distance au soleil / moyenne
1717	endif
1718	if(prt_level.ge.1) &
1719	write(lunout,*)'Longitude solaire ',zlongi,solarlong0,dist
1720
1721
1722	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1723	! Calcul de l'ensoleillement :
1724	! ============================
1725	! Pour une solarlong0=1000., on calcule un ensoleillement moyen sur
1726	! l'annee a partir d'une formule analytique.
1727	! Cet ensoleillement est sym\'etrique autour de l'\'equateur et
1728	! non nul aux poles.
1729	IF (abs(solarlong0-1000.)<1.e-4) then
1730	call zenang_an(cycle_diurne,jH_cur,rlat,rlon,rmu0,fract)
1731	ELSE
1732	! Avec ou sans cycle diurne
1733	IF (cycle_diurne) THEN
1734	zdtime=dtime*REAL(radpas) ! pas de temps du rayonnement (s)
1735	CALL zenang(zlongi,jH_cur,zdtime,rlat,rlon,rmu0,fract)
1736	ELSE
1737	CALL angle(zlongi, rlat, fract, rmu0)
1738	ENDIF
1739	ENDIF
1740
1741	! AI Janv 2014
1742	do i = 1, klon
1743	if (fract(i).le.0.) then
1744	JrNt(i)=0.
1745	else
1746	JrNt(i)=1.
1747	endif
1748	enddo
1749
1750	if (mydebug) then
1751	call writefield_phy('u_seri',u_seri,llm)
1752	call writefield_phy('v_seri',v_seri,llm)
1753	call writefield_phy('t_seri',t_seri,llm)
1754	call writefield_phy('q_seri',q_seri,llm)
1755	endif
1756
1757	!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1758	! Appel au pbl_surface : Planetary Boudary Layer et Surface
1759	! Cela implique tous les interactions des sous-surfaces et la partie diffusion
1760	! turbulent du couche limit.
1761	!
1762	! Certains varibales de sorties de pbl_surface sont utiliser que pour
1763	! ecriture des fihiers hist_XXXX.nc, ces sont :
1764	! qsol, zq2m, s_pblh, s_lcl,
1765	! s_capCL, s_oliqCL, s_cteiCL,s_pblT,
1766	! s_therm, s_trmb1, s_trmb2, s_trmb3,
1767	! zxrugs, zu10m, zv10m, fder,
1768	! zxqsurf, rh2m, zxfluxu, zxfluxv,
1769	! frugs, agesno, fsollw, fsolsw,
1770	! d_ts, fevap, fluxlat, t2m,
1771	! wfbils, wfbilo, fluxt, fluxu, fluxv,
1772	!
1773	! Certains ne sont pas utiliser du tout :
1774	! dsens, devap, zxsnow, zxfluxt, zxfluxq, q2m, fluxq
1775	!
1776
1777	! Calcul de l'humidite de saturation au niveau du sol
1778
1779
1780
1781	if (iflag_pbl/=0) then
1782
1783	CALL pbl_surface( &
1784	dtime, date0, itap, days_elapsed+1, &
1785	debut, lafin, &
1786	rlon, rlat, rugoro, rmu0, &
1787	zsig, sollwdown, pphi, cldt, &
1788	rain_fall, snow_fall, solsw, sollw, &
1789	t_seri, q_seri, u_seri, v_seri, &
1790	pplay, paprs, pctsrf, &
1791	ftsol,falb1,falb2,ustar,u10m,v10m,wstar, &
1792	sollwdown, cdragh, cdragm, u1, v1, &
1793	albsol1, albsol2, sens, evap, &
1794	albsol3_lic,runoff, snowhgt, qsnow, to_ice, sissnow, &
1795	zxtsol, zxfluxlat, zt2m, qsat2m, &
1796	d_t_vdf, d_q_vdf, d_u_vdf, d_v_vdf, d_t_diss, &
1797	coefh(1:klon,1:klev,1:nbsrf+1), coefm(1:klon,1:klev,1:nbsrf+1), &
1798	slab_wfbils, &
1799	qsol, zq2m, s_pblh, s_lcl, &
1800	s_capCL, s_oliqCL, s_cteiCL,s_pblT, &
1801	s_therm, s_trmb1, s_trmb2, s_trmb3, &
1802	zxrugs, zustar, zu10m, zv10m, fder, &
1803	zxqsurf, rh2m, zxfluxu, zxfluxv, &
1804	frugs, agesno, fsollw, fsolsw, &
1805	d_ts, fevap, fluxlat, t2m, &
1806	wfbils, wfbilo, fluxt, fluxu, fluxv, &
1807	dsens, devap, zxsnow, &
1808	zxfluxt, zxfluxq, q2m, fluxq, pbl_tke )
1809
1810
1811	!---------------------------------------------------------------------
1812	! ajout des tendances de la diffusion turbulente
1813	IF (klon_glo==1) THEN
1814	CALL add_pbl_tend &
1815	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
1816	ELSE
1817	CALL add_phys_tend &
1818	(d_u_vdf,d_v_vdf,d_t_vdf+d_t_diss,d_q_vdf,dql0,dqi0,paprs,'vdf')
1819	ENDIF
1820	!--------------------------------------------------------------------
1821
1822	if (mydebug) then
1823	call writefield_phy('u_seri',u_seri,llm)
1824	call writefield_phy('v_seri',v_seri,llm)
1825	call writefield_phy('t_seri',t_seri,llm)
1826	call writefield_phy('q_seri',q_seri,llm)
1827	endif
1828
1829	CALL evappot(klon,nbsrf,ftsol,pplay(:,1),cdragh, &
1830	t_seri(:,1),q_seri(:,1),u_seri(:,1),v_seri(:,1),evap_pot)
1831
1832
1833	IF (ip_ebil_phy.ge.2) THEN
1834	ztit='after surface_main'
1835	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
1836	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
1837	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1838	call diagphy(airephy,ztit,ip_ebil_phy &
1839	, zero_v, zero_v, zero_v, zero_v, sens &
1840	, evap , zero_v, zero_v, ztsol &
1841	, d_h_vcol, d_qt, d_ec &
1842	, fs_bound, fq_bound )
1843	END IF
1844
1845	ENDIF
1846	! =================================================================== c
1847	! Calcul de Qsat
1848
1849	DO k = 1, klev
1850	DO i = 1, klon
1851	zx_t = t_seri(i,k)
1852	IF (thermcep) THEN
1853	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1854	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1855	zx_qs = MIN(0.5,zx_qs)
1856	zcor = 1./(1.-retv*zx_qs)
1857	zx_qs = zx_qs*zcor
1858	ELSE
1859	IF (zx_t.LT.t_coup) THEN
1860	zx_qs = qsats(zx_t)/pplay(i,k)
1861	ELSE
1862	zx_qs = qsatl(zx_t)/pplay(i,k)
1863	ENDIF
1864	ENDIF
1865	zqsat(i,k)=zx_qs
1866	ENDDO
1867	ENDDO
1868
1869	if (prt_level.ge.1) then
1870	write(lunout,*) 'L qsat (g/kg) avant clouds_gno'
1871	write(lunout,'(i4,f15.4)') (k,1000.*zqsat(igout,k),k=1,klev)
1872	endif
1873	!
1874	! Appeler la convection (au choix)
1875	!
1876	DO k = 1, klev
1877	DO i = 1, klon
1878	conv_q(i,k) = d_q_dyn(i,k) &
1879	+ d_q_vdf(i,k)/dtime
1880	conv_t(i,k) = d_t_dyn(i,k) &
1881	+ d_t_vdf(i,k)/dtime
1882	ENDDO
1883	ENDDO
1884	IF (check) THEN
1885	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
1886	WRITE(lunout,*) "avantcon=", za
1887	ENDIF
1888	zx_ajustq = .FALSE.
1889	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1890	IF (zx_ajustq) THEN
1891	DO i = 1, klon
1892	z_avant(i) = 0.0
1893	ENDDO
1894	DO k = 1, klev
1895	DO i = 1, klon
1896	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k)) &
1897	*(paprs(i,k)-paprs(i,k+1))/RG
1898	ENDDO
1899	ENDDO
1900	ENDIF
1901
1902	! Calcule de vitesse verticale a partir de flux de masse verticale
1903	DO k = 1, klev
1904	DO i = 1, klon
1905	omega(i,k) = RG*flxmass_w(i,k) / airephy(i)
1906	END DO
1907	END DO
1908	if (prt_level.ge.1) write(lunout,*) 'omega(igout, :) = ', &
1909	omega(igout, :)
1910
1911	IF (iflag_con.EQ.1) THEN
1912	abort_message ='reactiver le call conlmd dans physiq.F'
1913	CALL abort_gcm (modname,abort_message,1)
1914	! CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1915	! . d_t_con, d_q_con,
1916	! . rain_con, snow_con, ibas_con, itop_con)
1917	ELSE IF (iflag_con.EQ.2) THEN
1918	CALL conflx(dtime, paprs, pplay, t_seri, q_seri, &
1919	conv_t, conv_q, -evap, omega, &
1920	d_t_con, d_q_con, rain_con, snow_con, &
1921	pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
1922	kcbot, kctop, kdtop, pmflxr, pmflxs)
1923	d_u_con = 0.
1924	d_v_con = 0.
1925
1926	WHERE (rain_con < 0.) rain_con = 0.
1927	WHERE (snow_con < 0.) snow_con = 0.
1928	DO i = 1, klon
1929	ibas_con(i) = klev+1 - kcbot(i)
1930	itop_con(i) = klev+1 - kctop(i)
1931	ENDDO
1932	ELSE IF (iflag_con.GE.3) THEN
1933	! nb of tracers for the KE convection:
1934	! MAF la partie traceurs est faite dans phytrac
1935	! on met ntra=1 pour limiter les appels mais on peut
1936	! supprimer les calculs / ftra.
1937	ntra = 1
1938
1939	!=========================================================================
1940	!ajout pour la parametrisation des poches froides: calcul de
1941	!t_wake et t_undi: si pas de poches froides, t_wake=t_undi=t_seri
1942	do k=1,klev
1943	do i=1,klon
1944	if (iflag_wake>=1) then
1945	t_wake(i,k) = t_seri(i,k) &
1946	+(1-wake_s(i))*wake_deltat(i,k)
1947	q_wake(i,k) = q_seri(i,k) &
1948	+(1-wake_s(i))*wake_deltaq(i,k)
1949	t_undi(i,k) = t_seri(i,k) &
1950	-wake_s(i)*wake_deltat(i,k)
1951	q_undi(i,k) = q_seri(i,k) &
1952	-wake_s(i)*wake_deltaq(i,k)
1953	else
1954	t_wake(i,k) = t_seri(i,k)
1955	q_wake(i,k) = q_seri(i,k)
1956	t_undi(i,k) = t_seri(i,k)
1957	q_undi(i,k) = q_seri(i,k)
1958	endif
1959	enddo
1960	enddo
1961
1962	! Calcul de l'energie disponible ALE (J/kg) et de la puissance
1963	! disponible ALP (W/m2) pour le soulevement des particules dans
1964	! le modele convectif
1965	!
1966	do i = 1,klon
1967	ALE(i) = 0.
1968	ALP(i) = 0.
1969	enddo
1970	!
1971	!calcul de ale_wake et alp_wake
1972	if (iflag_wake>=1) then
1973	if (itap .le. it_wape_prescr) then
1974	do i = 1,klon
1975	ale_wake(i) = wape_prescr
1976	alp_wake(i) = fip_prescr
1977	enddo
1978	else
1979	do i = 1,klon
1980	!jyg ALE=WAPE au lieu de ALE = 1/2 Cstar**2
1981	!cc ale_wake(i) = 0.5wake_cstar(i)*2
1982	ale_wake(i) = wake_pe(i)
1983	alp_wake(i) = wake_fip(i)
1984	enddo
1985	endif
1986	else
1987	do i = 1,klon
1988	ale_wake(i) = 0.
1989	alp_wake(i) = 0.
1990	enddo
1991	endif
1992	!combinaison avec ale et alp de couche limite: constantes si pas
1993	!de couplage, valeurs calculees dans le thermique sinon
1994	if (iflag_coupl.eq.0) then
1995	if (debut.and.prt_level.gt.9) &
1996	WRITE(lunout,*)'ALE et ALP imposes'
1997	do i = 1,klon
1998	!on ne couple que ale
1999	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2000	ALE(i) = max(ale_wake(i),ale_bl_prescr)
2001	!on ne couple que alp
2002	! ALP(i) = alp_wake(i) + Alp_bl(i)
2003	ALP(i) = alp_wake(i) + alp_bl_prescr
2004	enddo
2005	else
2006	IF(prt_level>9)WRITE(lunout,*)'ALE et ALP couples au thermique'
2007	! do i = 1,klon
2008	! ALE(i) = max(ale_wake(i),Ale_bl(i))
2009	! avant ALP(i) = alp_wake(i) + Alp_bl(i)
2010	! ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2011	! write(20,*)'ALE',ALE(i),Ale_bl(i),ale_wake(i)
2012	! write(21,*)'ALP',ALP(i),Alp_bl(i),alp_wake(i)
2013	! enddo
2014
2015	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2016	! Modif FH 2010/04/27. Sans doute temporaire.
2017	! Deux options pour le alp_offset : constant si >?? 0 ou
2018	! proportionnel ??a w si <0
2019	! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2020	! Estimation d'une vitesse verticale effective pour ALP
2021	www(1:klon)=0.
2022	do k=2,klev-1
2023	do i=1,klon
2024	www(i)=max(www(i),-omega(i,k)RDt_seri(i,k)/(RG*paprs(i,k)) &
2025	& zw2(i,k)zw2(i,k))
2026	! if (paprs(i,k)>pbase(i)) then
2027	! calcul approche de la vitesse verticale en m/s
2028	! www(i)=max(www(i),-omega(i,k)RDtemp(i,k)/(RG*paprs(i,k))
2029	! endif
2030	! Le 0.1 est en gros H / ps = 1e5 / 1e4
2031	enddo
2032	enddo
2033	do i=1,klon
2034	if (www(i)>0. .and. ale_bl(i)>0. ) www(i)=www(i)/ale_bl(i)
2035	enddo
2036
2037
2038	do i = 1,klon
2039	ALE(i) = max(ale_wake(i),Ale_bl(i))
2040	!cc nrlmd le 10/04/2012----------Stochastic triggering--------------
2041	if (iflag_trig_bl.ge.1) then
2042	ALE(i) = max(ale_wake(i),Ale_bl_trig(i))
2043	endif
2044	!cc fin nrlmd le 10/04/2012
2045	if (alp_offset>=0.) then
2046	ALP(i) = alp_wake(i) + Alp_bl(i) + alp_offset ! modif sb
2047	else
2048
2049	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2050	! _ _
2051	! Ajout d'une composante 3 * A * w w'2 a w'3 avec w=www : w max sous pbase
2052	! ou A est la fraction couverte par les ascendances w'
2053	! on utilise le fait que A * w'3 = ALP
2054	! et donc A * w'2 ~ ALP / sqrt(ALE) (on ajoute 0.1 pour les
2055	! singularites)
2056	ALP(i)=alp_wake(i)(1.+3.www(i)/( sqrt(ale_wake(i))+0.1) ) &
2057	& +alp_bl(i) (1.+3.www(i)/( sqrt(ale_bl(i)) +0.1) )
2058	! ALP(i)=alp_wake(i)+Alp_bl(i)+alp_offset*min(omega(i,6),0.)
2059	! if (alp(i)<0.) then
2060	! print*,'ALP ',alp(i),alp_wake(i) &
2061	! ,Alp_bl(i),alp_offset*min(omega(i,6),0.)
2062	! endif
2063	endif
2064	enddo
2065	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2066
2067	endif
2068	do i=1,klon
2069	if (alp(i)>alp_max) then
2070	IF(prt_level>9)WRITE(lunout,*) &
2071	'WARNING SUPER ALP (seuil=',alp_max, &
2072	'): i, alp, alp_wake,ale',i,alp(i),alp_wake(i),ale(i)
2073	alp(i)=alp_max
2074	endif
2075	if (ale(i)>ale_max) then
2076	IF(prt_level>9)WRITE(lunout,*) &
2077	'WARNING SUPER ALE (seuil=',ale_max, &
2078	'): i, alp, alp_wake,ale',i,ale(i),ale_wake(i),alp(i)
2079	ale(i)=ale_max
2080	endif
2081	enddo
2082
2083	!fin calcul ale et alp
2084	!=======================================================================
2085
2086
2087	! sb, oct02:
2088	! Schema de convection modularise et vectorise:
2089	! (driver commun aux versions 3 et 4)
2090	!
2091	IF (ok_cvl) THEN ! new driver for convectL
2092
2093	IF (type_trac == 'repr') THEN
2094	nbtr_tmp=ntra
2095	ELSE
2096	nbtr_tmp=nbtr
2097	END IF
2098	!jyg iflag_con est dans clesphys
2099	!c CALL concvl (iflag_con,iflag_clos,
2100	CALL concvl (iflag_clos, &
2101	dtime,paprs,pplay,t_undi,q_undi, &
2102	t_wake,q_wake,wake_s, &
2103	u_seri,v_seri,tr_seri,nbtr_tmp, &
2104	ALE,ALP, &
2105	sig1,w01, &
2106	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2107	rain_con, snow_con, ibas_con, itop_con, sigd, &
2108	ema_cbmf,plcl,plfc,wbeff,upwd,dnwd,dnwd0, &
2109	Ma,mip,Vprecip,cape,cin,tvp,Tconv,iflagctrl, &
2110	pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd, &
2111	! RomP >>>
2112	!! . pmflxr,pmflxs,da,phi,mp,
2113	!! . ftd,fqd,lalim_conv,wght_th)
2114	pmflxr,pmflxs,da,phi,mp,phi2,d1a,dam,sij,clw,elij, &
2115	ftd,fqd,lalim_conv,wght_th, &
2116	ev, ep,epmlmMm,eplaMm, &
2117	wdtrainA,wdtrainM,wght_cvfd)
2118	! RomP <<<
2119
2120	!IM begin
2121	! print*,'physiq: cin pbase dnwd0 ftd fqd ',cin(1),pbase(1),
2122	! .dnwd0(1,1),ftd(1,1),fqd(1,1)
2123	!IM end
2124	!IM cf. FH
2125	clwcon0=qcondc
2126	pmfu(:,:)=upwd(:,:)+dnwd(:,:)
2127
2128	do i = 1, klon
2129	if (iflagctrl(i).le.1) itau_con(i)=itau_con(i)+1
2130	enddo
2131
2132	ELSE ! ok_cvl
2133
2134	! MAF conema3 ne contient pas les traceurs
2135	CALL conema3 (dtime, &
2136	paprs,pplay,t_seri,q_seri, &
2137	u_seri,v_seri,tr_seri,ntra, &
2138	sig1,w01, &
2139	d_t_con,d_q_con,d_u_con,d_v_con,d_tr, &
2140	rain_con, snow_con, ibas_con, itop_con, &
2141	upwd,dnwd,dnwd0,bas,top, &
2142	Ma,cape,tvp,rflag, &
2143	pbase &
2144	,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr &
2145	,clwcon0)
2146
2147	ENDIF ! ok_cvl
2148
2149	!
2150	! Correction precip
2151	rain_con = rain_con * cvl_corr
2152	snow_con = snow_con * cvl_corr
2153	!
2154
2155	IF (.NOT. ok_gust) THEN
2156	do i = 1, klon
2157	wd(i)=0.0
2158	enddo
2159	ENDIF
2160
2161	! =================================================================== c
2162	! Calcul des proprietes des nuages convectifs
2163	!
2164
2165	! calcul des proprietes des nuages convectifs
2166	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
2167	call clouds_gno &
2168	(klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
2169
2170	! =================================================================== c
2171
2172	DO i = 1, klon
2173	itop_con(i) = min(max(itop_con(i),1),klev)
2174	ibas_con(i) = min(max(ibas_con(i),1),itop_con(i))
2175	ENDDO
2176
2177	DO i = 1, klon
2178	ema_pcb(i) = paprs(i,ibas_con(i))
2179	ENDDO
2180	DO i = 1, klon
2181	! L'idicage de itop_con peut cacher un pb potentiel
2182	! FH sous la dictee de JYG, CR
2183	ema_pct(i) = paprs(i,itop_con(i)+1)
2184
2185	if (itop_con(i).gt.klev-3) then
2186	if(prt_level >= 9) then
2187	write(lunout,*)'La convection monte trop haut '
2188	write(lunout,*)'itop_con(,',i,',)=',itop_con(i)
2189	endif
2190	endif
2191	ENDDO
2192	ELSE IF (iflag_con.eq.0) THEN
2193	write(lunout,*) 'On n appelle pas la convection'
2194	clwcon0=0.
2195	rnebcon0=0.
2196	d_t_con=0.
2197	d_q_con=0.
2198	d_u_con=0.
2199	d_v_con=0.
2200	rain_con=0.
2201	snow_con=0.
2202	bas=1
2203	top=1
2204	ELSE
2205	WRITE(lunout,*) "iflag_con non-prevu", iflag_con
2206	call abort_gcm("physiq", "", 1)
2207	ENDIF
2208
2209	! CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
2210	! . d_u_con, d_v_con)
2211
2212	CALL add_phys_tend(d_u_con, d_v_con, d_t_con, d_q_con, dql0, dqi0, paprs, &
2213	'convection')
2214	!----------------------------------------------------------------------------
2215
2216	if (mydebug) then
2217	call writefield_phy('u_seri',u_seri,llm)
2218	call writefield_phy('v_seri',v_seri,llm)
2219	call writefield_phy('t_seri',t_seri,llm)
2220	call writefield_phy('q_seri',q_seri,llm)
2221	endif
2222
2223	!IM
2224	IF (ip_ebil_phy.ge.2) THEN
2225	ztit='after convect'
2226	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2227	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2228	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2229	call diagphy(airephy,ztit,ip_ebil_phy &
2230	, zero_v, zero_v, zero_v, zero_v, zero_v &
2231	, zero_v, rain_con, snow_con, ztsol &
2232	, d_h_vcol, d_qt, d_ec &
2233	, fs_bound, fq_bound )
2234	END IF
2235	!
2236	IF (check) THEN
2237	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2238	WRITE(lunout,*)"aprescon=", za
2239	zx_t = 0.0
2240	za = 0.0
2241	DO i = 1, klon
2242	za = za + airephy(i)/REAL(klon)
2243	zx_t = zx_t + (rain_con(i)+ &
2244	snow_con(i))*airephy(i)/REAL(klon)
2245	ENDDO
2246	zx_t = zx_t/za*dtime
2247	WRITE(lunout,*)"Precip=", zx_t
2248	ENDIF
2249	IF (zx_ajustq) THEN
2250	DO i = 1, klon
2251	z_apres(i) = 0.0
2252	ENDDO
2253	DO k = 1, klev
2254	DO i = 1, klon
2255	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k)) &
2256	*(paprs(i,k)-paprs(i,k+1))/RG
2257	ENDDO
2258	ENDDO
2259	DO i = 1, klon
2260	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime) &
2261	/z_apres(i)
2262	ENDDO
2263	DO k = 1, klev
2264	DO i = 1, klon
2265	IF (z_factor(i).GT.(1.0+1.0E-08) .OR. &
2266	z_factor(i).LT.(1.0-1.0E-08)) THEN
2267	q_seri(i,k) = q_seri(i,k) * z_factor(i)
2268	ENDIF
2269	ENDDO
2270	ENDDO
2271	ENDIF
2272	zx_ajustq=.FALSE.
2273
2274	!
2275	!=============================================================================
2276	!RR:Evolution de la poche froide: on ne fait pas de separation wake/env
2277	!pour la couche limite diffuse pour l instant
2278	!
2279	if (iflag_wake>=1) then
2280	DO k=1,klev
2281	DO i=1,klon
2282	dt_dwn(i,k) = ftd(i,k)
2283	wdt_PBL(i,k) = 0.
2284	dq_dwn(i,k) = fqd(i,k)
2285	wdq_PBL(i,k) = 0.
2286	M_dwn(i,k) = dnwd0(i,k)
2287	M_up(i,k) = upwd(i,k)
2288	dt_a(i,k) = d_t_con(i,k)/dtime - ftd(i,k)
2289	udt_PBL(i,k) = 0.
2290	dq_a(i,k) = d_q_con(i,k)/dtime - fqd(i,k)
2291	udq_PBL(i,k) = 0.
2292	ENDDO
2293	ENDDO
2294
2295	IF (iflag_wake==2) THEN
2296	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2297	DO k = 1,klev
2298	dt_dwn(:,k)= dt_dwn(:,k)+ &
2299	ok_wk_lsp(:)*(d_t_eva(:,k)+d_t_lsc(:,k))/dtime
2300	dq_dwn(:,k)= dq_dwn(:,k)+ &
2301	ok_wk_lsp(:)*(d_q_eva(:,k)+d_q_lsc(:,k))/dtime
2302	ENDDO
2303	ELSEIF (iflag_wake==3) THEN
2304	ok_wk_lsp(:)=max(sign(1.,wake_s(:)-wake_s_min_lsp),0.)
2305	DO k = 1,klev
2306	DO i=1,klon
2307	IF (rneb(i,k)==0.) THEN
2308	! On ne tient compte des tendances qu'en dehors des nuages (c'est �\| dire
2309	! a priri dans une region ou l'eau se reevapore).
2310	dt_dwn(i,k)= dt_dwn(i,k)+ &
2311	ok_wk_lsp(i)*d_t_lsc(i,k)/dtime
2312	dq_dwn(i,k)= dq_dwn(i,k)+ &
2313	ok_wk_lsp(i)*d_q_lsc(i,k)/dtime
2314	ENDIF
2315	ENDDO
2316	ENDDO
2317	ENDIF
2318
2319	!
2320	!calcul caracteristiques de la poche froide
2321	call calWAKE (paprs,pplay,dtime &
2322	,t_seri,q_seri,omega &
2323	,dt_dwn,dq_dwn,M_dwn,M_up &
2324	,dt_a,dq_a,sigd &
2325	,wdt_PBL,wdq_PBL &
2326	,udt_PBL,udq_PBL &
2327	,wake_deltat,wake_deltaq,wake_dth &
2328	,wake_h,wake_s,wake_dens &
2329	,wake_pe,wake_fip,wake_gfl &
2330	,dt_wake,dq_wake &
2331	,wake_k, t_undi,q_undi &
2332	,wake_omgbdth,wake_dp_omgb &
2333	,wake_dtKE,wake_dqKE &
2334	,wake_dtPBL,wake_dqPBL &
2335	,wake_omg,wake_dp_deltomg &
2336	,wake_spread,wake_Cstar,wake_d_deltat_gw &
2337	,wake_ddeltat,wake_ddeltaq)
2338	!
2339	!-------------------------------------------------------------------------
2340	! ajout des tendances des poches froides
2341	! Faire rapidement disparaitre l'ancien dt_wake pour garder un d_t_wake
2342	! coherent avec les autres d_t_...
2343	d_t_wake(:,:)=dt_wake(:,:)*dtime
2344	d_q_wake(:,:)=dq_wake(:,:)*dtime
2345	CALL add_phys_tend(du0,dv0,d_t_wake,d_q_wake,dql0,dqi0,paprs,'wake')
2346	!------------------------------------------------------------------------
2347
2348	endif
2349	!
2350	!===================================================================
2351	!JYG
2352	IF (ip_ebil_phy.ge.2) THEN
2353	ztit='after wake'
2354	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2355	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2356	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2357	call diagphy(airephy,ztit,ip_ebil_phy &
2358	, zero_v, zero_v, zero_v, zero_v, zero_v &
2359	, zero_v, zero_v, zero_v, ztsol &
2360	, d_h_vcol, d_qt, d_ec &
2361	, fs_bound, fq_bound )
2362	END IF
2363
2364	! print*,'apres callwake iflag_cldcon=', iflag_cldcon
2365	!
2366	!===================================================================
2367	! Convection seche (thermiques ou ajustement)
2368	!===================================================================
2369	!
2370	call stratocu_if(klon,klev,pctsrf,paprs, pplay,t_seri &
2371	,seuil_inversion,weak_inversion,dthmin)
2372
2373
2374
2375	d_t_ajsb(:,:)=0.
2376	d_q_ajsb(:,:)=0.
2377	d_t_ajs(:,:)=0.
2378	d_u_ajs(:,:)=0.
2379	d_v_ajs(:,:)=0.
2380	d_q_ajs(:,:)=0.
2381	clwcon0th(:,:)=0.
2382	!
2383	! fm_therm(:,:)=0.
2384	! entr_therm(:,:)=0.
2385	! detr_therm(:,:)=0.
2386	!
2387	IF(prt_level>9)WRITE(lunout,*) &
2388	'AVANT LA CONVECTION SECHE , iflag_thermals=' &
2389	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2390	if(iflag_thermals<0) then
2391	! Rien
2392	! ====
2393	IF(prt_level>9)WRITE(lunout,*)'pas de convection seche'
2394
2395
2396	else
2397
2398	! Thermiques
2399	! ==========
2400	IF(prt_level>9)WRITE(lunout,*)'JUSTE AVANT , iflag_thermals=' &
2401	,iflag_thermals,' nsplit_thermals=',nsplit_thermals
2402
2403
2404	!cc nrlmd le 10/04/2012
2405	DO k=1,klev+1
2406	DO i=1,klon
2407	pbl_tke_input(i,k,is_oce)=pbl_tke(i,k,is_oce)
2408	pbl_tke_input(i,k,is_ter)=pbl_tke(i,k,is_ter)
2409	pbl_tke_input(i,k,is_lic)=pbl_tke(i,k,is_lic)
2410	pbl_tke_input(i,k,is_sic)=pbl_tke(i,k,is_sic)
2411	ENDDO
2412	ENDDO
2413	!cc fin nrlmd le 10/04/2012
2414
2415	if (iflag_thermals>=1) then
2416	call calltherm(pdtphys &
2417	,pplay,paprs,pphi,weak_inversion &
2418	,u_seri,v_seri,t_seri,q_seri,zqsat,debut &
2419	,d_u_ajs,d_v_ajs,d_t_ajs,d_q_ajs &
2420	,fm_therm,entr_therm,detr_therm &
2421	,zqasc,clwcon0th,lmax_th,ratqscth &
2422	,ratqsdiff,zqsatth &
2423	!on rajoute ale et alp, et les caracteristiques de la couche alim
2424	,Ale_bl,Alp_bl,lalim_conv,wght_th, zmax0, f0, zw2,fraca &
2425	,ztv,zpspsk,ztla,zthl &
2426	!cc nrlmd le 10/04/2012
2427	,pbl_tke_input,pctsrf,omega,airephy &
2428	,zlcl_th,fraca0,w0,w_conv,therm_tke_max0,env_tke_max0 &
2429	,n2,s2,ale_bl_stat &
2430	,therm_tke_max,env_tke_max &
2431	,alp_bl_det,alp_bl_fluct_m,alp_bl_fluct_tke &
2432	,alp_bl_conv,alp_bl_stat &
2433	!cc fin nrlmd le 10/04/2012
2434	,zqla,ztva )
2435
2436	!cc nrlmd le 10/04/2012
2437	!-----------Stochastic triggering-----------
2438	if (iflag_trig_bl.ge.1) then
2439	!
2440	IF (prt_level .GE. 10) THEN
2441	print *,'cin, ale_bl_stat, alp_bl_stat ', &
2442	cin, ale_bl_stat, alp_bl_stat
2443	ENDIF
2444
2445
2446	!----Initialisations
2447	do i=1,klon
2448	proba_notrig(i)=1.
2449	random_notrig(i)=1e6ale_bl_stat(i)-int(1e6ale_bl_stat(i))
2450	if ( ale_bl_trig(i) .lt. abs(cin(i))+1.e-10 ) then
2451	tau_trig(i)=tau_trig_shallow
2452	else
2453	tau_trig(i)=tau_trig_deep
2454	endif
2455	enddo
2456	!
2457	IF (prt_level .GE. 10) THEN
2458	print *,'random_notrig, tau_trig ', &
2459	random_notrig, tau_trig
2460	print *,'s_trig,s2,n2 ', &
2461	s_trig,s2,n2
2462	ENDIF
2463
2464	!Option pour re-activer l'ancien calcul de Ale_bl (iflag_trig_bl=2)
2465	IF (iflag_trig_bl.eq.1) then
2466
2467	!----Tirage al\'eatoire et calcul de ale_bl_trig
2468	do i=1,klon
2469	if ( (ale_bl_stat(i) .gt. abs(cin(i))+1.e-10) ) then
2470	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2471	(n2(i)*dtime/tau_trig(i))
2472	! print *, 'proba_notrig(i) ',proba_notrig(i)
2473	if (random_notrig(i) .ge. proba_notrig(i)) then
2474	ale_bl_trig(i)=ale_bl_stat(i)
2475	else
2476	ale_bl_trig(i)=0.
2477	endif
2478	else
2479	proba_notrig(i)=1.
2480	random_notrig(i)=0.
2481	ale_bl_trig(i)=0.
2482	endif
2483	enddo
2484
2485	ELSE IF (iflag_trig_bl.eq.2) then
2486
2487	do i=1,klon
2488	if ( (Ale_bl(i) .gt. abs(cin(i))+1.e-10) ) then
2489	proba_notrig(i)=(1.-exp(-s_trig/s2(i)))** &
2490	(n2(i)*dtime/tau_trig(i))
2491	! print *, 'proba_notrig(i) ',proba_notrig(i)
2492	if (random_notrig(i) .ge. proba_notrig(i)) then
2493	ale_bl_trig(i)=Ale_bl(i)
2494	else
2495	ale_bl_trig(i)=0.
2496	endif
2497	else
2498	proba_notrig(i)=1.
2499	random_notrig(i)=0.
2500	ale_bl_trig(i)=0.
2501	endif
2502	enddo
2503
2504	ENDIF
2505
2506	!
2507	IF (prt_level .GE. 10) THEN
2508	print *,'proba_notrig, ale_bl_trig ', &
2509	proba_notrig, ale_bl_trig
2510	ENDIF
2511
2512	endif !(iflag_trig_bl)
2513
2514	!-----------Statistical closure-----------
2515	if (iflag_clos_bl.eq.1) then
2516
2517	do i=1,klon
2518	!CR: alp probabiliste
2519	if (ale_bl_trig(i).gt.0.) then
2520	alp_bl(i)=alp_bl(i)/(1.-min(proba_notrig(i),0.999))
2521	endif
2522	enddo
2523
2524	else if (iflag_clos_bl.eq.2) then
2525
2526	!CR: alp calculee dans thermcell_main
2527	do i=1,klon
2528	alp_bl(i)=alp_bl_stat(i)
2529	enddo
2530
2531	else
2532
2533	alp_bl_stat(:)=0.
2534
2535	endif !(iflag_clos_bl)
2536
2537	IF (prt_level .GE. 10) THEN
2538	print *,'ale_bl_trig, alp_bl_stat ',ale_bl_trig, alp_bl_stat
2539	ENDIF
2540
2541	!cc fin nrlmd le 10/04/2012
2542
2543	! ----------------------------------------------------------------------
2544	! Transport de la TKE par les panaches thermiques.
2545	! FH : 2010/02/01
2546	! if (iflag_pbl.eq.10) then
2547	! call thermcell_dtke(klon,klev,nbsrf,pdtphys,fm_therm,entr_therm,
2548	! s rg,paprs,pbl_tke)
2549	! endif
2550	! ----------------------------------------------------------------------
2551	!IM/FH: 2011/02/23
2552	! Couplage Thermiques/Emanuel seulement si T<0
2553	if (iflag_coupl==2) then
2554	print*,'Couplage Thermiques/Emanuel seulement si T<0'
2555	do i=1,klon
2556	if (t_seri(i,lmax_th(i))>273.) then
2557	Ale_bl(i)=0.
2558	endif
2559	enddo
2560	endif
2561
2562	do i=1,klon
2563	! zmax_th(i)=pphi(i,lmax_th(i))/rg
2564	!CR:04/05/12:correction calcul zmax
2565	zmax_th(i)=zmax0(i)
2566	enddo
2567
2568	endif
2569
2570
2571	! Ajustement sec
2572	! ==============
2573
2574	! Dans le cas o\`u on active les thermiques, on fait partir l'ajustement
2575	! a partir du sommet des thermiques.
2576	! Dans le cas contraire, on demarre au niveau 1.
2577
2578	if (iflag_thermals>=13.or.iflag_thermals<=0) then
2579
2580	if(iflag_thermals.eq.0) then
2581	IF(prt_level>9)WRITE(lunout,*)'ajsec'
2582	limbas(:)=1
2583	else
2584	limbas(:)=lmax_th(:)
2585	endif
2586
2587	! Attention : le call ajsec_convV2 n'est maintenu que momentanneement
2588	! pour des test de convergence numerique.
2589	! Le nouveau ajsec est a priori mieux, meme pour le cas
2590	! iflag_thermals = 0 (l'ancienne version peut faire des tendances
2591	! non nulles numeriquement pour des mailles non concernees.
2592
2593	if (iflag_thermals==0) then
2594	! Calling adjustment alone (but not the thermal plume model)
2595	CALL ajsec_convV2(paprs, pplay, t_seri,q_seri &
2596	, d_t_ajsb, d_q_ajsb)
2597	else if (iflag_thermals>0) then
2598	! Calling adjustment above the top of thermal plumes
2599	CALL ajsec(paprs, pplay, t_seri,q_seri,limbas &
2600	, d_t_ajsb, d_q_ajsb)
2601	endif
2602
2603	!-----------------------------------------------------------------------
2604	! ajout des tendances de l'ajustement sec ou des thermiques
2605	CALL add_phys_tend(du0,dv0,d_t_ajsb,d_q_ajsb,dql0,dqi0,paprs,'ajsb')
2606	d_t_ajs(:,:)=d_t_ajs(:,:)+d_t_ajsb(:,:)
2607	d_q_ajs(:,:)=d_q_ajs(:,:)+d_q_ajsb(:,:)
2608
2609	!---------------------------------------------------------------------
2610
2611	endif
2612
2613	endif
2614	!
2615	!===================================================================
2616	!IM
2617	IF (ip_ebil_phy.ge.2) THEN
2618	ztit='after dry_adjust'
2619	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2620	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2621	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2622	call diagphy(airephy,ztit,ip_ebil_phy &
2623	, zero_v, zero_v, zero_v, zero_v, zero_v &
2624	, zero_v, zero_v, zero_v, ztsol &
2625	, d_h_vcol, d_qt, d_ec &
2626	, fs_bound, fq_bound )
2627	END IF
2628
2629
2630	!-------------------------------------------------------------------------
2631	! Computation of ratqs, the width (normalized) of the subrid scale
2632	! water distribution
2633	CALL calcratqs(klon,klev,prt_level,lunout, &
2634	iflag_ratqs,iflag_con,iflag_cldcon,pdtphys, &
2635	ratqsbas,ratqshaut,tau_ratqs,fact_cldcon, &
2636	ptconv,ptconvth,clwcon0th, rnebcon0th, &
2637	paprs,pplay,q_seri,zqsat,fm_therm, &
2638	ratqs,ratqsc)
2639
2640
2641	!
2642	! Appeler le processus de condensation a grande echelle
2643	! et le processus de precipitation
2644	!-------------------------------------------------------------------------
2645	IF (prt_level .GE.10) THEN
2646	print *,' ->fisrtilp '
2647	ENDIF
2648	!-------------------------------------------------------------------------
2649	CALL fisrtilp(dtime,paprs,pplay, &
2650	t_seri, q_seri,ptconv,ratqs, &
2651	d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, cldliq, &
2652	rain_lsc, snow_lsc, &
2653	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
2654	frac_impa, frac_nucl, beta_prec_fisrt, &
2655	prfl, psfl, rhcl, &
2656	zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cldcon, &
2657	iflag_ice_thermo)
2658
2659	WHERE (rain_lsc < 0) rain_lsc = 0.
2660	WHERE (snow_lsc < 0) snow_lsc = 0.
2661
2662	CALL add_phys_tend(du0,dv0,d_t_lsc,d_q_lsc,d_ql_lsc,d_qi_lsc,paprs,'lsc')
2663	!---------------------------------------------------------------------------
2664	DO k = 1, klev
2665	DO i = 1, klon
2666	cldfra(i,k) = rneb(i,k)
2667	!CR: a quoi ca sert? Faut-il ajouter qs_seri?
2668	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
2669	ENDDO
2670	ENDDO
2671	IF (check) THEN
2672	za = qcheck(klon,klev,paprs,q_seri,ql_seri,airephy)
2673	WRITE(lunout,*)"apresilp=", za
2674	zx_t = 0.0
2675	za = 0.0
2676	DO i = 1, klon
2677	za = za + airephy(i)/REAL(klon)
2678	zx_t = zx_t + (rain_lsc(i) &
2679	+ snow_lsc(i))*airephy(i)/REAL(klon)
2680	ENDDO
2681	zx_t = zx_t/za*dtime
2682	WRITE(lunout,*)"Precip=", zx_t
2683	ENDIF
2684	!IM
2685	IF (ip_ebil_phy.ge.2) THEN
2686	ztit='after fisrt'
2687	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2688	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2689	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2690	call diagphy(airephy,ztit,ip_ebil_phy &
2691	, zero_v, zero_v, zero_v, zero_v, zero_v &
2692	, zero_v, rain_lsc, snow_lsc, ztsol &
2693	, d_h_vcol, d_qt, d_ec &
2694	, fs_bound, fq_bound )
2695	END IF
2696
2697	if (mydebug) then
2698	call writefield_phy('u_seri',u_seri,llm)
2699	call writefield_phy('v_seri',v_seri,llm)
2700	call writefield_phy('t_seri',t_seri,llm)
2701	call writefield_phy('q_seri',q_seri,llm)
2702	endif
2703
2704	!
2705	!-------------------------------------------------------------------
2706	! PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
2707	!-------------------------------------------------------------------
2708
2709	! 1. NUAGES CONVECTIFS
2710	!
2711	!IM cf FH
2712	! IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
2713	IF (iflag_cldcon.le.-1) THEN ! seulement pour Tiedtke
2714	snow_tiedtke=0.
2715	! print*,'avant calcul de la pseudo precip '
2716	! print*,'iflag_cldcon',iflag_cldcon
2717	if (iflag_cldcon.eq.-1) then
2718	rain_tiedtke=rain_con
2719	else
2720	! print*,'calcul de la pseudo precip '
2721	rain_tiedtke=0.
2722	! print*,'calcul de la pseudo precip 0'
2723	do k=1,klev
2724	do i=1,klon
2725	if (d_q_con(i,k).lt.0.) then
2726	rain_tiedtke(i)=rain_tiedtke(i)-d_q_con(i,k)/pdtphys &
2727	*(paprs(i,k)-paprs(i,k+1))/rg
2728	endif
2729	enddo
2730	enddo
2731	endif
2732	!
2733	! call dump2d(iim,jjm,rain_tiedtke(2:klon-1),'PSEUDO PRECIP ')
2734	!
2735
2736	! Nuages diagnostiques pour Tiedtke
2737	CALL diagcld1(paprs,pplay, &
2738	!IM cf FH . rain_con,snow_con,ibas_con,itop_con,
2739	rain_tiedtke,snow_tiedtke,ibas_con,itop_con, &
2740	diafra,dialiq)
2741	DO k = 1, klev
2742	DO i = 1, klon
2743	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2744	cldliq(i,k) = dialiq(i,k)
2745	cldfra(i,k) = diafra(i,k)
2746	ENDIF
2747	ENDDO
2748	ENDDO
2749
2750	ELSE IF (iflag_cldcon.ge.3) THEN
2751	! On prend pour les nuages convectifs le max du calcul de la
2752	! convection et du calcul du pas de temps precedent diminue d'un facteur
2753	! facttemps
2754	facteur = pdtphys *facttemps
2755	do k=1,klev
2756	do i=1,klon
2757	rnebcon(i,k)=rnebcon(i,k)*facteur
2758	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k)) &
2759	then
2760	rnebcon(i,k)=rnebcon0(i,k)
2761	clwcon(i,k)=clwcon0(i,k)
2762	endif
2763	enddo
2764	enddo
2765
2766	!
2767	!jq - introduce the aerosol direct and first indirect radiative forcings
2768	!jq - Johannes Quaas, 27/11/2003 (quaas@lmd.jussieu.fr)
2769	IF (flag_aerosol .gt. 0) THEN
2770	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
2771	IF (.NOT. aerosol_couple) THEN
2772	!
2773	CALL readaerosol_optic( &
2774	debut, new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
2775	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
2776	mass_solu_aero, mass_solu_aero_pi, &
2777	tau_aero, piz_aero, cg_aero, &
2778	tausum_aero, tau3d_aero)
2779	ENDIF
2780	ELSE ! RRTM radiation
2781	IF (aerosol_couple .AND. config_inca == 'aero' ) THEN
2782	abort_message='config_inca=aero et rrtm=1 impossible'
2783	call abort_gcm(modname,abort_message,1)
2784	ELSE
2785	!
2786	#ifdef CPP_RRTM
2787	CALL readaerosol_optic_rrtm( debut, aerosol_couple, &
2788	new_aod, flag_aerosol, itap, jD_cur-jD_ref, &
2789	pdtphys, pplay, paprs, t_seri, rhcl, presnivs, &
2790	tr_seri, mass_solu_aero, mass_solu_aero_pi, &
2791	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm, &
2792	tausum_aero, tau3d_aero)
2793	#else
2794
2795	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
2796	call abort_gcm(modname,abort_message,1)
2797	#endif
2798	!
2799	ENDIF
2800	ENDIF
2801	ELSE
2802	tausum_aero(:,:,:) = 0.
2803	IF (iflag_rrtm .EQ. 0) THEN !--old radiation
2804	tau_aero(:,:,:,:) = 0.
2805	piz_aero(:,:,:,:) = 0.
2806	cg_aero(:,:,:,:) = 0.
2807	ELSE
2808	tau_aero_sw_rrtm(:,:,:,:)=0.0
2809	piz_aero_sw_rrtm(:,:,:,:)=0.0
2810	cg_aero_sw_rrtm(:,:,:,:)=0.0
2811	ENDIF
2812	ENDIF
2813	!
2814	!--STRAT AEROSOL
2815	!--updates tausum_aero,tau_aero,piz_aero,cg_aero
2816	IF (flag_aerosol_strat) THEN
2817	PRINT *,'appel a readaerosolstrat', mth_cur
2818	IF (iflag_rrtm.EQ.0) THEN
2819	CALL readaerosolstrato(debut)
2820	ELSE
2821	#ifdef CPP_RRTM
2822	CALL readaerosolstrato_rrtm(debut)
2823	#else
2824
2825	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
2826	call abort_gcm(modname,abort_message,1)
2827	#endif
2828	ENDIF
2829	ENDIF
2830	!--fin STRAT AEROSOL
2831
2832
2833	! On prend la somme des fractions nuageuses et des contenus en eau
2834
2835	if (iflag_cldcon>=5) then
2836
2837	do k=1,klev
2838	ptconvth(:,k)=fm_therm(:,k+1)>0.
2839	enddo
2840
2841	if (iflag_coupl==4) then
2842
2843	! Dans le cas iflag_coupl==4, on prend la somme des convertures
2844	! convectives et lsc dans la partie des thermiques
2845	! Le controle par iflag_coupl est peut etre provisoire.
2846	do k=1,klev
2847	do i=1,klon
2848	if (ptconv(i,k).and.ptconvth(i,k)) then
2849	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
2850	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
2851	else if (ptconv(i,k)) then
2852	cldfra(i,k)=rnebcon(i,k)
2853	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
2854	endif
2855	enddo
2856	enddo
2857
2858	else if (iflag_coupl==5) then
2859	do k=1,klev
2860	do i=1,klon
2861	cldfra(i,k)=min(cldfra(i,k)+rnebcon(i,k),1.)
2862	cldliq(i,k)=cldliq(i,k)+rnebcon(i,k)*clwcon(i,k)
2863	enddo
2864	enddo
2865
2866	else
2867
2868	! Si on est sur un point touche par la convection profonde et pas
2869	! par les thermiques, on prend la couverture nuageuse et l'eau nuageuse
2870	! de la convection profonde.
2871
2872	!IM/FH: 2011/02/23
2873	! definition des points sur lesquels ls thermiques sont actifs
2874
2875	do k=1,klev
2876	do i=1,klon
2877	if (ptconv(i,k).and. .not. ptconvth(i,k)) then
2878	cldfra(i,k)=rnebcon(i,k)
2879	cldliq(i,k)=rnebcon(i,k)*clwcon(i,k)
2880	endif
2881	enddo
2882	enddo
2883
2884	endif
2885
2886	else
2887
2888	! Ancienne version
2889	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2890	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2891	endif
2892
2893	ENDIF
2894
2895	! plulsc(:)=0.
2896	! do k=1,klev,-1
2897	! do i=1,klon
2898	! zzz=prfl(:,k)+psfl(:,k)
2899	! if (.not.ptconvth.zzz.gt.0.)
2900	! enddo prfl, psfl,
2901	! enddo
2902	!
2903	! 2. NUAGES STARTIFORMES
2904	!
2905	IF (ok_stratus) THEN
2906	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2907	DO k = 1, klev
2908	DO i = 1, klon
2909	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2910	cldliq(i,k) = dialiq(i,k)
2911	cldfra(i,k) = diafra(i,k)
2912	ENDIF
2913	ENDDO
2914	ENDDO
2915	ENDIF
2916	!
2917	! Precipitation totale
2918	!
2919	DO i = 1, klon
2920	rain_fall(i) = rain_con(i) + rain_lsc(i)
2921	snow_fall(i) = snow_con(i) + snow_lsc(i)
2922	ENDDO
2923	!IM
2924	IF (ip_ebil_phy.ge.2) THEN
2925	ztit="after diagcld"
2926	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
2927	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
2928	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2929	call diagphy(airephy,ztit,ip_ebil_phy &
2930	, zero_v, zero_v, zero_v, zero_v, zero_v &
2931	, zero_v, zero_v, zero_v, ztsol &
2932	, d_h_vcol, d_qt, d_ec &
2933	, fs_bound, fq_bound )
2934	END IF
2935	!
2936	! Calculer l'humidite relative pour diagnostique
2937	!
2938	DO k = 1, klev
2939	DO i = 1, klon
2940	zx_t = t_seri(i,k)
2941	IF (thermcep) THEN
2942	if (iflag_ice_thermo.eq.0) then
2943	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2944	else
2945	zdelta = MAX(0.,SIGN(1.,t_glace_min-zx_t))
2946	endif
2947	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2948	zx_qs = MIN(0.5,zx_qs)
2949	zcor = 1./(1.-retv*zx_qs)
2950	zx_qs = zx_qs*zcor
2951	ELSE
2952	IF (zx_t.LT.t_coup) THEN
2953	zx_qs = qsats(zx_t)/pplay(i,k)
2954	ELSE
2955	zx_qs = qsatl(zx_t)/pplay(i,k)
2956	ENDIF
2957	ENDIF
2958	zx_rh(i,k) = q_seri(i,k)/zx_qs
2959	zqsat(i,k)=zx_qs
2960	ENDDO
2961	ENDDO
2962
2963	!IM Calcul temp.potentielle a 2m (tpot) et temp. potentielle
2964	! equivalente a 2m (tpote) pour diagnostique
2965	!
2966	DO i = 1, klon
2967	tpot(i)=zt2m(i)(100000./paprs(i,1))*RKAPPA
2968	IF (thermcep) THEN
2969	IF(zt2m(i).LT.RTT) then
2970	Lheat=RLSTT
2971	ELSE
2972	Lheat=RLVTT
2973	ENDIF
2974	ELSE
2975	IF (zt2m(i).LT.RTT) THEN
2976	Lheat=RLSTT
2977	ELSE
2978	Lheat=RLVTT
2979	ENDIF
2980	ENDIF
2981	tpote(i) = tpot(i)* &
2982	EXP((Lheat qsat2m(i))/(RCPDzt2m(i)))
2983	ENDDO
2984
2985	IF (type_trac == 'inca') THEN
2986	#ifdef INCA
2987	CALL VTe(VTphysiq)
2988	CALL VTb(VTinca)
2989	calday = REAL(days_elapsed + 1) + jH_cur
2990
2991	call chemtime(itap+itau_phy-1, date0, dtime)
2992	IF (config_inca == 'aero' .OR. config_inca == 'aeNP') THEN
2993	CALL AEROSOL_METEO_CALC( &
2994	calday,pdtphys,pplay,paprs,t,pmflxr,pmflxs, &
2995	prfl,psfl,pctsrf,airephy,rlat,rlon,u10m,v10m)
2996	END IF
2997
2998	zxsnow_dummy(:) = 0.0
2999
3000	CALL chemhook_begin (calday, &
3001	days_elapsed+1, &
3002	jH_cur, &
3003	pctsrf(1,1), &
3004	rlat, &
3005	rlon, &
3006	airephy, &
3007	paprs, &
3008	pplay, &
3009	coefh(1:klon,1:klev,is_ave), &
3010	pphi, &
3011	t_seri, &
3012	u, &
3013	v, &
3014	wo(:, :, 1), &
3015	q_seri, &
3016	zxtsol, &
3017	zxsnow_dummy, &
3018	solsw, &
3019	albsol1, &
3020	rain_fall, &
3021	snow_fall, &
3022	itop_con, &
3023	ibas_con, &
3024	cldfra, &
3025	iim, &
3026	jjm, &
3027	tr_seri, &
3028	ftsol, &
3029	paprs, &
3030	cdragh, &
3031	cdragm, &
3032	pctsrf, &
3033	pdtphys, &
3034	itap)
3035
3036	CALL VTe(VTinca)
3037	CALL VTb(VTphysiq)
3038	#endif
3039	END IF !type_trac = inca
3040	!
3041	! Calculer les parametres optiques des nuages et quelques
3042	! parametres pour diagnostiques:
3043	!
3044
3045	IF (aerosol_couple) THEN
3046	mass_solu_aero(:,:) = ccm(:,:,1)
3047	mass_solu_aero_pi(:,:) = ccm(:,:,2)
3048	END IF
3049
3050	if (ok_newmicro) then
3051	IF (iflag_rrtm.NE.0) THEN
3052	#ifdef CPP_RRTM
3053	IF (ok_cdnc.AND.NRADLP.NE.3) THEN
3054	abort_message='RRTM choix incoherent NRADLP doit etre egal a 3 pour ok_cdnc'
3055	call abort_gcm(modname,abort_message,1)
3056	endif
3057	#else
3058
3059	abort_message='You should compile with -rrtm if running with iflag_rrtm=1'
3060	call abort_gcm(modname,abort_message,1)
3061	#endif
3062	ENDIF
3063	CALL newmicro (ok_cdnc, bl95_b0, bl95_b1, &
3064	paprs, pplay, t_seri, cldliq, cldfra, &
3065	cldtau, cldemi, cldh, cldl, cldm, cldt, cldq, &
3066	flwp, fiwp, flwc, fiwc, &
3067	mass_solu_aero, mass_solu_aero_pi, &
3068	cldtaupi, re, fl, ref_liq, ref_ice, &
3069	ref_liq_pi, ref_ice_pi)
3070	else
3071	CALL nuage (paprs, pplay, &
3072	t_seri, cldliq, cldfra, cldtau, cldemi, &
3073	cldh, cldl, cldm, cldt, cldq, &
3074	ok_aie, &
3075	mass_solu_aero, mass_solu_aero_pi, &
3076	bl95_b0, bl95_b1, &
3077	cldtaupi, re, fl)
3078	endif
3079	!
3080	!IM betaCRF
3081	!
3082	cldtaurad = cldtau
3083	cldtaupirad = cldtaupi
3084	cldemirad = cldemi
3085
3086	!
3087	if(lon1_beta.EQ.-180..AND.lon2_beta.EQ.180..AND. &
3088	lat1_beta.EQ.90..AND.lat2_beta.EQ.-90.) THEN
3089	!
3090	! global
3091	!
3092	DO k=1, klev
3093	DO i=1, klon
3094	if (pplay(i,k).GE.pfree) THEN
3095	beta(i,k) = beta_pbl
3096	else
3097	beta(i,k) = beta_free
3098	endif
3099	if (mskocean_beta) THEN
3100	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3101	endif
3102	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3103	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3104	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3105	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3106	ENDDO
3107	ENDDO
3108	!
3109	else
3110	!
3111	! regional
3112	!
3113	DO k=1, klev
3114	DO i=1,klon
3115	!
3116	if (rlon(i).ge.lon1_beta.AND.rlon(i).le.lon2_beta.AND. &
3117	rlat(i).le.lat1_beta.AND.rlat(i).ge.lat2_beta) THEN
3118	if (pplay(i,k).GE.pfree) THEN
3119	beta(i,k) = beta_pbl
3120	else
3121	beta(i,k) = beta_free
3122	endif
3123	if (mskocean_beta) THEN
3124	beta(i,k) = beta(i,k) * pctsrf(i,is_oce)
3125	endif
3126	cldtaurad(i,k) = cldtau(i,k) * beta(i,k)
3127	cldtaupirad(i,k) = cldtaupi(i,k) * beta(i,k)
3128	cldemirad(i,k) = cldemi(i,k) * beta(i,k)
3129	cldfrarad(i,k) = cldfra(i,k) * beta(i,k)
3130	endif
3131	!
3132	ENDDO
3133	ENDDO
3134	!
3135	endif
3136	!
3137	! Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
3138	!
3139	IF (MOD(itaprad,radpas).EQ.0) THEN
3140
3141	DO i = 1, klon
3142	albsol1(i) = falb1(i,is_oce) * pctsrf(i,is_oce) &
3143	+ falb1(i,is_lic) * pctsrf(i,is_lic) &
3144	+ falb1(i,is_ter) * pctsrf(i,is_ter) &
3145	+ falb1(i,is_sic) * pctsrf(i,is_sic)
3146	albsol2(i) = falb2(i,is_oce) * pctsrf(i,is_oce) &
3147	+ falb2(i,is_lic) * pctsrf(i,is_lic) &
3148	+ falb2(i,is_ter) * pctsrf(i,is_ter) &
3149	+ falb2(i,is_sic) * pctsrf(i,is_sic)
3150	ENDDO
3151
3152	if (mydebug) then
3153	call writefield_phy('u_seri',u_seri,llm)
3154	call writefield_phy('v_seri',v_seri,llm)
3155	call writefield_phy('t_seri',t_seri,llm)
3156	call writefield_phy('q_seri',q_seri,llm)
3157	endif
3158
3159	IF (aerosol_couple) THEN
3160	#ifdef INCA
3161	CALL radlwsw_inca &
3162	(kdlon,kflev,dist, rmu0, fract, solaire, &
3163	paprs, pplay,zxtsol,albsol1, albsol2, t_seri,q_seri, &
3164	wo(:, :, 1), &
3165	cldfrarad, cldemirad, cldtaurad, &
3166	heat,heat0,cool,cool0,radsol,albpla, &
3167	topsw,toplw,solsw,sollw, &
3168	sollwdown, &
3169	topsw0,toplw0,solsw0,sollw0, &
3170	lwdn0, lwdn, lwup0, lwup, &
3171	swdn0, swdn, swup0, swup, &
3172	ok_ade, ok_aie, &
3173	tau_aero, piz_aero, cg_aero, &
3174	topswad_aero, solswad_aero, &
3175	topswad0_aero, solswad0_aero, &
3176	topsw_aero, topsw0_aero, &
3177	solsw_aero, solsw0_aero, &
3178	cldtaupirad, &
3179	topswai_aero, solswai_aero)
3180
3181	#endif
3182	ELSE
3183	!
3184	!IM calcul radiatif pour le cas actuel
3185	!
3186	RCO2 = RCO2_act
3187	RCH4 = RCH4_act
3188	RN2O = RN2O_act
3189	RCFC11 = RCFC11_act
3190	RCFC12 = RCFC12_act
3191	!
3192	IF (prt_level .GE.10) THEN
3193	print *,' ->radlwsw, number 1 '
3194	ENDIF
3195	!
3196	CALL radlwsw &
3197	(dist, rmu0, fract, &
3198	paprs, pplay,zxtsol,albsol1, albsol2, &
3199	t_seri,q_seri,wo, &
3200	cldfrarad, cldemirad, cldtaurad, &
3201	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3202	flag_aerosol_strat, &
3203	tau_aero, piz_aero, cg_aero, &
3204	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! Rajoute par OB pour RRTM
3205	tau_aero_lw_rrtm, &
3206	cldtaupirad,new_aod, &
3207	zqsat, flwc, fiwc, &
3208	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3209	heat,heat0,cool,cool0,radsol,albpla, &
3210	topsw,toplw,solsw,sollw, &
3211	sollwdown, &
3212	topsw0,toplw0,solsw0,sollw0, &
3213	lwdn0, lwdn, lwup0, lwup, &
3214	swdn0, swdn, swup0, swup, &
3215	topswad_aero, solswad_aero, &
3216	topswai_aero, solswai_aero, &
3217	topswad0_aero, solswad0_aero, &
3218	topsw_aero, topsw0_aero, &
3219	solsw_aero, solsw0_aero, &
3220	topswcf_aero, solswcf_aero, &
3221	!-C. Kleinschmitt for LW diagnostics
3222	toplwad_aero, sollwad_aero,&
3223	toplwai_aero, sollwai_aero, &
3224	toplwad0_aero, sollwad0_aero,&
3225	!-end
3226	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3227	ZSWFT0_i, ZFSDN0, ZFSUP0)
3228
3229	!
3230	!IM 2eme calcul radiatif pour le cas perturbe ou au moins un
3231	!IM des taux doit etre different du taux actuel
3232	!IM Par defaut on a les taux perturbes egaux aux taux actuels
3233	!
3234	if (ok_4xCO2atm) then
3235	if (RCO2_per.NE.RCO2_act.OR.RCH4_per.NE.RCH4_act.OR. &
3236	RN2O_per.NE.RN2O_act.OR.RCFC11_per.NE.RCFC11_act.OR. &
3237	RCFC12_per.NE.RCFC12_act) THEN
3238	!
3239	RCO2 = RCO2_per
3240	RCH4 = RCH4_per
3241	RN2O = RN2O_per
3242	RCFC11 = RCFC11_per
3243	RCFC12 = RCFC12_per
3244	!
3245	IF (prt_level .GE.10) THEN
3246	print *,' ->radlwsw, number 2 '
3247	ENDIF
3248	!
3249	CALL radlwsw &
3250	(dist, rmu0, fract, &
3251	paprs, pplay,zxtsol,albsol1, albsol2, &
3252	t_seri,q_seri,wo, &
3253	cldfra, cldemi, cldtau, &
3254	ok_ade.OR.flag_aerosol_strat, ok_aie, flag_aerosol, &
3255	flag_aerosol_strat, &
3256	tau_aero, piz_aero, cg_aero, &
3257	tau_aero_sw_rrtm, piz_aero_sw_rrtm, cg_aero_sw_rrtm,& ! Rajoute par OB pour RRTM
3258	tau_aero_lw_rrtm, &
3259	cldtaupi,new_aod, &
3260	zqsat, flwc, fiwc, &
3261	ref_liq, ref_ice, ref_liq_pi, ref_ice_pi, &
3262	heatp,heat0p,coolp,cool0p,radsolp,albplap, &
3263	topswp,toplwp,solswp,sollwp, &
3264	sollwdownp, &
3265	topsw0p,toplw0p,solsw0p,sollw0p, &
3266	lwdn0p, lwdnp, lwup0p, lwupp, &
3267	swdn0p, swdnp, swup0p, swupp, &
3268	topswad_aerop, solswad_aerop, &
3269	topswai_aerop, solswai_aerop, &
3270	topswad0_aerop, solswad0_aerop, &
3271	topsw_aerop, topsw0_aerop, &
3272	solsw_aerop, solsw0_aerop, &
3273	topswcf_aerop, solswcf_aerop, &
3274	!-C. Kleinschmitt for LW diagnostics
3275	toplwad_aerop, sollwad_aerop,&
3276	toplwai_aerop, sollwai_aerop, &
3277	toplwad0_aerop, sollwad0_aerop,&
3278	!-end
3279	ZLWFT0_i, ZFLDN0, ZFLUP0, &
3280	ZSWFT0_i, ZFSDN0, ZFSUP0)
3281	endif
3282	endif
3283	!
3284	ENDIF ! aerosol_couple
3285	itaprad = 0
3286	ENDIF ! MOD(itaprad,radpas)
3287	itaprad = itaprad + 1
3288
3289	IF (iflag_radia.eq.0) THEN
3290	IF (prt_level.ge.9) THEN
3291	PRINT *,'--------------------------------------------------'
3292	PRINT *,'>>>> ATTENTION rayonnement desactive pour ce cas'
3293	PRINT *,'>>>> heat et cool mis a zero '
3294	PRINT *,'--------------------------------------------------'
3295	END IF
3296	heat=0.
3297	cool=0.
3298	sollw=0. ! MPL 01032011
3299	solsw=0.
3300	radsol=0.
3301	swup=0. ! MPL 27102011 pour les fichiers AMMA_profiles et AMMA_scalars
3302	swup0=0.
3303	swdn=0.
3304	swdn0=0.
3305	lwup=0.
3306	lwup0=0.
3307	lwdn=0.
3308	lwdn0=0.
3309	END IF
3310
3311	!
3312	! Ajouter la tendance des rayonnements (tous les pas)
3313	!
3314	DO k = 1, klev
3315	DO i = 1, klon
3316	t_seri(i,k) = t_seri(i,k) &
3317	+ (heat(i,k)-cool(i,k)) * dtime/RDAY
3318	ENDDO
3319	ENDDO
3320	!
3321	if (mydebug) then
3322	call writefield_phy('u_seri',u_seri,llm)
3323	call writefield_phy('v_seri',v_seri,llm)
3324	call writefield_phy('t_seri',t_seri,llm)
3325	call writefield_phy('q_seri',q_seri,llm)
3326	endif
3327
3328	!IM
3329	IF (ip_ebil_phy.ge.2) THEN
3330	ztit='after rad'
3331	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3332	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3333	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3334	call diagphy(airephy,ztit,ip_ebil_phy &
3335	, topsw, toplw, solsw, sollw, zero_v &
3336	, zero_v, zero_v, zero_v, ztsol &
3337	, d_h_vcol, d_qt, d_ec &
3338	, fs_bound, fq_bound )
3339	END IF
3340	!
3341	!
3342	! Calculer l'hydrologie de la surface
3343	!
3344	! CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
3345	! . agesno, ftsol,fqsurf,fsnow, ruis)
3346	!
3347
3348	!
3349	! Calculer le bilan du sol et la derive de temperature (couplage)
3350	!
3351	DO i = 1, klon
3352	! bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
3353	! a la demande de JLD
3354	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
3355	ENDDO
3356	!
3357	!moddeblott(jan95)
3358	! Appeler le programme de parametrisation de l'orographie
3359	! a l'echelle sous-maille:
3360	!
3361	IF (prt_level .GE.10) THEN
3362	print *,' call orography ? ', ok_orodr
3363	ENDIF
3364	!
3365	IF (ok_orodr) THEN
3366	!
3367	! selection des points pour lesquels le shema est actif:
3368	igwd=0
3369	DO i=1,klon
3370	itest(i)=0
3371	! IF ((zstd(i).gt.10.0)) THEN
3372	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
3373	itest(i)=1
3374	igwd=igwd+1
3375	idx(igwd)=i
3376	ENDIF
3377	ENDDO
3378	! igwdim=MAX(1,igwd)
3379	!
3380	IF (ok_strato) THEN
3381
3382	CALL drag_noro_strato(klon,klev,dtime,paprs,pplay, &
3383	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3384	igwd,idx,itest, &
3385	t_seri, u_seri, v_seri, &
3386	zulow, zvlow, zustrdr, zvstrdr, &
3387	d_t_oro, d_u_oro, d_v_oro)
3388
3389	ELSE
3390	CALL drag_noro(klon,klev,dtime,paprs,pplay, &
3391	zmea,zstd, zsig, zgam, zthe,zpic,zval, &
3392	igwd,idx,itest, &
3393	t_seri, u_seri, v_seri, &
3394	zulow, zvlow, zustrdr, zvstrdr, &
3395	d_t_oro, d_u_oro, d_v_oro)
3396	ENDIF
3397	!
3398	! ajout des tendances
3399	!-----------------------------------------------------------------------------------------
3400	! ajout des tendances de la trainee de l'orographie
3401	CALL add_phys_tend(d_u_oro,d_v_oro,d_t_oro,dq0,dql0,dqi0,paprs,'oro')
3402	!-----------------------------------------------------------------------------------------
3403	!
3404	ENDIF ! fin de test sur ok_orodr
3405	!
3406	if (mydebug) then
3407	call writefield_phy('u_seri',u_seri,llm)
3408	call writefield_phy('v_seri',v_seri,llm)
3409	call writefield_phy('t_seri',t_seri,llm)
3410	call writefield_phy('q_seri',q_seri,llm)
3411	endif
3412
3413	IF (ok_orolf) THEN
3414	!
3415	! selection des points pour lesquels le shema est actif:
3416	igwd=0
3417	DO i=1,klon
3418	itest(i)=0
3419	IF ((zpic(i)-zmea(i)).GT.100.) THEN
3420	itest(i)=1
3421	igwd=igwd+1
3422	idx(igwd)=i
3423	ENDIF
3424	ENDDO
3425	! igwdim=MAX(1,igwd)
3426	!
3427	IF (ok_strato) THEN
3428
3429	CALL lift_noro_strato(klon,klev,dtime,paprs,pplay, &
3430	rlat,zmea,zstd,zpic,zgam,zthe,zpic,zval, &
3431	igwd,idx,itest, &
3432	t_seri, u_seri, v_seri, &
3433	zulow, zvlow, zustrli, zvstrli, &
3434	d_t_lif, d_u_lif, d_v_lif )
3435
3436	ELSE
3437	CALL lift_noro(klon,klev,dtime,paprs,pplay, &
3438	rlat,zmea,zstd,zpic, &
3439	itest, &
3440	t_seri, u_seri, v_seri, &
3441	zulow, zvlow, zustrli, zvstrli, &
3442	d_t_lif, d_u_lif, d_v_lif)
3443	ENDIF
3444	!
3445	!-----------------------------------------------------------------------------------------
3446	! ajout des tendances de la portance de l'orographie
3447	CALL add_phys_tend(d_u_lif,d_v_lif,d_t_lif,dq0,dql0,dqi0,paprs,'lif')
3448	!-----------------------------------------------------------------------------------------
3449	!
3450	ENDIF ! fin de test sur ok_orolf
3451	! HINES GWD PARAMETRIZATION
3452
3453	IF (ok_hines) then
3454
3455	CALL hines_gwd(klon,klev,dtime,paprs,pplay, &
3456	rlat,t_seri,u_seri,v_seri, &
3457	zustrhi,zvstrhi, &
3458	d_t_hin, d_u_hin, d_v_hin)
3459	!
3460	! ajout des tendances
3461	CALL add_phys_tend(d_u_hin,d_v_hin,d_t_hin,dq0,dql0,dqi0,paprs,'hin')
3462
3463	ENDIF
3464
3465	if (ok_gwd_rando) then
3466	call FLOTT_GWD_rando(DTIME, pplay, t_seri, u_seri, v_seri, &
3467	rain_fall + snow_fall, zustr_gwd_rando, zvstr_gwd_rando, &
3468	du_gwd_rando, dv_gwd_rando)
3469	CALL add_phys_tend(du_gwd_rando, dv_gwd_rando, dt0, dq0, dql0,dqi0,paprs, &
3470	'flott_gwd_rando')
3471	end if
3472
3473	! STRESS NECESSAIRES: TOUTE LA PHYSIQUE
3474
3475	if (mydebug) then
3476	call writefield_phy('u_seri',u_seri,llm)
3477	call writefield_phy('v_seri',v_seri,llm)
3478	call writefield_phy('t_seri',t_seri,llm)
3479	call writefield_phy('q_seri',q_seri,llm)
3480	endif
3481
3482	DO i = 1, klon
3483	zustrph(i)=0.
3484	zvstrph(i)=0.
3485	ENDDO
3486	DO k = 1, klev
3487	DO i = 1, klon
3488	zustrph(i)=zustrph(i)+(u_seri(i,k)-u(i,k))/dtime* &
3489	(paprs(i,k)-paprs(i,k+1))/rg
3490	zvstrph(i)=zvstrph(i)+(v_seri(i,k)-v(i,k))/dtime* &
3491	(paprs(i,k)-paprs(i,k+1))/rg
3492	ENDDO
3493	ENDDO
3494	!
3495	!IM calcul composantes axiales du moment angulaire et couple des montagnes
3496	!
3497	IF (is_sequential .and. ok_orodr) THEN
3498	CALL aaam_bud (27,klon,klev,jD_cur-jD_ref,jH_cur, &
3499	ra,rg,romega, &
3500	rlat,rlon,pphis, &
3501	zustrdr,zustrli,zustrph, &
3502	zvstrdr,zvstrli,zvstrph, &
3503	paprs,u,v, &
3504	aam, torsfc)
3505	ENDIF
3506	!IM cf. FLott END
3507	!IM
3508	IF (ip_ebil_phy.ge.2) THEN
3509	ztit='after orography'
3510	CALL diagetpq(airephy,ztit,ip_ebil_phy,2,2,dtime &
3511	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3512	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3513	call diagphy(airephy,ztit,ip_ebil_phy &
3514	, zero_v, zero_v, zero_v, zero_v, zero_v &
3515	, zero_v, zero_v, zero_v, ztsol &
3516	, d_h_vcol, d_qt, d_ec &
3517	, fs_bound, fq_bound )
3518	END IF
3519
3520	!DC Calcul de la tendance due au methane
3521	IF(ok_qch4) THEN
3522	CALL METHOX(1,klon,klon,klev,q_seri,d_q_ch4,pplay)
3523	! ajout de la tendance d'humidite due au methane
3524	CALL add_phys_tend(du0,dv0,dt0,d_q_ch4*dtime,dql0,'q_ch4')
3525	END IF
3526	!
3527	!
3528	!====================================================================
3529	! Interface Simulateur COSP (Calipso, ISCCP, MISR, ..)
3530	!====================================================================
3531	! Abderrahmane 24.08.09
3532
3533	IF (ok_cosp) THEN
3534	! adeclarer
3535	#ifdef CPP_COSP
3536	IF (itap.eq.1.or.MOD(itap,NINT(freq_cosp/dtime)).EQ.0) THEN
3537
3538	print*,'freq_cosp',freq_cosp
3539	mr_ozone=wo(:, :, 1) * dobson_u * 1e3 / zmasse
3540	! print*,'Dans physiq.F avant appel cosp ref_liq,ref_ice=',
3541	! s ref_liq,ref_ice
3542	call phys_cosp(itap,dtime,freq_cosp, &
3543	ok_mensuelCOSP,ok_journeCOSP,ok_hfCOSP, &
3544	ecrit_mth,ecrit_day,ecrit_hf, ok_all_xml, &
3545	klon,klev,rlon,rlat,presnivs,overlap, &
3546	fract,ref_liq,ref_ice, &
3547	pctsrf(:,is_ter)+pctsrf(:,is_lic), &
3548	zu10m,zv10m,pphis, &
3549	zphi,paprs(:,1:klev),pplay,zxtsol,t_seri, &
3550	qx(:,:,ivap),zx_rh,cldfra,rnebcon,flwc,fiwc, &
3551	prfl(:,1:klev),psfl(:,1:klev), &
3552	pmflxr(:,1:klev),pmflxs(:,1:klev), &
3553	mr_ozone,cldtau, cldemi)
3554
3555	! L calipso2D,calipso3D,cfadlidar,parasolrefl,atb,betamol,
3556	! L cfaddbze,clcalipso2,dbze,cltlidarradar,
3557	! M clMISR,
3558	! R clisccp2,boxtauisccp,boxptopisccp,tclisccp,ctpisccp,
3559	! I tauisccp,albisccp,meantbisccp,meantbclrisccp)
3560
3561	ENDIF
3562
3563	#endif
3564	ENDIF !ok_cosp
3565	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3566	!AA
3567	!AA Installation de l'interface online-offline pour traceurs
3568	!AA
3569	!====================================================================
3570	! Calcul des tendances traceurs
3571	!====================================================================
3572	!
3573
3574	IF (type_trac=='repr') THEN
3575	sh_in(:,:) = q_seri(:,:)
3576	ELSE
3577	sh_in(:,:) = qx(:,:,ivap)
3578	END IF
3579
3580	call phytrac ( &
3581	itap, days_elapsed+1, jH_cur, debut, &
3582	lafin, dtime, u, v, t, &
3583	paprs, pplay, pmfu, pmfd, &
3584	pen_u, pde_u, pen_d, pde_d, &
3585	cdragh, coefh(1:klon,1:klev,is_ave), fm_therm, entr_therm, &
3586	u1, v1, ftsol, pctsrf, &
3587	zustar, zu10m, zv10m, &
3588	wstar(:,is_ave), ale_bl, ale_wake, &
3589	rlat, rlon, &
3590	frac_impa,frac_nucl, beta_prec_fisrt,beta_prec, &
3591	presnivs, pphis, pphi, albsol1, &
3592	sh_in, rhcl, cldfra, rneb, &
3593	diafra, cldliq, itop_con, ibas_con, &
3594	pmflxr, pmflxs, prfl, psfl, &
3595	da, phi, mp, upwd, &
3596	phi2, d1a, dam, sij, wght_cvfd, & !<<RomP+RL
3597	wdtrainA, wdtrainM, sigd, clw,elij, & !<<RomP
3598	ev, ep, epmlmMm, eplaMm, & !<<RomP
3599	dnwd, aerosol_couple, flxmass_w, &
3600	tau_aero, piz_aero, cg_aero, ccm, &
3601	rfname, &
3602	d_tr_dyn, & !<<RomP
3603	tr_seri)
3604
3605	IF (offline) THEN
3606
3607	IF (prt_level.ge.9) &
3608	print*,'Attention on met a 0 les thermiques pour phystoke'
3609	call phystokenc ( &
3610	nlon,klev,pdtphys,rlon,rlat, &
3611	t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
3612	fm_therm,entr_therm, &
3613	cdragh,coefh(1:klon,1:klev,is_ave),u1,v1,ftsol,pctsrf, &
3614	frac_impa, frac_nucl, &
3615	pphis,airephy,dtime,itap, &
3616	qx(:,:,ivap),da,phi,mp,upwd,dnwd)
3617
3618
3619	ENDIF
3620
3621	!
3622	! Calculer le transport de l'eau et de l'energie (diagnostique)
3623	!
3624	CALL transp (paprs,zxtsol, &
3625	t_seri, q_seri, u_seri, v_seri, zphi, &
3626	ve, vq, ue, uq)
3627	!
3628	!IM global posePB BEG
3629	IF(1.EQ.0) THEN
3630	!
3631	CALL transp_lay (paprs,zxtsol, &
3632	t_seri, q_seri, u_seri, v_seri, zphi, &
3633	ve_lay, vq_lay, ue_lay, uq_lay)
3634	!
3635	ENDIF !(1.EQ.0) THEN
3636	!IM global posePB END
3637	! Accumuler les variables a stocker dans les fichiers histoire:
3638	!
3639
3640	!================================================================
3641	! Conversion of kinetic and potential energy into heat, for
3642	! parameterisation of subgrid-scale motions
3643	!================================================================
3644
3645	d_t_ec(:,:)=0.
3646	forall (k=1: llm) exner(:, k) = (pplay(:, k)/paprs(:,1))**RKAPPA
3647	CALL ener_conserv(klon,klev,pdtphys,u,v,t,qx(:,:,ivap), &
3648	u_seri,v_seri,t_seri,q_seri,pbl_tke(:,:,is_ave)-tke0(:,:), &
3649	zmasse,exner,d_t_ec)
3650	t_seri(:,:)=t_seri(:,:)+d_t_ec(:,:)
3651
3652	!IM
3653	IF (ip_ebil_phy.ge.1) THEN
3654	ztit='after physic'
3655	CALL diagetpq(airephy,ztit,ip_ebil_phy,1,1,dtime &
3656	, t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay &
3657	, d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
3658	! Comme les tendances de la physique sont ajoute dans la dynamique,
3659	! on devrait avoir que la variation d'entalpie par la dynamique
3660	! est egale a la variation de la physique au pas de temps precedent.
3661	! Donc la somme de ces 2 variations devrait etre nulle.
3662
3663	call diagphy(airephy,ztit,ip_ebil_phy &
3664	, topsw, toplw, solsw, sollw, sens &
3665	, evap, rain_fall, snow_fall, ztsol &
3666	, d_h_vcol, d_qt, d_ec &
3667	, fs_bound, fq_bound )
3668	!
3669	d_h_vcol_phy=d_h_vcol
3670	!
3671	END IF
3672	!
3673	!=======================================================================
3674	! SORTIES
3675	!=======================================================================
3676	!
3677	!IM initialisation + calculs divers diag AMIP2
3678	!
3679	include "calcul_divers.h"
3680	!
3681	!IM Interpolation sur les niveaux de pression du NMC
3682	! -------------------------------------------------
3683	#ifdef CPP_XIOS
3684	!$OMP MASTER
3685	!On recupere la valeur de la missing value donnee dans le xml
3686	CALL xios_get_field_attr("t850",default_value=missing_val_omp)
3687	! PRINT *,"ARNAUD value missing ",missing_val_omp
3688	!$OMP END MASTER
3689	!$OMP BARRIER
3690	missing_val=missing_val_omp
3691	#endif
3692	#ifndef CPP_XIOS
3693	missing_val=missing_val_nf90
3694	#endif
3695	!
3696	include "calcul_STDlev.h"
3697	!
3698	! slp sea level pressure
3699	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
3700	!
3701	!cc prw = eau precipitable
3702	DO i = 1, klon
3703	prw(i) = 0.
3704	DO k = 1, klev
3705	prw(i) = prw(i) + &
3706	q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
3707	ENDDO
3708	ENDDO
3709	!
3710	IF (type_trac == 'inca') THEN
3711	#ifdef INCA
3712	CALL VTe(VTphysiq)
3713	CALL VTb(VTinca)
3714
3715	CALL chemhook_end ( &
3716	dtime, &
3717	pplay, &
3718	t_seri, &
3719	tr_seri, &
3720	nbtr, &
3721	paprs, &
3722	q_seri, &
3723	airephy, &
3724	pphi, &
3725	pphis, &
3726	zx_rh)
3727
3728	CALL VTe(VTinca)
3729	CALL VTb(VTphysiq)
3730	#endif
3731	END IF
3732
3733
3734	!
3735	! Convertir les incrementations en tendances
3736	!
3737	IF (prt_level .GE.10) THEN
3738	print *,'Convertir les incrementations en tendances '
3739	ENDIF
3740	!
3741	if (mydebug) then
3742	call writefield_phy('u_seri',u_seri,llm)
3743	call writefield_phy('v_seri',v_seri,llm)
3744	call writefield_phy('t_seri',t_seri,llm)
3745	call writefield_phy('q_seri',q_seri,llm)
3746	endif
3747
3748	DO k = 1, klev
3749	DO i = 1, klon
3750	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
3751	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
3752	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
3753	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
3754	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
3755	!CR: on ajoute le contenu en glace
3756	if (nqo.eq.3) then
3757	d_qx(i,k,isol) = ( qs_seri(i,k) - qx(i,k,isol) ) / dtime
3758	endif
3759	ENDDO
3760	ENDDO
3761	!
3762	!CR: nb de traceurs eau: nqo
3763	! IF (nqtot.GE.3) THEN
3764	IF (nqtot.GE.(nqo+1)) THEN
3765	! DO iq = 3, nqtot
3766	DO iq = nqo+1, nqtot
3767	DO k = 1, klev
3768	DO i = 1, klon
3769	! d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
3770	d_qx(i,k,iq) = ( tr_seri(i,k,iq-nqo) - qx(i,k,iq) ) / dtime
3771	ENDDO
3772	ENDDO
3773	ENDDO
3774	ENDIF
3775	!
3776	!IM rajout diagnostiques bilan KP pour analyse MJO par Jun-Ichi Yano
3777	!IM global posePB include "write_bilKP_ins.h"
3778	!IM global posePB include "write_bilKP_ave.h"
3779	!
3780
3781	! Sauvegarder les valeurs de t et q a la fin de la physique:
3782	!
3783	DO k = 1, klev
3784	DO i = 1, klon
3785	u_ancien(i,k) = u_seri(i,k)
3786	v_ancien(i,k) = v_seri(i,k)
3787	t_ancien(i,k) = t_seri(i,k)
3788	q_ancien(i,k) = q_seri(i,k)
3789	ENDDO
3790	ENDDO
3791
3792	!!! RomP >>>
3793	!CR: nb de traceurs eau: nqo
3794	! IF (nqtot.GE.3) THEN
3795	IF (nqtot.GE.(nqo+1)) THEN
3796	! DO iq = 3, nqtot
3797	DO iq = nqo+1, nqtot
3798	DO k = 1, klev
3799	DO i = 1, klon
3800	! tr_ancien(i,k,iq-2) = tr_seri(i,k,iq-2)
3801	tr_ancien(i,k,iq-nqo) = tr_seri(i,k,iq-nqo)
3802	ENDDO
3803	ENDDO
3804	ENDDO
3805	ENDIF
3806	!!! RomP <<<
3807	!==========================================================================
3808	! Sorties des tendances pour un point particulier
3809	! a utiliser en 1D, avec igout=1 ou en 3D sur un point particulier
3810	! pour le debug
3811	! La valeur de igout est attribuee plus haut dans le programme
3812	!==========================================================================
3813
3814	if (prt_level.ge.1) then
3815	write(lunout,*) 'FIN DE PHYSIQ !!!!!!!!!!!!!!!!!!!!'
3816	write(lunout,*) &
3817	'nlon,klev,nqtot,debut,lafin,jD_cur, jH_cur, pdtphys pct tlos'
3818	write(lunout,*) &
3819	nlon,klev,nqtot,debut,lafin, jD_cur, jH_cur ,pdtphys, &
3820	pctsrf(igout,is_ter), pctsrf(igout,is_lic),pctsrf(igout,is_oce), &
3821	pctsrf(igout,is_sic)
3822	write(lunout,*) 'd_t_dyn,d_t_con,d_t_lsc,d_t_ajsb,d_t_ajs,d_t_eva'
3823	do k=1,klev
3824	write(lunout,*) d_t_dyn(igout,k),d_t_con(igout,k), &
3825	d_t_lsc(igout,k),d_t_ajsb(igout,k),d_t_ajs(igout,k), &
3826	d_t_eva(igout,k)
3827	enddo
3828	write(lunout,*) 'cool,heat'
3829	do k=1,klev
3830	write(lunout,*) cool(igout,k),heat(igout,k)
3831	enddo
3832
3833	write(lunout,*) 'd_t_oli,d_t_vdf,d_t_oro,d_t_lif,d_t_ec'
3834	do k=1,klev
3835	write(lunout,*) d_t_oli(igout,k),d_t_vdf(igout,k), &
3836	d_t_oro(igout,k),d_t_lif(igout,k),d_t_ec(igout,k)
3837	enddo
3838
3839	write(lunout,*) 'd_ps ',d_ps(igout)
3840	write(lunout,*) 'd_u, d_v, d_t, d_qx1, d_qx2 '
3841	do k=1,klev
3842	write(lunout,*) d_u(igout,k),d_v(igout,k),d_t(igout,k), &
3843	d_qx(igout,k,1),d_qx(igout,k,2)
3844	enddo
3845	endif
3846
3847	!==========================================================================
3848
3849	!============================================================
3850	! Calcul de la temperature potentielle
3851	!============================================================
3852	DO k = 1, klev
3853	DO i = 1, klon
3854	!JYG/IM theta en debut du pas de temps
3855	!JYG/IM theta(i,k)=t(i,k)(100000./pplay(i,k))*(RD/RCPD)
3856	!JYG/IM theta en fin de pas de temps de physique
3857	theta(i,k)=t_seri(i,k)(100000./pplay(i,k))*(RD/RCPD)
3858	! thetal: 2 lignes suivantes a decommenter si vous avez les fichiers MPL 20130625
3859	! fth_fonctions.F90 et parkind1.F90
3860	! sinon thetal=theta
3861	! thetal(i,k)=fth_thetal(pplay(i,k),t_seri(i,k),q_seri(i,k),
3862	! : ql_seri(i,k))
3863	thetal(i,k)=theta(i,k)
3864	ENDDO
3865	ENDDO
3866	!
3867
3868	! 22.03.04 BEG
3869	!=============================================================
3870	! Ecriture des sorties
3871	!=============================================================
3872	#ifdef CPP_IOIPSL
3873
3874	! Recupere des varibles calcule dans differents modules
3875	! pour ecriture dans histxxx.nc
3876
3877	! Get some variables from module fonte_neige_mod
3878	CALL fonte_neige_get_vars(pctsrf, &
3879	zxfqcalving, zxfqfonte, zxffonte)
3880
3881
3882
3883
3884	!=============================================================
3885	! Separation entre thermiques et non thermiques dans les sorties
3886	! de fisrtilp
3887	!=============================================================
3888
3889	if (iflag_thermals>=1) then
3890	d_t_lscth=0.
3891	d_t_lscst=0.
3892	d_q_lscth=0.
3893	d_q_lscst=0.
3894	do k=1,klev
3895	do i=1,klon
3896	if (ptconvth(i,k)) then
3897	d_t_lscth(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
3898	d_q_lscth(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
3899	else
3900	d_t_lscst(i,k)=d_t_eva(i,k)+d_t_lsc(i,k)
3901	d_q_lscst(i,k)=d_q_eva(i,k)+d_q_lsc(i,k)
3902	endif
3903	enddo
3904	enddo
3905
3906	do i=1,klon
3907	plul_st(i)=prfl(i,lmax_th(i)+1)+psfl(i,lmax_th(i)+1)
3908	plul_th(i)=prfl(i,1)+psfl(i,1)
3909	enddo
3910	endif
3911
3912
3913	!On effectue les sorties:
3914
3915	CALL phys_output_write(itap, pdtphys, paprs, pphis, &
3916	pplay, lmax_th, aerosol_couple, &
3917	ok_ade, ok_aie, ivap, new_aod, ok_sync, &
3918	ptconv, read_climoz, clevSTD, &
3919	ptconvth, d_t, qx, d_qx, zmasse, &
3920	flag_aerosol, flag_aerosol_strat, ok_cdnc)
3921
3922
3923
3924
3925	include "write_histday_seri.h"
3926
3927	include "write_paramLMDZ_phy.h"
3928
3929	#endif
3930
3931	! 22.03.04 END
3932	!
3933	!====================================================================
3934	! Si c'est la fin, il faut conserver l'etat de redemarrage
3935	!====================================================================
3936	!
3937
3938	! -----------------------------------------------------------------
3939	! WSTATS: Saving statistics
3940	! -----------------------------------------------------------------
3941	! ("stats" stores and accumulates 8 key variables in file "stats.nc"
3942	! which can later be used to make the statistic files of the run:
3943	! "stats") only possible in 3D runs !
3944
3945
3946	IF (callstats) THEN
3947
3948	call wstats(klon,o_psol%name,"Surface pressure","Pa" &
3949	,2,paprs(:,1))
3950	call wstats(klon,o_tsol%name,"Surface temperature","K", &
3951	2,zxtsol)
3952	zx_tmp_fi2d(:) = rain_fall(:) + snow_fall(:)
3953	call wstats(klon,o_precip%name,"Precip Totale liq+sol", &
3954	"kg/(s*m2)",2,zx_tmp_fi2d)
3955	zx_tmp_fi2d(:) = rain_lsc(:) + snow_lsc(:)
3956	call wstats(klon,o_plul%name,"Large-scale Precip", &
3957	"kg/(s*m2)",2,zx_tmp_fi2d)
3958	zx_tmp_fi2d(:) = rain_con(:) + snow_con(:)
3959	call wstats(klon,o_pluc%name,"Convective Precip", &
3960	"kg/(s*m2)",2,zx_tmp_fi2d)
3961	call wstats(klon,o_sols%name,"Solar rad. at surf.", &
3962	"W/m2",2,solsw)
3963	call wstats(klon,o_soll%name,"IR rad. at surf.", &
3964	"W/m2",2,sollw)
3965	zx_tmp_fi2d(:) = topsw(:)-toplw(:)
3966	call wstats(klon,o_nettop%name,"Net dn radiatif flux at TOA", &
3967	"W/m2",2,zx_tmp_fi2d)
3968
3969
3970
3971	call wstats(klon,o_temp%name,"Air temperature","K", &
3972	3,t_seri)
3973	call wstats(klon,o_vitu%name,"Zonal wind","m.s-1", &
3974	3,u_seri)
3975	call wstats(klon,o_vitv%name,"Meridional wind", &
3976	"m.s-1",3,v_seri)
3977	call wstats(klon,o_vitw%name,"Vertical wind", &
3978	"m.s-1",3,omega)
3979	call wstats(klon,o_ovap%name,"Specific humidity", "kg/kg", &
3980	3,q_seri)
3981
3982
3983
3984	IF(lafin) THEN
3985	write (,) "Writing stats..."
3986	call mkstats(ierr)
3987	ENDIF
3988
3989	ENDIF !if callstats
3990
3991	IF (lafin) THEN
3992	itau_phy = itau_phy + itap
3993	CALL phyredem ("restartphy.nc")
3994	! open(97,form="unformatted",file="finbin")
3995	! write(97) u_seri,v_seri,t_seri,q_seri
3996	! close(97)
3997	!$OMP MASTER
3998	if (read_climoz >= 1) then
3999	if (is_mpi_root) then
4000	call nf95_close(ncid_climoz)
4001	end if
4002	deallocate(press_climoz) ! pointer
4003	end if
4004	!$OMP END MASTER
4005	ENDIF
4006
4007	! first=.false.
4008
4009	RETURN
4010	END SUBROUTINE physiq
4011	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
4012	IMPLICIT none
4013	!
4014	! Calculer et imprimer l'eau totale. A utiliser pour verifier
4015	! la conservation de l'eau
4016	!
4017	include "YOMCST.h"
4018	INTEGER klon,klev
4019	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
4020	REAL aire(klon)
4021	REAL qtotal, zx, qcheck
4022	INTEGER i, k
4023	!
4024	zx = 0.0
4025	DO i = 1, klon
4026	zx = zx + aire(i)
4027	ENDDO
4028	qtotal = 0.0
4029	DO k = 1, klev
4030	DO i = 1, klon
4031	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i) &
4032	*(paprs(i,k)-paprs(i,k+1))/RG
4033	ENDDO
4034	ENDDO
4035	!
4036	qcheck = qtotal/zx
4037	!
4038	RETURN
4039	END FUNCTION qcheck
4040	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
4041	IMPLICIT none
4042	!
4043	! Tranformer une variable de la grille physique a
4044	! la grille d'ecriture
4045	!
4046	INTEGER nfield,nlon,iim,jjmp1, jjm
4047	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
4048	!
4049	INTEGER i, n, ig
4050	!
4051	jjm = jjmp1 - 1
4052	DO n = 1, nfield
4053	DO i=1,iim
4054	ecrit(i,n) = fi(1,n)
4055	ecrit(i+jjm*iim,n) = fi(nlon,n)
4056	ENDDO
4057	DO ig = 1, nlon - 2
4058	ecrit(iim+ig,n) = fi(1+ig,n)
4059	ENDDO
4060	ENDDO
4061	RETURN
4062	END SUBROUTINE gr_fi_ecrit

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

Download in other formats:

Original Format