source: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_emiss.F @ 3722

Last change on this file since 3722 was 3714, checked in by slebonnois, 8 months ago

SL: small adjustment in radlwsw.F
File size: 8.8 KB
Line 
1!
2! $Header: /home/cvsroot/LMDZ4/libf/phylmd/phytrac.F,v 1.16 2006/03/24 15:06:23 lmdzadmin Exp $
3!
4c
5c
6      SUBROUTINE phytrac_emiss (debutphy,
7     I                    lafin,
8     I                    nqmax,
9     I                    nlon,
10     I                    nlev,
11     I                    pdtphys,
12     I                    paprs,
13     I                    xlat,xlon,
14     O                    tr_seri)
15     
16     
17
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
25c
26c SL: Janvier 2014
27c Version developed for surface emission
28c Maybe could be used just to compute the 'source' variable from physiq
29c
30c======================================================================
31      use infotrac_phy, only: tname, nqtot
32#ifdef CPP_XIOS     
33      use xios_output_mod, only:  send_xios_field
34#endif
35      use dimphy
36      USE geometry_mod, only: cell_area
37      USE chemparam_mod,only:M_tr,type_tr
38      USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
39      IMPLICIT none
40#include "YOMCST.h"
41#include "clesphys.h"
42     
43c======================================================================
44
45c Arguments:
46
47c   EN ENTREE:
48c   ==========
49
50      logical debutphy       ! le flag de l'initialisation de la physique
51      logical lafin          ! le flag de la fin de la physique
52      integer nqmax ! nombre de traceurs auxquels on applique la physique
53      integer nlon  ! nombre de points horizontaux
54      integer nlev  ! nombre de couches verticales
55      real pdtphys  ! pas d'integration pour la physique (seconde)
56      real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
57      REAL xlat(nlon)       ! latitudes pour chaque point
58      REAL xlon(nlon)       ! longitudes pour chaque point
59
60c   EN ENTREE/SORTIE:
61c   =================
62
63      real tr_seri(nlon,nlev,nqmax) ! traceur 
64
65cAA ----------------------------
66cAA  VARIABLES LOCALES TRACEURS
67cAA ----------------------------
68
69c pour emission type effusion
70      real :: deltatr(klon,klev,nqtot)
71
72c several identical tracers emitted at different places, with different
73c fluxes, at the same time :
74c same_tracer=.true. and various number of areas and fluxes
75c      logical,parameter :: same_tracer=.true.
76c      integer,parameter :: nblat=3,nblon=3,nbaire=3,nbflux=2,maxcell=16
77c Otherwise, emission of one of the tracers only, in one area only
78      logical,parameter :: same_tracer=.false.
79      integer,parameter :: nblat=1,nblon=1,nbaire=1,nbflux=1
80      integer,save :: iq_co2,iq_n2
81
82      integer,parameter :: maxcell=10000
83      integer,parameter :: Nheight=3 ! layer emission (150m)
84      integer,save :: Ncell(nbaire)
85      real,save :: flux_surface_co2(nbflux) ! flux de CO2 emis (kg/m2/s)
86      real :: flux(nlon,nqtot)
87      real,save :: lat_zone(nblat),lon_zone(nblon)
88      integer,save :: ig_zone(nblat,nblon,nbaire,nbflux,maxcell)
89      integer,save :: numcell(nblat,nblon,nbaire,nbflux)
90
91      INTEGER i, k, it
92      integer ilat,ilon,iaire,iflux,ipos
93      real    deltalat,deltalon
94     
95
96
97c======================================================================
98
99c EMISSIONS TRACEURS
100
101c---------
102c debutphy
103c---------     
104      if (debutphy) then
105
106        print*,"DEBUT PHYTRAC"
107        print*,"PHYTRAC: EMISSION"
108
109        ALLOCATE(M_tr(nqtot))
110        ALLOCATE(type_tr(nqtot))
111        do i=1,nqtot
112           if (tname(i)(1:3).eq.'co2') then   ! CO2
113               M_tr(i)=44.
114               iq_co2 =i
115               type_tr(i)=1
116           endif
117           if (tname(i)(1:2).eq.'n2')  then   ! N2
118               M_tr(i)=28.
119               iq_n2 =i
120               type_tr(i)=1
121           endif
122        enddo
123
124C=========================================================================
125c Caracteristiques des traceurs emis:
126C=========================================================================
127
128
129c nombre total de traceur
130         if (nblat*nblon*nbaire*nbflux+1 .gt. nqtot) then
131            print*, nblat*nblon*nbaire*nbflux+1, nqtot
132            write(*,*) "Attention, pas assez de traceurs"
133            write(*,*) "le dernier sera bien le dernier"
134         endif
135                 
136c flux de CO2 (kg/s/m2)
137c 4.5e-7 correspond a 1% (en masse) d une couche de 2E6 Pa en 500 jV
138c donc avec 5E-7, on vide le N2 en 1000 jV
139c         flux_surface_co2(1) = 5.E-5
140         flux_surface_co2(1) = 0.
141
142c nombre de cellule pour le cote du carre d'aire
143         Ncell(1)= 9999
144
145c localisation zone emission
146         lat_zone(1) = 0.
147         lon_zone(1) = 0.
148 
149         if ((nbp_lon*nbp_lat)==1) then ! running a 1D simulation
150           deltalat=180.
151           deltalon=360.
152         else
153           deltalat = 180./(nbp_lat-1)
154           deltalon = 360./nbp_lon
155         endif
156
157         numcell(:,:,:,:)=0
158         ig_zone(:,:,:,:,:)=0
159         do i=1,nlon
160          do ilat=1,nblat
161           do ilon=1,nblon
162            do iaire=1,nbaire
163
164            if ( 
165     &       (Ncell(iaire).ne.9999)
166c emission on some areas only
167     &      .and.((xlat(i).ge.lat_zone(ilat))
168     &      .and.((xlat(i)-Ncell(iaire)*deltalat)
169     &      .lt.lat_zone(ilat))
170     &      .and.(xlon(i).le.lon_zone(ilon))
171     &      .and.((xlon(i)+Ncell(iaire)*deltalon)
172     &      .gt.lon_zone(ilon))) ) then
173
174              do iflux=1,nbflux
175         numcell(ilat,ilon,iaire,iflux)=numcell(ilat,ilon,iaire,iflux)+1
176        ig_zone(ilat,ilon,iaire,iflux,numcell(ilat,ilon,iaire,iflux))= i
177             print*,"Lat,lon,naire,nflux,nlon=",ilat,ilon,iaire,iflux
178     &        ,i," OK"
179              end do
180
181             end if
182
183            end do
184           end do
185          end do
186         end do
187
188c Reinit des traceurs si necessaire
189         if (reinit_trac) then
190          tr_seri(:,:,:)=0.
191
192           do i=1,klon
193            do k=1,klev
194c Profile for VeGa2
195c              tr_seri(i,k,iq_n2) = M_tr(i_n2)/43.44*max(min(0.035,
196c     &           0.035*(1.-log(paprs(i,k)/6.e6)/log(9.e6/6.e6))),0.)
197c              tr_seri(i,k,iq_co2)= 1. - tr_seri(i,k,iq_n2)
198
199c Uniform initialization, yN2=3.5%, regular mean mol mass = 43.44 g/mol
200              tr_seri(i,k,iq_n2) = 0.035*M_tr(iq_n2)/43.44
201              tr_seri(i,k,iq_co2)= 1. - tr_seri(i,k,iq_n2)
202            end do
203           end do
204         endif
205 
206C=========================================================================
207C=========================================================================
208      ENDIF  ! fin debutphy
209c-------------
210c fin debutphy
211c-------------
212
213c======================================================================
214c Emission continue d'un traceur
215c necessite raz_date=1 dans run.def
216c et reinit_trac=y
217c======================================================================
218       deltatr(:,:,:) = 0.
219       flux(:,:)=0.
220
221        if (same_tracer) then
222c emet les traceurs qui sont presents sur la grille
223        do ilat  = 1,nblat
224        do ilon  = 1,nblon
225        do iaire = 1,nbaire
226        do iflux = 1,nbflux
227     
228              it=min( (ilat-1)*nblon*nbflux*nbaire+(iaire-1)*nbflux
229     &         +(ilon-1)*nbaire*nbflux+iflux , nqtot )   
230 
231         if (Ncell(iaire).ne.9999) then
232c column injection
233          do ipos=1,maxcell
234 
235           if (ig_zone(ilat,ilon,iaire,iflux,ipos).ne.0) then
236
237             i=ig_zone(ilat,ilon,iaire,iflux,ipos)         
238             flux(i,it)=flux_surface_co2(iflux) 
239             
240            do k=1,Nheight
241             deltatr(i,k,it) = flux_surface_co2(iflux) ! kg/s/m2
242     $         *RG/(paprs(i,1)-paprs(i,Nheight+1))    ! / (kg/m2) (masse colonne)
243     
244             tr_seri(i,k,it) = tr_seri(i,k,it)+deltatr(i,k,it)*pdtphys
245            end do
246
247           end if
248          end do
249         endif
250        end do
251        end do
252        end do
253        end do
254
255        else  ! same_tracer=.false.
256
257c column injection !! with constant mass !!
258          do i=1,nlon
259 
260             flux(i,iq_co2)=flux_surface_co2(1) 
261             
262            do k=1,Nheight
263          deltatr(i,k,iq_co2) = min((1.-tr_seri(i,k,iq_co2))/pdtphys,
264     $                                flux(i,iq_co2)) ! kg/s/m2
265     $         *RG/(paprs(i,1)-paprs(i,Nheight+1))    ! / (kg/m2) (masse colonne)
266             deltatr(i,k,iq_n2)  = -deltatr(i,k,iq_co2)
267             tr_seri(i,k,iq_co2) = tr_seri(i,k,iq_co2)
268     $                            +deltatr(i,k,iq_co2)*pdtphys
269             tr_seri(i,k,iq_n2)  = tr_seri(i,k,iq_n2)
270     $                            +deltatr(i,k,iq_n2)*pdtphys
271            end do
272
273          end do
274
275        endif ! same_tracer
276 
277       
278c======================================================================
279c======================================================================
280
281#ifdef CPP_XIOS     
282!      do it=1,nqtot
283!      CALL  send_xios_field("flux_"//tname(it),
284!    &                     flux(:,it))
285!     end do
286#endif
287
288      RETURN
289      END
Note: See TracBrowser for help on using the repository browser.