source: trunk/libf/phytitan/phytrac.F @ 21

Last change on this file since 21 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: 11.8 KB
RevLine 
[3]1      SUBROUTINE phytrac (firstcall,lastcall,
2     .                   nqmax,nmicro,ptimestep,appkim,dtkim,
3     .                   pplev,pplay,delp,ptemp,pmu0,pfract,pdecli,
4     .                   lonsol,tr_seri,d_tr_mph,d_tr_kim)
5
6      IMPLICIT none
7
8c======================================================================
9c    S. Lebonnois, mai 2008
10c
11c  Arguments:
12c
13c firstcall----input-L-variable logique indiquant le premier passage
14c lastcall-----input-L-variable logique indiquant le dernier passage
15c nqmax--------input-I-nombre de traceurs (total)
16c nmicro-------input-I-nombre de traceurs microphysiques !! doivent etre toujours en premiers!!
17c ptimestep----input-R-pas d'integration pour la physique (seconde)
18c appkim-------input-I-appel a la chimie
19c dtkim--------input-R-pas de temps chimique (seconde)
20c pplev--------input-R-pression pour chaque inter-couche (en Pa)
21c pplay--------input-R-pression pour chaque couche (en Pa)
22c delp---------input-R-epaisseur d'une couche (en Pa)
23c ptemp--------input-R-temperature (K)
24c pmu0---------input-R-cos angle zenithal
25c pfract-------input-R-fractional day
26c pdecli-------input-R-declinaison en radian
27c lonsol-------input-R-longitude solaire en radian
28c tr_seri------input-R-mass mixing ratio traceurs (kg/kg)
29c d_tr_mph----output-R-tendance microphysique de "qx" (kg/kg/s)
30c d_tr_kim----output-R-tendance chimique de "qx" (kg/kg/s)
31c======================================================================
32#include "dimensions.h"
33#include "dimphy.h"
34#include "clesphys.h"
35#include "YOMCST.h"
36#include "advtrac.h"
37
38c======================================================================
39c Variables argument:
40c
41      LOGICAL firstcall,lastcall
42      INTEGER nqmax,nmicro,nlat,appkim
43      REAL ptimestep,dtkim
44      REAL pplev(klon,klev+1),pplay(klon,klev+1),delp(klon,klev)
45      REAL ptemp(klon,klev)
46      REAL pmu0(klon), pfract(klon), pdecli, lonsol
47      REAL tr_seri(klon,klev,nqmax)
48      REAL d_tr_mph(klon,klev,nqmax),d_tr_kim(klon,klev,nqmax)
49
50c======================================================================
51c Local variables
52
53* common relatifs aux aerosols
54      REAL qaer(klon,klev,nqmx)
55      common/traceurs/qaer
56
57c grandeurs en moyennes zonales
58      REAL zplev(jjm+1,klev+1),zplay(jjm+1,klev)
59      REAL ztemp(jjm+1,klev),zmu0(jjm+1),zfract(jjm+1)
60      real temp_eq(klev),press_eq(klev)
61      REAL zqaer(jjm+1,klev,nqmax)    ! et non nmicro... Permet nmicro=0.
62      REAL zqaer0(jjm+1,klev,nqmax)
63      REAL pdqmfi(jjm+1,klev,nqmax)
64      REAL ychim(jjm+1,klev,nqmax-nmicro)
65      REAL qysat(klev,nqmx)        ! dim nqmx, mais en fait nqmax-nmicro (save...)
66      REAL pdyfi(jjm+1,klev,nqmx)  ! dim nqmx, mais en fait nqmax-nmicro (save...)
67      character*10 nomqy(nqmax-nmicro+1)
68      integer      i,j,l,iq,ig0
69
70c La saturation n est calculee qu une seule fois: sauvegarde qysat
71c La chimie n est pas calculee tous les pas, il faut donc
72c                      sauvegarder les sorties de la chimie
73
74      SAVE pdyfi,qysat 
75     
76c======================================================================
77c======================================================================
78
79c-----------------------------------------------------------------------
80c   convertion moyennes zonales et changement d unites pour microphy
81c   ---------------------------------
82
83c     print*,'CONVERSION 2D ET CHANGEMENT UNITES (PHYTRAC)'
84
85      zplev = 0.0
86      zplay = 0.0
87      ztemp = 0.0
88      zqaer = 0.0
89      ychim = 0.0
90      zmu0  = 0.0
91      zfract= 0.0
92     
93       do l=1,llm+1
94         zplev(1,l) = pplev(1,l)
95         do j=2,jjm
96            ig0=1+(j-2)*iim
97            do i=1,iim
98               zplev(j,l) = zplev(j,l) + pplev(ig0+i,l)/iim
99            enddo
100         enddo
101         zplev(jjm+1,l) = pplev(klon,l)
102       enddo
103
104       do l=1,llm
105         ztemp(1,l) = ptemp(1,l)
106         zplay(1,l) = pplay(1,l)
107         do j=2,jjm
108            ig0=1+(j-2)*iim
109            do i=1,iim
110               ztemp(j,l) = ztemp(j,l) + ptemp(ig0+i,l)/iim
111               zplay(j,l) = zplay(j,l) + pplay(ig0+i,l)/iim
112            enddo
113         enddo
114         ztemp(jjm+1,l) = ptemp(klon,l)
115         zplay(jjm+1,l) = pplay(klon,l)
116         temp_eq  = ztemp((jjm+1)/2,:)
117         press_eq = zplay((jjm+1)/2,:)/100. ! en mbar
118       enddo
119
120         zmu0(1)   = pmu0(1)
121         zfract(1) = pfract(1)
122         do j=2,jjm
123            ig0=1+(j-2)*iim
124            do i=1,iim
125               zmu0(j)   = zmu0(j)   + pmu0(ig0+i)/iim
126               zfract(j) = zfract(j) + pfract(ig0+i)/iim
127            enddo
128         enddo
129         zmu0(jjm+1)   = pmu0(klon)
130         zfract(jjm+1) = pfract(klon)
131
132c TRACEURS MICROPHYSIQUES
133
134      if (microfi.eq.1) then
135
136      do iq=1,nmicro
137c        print*,tname(iq)
138
139c       Traceurs microphysiques: passage en extensif: n/kg --> n/m^2
140         DO l=1,llm
141          DO i = 1, klon
142            qaer(i,l,iq) = tr_seri(i,l,iq)*delp(i,l)/RG
143          ENDDO
144         ENDDO
145
146       do l=1,llm
147         zqaer(1,l,iq) = qaer(1,l,iq)
148         do j=2,jjm
149            ig0=1+(j-2)*iim
150            do i=1,iim
151               zqaer(j,l,iq) = zqaer(j,l,iq) + qaer(ig0+i,l,iq)/iim
152            enddo
153         enddo
154         zqaer(jjm+1,l,iq) = qaer(klon,l,iq)
155       enddo
156      enddo
157
158      endif   ! microfi
159
160c AUTRES TRACEURS
161     
162      if (nqmax.gt.nmicro) then
163      do iq=nmicro+1,nqmax
164       do l=1,llm
165         ychim(1,l,iq-nmicro) = tr_seri(1,l,iq)
166         do j=2,jjm
167            ig0=1+(j-2)*iim
168            do i=1,iim
169               ychim(j,l,iq-nmicro) = ychim(j,l,iq-nmicro)
170     .                              + tr_seri(ig0+i,l,iq)/iim
171            enddo
172         enddo
173         ychim(jjm+1,l,iq-nmicro) = tr_seri(klon,l,iq)
174       enddo
175         nomqy(iq-nmicro) = tname(iq)
176
177c       print*,iq-nmicro,nomqy(iq-nmicro)
178
179      enddo
180      nomqy(nqmax-nmicro+1) = "HV"
181      endif
182
183c-----------------------------------------------------------------------
184c   initialisation des qysat au premier appel:
185c   ---------------------------------
186
187c!! ATTENTION, qysat pris uniquement a l'equateur
188c!!  justifie puisque dans cette region, les var de t et p sont faibles...
189
190      if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
191           call inicondens(nqmax-nmicro,press_eq,temp_eq,nomqy,qysat)
192      endif
193
194c-----------------------------------------------------------------------
195c     Appel de la microphysique
196c    --------------------------
197
198       pdqmfi = 0.
199       zqaer0 = zqaer
200       
201       IF(firstcall) THEN
202        print*,'MICROPHYSIQUE ',MICROFI
203       ENDIF
204
205       IF(MICROFI.eq.1) THEN
206
207         IF(firstcall) THEN
208          print*,'OH! On passe dans la microphysique'
209         ENDIF
210
211         CALL pg2(zplev,ztemp,zqaer,pdqmfi    ! tendances 2D, /s
212     &   ,nmicro,ptimestep,zmu0,zfract,lastcall)
213
214       ELSE
215
216        IF(firstcall) THEN
217          print*,'MICROPHYSIQUE OFF-LINE',MICROFI
218          if (nmicro.gt.0) then
219            CALL pg2(zplev,ztemp,zqaer,pdqmfi
220     &    ,nmicro,ptimestep,zmu0,zfract,lastcall)
221          endif
222        ENDIF
223
224       ENDIF
225
226       zqaer  = zqaer+pdqmfi*ptimestep
227       pdqmfi = (zqaer-zqaer0)/ptimestep
228       
229c-----------------------------------------------------------------------
230c     Condensation
231c    -------------
232
233      if ((chimi).and.(nqmax.gt.nmicro)) then
234
235c   tendance (en /s) passee sur pdqmfi(nmicro+1 a nqmax)
236c        print*,'Condensation'
237
238        do iq=1,nqmax-nmicro
239           do l=1,llm
240              do j=1,jjm+1
241                 if (ychim(j,l,iq).gt.qysat(l,iq)) then
242           pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+qysat(l,iq)) !delta y
243     .                            / ptimestep                  ! / dt
244                 endif
245              enddo
246           enddo
247        enddo
248
249      ENDIF
250
251c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
252c eventuellement, modif initiale de la compo
253c
254c   tendance (en /s) passee sur pdqmfi(nmicro+1 a nqmax)
255c
256c     if (firstcall .and. chimi .and.(nqmax.gt.nmicro)) then
257c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
258c!!!remise de CH4 a 1.5%!!!!!!!!!!!!!!!!!!!!!!
259c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
260c       do iq=1,nqmax-nmicro
261c         if (nomqy(iq).eq."CH4") then
262c          do l=1,llm
263c             do j=1,jjm+1
264c                if (ychim(j,l,iq).le.0.015) then
265c          pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+0.015) !delta y
266c    .                            / ptimestep                  ! / dt
267c                endif
268c             enddo
269c          enddo
270c         endif
271c       enddo
272c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
273c         
274c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
275c!!!remise de C2H2 a 1.e-5 max !!!!!!!!!!!!!
276c!!!remise de C2H6 a 3.e-5 max !!!!!!!!!!!!!
277c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
278c       do iq=1,nqmax-nmicro
279c         if (nomqy(iq).eq."C2H2") then
280c          do l=1,llm
281c             do j=1,jjm+1
282c                if (ychim(j,l,iq).gt.1.e-5) then
283c          pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+1.e-5) !delta y
284c    .                            / ptimestep                  ! / dt
285c                endif
286c             enddo
287c          enddo
288c         endif
289c         if (nomqy(iq).eq."C2H6") then
290c          do l=1,llm
291c             do j=1,jjm+1
292c                if (ychim(j,l,iq).gt.3.e-5) then
293c          pdqmfi(j,l,nmicro+iq)= (-ychim(j,l,iq)+3.e-5) !delta y
294c    .                            / ptimestep                  ! / dt
295c                endif
296c             enddo
297c          enddo
298c         endif
299c       enddo
300c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
301c     endif
302c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
303         
304c-----------------------------------------------------------------------
305c     Appel de la chimie
306c    --------------------------
307
308      if((appkim.eq.1).and.(chimi)) then
309        print*,'On passe dans la CHIMIE'
310
311c       do iq=1,nqmax-nmicro
312c         if (nomqy(iq).eq."C2H2") then
313c           print*,"C2H2top=",ychim(:,klev,iq)
314c         endif
315c       enddo
316
317c Appel Chimie
318c ------------
319       CALL calchim(nqmax-nmicro,ychim,nomqy,pdecli,lonsol,dtkim,
320     .              ztemp,zplay,zplev,
321     .              pdyfi)   
322c ychim ne doit pas etre modifie, pdyfi en /s
323         
324      endif
325     
326c-----------------------------------------------------------------------
327c   retour des tendances vers 3D
328c   ---------------------------------
329
330c TRACEURS MICROPHYSIQUES
331
332      if (microfi.eq.1) then
333
334      DO iq=1,nmicro
335         DO l=1,llm
336            d_tr_mph(1,l,iq) = pdqmfi(1,l,iq)
337            ig0          = 2
338            DO j=2,jjm
339               DO i = 1, iim
340                  d_tr_mph(ig0,l,iq)  = pdqmfi(j,l,iq)
341     &                 *qaer(ig0,l,iq)/zqaer0(j,l,iq)
342                  ig0             = ig0 + 1
343               ENDDO
344            ENDDO
345            d_tr_mph(ig0,l,iq) = pdqmfi(jjm+1,l,iq)
346         ENDDO
347      ENDDO
348
349      do iq=1,nmicro
350         DO l=1,llm
351          DO i = 1, klon
352c incrementation de la tendance sur qaer (pour sorties dans physiq.F)
353            qaer(i,l,iq) = qaer(i,l,iq) + d_tr_mph(i,l,iq)*ptimestep
354c  Traceurs microphysiques: passage en intensif: n/m^2 --> n/kg
355            d_tr_mph(i,l,iq) = d_tr_mph(i,l,iq)*RG/delp(i,l)
356          ENDDO
357         ENDDO
358      enddo
359
360      endif   ! microfi
361
362c AUTRES TRACEURS
363
364      if ((chimi).and.(nqmax.gt.nmicro)) then
365c on passe de pdyfi (tendance chimique en /s calculee quand chimie appelee)
366c          a  d_tr_kim (tendance chimique 3D en /s, passee a physiq)
367c  et de pdqmfi a d_tr_mph (tendance condensation 3D en /s passee a physiq)
368
369      DO iq=nmicro+1,nqmax
370         DO l=1,llm
371            d_tr_kim(1,l,iq) = pdyfi(1,l,iq-nmicro)
372            d_tr_mph(1,l,iq) = pdqmfi(1,l,iq)
373            ig0          = 2
374            DO j=2,jjm
375               DO i = 1, iim
376                  d_tr_kim(ig0,l,iq)  = pdyfi(j,l,iq-nmicro)
377     &                 *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
378                  d_tr_mph(ig0,l,iq)  = pdqmfi(j,l,iq)
379     &                 *tr_seri(ig0,l,iq)/ychim(j,l,iq-nmicro)
380                  ig0             = ig0 + 1
381               ENDDO
382            ENDDO
383            d_tr_kim(ig0,l,iq) = pdyfi(jjm+1,l,iq-nmicro)
384            d_tr_mph(ig0,l,iq) = pdqmfi(jjm+1,l,iq)
385         ENDDO
386      ENDDO
387
388      endif   ! chimi
389
390      RETURN
391      END
Note: See TracBrowser for help on using the repository browser.