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

Last change on this file since 3567 was 1917, checked in by emillour, 7 years ago

Mars GCM:
Code cleanup: get rid of routine "zerophys".
EM

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