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

Last change on this file since 1112 was 1047, checked in by emillour, 11 years ago

Mars GCM:

  • IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
  • Also updated EOF routines, everything is now in eofdump_mod.F90
  • Removed unused routine lectfux.F (in dyn3d)

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