MODULE conduction_mod IMPLICIT NONE CONTAINS SUBROUTINE conduction(ngrid,nlayer,nq,ptimestep,pplay,pplev,zt, & tsurf,zzlev,zzlay,muvar,qvar,firstcall,zdtconduc) use comcstfi_mod, only: r, cpp, mugaz use callkeys_mod, only: phitop,zztop,a_coeff,s_coeff,force_conduction use gases_h !======================================================================= ! ! Molecular thermal conduction ! ! N. Descamp, F. Forget 05/1999 ! !======================================================================= !----------------------------------------------------------------------- ! Declarations: !----------------------------------------------------------------------- ! Arguments: ! ---------- INTEGER,intent(in) :: ngrid ! number of atmospheric columns INTEGER,intent(in) :: nlayer ! number of atmospheric layers INTEGER,intent(in) :: nq ! number of tracers REAL,intent(in) :: ptimestep REAL,intent(in) :: pplay(ngrid,nlayer) ! pressure at middle of layers (Pa) REAL,intent(in) :: pplev(ngrid,nlayer+1) ! (Pa) REAL,intent(in) :: zzlay(ngrid,nlayer) ! (m) REAL,intent(in) :: zzlev(ngrid,nlayer+1) ! (m) REAL,intent(in) :: zt(ngrid,nlayer) ! Temperature [K] REAL,intent(in) :: tsurf(ngrid) ! Surface temperature [K] REAL,intent(in) :: muvar(ngrid,nlayer+1) ! Molar mass (kg.mol-1) REAL,intent(in) :: qvar(ngrid,nlayer,nq) ! Tracers (kg/kg) LOGICAL,intent(in) :: firstcall ! Signals first call to physics. REAL,intent(out) :: zdtconduc(ngrid,nlayer) ! [K.s-1] ! Local: ! ------ INTEGER :: i,ig,l,igas,kgas REAL :: alpha(ngrid,nlayer) REAL :: lambda(ngrid,nlayer) REAL :: muvol(ngrid,nlayer) ! kg.m-3 REAL :: C(ngrid,nlayer) REAL :: D(ngrid,nlayer) REAL :: den(ngrid,nlayer) REAL :: pdt(ngrid,nlayer) REAL :: zlev(ngrid,nlayer+1) REAL,ALLOCATABLE,SAVE :: akk(:) ! Akk conductivity coefficient for each species REAL,ALLOCATABLE,SAVE :: skk(:) ! skk conductivity coefficient for each species REAL,ALLOCATABLE,SAVE :: molar_mass(:) ! molar mass of each species (kg.m-3) REAL,ALLOCATABLE,SAVE :: akk_visc(:) ! Akk viscosity coefficient for each species REAL,ALLOCATABLE,SAVE :: skk_visc(:) ! skk viscosity coefficient for each species REAL,ALLOCATABLE,SAVE :: molar_frac(:,:,:) ! Molar fraction of each species (mol/mol) REAL :: mass_frac(ngrid,nlayer,ngasmx) ! Mass fraction of each species (kg/kg) REAL :: lambda_i(ngrid,nlayer,ngasmx) ! Conductivity of each species (W.m-1.K-1) REAL :: nu_i(ngrid,nlayer,ngasmx) ! Viscosity of each species (Pa.s) REAL :: G_ik(ngrid,nlayer) REAL :: somme(ngrid,nlayer) LOGICAL,ALLOCATABLE,SAVE :: here(:) !----------------------------------------------------------------------- ! The atmospheric conductivity is a function of temperature T : ! conductivity = Akk* T**skk !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! ! Calculation of alpha and lambda coefficients ! !----------------------------------------------------------------------- zlev(:,:) = zzlev(:,:) zlev(:,nlayer+1)= zzlev(:,nlayer)+zztop if(firstcall) then if(.not.force_conduction) then allocate(akk(ngasmx),skk(ngasmx),molar_mass(ngasmx),here(ngasmx)) allocate(akk_visc(ngasmx),skk_visc(ngasmx),molar_frac(ngrid,nlayer,ngasmx)) do igas=1,ngasmx if(igas.eq.igas_H2) then !Interpolated from Mehl et al., (2010) !valid between 20 and 1000 K (max 7 percent of error) !Equilibrium at local temperature (ortho/para effect takes into account) akk(igas) = 0.001539 skk(igas) = 0.8265 molar_mass(igas) = 2.01e-3 akk_visc(igas) = 1.2613e-7 skk_visc(igas) = 0.7445 molar_frac(:,:,igas) = gfrac(igas) here(igas) = .true. elseif(igas.eq.igas_He) then !Interpolated from Hurly et al., (2007) !valid between 20 and 1000 K (max 3 percent of error) akk(igas) = 0.003530 skk(igas) = 0.6655 molar_mass(igas) = 4.003e-3 akk_visc(igas) = 4.5054e-7 skk_visc(igas) = 0.6658 molar_frac(:,:,igas) = gfrac(igas) here(igas) = .true. ! Add more molecules here and the reference PLEASE! else akk(igas) = 0.0 skk(igas) = 0.0 molar_mass(igas) = 0.0 akk_visc(igas) = 0.0 skk_visc(igas) = 0.0 molar_frac(:,:,igas) = 0.0 here(igas) = .false. endif enddo endif endif if(force_conduction) then lambda(:,1) = a_coeff*tsurf(:)**s_coeff / zzlay(:,1) DO i = 2 , nlayer lambda(:,i) = a_coeff*zt(:,i)**s_coeff / (zzlay(:,i)-zzlay(:,i-1)) ENDDO else ! Total conductivity is not equal to sum of conductivities of different species ! We use the semi-empirical formulation from Mason et al., (1959) mass_frac = 0.0 do igas=1,ngasmx if(gfrac(igas).eq.-1) then mass_frac(:,:,igas) = qvar(:,:,igas) else mass_frac(:,:,igas) = molar_frac(:,:,igas)*molar_mass(igas)/muvar(:,:) endif enddo lambda_i(:,:,:) = 0.0 nu_i(:,:,:) = 0.0 somme(:,:) = 0.0 do igas=1,ngasmx if(.not.(here(igas))) cycle lambda_i(:,:,igas) = akk(igas)*zt(:,:)**skk(igas) lambda_i(:,1,igas) = akk(igas)*tsurf(:)**skk(igas) nu_i(:,:,igas) = akk_visc(igas)*zt(:,:)**skk_visc(igas) nu_i(:,1,igas) = akk_visc(igas)*tsurf(:)**skk_visc(igas) enddo do igas=1,ngasmx if(.not.(here(igas))) cycle G_ik(:,:) = 0.0 do kgas=1,ngasmx if(.not.(here(kgas))) cycle if(kgas.ne.igas) then G_ik(:,:) = G_ik(:,:) + (mass_frac(:,:,kgas)/mass_frac(:,:,igas))*(1.065/(2.*SQRT(2.)))*(1.+molar_mass(igas)/molar_mass(kgas))**(-1./2.)* & (1.+((nu_i(:,:,igas)*molar_mass(kgas)/(nu_i(:,:,kgas)*molar_mass(igas)))**(1./2.))*(molar_mass(igas)/molar_mass(kgas))**(1./4.))**2 endif enddo somme(:,:) = somme(:,:) + lambda_i(:,:,igas)/(1+G_ik(:,:)) enddo lambda(:,1) = somme(:,1) / zzlay(:,1) DO i = 2 , nlayer lambda(:,i) = somme(:,i) / (zzlay(:,i)-zzlay(:,i-1)) ENDDO endif DO i=1,nlayer-1 muvol(:,i)=pplay(:,i)/(r*zt(:,i)) alpha(:,i)=cpp*(muvol(:,i)/ptimestep) & *(zlev(:,i+1)-zlev(:,i)) ENDDO muvol(:,nlayer)=pplay(:,nlayer)/(r*zt(:,nlayer)) alpha(:,nlayer)=cpp*(muvol(:,nlayer)/ptimestep) & *(zlev(:,nlayer+1)-zlev(:,nlayer)) !-------------------------------------------------------------------- ! ! Calculation of C and D coefficients ! !------------------------------------------------------------------- den(:,1)=alpha(:,1)+lambda(:,2)+lambda(:,1) C(:,1) =lambda(:,1)*(tsurf(:)-zt(:,1)) & +lambda(:,2)*(zt(:,2)-zt(:,1)) C(:,1) =C(:,1)/den(:,1) D(:,1) =lambda(:,2)/den(:,1) DO i = 2,nlayer-1 den(:,i)=alpha(:,i)+lambda(:,i+1) den(:,i)=den(:,i)+lambda(:,i)*(1-D(:,i-1)) C(:,i) =lambda(:,i+1)*(zt(:,i+1)-zt(:,i)) & +lambda(:,i)*(zt(:,i-1)-zt(:,i)+C(:,i-1)) C(:,i) =C(:,i)/den(:,i) D(:,i) =lambda(:,i+1) / den(:,i) ENDDO den(:,nlayer)=alpha(:,nlayer) + lambda(:,nlayer) & * (1-D(:,nlayer-1)) C(:,nlayer) =C(:,nlayer-1)+zt(:,nlayer-1)-zt(:,nlayer) C(:,nlayer) =(C(:,nlayer)*lambda(:,nlayer)+phitop) & / den(:,nlayer) !---------------------------------------------------------------------- ! ! Calculation of new temperature pdt ! !---------------------------------------------------------------------- DO i=1,nlayer pdt(:,i)=0. ENDDO pdt(:,nlayer)=C(:,nlayer) DO i=nlayer-1,1,-1 pdt(:,i)=C(:,i)+D(:,i)*pdt(:,i+1) ENDDO !----------------------------------------------------------------------- ! ! Calculation of zdtconduc ! !----------------------------------------------------------------------- DO i=1,nlayer zdtconduc(:,i)=pdt(:,i)/ptimestep ENDDO END subroutine conduction END MODULE conduction_mod