source: trunk/WRF.COMMON/WRFV2/mars_lmd/libf/phymars/lwu.F

Last change on this file was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 8.8 KB
Line 
1      subroutine lwu (kdlon,kflev
2     &                ,dp,plev,tlay,aerosol
3     &                ,aer_t,co2_u,co2_up
4     &                ,tautotal,omegtotal,gtotal)
5
6c----------------------------------------------------------------------
7c     LWU   computes  - co2: longwave effective absorber amounts including
8c                      pressure and temperature effects
9c                     - aerosols: amounts for every band
10c                                transmission for bandes 1 and 2 of co2
11c----------------------------------------------------------------------
12
13c-----------------------------------------------------------------------
14c ATTENTION AUX UNITES:
15c le facteur 10*g fait passer des kg m-2 aux g cm-2
16c-----------------------------------------------------------------------
17c! modif diffusion
18c! on ne change rien a la bande CO2 : les quantites d'absorbant CO2
19c! sont multipliees par 1.66
20c! pview= 1/cos(teta0)=1.66
21c-----------------------------------------------------------------------
22 
23      implicit none
24
25#include "dimensions.h"
26#include "dimphys.h"
27#include "dimradmars.h"
28#include "comcstfi.h"
29
30#include "yomaer.h"
31#include "yomlw.h"
32
33#include "fisice.h"
34#include "callkeys.h"
35
36#include "aerice.h"
37 
38c----------------------------------------------------------------------
39c         0.1   arguments
40c               ---------                                   
41c                                                            inputs:
42c                                                            -------
43      integer kdlon        ! part of ngrid
44      integer kflev        ! part of nalyer
45
46      real dp (ndlo2,kflev)         ! layer pressure thickness (Pa)
47      real plev (ndlo2,kflev+1)     ! level pressure (Pa)
48      real tlay (ndlo2,kflev)       ! layer temperature (K)
49      real aerosol (ndlo2,kflev,naerkind) ! aerosol extinction optical depth
50c                         at reference wavelength "longrefvis" set
51c                         in dimradmars.h , in each layer, for one of
52c                         the "naerkind" kind of aerosol optical properties.
53
54
55c                                                            outputs:
56c                                                            --------
57      real aer_t (ndlo2,nuco2,kflev+1)   ! transmission (aer)
58      real co2_u (ndlo2,nuco2,kflev+1)   ! absorber amounts (co2)
59      real co2_up (ndlo2,nuco2,kflev+1)  ! idem scaled by the pressure (co2)
60
61      real tautotal(ndlo2,kflev,nir)  ! \   Total single scattering
62      real omegtotal(ndlo2,kflev,nir) !  >  properties (Addition of the
63      real gtotal(ndlo2,kflev,nir)    ! /   NAERKIND aerosols properties)
64
65c----------------------------------------------------------------------
66c         0.2   local arrays
67c               ------------
68
69      integer jl,jk,jkl,ja,n
70 
71      real aer_u (ndlon,nir,nflev+1)     ! absorber amounts (aer) extinction
72      real co2c           ! co2 concentration (pa/pa)
73      real pview          ! cosecant of viewing angle
74      real pref           ! reference pressure (1013 mb = 101325 Pa)
75      real tx,tx2
76      real phi (ndlon,nuco2)
77      real psi (ndlon,nuco2)
78      real plev2 (ndlon,nflev+1)
79      real zzz
80
81      real ray,coefsize,coefsizew,coefsizeg
82
83c************************************************************************
84c----------------------------------------------------------------------
85c         0.3  Initialisation 
86c               -------------
87
88      pview = 1.66
89      co2c = 0.95           
90      pref = 101325.
91
92      do jk=1,nlaylte+1
93        do jl=1,kdlon
94          plev2(jl,jk)=plev(jl,jk)*plev(jl,jk)
95        enddo
96      enddo
97
98c----------------------------------------------------------------------
99c  Computing TOTAL single scattering parameters by adding properties of
100c  all the NAERKIND kind of aerosols in each IR band
101
102      call zerophys(ndlon*kflev*nir,tautotal)
103      call zerophys(ndlon*kflev*nir,omegtotal)
104      call zerophys(ndlon*kflev*nir,gtotal)
105
106      do n=1,naerkind
107        do ja=1,nir
108          do jk=1,nlaylte
109            do jl = 1,kdlon
110c TEST : to account for the varying sol/ir optical depth of ice with varying crystal size     
111c      : and for varying w,g with varying crystal size     
112              if (activice.and.n.eq.naerkind) then
113        ray=min( max(rice(jl,jk)*1.e+6, 1.),10.)
114            if (ja.eq.1.or.ja.eq.2) then
115          coefsize=(-0.00382*ray**3.+0.0503*ray**2.+0.03531)
116     &             /QIRsQREF(ja,n)
117          coefsizew=(-0.011*ray**2.+0.1824*ray-0.1283)
118     &           /omegaIR(ja,n)
119          coefsizeg=(-0.00122*ray**3.+0.009161*ray**2.+ 
120     &          0.1182*ray-0.096877)/gIR(ja,n)
121                elseif (ja.eq.3) then
122          coefsize=(-0.00324*ray**3.+0.0419*ray**2.+0.0459)
123     &             /QIRsQREF(ja,n)
124          coefsizew=(-0.01292*ray**2.+0.1963*ray-0.06566)
125     &           /omegaIR(ja,n)
126          coefsizeg=(0.00271*ray**3.-0.05959*ray**2.+
127     &          0.4411*ray-0.2724)/gIR(ja,n)
128                elseif (ja.eq.4) then
129          coefsize=(-0.0003823*ray**3.+0.0104*ray**2.+0.005966)
130     &             /QIRsQREF(ja,n)
131              coefsizew=(-0.002*ray**3.+0.02623*ray**2.-0.014465)
132     &           /omegaIR(ja,n)
133          coefsizeg=(-0.0017192*ray**3.+0.0259*ray**2.-
134     &          0.027692*ray+0.016099)/gIR(ja,n)
135            endif
136        if (coefsize.le.0.or.coefsizew.le.0.or.coefsizeg.le.0)
137     &        stop'pb dans lwu avec prp opt glace'
138
139            tautotal(jl,jk,ja)=tautotal(jl,jk,ja) +
140     &          QIRsQREF(ja,n)*aerosol(jl,jk,n)
141     &          *coefsize
142                omegtotal(jl,jk,ja) =  omegtotal(jl,jk,ja) +
143     &           QIRsQREF(ja,n)*aerosol(jl,jk,n)*omegaIR(ja,n)
144     &          *coefsize*coefsizew
145                gtotal(jl,jk,ja) =  gtotal(jl,jk,ja) +
146     &           QIRsQREF(ja,n)*aerosol(jl,jk,n)*omegaIR(ja,n)*gIR(ja,n)
147     &          *coefsize*coefsizew*coefsizeg     
148              else
149              tautotal(jl,jk,ja)=tautotal(jl,jk,ja) +
150     &           QIRsQREF(ja,n)*aerosol(jl,jk,n)
151              omegtotal(jl,jk,ja) =  omegtotal(jl,jk,ja) +
152     &           QIRsQREF(ja,n)*aerosol(jl,jk,n)*omegaIR(ja,n)
153              gtotal(jl,jk,ja) =  gtotal(jl,jk,ja) +
154     &           QIRsQREF(ja,n)*aerosol(jl,jk,n)*omegaIR(ja,n)*gIR(ja,n)
155              endif
156            enddo
157          enddo
158        enddo
159      enddo
160      do ja=1,nir
161        do jk=1,nlaylte
162          do jl = 1,kdlon
163             gtotal(jl,jk,ja)=gtotal(jl,jk,ja)/omegtotal(jl,jk,ja)
164             omegtotal(jl,jk,ja)=omegtotal(jl,jk,ja)/tautotal(jl,jk,ja)
165          enddo
166        enddo
167      enddo
168
169c----------------------------------------------------------------------
170c         1.0   cumulative (aerosol) amounts (for every band)
171c               ----------------------------
172
173      jk=nlaylte+1
174      do ja=1,nir
175        do jl = 1 , kdlon
176          aer_u(jl,ja,jk)=0.
177        enddo
178      enddo
179
180      do jk=1,nlaylte
181        jkl=nlaylte+1-jk
182        do ja=1,nir
183          do jl=1,kdlon
184             aer_u(jl,ja,jkl)=aer_u(jl,ja,jkl+1)+ tautotal(jl,jkl,ja)
185          enddo
186        enddo
187      enddo                 
188
189c----------------------------------------------------------------------
190c         1.0   bands 1 and 2 of co2
191c               --------------------
192
193      jk=nlaylte+1
194      do ja=1,nuco2
195        do jl = 1 , kdlon
196          co2_u(jl,ja,jk)=0.
197          co2_up(jl,ja,jk)=0.
198          aer_t(jl,ja,jk)=1.
199        enddo
200      enddo
201
202      do jk=1,nlaylte
203                       jkl=nlaylte+1-jk
204        do ja=1,nuco2
205          do jl=1,kdlon
206
207c                 introduces temperature effects on absorber(co2) amounts
208c                 -------------------------------------------------------
209            tx = sign(min(abs(tlay(jl,jkl)-tref),70.)
210     .         ,tlay(jl,jkl)-tref)
211            tx2=tx*tx
212            phi(jl,ja)=at(1,ja)*tx+bt(1,ja)*tx2
213            psi(jl,ja)=at(2,ja)*tx+bt(2,ja)*tx2
214            phi(jl,ja)=exp(phi(jl,ja)/cst_voigt(2,ja))
215            psi(jl,ja)=exp(2.*psi(jl,ja))
216
217c                                        cumulative absorber(co2) amounts
218c                                        --------------------------------
219            co2_u(jl,ja,jkl)=co2_u(jl,ja,jkl+1)
220     .     +         pview/(10*g)*phi(jl,ja)*dp(jl,jkl)*co2c
221
222            co2_up(jl,ja,jkl)=co2_up(jl,ja,jkl+1)
223     .     +         pview/(10*g*2*pref)*psi(jl,ja)
224     .     *        (plev2(jl,jkl)-plev2(jl,jkl+1))*co2c
225
226c                                                  (aerosol) transmission
227c                                                  ----------------------
228c   on calcule directement les transmissions pour les aerosols.
229c   on multiplie le Qext  par 1-omega dans la bande du CO2.
230c   et pourquoi pas d'abord?  hourdin@lmd.ens.fr
231
232           zzz=pview*(1.-omegtotal(jl,jkl,ja))
233           aer_t(jl,ja,jkl)=exp(-zzz*aer_u(jl,ja,jkl))
234         
235          enddo
236        enddo
237      enddo                 
238
239c----------------------------------------------------------------------
240      return
241      end
Note: See TracBrowser for help on using the repository browser.