source: LMDZ4/branches/IPSL-CM4_IPCC_branch/libf/phylmd/phytrac.F @ 5427

Last change on this file since 5427 was 524, checked in by lmdzadmin, 21 years ago

Initial revision

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