source: LMDZ.3.3/trunk/libf/phylmd/phytrac.F @ 204

Last change on this file since 204 was 204, checked in by lmdz, 24 years ago

Debogage du guidage et de la version debranchee et abandon de la version
debranchee non-netcdf FH/MAF
LF

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