source: LMDZ.3.3/branches/rel-LF/libf/phylmd/phytrac.F @ 436

Last change on this file since 436 was 353, checked in by lmdzadmin, 23 years ago

2 changements pour les fichiers histoire:

  • utilisation de l'entree "rectilineaire" de IOIPSL pour ne plus avoir

a

lancer ncregular a chaque fois

  • le calendrier des fichiers histoire est maintenant base sur la date d'initialisation de la simulation plutot que sur la date de depart du

job

en cours

LF

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