Context Navigation

physiq.F90 @ 2076

Last change on this file since 2076 was 2076, checked in by Laurent Fairhead, 10 years ago

Pour éviter des fsync intempestifs sur le $SCRATCHDIR de curie

To prevent superfluous fsync on curie's disks

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

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

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

Download in other formats:

Original Format