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

Last change on this file since 103 was 102, checked in by slebonnois, 15 years ago

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