Context Navigation

physiq.F90 @ 2271

Last change on this file since 2271 was 2271, checked in by musat, 9 years ago
Corrections for standard level pressure outputs hist*NMC.nc by XIOS and IOIPSL. For XIOS one need to specify a default_value="1.e+20" in the field_def_lmdz.xml and a detect_missing_value=".true." in the file_def_histins_lmdz.xml, file_def_histday_lmdz.xml, etc files.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 138.5 KB

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

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

Download in other formats:

Original Format