Context Navigation

physiq.F90 @ 2262

Last change on this file since 2262 was 2262, checked in by jyg, 9 years ago

Correction of a bug concerning the number 'nqo' of
water phases transported by the dynamic : the
default value (= 2, corresponding to vapour and
liquid phases) was still explicitely present in
various places.

Modified files:

infotrac.F90,
physiq.F90

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

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

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

Download in other formats:

Original Format