Context Navigation

physiq.F @ 478

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

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

Download in other formats:

Original Format