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

Last change on this file since 662 was 662, checked in by lmdzadmin, 19 years ago

Suppression du logical "thermiques" remplace par iflag_thermals
MAF

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