source: LMDZ4/branches/V3_test/libf/phylmd/phytrac.F @ 4487

Last change on this file since 4487 was 735, checked in by Laurent Fairhead, 18 years ago

Modifications version parallele
YM/LF

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