source: trunk/libf/phyvenus/phytrac.old @ 86

Last change on this file since 86 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: 24.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
116
117      CHARACTER*2 itn
118C maf ioipsl
119      CHARACTER*2 str2
120      INTEGER nhori, nvert
121      REAL zsto, zout, zjulian
122      INTEGER nid_tra
123      SAVE nid_tra
124      INTEGER nid_tra2,nid_tra3
125      SAVE nid_tra2,nid_tra3
126      INTEGER ndex(1)
127      INTEGER ndex2d(iim*(jjm+1)),ndex3d(iim*(jjm+1)*klev)
128      REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
129      REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
130c
131      integer itau_w   ! pas de temps ecriture = nstep + itau_phy
132c
133
134C
135C Variables liees a l'ecriture de la bande histoire : phytrac.nc
136c
137c      INTEGER ecrit_tra
138c      SAVE ecrit_tra   
139      logical ok_sync
140      parameter (ok_sync = .true.)
141C
142C les traceurs
143C
144c===================
145c it--------indice de traceur
146c k,i---------indices long, vert
147c===================
148c Variables deja declarees dont on a besoin pour traceurs   
149c   k,i,it,tr_seri(klon,klev,nqmax),pplay(nlon,nlev),
150      real zprof(klon,klev,nqmx)
151c      real pzero,gamma
152c      parameter (pzero=85000.)
153c      parameter (gamma=5000.)
154      REAL alpha
155      real deltatr(klon,klev,nqmx) ! ecart au profil de ref zprof
156      real tau(klev,nqmx)          ! temps de relaxation vers le profil (s)
157      save zprof,tau
158c======================================================================
159c
160c Declaration des procedures appelees
161c
162c--modif convection tiedtke
163      INTEGER i, k, it
164      INTEGER iq, iiq
165      REAL delp(klon,klev)
166c--end modif
167c
168c Variables liees a l'ecriture de la bande histoire physique
169c
170c Variables locales pour effectuer les appels en serie
171c----------------------------------------------------
172c
173      REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
174      REAL d_tr_cl(klon,klev,nqmax) ! tendance de traceurs  couche limite
175      REAL d_tr_cv(klon,klev,nqmax) ! tendance de traceurs  conv pour chq traceur
176C
177      character*80 abort_message
178c
179c   Controles
180c-------------
181      logical first,couchelimite,convection
182      save first,couchelimite,convection
183c Olivia
184       data first,couchelimite,convection
185     s     /.true.,.false.,.false./
186
187c======================================================================
188         ps(:)=paprs(:,1)
189c---------
190c debutphy
191c---------
192         if (debutphy) then
193                 print*,"DEBUT PHYTRAC"
194C         
195c=============================================================
196c   Initialisation des traceurs
197c=============================================================
198c
199c===========
200c definition de traceurs idealises
201c==========
202c
203c I) Declaration directe du traceur a altitude fixee
204c
205c a) traceur en carre OK
206c
207c         do i=1,klon
208c         tr_seri(i,:,1)=0.
209c        if ((xlat(i)>=0.).and.(xlat(i)<=-30.)) then
210c          if ((xlon(i)>=0.).and.(xlon(i)<=40.)) then
211c              tr_seri(i,10,1)=1.
212c          endif
213c        endif
214c      end do
215c
216c a bis) 2 traceurs en carre lat/alt, uniforme en longitude OK
217c
218C entre 45-55 km
219c
220c         do i=1,klon
221c         do k=1,klev+1
222cc         tr_seri(i,k,1)=0.
223c           if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
224c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
225c           if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=4.e5)) then
226c               tr_seri(i,k,1)=1.
227c           endif
228c           endif
229c           endif
230c         else
231c            tr_seri(i,k,1)=0.
232c         end do
233c         end do
234cc
235c         do i=1,klon
236c         do k=1,klev+1
237cc         tr_seri(i,k,2)=0.
238c           if ((xlat(i)>=-60.).and.(xlat(i)<=-80.)) then
239c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
240c           if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=4.e5)) then
241c               tr_seri(i,k,2)=1.
242c           endif
243c           endif
244c           endif
245c         else
246c            tr_seri(i,k,2)=0.
247c         end do
248c         end do
249cc
250c         do i=1,klon
251c         do k=1,klev+1
252cc         tr_seri(i,k,3)=0.
253c           if ((xlat(i)>=40.).and.(xlat(i)<=60.)) then
254c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
255c           if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=4.e5)) then
256c               tr_seri(i,k,3)=1.
257c           endif
258c           endif
259c           endif
260c         else
261c            tr_seri(i,k,3)=0.
262c         end do
263c         end do
264cc
265c         do i=1,klon
266c         do k=1,klev+1
267cc         tr_seri(i,k,4)=0.
268c           if ((xlat(i)>=-40.).and.(xlat(i)<=-60.)) then
269c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
270c           if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=4.e5)) then
271c               tr_seri(i,k,4)=1.
272c           endif
273c           endif
274c           endif
275c         else
276c            tr_seri(i,k,4)=0.
277c         end do
278c         end do
279cc
280c         do i=1,klon
281c         do k=1,klev+1
282cc         tr_seri(i,k,5)=0.
283c           if ((xlat(i)>=-20.).and.(xlat(i)<=20.)) then
284c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
285c           if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=4.e5)) then
286c              tr_seri(i,k,5)=1.
287c           endif
288c           endif
289c           endif
290c         else
291c            tr_seri(i,k,5)=0.
292c         end do
293c         end do
294c
295c entre 35-45 km
296c
297c         do i=1,klon
298c         do k=1,klev+1
299cc         tr_seri(i,k,6)=0.
300c           if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
301c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
302c           if ((pplay(klon/2,k)>=4.e5).and.(pplay(klon/2,k)<=8.e6)) then
303c               tr_seri(i,k,6)=1.
304c           endif
305c           endif
306c           endif
307c         else
308c            tr_seri(i,k,6)=0.
309c         end do
310c         end do
311c
312c         do i=1,klon
313c         do k=1,klev+1
314cc         tr_seri(i,k,7)=0.
315c           if ((xlat(i)>=-60.).and.(xlat(i)<=-80.)) then
316c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
317c           if ((pplay(klon/2,k)>=4.e5).and.(pplay(klon/2,k)<=8.e6)) then
318c               tr_seri(i,k,7)=1.
319c           endif
320c           endif
321c           endif
322c         else
323c            tr_seri(i,k,7)=0.
324c         end do
325c         end do
326c
327C entre 50-60 km
328c
329c         do i=1,klon
330c         do k=1,klev+1
331cc         tr_seri(i,k,8)=0.
332c           if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
333c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
334c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
335c               tr_seri(i,k,8)=1.
336c           endif
337c           endif
338c          endif
339c         else
340c            tr_seri(i,k,8)=0.
341c         end do
342c         end do
343c
344c         do i=1,klon
345c         do k=1,klev+1
346cc         tr_seri(i,k,9)=0.
347c           if ((xlat(i)>=-80.).and.(xlat(i)<=-60.)) then
348c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
349c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
350c               tr_seri(i,k,9)=1.
351c           endif
352c           endif
353c           endif
354c         else
355c            tr_seri(i,k,9)=0.
356c         end do
357c         end do
358c
359c         do i=1,klon
360c         do k=1,klev+1
361cc         tr_seri(i,k,10)=0.
362c           if ((xlat(i)>=40.).and.(xlat(i)<=60.)) then
363c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
364c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
365c               tr_seri(i,k,10)=1.
366c           endif
367c           endif
368c           endif
369c         else
370c            tr_seri(i,k,10)=0.
371c         end do
372c         end do
373c
374c         do i=1,klon
375c         do k=1,klev+1
376cc         tr_seri(i,k,11)=0.
377c           if ((xlat(i)>=-60.).and.(xlat(i)<=-40.)) then
378c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
379c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
380c               tr_seri(i,k,11)=1.
381c           endif
382c           endif
383c           endif
384c         else
385c            tr_seri(i,k,11)=0.
386c         end do
387c         end do
388c
389c         do i=1,klon
390c         do k=1,klev+1
391cc         tr_seri(i,k,12)=0.
392c           if ((xlat(i)>=-20.).and.(xlat(i)<=20.)) then
393c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
394c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e5)) then
395c               tr_seri(i,k,12)=1.
396c           endif
397c           endif
398c           endif
399c         else
400c            tr_seri(i,k,12)=0.
401c         end do
402c         end do
403c
404c entre 20-30 km
405c
406c         do i=1,klon
407c         do k=1,klev+1
408cc         tr_seri(i,k,13)=0.
409c           if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
410c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
411c           if ((pplay(klon/2,k)>=1.e6).and.(pplay(klon/2,k)<=2.e6)) then
412c               tr_seri(i,k,13)=1.
413c           endif
414c           endif
415c           endif
416c         else
417c            tr_seri(i,k,13)=0.
418c         end do
419c         end do
420c
421c         do i=1,klon
422c         do k=1,klev+1
423cc         tr_seri(i,k,14)=0.
424c           if ((xlat(i)>=-80.).and.(xlat(i)<=-60.)) then
425c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
426c           if ((pplay(klon/2,k)>=1.e6).and.(pplay(klon/2,k)<=2.e6)) then
427c               tr_seri(i,k,14)=1.
428c           endif
429c           endif
430c           endif
431c         else
432c            tr_seri(i,k,14)=0.
433c         end do
434c         end do
435c
436c         do i=1,klon
437c         do k=1,klev+1
438cc         tr_seri(i,k,15)=0.
439c           if ((xlat(i)>=-20.).and.(xlat(i)<=20.)) then
440c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
441c           if ((pplay(klon/2,k)>=1.e6).and.(pplay(klon/2,k)<=2.e6)) then
442c               tr_seri(i,k,15)=1.
443c           endif
444c           endif
445c           endif
446c         else
447c            tr_seri(i,k,15)=0.
448c         end do
449c         end do
450c
451c entre 55-65 km
452c
453c         do i=1,klon
454c         do k=1,klev+1
455cc         tr_seri(i,k,16)=0.
456c           if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
457c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
458c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=5.e4)) then
459c               tr_seri(i,k,16)=1.
460c           endif
461c           endif
462c           endif
463c           endif
464c         else
465c            tr_seri(i,k,16)=0.
466c         end do
467c         end do
468c
469c         do i=1,klon
470c         do k=1,klev+1
471cc         tr_seri(i,k,17)=0.
472c           if ((xlat(i)>=-80.).and.(xlat(i)<=-60.)) then
473c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
474c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=5.e4)) then
475c               tr_seri(i,k,17)=1.
476c          endif
477c          endif
478c           endif
479c           endif
480c         else
481c            tr_seri(i,k,17)=0.
482c         end do
483c         end do
484c
485c         do i=1,klon
486c         do k=1,klev+1
487cc         tr_seri(i,k,18)=0.
488c           if ((xlat(i)>=-20.).and.(xlat(i)<=20.)) then
489c           if ((xlon(i)>=-180.).and.(xlon(i)<=180.)) then
490c           if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=5.e4)) then
491c               tr_seri(i,k,18)=1.
492c           endif
493c           endif
494c           endif
495c           endif
496c         else
497c            tr_seri(i,k,18)=0.
498c         end do
499c         end do
500c
501c b) traceur a une bande en latitudeOK
502c
503c a 65km
504c
505c        do i=1,klon
506c         tr_seri(i,:,1)=0.
507c        if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
508c              tr_seri(i,20,1)=1.
509c        endif
510c      end do 
511c
512c        do i=1,klon
513c         tr_seri(i,:,2)=0.
514c        if ((xlat(i)>=40.).and.(xlat(i)<=60.)) then
515c              tr_seri(i,20,2)=1.
516c        endif
517c      end do 
518c
519c        do i=1,klon
520c         tr_seri(i,:,3)=0.
521c        if ((xlat(i)>=20.).and.(xlat(i)<=40.)) then
522c              tr_seri(i,20,3)=1.
523c        endif
524c      end do 
525c
526c        do i=1,klon
527c         tr_seri(i,:,4)=0.
528c        if ((xlat(i)>=0.).and.(xlat(i)<=20.)) then
529c              tr_seri(i,20,4)=1.
530c        endif
531c      end do 
532c
533c        do i=1,klon
534c         tr_seri(i,:,5)=0.
535c        if ((xlat(i)>=-20.).and.(xlat(i)<=0.)) then
536c              tr_seri(i,20,5)=1.
537c        endif
538c      end do 
539c
540c        do i=1,klon
541c         tr_seri(i,:,6)=0.
542c        if ((xlat(i)>=-40.).and.(xlat(i)<=-20.)) then
543c              tr_seri(i,20,6)=1.
544c        endif
545c      end do 
546c
547c        do i=1,klon
548c         tr_seri(i,:,7)=0.
549c        if ((xlat(i)>=-60.).and.(xlat(i)<=-40.)) then
550c              tr_seri(i,20,7)=1.
551c        endif
552c      end do 
553c
554c        do i=1,klon
555c         tr_seri(i,:,8)=0.
556c        if ((xlat(i)>=-80.).and.(xlat(i)<=-60.)) then
557c              tr_seri(i,20,8)=1.
558c        endif
559c      end do 
560c
561c a 50km
562c
563c        do i=1,klon
564c        tr_seri(i,:,1)=0.
565c        if ((xlat(i)>=40.).and.(xlat(i)<=60.)) then
566c              tr_seri(i,27,1)=1.
567c        endif
568c      end do 
569c
570c        do i=1,klon
571c         tr_seri(i,:,2)=0.
572c        if ((xlat(i)>=60.).and.(xlat(i)<=80.)) then
573c              tr_seri(i,27,2)=1.
574c        endif
575c      end do 
576c
577c        do i=1,klon
578c         tr_seri(i,:,3)=0.
579c        if ((xlat(i)>=20.).and.(xlat(i)<=40.)) then
580c              tr_seri(i,27,3)=1.
581c        endif
582c      end do 
583c
584c        do i=1,klon
585c         tr_seri(i,:4)=0.
586c        if ((xlat(i)>=0.).and.(xlat(i)<=20.)) then
587c              tr_seri(i,27,4)=1.
588c       endif
589c      end do 
590c
591c        do i=1,klon
592c         tr_seri(i,:,5)=0.
593c        if ((xlat(i)>=-20.).and.(xlat(i)<=0.)) then
594c              tr_seri(i,27,5)=1.
595c        endif
596c      end do 
597c
598c        do i=1,klon
599c         tr_seri(i,:,6)=0.
600c        if ((xlat(i)>=-40.).and.(xlat(i)<=-20.)) then
601c              tr_seri(i,27,6)=1.
602c        endif
603c      end do 
604c
605c        do i=1,klon
606c         tr_seri(i,:,7)=0.
607c        if ((xlat(i)>=-60.).and.(xlat(i)<=-40.)) then
608c              tr_seri(i,27,7)=1.
609c        endif
610c      end do 
611c
612c        do i=1,klon
613c         tr_seri(i,:,8)=0.
614c        if ((xlat(i)>=-80.).and.(xlat(i)<=-60.)) then
615c              tr_seri(i,27,8)=1.
616c        endif
617c      end do 
618c
619c c) traceur a plusieurs bandes en latitude OK
620c
621c         do i=1,klon
622c        tr_seri(i,:,2)=0.
623c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
624c             tr_seri(i,10,2)=1.
625c        endif
626c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
627c              tr_seri(i,10,2)=1.
628c        endif
629c
630c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
631c              tr_seri(i,10,2)=1.
632c        endif         
633c      end do
634c
635c d) traceur a une bande en altitude OK
636c
637c       do k=1,klev+1
638c         tr_seri(:,k,1)=0.
639c         if ((pplay(klon/2,k)>=1.e5).and.(pplay(klon/2,k)<=1.e6)) then
640c             tr_seri(:,k,1)=1.
641c         endif
642c       end do
643c
644c dbis) plusieurs traceurs a une bande en altitude OK
645c
646c bande tres basse tropo
647c      do k=1,klev
648c        tr_seri(:,k,1)=0.
649c        if ((pplay(klon/2,k)>=5.e5).and.(pplay(klon/2,k)<=5.e6)) then
650c            tr_seri(:,k,1)=1.
651c        endif
652c      end do
653c bande dans les nuages et un peu en-dessous   
654c      do k=1,klev
655c         tr_seri(:,k,2)=0.
656c         if ((pplay(klon/2,k)>=5.e4).and.(pplay(klon/2,k)<=5.e5)) then
657c             tr_seri(:,k,2)=1.
658c         endif
659c       end do
660cune grosse epaisseur: inclue toute la circulation meridienne
661c      do k=1,klev
662c        tr_seri(:,k,1)=0.
663c        if ((pplay(klon/2,k)>=1.e4).and.(pplay(klon/2,k)<=1.e6)) then
664c            tr_seri(:,k,1)=1.
665c        endif
666c      end do
667cune grosse epaisseur: inclue la mesosphere
668c      do k=1,klev
669c         tr_seri(:,k,2)=0.
670c         if ((pplay(klon/2,k)>=2.e2).and.(pplay(klon/2,k)<=1.e4)) then
671c             tr_seri(:,k,2)=1.
672c         endif
673c       end do
674c
675c      do k=1,klev
676c        tr_seri(:,k,3)=0.
677c        if ((pplay(klon/2,k)>=5.e1).and.(pplay(klon/2,k)<=5.e2)) then
678c            tr_seri(:,k,3)=1.
679c        endif
680c      end do
681c
682c e) plusieurs couches verticales de traceurs, a plusieurs bandes en latitude???
683c       
684c au sol
685c         do i=1,klon
686c        tr_seri(i,:,1)=0.
687c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
688c             tr_seri(i,5,1)=1.
689c        endif
690c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
691c              tr_seri(i,5,1)=1.
692c        endif
693c
694c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
695c              tr_seri(i,5,1)=1.
696c        endif         
697c      end do
698c
699c         do i=1,klon
700c        tr_seri(i,:,2)=0.
701c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
702c             tr_seri(i,10,2)=1.
703c        endif
704c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
705c              tr_seri(i,10,2)=1.
706c      endif
707c
708c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
709c              tr_seri(i,10,2)=1.
710c        endif         
711c      end do
712c
713c         do i=1,klon
714c        tr_seri(i,:,3)=0.
715c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
716c             tr_seri(i,30,3)=1.
717c        endif
718c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
719c              tr_seri(i,30,3)=1.
720c        endif
721c
722c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
723c              tr_seri(i,30,3)=1.
724c        endif         
725c      end do
726c       
727c        do i=1,klon
728c        tr_seri(i,:,4)=0.
729c        if ((xlat(i)>=50.).and.(xlat(i)<=70.)) then
730c             tr_seri(i,45,4)=1.
731c        endif
732c        if ((xlat(i)>=-10.).and.(xlat(i)<=10.)) then
733c              tr_seri(i,45,4)=1.
734c        endif
735c
736c        if ((xlat(i)>=-70.).and.(xlat(i)<=-50.)) then
737c              tr_seri(i,45,4)=1.
738c        endif         
739c      end do
740c     
741
742c II) Declaration d'un profil vertical de traceur OK
743c
744c tr_seri = profil initial
745c zprof   = profil de rappel
746c           (identiques a l'initialisation)
747c
748c 1 -> CO ; 2 -> OCS
749c def des profils en log(a) = a * log(P) + b par morceaux, cf. pollack et al
750c tr_seri en ppm
751
752c Constantes de rappel:
753
754       do k=1,klev
755         tau(k,1)=1.e6
756         tau(k,2)=1.e7
757         tau(k,3)=1.e8
758         tau(k,4)=1.e6
759         tau(k,5)=1.e7
760         tau(k,6)=1.e8
761       enddo
762
763c CO
764
765      do it=1,3
766       do k=1,klev
767         tr_seri(:,k,it)=0.
768c pour l'instant, tau fixe, mais possibilite de le faire varier avec z
769        if (pplay(klon/2,k) >= 1.9e6) then
770           tr_seri(:,k,it)=20.
771        endif
772        if ((pplay(klon/2,k)<=1.9e6).and.(pplay(klon/2,k)>=2.73e5)) then
773           alpha=(log(pplay(klon/2,k))-log(2.73e5))/
774     .     (log(1.9e6)-log(2.73e5))
775           tr_seri(:,k,it)=30.*(20./30.)**alpha
776        endif
777        if ((pplay(klon/2,k)<=2.73e5).and.(pplay(klon/2,k)>=1.1e4)) then
778           alpha=(log(pplay(klon/2,k))-log(1.1e4))/
779     .     (log(2.73e5)-log(1.1e4))
780           tr_seri(:,k,it)=50.*(30./50.)**alpha
781        endif
782        if ((pplay(klon/2,k)<=1.1e4).and.(pplay(klon/2,k)>=1.3e3)) then
783           alpha=(log(pplay(klon/2,k))-log(1.3e3))/
784     .     (log(1.1e4)-log(1.3e3))
785           tr_seri(:,k,it)=2.*(50./2.)**alpha
786        endif
787        if ((pplay(klon/2,k)<=1.3e3).and.(pplay(klon/2,k)>=2.4)) then
788           alpha=(log(pplay(klon/2,k))-log(2.4))/
789     .     (log(1.3e3)-log(2.4))
790           tr_seri(:,k,it)=1000.*(2./1000.)**alpha
791        endif
792        if (pplay(klon/2,k) <= 2.4) then
793           tr_seri(:,k,it)=1000.
794        endif
795       enddo
796 
797c OCS
798       do k=1,klev
799         tr_seri(:,k,it+3)=0.
800         if (pplay(klon/2,k) >= 9.45e5) then
801           tr_seri(:,k,it+3)=16.
802         endif
803         if ((pplay(klon/2,k)<=9.45e5).and.(pplay(klon/2,k)>=4.724e5))
804     *   then
805           alpha=(log(pplay(klon/2,k))-log(4.724e5))/
806     *     (log(9.45e5)-log(4.724e5))
807           tr_seri(:,k,it+3)=0.5*(16/0.5)**alpha
808         endif
809         if ((pplay(klon/2,k)<=4.724e5).and.(pplay(klon/2,k)>=1.1e4))
810     *   then
811           alpha=(log(pplay(klon/2,k))-log(1.1e4))/
812     *     (log(4.724e5)-log(1.1e4))
813           tr_seri(:,k,it+3)=0.005*(0.5/0.005)**alpha
814         endif
815         if (pplay(klon/2,k)<=1.1e4) then
816           tr_seri(:,k,it+3)=0.
817         endif
818       end do
819      enddo
820
821       do it=1,nqmax
822         do k=1,klev
823           do i=1,klon
824             zprof(i,k,it) = tr_seri(i,k,it)
825           enddo
826         enddo
827       enddo
828
829c-------------
830c fin debutphy
831c-------------
832      ENDIF  ! fin debutphy
833
834c======================================================================
835c Rappel vers un profil
836c======================================================================
837         do it=1,nqmax
838           do k=1,klev
839             do i=1,klon
840c VERIF
841           if (tr_seri(i,k,it).lt.0) then
842             print*,"Traceur negatif AVANT rappel:",i,k,it
843             stop
844           endif
845c FIN VERIF
846
847           deltatr(i,k,it) = (-tr_seri(i,k,it)+zprof(i,k,it))/tau(k,it)
848           tr_seri(i,k,it) =  tr_seri(i,k,it) + deltatr(i,k,it)*pdtphys
849
850c VERIF
851           if (tr_seri(i,k,it).lt.0) then
852             print*,"APRES rappel:",i,k,it,
853     .  deltatr(i,k,it),zprof(i,k,it),tr_seri(i,k,it),pdtphys/tau(k,it)
854             stop
855           endif
856c FIN VERIF
857             enddo
858           enddo
859         enddo
860
861c======================================================================
862c   Calcul de l'effet de la couche limite remis directement dans physiq
863c======================================================================
864
865
866      RETURN
867      END
Note: See TracBrowser for help on using the repository browser.