Context Navigation

physiq.F @ 544

Last change on this file since 544 was 544, checked in by lmdzadmin, 21 years ago
Incorporation des modifications necessaires a l'utilisation de la librairie Psmile/PRISM, et creation d'un tag IPSL-CM4_PSMILE, selon M.-E. Demory LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 86.5 KB

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

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

Download in other formats:

Original Format