source: LMDZ.3.3/branches/rel-1-0-patch/libf/phylmd/phytrac.F @ 5235

Last change on this file since 5235 was 349, checked in by lmdz, 23 years ago

Regle le probleme de decalage de 1 jour au debut de chaque simulation
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 26.0 KB
Line 
1c
2c $Header$
3c
4      SUBROUTINE phytrac (rnpb,
5     I                   debutphy,lafin,
6     I                   nqmax,
7     I                   nlon,nlev,pdtphys,
8     I                   t_seri,paprs,pplay,
9     I                   pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
10     I                   coefh,yu1,yv1,ftsol,pctsrf,xlat,
11     I                   frac_impa, frac_nucl,
12     I                   xlon,presnivs,paire,pphis,
13     O                   tr_seri)
14      USE ioipsl
15      IMPLICIT none
16c======================================================================
17c Auteur(s) FH
18c Objet: Moniteur general des tendances traceurs
19c
20cAA Remarques en vrac:
21cAA--------------------
22cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
23cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
24cAA 2/ Le choix du radon et du pb se fait juste avec un data
25cAA    (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
26cAA    une variable "type de traceur"
27c======================================================================
28#include "YOMCST.h"
29#include "dimensions.h"
30#include "dimphy.h"
31#include "indicesol.h"
32#include "temps.h"
33#include "control.h"
34c======================================================================
35
36c Arguments:
37c
38c   EN ENTREE:
39c   ==========
40c
41c   divers:
42c   -------
43c
44      integer nlon  ! nombre de points horizontaux
45      integer nlev  ! nombre de couches verticales
46      integer nqmax ! nombre de traceurs auxquels on applique la physique
47      real pdtphys  ! pas d'integration pour la physique (seconde)
48      real t_seri(nlon,nlev) ! temperature
49      real tr_seri(nlon,nlev,nbtr) ! traceur 
50      real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
51      real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
52      real presnivs(klev) ! pressions approximat. des milieux couches ( en PA)
53      real paire(klon)
54      real pphis(klon)
55      logical debutphy       ! le flag de l'initialisation de la physique
56      logical lafin          ! le flag de la fin de la physique
57
58      integer ll
59c
60cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
61c
62c   convection:
63c   -----------
64c
65      REAL pmfu(nlon,nlev)  ! flux de masse dans le panache montant
66      REAL pmfd(nlon,nlev)  ! flux de masse dans le panache descendant
67      REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
68      REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
69      REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
70      REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
71c
72c   Couche limite:
73c   --------------
74c
75      REAL coefh(nlon,nlev) ! coeff melange CL
76      REAL yu1(nlon)        ! vents au premier niveau
77      REAL yv1(nlon)        ! vents au premier niveau
78      REAL xlat(nlon)       ! latitudes pour chaque point
79      REAL xlon(nlon)       ! longitudes pour chaque point
80
81c
82c   Lessivage:
83c   ----------
84c
85c pour le ON-LINE
86c
87      REAL frac_impa(nlon,nlev)  ! fraction d'aerosols impactes
88      REAL frac_nucl(nlon,nlev)  ! fraction d'aerosols nuclees
89c
90cAA
91cAA Arguments necessaires pour les sources et puits de traceur:
92cAA ----------------
93cAA
94      real ftsol(nlon,nbsrf)  ! Temperature du sol (surf)(Kelvin)
95      real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
96c abder
97      real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon)
98      real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon)
99c fin
100cAA ----------------------------
101cAA  VARIABLES LOCALES TRACEURS
102cAA ----------------------------
103cAA
104cAA Sources et puits des traceurs:
105cAA ------------------------------
106cAA
107cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
108
109      REAL source(klon)       ! a voir lorsque le flux est prescrit
110      REAL topuits(klev,nbtr)  ! a voir
111cAA
112cAA Pour la source de radon et son reservoir de sol
113cAA ................................................
114 
115      REAL trs(klon,nbtr)    ! Conc. radon ds le sol
116      SAVE trs
117
118      REAL masktr(klon,nbtr) ! Masque reservoir de sol traceur
119c                            Masque de l'echange avec la surface
120c                           (1 = reservoir) ou (possible => 1 )
121      SAVE masktr
122      REAL fshtr(klon,nbtr)  ! Flux surfacique dans le reservoir de sol
123      SAVE fshtr
124      REAL hsoltr(nbtr)      ! Epaisseur equivalente du reservoir de sol
125      SAVE hsoltr
126      REAL tautr(nbtr)       ! Constante de decroissance radioactive
127      SAVE tautr
128      REAL vdeptr(nbtr)      ! Vitesse de depot sec dans la couche Brownienne
129      SAVE vdeptr
130      REAL scavtr(nbtr)      ! Coefficient de lessivage
131      SAVE scavtr
132cAA
133      CHARACTER*1 itn
134C maf ioipsl
135      CHARACTER*2 str2
136      INTEGER nhori, nvert
137      REAL zsto, zout, zjulian
138      integer idayref
139      INTEGER nid_tra
140      SAVE nid_tra
141c     REAL x(klon,klev,nbtr+2) ! traceurs
142      INTEGER ndex(1)
143      REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
144      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
145      INTEGER itra
146      SAVE itra                   ! compteur pour la physique
147C
148C Variables liees a l'ecriture de la bande histoire : phytrac.nc
149c
150      INTEGER ecrit_tra
151      SAVE ecrit_tra   
152      logical ok_sync
153      parameter (ok_sync = .true.)
154C
155C nature du traceur
156c
157      logical aerosol(nbtr)  ! Nature du traceur
158c                            ! aerosol(it) = true  => aerosol
159c                            ! aerosol(it) = false => gaz
160c                            ! nat_trac(it) = 1. aerosolc
161      logical clsol(nbtr)    ! clsol(it) = true => CL sol calculee
162      logical radio(nbtr)    ! radio(it)=true => decroisssance radioactive
163      save aerosol,clsol,radio
164C
165c======================================================================
166c
167c Declaration des procedures appelees
168c
169c--modif convection tiedtke
170      INTEGER i, k, it,itap
171        save itap
172      REAL delp(klon,klev)
173c--end modif
174c
175c Variables liees a l'ecriture de la bande histoire physique
176c
177c Variables locales pour effectuer les appels en serie
178c----------------------------------------------------
179c
180      REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
181      REAL d_tr_cl(klon,klev) ! tendance de traceurs  couche limite
182      REAL d_tr_cv(klon,klev) ! tendance de traceurs  convection
183      REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
184c                                   ! radioactive du rn - > pb
185      REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
186c                                          ! par impaction
187      REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
188c                                          ! par nucleation
189      REAL fluxrn(klon,klev)
190      REAL fluxpb(klon,klev)
191      REAL pbimpa(klon,klev)
192      REAL pbnucl(klon,klev)
193      REAL rn(klon,klev)
194      REAL pb(klon,klev)
195      REAL flestottr(klon,klev,nbtr) ! flux de lessivage
196c                                    ! dans chaque couche
197
198C
199      character*20 modname
200      character*80 abort_message
201c
202c   Controles
203c-------------
204      logical first,couchelimite,convection,lessivage,sorties,
205     s        rnpb,inirnpb
206      save first,couchelimite,convection,lessivage,sorties,
207     s     inirnpb
208      data first,couchelimite,convection,lessivage,sorties
209     s     /.true.,.true.,.true.,.true.,.true./
210c
211c======================================================================
212c         ecrit_tra = NINT(86400./pdtphys *1.0)  ! tous les jours
213         modname='phytrac'
214
215c        print*,'DANS PHYTRAC debutphy=',debutphy
216
217         if (debutphy) then
218
219          print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
220          ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H
221c         ecrit_tra = NINT(86400./pdtphys) ! tous les 24H
222
223         if(nbtr.lt.nqmax) then
224           print*,'NQMAX=',nqmax
225           print*,'NBTR=',nbtr
226           abort_message='See above'
227           call abort_gcm(modname,abort_message,1)
228         endif
229
230         inirnpb=rnpb
231         PRINT*, 'La frequence de sortie traceurs est  ', ecrit_tra
232         itra=0
233         itap=0
234C         
235         idayref = day_ini
236         CALL ymds2ju(anne_ini, 1, idayref, 0.0, zjulian)
237c
238         CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon)
239         DO i = 1, iim
240            zx_lon(i,1) = xlon(i+1)
241            zx_lon(i,jjm+1) = xlon(i+1)
242         ENDDO
243c         DO ll=1,klev
244c            znivsig(ll)=float(ll)
245c         ENDDO
246         CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
247         CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
248     .                 1,iim,1,jjm+1, 0, zjulian, pdtphys,
249     .                 nhori, nid_tra)
250         CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
251     .                 klev, presnivs, nvert)
252         zsto = pdtphys
253         zout = pdtphys * FLOAT(ecrit_tra)
254c
255         CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
256     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
257     .                "once",  zsto,zout)
258c
259         CALL histdef(nid_tra, "aire", "Grid area", "-",
260     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
261     .                "once",  zsto,zout)
262
263        goto 666
264         CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-",
265     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
266     .                "inst(X)",  zsto,zout)
267
268         CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-",
269     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
270     .                "inst(X)",  zsto,zout)
271         CALL histdef(nid_tra, "psrf1", "nature sol", "-",
272     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
273     .                "inst(X)",  zsto,zout)
274         CALL histdef(nid_tra, "psrf2", "nature sol", "-",
275     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
276     .                "inst(X)",  zsto,zout)
277         CALL histdef(nid_tra, "psrf3", "nature sol", "-",
278     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
279     .                "inst(X)",  zsto,zout)
280         CALL histdef(nid_tra, "psrf4", "nature sol", "-",
281     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
282     .                "inst(X)",  zsto,zout)
283         CALL histdef(nid_tra, "ftsol1", "temper sol", "-",
284     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
285     .                "inst(X)",  zsto,zout)
286         CALL histdef(nid_tra, "ftsol2", "temper sol", "-",
287     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
288     .                "inst(X)",  zsto,zout)
289         CALL histdef(nid_tra, "ftsol3", "temper sol", "-",
290     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
291     .                "inst",  zsto,zout)
292         CALL histdef(nid_tra, "ftsol4", "temper sol", "-",
293     .                iim,jjm+1,nhori, 1,1,1, -99, 32,
294     .                "inst(X)",  zsto,zout)
295         CALL histdef(nid_tra, "pplay", "flux u mont","-",
296     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
297     .                "inst(X)", zsto,zout)
298         CALL histdef(nid_tra, "t", "flux u mont","-",
299     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
300     .                "inst(X)", zsto,zout)
301         CALL histdef(nid_tra, "mfu", "flux u mont","-",
302     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
303     .                "ave(X)", zsto,zout)
304         CALL histdef(nid_tra, "mfd", "flux u decen","-",
305     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
306     .                "ave(X)", zsto,zout)
307         CALL histdef(nid_tra, "en_u", "flux u mont","-",
308     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
309     .                "ave(X)", zsto,zout)
310         CALL histdef(nid_tra, "en_d", "flux u mont","-",
311     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
312     .                "ave(X)", zsto,zout)
313         CALL histdef(nid_tra, "de_u", "flux u mont","-",
314     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
315     .                "ave(X)", zsto,zout)
316         CALL histdef(nid_tra, "de_d", "flux u mont","-",
317     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
318     .                "ave(X)", zsto,zout)
319         CALL histdef(nid_tra, "coefh", "turbulent coef","-",
320     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
321     .                "ave(X)", zsto,zout)
322
323666     continue
324c
325         DO it=1,nqmax
326         IF (it.LE.99) THEN
327         WRITE(str2,'(i2.2)') it
328C champ 3D
329         CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga",
330     .                iim,jjm+1,nhori, klev,1,klev,nvert, 32,
331     .                "ave(X)", zsto,zout)
332         IF (lessivage) THEN
333         CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2,
334     .                "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32,
335     .                "ave(X)", zsto,zout)
336         ENDIF
337         ELSE
338         PRINT*, "Trop de traceurs"
339         CALL abort
340         ENDIF
341         ENDDO
342         CALL histend(nid_tra)
343         ndex(1) = 0
344c
345         i = NINT(zout/zsto)
346         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
347         CALL histwrite(nid_tra,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex)
348C
349         i = NINT(zout/zsto)
350         CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
351         CALL histwrite(nid_tra,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex)
352
353c======================================================================
354c   Initialisation de certaines variables pour le Rn et le Pb
355c======================================================================
356
357c Initialisation du traceur dans le sol (couche limite radonique)
358c
359c        print*,'valeur de debut dans phytrac :',debutphy
360         do it=1,nqmax
361            do i=1,klon
362               trs(i,it) = 0.
363            enddo
364         END DO
365
366         open (99,file='starttrac',status='old',
367     .         err=999,form='formatted')
368         read(99,*) (trs(i,1),i=1,klon)
369999      close(99)
370         print*, 'apres starttrac'
371
372c Initialisation de la fraction d'aerosols lessivee
373c
374         DO it=1,nqmax
375           DO k = 1, nlev
376              DO i = 1, klon
377                 d_tr_lessi_impa(i,k,it) =  0.0
378                 d_tr_lessi_nucl(i,k,it) =  0.0
379             ENDDO
380           ENDDO
381         ENDDO
382c
383c Initialisation de la nature des traceurs
384c
385         DO it = 1, nqmax
386            aerosol(it) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
387            radio(it) = .FALSE.    ! Par defaut pas de passage par radiornpb
388            clsol(it) = .FALSE.  ! Par defaut couche limite avec flux prescrit
389         ENDDO
390c
391      ENDIF  ! fin debutphy
392c Initialisation du traceur dans le sol (couche limite radonique)
393      if(inirnpb) THEN
394c
395         radio(1)= .true.
396         radio(2)= .true.
397         clsol(1)= .true.
398         clsol(2)= .true.
399         aerosol(2) = .TRUE. ! le Pb est un aerosol
400c
401         call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
402     .                   ,vdeptr,scavtr)
403         inirnpb=.false.
404      endif
405      if(nqmax.gt.2) aerosol(3)=.true.
406
407
408c  abder
409        goto 777
410            do i=1,nlon
411               pftsol1(i) = ftsol(i,1)
412               pftsol2(i) = ftsol(i,2)
413               pftsol3(i) = ftsol(i,3)
414               pftsol4(i) = ftsol(i,4)
415
416               ppsrf1(i) = pctsrf(i,1)
417               ppsrf2(i) = pctsrf(i,2)
418               ppsrf3(i) = pctsrf(i,3)
419               ppsrf4(i) = pctsrf(i,4)
420
421            enddo
422         ndex(1)=0
423         itap=itap+1
424         CALL gr_fi_ecrit(1,klon,iim,jjm+1,yu1,zx_tmp_2d)
425         CALL histwrite(nid_tra,"pyu1",itap,zx_tmp_2d,
426     s                                  iim*(jjm+1),ndex)
427         
428         CALL gr_fi_ecrit(1,klon,iim,jjm+1,yv1,zx_tmp_2d)
429         CALL histwrite(nid_tra,"pyv1",itap,zx_tmp_2d,
430     s                                  iim*(jjm+1),ndex)
431
432         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol1,zx_tmp_2d)
433         CALL histwrite(nid_tra,"ftsol1",itap,zx_tmp_2d,
434     s                                       iim*(jjm+1),ndex)
435
436         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol2,zx_tmp_2d)
437         CALL histwrite(nid_tra,"ftsol2",itap,zx_tmp_2d,
438     s                                       iim*(jjm+1),ndex)
439
440         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol3,zx_tmp_2d)
441         CALL histwrite(nid_tra,"ftsol3",itap,zx_tmp_2d,
442     s                                      iim*(jjm+1),ndex)
443
444         CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol4,zx_tmp_2d)
445         CALL histwrite(nid_tra,"ftsol4",itap,zx_tmp_2d,
446     s                                      iim*(jjm+1),ndex)
447
448         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf1,zx_tmp_2d)
449         CALL histwrite(nid_tra,"psrf1",itap,zx_tmp_2d,
450     s                                     iim*(jjm+1),ndex)
451
452         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf2,zx_tmp_2d)
453         CALL histwrite(nid_tra,"psrf2",itap,zx_tmp_2d,
454     s                                     iim*(jjm+1),ndex)
455
456         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf3,zx_tmp_2d)
457         CALL histwrite(nid_tra,"psrf3",itap,zx_tmp_2d,
458     s                                     iim*(jjm+1),ndex)
459
460         CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf4,zx_tmp_2d)
461         CALL histwrite(nid_tra,"psrf4",itap,zx_tmp_2d,
462     s                                     iim*(jjm+1),ndex)
463777     continue
464c======================================================================
465c   Calcul de l'effet de la convection
466c======================================================================
467        print*,'Avant convection'
468      do it=1,nqmax
469         WRITE(itn,'(i1)') it
470c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
471      enddo
472
473      if (convection) then
474
475      print*,'Pas de temps dans phytrac : ',pdtphys
476      DO it=1, nqmax
477      CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
478     .            pplay, paprs, tr_seri(1,1,it), d_tr_cv)
479      DO k = 1, nlev
480      DO i = 1, klon
481         tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
482      ENDDO
483      ENDDO
484c      WRITE(itn,'(i1)') it
485c      CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
486      ENDDO
487c      print*,'apres nflxtr'
488
489
490      endif ! convection
491c        print*,'Apres convection'
492c      do it=1,nqmax
493c         WRITE(itn,'(i1)') it
494c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
495c      enddo
496
497c======================================================================
498c   Calcul de l'effet de la couche limite
499c======================================================================
500c       print *,'Avant couchelimite'
501c      do it=1,nqmax
502c         WRITE(itn,'(i1)') it
503c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
504c      enddo
505
506      if (couchelimite) then
507
508      DO k = 1, nlev
509      DO i = 1, klon
510         delp(i,k) = paprs(i,k)-paprs(i,k+1)
511      ENDDO
512      ENDDO
513
514C maf modif pour tenir compte du cas rnpb + traceur
515      DO it=1, nqmax
516c     print *,'it',it,clsol(it)
517      if (clsol(it)) then  ! couche limite avec quantite dans le sol calculee
518          CALL cltracrn(it, pdtphys, yu1, yv1,
519     e                    coefh,t_seri,ftsol,pctsrf,
520     e                    tr_seri(1,1,it),trs(1,it),
521     e                    paprs, pplay, delp,
522     e                    masktr(1,it),fshtr(1,it),hsoltr(it),
523     e                    tautr(it),vdeptr(it),
524     e                    xlat,
525     s                    d_tr_cl,d_trs)
526          DO k = 1, nlev
527            DO i = 1, klon
528              tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
529            ENDDO
530          ENDDO
531c
532c Traceur ds sol
533c
534          DO i = 1, klon
535            trs(i,it) = trs(i,it) + d_trs(i)
536          END DO
537C
538C maf provisoire suppression des prints
539C         WRITE(itn,'(i1)') it
540C         CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
541      else ! couche limite avec flux prescrit
542Cmaf provisoire source / traceur a creer
543        DO i=1, klon
544          source(i) = 0.0 ! pas de source, pour l'instant
545        ENDDO
546C
547          CALL cltrac(pdtphys, coefh,t_seri,
548     s               tr_seri(1,1,it), source,
549     e               paprs, pplay, delp,
550     s               d_tr )
551            DO k = 1, nlev
552               DO i = 1, klon
553                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
554               ENDDO
555            ENDDO
556Cmaf provisoire suppression des prints
557Cmaf          WRITE(itn,'(i1)') it
558cmaf          CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
559      endif
560      ENDDO
561c
562      endif ! couche limite
563
564c      print*,'Apres couchelimite'
565c      do it=1,nqmax
566c         WRITE(itn,'(i1)') it
567c        call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL  '//itn)
568c      enddo
569
570c======================================================================
571c   Calcul de l'effet du puits radioactif
572c======================================================================
573
574C MAF il faudrait faire une modification pour passer dans radiornpb
575c si radio=true mais pour l'instant radiornpb propre au cas rnpb
576      if(rnpb) then
577        call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
578C
579        DO it=1,nqmax
580            if(radio(it)) then
581            DO k = 1, nlev
582               DO i = 1, klon
583                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
584               ENDDO
585            ENDDO
586            WRITE(itn,'(i1)') it
587            CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
588            endif
589        ENDDO
590c
591      endif ! rnpb decroissance  radioactive
592C
593c======================================================================
594c   Calcul de l'effet de la precipitation
595c======================================================================
596
597      print*,'LESSIVAGE =',lessivage
598      IF (lessivage) THEN
599
600c     print*,'avant lessivage'
601
602       DO it = 1, nqmax
603           DO k = 1, nlev
604              DO i = 1, klon
605               d_tr_lessi_nucl(i,k,it) = 0.0
606               d_tr_lessi_impa(i,k,it) = 0.0
607               flestottr(i,k,it) = 0.0
608              ENDDO
609           ENDDO
610       ENDDO
611c
612c tendance des aerosols nuclees et impactes
613c
614       DO it = 1, nqmax
615         IF (aerosol(it)) THEN
616           DO k = 1, nlev
617              DO i = 1, klon
618               d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
619     s                  ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
620               d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
621     s                  ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
622              ENDDO
623           ENDDO
624         ENDIF
625       ENDDO
626c
627c Mises a jour des traceurs + calcul des flux de lessivage
628c Mise a jour due a l'impaction et a la nucleation
629c
630c      call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
631c      call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
632c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
633       DO it = 1, nqmax
634c         print*,'IT=',it,aerosol(it)
635         IF (aerosol(it)) THEN
636c           print*,'IT=',it,' On lessive'
637           DO k = 1, nlev
638              DO i = 1, klon
639               tr_seri(i,k,it)=tr_seri(i,k,it)
640     s         *frac_impa(i,k)*frac_nucl(i,k)
641              ENDDO
642           ENDDO
643         ENDIF
644       ENDDO
645c      call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
646c
647c Flux lessivage total
648c
649      DO it = 1, nqmax
650           DO k = 1, nlev
651            DO i = 1, klon
652               flestottr(i,k,it) = flestottr(i,k,it) -
653     s                   ( d_tr_lessi_nucl(i,k,it)   +
654     s                     d_tr_lessi_impa(i,k,it) ) *
655     s                   ( paprs(i,k)-paprs(i,k+1) ) /
656     s                   (RG * pdtphys)
657            ENDDO
658           ENDDO
659c
660Cmaf suppression provisoire
661Cmaf        WRITE(itn,'(i1)') it
662Cmaf    CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
663      ENDDO
664c
665c     print*,'apres lessivage'
666      ENDIF
667Cc
668      DO k = 1, nlev
669         DO i = 1, klon
670            fluxrn(i,k) = flestottr(i,k,1)
671            fluxpb(i,k) = flestottr(i,k,2)
672            rn(i,k) = tr_seri(i,k,1)
673            pb(i,k) = tr_seri(i,k,2)
674            pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
675            pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
676         ENDDO
677      ENDDO
678      itra=itra+1
679      ndex(1) = 0
680      DO it=1,nqmax
681      IF (it.LE.99) THEN
682       WRITE(str2,'(i2.2)') it
683C champs 3D
684       CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
685       CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
686     .                                   iim*(jjm+1)*klev,ndex)
687c      IF (lessivage) THEN
688c      CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
689c      CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
690c    .                                   iim*(jjm+1)*klev,ndex)
691c     ENDIF
692      ELSE
693         PRINT*, "Trop de traceurs"
694         CALL abort
695      ENDIF
696      ENDDO
697
698        goto 888
699        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pplay,zx_tmp_3d)
700        CALL histwrite(nid_tra,"pplay",itra,zx_tmp_3d,
701     .                  iim*(jjm+1)*klev,ndex)
702
703        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,t_seri,zx_tmp_3d)
704        CALL histwrite(nid_tra,"t",itra,zx_tmp_3d,
705     .                  iim*(jjm+1)*klev,ndex)
706        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfu,zx_tmp_3d)
707        CALL histwrite(nid_tra,"mfu",itra,zx_tmp_3d,
708     .                  iim*(jjm+1)*klev,ndex)
709        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfd,zx_tmp_3d)
710        CALL histwrite(nid_tra,"mfd",itra,zx_tmp_3d,
711     .                  iim*(jjm+1)*klev,ndex)
712        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_u,zx_tmp_3d)
713        CALL histwrite(nid_tra,"en_u",itra,zx_tmp_3d,
714     .                  iim*(jjm+1)*klev,ndex)
715        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_d,zx_tmp_3d)
716        CALL histwrite(nid_tra,"en_d",itra,zx_tmp_3d,
717     .                  iim*(jjm+1)*klev,ndex)
718        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_d,zx_tmp_3d)
719        CALL histwrite(nid_tra,"de_d",itra,zx_tmp_3d,
720     .                  iim*(jjm+1)*klev,ndex)
721        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_u,zx_tmp_3d)
722        CALL histwrite(nid_tra,"de_u",itra,zx_tmp_3d,
723     .                  iim*(jjm+1)*klev,ndex)
724        CALL gr_fi_ecrit(klev,klon,iim,jjm+1,coefh,zx_tmp_3d)
725        CALL histwrite(nid_tra,"coefh",itra,zx_tmp_3d,
726     .                  iim*(jjm+1)*klev,ndex)
727
728888     continue
729
730c       print*,'Sortie phytrac'
731c      do it=1,nqmax
732c         WRITE(itn,'(i1)') it
733c        call diagtracphy(tr_seri(:,:,it),paprs,'Fin Phys  '//itn)
734c      enddo
735
736      if (lafin) then
737         print*, 'c est la fin de la physique'
738         open (99,file='restarttrac',  form='formatted')
739         do i=1,klon
740             write(99,*) trs(i,1)
741         enddo
742         PRINT*, 'Ecriture du fichier restarttrac'
743         close(99)
744      else
745         print*, 'physique pas fini'
746      endif
747
748
749      RETURN
750      END
Note: See TracBrowser for help on using the repository browser.