Context Navigation

physiq.F90 @ 2325

Last change on this file since 2325 was 2325, checked in by oboucher, 9 years ago
Enabling aerosols for the case RRTM + NSW=2 In this case we use the old aerosol properties
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: 139.8 KB

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

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

Download in other formats:

Original Format