Context Navigation

physiq.F @ 476

Last change on this file since 476 was 476, checked in by lmdzadmin, 21 years ago
Erreur sur le passage du taux de CO2 + remise à zero de la date d'origine du calendrier de la simulation 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.2 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	c
908	C on remet le calendrier à zero
909	c
910	IF (raz_date .eq. 1) THEN
911	itau_phy = 0
912	ENDIF
913
914	c
915	CALL printflag( tabcntr0,radpas,ok_ocean,ok_oasis ,ok_journe,
916	, ok_instan, ok_region )
917	c
918	IF (ABS(dtime-pdtphys).GT.0.001) THEN
919	PRINT*, 'Pas physique n est pas correcte',dtime,pdtphys
920	abort_message=' See above '
921	call abort_gcm(modname,abort_message,1)
922	ENDIF
923	IF (nlon .NE. klon) THEN
924	PRINT*, 'nlon et klon ne sont pas coherents', nlon, klon
925	abort_message=' See above '
926	call abort_gcm(modname,abort_message,1)
927	ENDIF
928	IF (nlev .NE. klev) THEN
929	PRINT*, 'nlev et klev ne sont pas coherents', nlev, klev
930	abort_message=' See above '
931	call abort_gcm(modname,abort_message,1)
932	ENDIF
933	c
934	IF (dtime*FLOAT(radpas).GT.21600..AND.cycle_diurne) THEN
935	PRINT*, 'Nbre d appels au rayonnement insuffisant'
936	PRINT*, "Au minimum 4 appels par jour si cycle diurne"
937	abort_message=' See above '
938	call abort_gcm(modname,abort_message,1)
939	ENDIF
940	PRINT*, "Clef pour la convection, iflag_con=", iflag_con
941	PRINT*, "Clef pour le driver de la convection, ok_cvl=", ok_cvl
942	c
943	cKE43
944	c Initialisation pour la convection de K.E. (sb):
945	IF (iflag_con.GE.3) THEN
946
947	PRINT, "** Convection de Kerry Emanuel 4.3 "
948	PRINT*, "On va utiliser le melange convectif des traceurs qui"
949	PRINT*, "est calcule dans convect4.3"
950	PRINT*, " !!! penser aux logical flags de phytrac"
951
952	DO i = 1, klon
953	ema_cbmf(i) = 0.
954	ema_pcb(i) = 0.
955	ema_pct(i) = 0.
956	ema_workcbmf(i) = 0.
957	ENDDO
958
959	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>BEG
960	DO i = 1, klon
961	ibas_con(i) = 1
962	itop_con(i) = klev+1
963	ENDDO
964	cIM15/11/02 rajout initialisation ibas_con,itop_con cf. SB =>END
965
966	ENDIF
967
968	c34EK
969	IF (ok_orodr) THEN
970	DO i=1,klon
971	rugoro(i) = MAX(1.0e-05, zstd(i)*zsig(i)/2.0)
972	ENDDO
973	CALL SUGWD(klon,klev,paprs,pplay)
974	DO i=1,klon
975	zuthe(i)=0.
976	zvthe(i)=0.
977	if(zstd(i).gt.10.)then
978	zuthe(i)=(1.-zgam(i))*cos(zthe(i))
979	zvthe(i)=(1.-zgam(i))*sin(zthe(i))
980	endif
981	ENDDO
982	ENDIF
983	c
984	c
985	lmt_pas = NINT(86400./dtime * 1.0) ! tous les jours
986	PRINT*,'La frequence de lecture surface est de ', lmt_pas
987	c
988	ecrit_mth = NINT(86400./dtime *ecritphy) ! tous les ecritphy jours
989	IF (ok_mensuel) THEN
990	PRINT*, 'La frequence de sortie mensuelle est de ', ecrit_mth
991	ENDIF
992	ecrit_day = NINT(86400./dtime *1.0) ! tous les jours
993	IF (ok_journe) THEN
994	PRINT*, 'La frequence de sortie journaliere est de ',ecrit_day
995	ENDIF
996	ccc ecrit_ins = NINT(86400./dtime *0.5) ! 2 fois par jour
997	ccc ecrit_ins = NINT(86400./dtime *0.25) ! 4 fois par jour
998	ecrit_ins = NINT(86400./dtime/48.) ! a chaque pas de temps ==> PB. dans time_counter pour 1mois
999	ecrit_ins = NINT(86400./dtime/12.) ! toutes les deux heures
1000	IF (ok_instan) THEN
1001	PRINT*, 'La frequence de sortie instant. est de ', ecrit_ins
1002	ENDIF
1003	ecrit_reg = NINT(86400./dtime *0.25) ! 4 fois par jour
1004	IF (ok_region) THEN
1005	PRINT*, 'La frequence de sortie region est de ', ecrit_reg
1006	ENDIF
1007
1008	c
1009	c Initialiser le couplage si necessaire
1010	c
1011	npas = 0
1012	nexca = 0
1013	if (ocean == 'couple') then
1014	npas = itaufin/ iphysiq
1015	nexca = 86400 / dtime
1016	write(,)' ##### Ocean couple #####'
1017	write(,)' Valeurs des pas de temps'
1018	write(,)' npas = ', npas
1019	write(,)' nexca = ', nexca
1020	endif
1021	c
1022	c
1023	cccIM
1024	capemaxcels = 't_max(X)'
1025	t2mincels = 't_min(X)'
1026	t2maxcels = 't_max(X)'
1027
1028	cccIM cf. FH
1029	c
1030	c=============================================================
1031	c Initialisation des sorties
1032	c=============================================================
1033
1034	#ifdef histISCCP
1035	#include "ini_histISCCP.h"
1036	#endif
1037
1038	#ifdef histhf
1039	#include "ini_histhf.h"
1040	#endif
1041
1042	#include "ini_histday.h"
1043	#include "ini_histmth.h"
1044	#include "ini_histins.h"
1045
1046	cXXXPB Positionner date0 pour initialisation de ORCHIDEE
1047	date0 = zjulian
1048	C date0 = day_ini
1049	WRITE(,) 'physiq date0 : ',date0
1050	c
1051	c
1052	c
1053	c Prescrire l'ozone dans l'atmosphere
1054	c
1055	c
1056	cc DO i = 1, klon
1057	cc DO k = 1, klev
1058	cc CALL o3cm (paprs(i,k)/100.,paprs(i,k+1)/100., wo(i,k),20)
1059	cc ENDDO
1060	cc ENDDO
1061	c
1062	c
1063	ENDIF
1064	c
1065	c ************** Fin de IF ( debut ) *************
1066	c
1067	c
1068	c Mettre a zero des variables de sortie (pour securite)
1069	c
1070	DO i = 1, klon
1071	d_ps(i) = 0.0
1072	ENDDO
1073	DO k = 1, klev
1074	DO i = 1, klon
1075	d_t(i,k) = 0.0
1076	d_u(i,k) = 0.0
1077	d_v(i,k) = 0.0
1078	ENDDO
1079	ENDDO
1080	DO iq = 1, nqmax
1081	DO k = 1, klev
1082	DO i = 1, klon
1083	d_qx(i,k,iq) = 0.0
1084	ENDDO
1085	ENDDO
1086	ENDDO
1087	c
1088	c Ne pas affecter les valeurs entrees de u, v, h, et q
1089	c
1090	DO k = 1, klev
1091	DO i = 1, klon
1092	t_seri(i,k) = t(i,k)
1093	u_seri(i,k) = u(i,k)
1094	v_seri(i,k) = v(i,k)
1095	q_seri(i,k) = qx(i,k,ivap)
1096	ql_seri(i,k) = qx(i,k,iliq)
1097	qs_seri(i,k) = 0.
1098	ENDDO
1099	ENDDO
1100	IF (nqmax.GE.3) THEN
1101	DO iq = 3, nqmax
1102	DO k = 1, klev
1103	DO i = 1, klon
1104	tr_seri(i,k,iq-2) = qx(i,k,iq)
1105	ENDDO
1106	ENDDO
1107	ENDDO
1108	ELSE
1109	DO k = 1, klev
1110	DO i = 1, klon
1111	tr_seri(i,k,1) = 0.0
1112	ENDDO
1113	ENDDO
1114	ENDIF
1115	C
1116	DO i = 1, klon
1117	ztsol(i) = 0.
1118	ENDDO
1119	DO nsrf = 1, nbsrf
1120	DO i = 1, klon
1121	ztsol(i) = ztsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1122	ENDDO
1123	ENDDO
1124	C
1125	IF (if_ebil.ge.1) THEN
1126	ztit='after dynamic'
1127	CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
1128	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1129	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1130	C Comme les tendances de la physique sont ajoute dans la dynamique,
1131	C on devrait avoir que la variation d'entalpie par la dynamique
1132	C est egale a la variation de la physique au pas de temps precedent.
1133	C Donc la somme de ces 2 variations devrait etre nulle.
1134	call diagphy(paire,ztit,ip_ebil
1135	e , zero_v, zero_v, zero_v, zero_v, zero_v
1136	e , zero_v, zero_v, zero_v, ztsol
1137	e , d_h_vcol+d_h_vcol_phy, d_qt, 0.
1138	s , fs_bound, fq_bound )
1139	END IF
1140
1141	c Diagnostiquer la tendance dynamique
1142	c
1143	IF (ancien_ok) THEN
1144	DO k = 1, klev
1145	DO i = 1, klon
1146	d_t_dyn(i,k) = (t_seri(i,k)-t_ancien(i,k))/dtime
1147	d_q_dyn(i,k) = (q_seri(i,k)-q_ancien(i,k))/dtime
1148	ENDDO
1149	ENDDO
1150	ELSE
1151	DO k = 1, klev
1152	DO i = 1, klon
1153	d_t_dyn(i,k) = 0.0
1154	d_q_dyn(i,k) = 0.0
1155	ENDDO
1156	ENDDO
1157	ancien_ok = .TRUE.
1158	ENDIF
1159	c
1160	c Ajouter le geopotentiel du sol:
1161	c
1162	DO k = 1, klev
1163	DO i = 1, klon
1164	zphi(i,k) = pphi(i,k) + pphis(i)
1165	ENDDO
1166	ENDDO
1167	c
1168	c Verifier les temperatures
1169	c
1170	CALL hgardfou(t_seri,ftsol,'debutphy')
1171	c
1172	c Incrementer le compteur de la physique
1173	c
1174	itap = itap + 1
1175	julien = MOD(NINT(xjour),360)
1176	c
1177	c Mettre en action les conditions aux limites (albedo, sst, etc.).
1178	c Prescrire l'ozone et calculer l'albedo sur l'ocean.
1179	c
1180	IF (MOD(itap-1,lmt_pas) .EQ. 0) THEN
1181	PRINT *,' PHYS cond julien ',julien
1182	CALL ozonecm( FLOAT(julien), rlat, paprs, wo)
1183	ENDIF
1184	c
1185	c Re-evaporer l'eau liquide nuageuse
1186	c
1187	DO k = 1, klev ! re-evaporation de l'eau liquide nuageuse
1188	DO i = 1, klon
1189	zlvdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1190	c zlsdcp=RLSTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1191	zlsdcp=RLVTT/RCPD/(1.0+RVTMP2*q_seri(i,k))
1192	zdelta = MAX(0.,SIGN(1.,RTT-t_seri(i,k)))
1193	zb = MAX(0.0,ql_seri(i,k))
1194	za = - MAX(0.0,ql_seri(i,k))
1195	. * (zlvdcp(1.-zdelta)+zlsdcpzdelta)
1196	t_seri(i,k) = t_seri(i,k) + za
1197	q_seri(i,k) = q_seri(i,k) + zb
1198	ql_seri(i,k) = 0.0
1199	d_t_eva(i,k) = za
1200	d_q_eva(i,k) = zb
1201	ENDDO
1202	ENDDO
1203	c
1204	IF (if_ebil.ge.2) THEN
1205	ztit='after reevap'
1206	CALL diagetpq(paire,ztit,ip_ebil,2,1,dtime
1207	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1208	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1209	call diagphy(paire,ztit,ip_ebil
1210	e , zero_v, zero_v, zero_v, zero_v, zero_v
1211	e , zero_v, zero_v, zero_v, ztsol
1212	e , d_h_vcol, d_qt, d_ec
1213	s , fs_bound, fq_bound )
1214	C
1215	END IF
1216	C
1217	c
1218	c Appeler la diffusion verticale (programme de couche limite)
1219	c
1220	DO i = 1, klon
1221	c if (.not. ok_veget) then
1222	c frugs(i,is_ter) = SQRT(frugs(i,is_ter)2+rugoro(i)2)
1223	c endif
1224	c frugs(i,is_lic) = rugoro(i)
1225	c frugs(i,is_oce) = rugmer(i)
1226	c frugs(i,is_sic) = 0.001
1227	zxrugs(i) = 0.0
1228	ENDDO
1229	DO nsrf = 1, nbsrf
1230	DO i = 1, klon
1231	c frugs(i,nsrf) = MAX(frugs(i,nsrf),0.001)
1232	frugs(i,nsrf) = MAX(frugs(i,nsrf),0.000015)
1233	ENDDO
1234	ENDDO
1235	DO nsrf = 1, nbsrf
1236	DO i = 1, klon
1237	zxrugs(i) = zxrugs(i) + frugs(i,nsrf)*pctsrf(i,nsrf)
1238	ENDDO
1239	ENDDO
1240	c
1241	C calculs necessaires au calcul de l'albedo dans l'interface
1242	c
1243	CALL orbite(FLOAT(julien),zlongi,dist)
1244	IF (cycle_diurne) THEN
1245	zdtime=dtime*FLOAT(radpas) ! pas de temps du rayonnement (s)
1246	CALL zenang(zlongi,gmtime,zdtime,rlat,rlon,rmu0,fract)
1247	ELSE
1248	rmu0 = -999.999
1249	ENDIF
1250	cIM BEG
1251	DO i=1, klon
1252	sunlit(i)=1
1253	IF(rmu0(i).EQ.0.) sunlit(i)=0
1254	c IF(rmu0(i).EQ.0.) THEN
1255	c sunlit(i)=0
1256	c PRINT*,' il fait nuit ',i,rlat(i),rlon(i)
1257	c ENDIF
1258	ENDDO
1259	cIM END
1260	C Calcul de l'abedo moyen par maille
1261	albsol(:)=0.
1262	albsollw(:)=0.
1263	DO nsrf = 1, nbsrf
1264	DO i = 1, klon
1265	albsol(i) = albsol(i) + falbe(i,nsrf) * pctsrf(i,nsrf)
1266	albsollw(i) = albsollw(i) + falblw(i,nsrf) * pctsrf(i,nsrf)
1267	ENDDO
1268	ENDDO
1269	C
1270	C Repartition sous maille des flux LW et SW
1271	C Modif OM+PASB+JLD
1272	C Repartition du longwave par sous-surface linearisee
1273	Cn
1274	DO nsrf = 1, nbsrf
1275	DO i = 1, klon
1276	c$$$ fsollw(i,nsrf) = sollwdown(i) - RSIGMAftsol(i,nsrf)*4
1277	c$$$ fsollw(i,nsrf) = sollw(i)
1278	fsollw(i,nsrf) = sollw(i)
1279	$ + 4.0RSIGMAztsol(i)*3 (ztsol(i)-ftsol(i,nsrf))
1280	fsolsw(i,nsrf) = solsw(i)*(1.-falbe(i,nsrf))/(1.-albsol(i))
1281	ENDDO
1282	ENDDO
1283
1284	fder = dlw
1285
1286	CALL clmain(dtime,itap,date0,pctsrf,
1287	e t_seri,q_seri,u_seri,v_seri,
1288	e julien, rmu0, co2_ppm,
1289	e ok_veget, ocean, npas, nexca, ftsol,
1290	$ soil_model,cdmmax, cdhmax, ftsoil, qsol,
1291	$ paprs,pplay,radsol, fsnow,fqsurf,fevap,falbe,falblw,
1292	$ fluxlat,
1293	cIM cf. JLD e rain_fall, snow_fall, solsw, sollw, sollwdown, fder,
1294	e rain_fall, snow_fall, fsolsw, fsollw, sollwdown, fder,
1295	e rlon, rlat, cufi, cvfi, frugs,
1296	e debut, lafin, agesno,rugoro ,
1297	s d_t_vdf,d_q_vdf,d_u_vdf,d_v_vdf,d_ts,
1298	s fluxt,fluxq,fluxu,fluxv,cdragh,cdragm,
1299	s dsens, devap,
1300	cIM cf JLD s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m)
1301	s ycoefh,yu1,yv1, t2m, q2m, u10m, v10m,
1302	s fqcalving,ffonte)
1303	c
1304	CXXX PB
1305	CXXX Incrementation des flux
1306	CXXX
1307	zxfluxt=0.
1308	zxfluxq=0.
1309	zxfluxu=0.
1310	zxfluxv=0.
1311	DO nsrf = 1, nbsrf
1312	DO k = 1, klev
1313	DO i = 1, klon
1314	zxfluxt(i,k) = zxfluxt(i,k) +
1315	$ fluxt(i,k,nsrf) * pctsrf( i, nsrf)
1316	zxfluxq(i,k) = zxfluxq(i,k) +
1317	$ fluxq(i,k,nsrf) * pctsrf( i, nsrf)
1318	zxfluxu(i,k) = zxfluxu(i,k) +
1319	$ fluxu(i,k,nsrf) * pctsrf( i, nsrf)
1320	zxfluxv(i,k) = zxfluxv(i,k) +
1321	$ fluxv(i,k,nsrf) * pctsrf( i, nsrf)
1322	END DO
1323	END DO
1324	END DO
1325	DO i = 1, klon
1326	sens(i) = - zxfluxt(i,1) ! flux de chaleur sensible au sol
1327	c evap(i) = - fluxq(i,1) ! flux d'evaporation au sol
1328	evap(i) = - zxfluxq(i,1) ! flux d'evaporation au sol
1329	fder(i) = dlw(i) + dsens(i) + devap(i)
1330	ENDDO
1331
1332
1333	DO k = 1, klev
1334	DO i = 1, klon
1335	t_seri(i,k) = t_seri(i,k) + d_t_vdf(i,k)
1336	q_seri(i,k) = q_seri(i,k) + d_q_vdf(i,k)
1337	u_seri(i,k) = u_seri(i,k) + d_u_vdf(i,k)
1338	v_seri(i,k) = v_seri(i,k) + d_v_vdf(i,k)
1339	ENDDO
1340	ENDDO
1341	c
1342	IF (if_ebil.ge.2) THEN
1343	ztit='after clmain'
1344	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1345	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1346	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1347	call diagphy(paire,ztit,ip_ebil
1348	e , zero_v, zero_v, zero_v, zero_v, sens
1349	e , evap , zero_v, zero_v, ztsol
1350	e , d_h_vcol, d_qt, d_ec
1351	s , fs_bound, fq_bound )
1352	END IF
1353	C
1354	c
1355	c Incrementer la temperature du sol
1356	c
1357	DO i = 1, klon
1358	zxtsol(i) = 0.0
1359	zxfluxlat(i) = 0.0
1360	cccIM
1361	zt2m(i) = 0.0
1362	zq2m(i) = 0.0
1363	zu10m(i) = 0.0
1364	zv10m(i) = 0.0
1365	cIM cf JLD ??
1366	zxffonte(i) = 0.0
1367	zxfqcalving(i) = 0.0
1368	c
1369	IF ( abs( pctsrf(i, is_ter) + pctsrf(i, is_lic) +
1370	$ pctsrf(i, is_oce) + pctsrf(i, is_sic) - 1.) .GT. EPSFRA)
1371	$ THEN
1372	WRITE(,) 'physiq : pb sous surface au point ', i,
1373	$ pctsrf(i, 1 : nbsrf)
1374	ENDIF
1375	ENDDO
1376	DO nsrf = 1, nbsrf
1377	DO i = 1, klon
1378	c IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1379	ftsol(i,nsrf) = ftsol(i,nsrf) + d_ts(i,nsrf)
1380	cIM cf. JLD
1381	wfbils(i,nsrf) = ( fsolsw(i,nsrf) + fsollw(i,nsrf)
1382	$ + fluxt(i,1,nsrf) + fluxlat(i,nsrf) ) * pctsrf(i,nsrf)
1383	zxtsol(i) = zxtsol(i) + ftsol(i,nsrf)*pctsrf(i,nsrf)
1384	zxfluxlat(i) = zxfluxlat(i) + fluxlat(i,nsrf)*pctsrf(i,nsrf)
1385	cccIM
1386	zt2m(i) = zt2m(i) + t2m(i,nsrf)*pctsrf(i,nsrf)
1387	zq2m(i) = zq2m(i) + q2m(i,nsrf)*pctsrf(i,nsrf)
1388	zu10m(i) = zu10m(i) + u10m(i,nsrf)*pctsrf(i,nsrf)
1389	zv10m(i) = zv10m(i) + v10m(i,nsrf)*pctsrf(i,nsrf)
1390	cIM cf JLD ??
1391	zxffonte(i) = zxffonte(i) + ffonte(i,nsrf)*pctsrf(i,nsrf)
1392	zxfqcalving(i) = zxfqcalving(i) +
1393	. fqcalving(i,nsrf)*pctsrf(i,nsrf)
1394	c ENDIF
1395	ENDDO
1396	ENDDO
1397
1398	c
1399	c Si une sous-fraction n'existe pas, elle prend la temp. moyenne
1400	c
1401	DO nsrf = 1, nbsrf
1402	DO i = 1, klon
1403	IF (pctsrf(i,nsrf) .LT. epsfra) ftsol(i,nsrf) = zxtsol(i)
1404	cccIM
1405	IF (pctsrf(i,nsrf) .LT. epsfra) t2m(i,nsrf) = zt2m(i)
1406	IF (pctsrf(i,nsrf) .LT. epsfra) q2m(i,nsrf) = zq2m(i)
1407	IF (pctsrf(i,nsrf) .LT. epsfra) u10m(i,nsrf) = zu10m(i)
1408	IF (pctsrf(i,nsrf) .LT. epsfra) v10m(i,nsrf) = zv10m(i)
1409	cIM cf JLD ??
1410	IF (pctsrf(i,nsrf) .LT. epsfra) ffonte(i,nsrf) = zxffonte(i)
1411	IF (pctsrf(i,nsrf) .LT. epsfra)
1412	. fqcalving(i,nsrf) = zxfqcalving(i)
1413	ENDDO
1414	ENDDO
1415	c
1416	c
1417	c Calculer la derive du flux infrarouge
1418	c
1419	cXXX DO nsrf = 1, nbsrf
1420	DO i = 1, klon
1421	cXXX IF (pctsrf(i,nsrf) .GE. EPSFRA) THEN
1422	dlw(i) = - 4.0RSIGMAzxtsol(i)**3
1423	cXXX . *(ftsol(i,nsrf)-zxtsol(i))
1424	cXXX . *pctsrf(i,nsrf)
1425	cXXX ENDIF
1426	cXXX ENDDO
1427	ENDDO
1428	c
1429	c Appeler la convection (au choix)
1430	c
1431	DO k = 1, klev
1432	DO i = 1, klon
1433	conv_q(i,k) = d_q_dyn(i,k)
1434	. + d_q_vdf(i,k)/dtime
1435	conv_t(i,k) = d_t_dyn(i,k)
1436	. + d_t_vdf(i,k)/dtime
1437	ENDDO
1438	ENDDO
1439	IF (check) THEN
1440	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1441	PRINT*, "avantcon=", za
1442	ENDIF
1443	zx_ajustq = .FALSE.
1444	IF (iflag_con.EQ.2) zx_ajustq=.TRUE.
1445	IF (zx_ajustq) THEN
1446	DO i = 1, klon
1447	z_avant(i) = 0.0
1448	ENDDO
1449	DO k = 1, klev
1450	DO i = 1, klon
1451	z_avant(i) = z_avant(i) + (q_seri(i,k)+ql_seri(i,k))
1452	. *(paprs(i,k)-paprs(i,k+1))/RG
1453	ENDDO
1454	ENDDO
1455	ENDIF
1456	IF (iflag_con.EQ.1) THEN
1457	stop'reactiver le call conlmd dans physiq.F'
1458	c CALL conlmd (dtime, paprs, pplay, t_seri, q_seri, conv_q,
1459	c . d_t_con, d_q_con,
1460	c . rain_con, snow_con, ibas_con, itop_con)
1461	ELSE IF (iflag_con.EQ.2) THEN
1462	CALL conflx(dtime, paprs, pplay, t_seri, q_seri,
1463	e conv_t, conv_q, zxfluxq(1,1), omega,
1464	s d_t_con, d_q_con, rain_con, snow_con,
1465	s pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
1466	s kcbot, kctop, kdtop, pmflxr, pmflxs)
1467	WHERE (rain_con < 0.) rain_con = 0.
1468	WHERE (snow_con < 0.) snow_con = 0.
1469	DO i = 1, klon
1470	ibas_con(i) = klev+1 - kcbot(i)
1471	itop_con(i) = klev+1 - kctop(i)
1472	ENDDO
1473	ELSE IF (iflag_con.GE.3) THEN
1474	c nb of tracers for the KE convection:
1475	if (nqmax .GE. 4) then
1476	ntra = nbtr
1477	else
1478	ntra = 1
1479	endif
1480	c
1481	c sb, oct02:
1482	c Schema de convection modularise et vectorise:
1483	c (driver commun aux versions 3 et 4)
1484	c
1485	IF (ok_cvl) THEN ! new driver for convectL
1486
1487	CALL concvl (iflag_con,
1488	. dtime,paprs,pplay,t_seri,q_seri,
1489	. u_seri,v_seri,tr_seri,nbtr,
1490	. ema_work1,ema_work2,
1491	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1492	. rain_con, snow_con, ibas_con, itop_con,
1493	. upwd,dnwd,dnwd0,
1494	. Ma,cape,tvp,iflagctrl,
1495	. pbase,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr,qcondc,wd)
1496	cIM cf. FH
1497	clwcon0=qcondc
1498
1499	ELSE ! ok_cvl
1500
1501	c print*,'Avant conema OUI'
1502	CALL conema3 (dtime,
1503	. paprs,pplay,t_seri,q_seri,
1504	. u_seri,v_seri,tr_seri,nbtr,
1505	. ema_work1,ema_work2,
1506	. d_t_con,d_q_con,d_u_con,d_v_con,d_tr,
1507	. rain_con, snow_con, ibas_con, itop_con,
1508	. upwd,dnwd,dnwd0,bas,top,
1509	. Ma,cape,tvp,rflag,
1510	. pbase
1511	. ,bbase,dtvpdt1,dtvpdq1,dplcldt,dplcldr
1512	. ,clwcon0)
1513	print*,'Apres conema3 '
1514
1515	ENDIF ! ok_cvl
1516
1517	IF (.NOT. ok_gust) THEN
1518	do i = 1, klon
1519	wd(i)=0.0
1520	enddo
1521	ENDIF
1522
1523	c =================================================================== c
1524	c Calcul des proprietes des nuages convectifs
1525	c
1526	DO k = 1, klev
1527	DO i = 1, klon
1528	zx_t = t_seri(i,k)
1529	IF (thermcep) THEN
1530	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
1531	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
1532	zx_qs = MIN(0.5,zx_qs)
1533	zcor = 1./(1.-retv*zx_qs)
1534	zx_qs = zx_qs*zcor
1535	ELSE
1536	IF (zx_t.LT.t_coup) THEN
1537	zx_qs = qsats(zx_t)/pplay(i,k)
1538	ELSE
1539	zx_qs = qsatl(zx_t)/pplay(i,k)
1540	ENDIF
1541	ENDIF
1542	zqsat(i,k)=zx_qs
1543	ENDDO
1544	ENDDO
1545
1546	c calcul des proprietes des nuages convectifs
1547	clwcon0(:,:)=fact_cldcon*clwcon0(:,:)
1548	call clouds_gno
1549	s (klon,klev,q_seri,zqsat,clwcon0,ptconv,ratqsc,rnebcon0)
1550
1551	c =================================================================== c
1552
1553	DO i = 1, klon
1554	ema_pcb(i) = pbase(i)
1555	ENDDO
1556	DO i = 1, klon
1557	ema_pct(i) = paprs(i,itop_con(i))
1558	ENDDO
1559	DO i = 1, klon
1560	ema_cbmf(i) = ema_workcbmf(i)
1561	ENDDO
1562	ELSE
1563	PRINT*, "iflag_con non-prevu", iflag_con
1564	CALL abort
1565	ENDIF
1566
1567	c CALL homogene(paprs, q_seri, d_q_con, u_seri,v_seri,
1568	c . d_u_con, d_v_con)
1569	DO k = 1, klev
1570	DO i = 1, klon
1571	t_seri(i,k) = t_seri(i,k) + d_t_con(i,k)
1572	q_seri(i,k) = q_seri(i,k) + d_q_con(i,k)
1573	u_seri(i,k) = u_seri(i,k) + d_u_con(i,k)
1574	v_seri(i,k) = v_seri(i,k) + d_v_con(i,k)
1575	ENDDO
1576	ENDDO
1577	c
1578	IF (if_ebil.ge.2) THEN
1579	ztit='after convect'
1580	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1581	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1582	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1583	call diagphy(paire,ztit,ip_ebil
1584	e , zero_v, zero_v, zero_v, zero_v, zero_v
1585	e , zero_v, rain_con, snow_con, ztsol
1586	e , d_h_vcol, d_qt, d_ec
1587	s , fs_bound, fq_bound )
1588	END IF
1589	C
1590	IF (check) THEN
1591	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1592	PRINT*, "aprescon=", za
1593	zx_t = 0.0
1594	za = 0.0
1595	DO i = 1, klon
1596	za = za + paire(i)/FLOAT(klon)
1597	zx_t = zx_t + (rain_con(i)+snow_con(i))*paire(i)/FLOAT(klon)
1598	ENDDO
1599	zx_t = zx_t/za*dtime
1600	PRINT*, "Precip=", zx_t
1601	ENDIF
1602	IF (zx_ajustq) THEN
1603	DO i = 1, klon
1604	z_apres(i) = 0.0
1605	ENDDO
1606	DO k = 1, klev
1607	DO i = 1, klon
1608	z_apres(i) = z_apres(i) + (q_seri(i,k)+ql_seri(i,k))
1609	. *(paprs(i,k)-paprs(i,k+1))/RG
1610	ENDDO
1611	ENDDO
1612	DO i = 1, klon
1613	z_factor(i) = (z_avant(i)-(rain_con(i)+snow_con(i))*dtime)
1614	. /z_apres(i)
1615	ENDDO
1616	DO k = 1, klev
1617	DO i = 1, klon
1618	IF (z_factor(i).GT.(1.0+1.0E-08) .OR.
1619	. z_factor(i).LT.(1.0-1.0E-08)) THEN
1620	q_seri(i,k) = q_seri(i,k) * z_factor(i)
1621	ENDIF
1622	ENDDO
1623	ENDDO
1624	ENDIF
1625	zx_ajustq=.FALSE.
1626	c
1627	IF (nqmax.GT.2) THEN !--melange convectif de traceurs
1628	c
1629	IF (iflag_con .NE. 2 .AND. debut) THEN
1630	PRINT*, 'Pour l instant, seul conflx fonctionne ',
1631	$ 'avec traceurs', iflag_con
1632	PRINT*,' Mettre iflag_con',
1633	$ ' = 2 dans run.def et repasser'
1634	c CALL abort
1635	ENDIF
1636	c
1637	ENDIF !--nqmax.GT.2
1638	c
1639	c Appeler l'ajustement sec
1640	c
1641	CALL ajsec(paprs, pplay, t_seri, q_seri, d_t_ajs, d_q_ajs)
1642	DO k = 1, klev
1643	DO i = 1, klon
1644	t_seri(i,k) = t_seri(i,k) + d_t_ajs(i,k)
1645	q_seri(i,k) = q_seri(i,k) + d_q_ajs(i,k)
1646	ENDDO
1647	ENDDO
1648	c
1649	IF (if_ebil.ge.2) THEN
1650	ztit='after dry_adjust'
1651	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1652	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1653	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1654	END IF
1655
1656
1657	c-------------------------------------------------------------------------
1658	c Caclul des ratqs
1659	c-------------------------------------------------------------------------
1660
1661	c print*,'calcul des ratqs'
1662	c ratqs convectifs a l'ancienne en fonction de q(z=0)-q / q
1663	c ----------------
1664	c on ecrase le tableau ratqsc calcule par clouds_gno
1665	if (iflag_cldcon.eq.1) then
1666	do k=1,klev
1667	do i=1,klon
1668	if(ptconv(i,k)) then
1669	ratqsc(i,k)=ratqsbas
1670	s +fact_cldcon*(q_seri(i,1)-q_seri(i,k))/q_seri(i,k)
1671	else
1672	ratqsc(i,k)=0.
1673	endif
1674	enddo
1675	enddo
1676	endif
1677
1678	c ratqs stables
1679	c -------------
1680	do k=1,klev
1681	ratqss(:,k)=ratqsbas+(ratqshaut-ratqsbas)*
1682	s min((paprs(:,1)-pplay(:,k))/(paprs(:,1)-30000.),1.)
1683	enddo
1684
1685
1686	c ratqs final
1687	c -----------
1688	if (iflag_cldcon.eq.1 .or.iflag_cldcon.eq.2) then
1689	c les ratqs sont une conbinaison de ratqss et ratqsc
1690	c ratqs final
1691	c 1e4 (en gros 3 heures), en dur pour le moment, est le temps de
1692	c relaxation des ratqs
1693	c facttemps=exp(-pdtphys/1.e4)
1694	facteur=exp(-pdtphys*facttemps)
1695	ratqs(:,:)=max(ratqs(:,:)*facteur,ratqss(:,:))
1696	ratqs(:,:)=max(ratqs(:,:),ratqsc(:,:))
1697	c print*,'calcul des ratqs fini'
1698	else
1699	c on ne prend que le ratqs stable pour fisrtilp
1700	ratqs(:,:)=ratqss(:,:)
1701	endif
1702
1703
1704	c
1705	c Appeler le processus de condensation a grande echelle
1706	c et le processus de precipitation
1707	c-------------------------------------------------------------------------
1708	CALL fisrtilp(dtime,paprs,pplay,
1709	. t_seri, q_seri,ptconv,ratqs,
1710	. d_t_lsc, d_q_lsc, d_ql_lsc, rneb, cldliq,
1711	. rain_lsc, snow_lsc,
1712	. pfrac_impa, pfrac_nucl, pfrac_1nucl,
1713	. frac_impa, frac_nucl,
1714	. prfl, psfl, rhcl)
1715
1716	WHERE (rain_lsc < 0) rain_lsc = 0.
1717	WHERE (snow_lsc < 0) snow_lsc = 0.
1718	DO k = 1, klev
1719	DO i = 1, klon
1720	t_seri(i,k) = t_seri(i,k) + d_t_lsc(i,k)
1721	q_seri(i,k) = q_seri(i,k) + d_q_lsc(i,k)
1722	ql_seri(i,k) = ql_seri(i,k) + d_ql_lsc(i,k)
1723	cldfra(i,k) = rneb(i,k)
1724	IF (.NOT.new_oliq) cldliq(i,k) = ql_seri(i,k)
1725	ENDDO
1726	ENDDO
1727	IF (check) THEN
1728	za = qcheck(klon,klev,paprs,q_seri,ql_seri,paire)
1729	PRINT*, "apresilp=", za
1730	zx_t = 0.0
1731	za = 0.0
1732	DO i = 1, klon
1733	za = za + paire(i)/FLOAT(klon)
1734	zx_t = zx_t + (rain_lsc(i)+snow_lsc(i))*paire(i)/FLOAT(klon)
1735	ENDDO
1736	zx_t = zx_t/za*dtime
1737	PRINT*, "Precip=", zx_t
1738	ENDIF
1739	c
1740	IF (if_ebil.ge.2) THEN
1741	ztit='after fisrt'
1742	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
1743	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
1744	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
1745	call diagphy(paire,ztit,ip_ebil
1746	e , zero_v, zero_v, zero_v, zero_v, zero_v
1747	e , zero_v, rain_lsc, snow_lsc, ztsol
1748	e , d_h_vcol, d_qt, d_ec
1749	s , fs_bound, fq_bound )
1750	END IF
1751	c
1752	c-------------------------------------------------------------------
1753	c PRESCRIPTION DES NUAGES POUR LE RAYONNEMENT
1754	c-------------------------------------------------------------------
1755
1756	c 1. NUAGES CONVECTIFS
1757	c
1758	IF (iflag_cldcon.eq.-1) THEN ! seulement pour Tiedtke
1759
1760	c Nuages diagnostiques pour Tiedtke
1761	CALL diagcld1(paprs,pplay,
1762	. rain_con,snow_con,ibas_con,itop_con,
1763	. diafra,dialiq)
1764	DO k = 1, klev
1765	DO i = 1, klon
1766	IF (diafra(i,k).GT.cldfra(i,k)) THEN
1767	cldliq(i,k) = dialiq(i,k)
1768	cldfra(i,k) = diafra(i,k)
1769	ENDIF
1770	ENDDO
1771	ENDDO
1772
1773	ELSE IF (iflag_cldcon.eq.3) THEN
1774	c On prend pour les nuages convectifs le max du calcul de la
1775	c convection et du calcul du pas de temps précédent diminué d'un facteur
1776	c facttemps
1777	c facttemps=pdtphys/1.e4
1778	facteur = pdtphys *facttemps
1779	do k=1,klev
1780	do i=1,klon
1781	rnebcon(i,k)=rnebcon(i,k)*facteur
1782	if (rnebcon0(i,k)clwcon0(i,k).gt.rnebcon(i,k)clwcon(i,k))
1783	s then
1784	rnebcon(i,k)=rnebcon0(i,k)
1785	clwcon(i,k)=clwcon0(i,k)
1786	endif
1787	enddo
1788	enddo
1789
1790	cIM ISCCP simulator BEGIN
1791	IF (ok_isccp) THEN
1792	cIM calcul tau. emi nuages convectifs
1793	convfra(:,:)=rnebcon(:,:)
1794	convliq(:,:)=rnebcon(:,:)*clwcon(:,:)
1795	c CALL newmicro (paprs, pplay,ok_newmicro,
1796	c . t_seri, cldliq, cldfra, cldtau, cldemi,
1797	c . cldh, cldl, cldm, cldt, cldq)
1798	CALL newmicro (paprs, pplay,ok_newmicro,
1799	. t_seri, convliq, convfra, dtau_c, dem_c,
1800	. cldh_c, cldl_c, cldm_c, cldt_c, cldq_c)
1801
1802	cIM calcul tau. emi nuages startiformes
1803	CALL newmicro (paprs, pplay,ok_newmicro,
1804	. t_seri, cldliq, cldfra, dtau_s, dem_s,
1805	. cldh_s, cldl_s, cldm_s, cldt_s, cldq_s)
1806	cIM calcul diagramme (PC, tau) cf. ISCCP D
1807	c seed=50
1808	c seed=ran0(klon)
1809	cT1O3
1810	c top_height=1
1811	cT3O3
1812	c top_height=3
1813	c overlap=3
1814	cIM cf GCM
1815	cldtot(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
1816
1817	cIM inversion des niveaux de pression ==> de haut en bas
1818	DO k=1,klev
1819	kinv=klev-k+1
1820	DO i=1,klon
1821	pfull(i,k)=pplay(i,kinv)
1822	c on met toutes les variables de Haut 2 Bas
1823	qv(i,k)=q_seri(i,kinv)
1824	cc(i,k)=cldtot(i,kinv)
1825	conv(i,k)=rnebcon(i,kinv)
1826	dtau_sH2B(i,k)=dtau_s(i,kinv)
1827	dtau_cH2B(i,k)=dtau_c(i,kinv)
1828	at(i,k)=t_seri(i,kinv)
1829	dem_sH2B(i,k)=dem_s(i,kinv)
1830	dem_cH2B(i,k)=dem_c(i,kinv)
1831
1832	ENDDO
1833	ENDDO
1834
1835	DO k=1,klev+1
1836	kinv=klev-k+2
1837	DO i=1,klon
1838	phalf(i,k)=paprs(i,kinv)
1839	ENDDO
1840	ENDDO
1841
1842	c open(99,file='tautab.bin',access='sequential',
1843	c $ form='unformatted',status='old')
1844	c read(99) tautab
1845
1846	cIM210503
1847	IF (debut) THEN
1848	open(99,file='tautab.formatted', FORM='FORMATTED')
1849	read(99,'(f30.20)') tautab
1850	close(99)
1851	c
1852	open(99,file='invtau.formatted',form='FORMATTED')
1853	read(99,'(i10)') invtau
1854	close(99)
1855	c
1856	nsrf=3
1857	DO nreg=1, nbreg
1858	DO i=1, klon
1859
1860	c IF (debut) THEN
1861	IF(rlon(i).LT.0.) THEN
1862	rlonPOS(i)=rlon(i)+360.
1863	ELSE
1864	rlonPOS(i)=rlon(i)
1865	ENDIF
1866	c ENDIF
1867
1868	c pct_ocean(i,nreg)=.FALSE.
1869	pct_ocean(i,nreg)=0
1870
1871	c DO nsrf = 1, nbsrf
1872
1873	c test si c'est 1 point d'ocean
1874	IF(pctsrf(i,nsrf).EQ.1.) THEN
1875
1876	IF(nreg.EQ.1) THEN
1877
1878	c TROP
1879	IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
1880	c pct_ocean(i,nreg)=.TRUE.
1881	pct_ocean(i,nreg)=1
1882	ENDIF
1883
1884	c PACIFIQUE NORD
1885	ELSEIF(nreg.EQ.2) THEN
1886	IF(rlat(i).GE.40.AND.rlat(i).LE.60.) THEN
1887	IF(rlonPOS(i).GE.160..AND.rlonPOS(i).LE.235.) THEN
1888	c pct_ocean(i,nreg)=.TRUE.
1889	pct_ocean(i,nreg)=1
1890	ENDIF
1891	ENDIF
1892	c CALIFORNIE ST-CU
1893	ELSEIF(nreg.EQ.3) THEN
1894	IF(rlonPOS(i).GE.220..AND.rlonPOS(i).LE.250.) THEN
1895	IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
1896	c pct_ocean(i,nreg)=.TRUE.
1897	pct_ocean(i,nreg)=1
1898	ENDIF
1899	ENDIF
1900	c HAWAI
1901	ELSEIF(nreg.EQ.4) THEN
1902	IF(rlonPOS(i).GE.180..AND.rlonPOS(i).LE.220.) THEN
1903	IF(rlat(i).GE.15.AND.rlat(i).LE.35.) THEN
1904	c pct_ocean(i,nreg)=.TRUE.
1905	pct_ocean(i,nreg)=1
1906	ENDIF
1907	ENDIF
1908	c WARM POOL
1909	ELSEIF(nreg.EQ.5) THEN
1910	IF(rlonPOS(i).GE.70..AND.rlonPOS(i).LE.150.) THEN
1911	IF(rlat(i).GE.-5.AND.rlat(i).LE.20.) THEN
1912	c pct_ocean(i,nreg)=.TRUE.
1913	pct_ocean(i,nreg)=1
1914	ENDIF
1915	ENDIF
1916	ENDIF !nbreg
1917	c TROP
1918	c IF(rlat(i).GE.-30.AND.rlat(i).LE.30.) THEN
1919	c pct_ocean(i)=.TRUE.
1920	c WRITE(,) 'pct_ocean =',i, rlon(i), rlat(i)
1921	c ENDIF !lon
1922	c ENDIF !lat
1923
1924	ENDIF !pctsrf
1925	c ENDDO
1926	ENDDO !klon
1927	ENDDO !nbreg
1928
1929	cIM somme de toutes les nhistoW BEG
1930	DO nreg = 1, nbreg
1931	DO k = 1, kmaxm1
1932	DO l = 1, lmaxm1
1933	DO iw = 1, iwmax
1934	nhistoWt(k,l,iw,nreg)=0.
1935	ENDDO
1936	ENDDO
1937	ENDDO
1938	ENDDO
1939	cIM somme de toutes les nhistoW END
1940	ENDIF
1941
1942
1943	c CALL ISCCP_CLOUD_TYPES(nlev,ncol,seed,pfull,phalf,qv,
1944	c & cc,conv,dtau_s,dtau_c,top_height,overlap,
1945	c & tautab,invtau,skt,emsfc_lw,at,dem_s,dem_c,fq_isccp,
1946	c & totalcldarea,meanptop,meantaucld,boxtau,boxptop)
1947
1948	c DO i=1, klon
1949	c i=1
1950	c1011 CONTINUE
1951	c
1952	cIM on verifie les donnees de INPUT en dehors du simulateur ISCCP
1953	cIM 1D non-vectorise (!) pour qu'on gagne du temps ...
1954	cIM
1955	c BEGIN find unpermittable data.....
1956	! ---------------------------------------------------!
1957	! find unpermittable data.....
1958	!
1959	c do 13 k=1,klev
1960	c ca prend trop de temps ??
1961	c cldtot(:,:) = min(max(cldtot(:,:),0.),1.)
1962	c rnebcon(:,:) = min(max(rnebcon(:,:),0.),1.)
1963	c dtau_s(:,:) = max(dtau_s(:,:),0.)
1964	c dem_s(:,:) = min(max(dem_s(:,:),0.),1.)
1965	c dtau_c(:,:) = max(dtau_c(:,:),0.)
1966	c dem_c(:,:) = min(max(dem_c(:,:),0.),1.)
1967	c ca prend trop de temps ??
1968
1969	c if (cldtot(i,k) .lt. 0.) then
1970	c print *, ' error = cloud fraction less than zero'
1971	c STOP
1972	c end if
1973	c if (cldtot(i,k) .gt. 1.) then
1974	c print *, ' error = cloud fraction greater than 1'
1975	c STOP
1976	c end if
1977	c if (rnebcon(i,k) .lt. 0.) then
1978	c print *,
1979	c & ' error = convective cloud fraction less than zero'
1980	c STOP
1981	c end if
1982	c if (rnebcon(i,k) .gt. 1.) then
1983	c print *,
1984	c & ' error = convective cloud fraction greater than 1'
1985	c STOP
1986	c end if
1987
1988	c if (dtau_s(i,k) .lt. 0.) then
1989	c print *,
1990	c & ' error = stratiform cloud opt. depth less than zero'
1991	c STOP
1992	c end if
1993	c if (dem_s(i,k) .lt. 0.) then
1994	c print *,
1995	c & ' error = stratiform cloud emissivity less than zero'
1996	c STOP
1997	c end if
1998	c if (dem_s(i,k) .gt. 1.) then
1999	c print *,
2000	c & ' error = stratiform cloud emissivity greater than 1'
2001	c STOP
2002	c end if
2003
2004	c if (dtau_c(i,k) .lt. 0.) then
2005	c print *,
2006	c & ' error = convective cloud opt. depth less than zero'
2007	c STOP
2008	c end if
2009	c if (dem_c(i,k) .lt. 0.) then
2010	c print *,
2011	c & ' error = convective cloud emissivity less than zero'
2012	c STOP
2013	c end if
2014	c if (dem_c(i,k) .gt. 1.) then
2015	c print *,
2016	c & ' error = convective cloud emissivity greater than 1'
2017	c STOP
2018	c end if
2019	c13 continue
2020
2021	! ---------------------------------------------------!
2022	c
2023	c END find unpermittable data.....
2024	cv2.2.1.1 DO i=1, klon
2025	c i=1
2026	c seed=i
2027	c
2028	cv3.4
2029	if (debut) then
2030	DO i=1, klon
2031	seed(i)=i+100
2032	c seed(i)=i+50
2033	ENDDO
2034	endif
2035	c seed=aint(ran0(klon))
2036	c CALL ISCCP_CLOUD_TYPES(klev,ncol,seed,pfull(i,:),phalf(i,:)
2037	cv2.2.1.1
2038	c CALL ISCCP_CLOUD_TYPES(klev,ncol,seed(i),pfull(i,:),phalf(i,:)
2039	c & ,q_seri(i,:),
2040	c & cldtot(i,:),rnebcon(i,:),dtau_s(i,:),dtau_c(i,:),
2041	c & top_height,overlap,
2042	c & tautab,invtau,ztsol,emsfc_lw,t_seri(i,:),dem_s(i,:),
2043	c & dem_c(i,:),
2044	c & fq_isccp(i,:,:),
2045	c & totalcldarea(i),meanptop(i),meantaucld(i),
2046	c & boxtau(i,:),boxptop(i,:))
2047	cv2.2.1.1
2048	cv3.4
2049	debug=0
2050	debugcol=0
2051	cIM260503
2052	c o500 ==> distribution nuage ftion du regime dynamique
2053	DO i=1, klon
2054	o500(i)=omega(i,8)*864.
2055	c PRINT*,'pphi8 ',pphi(i,8),'zphi8,11,12',zphi(i,8),
2056	c & zphi(i,11),zphi(i,12)
2057	ENDDO
2058
2059	c axe vertical pour les differents niveaux des histogrammes
2060	c DO iw=1, iwmax
2061	c zx_o500(iw)=wmin+(iw-1./2.)*pas_w
2062	c ENDDO
2063	c PRINT*,' phys AVANT seed(3361)=',seed(3361)
2064	CALL ISCCP_CLOUD_TYPES(
2065	& debug,
2066	& debugcol,
2067	& klon,
2068	& sunlit,
2069	& klev,
2070	& ncol,
2071	& seed,
2072	& pfull,
2073	& phalf,
2074	c var de bas en haut ==> PB !
2075	c & q_seri,
2076	c & cldtot,
2077	c & rnebcon,
2078	c & dtau_s,
2079	c & dtau_c,
2080	c var de Haut en Bas BEG
2081	& qv, cc, conv, dtau_sH2B, dtau_cH2B,
2082	c var de Haut en Bas END
2083	& top_height,
2084	& overlap,
2085	& tautab,
2086	& invtau,
2087	& ztsol,
2088	& emsfc_lw,
2089	c var de bas en haut ==> PB !
2090	c & t_seri,
2091	c & dem_s,
2092	c & dem_c,
2093	c var de Haut en Bas BEG
2094	& at, dem_sH2B, dem_cH2B,
2095	cIM260503
2096	c & o500, pct_ocean,
2097	c var de Haut en Bas END
2098	& fq_isccp,
2099	& totalcldarea,
2100	& meanptop,
2101	& meantaucld,
2102	& boxtau,
2103	& boxptop)
2104	c & boxptop,
2105	cIM 260503
2106	c & histoW,
2107	c & nhistoW
2108	c &)
2109
2110	cIM 200603
2111	c PRINT*,'physiq fq_isccp(6,1,1)',fq_isccp(6,1,1)
2112
2113	cIM 200603
2114	cIM somme de toutes les nhistoW BEG
2115	c DO k = 1, kmaxm1
2116	c DO l = 1, lmaxm1
2117	c DO iw = 1, iwmax
2118	c nhistoWt(k,l,iw)=nhistoWt(k,l,iw)+nhistoW(k,l,iw)
2119	ccc IF(k.EQ.1.AND.l.EQ.1.AND.iw.EQ.1) then
2120	c IF(nhistoWt(k,l,iw).NE.0.) THEN
2121	c PRINT*,' physiq nWt', k,l,iw,nhistoWt(k,l,iw)
2122	c ENDIF
2123	c ENDDO
2124	c ENDDO
2125	c ENDDO
2126	cIM somme de toutes les nhistoW END
2127	c PRINT*,' phys APRES seed(3361)=',seed(3361)
2128	cv3.4
2129	c i=i+1
2130	c IF(i.LE.klon) THEN
2131	c GOTO 1011
2132	c ENDIF
2133	cv2.2.1.1 ENDDO
2134
2135	c passage de la grille (klon,7,7) a (iim,jjmp1,7,7)
2136	c minfq3d=100.
2137	c maxfq3d=0.
2138	cIM calcul des correspondances entre la grille physique et
2139	cIM la grille dynamique
2140	c DO i=1, klon
2141	c grid_phys(i)=i
2142	c PRINT*,'i, grid_phys',i,grid_phys(i)
2143	c ENDDO
2144	c CALL gr_fi_dyn(1,klon,iimp1,jjmp1,grid_phys,grid_dyn)
2145	c DO j=1, jjmp1
2146	c DO i=1, iimp1
2147	c PRINT*,'i,j grid_dyn ',i,j,grid_dyn(i,j)
2148	c ENDDO
2149	c ENDDO
2150	c
2151	DO l=1, lmax
2152	DO k=1, kmax
2153	cIM grille physique ==> grille ecriture 2D (iim,jjmp1)
2154	c
2155	DO i=1, iim
2156	fq4d(i,1,k,l)=fq_isccp(1,k,l)
2157	cc PRINT*,'first j=1 i =',i
2158	ENDDO
2159	DO j=2, jjm
2160	DO i=1, iim
2161	cERROR ?? ig=i+iim*(j-1)
2162	ig=i+1+(j-2)*iim
2163	cc PRINT*,'i =',i,'j =',j,'ig=',ig
2164	fq4d(i,j,k,l)=fq_isccp(ig,k,l)
2165	ENDDO
2166	ENDDO
2167	DO i=1, iim
2168	fq4d(i,jjmp1,k,l)=fq_isccp(klon,k,l)
2169	cc PRINT*,'last jjmp1 i =',i
2170	ENDDO
2171	IF(debut) THEN
2172	DO j=1, jjmp1
2173	DO i=1, iim
2174	IF(j.GE.2.AND.j.LE.jjm) THEN
2175	igfi2D(i,j)=i+1+(j-2)*iim
2176	c PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
2177	ELSEIF(j.EQ.1) THEN
2178	igfi2D(i,j)=1
2179	c PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
2180	ELSEIF(j.EQ.jjmp1) THEN
2181	igfi2D(i,j)=klon
2182	c PRINT*,'i=',i,'j=',j,'ig=',igfi2D(i,j)
2183	ENDIF
2184	ENDDO
2185	ENDDO
2186	ENDIF
2187	c STOP
2188	c
2189	c CALL gr_fi_ecrit(1,klon,iim,jjmp1,fq_isccp(:,k,l),
2190	c $ fq4d(:,:,k,l))
2191	ENDDO
2192	ENDDO
2193	DO l=1, lmax
2194	fq4d(:,:,8,l)=-1.e+10
2195	fq4d(:,:,l,8)=-1.e+10
2196	ENDDO
2197	DO l=1, lmax
2198	DO k=1, kmax
2199	DO j=1, jjmp1
2200	DO i=1, iim
2201
2202	c inversion TAU ?!
2203	c ni=(i-1)*lmax+l
2204	c nj=(j-1)*kmax+kmax-k+1
2205	c
2206	c210503 inversion en PC ==> pas besoin !!!
2207	c ni=(i-1)*lmax+lmax-l+1
2208	c nj=(j-1)*kmax+k
2209	c
2210	c210503
2211	ni=(i-1)*lmax+l
2212	nj=(j-1)*kmax+k
2213
2214	c210503 if(k.EQ.8) then
2215	c fq4d(i,j,8,l)=-1.e+10
2216	c endif
2217
2218	c210503 if(l.EQ.8) then
2219	c fq4d(i,j,k,8)=-1.e+10
2220	c endif
2221
2222	fq3d(ni,nj)=fq4d(i,j,k,l)
2223
2224	c if(fq3d(ni,nj).lt.0.) then
2225	c print*,' fq3d LT ZERO ',ni,nj,fq3d(ni,nj)
2226	c endif
2227	c if(fq3d(ni,nj).gt.100.) then
2228	c print*,' fq3d GT 100 ',ni,nj,fq3d(ni,nj)
2229	c endif
2230	c max & min fq3d
2231	c if(fq3d(ni,nj).gt.maxfq3d) maxfq3d=fq3d(ni,nj)
2232	c if(fq3d(ni,nj).lt.minfq3d) minfq3d=fq3d(ni,nj)
2233
2234	ENDDO
2235	ENDDO
2236	c fq4d(:,:,8,l)=-1.e+10
2237	c fq4d(:,:,k,8)=-1.e+10
2238	c k=k+1
2239	c if(k.LE.kmax) then
2240	c goto 1022
2241	c endif
2242	ENDDO
2243	c l=l+1
2244	c if(l.LE.lmax) then
2245	c goto 1021
2246	c endif
2247	ENDDO
2248
2249	c print*,' minfq3d=',minfq3d,' maxfq3d=',maxfq3d
2250	c
2251	c calculs statistiques distribution nuage ftion du regime dynamique
2252	c DO i=1, klon
2253	c! o500(i)=omega(i,9)*864.
2254	c! PRINT*,' o500=',o500(i),' pphi(9)=',pphi(i,9)
2255	c o500(i)=omega(i,8)*864.
2256	cc PRINT*,' pphi(8)',pphi(i,8),'pphi(11)',pphi(i,11),
2257	cc .'pphi(12)',pphi(i,12)
2258	cc PRINT*,' zphi8,11,12=',zphi(i,8),zphi(i,11),zphi(i,12)
2259	cc PRINT*,' o500',o500(i),' w500',w500(i)
2260	c ENDDO
2261
2262	c axe vertical pour les differents niveaux des histogrammes
2263	c DO iw=1, iwmax
2264	c zx_o500(iw)=wmin+(iw-1./2.)*pas_w
2265	c ENDDO
2266
2267
2268	c Ce calcul doit etre fait a partir de valeurs mensuelles ??
2269	cc CALL histo_o500_pctau(o500,fq4d,histoW)
2270	cc CALL histo_o500_pctau(paire,pctsrf,o500,fq4d,histoW)
2271	cc CALL histo_o500_pctau(pct_ocean,rlat,o500,fq4d,histoW)
2272	ccOK ??? CALL histo_o500_pctau(pct_ocean,o500,fq4d,histoW)
2273	c CALL histo_o500_pctau(klon,pct_ocean,o500,fq4d,histoW,nhistoW)
2274	c CALL histo_o500_pctau(klon,pct_ocean,o500,fq_isccp,
2275	CALL histo_o500_pctau(nbreg,pct_ocean,o500,fq_isccp,
2276	&histoW,nhistoW)
2277	c
2278	cIM somme de toutes les nhistoW BEG
2279	DO nreg=1, nbreg
2280	DO k = 1, kmaxm1
2281	DO l = 1, lmaxm1
2282	DO iw = 1, iwmax
2283	nhistoWt(k,l,iw,nreg)=nhistoWt(k,l,iw,nreg)+
2284	& nhistoW(k,l,iw,nreg)
2285	ccc IF(k.EQ.1.AND.l.EQ.1.AND.iw.EQ.1) then
2286	c IF(nhistoWt(k,l,iw).NE.0.) THEN
2287	c PRINT*,' physiq nWt', k,l,iw,nhistoWt(k,l,iw)
2288	c ENDIF
2289	ENDDO
2290	ENDDO
2291	ENDDO
2292	ENDDO
2293	cIM somme de toutes les nhistoW END
2294	c
2295	c IF(lafin) THEN
2296	c DO nreg=1, nbreg
2297	c DO iw=1, iwmax
2298	c DO l=1,lmaxm1
2299	c DO k=1,kmaxm1
2300	c IF(histoW(k,l,iw,nreg).NE.0.) then
2301	c PRINT*,'physiq H nH',k,l,iw,
2302	c & histoW(k,l,iw,nreg),
2303	c & nhistoW(k,l,iw,nreg),nhistoWt(k,l,iw,nreg)
2304	c ENDIF
2305	c ENDDO
2306	c ENDDO
2307	c ENDDO
2308	c ENDDO
2309	cIM verif fq_isccp, fq4d, fq3d
2310	c DO l=1, lmaxm1
2311	c DO k=1,kmaxm1
2312	c i=74
2313	c j=36
2314	c DO j=1, jjmp1
2315	c DO i=1, iim
2316	c DO l=1, lmaxm1
2317	c WRITE(*,'(a,3i4,7f10.4)')
2318	c & 'fq_isccp,j,i,l=',j,i,l,
2319	c & (fq_isccp(igfi2D(i,j),k,l),k=1,kmaxm1)
2320	c WRITE(*,'(a,3i4,7f10.4)')
2321	c & 'fq4d,j,i,l=',j,i,l,(fq4d(i,j,k,l),k=1,kmaxm1)
2322	c ENDDO
2323	c ENDDO
2324	c ENDDO
2325	c ni1=(i-1)*8+1
2326	c ni2=i*8
2327	c nj1=(j-1)*8+1
2328	c nj2=j*8
2329	c DO ni=ni1,ni2
2330	c WRITE(*,'(a,2i4,7f10.4)')
2331	c & 'fq3d, ni,nj=',ni,nj,
2332	c & (fq3d(ni,nj),nj=nj1,nj2)
2333	c ENDDO
2334	c ENDIF
2335
2336	c DO iw=1, iwmax
2337	c DO l=1,lmaxm1
2338	c DO k=1,kmaxm1
2339	c PRINT*,' iw,l,k,nhistoW=',iw,l,k,nhistoW(k,l,iw)
2340	c ENDDO
2341	c ENDDO
2342	c ENDDO
2343
2344	c DO iw=1, iwmax
2345	c DO l=1, lmaxm1
2346	c linv=lmaxm1-l+1
2347	c DO k=1, kmaxm1
2348	c histoWinv(k,l,iw)=histoW(iw,k,l)
2349	c ENDDO
2350	c ENDDO
2351	c ENDDO
2352	c
2353	c pb syncronisation ?? : 48 * 30 * 7 (jour1) + 48* 29 * 7 (jour suivant)
2354	c
2355
2356
2357	ENDIF !ok_isccp
2358	cIM ISCCP simulator END
2359
2360	c On prend la somme des fractions nuageuses et des contenus en eau
2361	cldfra(:,:)=min(max(cldfra(:,:),rnebcon(:,:)),1.)
2362	cldliq(:,:)=cldliq(:,:)+rnebcon(:,:)*clwcon(:,:)
2363
2364
2365	ENDIF
2366	c
2367	c 2. NUAGES STARTIFORMES
2368	c
2369	IF (ok_stratus) THEN
2370	CALL diagcld2(paprs,pplay,t_seri,q_seri, diafra,dialiq)
2371	DO k = 1, klev
2372	DO i = 1, klon
2373	IF (diafra(i,k).GT.cldfra(i,k)) THEN
2374	cldliq(i,k) = dialiq(i,k)
2375	cldfra(i,k) = diafra(i,k)
2376	ENDIF
2377	ENDDO
2378	ENDDO
2379	ENDIF
2380	c
2381	c Precipitation totale
2382	c
2383	DO i = 1, klon
2384	rain_fall(i) = rain_con(i) + rain_lsc(i)
2385	snow_fall(i) = snow_con(i) + snow_lsc(i)
2386	ENDDO
2387	c
2388	IF (if_ebil.ge.2) THEN
2389	ztit="after diagcld"
2390	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2391	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2392	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2393	END IF
2394	c
2395	c Calculer l'humidite relative pour diagnostique
2396	c
2397	DO k = 1, klev
2398	DO i = 1, klon
2399	zx_t = t_seri(i,k)
2400	IF (thermcep) THEN
2401	zdelta = MAX(0.,SIGN(1.,rtt-zx_t))
2402	zx_qs = r2es * FOEEW(zx_t,zdelta)/pplay(i,k)
2403	zx_qs = MIN(0.5,zx_qs)
2404	zcor = 1./(1.-retv*zx_qs)
2405	zx_qs = zx_qs*zcor
2406	ELSE
2407	IF (zx_t.LT.t_coup) THEN
2408	zx_qs = qsats(zx_t)/pplay(i,k)
2409	ELSE
2410	zx_qs = qsatl(zx_t)/pplay(i,k)
2411	ENDIF
2412	ENDIF
2413	zx_rh(i,k) = q_seri(i,k)/zx_qs
2414	zqsat(i,k)=zx_qs
2415	ENDDO
2416	ENDDO
2417	c
2418	c Calculer les parametres optiques des nuages et quelques
2419	c parametres pour diagnostiques:
2420	c
2421	if (ok_newmicro) then
2422	CALL newmicro (paprs, pplay,ok_newmicro,
2423	. t_seri, cldliq, cldfra, cldtau, cldemi,
2424	. cldh, cldl, cldm, cldt, cldq)
2425	else
2426	CALL nuage (paprs, pplay,
2427	. t_seri, cldliq, cldfra, cldtau, cldemi,
2428	. cldh, cldl, cldm, cldt, cldq)
2429	endif
2430	c
2431	c Appeler le rayonnement mais calculer tout d'abord l'albedo du sol.
2432	c
2433	IF (MOD(itaprad,radpas).EQ.0) THEN
2434	DO i = 1, klon
2435	albsol(i) = falbe(i,is_oce) * pctsrf(i,is_oce)
2436	. + falbe(i,is_lic) * pctsrf(i,is_lic)
2437	. + falbe(i,is_ter) * pctsrf(i,is_ter)
2438	. + falbe(i,is_sic) * pctsrf(i,is_sic)
2439	albsollw(i) = falblw(i,is_oce) * pctsrf(i,is_oce)
2440	. + falblw(i,is_lic) * pctsrf(i,is_lic)
2441	. + falblw(i,is_ter) * pctsrf(i,is_ter)
2442	. + falblw(i,is_sic) * pctsrf(i,is_sic)
2443	ENDDO
2444	! if (debut) then
2445	! albsol1 = albsol
2446	! albsollw1 = albsollw
2447	! endif
2448	! albsol = albsol1
2449	! albsollw = albsollw1
2450	CALL radlwsw ! nouveau rayonnement (compatible Arpege-IFS)
2451	cIM e (dist, rmu0, fract, co2_ppm, solaire,
2452	e (dist, rmu0, fract,
2453	e paprs, pplay,zxtsol,albsol, albsollw, t_seri,q_seri,
2454	e wo,
2455	e cldfra, cldemi, cldtau,
2456	s heat,heat0,cool,cool0,radsol,albpla,
2457	s topsw,toplw,solsw,sollw,
2458	s sollwdown,
2459	cccIMs topsw0,toplw0,solsw0,sollw0)
2460	s topsw0,toplw0,solsw0,sollw0,
2461	s swdn0, swdn, swup0, swup )
2462	itaprad = 0
2463	ENDIF
2464	itaprad = itaprad + 1
2465
2466	c
2467	c Ajouter la tendance des rayonnements (tous les pas)
2468	c
2469	DO k = 1, klev
2470	DO i = 1, klon
2471	t_seri(i,k) = t_seri(i,k)
2472	. + (heat(i,k)-cool(i,k)) * dtime/86400.
2473	ENDDO
2474	ENDDO
2475	c
2476	IF (if_ebil.ge.2) THEN
2477	ztit='after rad'
2478	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2479	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2480	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2481	call diagphy(paire,ztit,ip_ebil
2482	e , topsw, toplw, solsw, sollw, zero_v
2483	e , zero_v, zero_v, zero_v, ztsol
2484	e , d_h_vcol, d_qt, d_ec
2485	s , fs_bound, fq_bound )
2486	END IF
2487	c
2488	c
2489	c Calculer l'hydrologie de la surface
2490	c
2491	c CALL hydrol(dtime,pctsrf,rain_fall, snow_fall, zxevap,
2492	c . agesno, ftsol,fqsurf,fsnow, ruis)
2493	c
2494	DO i = 1, klon
2495	zxqsurf(i) = 0.0
2496	zxsnow(i) = 0.0
2497	ENDDO
2498	DO nsrf = 1, nbsrf
2499	DO i = 1, klon
2500	zxqsurf(i) = zxqsurf(i) + fqsurf(i,nsrf)*pctsrf(i,nsrf)
2501	zxsnow(i) = zxsnow(i) + fsnow(i,nsrf)*pctsrf(i,nsrf)
2502	ENDDO
2503	ENDDO
2504	c
2505	c Si une sous-fraction n'existe pas, elle prend la valeur moyenne
2506	c
2507	cXXX DO nsrf = 1, nbsrf
2508	cXXX DO i = 1, klon
2509	cXXX IF (pctsrf(i,nsrf).LT.epsfra) THEN
2510	cXXX fqsurf(i,nsrf) = zxqsurf(i)
2511	cXXX fsnow(i,nsrf) = zxsnow(i)
2512	cXXX ENDIF
2513	cXXX ENDDO
2514	cXXX ENDDO
2515	c
2516	c Calculer le bilan du sol et la derive de temperature (couplage)
2517	c
2518	DO i = 1, klon
2519	c bils(i) = radsol(i) - sens(i) - evap(i)*RLVTT
2520	c a la demande de JLD
2521	bils(i) = radsol(i) - sens(i) + zxfluxlat(i)
2522	ENDDO
2523	c
2524	cmoddeblott(jan95)
2525	c Appeler le programme de parametrisation de l'orographie
2526	c a l'echelle sous-maille:
2527	c
2528	IF (ok_orodr) THEN
2529	c
2530	c selection des points pour lesquels le shema est actif:
2531	igwd=0
2532	DO i=1,klon
2533	itest(i)=0
2534	c IF ((zstd(i).gt.10.0)) THEN
2535	IF (((zpic(i)-zmea(i)).GT.100.).AND.(zstd(i).GT.10.0)) THEN
2536	itest(i)=1
2537	igwd=igwd+1
2538	idx(igwd)=i
2539	ENDIF
2540	ENDDO
2541	c igwdim=MAX(1,igwd)
2542	c
2543	CALL drag_noro(klon,klev,dtime,paprs,pplay,
2544	e zmea,zstd, zsig, zgam, zthe,zpic,zval,
2545	e igwd,idx,itest,
2546	e t_seri, u_seri, v_seri,
2547	s zulow, zvlow, zustr, zvstr,
2548	s d_t_oro, d_u_oro, d_v_oro)
2549	c
2550	c ajout des tendances
2551	DO k = 1, klev
2552	DO i = 1, klon
2553	t_seri(i,k) = t_seri(i,k) + d_t_oro(i,k)
2554	u_seri(i,k) = u_seri(i,k) + d_u_oro(i,k)
2555	v_seri(i,k) = v_seri(i,k) + d_v_oro(i,k)
2556	ENDDO
2557	ENDDO
2558	c
2559	ENDIF ! fin de test sur ok_orodr
2560	c
2561	IF (ok_orolf) THEN
2562	c
2563	c selection des points pour lesquels le shema est actif:
2564	igwd=0
2565	DO i=1,klon
2566	itest(i)=0
2567	IF ((zpic(i)-zmea(i)).GT.100.) THEN
2568	itest(i)=1
2569	igwd=igwd+1
2570	idx(igwd)=i
2571	ENDIF
2572	ENDDO
2573	c igwdim=MAX(1,igwd)
2574	c
2575	CALL lift_noro(klon,klev,dtime,paprs,pplay,
2576	e rlat,zmea,zstd,zpic,
2577	e itest,
2578	e t_seri, u_seri, v_seri,
2579	s zulow, zvlow, zustr, zvstr,
2580	s d_t_lif, d_u_lif, d_v_lif)
2581	c
2582	c ajout des tendances
2583	DO k = 1, klev
2584	DO i = 1, klon
2585	t_seri(i,k) = t_seri(i,k) + d_t_lif(i,k)
2586	u_seri(i,k) = u_seri(i,k) + d_u_lif(i,k)
2587	v_seri(i,k) = v_seri(i,k) + d_v_lif(i,k)
2588	ENDDO
2589	ENDDO
2590	c
2591	ENDIF ! fin de test sur ok_orolf
2592	c
2593	IF (if_ebil.ge.2) THEN
2594	ztit='after orography'
2595	CALL diagetpq(paire,ztit,ip_ebil,2,2,dtime
2596	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2597	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2598	END IF
2599	c
2600	c
2601	cAA
2602	cAA Installation de l'interface online-offline pour traceurs
2603	cAA
2604	c====================================================================
2605	c Calcul des tendances traceurs
2606	c====================================================================
2607	C Pascale : il faut quand meme apeller phytrac car il gere les sorties
2608	cKE43 des traceurs => il faut donc mettre des flags a .false.
2609	IF (iflag_con.GE.3) THEN
2610	c on ajoute les tendances calculees par KE43
2611	cXXX OM on onhibe la convection sur les traceurs
2612	DO iq=1, nqmax-2 ! Sandrine a -3 ???
2613	cXXX OM on inhibe la convection sur les traceur
2614	cXXX DO k = 1, nlev
2615	cXXX DO i = 1, klon
2616	cXXX tr_seri(i,k,iq) = tr_seri(i,k,iq) + d_tr(i,k,iq)
2617	cXXX ENDDO
2618	cXXX ENDDO
2619	WRITE(iqn,'(i2.2)') iq
2620	CALL minmaxqfi(tr_seri(1,1,iq),0.,1.e33,'couche lim iq='//iqn)
2621	ENDDO
2622	CMAF modif pour garder info du nombre de traceurs auxquels
2623	C la physique s'applique
2624	ELSE
2625	CMAF modif pour garder info du nombre de traceurs auxquels
2626	C la physique s'applique
2627	C
2628	call phytrac (rnpb,
2629	I debut,lafin,
2630	I nqmax-2,
2631	I nlon,nlev,dtime,
2632	I t,paprs,pplay,
2633	I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2634	I ycoefh,yu1,yv1,ftsol,pctsrf,rlat,
2635	I frac_impa, frac_nucl,
2636	I rlon,presnivs,paire,pphis,
2637	O tr_seri)
2638	ENDIF
2639
2640	IF (offline) THEN
2641
2642	call phystokenc (
2643	I nlon,nlev,pdtphys,rlon,rlat,
2644	I t,pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
2645	I ycoefh,yu1,yv1,ftsol,pctsrf,
2646	I frac_impa, frac_nucl,
2647	I pphis,paire,dtime,itap)
2648
2649
2650	ENDIF
2651
2652	c
2653	c Calculer le transport de l'eau et de l'energie (diagnostique)
2654	c
2655	CALL transp (paprs,zxtsol,
2656	e t_seri, q_seri, u_seri, v_seri, zphi,
2657	s ve, vq, ue, uq)
2658	c
2659	c Accumuler les variables a stocker dans les fichiers histoire:
2660	c
2661	c
2662	c
2663	c+jld ec_conser
2664	DO k = 1, klev
2665	DO i = 1, klon
2666	ZRCPD = RCPD(1.0+RVTMP2q_seri(i,k))
2667	d_t_ec(i,k)=0.5/ZRCPD
2668	$ (u(i,k)2+v(i,k)2-u_seri(i,k)2-v_seri(i,k)*2)
2669	t_seri(i,k)=t_seri(i,k)+d_t_ec(i,k)
2670	d_t_ec(i,k) = d_t_ec(i,k)/dtime
2671	END DO
2672	END DO
2673	c-jld ec_conser
2674	IF (if_ebil.ge.1) THEN
2675	ztit='after physic'
2676	CALL diagetpq(paire,ztit,ip_ebil,1,1,dtime
2677	e , t_seri,q_seri,ql_seri,qs_seri,u_seri,v_seri,paprs,pplay
2678	s , d_h_vcol, d_qt, d_qw, d_ql, d_qs, d_ec)
2679	C Comme les tendances de la physique sont ajoute dans la dynamique,
2680	C on devrait avoir que la variation d'entalpie par la dynamique
2681	C est egale a la variation de la physique au pas de temps precedent.
2682	C Donc la somme de ces 2 variations devrait etre nulle.
2683	call diagphy(paire,ztit,ip_ebil
2684	e , topsw, toplw, solsw, sollw, sens
2685	e , evap, rain_fall, snow_fall, ztsol
2686	e , d_h_vcol, d_qt, d_ec
2687	s , fs_bound, fq_bound )
2688	C
2689	d_h_vcol_phy=d_h_vcol
2690	C
2691	END IF
2692	C
2693	cccIM cf. FH
2694	c=======================================================================
2695	c SORTIES
2696	c=======================================================================
2697
2698	c Interpollation sur quelques niveaux de pression
2699	c -----------------------------------------------
2700
2701	call plevel(klon,klev,.true. ,pplay,85000.,u_seri,u850)
2702	call plevel(klon,klev,.false.,pplay,85000.,v_seri,v850)
2703	call plevel(klon,klev,.true. ,pplay,50000.,u_seri,u500)
2704	call plevel(klon,klev,.false.,pplay,50000.,v_seri,v500)
2705	call plevel(klon,klev,.true. ,pplay,20000.,u_seri,u200)
2706	call plevel(klon,klev,.false.,pplay,20000.,v_seri,v200)
2707	call plevel(klon,klev,.true. ,pplay,50000.,zphi,phi500)
2708	call plevel(klon,klev,.true. ,paprs,50000.,omega,w500)
2709
2710	cIM cf. FH slp(:) = paprs(:,1)exp(pphis(:)/(289.t_seri(:,1)))
2711	slp(:) = paprs(:,1)exp(pphis(:)/(RDt_seri(:,1)))
2712	c PRINT*,' physiq slp ',slp(2185),paprs(2185,1),pphis(2185),
2713	c . RD,t_seri(2185,1)
2714	c
2715	ccc prw = eau precipitable
2716	DO i = 1, klon
2717	prw(i) = 0.
2718	DO k = 1, klev
2719	prw(i) = prw(i) +
2720	. q_seri(i,k)*(paprs(i,k)-paprs(i,k+1))/RG
2721	ENDDO
2722	c PRINT*,' i ',i,' prw',prw(i)
2723	ENDDO
2724	c
2725
2726	c=============================================================
2727	c Ecriture des sorties
2728	c=============================================================
2729
2730	#ifdef histISCCP
2731	#include "write_histISCCP.h"
2732	#endif
2733
2734	#ifdef histhf
2735	#include "write_histhf.h"
2736	#endif
2737
2738	#include "write_histday.h"
2739	#include "write_histmth.h"
2740	#include "write_histins.h"
2741
2742	c=============================================================
2743	c
2744	c Convertir les incrementations en tendances
2745	c
2746	DO k = 1, klev
2747	DO i = 1, klon
2748	d_u(i,k) = ( u_seri(i,k) - u(i,k) ) / dtime
2749	d_v(i,k) = ( v_seri(i,k) - v(i,k) ) / dtime
2750	d_t(i,k) = ( t_seri(i,k)-t(i,k) ) / dtime
2751	d_qx(i,k,ivap) = ( q_seri(i,k) - qx(i,k,ivap) ) / dtime
2752	d_qx(i,k,iliq) = ( ql_seri(i,k) - qx(i,k,iliq) ) / dtime
2753	ENDDO
2754	ENDDO
2755	c
2756	IF (nqmax.GE.3) THEN
2757	DO iq = 3, nqmax
2758	DO k = 1, klev
2759	DO i = 1, klon
2760	d_qx(i,k,iq) = ( tr_seri(i,k,iq-2) - qx(i,k,iq) ) / dtime
2761	ENDDO
2762	ENDDO
2763	ENDDO
2764	ENDIF
2765	c
2766	c Sauvegarder les valeurs de t et q a la fin de la physique:
2767	c
2768	DO k = 1, klev
2769	DO i = 1, klon
2770	t_ancien(i,k) = t_seri(i,k)
2771	q_ancien(i,k) = q_seri(i,k)
2772	ENDDO
2773	ENDDO
2774	c
2775	c====================================================================
2776	c Si c'est la fin, il faut conserver l'etat de redemarrage
2777	c====================================================================
2778	c
2779	IF (lafin) THEN
2780	itau_phy = itau_phy + itap
2781	ccc IF (ok_oasis) CALL quitcpl
2782	CALL phyredem ("restartphy.nc",dtime,radpas,
2783	. rlat, rlon, pctsrf, ftsol, ftsoil, deltat, fqsurf, qsol,
2784	. fsnow, falbe,falblw, fevap, rain_fall, snow_fall,
2785	. solsw, sollwdown,dlw,
2786	. radsol,frugs,agesno,
2787	. zmea,zstd,zsig,zgam,zthe,zpic,zval,rugoro,
2788	. t_ancien, q_ancien, rnebcon, ratqs, clwcon)
2789	ENDIF
2790
2791	RETURN
2792	END
2793	FUNCTION qcheck(klon,klev,paprs,q,ql,aire)
2794	IMPLICIT none
2795	c
2796	c Calculer et imprimer l'eau totale. A utiliser pour verifier
2797	c la conservation de l'eau
2798	c
2799	#include "YOMCST.h"
2800	INTEGER klon,klev
2801	REAL paprs(klon,klev+1), q(klon,klev), ql(klon,klev)
2802	REAL aire(klon)
2803	REAL qtotal, zx, qcheck
2804	INTEGER i, k
2805	c
2806	zx = 0.0
2807	DO i = 1, klon
2808	zx = zx + aire(i)
2809	ENDDO
2810	qtotal = 0.0
2811	DO k = 1, klev
2812	DO i = 1, klon
2813	qtotal = qtotal + (q(i,k)+ql(i,k)) * aire(i)
2814	. *(paprs(i,k)-paprs(i,k+1))/RG
2815	ENDDO
2816	ENDDO
2817	c
2818	qcheck = qtotal/zx
2819	c
2820	RETURN
2821	END
2822	SUBROUTINE gr_fi_ecrit(nfield,nlon,iim,jjmp1,fi,ecrit)
2823	IMPLICIT none
2824	c
2825	c Tranformer une variable de la grille physique a
2826	c la grille d'ecriture
2827	c
2828	INTEGER nfield,nlon,iim,jjmp1, jjm
2829	REAL fi(nlon,nfield), ecrit(iim*jjmp1,nfield)
2830	c
2831	INTEGER i, n, ig
2832	c
2833	jjm = jjmp1 - 1
2834	DO n = 1, nfield
2835	DO i=1,iim
2836	ecrit(i,n) = fi(1,n)
2837	ecrit(i+jjm*iim,n) = fi(nlon,n)
2838	ENDDO
2839	DO ig = 1, nlon - 2
2840	ecrit(iim+ig,n) = fi(1+ig,n)
2841	ENDDO
2842	ENDDO
2843	RETURN
2844	END
2845

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

Download in other formats:

Original Format