source: trunk/LMDZ.GENERIC/libf/phystd/conduction.F90

Last change on this file was 3235, checked in by gmilcareck, 9 months ago

Add thermal conduction from thermosphere.
File size: 8.8 KB
Line 
1MODULE conduction_mod
2
3IMPLICIT NONE
4
5CONTAINS
6
7       SUBROUTINE conduction(ngrid,nlayer,nq,ptimestep,pplay,pplev,zt, &
8                        tsurf,zzlev,zzlay,muvar,qvar,firstcall,zdtconduc)
9   
10      use comcstfi_mod, only: r, cpp, mugaz
11      use callkeys_mod, only: phitop,zztop,a_coeff,s_coeff,force_conduction
12      use gases_h
13
14!=======================================================================
15!
16!   Molecular thermal conduction
17!   
18!   N. Descamp, F. Forget 05/1999
19!
20!=======================================================================
21
22!-----------------------------------------------------------------------
23!   Declarations:
24!-----------------------------------------------------------------------
25
26!   Arguments:
27!   ----------
28
29      INTEGER,intent(in) :: ngrid  ! number of atmospheric columns
30      INTEGER,intent(in) :: nlayer ! number of atmospheric layers
31      INTEGER,intent(in) :: nq     ! number of tracers
32      REAL,intent(in) :: ptimestep
33      REAL,intent(in) :: pplay(ngrid,nlayer)   ! pressure at middle of layers (Pa)
34      REAL,intent(in) :: pplev(ngrid,nlayer+1) ! (Pa)
35      REAL,intent(in) :: zzlay(ngrid,nlayer)   ! (m)
36      REAL,intent(in) :: zzlev(ngrid,nlayer+1) ! (m)
37      REAL,intent(in) :: zt(ngrid,nlayer)      ! Temperature [K]
38      REAL,intent(in) :: tsurf(ngrid)          ! Surface temperature [K]
39      REAL,intent(in) :: muvar(ngrid,nlayer+1) ! Molar mass (kg.mol-1)
40      REAL,intent(in) :: qvar(ngrid,nlayer,nq) ! Tracers (kg/kg)
41      LOGICAL,intent(in) :: firstcall ! Signals first call to physics.
42
43      REAL,intent(out) :: zdtconduc(ngrid,nlayer) ! [K.s-1]
44
45!   Local:
46!   ------
47
48      INTEGER :: i,ig,l,igas,kgas
49      REAL :: alpha(ngrid,nlayer)
50      REAL :: lambda(ngrid,nlayer)
51      REAL :: muvol(ngrid,nlayer)   ! kg.m-3
52      REAL :: C(ngrid,nlayer)
53      REAL :: D(ngrid,nlayer)
54      REAL :: den(ngrid,nlayer)
55      REAL :: pdt(ngrid,nlayer)
56      REAL :: zlev(ngrid,nlayer+1)
57      REAL,ALLOCATABLE,SAVE :: akk(:)        ! Akk conductivity coefficient for each species
58      REAL,ALLOCATABLE,SAVE :: skk(:)        ! skk conductivity coefficient for each species
59      REAL,ALLOCATABLE,SAVE :: molar_mass(:) ! molar mass of each species (kg.m-3)
60      REAL,ALLOCATABLE,SAVE :: akk_visc(:)   ! Akk viscosity coefficient for each species
61      REAL,ALLOCATABLE,SAVE :: skk_visc(:)   ! skk viscosity coefficient for each species
62      REAL,ALLOCATABLE,SAVE :: molar_frac(:,:,:) ! Molar fraction of each species (mol/mol)
63      REAL :: mass_frac(ngrid,nlayer,ngasmx) ! Mass fraction of each species (kg/kg)
64      REAL :: lambda_i(ngrid,nlayer,ngasmx)  ! Conductivity of each species (W.m-1.K-1)
65      REAL :: nu_i(ngrid,nlayer,ngasmx)      ! Viscosity of each species (Pa.s)
66      REAL :: G_ik(ngrid,nlayer)
67      REAL :: somme(ngrid,nlayer)
68      LOGICAL,ALLOCATABLE,SAVE :: here(:)
69
70
71!-----------------------------------------------------------------------
72!     The atmospheric conductivity is a function of temperature T :
73!      conductivity = Akk* T**skk
74!-----------------------------------------------------------------------
75     
76
77!-----------------------------------------------------------------------
78!
79!   Calculation of alpha and lambda coefficients
80!
81!-----------------------------------------------------------------------
82
83     
84      zlev(:,:) = zzlev(:,:)
85      zlev(:,nlayer+1)= zzlev(:,nlayer)+zztop
86     
87      if(firstcall) then
88        if(.not.force_conduction) then
89          allocate(akk(ngasmx),skk(ngasmx),molar_mass(ngasmx),here(ngasmx))
90          allocate(akk_visc(ngasmx),skk_visc(ngasmx),molar_frac(ngrid,nlayer,ngasmx))
91          do igas=1,ngasmx
92            if(igas.eq.igas_H2) then
93              !Interpolated from Mehl et al., (2010)
94              !valid between 20 and 1000 K (max 7 percent of error)
95              !Equilibrium at local temperature (ortho/para effect takes into account)
96              akk(igas) = 0.001539
97              skk(igas) = 0.8265
98              molar_mass(igas) = 2.01e-3
99              akk_visc(igas) = 1.2613e-7
100              skk_visc(igas) = 0.7445
101              molar_frac(:,:,igas) = gfrac(igas)
102              here(igas) = .true.
103            elseif(igas.eq.igas_He) then
104              !Interpolated from Hurly et al., (2007)
105              !valid between 20 and 1000 K (max 3 percent of error)
106              akk(igas) = 0.003530
107              skk(igas) = 0.6655
108              molar_mass(igas) = 4.003e-3
109              akk_visc(igas) = 4.5054e-7
110              skk_visc(igas) = 0.6658
111              molar_frac(:,:,igas) = gfrac(igas)
112              here(igas) = .true.
113              ! Add more molecules here and the reference PLEASE!
114            else
115              akk(igas) = 0.0
116              skk(igas) = 0.0
117              molar_mass(igas) = 0.0
118              akk_visc(igas) = 0.0
119              skk_visc(igas) = 0.0
120              molar_frac(:,:,igas) = 0.0
121              here(igas) = .false.
122            endif
123          enddo
124        endif
125      endif
126     
127      if(force_conduction) then
128          lambda(:,1) = a_coeff*tsurf(:)**s_coeff / zzlay(:,1)
129          DO i = 2 , nlayer
130            lambda(:,i) = a_coeff*zt(:,i)**s_coeff / (zzlay(:,i)-zzlay(:,i-1))
131          ENDDO
132      else
133          ! Total conductivity is not equal to sum of conductivities of different species
134          ! We use the semi-empirical formulation from Mason et al., (1959)
135          mass_frac = 0.0
136          do igas=1,ngasmx
137            if(gfrac(igas).eq.-1) then
138              mass_frac(:,:,igas) = qvar(:,:,igas)
139            else
140              mass_frac(:,:,igas) = molar_frac(:,:,igas)*molar_mass(igas)/muvar(:,:)
141            endif
142          enddo
143          lambda_i(:,:,:) = 0.0
144          nu_i(:,:,:) = 0.0
145          somme(:,:) = 0.0
146          do igas=1,ngasmx
147            if(.not.(here(igas))) cycle
148            lambda_i(:,:,igas) = akk(igas)*zt(:,:)**skk(igas)
149            lambda_i(:,1,igas) = akk(igas)*tsurf(:)**skk(igas)
150            nu_i(:,:,igas) = akk_visc(igas)*zt(:,:)**skk_visc(igas)
151            nu_i(:,1,igas) = akk_visc(igas)*tsurf(:)**skk_visc(igas)
152          enddo
153          do igas=1,ngasmx
154            if(.not.(here(igas))) cycle
155            G_ik(:,:) = 0.0
156            do kgas=1,ngasmx
157              if(.not.(here(kgas))) cycle
158              if(kgas.ne.igas) then
159                G_ik(:,:) = G_ik(:,:) + (mass_frac(:,:,kgas)/mass_frac(:,:,igas))*(1.065/(2.*SQRT(2.)))*(1.+molar_mass(igas)/molar_mass(kgas))**(-1./2.)* &
160                      (1.+((nu_i(:,:,igas)*molar_mass(kgas)/(nu_i(:,:,kgas)*molar_mass(igas)))**(1./2.))*(molar_mass(igas)/molar_mass(kgas))**(1./4.))**2
161              endif
162            enddo
163            somme(:,:) = somme(:,:) + lambda_i(:,:,igas)/(1+G_ik(:,:))
164          enddo
165          lambda(:,1) = somme(:,1) / zzlay(:,1)
166          DO i = 2 , nlayer
167            lambda(:,i) = somme(:,i) / (zzlay(:,i)-zzlay(:,i-1))
168          ENDDO
169      endif
170
171      DO i=1,nlayer-1
172         muvol(:,i)=pplay(:,i)/(r*zt(:,i))
173         alpha(:,i)=cpp*(muvol(:,i)/ptimestep) &
174                             *(zlev(:,i+1)-zlev(:,i))
175      ENDDO
176
177      muvol(:,nlayer)=pplay(:,nlayer)/(r*zt(:,nlayer))
178      alpha(:,nlayer)=cpp*(muvol(:,nlayer)/ptimestep) &
179                             *(zlev(:,nlayer+1)-zlev(:,nlayer))
180
181!--------------------------------------------------------------------
182!
183!     Calculation of C and D coefficients
184!
185!-------------------------------------------------------------------
186
187      den(:,1)=alpha(:,1)+lambda(:,2)+lambda(:,1)
188      C(:,1)  =lambda(:,1)*(tsurf(:)-zt(:,1)) &
189              +lambda(:,2)*(zt(:,2)-zt(:,1))
190      C(:,1)  =C(:,1)/den(:,1)   
191      D(:,1)  =lambda(:,2)/den(:,1)           
192   
193      DO i = 2,nlayer-1
194         den(:,i)=alpha(:,i)+lambda(:,i+1)
195         den(:,i)=den(:,i)+lambda(:,i)*(1-D(:,i-1))
196           
197         C(:,i)  =lambda(:,i+1)*(zt(:,i+1)-zt(:,i)) &
198                 +lambda(:,i)*(zt(:,i-1)-zt(:,i)+C(:,i-1))   
199         C(:,i)  =C(:,i)/den(:,i)           
200
201         D(:,i)  =lambda(:,i+1) / den(:,i)
202      ENDDO
203
204      den(:,nlayer)=alpha(:,nlayer) + lambda(:,nlayer) &
205                   * (1-D(:,nlayer-1))
206      C(:,nlayer)  =C(:,nlayer-1)+zt(:,nlayer-1)-zt(:,nlayer)
207      C(:,nlayer)  =(C(:,nlayer)*lambda(:,nlayer)+phitop) &
208                   / den(:,nlayer)
209
210!----------------------------------------------------------------------
211!
212!     Calculation of new temperature pdt
213!
214!----------------------------------------------------------------------
215
216      DO i=1,nlayer
217         pdt(:,i)=0.
218      ENDDO
219      pdt(:,nlayer)=C(:,nlayer)
220      DO i=nlayer-1,1,-1
221         pdt(:,i)=C(:,i)+D(:,i)*pdt(:,i+1)
222      ENDDO
223!-----------------------------------------------------------------------
224!
225!     Calculation of zdtconduc
226!
227!-----------------------------------------------------------------------
228   
229      DO i=1,nlayer
230         zdtconduc(:,i)=pdt(:,i)/ptimestep
231      ENDDO
232
233      END subroutine conduction
234
235END MODULE conduction_mod
Note: See TracBrowser for help on using the repository browser.