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

Last change on this file since 4033 was 4033, checked in by mmaurice, 29 hours ago

Generic PCM

Introduce molecular viscosity module. Add callconduc and callmolvis
flags to call conduction and molecular viscosity, respectively. Change
phitop in conduction to phitop_conduc (there is also phitop_molvis now).
Move thermal conductivity lambda to conc_mod.F90.

MM
File size: 9.2 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_conduc,zztop,a_coeff,s_coeff,force_conduction
12      use conc_mod,     only: lambda
13      use gases_h
14
15!=======================================================================
16!
17!   Molecular thermal conduction
18!   
19!   N. Descamp, F. Forget 05/1999
20!
21!=======================================================================
22
23!-----------------------------------------------------------------------
24!   Declarations:
25!-----------------------------------------------------------------------
26
27!   Arguments:
28!   ----------
29
30      INTEGER,intent(in) :: ngrid  ! number of atmospheric columns
31      INTEGER,intent(in) :: nlayer ! number of atmospheric layers
32      INTEGER,intent(in) :: nq     ! number of tracers
33      REAL,intent(in) :: ptimestep
34      REAL,intent(in) :: pplay(ngrid,nlayer)   ! pressure at middle of layers (Pa)
35      REAL,intent(in) :: pplev(ngrid,nlayer+1) ! (Pa)
36      REAL,intent(in) :: zzlay(ngrid,nlayer)   ! (m)
37      REAL,intent(in) :: zzlev(ngrid,nlayer+1) ! (m)
38      REAL,intent(in) :: zt(ngrid,nlayer)      ! Temperature [K]
39      REAL,intent(in) :: tsurf(ngrid)          ! Surface temperature [K]
40      REAL,intent(in) :: muvar(ngrid,nlayer+1) ! Molar mass (kg.mol-1)
41      REAL,intent(in) :: qvar(ngrid,nlayer,nq) ! Tracers (kg/kg)
42      LOGICAL,intent(in) :: firstcall ! Signals first call to physics.
43
44      REAL,intent(out) :: zdtconduc(ngrid,nlayer) ! [K.s-1]
45
46!   Local:
47!   ------
48
49      INTEGER :: i,ig,l,igas,kgas
50      REAL :: alpha(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            elseif(igas.eq.igas_CO2) then
114              !Interpolated from Hurly et al., (2007)
115              !valid between 20 and 1000 K (max 3 percent of error)
116              akk(igas) = 3.072e-4
117              skk(igas) = 0.69
118              molar_mass(igas) = 44.010e-3
119              akk_visc(igas) = 4.5054e-7
120              skk_visc(igas) = 0.6658
121              molar_frac(:,:,igas) = gfrac(igas)
122              here(igas) = .true.
123              ! Add more molecules here and the reference PLEASE!
124            else
125              akk(igas) = 0.0
126              skk(igas) = 0.0
127              molar_mass(igas) = 0.0
128              akk_visc(igas) = 0.0
129              skk_visc(igas) = 0.0
130              molar_frac(:,:,igas) = 0.0
131              here(igas) = .false.
132            endif
133          enddo
134        endif
135      endif
136     
137      if(force_conduction) then
138          lambda(:,1) = a_coeff*tsurf(:)**s_coeff / zzlay(:,1)
139          DO i = 2 , nlayer
140            lambda(:,i) = a_coeff*zt(:,i)**s_coeff / (zzlay(:,i)-zzlay(:,i-1))
141          ENDDO
142      else
143          ! Total conductivity is not equal to sum of conductivities of different species
144          ! We use the semi-empirical formulation from Mason et al., (1959)
145          mass_frac = 0.0
146          do igas=1,ngasmx
147            if(gfrac(igas).eq.-1) then
148              mass_frac(:,:,igas) = qvar(:,:,igas)
149            else
150              mass_frac(:,:,igas) = molar_frac(:,:,igas)*molar_mass(igas)/muvar(:,:)
151            endif
152          enddo
153          lambda_i(:,:,:) = 0.0
154          nu_i(:,:,:) = 0.0
155          somme(:,:) = 0.0
156          do igas=1,ngasmx
157            if(.not.(here(igas))) cycle
158            lambda_i(:,:,igas) = akk(igas)*zt(:,:)**skk(igas)
159            lambda_i(:,1,igas) = akk(igas)*tsurf(:)**skk(igas)
160            nu_i(:,:,igas) = akk_visc(igas)*zt(:,:)**skk_visc(igas)
161            nu_i(:,1,igas) = akk_visc(igas)*tsurf(:)**skk_visc(igas)
162          enddo
163          do igas=1,ngasmx
164            if(.not.(here(igas))) cycle
165            G_ik(:,:) = 0.0
166            do kgas=1,ngasmx
167              if(.not.(here(kgas))) cycle
168              if(kgas.ne.igas) then
169                G_ik(:,:) = G_ik(:,:) + (mass_frac(:,:,kgas)/mass_frac(:,:,igas))*(1.065/(2.*SQRT(2.)))*(1.+molar_mass(igas)/molar_mass(kgas))**(-1./2.)* &
170                      (1.+((nu_i(:,:,igas)*molar_mass(kgas)/(nu_i(:,:,kgas)*molar_mass(igas)))**(1./2.))*(molar_mass(igas)/molar_mass(kgas))**(1./4.))**2
171              endif
172            enddo
173            somme(:,:) = somme(:,:) + lambda_i(:,:,igas)/(1+G_ik(:,:))
174          enddo
175          lambda(:,1) = somme(:,1) / zzlay(:,1)
176          DO i = 2 , nlayer
177            lambda(:,i) = somme(:,i) / (zzlay(:,i)-zzlay(:,i-1))
178          ENDDO
179      endif
180
181      DO i=1,nlayer-1
182         muvol(:,i)=pplay(:,i)/(r*zt(:,i))
183         alpha(:,i)=cpp*(muvol(:,i)/ptimestep) &
184                             *(zlev(:,i+1)-zlev(:,i))
185      ENDDO
186
187      muvol(:,nlayer)=pplay(:,nlayer)/(r*zt(:,nlayer))
188      alpha(:,nlayer)=cpp*(muvol(:,nlayer)/ptimestep) &
189                             *(zlev(:,nlayer+1)-zlev(:,nlayer))
190
191!--------------------------------------------------------------------
192!
193!     Calculation of C and D coefficients
194!
195!-------------------------------------------------------------------
196
197      den(:,1)=alpha(:,1)+lambda(:,2)+lambda(:,1)
198      C(:,1)  =lambda(:,1)*(tsurf(:)-zt(:,1)) &
199              +lambda(:,2)*(zt(:,2)-zt(:,1))
200      C(:,1)  =C(:,1)/den(:,1)   
201      D(:,1)  =lambda(:,2)/den(:,1)           
202   
203      DO i = 2,nlayer-1
204         den(:,i)=alpha(:,i)+lambda(:,i+1)
205         den(:,i)=den(:,i)+lambda(:,i)*(1-D(:,i-1))
206           
207         C(:,i)  =lambda(:,i+1)*(zt(:,i+1)-zt(:,i)) &
208                 +lambda(:,i)*(zt(:,i-1)-zt(:,i)+C(:,i-1))   
209         C(:,i)  =C(:,i)/den(:,i)           
210
211         D(:,i)  =lambda(:,i+1) / den(:,i)
212      ENDDO
213
214      den(:,nlayer)=alpha(:,nlayer) + lambda(:,nlayer) &
215                   * (1-D(:,nlayer-1))
216      C(:,nlayer)  =C(:,nlayer-1)+zt(:,nlayer-1)-zt(:,nlayer)
217      C(:,nlayer)  =(C(:,nlayer)*lambda(:,nlayer)+phitop_conduc) &
218                   / den(:,nlayer)
219
220!----------------------------------------------------------------------
221!
222!     Calculation of new temperature pdt
223!
224!----------------------------------------------------------------------
225
226      DO i=1,nlayer
227         pdt(:,i)=0.
228      ENDDO
229      pdt(:,nlayer)=C(:,nlayer)
230      DO i=nlayer-1,1,-1
231         pdt(:,i)=C(:,i)+D(:,i)*pdt(:,i+1)
232      ENDDO
233!-----------------------------------------------------------------------
234!
235!     Calculation of zdtconduc
236!
237!-----------------------------------------------------------------------
238   
239      DO i=1,nlayer
240         zdtconduc(:,i)=pdt(:,i)/ptimestep
241      ENDDO
242
243      END subroutine conduction
244
245END MODULE conduction_mod
Note: See TracBrowser for help on using the repository browser.