source: trunk/LMDZ.MARS/libf/phymars/lwu.F @ 308

Last change on this file since 308 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM

File size: 7.6 KB
Line 
1      subroutine lwu (kdlon,kflev
2     &                ,dp,plev,tlay,aerosol
3     &                ,QIRsQREF3d,omegaIR3d,gIR3d
4     &                ,aer_t,co2_u,co2_up
5     &                ,tautotal,omegtotal,gtotal)
6
7c----------------------------------------------------------------------
8c     LWU   computes  - co2: longwave effective absorber amounts including
9c                      pressure and temperature effects
10c                     - aerosols: amounts for every band
11c                                transmission for bandes 1 and 2 of co2
12c----------------------------------------------------------------------
13
14c-----------------------------------------------------------------------
15c ATTENTION AUX UNITES:
16c le facteur 10*g fait passer des kg m-2 aux g cm-2
17c-----------------------------------------------------------------------
18c! modif diffusion
19c! on ne change rien a la bande CO2 : les quantites d'absorbant CO2
20c! sont multipliees par 1.66
21c! pview= 1/cos(teta0)=1.66
22c
23c Modif J.-B. Madeleine: Computing optical properties of dust and
24c   water-ice crystals in each gridbox. Optical parameters of
25c   water-ice clouds are convolved to crystal sizes predicted by
26c   the microphysical scheme.
27c
28c MODIF : FF : removing the monster bug on water ice clouds 11/2010
29c-----------------------------------------------------------------------
30 
31      implicit none
32
33#include "dimensions.h"
34#include "dimphys.h"
35#include "dimradmars.h"
36#include "comcstfi.h"
37
38#include "yomaer.h"
39#include "yomlw.h"
40
41#include "callkeys.h"
42
43c----------------------------------------------------------------------
44c         0.1   arguments
45c               ---------                                   
46c                                                            inputs:
47c                                                            -------
48      integer kdlon        ! part of ngrid
49      integer kflev        ! part of nalyer
50
51      real dp (ndlo2,kflev)         ! layer pressure thickness (Pa)
52      real plev (ndlo2,kflev+1)     ! level pressure (Pa)
53      real tlay (ndlo2,kflev)       ! layer temperature (K)
54      real aerosol (ndlo2,kflev,naerkind) ! aerosol extinction optical depth
55c                         at reference wavelength "longrefvis" set
56c                         in dimradmars.h , in each layer, for one of
57c                         the "naerkind" kind of aerosol optical properties.
58      REAL QIRsQREF3d(ndlo2,kflev,nir,naerkind)  ! 3d ext. coef.
59      REAL omegaIR3d(ndlo2,kflev,nir,naerkind)   ! 3d ssa
60      REAL gIR3d(ndlo2,kflev,nir,naerkind)       ! 3d assym. param.
61
62c                                                            outputs:
63c                                                            --------
64      real aer_t (ndlo2,nuco2,kflev+1)   ! transmission (aer)
65      real co2_u (ndlo2,nuco2,kflev+1)   ! absorber amounts (co2)
66      real co2_up (ndlo2,nuco2,kflev+1)  ! idem scaled by the pressure (co2)
67
68      real tautotal(ndlo2,kflev,nir)  ! \   Total single scattering
69      real omegtotal(ndlo2,kflev,nir) !  >  properties (Addition of the
70      real gtotal(ndlo2,kflev,nir)    ! /   NAERKIND aerosols properties)
71
72c----------------------------------------------------------------------
73c         0.2   local arrays
74c               ------------
75
76      integer jl,jk,jkl,ja,n
77 
78      real aer_a (ndlon,nir,nflev+1) ! absorber amounts (aer) ABSORPTION
79      real co2c           ! co2 concentration (pa/pa)
80      real pview          ! cosecant of viewing angle
81      real pref           ! reference pressure (1013 mb = 101325 Pa)
82      real tx,tx2
83      real phi (ndlon,nuco2)
84      real psi (ndlon,nuco2)
85      real plev2 (ndlon,nflev+1)
86      real zzz
87
88      real ray,coefsize,coefsizew,coefsizeg
89
90c************************************************************************
91c----------------------------------------------------------------------
92c         0.3  Initialisation 
93c               -------------
94
95      pview = 1.66
96      co2c = 0.95           
97      pref = 101325.
98
99      do jk=1,nlaylte+1
100        do jl=1,kdlon
101          plev2(jl,jk)=plev(jl,jk)*plev(jl,jk)
102        enddo
103      enddo
104
105c----------------------------------------------------------------------
106c  Computing TOTAL single scattering parameters by adding properties of
107c  all the NAERKIND kind of aerosols in each IR band
108
109      call zerophys(ndlon*kflev*nir,tautotal)
110      call zerophys(ndlon*kflev*nir,omegtotal)
111      call zerophys(ndlon*kflev*nir,gtotal)
112
113      do n=1,naerkind
114        do ja=1,nir
115          do jk=1,nlaylte
116            do jl = 1,kdlon
117              tautotal(jl,jk,ja)=tautotal(jl,jk,ja) +
118     &           QIRsQREF3d(jl,jk,ja,n)*aerosol(jl,jk,n)
119              omegtotal(jl,jk,ja) =  omegtotal(jl,jk,ja) +
120     &           QIRsQREF3d(jl,jk,ja,n)*aerosol(jl,jk,n)*
121     &           omegaIR3d(jl,jk,ja,n)
122              gtotal(jl,jk,ja) =  gtotal(jl,jk,ja) +
123     &           QIRsQREF3d(jl,jk,ja,n)*aerosol(jl,jk,n)*
124     &           omegaIR3d(jl,jk,ja,n)*gIR3d(jl,jk,ja,n)
125            enddo
126          enddo
127        enddo
128      enddo
129      do ja=1,nir
130        do jk=1,nlaylte
131          do jl = 1,kdlon
132             gtotal(jl,jk,ja)=gtotal(jl,jk,ja)/omegtotal(jl,jk,ja)
133             omegtotal(jl,jk,ja)=omegtotal(jl,jk,ja)/tautotal(jl,jk,ja)
134          enddo
135        enddo
136      enddo
137
138c----------------------------------------------------------------------
139c         1.0   cumulative (aerosol) amounts (for every band)
140c               ----------------------------
141
142      jk=nlaylte+1
143      do ja=1,nir
144        do jl = 1 , kdlon
145          aer_a(jl,ja,jk)=0.
146        enddo
147      enddo
148
149      do jk=1,nlaylte
150        jkl=nlaylte+1-jk
151        do ja=1,nir
152          do jl=1,kdlon
153             aer_a(jl,ja,jkl)=aer_a(jl,ja,jkl+1)+
154     &           tautotal(jl,jkl,ja)*(1.-omegtotal(jl,jkl,ja))
155          enddo
156        enddo
157      enddo                 
158
159c----------------------------------------------------------------------
160c         1.0   bands 1 and 2 of co2
161c               --------------------
162
163      jk=nlaylte+1
164      do ja=1,nuco2
165        do jl = 1 , kdlon
166          co2_u(jl,ja,jk)=0.
167          co2_up(jl,ja,jk)=0.
168          aer_t(jl,ja,jk)=1.
169        enddo
170      enddo
171
172      do jk=1,nlaylte
173                       jkl=nlaylte+1-jk
174        do ja=1,nuco2
175          do jl=1,kdlon
176
177c                 introduces temperature effects on absorber(co2) amounts
178c                 -------------------------------------------------------
179            tx = sign(min(abs(tlay(jl,jkl)-tref),70.)
180     .         ,tlay(jl,jkl)-tref)
181            tx2=tx*tx
182            phi(jl,ja)=at(1,ja)*tx+bt(1,ja)*tx2
183            psi(jl,ja)=at(2,ja)*tx+bt(2,ja)*tx2
184            phi(jl,ja)=exp(phi(jl,ja)/cst_voigt(2,ja))
185            psi(jl,ja)=exp(2.*psi(jl,ja))
186
187c                                        cumulative absorber(co2) amounts
188c                                        --------------------------------
189            co2_u(jl,ja,jkl)=co2_u(jl,ja,jkl+1)
190     .     +         pview/(10*g)*phi(jl,ja)*dp(jl,jkl)*co2c
191
192            co2_up(jl,ja,jkl)=co2_up(jl,ja,jkl+1)
193     .     +         pview/(10*g*2*pref)*psi(jl,ja)
194     .     *        (plev2(jl,jkl)-plev2(jl,jkl+1))*co2c
195
196
197c                                                  (aerosol) transmission
198c                                                  ----------------------
199c   on calcule directement les transmissions pour les aerosols.
200c   on multiplie le Qext  par 1-omega dans la bande du CO2.
201c   et pourquoi pas d'abord?  hourdin@lmd.ens.fr
202
203           aer_t(jl,ja,jkl)=exp(-pview*aer_a(jl,ja,jkl))
204         
205          enddo
206        enddo
207      enddo                 
208
209c----------------------------------------------------------------------
210      return
211      end
Note: See TracBrowser for help on using the repository browser.