source: LMDZ4/trunk/libf/phylmd/phytrac.F @ 541

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

Convergence avec la version d'Olivia Coindreau incluant:

  • le offline
  • les thermiques
  • mellor & yamada dans la couche limite

LF

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