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

Last change on this file since 3757 was 3757, checked in by emillour, 4 weeks ago

Mars PCM:
More code tidying: puting routines in modules and modernizing some old
constructs. Tested to not change results with respect to previous version.
EM

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