Context Navigation

physiq.F @ 517

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

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

Download in other formats:

Original Format