Context Navigation

physiq.F @ 523

Last change on this file since 523 was 523, checked in by lmdzadmin, 20 years ago

Dernieres modifs (ha, ha) sur le tag IPSL-CM4_IPCC, MED & JLD:

avant 1930 on force la lecture des sulfates naturels
bug de debordement de tableau dans physiq.F
taucalv passe à 10 ans

Property svn:eol-style set to native
Property svn:executable set to *
Property svn:keywords set to Author Date Id Revision

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

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

Download in other formats:

Original Format