Context Navigation

physiq.F90 @ 2086

Last change on this file since 2086 was 2086, checked in by fhourdin, 10 years ago
nclusion de la thermodynamique de la glace Ice thermodynamics (Catherine Rio)
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: 125.6 KB

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

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

Download in other formats:

Original Format