Context Navigation

physiq.F @ 467

Last change on this file since 467 was 467, checked in by lmdzadmin, 21 years ago
Modifs sur les seuils (cdrag etc...), inclusion des diagnostics ISCCP par Ionela LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 85.1 KB

Line
1	c
2	c $Header$
3	c
4	SUBROUTINE physiq (nlon,nlev,nqmax,
5	. debut,lafin,rjourvrai,rjour_ecri,gmtime,pdtphys,
6	. paprs,pplay,pphi,pphis,paire,presnivs,clesphy0,
7	. u,v,t,qx,
8	. omega, cufi, cvfi,
9	. d_u, d_v, d_t, d_qx, d_ps)
10	USE ioipsl
11	USE histcom
12	USE writephys
13
14	IMPLICIT none
15	c======================================================================
16	c
17	c Auteur(s) Z.X. Li (LMD/CNRS) date: 19930818
18	c
19	c Objet: Moniteur general de la physique du modele
20	cAA Modifications quant aux traceurs :
21	cAA - uniformisation des parametrisations ds phytrac
22	cAA - stockage des moyennes des champs necessaires
23	cAA en mode traceur off-line
24	c======================================================================
25	c modif ( P. Le Van , 12/10/98 )
26	c
27	c Arguments:
28	c
29	c nlon----input-I-nombre de points horizontaux
30	c nlev----input-I-nombre de couches verticales
31	c nqmax---input-I-nombre de traceurs (y compris vapeur d'eau) = 1
32	c debut---input-L-variable logique indiquant le premier passage
33	c lafin---input-L-variable logique indiquant le dernier passage
34	c rjour---input-R-numero du jour de l'experience
35	c gmtime--input-R-temps universel dans la journee (0 a 86400 s)
36	c pdtphys-input-R-pas d'integration pour la physique (seconde)
37	c paprs---input-R-pression pour chaque inter-couche (en Pa)
38	c pplay---input-R-pression pour le mileu de chaque couche (en Pa)
39	c pphi----input-R-geopotentiel de chaque couche (g z) (reference sol)
40	c pphis---input-R-geopotentiel du sol
41	c paire---input-R-aire de chaque maille
42	c presnivs-input_R_pressions approximat. des milieux couches ( en PA)
43	c u-------input-R-vitesse dans la direction X (de O a E) en m/s
44	c v-------input-R-vitesse Y (de S a N) en m/s
45	c t-------input-R-temperature (K)
46	c qx------input-R-humidite specifique (kg/kg) et d'autres traceurs
47	c d_t_dyn-input-R-tendance dynamique pour "t" (K/s)
48	c d_q_dyn-input-R-tendance dynamique pour "q" (kg/kg/s)
49	c omega---input-R-vitesse verticale en Pa/s
50	c cufi----input-R-resolution des mailles en x (m)
51	c cvfi----input-R-resolution des mailles en y (m)
52	c
53	c d_u-----output-R-tendance physique de "u" (m/s/s)
54	c d_v-----output-R-tendance physique de "v" (m/s/s)
55	c d_t-----output-R-tendance physique de "t" (K/s)
56	c d_qx----output-R-tendance physique de "qx" (kg/kg/s)
57	c d_ps----output-R-tendance physique de la pression au sol
58	c======================================================================
59	#include "dimensions.h"
60	integer jjmp1
61	parameter (jjmp1=jjm+1-1/jjm)
62	#include "dimphy.h"
63	#include "regdim.h"
64	#include "indicesol.h"
65	#include "dimsoil.h"
66	#include "clesphys.h"
67	#include "control.h"
68	#include "temps.h"
69	c======================================================================
70	LOGICAL ok_cvl ! pour activer le nouveau driver pour convection KE
71	PARAMETER (ok_cvl=.TRUE.)
72	LOGICAL ok_gust ! pour activer l'effet des gust sur flux surface
73	PARAMETER (ok_gust=.FALSE.)
74	c======================================================================
75	LOGICAL check ! Verifier la conservation du modele en eau
76	PARAMETER (check=.FALSE.)
77	LOGICAL ok_stratus ! Ajouter artificiellement les stratus
78	PARAMETER (ok_stratus=.FALSE.)
79	c======================================================================
80	c Parametres lies au coupleur OASIS:
81	#include "oasis.h"
82	INTEGER,SAVE :: npas, nexca
83	logical rnpb
84	parameter(rnpb=.true.)
85	c PARAMETER (npas=1440)
86	c PARAMETER (nexca=48)
87	EXTERNAL fromcpl, intocpl, inicma
88	c ocean = type de modele ocean a utiliser: force, slab, couple
89	character*6 ocean
90	SAVE ocean
91
92	c parameter (ocean = 'force ')
93	c parameter (ocean = 'couple')
94	logical ok_ocean
95	c======================================================================
96	c Clef controlant l'activation du cycle diurne:
97	ccc LOGICAL cycle_diurne
98	ccc PARAMETER (cycle_diurne=.FALSE.)
99	c======================================================================
100	c Modele thermique du sol, a activer pour le cycle diurne:
101	ccc LOGICAL soil_model
102	ccc PARAMETER (soil_model=.FALSE.)
103	logical ok_veget
104	save ok_veget
105	c parameter (ok_veget = .true.)
106	c parameter (ok_veget = .false.)
107	c======================================================================
108	c Dans les versions precedentes, l'eau liquide nuageuse utilisee dans
109	c le calcul du rayonnement est celle apres la precipitation des nuages.
110	c Si cette cle new_oliq est activee, ce sera une valeur moyenne entre
111	c la condensation et la precipitation. Cette cle augmente les impacts
112	c radiatifs des nuages.
113	ccc LOGICAL new_oliq
114	ccc PARAMETER (new_oliq=.FALSE.)
115	c======================================================================
116	c Clefs controlant deux parametrisations de l'orographie:
117	cc LOGICAL ok_orodr
118	ccc PARAMETER (ok_orodr=.FALSE.)
119	ccc LOGICAL ok_orolf
120	ccc PARAMETER (ok_orolf=.FALSE.)
121	c======================================================================
122	LOGICAL ok_journe ! sortir le fichier journalier
123	save ok_journe
124	c PARAMETER (ok_journe=.true.)
125
126	integer lev_histday
127	save lev_histday
128	data lev_histday/1/
129	c
130	LOGICAL ok_mensuel ! sortir le fichier mensuel
131	save ok_mensuel
132	c PARAMETER (ok_mensuel=.true.)
133	c
134	LOGICAL ok_instan ! sortir le fichier instantane
135	save ok_instan
136	c PARAMETER (ok_instan=.true.)
137	c
138	LOGICAL ok_region ! sortir le fichier regional
139	PARAMETER (ok_region=.FALSE.)
140	c======================================================================
141	c
142	INTEGER ivap ! indice de traceurs pour vapeur d'eau
143	PARAMETER (ivap=1)
144	INTEGER iliq ! indice de traceurs pour eau liquide
145	PARAMETER (iliq=2)
146
147	c
148	c
149	c Variables argument:
150	c
151	INTEGER nlon
152	INTEGER nlev
153	INTEGER nqmax
154	REAL rjourvrai, rjour_ecri
155	REAL gmtime
156	REAL pdtphys
157	LOGICAL debut, lafin
158	REAL paprs(klon,klev+1)
159	REAL pplay(klon,klev)
160	REAL pphi(klon,klev)
161	REAL pphis(klon)
162	REAL paire(klon)
163	REAL presnivs(klev)
164	REAL znivsig(klev)
165	REAL zsurf(nbsrf)
166	real cufi(klon), cvfi(klon)
167
168	REAL u(klon,klev)
169	REAL v(klon,klev)
170	REAL t(klon,klev)
171	REAL qx(klon,klev,nqmax)
172
173	REAL t_ancien(klon,klev), q_ancien(klon,klev)
174	SAVE t_ancien, q_ancien
175	LOGICAL ancien_ok
176	SAVE ancien_ok
177
178	REAL d_t_dyn(klon,klev)
179	REAL d_q_dyn(klon,klev)
180
181	REAL omega(klon,klev)
182
183	REAL d_u(klon,klev)
184	REAL d_v(klon,klev)
185	REAL d_t(klon,klev)
186	REAL d_qx(klon,klev,nqmax)
187	REAL d_ps(klon)
188
189	cccIM
190	INTEGER klevp1
191	PARAMETER(klevp1=klev+1)
192	#include "raddim.h"
193	cc REAL*8 ZFSUP(KDLON,KFLEV+1)
194	cc REAL*8 ZFSDN(KDLON,KFLEV+1)
195	cc REAL*8 ZFSUP0(KDLON,KFLEV+1)
196	cc REAL*8 ZFSDN0(KDLON,KFLEV+1)
197	c
198	REAL swdn0(klon,2), swdn(klon,2), swup0(klon,2), swup(klon,2)
199	SAVE swdn0 , swdn, swup0, swup
200
201	cccIM cf. FH
202	real u850(klon),v850(klon),u200(klon),v200(klon)
203	real u500(klon),v500(klon),phi500(klon),w500(klon)
204	cIM
205	real prw(klon)
206
207	cIM ISCCP - proprietes microphysiques des nuages convectifs
208	REAL convliq(klon,klev) ! eau liquide nuageuse convective
209	REAL convfra(klon,klev) ! fraction nuageuse convective
210
211	REAL cldl_c(klon),cldm_c(klon),cldh_c(klon) !nuages bas, moyen et haut
212	REAL cldt_c(klon),cldq_c(klon) !nuage total, eau liquide integree
213	REAL cldl_s(klon),cldm_s(klon),cldh_s(klon) !nuages bas, moyen et haut
214	REAL cldt_s(klon),cldq_s(klon) !nuage total, eau liquide integree
215
216	INTEGER kinv, linv
217
218	cIM ISCCP simulator BEGIN
219	INTEGER igfi2D(iim,jjmp1)
220	cv3.4
221	INTEGER debug, debugcol
222	INTEGER npoints
223	PARAMETER(npoints=klon)
224	INTEGER sunlit(klon)
225
226	INTEGER ncol, seed(klon)
227
228	cIM dans clesphys.h top_height, overlap
229	c PARAMETER(ncol=100)
230	c PARAMETER(ncol=625)
231	PARAMETER(ncol=10)
232	REAL tautab(0:255)
233	INTEGER invtau(-20:45000)
234	REAL emsfc_lw
235	PARAMETER(emsfc_lw=0.99)
236	REAL ran0 ! type for random number fuction
237
238	REAL pfull(klon,klev)
239	REAL phalf(klon,klev+1)
240	REAL cldtot(klon,klev)
241	REAL dtau_s(klon,klev)
242	REAL dtau_c(klon,klev)
243	REAL dem_s(klon,klev)
244	REAL dem_c(klon,klev)
245	cPLUS : variables de haut en bas pour le simulateur ISCCP
246	REAL qv(klon,klev)
247	REAL cc(klon,klev)
248	REAL conv(klon,klev)
249	REAL dtau_sH2B(klon,klev)
250	REAL dtau_cH2B(klon,klev)
251	REAL at(klon,klev)
252	REAL dem_sH2B(klon,klev)
253	REAL dem_cH2B(klon,klev)
254
255	c output from ISCCP
256	REAL fq_isccp(klon,7,7)
257	REAL totalcldarea(klon)
258	REAL meanptop(klon)
259	REAL meantaucld(klon)
260	REAL boxtau(klon,ncol)
261	REAL boxptop(klon,ncol)
262
263	c grille 4d physique
264	INTEGER l, ni, nj, kmax, lmax, nrec
265	INTEGER ni1, ni2, nj1, nj2
266	c PARAMETER(kmax=7, lmax=7)
267	PARAMETER(kmax=8, lmax=8)
268	INTEGER kmaxm1, lmaxm1
269	PARAMETER(kmaxm1=kmax-1, lmaxm1=lmax-1)
270	c INTEGER iimx7, jjmx7, jjmp1x7
271	c PARAMETER(iimx7=iim7, jjmx7=jjm7, jjmp1x7=jjmp1*7)
272	c REAL fq4d(iim,jjmp1,7,7)
273	c REAL fq3d(iimx7, jjmp1x7)
274	INTEGER iimx8, jjmx8, jjmp1x8
275	PARAMETER(iimx8=iim8, jjmx8=jjm8, jjmp1x8=jjmp1*8)
276	REAL fq4d(iim,jjmp1,8,8)
277	REAL fq3d(iimx8, jjmp1x8)
278	cIM180603 SAVE fq3d
279
280	c REAL maxfq3d, minfq3d
281	c
282	INTEGER iw, iwmax
283	REAL wmin, pas_w
284	c PARAMETER(wmin=-100.,pas_w=10.,iwmax=30)
285	PARAMETER(wmin=-200.,pas_w=10.,iwmax=40)
286	REAL o500(klon)
287	INTEGER nreg, nbreg
288	PARAMETER(nbreg=5)
289	c REAL histoW(iwmax,kmaxm1,lmaxm1)
290	REAL histoW(kmaxm1,lmaxm1,iwmax,nbreg)
291	REAL nhistoW(kmaxm1,lmaxm1,iwmax,nbreg)
292	cIM180603
293	c SAVE histoW, nhistoW
294	c SAVE nhistoW
295	REAL nhistoWt(kmaxm1,lmaxm1,iwmax,nbreg)
296	SAVE nhistoWt
297
298	c REAL histoWinv(kmaxm1,lmaxm1,iwmax)
299	c REAL nhistoW(kmaxm1,lmaxm1,iwmax)
300	INTEGER linv
301	c LOGICAL pct_ocean(klon,nbreg)
302	INTEGER pct_ocean(klon,nbreg)
303	REAL rlonPOS(klon)
304	c CHARACTER*4 pdirect
305
306	c sorties ISCCP
307
308	logical ok_isccp
309	real ecrit_isccp
310	integer nid_isccp
311	save ok_isccp, ecrit_isccp, nid_isccp
312
313	#define histISCCP
314	#undef histISCCP
315	#ifdef histISCCP
316	c data ok_isccp,ecrit_isccp/.true.,0.125/
317	c data ok_isccp,ecrit_isccp/.true.,1./
318	data ok_isccp/.true./
319	#else
320	data ok_isccp/.false./
321	#endif
322
323	REAL zx_tau(kmaxm1), zx_pc(lmaxm1), zx_o500(iwmax)
324	c DATA zx_tau/0.1, 1.3, 3.6, 9.4, 23., 60./
325	c DATA zx_pc/50., 180., 310., 440., 560., 680., 800., 1015./
326	c DATA zx_pc/50., 180., 310., 440., 560., 680., 800./
327	cOK DATA zx_tau/0.0, 0.1, 1.3, 3.6, 9.4, 23., 60./
328	cOK DATA zx_pc/800., 680., 560., 440., 310., 180., 50./
329
330	c tester l'alure
331	DATA zx_tau/1., 2., 3., 4., 5., 6., 7./
332	c DATA zx_pc/1., 2., 3., 4., 5., 6., 7./
333	DATA zx_pc/7., 6., 5., 4., 3., 2., 1./
334
335	INTEGER komega, nhoriRD
336
337	c statistiques regime dynamique END
338
339	c REAL del_lon(iim), del_lat(jjmp1)
340	REAL del_lon, del_lat
341	c REAL zx_lonx7(iimx7), zx_latx7(jjmp1x7)
342	REAL zx_lonx8(iimx8), zx_latx8(jjmp1x8)
343	c INTEGER nhorix7
344	INTEGER nhorix8
345
346	cIM ISCCP simulator END
347
348	logical ok_hf
349	real ecrit_hf
350	integer nid_hf
351	save ok_hf, ecrit_hf, nid_hf
352
353	c QUESTION : noms de variables ?
354
355	#define histhf
356	#ifdef histhf
357	data ok_hf,ecrit_hf/.true.,0.25/
358	#else
359	data ok_hf/.false./
360	#endif
361
362	INTEGER longcles
363	PARAMETER ( longcles = 20 )
364	REAL clesphy0( longcles )
365	c
366	c Variables quasi-arguments
367	c
368	REAL xjour
369	SAVE xjour
370	c
371	c
372	c Variables propres a la physique
373	c
374	REAL dtime
375	SAVE dtime ! pas temporel de la physique
376	c
377	INTEGER radpas
378	SAVE radpas ! frequence d'appel rayonnement
379	c
380	REAL radsol(klon)
381	SAVE radsol ! bilan radiatif au sol calcule par code radiatif
382	c
383	REAL rlat(klon)
384	SAVE rlat ! latitude pour chaque point
385	c
386	REAL rlon(klon)
387	SAVE rlon ! longitude pour chaque point
388	c
389	cc INTEGER iflag_con
390	cc SAVE iflag_con ! indicateur de la convection
391	c
392	INTEGER itap
393	SAVE itap ! compteur pour la physique
394	c
395	REAL co2_ppm_etat0
396	c
397	REAL solaire_etat0
398	c
399	real slp(klon) ! sea level pressure
400
401	REAL ftsol(klon,nbsrf)
402	SAVE ftsol ! temperature du sol
403	c
404	REAL ftsoil(klon,nsoilmx,nbsrf)
405	SAVE ftsoil ! temperature dans le sol
406	c
407	REAL fevap(klon,nbsrf)
408	SAVE fevap ! evaporation
409	REAL fluxlat(klon,nbsrf)
410	SAVE fluxlat
411	c
412	REAL deltat(klon)
413	SAVE deltat ! ecart avec la SST de reference
414	c
415	REAL fqsurf(klon,nbsrf)
416	SAVE fqsurf ! humidite de l'air au contact de la surface
417	c
418	REAL qsol(klon)
419	SAVE qsol ! hauteur d'eau dans le sol
420	c
421	REAL fsnow(klon,nbsrf)
422	SAVE fsnow ! epaisseur neigeuse
423	c
424	REAL falbe(klon,nbsrf)
425	SAVE falbe ! albedo par type de surface
426	REAL falblw(klon,nbsrf)
427	SAVE falblw ! albedo par type de surface
428
429	c
430	c
431	c Parametres de l'Orographie a l'Echelle Sous-Maille (OESM):
432	c
433	REAL zmea(klon)
434	SAVE zmea ! orographie moyenne
435	c
436	REAL zstd(klon)
437	SAVE zstd ! deviation standard de l'OESM
438	c
439	REAL zsig(klon)
440	SAVE zsig ! pente de l'OESM
441	c
442	REAL zgam(klon)
443	save zgam ! anisotropie de l'OESM
444	c
445	REAL zthe(klon)
446	SAVE zthe ! orientation de l'OESM
447	c
448	REAL zpic(klon)
449	SAVE zpic ! Maximum de l'OESM
450	c
451	REAL zval(klon)
452	SAVE zval ! Minimum de l'OESM
453	c
454	REAL rugoro(klon)
455	SAVE rugoro ! longueur de rugosite de l'OESM
456	c
457	REAL zulow(klon),zvlow(klon),zustr(klon), zvstr(klon)
458	c
459	REAL zuthe(klon),zvthe(klon)
460	SAVE zuthe
461	SAVE zvthe
462	INTEGER igwd,idx(klon),itest(klon)
463	c
464	REAL agesno(klon,nbsrf)
465	SAVE agesno ! age de la neige
466	c
467	REAL alb_neig(klon)
468	SAVE alb_neig ! albedo de la neige
469	cKE43
470	c Variables liees a la convection de K. Emanuel (sb):
471	c
472	REAL ema_workcbmf(klon) ! cloud base mass flux
473	SAVE ema_workcbmf
474
475	REAL ema_cbmf(klon) ! cloud base mass flux
476	SAVE ema_cbmf
477
478	REAL ema_pcb(klon) ! cloud base pressure
479	SAVE ema_pcb
480
481	REAL ema_pct(klon) ! cloud top pressure
482	SAVE ema_pct
483
484	REAL bas, top ! cloud base and top levels
485	SAVE bas
486	SAVE top
487
488	REAL Ma(klon,klev) ! undilute upward mass flux
489	SAVE Ma
490	REAL qcondc(klon,klev) ! in-cld water content from convect
491	SAVE qcondc
492	REAL ema_work1(klon, klev), ema_work2(klon, klev)
493	SAVE ema_work1, ema_work2
494	REAL wdn(klon), tdn(klon), qdn(klon)
495
496	REAL wd(klon) ! sb
497	SAVE wd ! sb
498
499	c Variables locales pour la couche limite (al1):
500	c
501	cAl1 REAL pblh(klon) ! Hauteur de couche limite
502	cAl1 SAVE pblh
503	c34EK
504	c
505	c Variables locales:
506	c
507	REAL cdragh(klon) ! drag coefficient pour T and Q
508	REAL cdragm(klon) ! drag coefficient pour vent
509	cAA
510	cAA Pour phytrac
511	cAA
512	REAL ycoefh(klon,klev) ! coef d'echange pour phytrac
513	REAL yu1(klon) ! vents dans la premiere couche U
514	REAL yv1(klon) ! vents dans la premiere couche V
515	cIM cf JLD
516	REAL ffonte(klon,nbsrf) !Flux thermique utilise pour fondre la neige
517	REAL fqcalving(klon,nbsrf) !Flux d'eau "perdue" par la surface
518	c !et necessaire pour limiter la
519	c !hauteur de neige, en kg/m2/s
520	REAL zxffonte(klon), zxfqcalving(klon)
521
522	LOGICAL offline ! Controle du stockage ds "physique"
523	PARAMETER (offline=.false.)
524	INTEGER physid
525	REAL pfrac_impa(klon,klev)! Produits des coefs lessivage impaction
526	save pfrac_impa
527	REAL pfrac_nucl(klon,klev)! Produits des coefs lessivage nucleation
528	save pfrac_nucl
529	REAL pfrac_1nucl(klon,klev)! Produits des coefs lessi nucl (alpha = 1)
530	save pfrac_1nucl
531	REAL frac_impa(klon,klev) ! fractions d'aerosols lessivees (impaction)
532	REAL frac_nucl(klon,klev) ! idem (nucleation)
533	cAA
534	REAL rain_fall(klon) ! pluie
535	REAL snow_fall(klon) ! neige
536	save snow_fall, rain_fall
537	REAL evap(klon), devap(klon) ! evaporation et sa derivee
538	REAL sens(klon), dsens(klon) ! chaleur sensible et sa derivee
539	REAL dlw(klon) ! derivee infra rouge
540	REAL bils(klon) ! bilan de chaleur au sol
541	cIM cf. JLD
542	REAL wfbils(klon,nbsrf) ! bilan de chaleur au sol, pour chaque
543	C type de sous-surface et pondere par la fraction
544	REAL fder(klon) ! Derive de flux (sensible et latente)
545	save fder
546	REAL ve(klon) ! integr. verticale du transport meri. de l'energie
547	REAL vq(klon) ! integr. verticale du transport meri. de l'eau
548	REAL ue(klon) ! integr. verticale du transport zonal de l'energie
549	REAL uq(klon) ! integr. verticale du transport zonal de l'eau
550	c
551	REAL frugs(klon,nbsrf) ! longueur de rugosite
552	save frugs
553	REAL zxrugs(klon) ! longueur de rugosite
554	c
555	c Conditions aux limites
556	c
557	INTEGER julien
558	c
559	INTEGER lmt_pas
560	SAVE lmt_pas ! frequence de mise a jour
561	REAL pctsrf(klon,nbsrf)
562	SAVE pctsrf ! sous-fraction du sol
563	REAL albsol(klon)
564	SAVE albsol ! albedo du sol total
565	REAL albsollw(klon)
566	SAVE albsollw ! albedo du sol total
567
568	REAL wo(klon,klev)
569	SAVE wo ! ozone
570	c======================================================================
571	c
572	c Declaration des procedures appelees
573	c
574	EXTERNAL angle ! calculer angle zenithal du soleil
575	EXTERNAL alboc ! calculer l'albedo sur ocean
576	EXTERNAL albsno ! calculer albedo sur neige
577	EXTERNAL ajsec ! ajustement sec
578	EXTERNAL clmain ! couche limite
579	EXTERNAL condsurf ! lire les conditions aux limites
580	EXTERNAL conlmd ! convection (schema LMD)
581	cKE43
582	EXTERNAL conema3 ! convect4.3
583	EXTERNAL fisrtilp ! schema de condensation a grande echelle (pluie)
584	cAA
585	EXTERNAL fisrtilp_tr ! schema de condensation a grande echelle (pluie)
586	c ! stockage des coefficients necessaires au
587	c ! lessivage OFF-LINE et ON-LINE
588	cAA
589	EXTERNAL hgardfou ! verifier les temperatures
590	EXTERNAL nuage ! calculer les proprietes radiatives
591	EXTERNAL o3cm ! initialiser l'ozone
592	EXTERNAL orbite ! calculer l'orbite terrestre
593	EXTERNAL ozonecm ! prescrire l'ozone
594	EXTERNAL phyetat0 ! lire l'etat initial de la physique
595	EXTERNAL phyredem ! ecrire l'etat de redemarrage de la physique
596	EXTERNAL radlwsw ! rayonnements solaire et infrarouge
597	EXTERNAL suphec ! initialiser certaines constantes
598	EXTERNAL transp ! transport total de l'eau et de l'energie
599	EXTERNAL ecribina ! ecrire le fichier binaire global
600	EXTERNAL ecribins ! ecrire le fichier binaire global
601	EXTERNAL ecrirega ! ecrire le fichier binaire regional
602	EXTERNAL ecriregs ! ecrire le fichier binaire regional
603	c
604	c Variables locales
605	c
606	real clwcon(klon,klev),rnebcon(klon,klev)
607	real clwcon0(klon,klev),rnebcon0(klon,klev)
608	save rnebcon, clwcon
609
610	REAL rhcl(klon,klev) ! humiditi relative ciel clair
611	REAL dialiq(klon,klev) ! eau liquide nuageuse
612	REAL diafra(klon,klev) ! fraction nuageuse
613	REAL cldliq(klon,klev) ! eau liquide nuageuse
614	REAL cldfra(klon,klev) ! fraction nuageuse
615	REAL cldtau(klon,klev) ! epaisseur optique
616	REAL cldemi(klon,klev) ! emissivite infrarouge
617	c
618	CXXX PB
619	REAL fluxq(klon,klev, nbsrf) ! flux turbulent d'humidite
620	REAL fluxt(klon,klev, nbsrf) ! flux turbulent de chaleur
621	REAL fluxu(klon,klev, nbsrf) ! flux turbulent de vitesse u
622	REAL fluxv(klon,klev, nbsrf) ! flux turbulent de vitesse v
623	c
624	REAL zxfluxt(klon, klev)
625	REAL zxfluxq(klon, klev)
626	REAL zxfluxu(klon, klev)
627	REAL zxfluxv(klon, klev)
628	CXXX
629	REAL heat(klon,klev) ! chauffage solaire
630	REAL heat0(klon,klev) ! chauffage solaire ciel clair
631	REAL cool(klon,klev) ! refroidissement infrarouge
632	REAL cool0(klon,klev) ! refroidissement infrarouge ciel clair
633	REAL topsw(klon), toplw(klon), solsw(klon), sollw(klon)
634	real sollwdown(klon) ! downward LW flux at surface
635	REAL topsw0(klon), toplw0(klon), solsw0(klon), sollw0(klon)
636	REAL albpla(klon)
637	cIM cf. JLD
638	REAL fsollw(klon, nbsrf) ! bilan flux IR pour chaque sous surface
639	REAL fsolsw(klon, nbsrf) ! flux solaire absorb. pour chaque sous surface
640	c Le rayonnement n'est pas calcule tous les pas, il faut donc
641	c sauvegarder les sorties du rayonnement
642	SAVE heat,cool,albpla,topsw,toplw,solsw,sollw,sollwdown
643	SAVE topsw0,toplw0,solsw0,sollw0, heat0, cool0
644	cccIM
645
646	INTEGER itaprad
647	SAVE itaprad
648	c
649	REAL conv_q(klon,klev) ! convergence de l'humidite (kg/kg/s)
650	REAL conv_t(klon,klev) ! convergence de la temperature(K/s)
651	c
652	REAL cldl(klon),cldm(klon),cldh(klon) !nuages bas, moyen et haut
653	REAL cldt(klon),cldq(klon) !nuage total, eau liquide integree
654	c
655	REAL zxtsol(klon), zxqsurf(klon), zxsnow(klon), zxfluxlat(klon)
656	c
657	REAL dist, rmu0(klon), fract(klon)
658	REAL zdtime, zlongi
659	c
660	CHARACTER*2 str2
661	CHARACTER*2 iqn
662	c
663	REAL qcheck
664	REAL z_avant(klon), z_apres(klon), z_factor(klon)
665	LOGICAL zx_ajustq
666	c
667	REAL za, zb
668	REAL zx_t, zx_qs, zdelta, zcor, zfra, zlvdcp, zlsdcp
669	real zqsat(klon,klev)
670	INTEGER i, k, iq, ig, j, nsrf, ll
671	REAL t_coup
672	PARAMETER (t_coup=234.0)
673	c
674	REAL zphi(klon,klev)
675	REAL zx_tmp_x(iim), zx_tmp_yjjmp1
676	REAL zx_relief(iim,jjmp1)
677	REAL zx_aire(iim,jjmp1)
678	cKE43
679	c Variables locales pour la convection de K. Emanuel (sb):
680	c
681	REAL upwd(klon,klev) ! saturated updraft mass flux
682	REAL dnwd(klon,klev) ! saturated downdraft mass flux
683	REAL dnwd0(klon,klev) ! unsaturated downdraft mass flux
684	REAL tvp(klon,klev) ! virtual temp of lifted parcel
685	REAL cape(klon) ! CAPE
686	SAVE cape
687	cccIM
688	CHARACTER*40 capemaxcels
689
690	REAL pbase(klon) ! cloud base pressure
691	SAVE pbase
692	REAL bbase(klon) ! cloud base buoyancy
693	SAVE bbase
694	REAL rflag(klon) ! flag fonctionnement de convect
695	INTEGER iflagctrl(klon) ! flag fonctionnement de convect
696	c -- convect43:
697	INTEGER ntra ! nb traceurs pour convect4.3
698	REAL pori_con(klon) ! pressure at the origin level of lifted parcel
699	REAL plcl_con(klon),dtma_con(klon),dtlcl_con(klon)
700	REAL dtvpdt1(klon,klev), dtvpdq1(klon,klev)
701	REAL dplcldt(klon), dplcldr(klon)
702	c? . condm_con(klon,klev),conda_con(klon,klev),
703	c? . mr_con(klon,klev),ep_con(klon,klev)
704	c? . ,sadiab(klon,klev),wadiab(klon,klev)
705	c --
706	c34EK
707	c
708	c Variables du changement
709	c
710	c con: convection
711	c lsc: condensation a grande echelle (Large-Scale-Condensation)
712	c ajs: ajustement sec
713	c eva: evaporation de l'eau liquide nuageuse
714	c vdf: couche limite (Vertical DiFfusion)
715	REAL d_t_con(klon,klev),d_q_con(klon,klev)
716	REAL d_u_con(klon,klev),d_v_con(klon,klev)
717	REAL d_t_lsc(klon,klev),d_q_lsc(klon,klev),d_ql_lsc(klon,klev)
718	REAL d_t_ajs(klon,klev), d_q_ajs(klon,klev)
719	REAL d_t_eva(klon,klev),d_q_eva(klon,klev)
720	REAL rneb(klon,klev)
721	c
722	REAL pmfu(klon,klev), pmfd(klon,klev)
723	REAL pen_u(klon,klev), pen_d(klon,klev)
724	REAL pde_u(klon,klev), pde_d(klon,klev)
725	INTEGER kcbot(klon), kctop(klon), kdtop(klon)
726	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1)
727	REAL prfl(klon,klev+1), psfl(klon,klev+1)
728	c
729	INTEGER ibas_con(klon), itop_con(klon)
730	REAL rain_con(klon), rain_lsc(klon)
731	REAL snow_con(klon), snow_lsc(klon)
732	REAL d_ts(klon,nbsrf)
733	c
734	REAL d_u_vdf(klon,klev), d_v_vdf(klon,klev)
735	REAL d_t_vdf(klon,klev), d_q_vdf(klon,klev)
736	c
737	REAL d_u_oro(klon,klev), d_v_oro(klon,klev)
738	REAL d_t_oro(klon,klev)
739	REAL d_u_lif(klon,klev), d_v_lif(klon,klev)
740	REAL d_t_lif(klon,klev)
741
742	REAL ratqs(klon,klev),ratqss(klon,klev),ratqsc(klon,klev)
743	real ratqsbas,ratqshaut
744	save ratqsbas,ratqshaut, ratqs
745	real zpt_conv(klon,klev)
746
747	c Parametres lies au nouveau schema de nuages (SB, PDF)
748	real fact_cldcon
749	real facttemps
750	logical ok_newmicro
751	save ok_newmicro
752	save fact_cldcon,facttemps
753	real facteur
754
755	integer iflag_cldcon
756	save iflag_cldcon
757
758	logical ptconv(klon,klev)
759
760	c
761	c Variables liees a l'ecriture de la bande histoire physique
762	c
763	INTEGER ecrit_mth
764	SAVE ecrit_mth ! frequence d'ecriture (fichier mensuel)
765	c
766	INTEGER ecrit_day
767	SAVE ecrit_day ! frequence d'ecriture (fichier journalier)
768	c
769	INTEGER ecrit_ins
770	SAVE ecrit_ins ! frequence d'ecriture (fichier instantane)
771	c
772	INTEGER ecrit_reg
773	SAVE ecrit_reg ! frequence d'ecriture
774	c
775	integer itau_w ! pas de temps ecriture = itap + itau_phy
776	c
777	c
778	c Variables locales pour effectuer les appels en serie
779	c
780	REAL t_seri(klon,klev), q_seri(klon,klev)
781	REAL ql_seri(klon,klev),qs_seri(klon,klev)
782	REAL u_seri(klon,klev), v_seri(klon,klev)
783	c
784	REAL tr_seri(klon,klev,nbtr)
785	REAL d_tr(klon,klev,nbtr)
786
787	REAL zx_rh(klon,klev)
788
789	INTEGER length
790	PARAMETER ( length = 100 )
791	REAL tabcntr0( length )
792	c
793	INTEGER ndex2d(iimjjmp1),ndex3d(iimjjmp1*klev)
794	REAL zx_tmp_fi2d(klon)
795	REAL zx_tmp_2d(iim,jjmp1), zx_tmp_3d(iim,jjmp1,klev)
796	REAL zx_lon(iim,jjmp1), zx_lat(iim,jjmp1)
797	c
798	INTEGER nid_day, nid_mth, nid_ins
799	SAVE nid_day, nid_mth, nid_ins
800	c
801	INTEGER nhori, nvert
802	REAL zsto, zout
803	real zjulian
804	save zjulian
805
806	character*20 modname
807	character*80 abort_message
808	logical ok_sync
809	real date0
810	integer idayref
811
812	C essai writephys
813	integer fid_day, fid_mth, fid_ins
814	parameter (fid_ins = 1, fid_day = 2, fid_mth = 3)
815	integer prof2d_on, prof3d_on, prof2d_av, prof3d_av
816	parameter (prof2d_on = 1, prof3d_on = 2,
817	. prof2d_av = 3, prof3d_av = 4)
818	character*30 nom_fichier
819	character*10 varname
820	character*40 vartitle
821	character*20 varunits
822	C Variables liees au bilan d'energie et d'enthalpi
823	REAL ztsol(klon)
824	REAL h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
825	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
826	SAVE h_vcol_tot, h_dair_tot, h_qw_tot, h_ql_tot
827	$ , h_qs_tot, qw_tot, ql_tot, qs_tot , ec_tot
828	REAL d_h_vcol, d_h_dair, d_qt, d_qw, d_ql, d_qs, d_ec
829	REAL d_h_vcol_phy
830	REAL fs_bound, fq_bound
831	SAVE d_h_vcol_phy
832	REAL zero_v(klon)
833	CHARACTER*15 ztit
834	INTEGER ip_ebil ! PRINT level for energy conserv. diag.
835	SAVE ip_ebil
836	DATA ip_ebil/2/
837	INTEGER if_ebil ! level for energy conserv. dignostics
838	SAVE if_ebil
839	c+jld ec_conser
840	REAL d_t_ec(klon,klev) ! tendance du a la conersion Ec -> E thermique
841	REAL ZRCPD
842	c-jld ec_conser
843	cIM
844	REAL t2m(klon,nbsrf), q2m(klon,nbsrf)
845	REAL u10m(klon,nbsrf), v10m(klon,nbsrf)
846	REAL zt2m(klon), zq2m(klon)
847	REAL zu10m(klon), zv10m(klon)
848	CHARACTER*40 t2mincels, t2maxcels
849	c
850	c Declaration des constantes et des fonctions thermodynamiques
851	c
852	#include "YOMCST.h"
853	#include "YOETHF.h"
854	#include "FCTTRE.h"
855	c======================================================================
856	modname = 'physiq'
857	IF (if_ebil.ge.1) THEN
858	DO i=1,klon
859	zero_v(i)=0.
860	END DO
861	END IF
862	ok_sync=.TRUE.
863	IF (nqmax .LT. 2) THEN
864	PRINT*, 'eaux vapeur et liquide sont indispensables'
865	CALL ABORT
866	ENDIF
867	IF (debut) THEN
868	CALL suphec ! initialiser constantes et parametres phys.
869	ENDIF
870
871
872	c======================================================================
873	xjour = rjourvrai
874	c
875	c Si c'est le debut, il faut initialiser plusieurs choses
876	c ********
877	c
878	IF (debut) THEN
879	C
880	IF (if_ebil.ge.1) d_h_vcol_phy=0.
881	c
882	c appel a la lecture du run.def physique
883	c
884	call conf_phys(ocean, ok_veget, ok_journe, ok_mensuel,
885	. ok_instan, fact_cldcon, facttemps,ok_newmicro,
886	. iflag_cldcon,ratqsbas,ratqshaut, if_ebil)
887	cIM . , RI0)
888
889	c
890	c
891	c Initialiser les compteurs:
892	c
893
894	frugs = 0.
895	itap = 0
896	itaprad = 0
897	c
898	CALL phyetat0 ("startphy.nc",dtime,co2_ppm_etat0,solaire_etat0,
899	. rlat,rlon,pctsrf, ftsol,ftsoil,deltat,fqsurf,qsol,fsnow,
900	. falbe, falblw, fevap, rain_fall,snow_fall,solsw, sollwdown,
901	. dlw,radsol,frugs,agesno,clesphy0,
902	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,tabcntr0,
903	. t_ancien, q_ancien, ancien_ok, rnebcon, ratqs,clwcon )
904
905	c
906	radpas = NINT( 86400./dtime/nbapp_rad)
907
908	c
909	CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
910	, ok_instan, ok_region )
911	c
912	IF (ABS(dtime-pdtphys).GT.0.001) THEN
913	PRINT*, 'Pas physique n est pas correcte',dtime,pdtphys
914	abort_message=' See above '
915	call abort_gcm(modname,abort_message,1)
916	ENDIF
917	IF (nlon .NE. klon) THEN
918	PRINT*, 'nlon et klon ne sont pas coherents', nlon, klon
919	abort_message=' See above '
920	call abort_gcm(modname,abort_message,1)
921	ENDIF
922	IF (nlev .NE. klev) THEN
923	PRINT*, 'nlev et klev ne sont pas coherents', nlev, klev
924	abort_message=' See above '
925	call abort_gcm(modname,abort_message,1)
926	ENDIF
927	c
928	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
929	PRINT*, 'Nbre d appels au rayonnement insuffisant'
930	PRINT*, "Au minimum 4 appels par jour si cycle diurne"
931	abort_message=' See above '
932	call abort_gcm(modname,abort_message,1)
933	ENDIF
934	PRINT*, "Clef pour la convection, iflag_con=", iflag_con
935	PRINT*, "Clef pour le driver de la convection, ok_cvl=", ok_cvl
936	c
937	cKE43
938	c Initialisation pour la convection de K.E. (sb):
939	IF (iflag_con.GE.3) THEN
940
941	PRINT, "** Convection de Kerry Emanuel 4.3 "
942	PRINT*, "On va utiliser le melange convectif des traceurs qui"
943	PRINT*, "est calcule dans convect4.3"
944	PRINT*, " !!! penser aux logical flags de phytrac"
945
946	DO i = 1, klon
947	ema_cbmf(i) = 0.
948	ema_pcb(i) = 0.
949	ema_pct(i) = 0.
950	ema_workcbmf(i) = 0.
951	ENDDO
952
953	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
954	DO i = 1, klon
955	ibas_con(i) = 1
956	itop_con(i) = klev+1
957	ENDDO
958	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
959
960	ENDIF
961
962	c34EK
963	IF (ok_orodr) THEN
964	DO i=1,klon
965	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
966	ENDDO
967	CALL SUGWD(klon,klev,paprs,pplay)
968	DO i=1,klon
969	zuthe(i)=0.
970	zvthe(i)=0.
971	if(zstd(i).gt.10.)then
972	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
973	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
974	endif
975	ENDDO
976	ENDIF
977	c
978	c
979	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
980	PRINT*,'La frequence de lecture surface est de ', lmt_pas
981	c
982	ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
983	IF (ok_mensuel) THEN
984	PRINT*, 'La frequence de sortie mensuelle est de ', ecrit_mth
985	ENDIF
986	ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
987	IF (ok_journe) THEN
988	PRINT*, 'La frequence de sortie journaliere est de ',ecrit_day
989	ENDIF
990	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
991	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
992	ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
993	ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
994	IF (ok_instan) THEN
995	PRINT*, 'La frequence de sortie instant. est de ', ecrit_ins
996	ENDIF
997	ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
998	IF (ok_region) THEN
999	PRINT*, 'La frequence de sortie region est de ', ecrit_reg
1000	ENDIF
1001
1002	c
1003	c Initialiser le couplage si necessaire
1004	c
1005	npas = 0
1006	nexca = 0
1007	if (ocean == 'couple') then
1008	npas = itaufin/ iphysiq
1009	nexca = 86400 / dtime
1010	write(,)' ##### Ocean couple #####'
1011	write(,)' Valeurs des pas de temps'
1012	write(,)' npas = ', npas
1013	write(,)' nexca = ', nexca
1014	endif
1015	c
1016	c
1017	cccIM
1018	capemaxcels = 't_max(X)'
1019	t2mincels = 't_min(X)'
1020	t2maxcels = 't_max(X)'
1021
1022	cccIM cf. FH
1023	c
1024	c=============================================================
1025	c Initialisation des sorties
1026	c=============================================================
1027
1028	#ifdef histISCCP
1029	#include "ini_histISCCP.h"
1030	#endif
1031
1032	#ifdef histhf
1033	#include "ini_histhf.h"
1034	#endif
1035
1036	#include "ini_histday.h"
1037	#include "ini_histmth.h"
1038	#include "ini_histins.h"
1039
1040	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1041	date0 = zjulian
1042	C date0 = day_ini
1043	WRITE(,) 'physiq date0 : ',date0
1044	c
1045	c
1046	c
1047	c Prescrire l'ozone dans l'atmosphere
1048	c
1049	c
1050	cc DO i = 1, klon
1051	cc DO k = 1, klev
1052	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1053	cc ENDDO
1054	cc ENDDO
1055	c
1056	c
1057	ENDIF
1058	c
1059	c ************** Fin de IF ( debut ) *************
1060	c
1061	c
1062	c Mettre a zero des variables de sortie (pour securite)
1063	c
1064	DO i = 1, klon
1065	d_ps(i) = 0.0
1066	ENDDO
1067	DO k = 1, klev
1068	DO i = 1, klon
1069	d_t(i,k) = 0.0
1070	d_u(i,k) = 0.0
1071	d_v(i,k) = 0.0
1072	ENDDO
1073	ENDDO
1074	DO iq = 1, nqmax
1075	DO k = 1, klev
1076	DO i = 1, klon
1077	d_qx(i,k,iq) = 0.0
1078	ENDDO
1079	ENDDO
1080	ENDDO
1081	c
1082	c Ne pas affecter les valeurs entrees de u, v, h, et q
1083	c
1084	DO k = 1, klev
1085	DO i = 1, klon
1086	t_seri(i,k) = t(i,k)
1087	u_seri(i,k) = u(i,k)
1088	v_seri(i,k) = v(i,k)
1089	q_seri(i,k) = qx(i,k,ivap)
1090	ql_seri(i,k) = qx(i,k,iliq)
1091	qs_seri(i,k) = 0.
1092	ENDDO
1093	ENDDO
1094	IF (nqmax.GE.3) THEN
1095	DO iq = 3, nqmax
1096	DO k = 1, klev
1097	DO i = 1, klon
1098	tr_seri(i,k,iq-2) = qx(i,k,iq)
1099	ENDDO
1100	ENDDO
1101	ENDDO
1102	ELSE
1103	DO k = 1, klev
1104	DO i = 1, klon
1105	tr_seri(i,k,1) = 0.0
1106	ENDDO
1107	ENDDO
1108	ENDIF
1109	C
1110	DO i = 1, klon
1111	ztsol(i) = 0.
1112	ENDDO
1113	DO nsrf = 1, nbsrf
1114	DO i = 1, klon
1115	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1116	ENDDO
1117	ENDDO
1118	C
1119	IF (if_ebil.ge.1) THEN
1120	ztit='after dynamic'
1121	CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
1122	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1123	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1124	C Comme les tendances de la physique sont ajoute dans la dynamique,
1125	C on devrait avoir que la variation d'entalpie par la dynamique
1126	C est egale a la variation de la physique au pas de temps precedent.
1127	C Donc la somme de ces 2 variations devrait etre nulle.
1128	call diagphy(paire,ztit,ip_ebil
1129	e , zero_v, zero_v, zero_v, zero_v, zero_v
1130	e , zero_v, zero_v, zero_v, ztsol
1131	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1132	s , fs_bound, fq_bound )
1133	END IF
1134
1135	c Diagnostiquer la tendance dynamique
1136	c
1137	IF (ancien_ok) THEN
1138	DO k = 1, klev
1139	DO i = 1, klon
1140	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1141	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1142	ENDDO
1143	ENDDO
1144	ELSE
1145	DO k = 1, klev
1146	DO i = 1, klon
1147	d_t_dyn(i,k) = 0.0
1148	d_q_dyn(i,k) = 0.0
1149	ENDDO
1150	ENDDO
1151	ancien_ok = .TRUE.
1152	ENDIF
1153	c
1154	c Ajouter le geopotentiel du sol:
1155	c
1156	DO k = 1, klev
1157	DO i = 1, klon
1158	zphi(i,k) = pphi(i,k) + pphis(i)
1159	ENDDO
1160	ENDDO
1161	c
1162	c Verifier les temperatures
1163	c
1164	CALL hgardfou(t_seri,ftsol,'debutphy')
1165	c
1166	c Incrementer le compteur de la physique
1167	c
1168	itap = itap + 1
1169	julien = MOD(NINT(xjour),360)
1170	c
1171	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1172	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1173	c
1174	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1175	PRINT *,' PHYS cond julien ',julien
1176	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1177	ENDIF
1178	c
1179	c Re-evaporer l'eau liquide nuageuse
1180	c
1181	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1182	DO i = 1, klon
1183	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1184	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1185	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1186	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1187	zb = MAX(0.0,ql_seri(i,k))
1188	za = - MAX(0.0,ql_seri(i,k))
1189	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1190	t_seri(i,k) = t_seri(i,k) + za
1191	q_seri(i,k) = q_seri(i,k) + zb
1192	ql_seri(i,k) = 0.0
1193	d_t_eva(i,k) = za
1194	d_q_eva(i,k) = zb
1195	ENDDO
1196	ENDDO
1197	c
1198	IF (if_ebil.ge.2) THEN
1199	ztit='after reevap'
1200	CALL diagetpq(paire,ztit,ip_ebil,2,1,dtime
1201	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1202	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1203	call diagphy(paire,ztit,ip_ebil
1204	e , zero_v, zero_v, zero_v, zero_v, zero_v
1205	e , zero_v, zero_v, zero_v, ztsol
1206	e , d_h_vcol, d_qt, d_ec
1207	s , fs_bound, fq_bound )
1208	C
1209	END IF
1210	C
1211	c
1212	c Appeler la diffusion verticale (programme de couche limite)
1213	c
1214	DO i = 1, klon
1215	c if (.not. ok_veget) then
1216	c frugs(i,is_ter) = SQRT(frugs(i,is_ter)2+rugoro(i)2)
1217	c endif
1218	c frugs(i,is_lic) = rugoro(i)
1219	c frugs(i,is_oce) = rugmer(i)
1220	c frugs(i,is_sic) = 0.001
1221	zxrugs(i) = 0.0
1222	ENDDO
1223	DO nsrf = 1, nbsrf
1224	DO i = 1, klon
1225	c frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1226	frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
1227	ENDDO
1228	ENDDO
1229	DO nsrf = 1, nbsrf
1230	DO i = 1, klon
1231	zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1232	ENDDO
1233	ENDDO
1234	c
1235	C calculs necessaires au calcul de l'albedo dans l'interface
1236	c
1237	CALL orbite(FLOAT(julien),zlongi,dist)
1238	IF (cycle_diurne) THEN
1239	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1240	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1241	ELSE
1242	rmu0 = -999.999
1243	ENDIF
1244	cIM BEG
1245	DO i=1, klon
1246	sunlit(i)=1
1247	IF(rmu0(i).EQ.0.) sunlit(i)=0
1248	c IF(rmu0(i).EQ.0.) THEN
1249	c sunlit(i)=0
1250	c PRINT*,' il fait nuit ',i,rlat(i),rlon(i)
1251	c ENDIF
1252	ENDDO
1253	cIM END
1254	C Calcul de l'abedo moyen par maille
1255	albsol(:)=0.
1256	albsollw(:)=0.
1257	DO nsrf = 1, nbsrf
1258	DO i = 1, klon
1259	albsol(i) = albsol(i) + falbe(i,nsrf) * pctsrf(i,nsrf)
1260	albsollw(i) = albsollw(i) + falblw(i,nsrf) * pctsrf(i,nsrf)
1261	ENDDO
1262	ENDDO
1263	C
1264	C Repartition sous maille des flux LW et SW
1265	C Modif OM+PASB+JLD
1266	C Repartition du longwave par sous-surface linearisee
1267	Cn
1268	DO nsrf = 1, nbsrf
1269	DO i = 1, klon
1270	c$$$ fsollw(i,nsrf) = sollwdown(i) - RSIGMAftsol(i,nsrf)*4
1271	c$$$ fsollw(i,nsrf) = sollw(i)
1272	fsollw(i,nsrf) = sollw(i)
1273	$ + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ftsol(i,nsrf))
1274	fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
1275	ENDDO
1276	ENDDO
1277
1278	fder = dlw
1279
1280	CALL clmain(dtime,itap,date0,pctsrf,
1281	e t_seri,q_seri,u_seri,v_seri,
1282	e julien, rmu0,
1283	e ok_veget, ocean, npas, nexca, ftsol,
1284	$ soil_model,cdmmax, cdhmax, ftsoil, qsol,
1285	$ paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
1286	$ fluxlat,
1287	cIM cf. JLD e rain_fall, snow_fall, solsw, sollw, sollwdown, fder,
1288	e rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder,
1289	e rlon, rlat, cufi, cvfi, frugs,
1290	e debut, lafin, agesno,rugoro ,
1291	s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1292	s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
1293	s dsens, devap,
1294	cIM cf JLD s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m)
1295	s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m,
1296	s fqcalving,ffonte)
1297	c
1298	CXXX PB
1299	CXXX Incrementation des flux
1300	CXXX
1301	zxfluxt=0.
1302	zxfluxq=0.
1303	zxfluxu=0.
1304	zxfluxv=0.
1305	DO nsrf = 1, nbsrf
1306	DO k = 1, klev
1307	DO i = 1, klon
1308	zxfluxt(i,k) = zxfluxt(i,k) +
1309	$ fluxt(i,k,nsrf) * pctsrf( i, nsrf)
1310	zxfluxq(i,k) = zxfluxq(i,k) +
1311	$ fluxq(i,k,nsrf) * pctsrf( i, nsrf)
1312	zxfluxu(i,k) = zxfluxu(i,k) +
1313	$ fluxu(i,k,nsrf) * pctsrf( i, nsrf)
1314	zxfluxv(i,k) = zxfluxv(i,k) +
1315	$ fluxv(i,k,nsrf) * pctsrf( i, nsrf)
1316	END DO
1317	END DO
1318	END DO
1319	DO i = 1, klon
1320	sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1321	c evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1322	evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1323	fder(i) = dlw(i) + dsens(i) + devap(i)
1324	ENDDO
1325
1326
1327	DO k = 1, klev
1328	DO i = 1, klon
1329	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1330	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1331	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1332	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1333	ENDDO
1334	ENDDO
1335	c
1336	IF (if_ebil.ge.2) THEN
1337	ztit='after clmain'
1338	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1339	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1340	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1341	call diagphy(paire,ztit,ip_ebil
1342	e , zero_v, zero_v, zero_v, zero_v, sens
1343	e , evap , zero_v, zero_v, ztsol
1344	e , d_h_vcol, d_qt, d_ec
1345	s , fs_bound, fq_bound )
1346	END IF
1347	C
1348	c
1349	c Incrementer la temperature du sol
1350	c
1351	DO i = 1, klon
1352	zxtsol(i) = 0.0
1353	zxfluxlat(i) = 0.0
1354	cccIM
1355	zt2m(i) = 0.0
1356	zq2m(i) = 0.0
1357	zu10m(i) = 0.0
1358	zv10m(i) = 0.0
1359	cIM cf JLD ??
1360	zxffonte(i) = 0.0
1361	zxfqcalving(i) = 0.0
1362	c
1363	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1364	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
1365	$ THEN
1366	WRITE(,) 'physiq : pb sous surface au point ', i,
1367	$ pctsrf(i, 1 : nbsrf)
1368	ENDIF
1369	ENDDO
1370	DO nsrf = 1, nbsrf
1371	DO i = 1, klon
1372	c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1373	ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
1374	cIM cf. JLD
1375	wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
1376	$ + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1377	zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1378	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
1379	cccIM
1380	zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
1381	zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
1382	zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
1383	zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
1384	cIM cf JLD ??
1385	zxffonte(i) = zxffonte(i) + ffonte(i,nsrf)*pctsrf(i,nsrf)
1386	zxfqcalving(i) = zxfqcalving(i) +
1387	. fqcalving(i,nsrf)*pctsrf(i,nsrf)
1388	c ENDIF
1389	ENDDO
1390	ENDDO
1391
1392	c
1393	c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1394	c
1395	DO nsrf = 1, nbsrf
1396	DO i = 1, klon
1397	IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
1398	cccIM
1399	IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i)
1400	IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i)
1401	IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i)
1402	IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i)
1403	cIM cf JLD ??
1404	IF (pctsrf(i,nsrf) .LT. epsfra) ffonte(i,nsrf) = zxffonte(i)
1405	IF (pctsrf(i,nsrf) .LT. epsfra)
1406	. fqcalving(i,nsrf) = zxfqcalving(i)
1407	ENDDO
1408	ENDDO
1409	c
1410	c
1411	c Calculer la derive du flux infrarouge
1412	c
1413	cXXX DO nsrf = 1, nbsrf
1414	DO i = 1, klon
1415	cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1416	dlw(i) = - 4.0RSIGMAzxtsol(i)**3
1417	cXXX . *(ftsol(i,nsrf)-zxtsol(i))
1418	cXXX . *pctsrf(i,nsrf)
1419	cXXX ENDIF
1420	cXXX ENDDO
1421	ENDDO
1422	c
1423	c Appeler la convection (au choix)
1424	c
1425	DO k = 1, klev
1426	DO i = 1, klon
1427	conv_q(i,k) = d_q_dyn(i,k)
1428	. + d_q_vdf(i,k)/dtime
1429	conv_t(i,k) = d_t_dyn(i,k)
1430	. + d_t_vdf(i,k)/dtime
1431	ENDDO
1432	ENDDO
1433	IF (check) THEN
1434	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1435	PRINT*, "avantcon=", za
1436	ENDIF
1437	zx_ajustq = .FALSE.
1438	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1439	IF (zx_ajustq) THEN
1440	DO i = 1, klon
1441	z_avant(i) = 0.0
1442	ENDDO
1443	DO k = 1, klev
1444	DO i = 1, klon
1445	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1446	. *(paprs(i,k)-paprs(i,k+1))/RG
1447	ENDDO
1448	ENDDO
1449	ENDIF
1450	IF (iflag_con.EQ.1) THEN
1451	stop'reactiver le call conlmd dans physiq.F'
1452	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1453	c . d_t_con, d_q_con,
1454	c . rain_con, snow_con, ibas_con, itop_con)
1455	ELSE IF (iflag_con.EQ.2) THEN
1456	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1457	e conv_t, conv_q, zxfluxq(1,1), omega,
1458	s d_t_con, d_q_con, rain_con, snow_con,
1459	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1460	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1461	WHERE (rain_con < 0.) rain_con = 0.
1462	WHERE (snow_con < 0.) snow_con = 0.
1463	DO i = 1, klon
1464	ibas_con(i) = klev+1 - kcbot(i)
1465	itop_con(i) = klev+1 - kctop(i)
1466	ENDDO
1467	ELSE IF (iflag_con.GE.3) THEN
1468	c nb of tracers for the KE convection:
1469	if (nqmax .GE. 4) then
1470	ntra = nbtr
1471	else
1472	ntra = 1
1473	endif
1474	c
1475	c sb, oct02:
1476	c Schema de convection modularise et vectorise:
1477	c (driver commun aux versions 3 et 4)
1478	c
1479	IF (ok_cvl) THEN ! new driver for convectL
1480
1481	CALL concvl (iflag_con,
1482	. dtime,paprs,pplay,t_seri,q_seri,
1483	. u_seri,v_seri,tr_seri,nbtr,
1484	. ema_work1,ema_work2,
1485	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1486	. rain_con, snow_con, ibas_con, itop_con,
1487	. upwd,dnwd,dnwd0,
1488	. Ma,cape,tvp,iflagctrl,
1489	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd)
1490	cIM cf. FH
1491	clwcon0=qcondc
1492
1493	ELSE ! ok_cvl
1494
1495	c print*,'Avant conema OUI'
1496	CALL conema3 (dtime,
1497	. paprs,pplay,t_seri,q_seri,
1498	. u_seri,v_seri,tr_seri,nbtr,
1499	. ema_work1,ema_work2,
1500	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1501	. rain_con, snow_con, ibas_con, itop_con,
1502	. upwd,dnwd,dnwd0,bas,top,
1503	. Ma,cape,tvp,rflag,
1504	. pbase
1505	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
1506	. ,clwcon0)
1507	print*,'Apres conema3 '
1508
1509	ENDIF ! ok_cvl
1510
1511	IF (.NOT. ok_gust) THEN
1512	do i = 1, klon
1513	wd(i)=0.0
1514	enddo
1515	ENDIF
1516
1517	c =================================================================== c
1518	c Calcul des proprietes des nuages convectifs
1519	c
1520	DO k = 1, klev
1521	DO i = 1, klon
1522	zx_t = t_seri(i,k)
1523	IF (thermcep) THEN
1524	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1525	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1526	zx_qs = MIN(0.5,zx_qs)
1527	zcor = 1./(1.-retv*zx_qs)
1528	zx_qs = zx_qs*zcor
1529	ELSE
1530	IF (zx_t.LT.t_coup) THEN
1531	zx_qs = qsats(zx_t)/pplay(i,k)
1532	ELSE
1533	zx_qs = qsatl(zx_t)/pplay(i,k)
1534	ENDIF
1535	ENDIF
1536	zqsat(i,k)=zx_qs
1537	ENDDO
1538	ENDDO
1539
1540	c calcul des proprietes des nuages convectifs
1541	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
1542	call clouds_gno
1543	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
1544
1545	c =================================================================== c
1546
1547	DO i = 1, klon
1548	ema_pcb(i) = pbase(i)
1549	ENDDO
1550	DO i = 1, klon
1551	ema_pct(i) = paprs(i,itop_con(i))
1552	ENDDO
1553	DO i = 1, klon
1554	ema_cbmf(i) = ema_workcbmf(i)
1555	ENDDO
1556	ELSE
1557	PRINT*, "iflag_con non-prevu", iflag_con
1558	CALL abort
1559	ENDIF
1560
1561	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
1562	c . d_u_con, d_v_con)
1563	DO k = 1, klev
1564	DO i = 1, klon
1565	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
1566	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
1567	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
1568	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
1569	ENDDO
1570	ENDDO
1571	c
1572	IF (if_ebil.ge.2) THEN
1573	ztit='after convect'
1574	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1575	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1576	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1577	call diagphy(paire,ztit,ip_ebil
1578	e , zero_v, zero_v, zero_v, zero_v, zero_v
1579	e , zero_v, rain_con, snow_con, ztsol
1580	e , d_h_vcol, d_qt, d_ec
1581	s , fs_bound, fq_bound )
1582	END IF
1583	C
1584	IF (check) THEN
1585	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1586	PRINT*, "aprescon=", za
1587	zx_t = 0.0
1588	za = 0.0
1589	DO i = 1, klon
1590	za = za + paire(i)/FLOAT(klon)
1591	zx_t = zx_t + (rain_con(i)+snow_con(i))*paire(i)/FLOAT(klon)
1592	ENDDO
1593	zx_t = zx_t/za*dtime
1594	PRINT*, "Precip=", zx_t
1595	ENDIF
1596	IF (zx_ajustq) THEN
1597	DO i = 1, klon
1598	z_apres(i) = 0.0
1599	ENDDO
1600	DO k = 1, klev
1601	DO i = 1, klon
1602	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
1603	. *(paprs(i,k)-paprs(i,k+1))/RG
1604	ENDDO
1605	ENDDO
1606	DO i = 1, klon
1607	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
1608	. /z_apres(i)
1609	ENDDO
1610	DO k = 1, klev
1611	DO i = 1, klon
1612	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
1613	. z_factor(i).LT.(1.0-1.0E-08)) THEN
1614	q_seri(i,k) = q_seri(i,k) * z_factor(i)
1615	ENDIF
1616	ENDDO
1617	ENDDO
1618	ENDIF
1619	zx_ajustq=.FALSE.
1620	c
1621	IF (nqmax.GT.2) THEN !--melange convectif de traceurs
1622	c
1623	IF (iflag_con .NE. 2 .AND. debut) THEN
1624	PRINT*, 'Pour l instant, seul conflx fonctionne ',
1625	$ 'avec traceurs', iflag_con
1626	PRINT*,' Mettre iflag_con',
1627	$ ' = 2 dans run.def et repasser'
1628	c CALL abort
1629	ENDIF
1630	c
1631	ENDIF !--nqmax.GT.2
1632	c
1633	c Appeler l'ajustement sec
1634	c
1635	CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
1636	DO k = 1, klev
1637	DO i = 1, klon
1638	t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
1639	q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
1640	ENDDO
1641	ENDDO
1642	c
1643	IF (if_ebil.ge.2) THEN
1644	ztit='after dry_adjust'
1645	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1646	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1647	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1648	END IF
1649
1650
1651	c-------------------------------------------------------------------------
1652	c Caclul des ratqs
1653	c-------------------------------------------------------------------------
1654
1655	c print*,'calcul des ratqs'
1656	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
1657	c ----------------
1658	c on ecrase le tableau ratqsc calcule par clouds_gno
1659	if (iflag_cldcon.eq.1) then
1660	do k=1,klev
1661	do i=1,klon
1662	if(ptconv(i,k)) then
1663	ratqsc(i,k)=ratqsbas
1664	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
1665	else
1666	ratqsc(i,k)=0.
1667	endif
1668	enddo
1669	enddo
1670	endif
1671
1672	c ratqs stables
1673	c -------------
1674	do k=1,klev
1675	ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
1676	s min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
1677	enddo
1678
1679
1680	c ratqs final
1681	c -----------
1682	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
1683	c les ratqs sont une conbinaison de ratqss et ratqsc
1684	c ratqs final
1685	c 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1686	c relaxation des ratqs
1687	c facttemps=exp(-pdtphys/1.e4)
1688	facteur=exp(-pdtphys*facttemps)
1689	ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
1690	ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
1691	c print*,'calcul des ratqs fini'
1692	else
1693	c on ne prend que le ratqs stable pour fisrtilp
1694	ratqs(:,:)=ratqss(:,:)
1695	endif
1696
1697
1698	c
1699	c Appeler le processus de condensation a grande echelle
1700	c et le processus de precipitation
1701	c-------------------------------------------------------------------------
1702	CALL fisrtilp(dtime,paprs,pplay,
1703	. t_seri, q_seri,ptconv,ratqs,
1704	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
1705	. rain_lsc, snow_lsc,
1706	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
1707	. frac_impa, frac_nucl,
1708	. prfl, psfl, rhcl)
1709
1710	WHERE (rain_lsc < 0) rain_lsc = 0.
1711	WHERE (snow_lsc < 0) snow_lsc = 0.
1712	DO k = 1, klev
1713	DO i = 1, klon
1714	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
1715	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
1716	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
1717	cldfra(i,k) = rneb(i,k)
1718	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
1719	ENDDO
1720	ENDDO
1721	IF (check) THEN
1722	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1723	PRINT*, "apresilp=", za
1724	zx_t = 0.0
1725	za = 0.0
1726	DO i = 1, klon
1727	za = za + paire(i)/FLOAT(klon)
1728	zx_t = zx_t + (rain_lsc(i)+snow_lsc(i))*paire(i)/FLOAT(klon)
1729	ENDDO
1730	zx_t = zx_t/za*dtime
1731	PRINT*, "Precip=", zx_t
1732	ENDIF
1733	c
1734	IF (if_ebil.ge.2) THEN
1735	ztit='after fisrt'
1736	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1737	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1738	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1739	call diagphy(paire,ztit,ip_ebil
1740	e , zero_v, zero_v, zero_v, zero_v, zero_v
1741	e , zero_v, rain_lsc, snow_lsc, ztsol
1742	e , d_h_vcol, d_qt, d_ec
1743	s , fs_bound, fq_bound )
1744	END IF
1745	c
1746	c-------------------------------------------------------------------
1747	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1748	c-------------------------------------------------------------------
1749
1750	c 1. NUAGES CONVECTIFS
1751	c
1752	IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
1753
1754	c Nuages diagnostiques pour Tiedtke
1755	CALL diagcld1(paprs,pplay,
1756	. rain_con,snow_con,ibas_con,itop_con,
1757	. diafra,dialiq)
1758	DO k = 1, klev
1759	DO i = 1, klon
1760	IF (diafra(i,k).GT.cldfra(i,k)) THEN
1761	cldliq(i,k) = dialiq(i,k)
1762	cldfra(i,k) = diafra(i,k)
1763	ENDIF
1764	ENDDO
1765	ENDDO
1766
1767	ELSE IF (iflag_cldcon.eq.3) THEN
1768	c On prend pour les nuages convectifs le max du calcul de la
1769	c convection et du calcul du pas de temps précédent diminué d'un facteur
1770	c facttemps
1771	c facttemps=pdtphys/1.e4
1772	facteur = pdtphys *facttemps
1773	do k=1,klev
1774	do i=1,klon
1775	rnebcon(i,k)=rnebcon(i,k)*facteur
1776	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
1777	s then
1778	rnebcon(i,k)=rnebcon0(i,k)
1779	clwcon(i,k)=clwcon0(i,k)
1780	endif
1781	enddo
1782	enddo
1783
1784	cIM ISCCP simulator BEGIN
1785	IF (ok_isccp) THEN
1786	cIM calcul tau. emi nuages convectifs
1787	convfra(:,:)=rnebcon(:,:)
1788	convliq(:,:)=rnebcon(:,:)*clwcon(:,:)
1789	c CALL newmicro (paprs, pplay,ok_newmicro,
1790	c . t_seri, cldliq, cldfra, cldtau, cldemi,
1791	c . cldh, cldl, cldm, cldt, cldq)
1792	CALL newmicro (paprs, pplay,ok_newmicro,
1793	. t_seri, convliq, convfra, dtau_c, dem_c,
1794	. cldh_c, cldl_c, cldm_c, cldt_c, cldq_c)
1795
1796	cIM calcul tau. emi nuages startiformes
1797	CALL newmicro (paprs, pplay,ok_newmicro,
1798	. t_seri, cldliq, cldfra, dtau_s, dem_s,
1799	. cldh_s, cldl_s, cldm_s, cldt_s, cldq_s)
1800	cIM calcul diagramme (PC, tau) cf. ISCCP D
1801	c seed=50
1802	c seed=ran0(klon)
1803	cT1O3
1804	c top_height=1
1805	cT3O3
1806	c top_height=3
1807	c overlap=3
1808	cIM cf GCM
1809	cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
1810
1811	cIM inversion des niveaux de pression ==> de haut en bas
1812	DO k=1,klev
1813	kinv=klev-k+1
1814	DO i=1,klon
1815	pfull(i,k)=pplay(i,kinv)
1816	c on met toutes les variables de Haut 2 Bas
1817	qv(i,k)=q_seri(i,kinv)
1818	cc(i,k)=cldtot(i,kinv)
1819	conv(i,k)=rnebcon(i,kinv)
1820	dtau_sH2B(i,k)=dtau_s(i,kinv)
1821	dtau_cH2B(i,k)=dtau_c(i,kinv)
1822	at(i,k)=t_seri(i,kinv)
1823	dem_sH2B(i,k)=dem_s(i,kinv)
1824	dem_cH2B(i,k)=dem_c(i,kinv)
1825
1826	ENDDO
1827	ENDDO
1828
1829	DO k=1,klev+1
1830	kinv=klev-k+2
1831	DO i=1,klon
1832	phalf(i,k)=paprs(i,kinv)
1833	ENDDO
1834	ENDDO
1835
1836	c open(99,file='tautab.bin',access='sequential',
1837	c $ form='unformatted',status='old')
1838	c read(99) tautab
1839
1840	cIM210503
1841	IF (debut) THEN
1842	open(99,file='tautab.formatted', FORM='FORMATTED')
1843	read(99,'(f30.20)') tautab
1844	close(99)
1845	c
1846	open(99,file='invtau.formatted',form='FORMATTED')
1847	read(99,'(i10)') invtau
1848	close(99)
1849	c
1850	nsrf=3
1851	DO nreg=1, nbreg
1852	DO i=1, klon
1853
1854	c IF (debut) THEN
1855	IF(rlon(i).LT.0.) THEN
1856	rlonPOS(i)=rlon(i)+360.
1857	ELSE
1858	rlonPOS(i)=rlon(i)
1859	ENDIF
1860	c ENDIF
1861
1862	c pct_ocean(i,nreg)=.FALSE.
1863	pct_ocean(i,nreg)=0
1864
1865	c DO nsrf = 1, nbsrf
1866
1867	c test si c'est 1 point d'ocean
1868	IF(pctsrf(i,nsrf).EQ.1.) THEN
1869
1870	IF(nreg.EQ.1) THEN
1871
1872	c TROP
1873	IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
1874	c pct_ocean(i,nreg)=.TRUE.
1875	pct_ocean(i,nreg)=1
1876	ENDIF
1877
1878	c PACIFIQUE NORD
1879	ELSEIF(nreg.EQ.2) THEN
1880	IF(rlat(i).GE.40.AND.rlat(i).LE.60.) THEN
1881	IF(rlonPOS(i).GE.160..AND.rlonPOS(i).LE.235.) THEN
1882	c pct_ocean(i,nreg)=.TRUE.
1883	pct_ocean(i,nreg)=1
1884	ENDIF
1885	ENDIF
1886	c CALIFORNIE ST-CU
1887	ELSEIF(nreg.EQ.3) THEN
1888	IF(rlonPOS(i).GE.220..AND.rlonPOS(i).LE.250.) THEN
1889	IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
1890	c pct_ocean(i,nreg)=.TRUE.
1891	pct_ocean(i,nreg)=1
1892	ENDIF
1893	ENDIF
1894	c HAWAI
1895	ELSEIF(nreg.EQ.4) THEN
1896	IF(rlonPOS(i).GE.180..AND.rlonPOS(i).LE.220.) THEN
1897	IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
1898	c pct_ocean(i,nreg)=.TRUE.
1899	pct_ocean(i,nreg)=1
1900	ENDIF
1901	ENDIF
1902	c WARM POOL
1903	ELSEIF(nreg.EQ.5) THEN
1904	IF(rlonPOS(i).GE.70..AND.rlonPOS(i).LE.150.) THEN
1905	IF(rlat(i).GE.-5.AND.rlat(i).LE.20.) THEN
1906	c pct_ocean(i,nreg)=.TRUE.
1907	pct_ocean(i,nreg)=1
1908	ENDIF
1909	ENDIF
1910	ENDIF !nbreg
1911	c TROP
1912	c IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
1913	c pct_ocean(i)=.TRUE.
1914	c WRITE(,) 'pct_ocean =',i, rlon(i), rlat(i)
1915	c ENDIF !lon
1916	c ENDIF !lat
1917
1918	ENDIF !pctsrf
1919	c ENDDO
1920	ENDDO !klon
1921	ENDDO !nbreg
1922
1923	cIM somme de toutes les nhistoW BEG
1924	DO nreg = 1, nbreg
1925	DO k = 1, kmaxm1
1926	DO l = 1, lmaxm1
1927	DO iw = 1, iwmax
1928	nhistoWt(k,l,iw,nreg)=0.
1929	ENDDO
1930	ENDDO
1931	ENDDO
1932	ENDDO
1933	cIM somme de toutes les nhistoW END
1934	ENDIF
1935
1936
1937	c CALL ISCCP_CLOUD_TYPES(nlev,ncol,seed,pfull,phalf,qv,
1938	c & cc,conv,dtau_s,dtau_c,top_height,overlap,
1939	c & tautab,invtau,skt,emsfc_lw,at,dem_s,dem_c,fq_isccp,
1940	c & totalcldarea,meanptop,meantaucld,boxtau,boxptop)
1941
1942	c DO i=1, klon
1943	c i=1
1944	c1011 CONTINUE
1945	c
1946	cIM on verifie les donnees de INPUT en dehors du simulateur ISCCP
1947	cIM 1D non-vectorise (!) pour qu'on gagne du temps ...
1948	cIM
1949	c BEGIN find unpermittable data.....
1950	! ---------------------------------------------------!
1951	! find unpermittable data.....
1952	!
1953	c do 13 k=1,klev
1954	c ca prend trop de temps ??
1955	c cldtot(:,:) = min(max(cldtot(:,:),0.),1.)
1956	c rnebcon(:,:) = min(max(rnebcon(:,:),0.),1.)
1957	c dtau_s(:,:) = max(dtau_s(:,:),0.)
1958	c dem_s(:,:) = min(max(dem_s(:,:),0.),1.)
1959	c dtau_c(:,:) = max(dtau_c(:,:),0.)
1960	c dem_c(:,:) = min(max(dem_c(:,:),0.),1.)
1961	c ca prend trop de temps ??
1962
1963	c if (cldtot(i,k) .lt. 0.) then
1964	c print *, ' error = cloud fraction less than zero'
1965	c STOP
1966	c end if
1967	c if (cldtot(i,k) .gt. 1.) then
1968	c print *, ' error = cloud fraction greater than 1'
1969	c STOP
1970	c end if
1971	c if (rnebcon(i,k) .lt. 0.) then
1972	c print *,
1973	c & ' error = convective cloud fraction less than zero'
1974	c STOP
1975	c end if
1976	c if (rnebcon(i,k) .gt. 1.) then
1977	c print *,
1978	c & ' error = convective cloud fraction greater than 1'
1979	c STOP
1980	c end if
1981
1982	c if (dtau_s(i,k) .lt. 0.) then
1983	c print *,
1984	c & ' error = stratiform cloud opt. depth less than zero'
1985	c STOP
1986	c end if
1987	c if (dem_s(i,k) .lt. 0.) then
1988	c print *,
1989	c & ' error = stratiform cloud emissivity less than zero'
1990	c STOP
1991	c end if
1992	c if (dem_s(i,k) .gt. 1.) then
1993	c print *,
1994	c & ' error = stratiform cloud emissivity greater than 1'
1995	c STOP
1996	c end if
1997
1998	c if (dtau_c(i,k) .lt. 0.) then
1999	c print *,
2000	c & ' error = convective cloud opt. depth less than zero'
2001	c STOP
2002	c end if
2003	c if (dem_c(i,k) .lt. 0.) then
2004	c print *,
2005	c & ' error = convective cloud emissivity less than zero'
2006	c STOP
2007	c end if
2008	c if (dem_c(i,k) .gt. 1.) then
2009	c print *,
2010	c & ' error = convective cloud emissivity greater than 1'
2011	c STOP
2012	c end if
2013	c13 continue
2014
2015	! ---------------------------------------------------!
2016	c
2017	c END find unpermittable data.....
2018	cv2.2.1.1 DO i=1, klon
2019	c i=1
2020	c seed=i
2021	c
2022	cv3.4
2023	if (debut) then
2024	DO i=1, klon
2025	seed(i)=i+100
2026	c seed(i)=i+50
2027	ENDDO
2028	endif
2029	c seed=aint(ran0(klon))
2030	c CALL ISCCP_CLOUD_TYPES(klev,ncol,seed,pfull(i,:),phalf(i,:)
2031	cv2.2.1.1
2032	c CALL ISCCP_CLOUD_TYPES(klev,ncol,seed(i),pfull(i,:),phalf(i,:)
2033	c & ,q_seri(i,:),
2034	c & cldtot(i,:),rnebcon(i,:),dtau_s(i,:),dtau_c(i,:),
2035	c & top_height,overlap,
2036	c & tautab,invtau,ztsol,emsfc_lw,t_seri(i,:),dem_s(i,:),
2037	c & dem_c(i,:),
2038	c & fq_isccp(i,:,:),
2039	c & totalcldarea(i),meanptop(i),meantaucld(i),
2040	c & boxtau(i,:),boxptop(i,:))
2041	cv2.2.1.1
2042	cv3.4
2043	debug=0
2044	debugcol=0
2045	cIM260503
2046	c o500 ==> distribution nuage ftion du regime dynamique
2047	DO i=1, klon
2048	o500(i)=omega(i,8)*864.
2049	c PRINT*,'pphi8 ',pphi(i,8),'zphi8,11,12',zphi(i,8),
2050	c & zphi(i,11),zphi(i,12)
2051	ENDDO
2052
2053	c axe vertical pour les differents niveaux des histogrammes
2054	c DO iw=1, iwmax
2055	c zx_o500(iw)=wmin+(iw-1./2.)*pas_w
2056	c ENDDO
2057	c PRINT*,' phys AVANT seed(3361)=',seed(3361)
2058	CALL ISCCP_CLOUD_TYPES(
2059	& debug,
2060	& debugcol,
2061	& klon,
2062	& sunlit,
2063	& klev,
2064	& ncol,
2065	& seed,
2066	& pfull,
2067	& phalf,
2068	c var de bas en haut ==> PB !
2069	c & q_seri,
2070	c & cldtot,
2071	c & rnebcon,
2072	c & dtau_s,
2073	c & dtau_c,
2074	c var de Haut en Bas BEG
2075	& qv, cc, conv, dtau_sH2B, dtau_cH2B,
2076	c var de Haut en Bas END
2077	& top_height,
2078	& overlap,
2079	& tautab,
2080	& invtau,
2081	& ztsol,
2082	& emsfc_lw,
2083	c var de bas en haut ==> PB !
2084	c & t_seri,
2085	c & dem_s,
2086	c & dem_c,
2087	c var de Haut en Bas BEG
2088	& at, dem_sH2B, dem_cH2B,
2089	cIM260503
2090	c & o500, pct_ocean,
2091	c var de Haut en Bas END
2092	& fq_isccp,
2093	& totalcldarea,
2094	& meanptop,
2095	& meantaucld,
2096	& boxtau,
2097	& boxptop)
2098	c & boxptop,
2099	cIM 260503
2100	c & histoW,
2101	c & nhistoW
2102	c &)
2103
2104	cIM 200603
2105	c PRINT*,'physiq fq_isccp(6,1,1)',fq_isccp(6,1,1)
2106
2107	cIM 200603
2108	cIM somme de toutes les nhistoW BEG
2109	c DO k = 1, kmaxm1
2110	c DO l = 1, lmaxm1
2111	c DO iw = 1, iwmax
2112	c nhistoWt(k,l,iw)=nhistoWt(k,l,iw)+nhistoW(k,l,iw)
2113	ccc IF(k.EQ.1.AND.l.EQ.1.AND.iw.EQ.1) then
2114	c IF(nhistoWt(k,l,iw).NE.0.) THEN
2115	c PRINT*,' physiq nWt', k,l,iw,nhistoWt(k,l,iw)
2116	c ENDIF
2117	c ENDDO
2118	c ENDDO
2119	c ENDDO
2120	cIM somme de toutes les nhistoW END
2121	c PRINT*,' phys APRES seed(3361)=',seed(3361)
2122	cv3.4
2123	c i=i+1
2124	c IF(i.LE.klon) THEN
2125	c GOTO 1011
2126	c ENDIF
2127	cv2.2.1.1 ENDDO
2128
2129	c passage de la grille (klon,7,7) a (iim,jjmp1,7,7)
2130	c minfq3d=100.
2131	c maxfq3d=0.
2132	cIM calcul des correspondances entre la grille physique et
2133	cIM la grille dynamique
2134	c DO i=1, klon
2135	c grid_phys(i)=i
2136	c PRINT*,'i, grid_phys',i,grid_phys(i)
2137	c ENDDO
2138	c CALL gr_fi_dyn(1,klon,iimp1,jjmp1,grid_phys,grid_dyn)
2139	c DO j=1, jjmp1
2140	c DO i=1, iimp1
2141	c PRINT*,'i,j grid_dyn ',i,j,grid_dyn(i,j)
2142	c ENDDO
2143	c ENDDO
2144	c
2145	DO l=1, lmax
2146	DO k=1, kmax
2147	cIM grille physique ==> grille ecriture 2D (iim,jjmp1)
2148	c
2149	DO i=1, iim
2150	fq4d(i,1,k,l)=fq_isccp(1,k,l)
2151	cc PRINT*,'first j=1 i =',i
2152	ENDDO
2153	DO j=2, jjm
2154	DO i=1, iim
2155	cERROR ?? ig=i+iim*(j-1)
2156	ig=i+1+(j-2)*iim
2157	cc PRINT*,'i =',i,'j =',j,'ig=',ig
2158	fq4d(i,j,k,l)=fq_isccp(ig,k,l)
2159	ENDDO
2160	ENDDO
2161	DO i=1, iim
2162	fq4d(i,jjmp1,k,l)=fq_isccp(klon,k,l)
2163	cc PRINT*,'last jjmp1 i =',i
2164	ENDDO
2165	IF(debut) THEN
2166	DO j=1, jjmp1
2167	DO i=1, iim
2168	IF(j.GE.2.AND.j.LE.jjm) THEN
2169	igfi2D(i,j)=i+1+(j-2)*iim
2170	c PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
2171	ELSEIF(j.EQ.1) THEN
2172	igfi2D(i,j)=1
2173	c PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
2174	ELSEIF(j.EQ.jjmp1) THEN
2175	igfi2D(i,j)=klon
2176	c PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
2177	ENDIF
2178	ENDDO
2179	ENDDO
2180	ENDIF
2181	c STOP
2182	c
2183	c CALL gr_fi_ecrit(1,klon,iim,jjmp1,fq_isccp(:,k,l),
2184	c $ fq4d(:,:,k,l))
2185	ENDDO
2186	ENDDO
2187	DO l=1, lmax
2188	fq4d(:,:,8,l)=-1.e+10
2189	fq4d(:,:,l,8)=-1.e+10
2190	ENDDO
2191	DO l=1, lmax
2192	DO k=1, kmax
2193	DO j=1, jjmp1
2194	DO i=1, iim
2195
2196	c inversion TAU ?!
2197	c ni=(i-1)*lmax+l
2198	c nj=(j-1)*kmax+kmax-k+1
2199	c
2200	c210503 inversion en PC ==> pas besoin !!!
2201	c ni=(i-1)*lmax+lmax-l+1
2202	c nj=(j-1)*kmax+k
2203	c
2204	c210503
2205	ni=(i-1)*lmax+l
2206	nj=(j-1)*kmax+k
2207
2208	c210503 if(k.EQ.8) then
2209	c fq4d(i,j,8,l)=-1.e+10
2210	c endif
2211
2212	c210503 if(l.EQ.8) then
2213	c fq4d(i,j,k,8)=-1.e+10
2214	c endif
2215
2216	fq3d(ni,nj)=fq4d(i,j,k,l)
2217
2218	c if(fq3d(ni,nj).lt.0.) then
2219	c print*,' fq3d LT ZERO ',ni,nj,fq3d(ni,nj)
2220	c endif
2221	c if(fq3d(ni,nj).gt.100.) then
2222	c print*,' fq3d GT 100 ',ni,nj,fq3d(ni,nj)
2223	c endif
2224	c max & min fq3d
2225	c if(fq3d(ni,nj).gt.maxfq3d) maxfq3d=fq3d(ni,nj)
2226	c if(fq3d(ni,nj).lt.minfq3d) minfq3d=fq3d(ni,nj)
2227
2228	ENDDO
2229	ENDDO
2230	c fq4d(:,:,8,l)=-1.e+10
2231	c fq4d(:,:,k,8)=-1.e+10
2232	c k=k+1
2233	c if(k.LE.kmax) then
2234	c goto 1022
2235	c endif
2236	ENDDO
2237	c l=l+1
2238	c if(l.LE.lmax) then
2239	c goto 1021
2240	c endif
2241	ENDDO
2242
2243	c print*,' minfq3d=',minfq3d,' maxfq3d=',maxfq3d
2244	c
2245	c calculs statistiques distribution nuage ftion du regime dynamique
2246	c DO i=1, klon
2247	c! o500(i)=omega(i,9)*864.
2248	c! PRINT*,' o500=',o500(i),' pphi(9)=',pphi(i,9)
2249	c o500(i)=omega(i,8)*864.
2250	cc PRINT*,' pphi(8)',pphi(i,8),'pphi(11)',pphi(i,11),
2251	cc .'pphi(12)',pphi(i,12)
2252	cc PRINT*,' zphi8,11,12=',zphi(i,8),zphi(i,11),zphi(i,12)
2253	cc PRINT*,' o500',o500(i),' w500',w500(i)
2254	c ENDDO
2255
2256	c axe vertical pour les differents niveaux des histogrammes
2257	c DO iw=1, iwmax
2258	c zx_o500(iw)=wmin+(iw-1./2.)*pas_w
2259	c ENDDO
2260
2261
2262	c Ce calcul doit etre fait a partir de valeurs mensuelles ??
2263	cc CALL histo_o500_pctau(o500,fq4d,histoW)
2264	cc CALL histo_o500_pctau(paire,pctsrf,o500,fq4d,histoW)
2265	cc CALL histo_o500_pctau(pct_ocean,rlat,o500,fq4d,histoW)
2266	ccOK ??? CALL histo_o500_pctau(pct_ocean,o500,fq4d,histoW)
2267	c CALL histo_o500_pctau(klon,pct_ocean,o500,fq4d,histoW,nhistoW)
2268	c CALL histo_o500_pctau(klon,pct_ocean,o500,fq_isccp,
2269	CALL histo_o500_pctau(nbreg,pct_ocean,o500,fq_isccp,
2270	&histoW,nhistoW)
2271	c
2272	cIM somme de toutes les nhistoW BEG
2273	DO nreg=1, nbreg
2274	DO k = 1, kmaxm1
2275	DO l = 1, lmaxm1
2276	DO iw = 1, iwmax
2277	nhistoWt(k,l,iw,nreg)=nhistoWt(k,l,iw,nreg)+
2278	& nhistoW(k,l,iw,nreg)
2279	ccc IF(k.EQ.1.AND.l.EQ.1.AND.iw.EQ.1) then
2280	c IF(nhistoWt(k,l,iw).NE.0.) THEN
2281	c PRINT*,' physiq nWt', k,l,iw,nhistoWt(k,l,iw)
2282	c ENDIF
2283	ENDDO
2284	ENDDO
2285	ENDDO
2286	ENDDO
2287	cIM somme de toutes les nhistoW END
2288	c
2289	c IF(lafin) THEN
2290	c DO nreg=1, nbreg
2291	c DO iw=1, iwmax
2292	c DO l=1,lmaxm1
2293	c DO k=1,kmaxm1
2294	c IF(histoW(k,l,iw,nreg).NE.0.) then
2295	c PRINT*,'physiq H nH',k,l,iw,
2296	c & histoW(k,l,iw,nreg),
2297	c & nhistoW(k,l,iw,nreg),nhistoWt(k,l,iw,nreg)
2298	c ENDIF
2299	c ENDDO
2300	c ENDDO
2301	c ENDDO
2302	c ENDDO
2303	cIM verif fq_isccp, fq4d, fq3d
2304	c DO l=1, lmaxm1
2305	c DO k=1,kmaxm1
2306	c i=74
2307	c j=36
2308	c DO j=1, jjmp1
2309	c DO i=1, iim
2310	c DO l=1, lmaxm1
2311	c WRITE(*,'(a,3i4,7f10.4)')
2312	c & 'fq_isccp,j,i,l=',j,i,l,
2313	c & (fq_isccp(igfi2D(i,j),k,l),k=1,kmaxm1)
2314	c WRITE(*,'(a,3i4,7f10.4)')
2315	c & 'fq4d,j,i,l=',j,i,l,(fq4d(i,j,k,l),k=1,kmaxm1)
2316	c ENDDO
2317	c ENDDO
2318	c ENDDO
2319	c ni1=(i-1)*8+1
2320	c ni2=i*8
2321	c nj1=(j-1)*8+1
2322	c nj2=j*8
2323	c DO ni=ni1,ni2
2324	c WRITE(*,'(a,2i4,7f10.4)')
2325	c & 'fq3d, ni,nj=',ni,nj,
2326	c & (fq3d(ni,nj),nj=nj1,nj2)
2327	c ENDDO
2328	c ENDIF
2329
2330	c DO iw=1, iwmax
2331	c DO l=1,lmaxm1
2332	c DO k=1,kmaxm1
2333	c PRINT*,' iw,l,k,nhistoW=',iw,l,k,nhistoW(k,l,iw)
2334	c ENDDO
2335	c ENDDO
2336	c ENDDO
2337
2338	c DO iw=1, iwmax
2339	c DO l=1, lmaxm1
2340	c linv=lmaxm1-l+1
2341	c DO k=1, kmaxm1
2342	c histoWinv(k,l,iw)=histoW(iw,k,l)
2343	c ENDDO
2344	c ENDDO
2345	c ENDDO
2346	c
2347	c pb syncronisation ?? : 48 * 30 * 7 (jour1) + 48* 29 * 7 (jour suivant)
2348	c
2349
2350
2351	ENDIF !ok_isccp
2352	cIM ISCCP simulator END
2353
2354	c On prend la somme des fractions nuageuses et des contenus en eau
2355	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2356	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2357
2358
2359	ENDIF
2360	c
2361	c 2. NUAGES STARTIFORMES
2362	c
2363	IF (ok_stratus) THEN
2364	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2365	DO k = 1, klev
2366	DO i = 1, klon
2367	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2368	cldliq(i,k) = dialiq(i,k)
2369	cldfra(i,k) = diafra(i,k)
2370	ENDIF
2371	ENDDO
2372	ENDDO
2373	ENDIF
2374	c
2375	c Precipitation totale
2376	c
2377	DO i = 1, klon
2378	rain_fall(i) = rain_con(i) + rain_lsc(i)
2379	snow_fall(i) = snow_con(i) + snow_lsc(i)
2380	ENDDO
2381	c
2382	IF (if_ebil.ge.2) THEN
2383	ztit="after diagcld"
2384	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2385	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2386	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2387	END IF
2388	c
2389	c Calculer l'humidite relative pour diagnostique
2390	c
2391	DO k = 1, klev
2392	DO i = 1, klon
2393	zx_t = t_seri(i,k)
2394	IF (thermcep) THEN
2395	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2396	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2397	zx_qs = MIN(0.5,zx_qs)
2398	zcor = 1./(1.-retv*zx_qs)
2399	zx_qs = zx_qs*zcor
2400	ELSE
2401	IF (zx_t.LT.t_coup) THEN
2402	zx_qs = qsats(zx_t)/pplay(i,k)
2403	ELSE
2404	zx_qs = qsatl(zx_t)/pplay(i,k)
2405	ENDIF
2406	ENDIF
2407	zx_rh(i,k) = q_seri(i,k)/zx_qs
2408	zqsat(i,k)=zx_qs
2409	ENDDO
2410	ENDDO
2411	c
2412	c Calculer les parametres optiques des nuages et quelques
2413	c parametres pour diagnostiques:
2414	c
2415	if (ok_newmicro) then
2416	CALL newmicro (paprs, pplay,ok_newmicro,
2417	. t_seri, cldliq, cldfra, cldtau, cldemi,
2418	. cldh, cldl, cldm, cldt, cldq)
2419	else
2420	CALL nuage (paprs, pplay,
2421	. t_seri, cldliq, cldfra, cldtau, cldemi,
2422	. cldh, cldl, cldm, cldt, cldq)
2423	endif
2424	c
2425	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2426	c
2427	IF (MOD(itaprad,radpas).EQ.0) THEN
2428	DO i = 1, klon
2429	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
2430	. + falbe(i,is_lic) * pctsrf(i,is_lic)
2431	. + falbe(i,is_ter) * pctsrf(i,is_ter)
2432	. + falbe(i,is_sic) * pctsrf(i,is_sic)
2433	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
2434	. + falblw(i,is_lic) * pctsrf(i,is_lic)
2435	. + falblw(i,is_ter) * pctsrf(i,is_ter)
2436	. + falblw(i,is_sic) * pctsrf(i,is_sic)
2437	ENDDO
2438	! if (debut) then
2439	! albsol1 = albsol
2440	! albsollw1 = albsollw
2441	! endif
2442	! albsol = albsol1
2443	! albsollw = albsollw1
2444	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2445	cIM e (dist, rmu0, fract, co2_ppm, solaire,
2446	e (dist, rmu0, fract,
2447	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2448	e wo,
2449	e cldfra, cldemi, cldtau,
2450	s heat,heat0,cool,cool0,radsol,albpla,
2451	s topsw,toplw,solsw,sollw,
2452	s sollwdown,
2453	cccIMs topsw0,toplw0,solsw0,sollw0)
2454	s topsw0,toplw0,solsw0,sollw0,
2455	s swdn0, swdn, swup0, swup )
2456	itaprad = 0
2457	ENDIF
2458	itaprad = itaprad + 1
2459
2460	c
2461	c Ajouter la tendance des rayonnements (tous les pas)
2462	c
2463	DO k = 1, klev
2464	DO i = 1, klon
2465	t_seri(i,k) = t_seri(i,k)
2466	. + (heat(i,k)-cool(i,k)) * dtime/86400.
2467	ENDDO
2468	ENDDO
2469	c
2470	IF (if_ebil.ge.2) THEN
2471	ztit='after rad'
2472	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2473	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2474	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2475	call diagphy(paire,ztit,ip_ebil
2476	e , topsw, toplw, solsw, sollw, zero_v
2477	e , zero_v, zero_v, zero_v, ztsol
2478	e , d_h_vcol, d_qt, d_ec
2479	s , fs_bound, fq_bound )
2480	END IF
2481	c
2482	c
2483	c Calculer l'hydrologie de la surface
2484	c
2485	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2486	c . agesno, ftsol,fqsurf,fsnow, ruis)
2487	c
2488	DO i = 1, klon
2489	zxqsurf(i) = 0.0
2490	zxsnow(i) = 0.0
2491	ENDDO
2492	DO nsrf = 1, nbsrf
2493	DO i = 1, klon
2494	zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf)
2495	zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
2496	ENDDO
2497	ENDDO
2498	c
2499	c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
2500	c
2501	cXXX DO nsrf = 1, nbsrf
2502	cXXX DO i = 1, klon
2503	cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN
2504	cXXX fqsurf(i,nsrf) = zxqsurf(i)
2505	cXXX fsnow(i,nsrf) = zxsnow(i)
2506	cXXX ENDIF
2507	cXXX ENDDO
2508	cXXX ENDDO
2509	c
2510	c Calculer le bilan du sol et la derive de temperature (couplage)
2511	c
2512	DO i = 1, klon
2513	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2514	c a la demande de JLD
2515	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2516	ENDDO
2517	c
2518	cmoddeblott(jan95)
2519	c Appeler le programme de parametrisation de l'orographie
2520	c a l'echelle sous-maille:
2521	c
2522	IF (ok_orodr) THEN
2523	c
2524	c selection des points pour lesquels le shema est actif:
2525	igwd=0
2526	DO i=1,klon
2527	itest(i)=0
2528	c IF ((zstd(i).gt.10.0)) THEN
2529	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2530	itest(i)=1
2531	igwd=igwd+1
2532	idx(igwd)=i
2533	ENDIF
2534	ENDDO
2535	c igwdim=MAX(1,igwd)
2536	c
2537	CALL drag_noro(klon,klev,dtime,paprs,pplay,
2538	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
2539	e igwd,idx,itest,
2540	e t_seri, u_seri, v_seri,
2541	s zulow, zvlow, zustr, zvstr,
2542	s d_t_oro, d_u_oro, d_v_oro)
2543	c
2544	c ajout des tendances
2545	DO k = 1, klev
2546	DO i = 1, klon
2547	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
2548	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
2549	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
2550	ENDDO
2551	ENDDO
2552	c
2553	ENDIF ! fin de test sur ok_orodr
2554	c
2555	IF (ok_orolf) THEN
2556	c
2557	c selection des points pour lesquels le shema est actif:
2558	igwd=0
2559	DO i=1,klon
2560	itest(i)=0
2561	IF ((zpic(i)-zmea(i)).GT.100.) THEN
2562	itest(i)=1
2563	igwd=igwd+1
2564	idx(igwd)=i
2565	ENDIF
2566	ENDDO
2567	c igwdim=MAX(1,igwd)
2568	c
2569	CALL lift_noro(klon,klev,dtime,paprs,pplay,
2570	e rlat,zmea,zstd,zpic,
2571	e itest,
2572	e t_seri, u_seri, v_seri,
2573	s zulow, zvlow, zustr, zvstr,
2574	s d_t_lif, d_u_lif, d_v_lif)
2575	c
2576	c ajout des tendances
2577	DO k = 1, klev
2578	DO i = 1, klon
2579	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
2580	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
2581	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
2582	ENDDO
2583	ENDDO
2584	c
2585	ENDIF ! fin de test sur ok_orolf
2586	c
2587	IF (if_ebil.ge.2) THEN
2588	ztit='after orography'
2589	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2590	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2591	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2592	END IF
2593	c
2594	c
2595	cAA
2596	cAA Installation de l'interface online-offline pour traceurs
2597	cAA
2598	c====================================================================
2599	c Calcul des tendances traceurs
2600	c====================================================================
2601	C Pascale : il faut quand meme apeller phytrac car il gere les sorties
2602	cKE43 des traceurs => il faut donc mettre des flags a .false.
2603	IF (iflag_con.GE.3) THEN
2604	c on ajoute les tendances calculees par KE43
2605	cXXX OM on onhibe la convection sur les traceurs
2606	DO iq=1, nqmax-2 ! Sandrine a -3 ???
2607	cXXX OM on inhibe la convection sur les traceur
2608	cXXX DO k = 1, nlev
2609	cXXX DO i = 1, klon
2610	cXXX tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq)
2611	cXXX ENDDO
2612	cXXX ENDDO
2613	WRITE(iqn,'(i2.2)') iq
2614	CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
2615	ENDDO
2616	CMAF modif pour garder info du nombre de traceurs auxquels
2617	C la physique s'applique
2618	ELSE
2619	CMAF modif pour garder info du nombre de traceurs auxquels
2620	C la physique s'applique
2621	C
2622	call phytrac (rnpb,
2623	I debut,lafin,
2624	I nqmax-2,
2625	I nlon,nlev,dtime,
2626	I t,paprs,pplay,
2627	I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2628	I ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
2629	I frac_impa, frac_nucl,
2630	I rlon,presnivs,paire,pphis,
2631	O tr_seri)
2632	ENDIF
2633
2634	IF (offline) THEN
2635
2636	call phystokenc (
2637	I nlon,nlev,pdtphys,rlon,rlat,
2638	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2639	I ycoefh,yu1,yv1,ftsol,pctsrf,
2640	I frac_impa, frac_nucl,
2641	I pphis,paire,dtime,itap)
2642
2643
2644	ENDIF
2645
2646	c
2647	c Calculer le transport de l'eau et de l'energie (diagnostique)
2648	c
2649	CALL transp (paprs,zxtsol,
2650	e t_seri, q_seri, u_seri, v_seri, zphi,
2651	s ve, vq, ue, uq)
2652	c
2653	c Accumuler les variables a stocker dans les fichiers histoire:
2654	c
2655	c
2656	c
2657	c+jld ec_conser
2658	DO k = 1, klev
2659	DO i = 1, klon
2660	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
2661	d_t_ec(i,k)=0.5/ZRCPD
2662	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
2663	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
2664	d_t_ec(i,k) = d_t_ec(i,k)/dtime
2665	END DO
2666	END DO
2667	c-jld ec_conser
2668	IF (if_ebil.ge.1) THEN
2669	ztit='after physic'
2670	CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
2671	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2672	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2673	C Comme les tendances de la physique sont ajoute dans la dynamique,
2674	C on devrait avoir que la variation d'entalpie par la dynamique
2675	C est egale a la variation de la physique au pas de temps precedent.
2676	C Donc la somme de ces 2 variations devrait etre nulle.
2677	call diagphy(paire,ztit,ip_ebil
2678	e , topsw, toplw, solsw, sollw, sens
2679	e , evap, rain_fall, snow_fall, ztsol
2680	e , d_h_vcol, d_qt, d_ec
2681	s , fs_bound, fq_bound )
2682	C
2683	d_h_vcol_phy=d_h_vcol
2684	C
2685	END IF
2686	C
2687	cccIM cf. FH
2688	c=======================================================================
2689	c SORTIES
2690	c=======================================================================
2691
2692	c Interpollation sur quelques niveaux de pression
2693	c -----------------------------------------------
2694
2695	call plevel(klon,klev,.true. ,pplay,85000.,u_seri,u850)
2696	call plevel(klon,klev,.false.,pplay,85000.,v_seri,v850)
2697	call plevel(klon,klev,.true. ,pplay,50000.,u_seri,u500)
2698	call plevel(klon,klev,.false.,pplay,50000.,v_seri,v500)
2699	call plevel(klon,klev,.true. ,pplay,20000.,u_seri,u200)
2700	call plevel(klon,klev,.false.,pplay,20000.,v_seri,v200)
2701	call plevel(klon,klev,.true. ,pplay,50000.,zphi,phi500)
2702	call plevel(klon,klev,.true. ,paprs,50000.,omega,w500)
2703
2704	cIM cf. FH slp(:) = paprs(:,1)exp(pphis(:)/(289.t_seri(:,1)))
2705	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
2706	c PRINT*,' physiq slp ',slp(2185),paprs(2185,1),pphis(2185),
2707	c . RD,t_seri(2185,1)
2708	c
2709	ccc prw = eau precipitable
2710	DO i = 1, klon
2711	prw(i) = 0.
2712	DO k = 1, klev
2713	prw(i) = prw(i) +
2714	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
2715	ENDDO
2716	c PRINT*,' i ',i,' prw',prw(i)
2717	ENDDO
2718	c
2719
2720	c=============================================================
2721	c Ecriture des sorties
2722	c=============================================================
2723
2724	#ifdef histISCCP
2725	#include "write_histISCCP.h"
2726	#endif
2727
2728	#ifdef histhf
2729	#include "write_histhf.h"
2730	#endif
2731
2732	#include "write_histday.h"
2733	#include "write_histmth.h"
2734	#include "write_histins.h"
2735
2736	c=============================================================
2737	c
2738	c Convertir les incrementations en tendances
2739	c
2740	DO k = 1, klev
2741	DO i = 1, klon
2742	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
2743	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
2744	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
2745	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
2746	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
2747	ENDDO
2748	ENDDO
2749	c
2750	IF (nqmax.GE.3) THEN
2751	DO iq = 3, nqmax
2752	DO k = 1, klev
2753	DO i = 1, klon
2754	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
2755	ENDDO
2756	ENDDO
2757	ENDDO
2758	ENDIF
2759	c
2760	c Sauvegarder les valeurs de t et q a la fin de la physique:
2761	c
2762	DO k = 1, klev
2763	DO i = 1, klon
2764	t_ancien(i,k) = t_seri(i,k)
2765	q_ancien(i,k) = q_seri(i,k)
2766	ENDDO
2767	ENDDO
2768	c
2769	c====================================================================
2770	c Si c'est la fin, il faut conserver l'etat de redemarrage
2771	c====================================================================
2772	c
2773	IF (lafin) THEN
2774	itau_phy = itau_phy + itap
2775	ccc IF (ok_oasis) CALL quitcpl
2776	CALL phyredem ("restartphy.nc",dtime,radpas,
2777	. rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsurf, qsol,
2778	. fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
2779	. solsw, sollwdown,dlw,
2780	. radsol,frugs,agesno,
2781	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
2782	. t_ancien, q_ancien, rnebcon, ratqs, clwcon)
2783	ENDIF
2784
2785	RETURN
2786	END
2787	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
2788	IMPLICIT none
2789	c
2790	c Calculer et imprimer l'eau totale. A utiliser pour verifier
2791	c la conservation de l'eau
2792	c
2793	#include "YOMCST.h"
2794	INTEGER klon,klev
2795	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
2796	REAL aire(klon)
2797	REAL qtotal, zx, qcheck
2798	INTEGER i, k
2799	c
2800	zx = 0.0
2801	DO i = 1, klon
2802	zx = zx + aire(i)
2803	ENDDO
2804	qtotal = 0.0
2805	DO k = 1, klev
2806	DO i = 1, klon
2807	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
2808	. *(paprs(i,k)-paprs(i,k+1))/RG
2809	ENDDO
2810	ENDDO
2811	c
2812	qcheck = qtotal/zx
2813	c
2814	RETURN
2815	END
2816	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
2817	IMPLICIT none
2818	c
2819	c Tranformer une variable de la grille physique a
2820	c la grille d'ecriture
2821	c
2822	INTEGER nfield,nlon,iim,jjmp1, jjm
2823	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
2824	c
2825	INTEGER i, n, ig
2826	c
2827	jjm = jjmp1 - 1
2828	DO n = 1, nfield
2829	DO i=1,iim
2830	ecrit(i,n) = fi(1,n)
2831	ecrit(i+jjm*iim,n) = fi(nlon,n)
2832	ENDDO
2833	DO ig = 1, nlon - 2
2834	ecrit(iim+ig,n) = fi(1+ig,n)
2835	ENDDO
2836	ENDDO
2837	RETURN
2838	END
2839

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

Download in other formats:

Original Format