Context Navigation

physiq.F90 @ 2237

Last change on this file since 2237 was 2236, checked in by fhourdin, 10 years ago

Changement du drapeau iflag_cldth en iflag_cld_th
et du défaut de iflag_cld_cv pour le choix entre schema
de nuages convectifs

lognormal = 0
bigaussien = 1

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

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

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

Download in other formats:

Original Format