source: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_relax.F @ 1443

Last change on this file since 1443 was 1443, checked in by emillour, 9 years ago

Titan and Venus GCMs:
Follow-up to the changes in dynamics/physics interface: ener.h, logic.h, serre.h and temps.h are now modules.
EM
File size: 8.2 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_relax (debutphy,lafin,
7     I                    nqmax,
8     I                    nlon,
9     I                    nlev,
10     I                    pdtphys,
11     I                    pplay,
12     O                    tr_seri)
13
14c======================================================================
15c Auteur(s) FH
16c Objet: Moniteur general des tendances traceurs
17c
18cAA Remarques en vrac:
19cAA--------------------
20cAA 1/ le call phytrac se fait avec nqmax
21c
22c SL: Janvier 2014
23c Version developed by E. Marcq for pseudo-chemistry relaxation
24c See Marcq&Lebonnois 2013.
25c
26c======================================================================
27      USE ioipsl
28      USE infotrac
29      USE control_mod
30      use dimphy
31      USE comgeomphy
32      IMPLICIT none
33#include "YOMCST.h"
34#include "dimensions.h"
35#include "clesphys.h"
36#include "paramet.h"
37c======================================================================
38
39c Arguments:
40
41c   EN ENTREE:
42c   ==========
43
44      logical debutphy       ! le flag de l'initialisation de la physique
45      logical lafin          ! le flag de la fin de la physique
46      integer nqmax ! nombre de traceurs auxquels on applique la physique
47      integer nlon  ! nombre de points horizontaux
48      integer nlev  ! nombre de couches verticales
49      real pdtphys  ! pas d'integration pour la physique (seconde)
50      real pplay(nlon,nlev)  ! pression pour le mileu de chaque couche (en Pa)
51
52c   EN ENTREE/SORTIE:
53c   =================
54
55      real tr_seri(nlon,nlev,nqmax) ! traceur 
56
57cAA ----------------------------
58cAA  VARIABLES LOCALES TRACEURS
59cAA ----------------------------
60
61C les traceurs
62
63c===================
64c it--------indice de traceur
65c k,i---------indices long, vert
66c===================
67c Variables deja declarees dont on a besoin pour traceurs   
68c   k,i,it,tr_seri(nlon,nlev,nqmax),pplay(nlon,nlev),
69      integer nqCO_OCS
70c      real pzero,gamma
71c      parameter (pzero=85000.)
72c      parameter (gamma=5000.)
73      REAL alpha
74      real deltatr(nlon,nlev,nqtot) ! ecart au profil de ref zprof
75      real,save,allocatable :: zprof(:,:)
76      real,save,allocatable :: tau(:,:) ! temps de relaxation vers le profil (s)
77c======================================================================
78
79      INTEGER i, k, it
80
81c Variables liees a l'ecriture de la bande histoire physique
82
83c Variables locales pour effectuer les appels en serie
84c----------------------------------------------------
85
86      REAL d_tr(nlon,nlev) ! tendances de traceurs
87
88      character*20 modname
89      character*80 abort_message
90
91c======================================================================
92
93      modname = 'phytrac_relax'
94c TRACEURS TYPE CO ET OCS
95      nqCO_OCS   = 6
96
97c---------
98c debutphy
99c---------
100      if (debutphy) then
101         print*,"DEBUT PHYTRAC"
102         print*,"PHYTRAC: RELAXATION"
103         allocate(zprof(nlev,nqtot),tau(nlev,nqtot))
104
105c=============================================================
106c=============================================================
107c=============================================================
108c   Initialisation des traceurs
109c=============================================================
110c=============================================================
111c=============================================================
112
113C=========================================================================
114C=========================================================================
115
116c II) Declaration d'un profil vertical de traceur OK
117c
118c zprof   = profil de rappel
119c
120c 1 -> CO ; 2 -> OCS
121c def des profils en log(a) = a * log(P) + b par morceaux, cf. pollack et al
122c tr_seri en ppm
123c (initialisation seulement si ceux-ci sont nuls)
124
125c ICI, ON UTILISE 3 CONSTANTES DE TEMPS DIFFERENTES POUR CHAQUE,
126c DONC TRACEURS 1 A 3 POUR CO ET 4 A 6 POUR OCS
127C=========================================================================
128
129
130c Constantes de rappel:
131
132       print*,"INIT TAU"
133       do k=1,nlev
134         tau(k,1)=1.e6
135         tau(k,2)=1.e7
136         tau(k,3)=1.e8
137         tau(k,4)=1.e6
138         tau(k,5)=1.e7
139         tau(k,6)=1.e8
140       enddo
141
142c CO
143
144      do it=1,3
145       print*,"INIT ZPROF ",tname(it)
146       do k=1,nlev
147         zprof(k,it)=0.
148c pour l'instant, tau fixe, mais possibilite de le faire varier avec z
149        if (pplay(nlon/2,k) >= 4.8e6) then
150           zprof(k,it)=14.
151        endif
152        if ((pplay(nlon/2,k)<=4.8e6).and.(pplay(nlon/2,k)>=1.9e6)) then
153           alpha=(log(pplay(nlon/2,k))-log(1.9e6))/
154     .     (log(4.8e6)-log(1.9e6))
155           zprof(k,it)=20.*(14./20.)**alpha
156        endif
157        if ((pplay(nlon/2,k)<=1.9e6).and.(pplay(nlon/2,k)>=1.5e5)) then
158           alpha=(log(pplay(nlon/2,k))-log(1.5e5))/
159     .     (log(1.9e6)-log(1.5e5))
160           zprof(k,it)=39.*(20./39.)**alpha
161        endif
162        if ((pplay(nlon/2,k)<=1.5e5).and.(pplay(nlon/2,k)>=1.1e4)) then
163           alpha=(log(pplay(nlon/2,k))-log(1.1e4))/
164     .     (log(2.73e5)-log(1.1e4))
165           zprof(k,it)=50.*(39./50.)**alpha
166        endif
167        if ((pplay(nlon/2,k)<=1.1e4).and.(pplay(nlon/2,k)>=1.3e3)) then
168           alpha=(log(pplay(nlon/2,k))-log(1.3e3))/
169     .     (log(1.1e4)-log(1.3e3))
170           zprof(k,it)=2.*(50./2.)**alpha
171        endif
172        if ((pplay(nlon/2,k)<=1.3e3).and.(pplay(nlon/2,k)>=2.4)) then
173           alpha=(log(pplay(nlon/2,k))-log(2.4))/
174     .     (log(1.3e3)-log(2.4))
175           zprof(k,it)=1000.*(2./1000.)**alpha
176        endif
177        if (pplay(nlon/2,k) <= 2.4) then
178           zprof(k,it)=1000.
179        endif
180       enddo
181       print*,zprof(:,it)
182 
183c OCS
184       print*,"INIT ZPROF ",tname(it+3)
185       do k=1,nlev
186         zprof(k,it+3)=0.
187         if (pplay(nlon/2,k) >= 4.8e6) then
188           zprof(k,it+3)=30.
189         endif
190         if ((pplay(nlon/2,k)<=4.8e6).and.(pplay(nlon/2,k)>=9.4e5))
191     *   then
192           alpha=(log(pplay(nlon/2,k))-log(9.4e5))/
193     *     (log(4.8e6)-log(9.4e5))
194           zprof(k,it+3)=20.*(30/20.)**alpha
195         endif
196         if ((pplay(nlon/2,k)<=9.4e5).and.(pplay(nlon/2,k)>=4.724e5))
197     *   then
198           alpha=(log(pplay(nlon/2,k))-log(4.724e5))/
199     *     (log(9.4e5)-log(4.724e5))
200           zprof(k,it+3)=0.5*(20/0.5)**alpha
201         endif
202         if ((pplay(nlon/2,k)<=4.724e5).and.(pplay(nlon/2,k)>=1.1e4))
203     *   then
204           alpha=(log(pplay(nlon/2,k))-log(1.1e4))/
205     *     (log(4.724e5)-log(1.1e4))
206           zprof(k,it+3)=0.005*(0.5/0.005)**alpha
207         endif
208         if (pplay(nlon/2,k)<=1.1e4) then
209           zprof(k,it+3)=0.
210         endif
211       end do
212       print*,zprof(:,it+3)
213      enddo
214
215c Initialisation du traceur s'il est nul:
216       do it=1,nqCO_OCS
217        if ((tr_seri(nlon/2,1,it).eq.0.).and.
218     .      (tr_seri(nlon/2,nlev/2,it).eq.0.).and.
219     .      (tr_seri(nlon/2,nlev,it).eq.0.)) then
220         print*,"INITIALISATION DE ",tname(it)
221         do k=1,nlev
222           do i=1,nlon
223             tr_seri(i,k,it) = zprof(k,it)
224           enddo
225         enddo
226        endif
227       enddo
228
229C=========================================================================
230C=========================================================================
231
232c-------------
233c fin debutphy
234c-------------
235      ENDIF  ! fin debutphy
236
237c======================================================================
238c Rappel vers un profil
239c======================================================================
240         do it=1,nqCO_OCS
241           do k=1,nlev
242             do i=1,nlon
243c VERIF
244           if (tr_seri(i,k,it).lt.0) then
245             print*,"Traceur negatif AVANT rappel:",i,k,it
246             stop
247           endif
248c FIN VERIF
249
250           deltatr(i,k,it) = (-tr_seri(i,k,it)+zprof(k,it))/tau(k,it)
251           tr_seri(i,k,it) =  tr_seri(i,k,it) + deltatr(i,k,it)*pdtphys
252
253c VERIF
254           if (tr_seri(i,k,it).lt.0) then
255             print*,"APRES rappel:",i,k,it,
256     .  deltatr(i,k,it),zprof(k,it),tr_seri(i,k,it),pdtphys/tau(k,it)
257             stop
258           endif
259c FIN VERIF
260             enddo
261           enddo
262         enddo
263
264c======================================================================
265c======================================================================
266
267
268      RETURN
269      END
270     
271     
Note: See TracBrowser for help on using the repository browser.