source: trunk/LMDZ.TITAN/libf/phytitan/phytrac.F @ 134

Last change on this file since 134 was 104, checked in by slebonnois, 14 years ago

SLebonnois: modification de makelmdz et create_make_gcm pour pouvoir
compiler la chimie titan. Pas de raison que ca gene les autres.
Dans cette version, les compilations de Venus et Titan fonctionnent.

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