source: trunk/libf/phyvenus/phytrac.older @ 4

Last change on this file since 4 was 3, checked in by slebonnois, 15 years ago

Creation de repertoires:

  • chantiers : pour communiquer sur nos projets de modifs
  • documentation : pour stocker les docs

Ajout de:

  • libf/phytitan : physique de Titan
  • libf/chimtitan: chimie de Titan
  • libf/phyvenus : physique de Venus
File size: 13.0 KB
RevLine 
[3]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 (nstep,
7     I                    gmtime,
8     I                    debutphy,
9     I                    lafin,
10     I                    nqmax,
11     I                    nlon,
12     I                    nlev,
13     I                    pdtphys,
14     I                    u,
15     I                    v,
16     I                    t_seri,
17     I                    paprs,
18     I                    pplay,
19     I                    coefh,
20     I                    yu1,
21     I                    yv1,
22     I                    ftsol,
23     I                    xlat,
24     I                    xlon,
25     I                    zlev,
26     I                    presnivs,
27     I                    pphis,
28     I                    pphi,
29     I                    albsol,
30     O                    tr_seri)
31
32      USE ioipsl
33
34      IMPLICIT none
35c======================================================================
36c Auteur(s) FH
37c Objet: Moniteur general des tendances traceurs
38c
39cAA Remarques en vrac:
40cAA--------------------
41cAA 1/ le call phytrac se fait avec nqmax
42c======================================================================
43#include "YOMCST.h"
44#include "dimensions.h"
45#include "dimphy.h"
46#include "clesphys.h" !///utile?
47#include "temps.h"
48#include "paramet.h"
49#include "control.h"
50#include "comgeomphy.h"
51#include "advtrac.h"
52c======================================================================
53
54c Arguments:
55c
56c   EN ENTREE:
57c   ==========
58c
59c   divers:
60c   -------
61c
62      integer nlon  ! nombre de points horizontaux
63      integer nlev  ! nombre de couches verticales
64      integer nqmax ! nombre de traceurs auxquels on applique la physique
65      integer nstep  ! appel physique
66      integer nseuil ! numero du premier traceur non CV
67c      integer julien !jour julien
68c      integer itop_con(nlon)
69c      integer ibas_con(nlon)
70      real gmtime
71      real pdtphys  ! pas d'integration pour la physique (seconde)
72      real t_seri(nlon,nlev) ! temperature
73      real tr_seri(nlon,nlev,nqmax) ! traceur 
74      real u(nlon,nlev)
75      real v(nlon,nlev)
76      real albsol(nlon)  ! albedo surface
77      real paprs(nlon,nlev+1)  ! pression pour chaque inter-couche (en Pa)
78      real ps(nlon)  ! pression surface
79      real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
80      real pphi(nlon,nlev) ! geopotentiel
81      real pphis(nlon)
82      REAL presnivs(nlev)
83      logical debutphy       ! le flag de l'initialisation de la physique
84      logical lafin          ! le flag de la fin de la physique
85c      REAL flxmass_w(nlon,nlev)
86
87cAA Rem : nqmax : nombre de vrais traceurs est defini dans dimphy.h
88c
89c   convection:
90c   -----------
91c
92      REAL pmfu(nlon,nlev)  ! flux de masse dans le panache montant
93      REAL pmfd(nlon,nlev)  ! flux de masse dans le panache descendant
94      REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
95
96c
97c   Couche limite:
98c   --------------
99c
100      REAL coefh(nlon,nlev) ! coeff melange CL
101      REAL yu1(nlon)        ! vents au premier niveau
102      REAL yv1(nlon)        ! vents au premier niveau
103      REAL xlat(nlon)       ! latitudes pour chaque point
104      REAL xlon(nlon)       ! longitudes pour chaque point
105      REAL zlev(nlon,nlev+1) ! altitude a chaque niveau (interface inferieure de la couche)
106cAA
107cAA Arguments necessaires pour les sources et puits de traceur:
108cAA ----------------
109cAA
110      real ftsol(nlon)  ! Temperature du sol (surf)(Kelvin)
111
112cAA ----------------------------
113cAA  VARIABLES LOCALES TRACEURS
114cAA ----------------------------
115cAA
116cAA Sources et puits des traceurs:
117cAA ------------------------------
118cAA
119      REAL source(klon)       ! a voir lorsque le flux est prescrit
120
121      CHARACTER*2 itn
122C maf ioipsl
123      CHARACTER*2 str2
124      INTEGER nhori, nvert
125      REAL zsto, zout, zjulian
126      INTEGER nid_tra
127      SAVE nid_tra
128      INTEGER nid_tra2,nid_tra3
129      SAVE nid_tra2,nid_tra3
130      INTEGER ndex(1)
131      INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
132      REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
133      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
134c
135      integer itau_w   ! pas de temps ecriture = nstep + itau_phy
136c
137
138C
139C Variables liees a l'ecriture de la bande histoire : phytrac.nc
140c
141c      INTEGER ecrit_tra
142c      SAVE ecrit_tra   
143      logical ok_sync
144      parameter (ok_sync = .true.)
145C
146C nature du traceur
147c
148      logical aerosol(nqmx)  ! Nature du traceur
149c                            ! aerosol(it) = true  => aerosol
150c                            ! aerosol(it) = false => gaz
151c                            ! nat_trac(it) = 1. aerosol
152      logical clsol(nqmx)    ! clsol(it) = true => CL sol calculee
153      save aerosol,clsol
154C
155C les traceurs
156C
157c===================
158c it--------indice de traceur
159c k,i---------indices long, vert
160c===================
161c Variables deja declarees dont on a besoin pour traceurs   
162c   k,i,it,tr_seri(klon,klev,nqmax),pplay(nlon,nlev),
163      real zprof(klon,klev)
164c      real pzero,gamma
165c      parameter (pzero=85000.)
166c      parameter (gamma=5000.)
167      real deltatr(klon,klev,nqmax) ! ecart au profil de ref zprof
168      real tau             ! temps de relaxation vers le profil
169c tau pourra dependre de it     
170      parameter (tau=1.) ! duree de vie du compose en secondes 3.e10
171c======================================================================
172c
173c Declaration des procedures appelees
174c
175c--modif convection tiedtke
176      INTEGER i, k, it
177      INTEGER iq, iiq
178      REAL delp(klon,klev)
179c--end modif
180c
181c Variables liees a l'ecriture de la bande histoire physique
182c
183c Variables locales pour effectuer les appels en serie
184c----------------------------------------------------
185c
186      REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
187      REAL d_tr_cl(klon,klev,nqmax) ! tendance de traceurs  couche limite
188      REAL d_tr_cv(klon,klev,nqmax) ! tendance de traceurs  conv pour chq traceur
189C
190      character*80 abort_message
191c
192c   Controles
193c-------------
194      logical first,couchelimite,convection
195      save first,couchelimite,convection
196c Olivia
197       data first,couchelimite,convection
198     s     /.true.,.false.,.false./
199
200c======================================================================
201         ps(:)=paprs(:,1)
202c---------
203c debutphy
204c---------
205         if (debutphy) then
206                 print*,"DEBUT PHYTRAC"
207C         
208c=============================================================
209c   Initialisation des traceurs
210c=============================================================
211c
212c Initialisation de la nature des traceurs
213c
214         DO it = 1, nqmax
215            aerosol(it) = .FALSE.  ! Tous les traceurs sont des gaz par defaut
216            clsol(it) = .FALSE.  ! Par defaut couche limite avec flux prescrit
217         ENDDO
218c
219c===========
220c definition de traceurs idealises
221c==========
222c
223c I) Declaration directe du traceur a altitude fixee
224c
225c a) traceur en carre OK
226c
227c         do i=1,klon
228c         tr_seri(i,:,1)=0.
229c        if ((xlat(i)>=0.).and.(xlat(i)<=-30.)) then
230c          if ((xlon(i)>=0.).and.(xlon(i)<=40.)) then
231c              tr_seri(i,10,1)=1.
232c          endif
233c        endif
234c      end do
235c
236c b) traceur a une bande en latitudeOK
237c
238c        do i=1,klon
239c         tr_seri(i,:,1)=0.
240c        if ((xlat(i)>=50.).and.(xlat(i)<=60.)) then
241c              tr_seri(i,25,1)=1.
242c        endif
243c      end do 
244c
245c c) traceur a plusieurs bandes en latitude OK
246c
247c         do i=1,klon
248c        tr_seri(i,:,2)=0.
249c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
250c             tr_seri(i,10,2)=1.
251c        endif
252c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
253c              tr_seri(i,10,2)=1.
254c        endif
255c
256c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
257c              tr_seri(i,10,2)=1.
258c        endif         
259c      end do
260c
261c d) traceur a une bande en altitude OK
262c
263c       do k=1,klev+1
264c         tr_seri(:,k,1)=0.
265c         if ((pplay(klon/2,k)>=1.e5).and.(pplay(klon/2,k)<=1.e6)) then
266c             tr_seri(:,k,1)=1.
267c         endif
268c       end do
269c
270c dbis) plusieurs traceurs a une bande en altitude OK
271c
272       do k=1,klev
273         tr_seri(:,k,1)=0.
274         if ((pplay(klon/2,k)>=9.e5).and.(pplay(klon/2,k)<=9.e6)) then
275             tr_seri(:,k,1)=1.
276         endif
277       end do
278
279       do k=1,klev
280         tr_seri(:,k,2)=0.
281         if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
282             tr_seri(:,k,2)=1.
283         endif
284       end do
285
286c      do k=1,klev
287c        tr_seri(:,k,3)=0.
288c        if ((pplay(klon/2,k)>=5.e1).and.(pplay(klon/2,k)<=5.e2)) then
289c            tr_seri(:,k,3)=1.
290c        endif
291c      end do
292c       
293c e) plusieurs couches verticales de traceurs, a plusieurs bandes en latitude???
294c       
295c au sol
296c         do i=1,klon
297c        tr_seri(i,:,1)=0.
298c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
299c             tr_seri(i,5,1)=1.
300c        endif
301c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
302c              tr_seri(i,5,1)=1.
303c        endif
304c
305c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
306c              tr_seri(i,5,1)=1.
307c        endif         
308c      end do
309c
310c         do i=1,klon
311c        tr_seri(i,:,2)=0.
312c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
313c             tr_seri(i,10,2)=1.
314c        endif
315c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
316c              tr_seri(i,10,2)=1.
317c      endif
318c
319c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
320c              tr_seri(i,10,2)=1.
321c        endif         
322c      end do
323c
324c         do i=1,klon
325c        tr_seri(i,:,3)=0.
326c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
327c             tr_seri(i,30,3)=1.
328c        endif
329c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
330c              tr_seri(i,30,3)=1.
331c        endif
332c
333c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
334c              tr_seri(i,30,3)=1.
335c        endif         
336c      end do
337c       
338c        do i=1,klon
339c        tr_seri(i,:,4)=0.
340c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
341c             tr_seri(i,45,4)=1.
342c        endif
343c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
344c              tr_seri(i,45,4)=1.
345c        endif
346c
347c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
348c              tr_seri(i,45,4)=1.
349c        endif         
350c      end do
351c     
352c II) Declaration d'un profil vertical de traceur OK
353c
354c      do i=1,klon
355c       do k=1,klev
356c         zprof(i,k)=0.5*(1.+tanh((pplay(i,k)-pzero)/gamma))
357c--------x=alog10P+b
358c        zprof(i,k)=-8.79e-6*(log10(pplay(i,k)))+7.6e-5
359c        zprof(i,k)=8.79e-6*(log10(pplay(i,k)))-7.6e-5
360c--------alog10x=alog10P+b       
361c        zprof(i,k)=-0.31*(log10(pplay(i,k)))+5.2
362c        zprof(i,k)=-0.31*(log10(pplay(i,k)))-2.68
363c       enddo 
364c      enddo
365c
366c      do i=1,klon
367c       do k=1,klev
368c         zprof(i,k)=0.5*(1.+tanh((pplay(i,k)-pzero)/gamma))
369c       enddo
370c      enddo
371c-------------
372c fin debutphy
373c-------------
374      ENDIF  ! fin debutphy
375
376c======================================================================
377c Rappel vers un profil
378c======================================================================
379c         do it=1,nqmax-2
380c           do k=1,klev
381c            do i=1,klon
382c             deltatr(i,k,1) = 0.
383c             tr_seri(i,k,1) = 0.
384c             deltatr(i,k,1) = (-tr_seri(i,k,1)+zprof(i,k))/tau
385c             tr_seri(i,k,1) =  tr_seri(i,k,1) + deltatr(i,k,1)
386c            enddo
387c            enddo
388c         enddo
389
390c======================================================================
391c   Calcul de l'effet de la couche limite
392c======================================================================
393
394      if (couchelimite) then
395
396      DO k = 1, klev
397      DO i = 1, klon
398         delp(i,k) = paprs(i,k)-paprs(i,k+1)
399      ENDDO
400      ENDDO
401
402C maf modif pour tenir compte du cas rnpb + traceur
403      DO it=1, nqmax
404
405c     print *,'it',it,clsol(it)
406      if (clsol(it)) then  ! couche limite avec quantite dans le sol calculee
407
408      else ! couche limite avec flux prescrit
409
410Cmaf provisoire source / traceur a creer
411        DO i=1, klon
412          source(i) = 0.0 ! pas de source, pour l'instant
413        ENDDO
414C
415
416          CALL cltrac(pdtphys, coefh,t_seri,
417     s               tr_seri(1,1,it), source,
418     e               paprs, pplay, delp,
419     s               d_tr_cl(1,1,it))
420            DO k = 1, nlev
421               DO i = 1, klon
422                  tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
423               ENDDO
424            ENDDO
425Cmaf          WRITE(itn,'(i1)') it
426cmaf          CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
427      endif
428      ENDDO
429c
430      endif ! couche limite
431
432c=============================================================
433
434      RETURN
435      END
Note: See TracBrowser for help on using the repository browser.