subroutine nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,pq,eps23,eps76) use comgeomfi_h implicit none !================================================================== ! Purpose ! ------- ! Calculation of epsilon - non LTE efficiency factor ! For CH4 heating and cooling rate calculation ! Cf Strobel et al 1996 ! ! Inputs ! ------ ! ngrid Number of vertical columns ! nlayer Number of layers ! pplay(ngrid,nlayer) Pressure layers ! pt ! ! Outputs ! ------- ! ! eps_sw ! efficiency factor for heating rate ! eps_lw ! efficiency factor for cooling rate ! ! Authors ! ------- ! Tanguy Bertrand (2015) ! !================================================================== #include "dimensions.h" #include "dimphys.h" #include "comcstfi.h" #include "surfdat.h" #include "comvert.h" #include "callkeys.h" #include "tracer.h" !----------------------------------------------------------------------- ! Arguments INTEGER ngrid, nlayer, nq REAL pplay(ngrid,nlayer) REAL pplev(ngrid,nlayer) REAL pt(ngrid,nlayer) REAL pq(ngrid,nlayer,nq) !----------------------------------------------------------------------- ! Local variables INTEGER l,ig REAL temper ! atmospheric temperature REAL cbol ! Boltzman constant REAL avog ! avogadro constant REAL beta23 REAL beta33 REAL nd(ngrid,nlayer) ! density number cm-3 REAL na(ngrid,nlayer) ! column density number cm-2 REAL rho(ngrid,nlayer) ! density of atmosphere REAL kn3,k3,k2n4,k3n4,k4 REAL An3,A3,A2n4,A3n4,A4 REAL dphin3,dphi3,dphi2n4,dphi3n4,dphi4 REAL phin3(ngrid,nlayer) REAL phi3(ngrid,nlayer) REAL phi2n4(ngrid,nlayer) REAL phi3n4(ngrid,nlayer) REAL phi4(ngrid,nlayer) ! DATA from Strobel et al 1996 REAL, DIMENSION(12) :: gam2 ! band escape fonction 2.3 REAL, DIMENSION(12) :: Naval2 ! column density number 2.3 REAL, DIMENSION(14) :: gam3 ! band escape fonction 3.3 REAL, DIMENSION(14) :: Naval3 ! column density number 3.3 REAL, DIMENSION(20) :: Navalcool REAL, DIMENSION(20) :: gamcool REAL gammafi23(nlayer),gammafi33(nlayer) ! Band espace fonction ! Output REAL eps23(ngrid,nlayer) REAL eps33(ngrid,nlayer) REAL eps76(ngrid,nlayer) ! Option simplified coefficient logical, save :: strobel=.true. !if false, coef are calculated from linear fit REAL klw, ksw, eps_sw1Pa !----------------------------------------------------------------------- IF (strobel) then ! Constantes cbol=1.3806488e-23 !Boltzman constant avog = 6.022141e23 temper=100. ! cte Temperature atmosphere K ! Calculation phi ! A2=25.8 !Einstein coefficient sec-1 ! P3=7.e-3 ! collisional probability ! band 3 k3=2.8e-11 !cm-3 A3=25.2 dphi3=k3/A3 ! band n3 kn3=2.8e-11 !cm-3 An3=27.3 dphin3=kn3/An3 ! band 3n4 k3n4=1.16e-11 !cm-3 A3n4=2.12 dphi3n4=k3n4/(3*A3n4) ! band 4 k4=4.5e-16 !cm-3 A4=2.12 dphi4=k4/A4 ! band 2n4 k2n4=1.16e-11 A2n4=2.12 dphi2n4=k2n4/(2*A2n4) ! 2.3 band : beta23=0.09 ! 3.3 band : beta23=0.13 ! Escape Function Gamma (Data Strobel et al) : ! 100 K 2.3 Naval2=(/1.e10,1.e16,8.e16,4.e17,2.e18,1.e19,5.e19,2.4e20,1.3e21,6.e21,1.e23,1.e24/) gam2=(/1.,1.,0.95,0.9,0.65,0.38,0.14,4.e-2,6.e-3,1.1e-3,6.e-5,5.e-6/) ! 40 K 2.3 ! gam2=(/1.,1.,0.95,0.7,0.47,0.16,0.05,7.e-3,1.2e-3,2.e-4,1.e-5,7.e-7/) ! 100 K 3.3 Naval3=(/1e10,1e15,1.4e15,3e15,1.5e16,7e16,4e17,2e18,1e19,5e19,2.7e20,1.3e21,7e21,1e24/) gam3=(/1.,1.,0.95,0.9,0.7,0.3,0.1,0.04,0.015,5.e-3,1.5e-3,4.e-4,8.e-5,1.e-6/) ! gam=(/1.,1.,1.,1.,0.95,0.9,0.8,0.7,0.6,0.5,0.2,0.12,0.07,0.03,0.014,0.006,0.0026,0.0007,1.e-4,1.e-5/) ! values of gamma ! Naval=(/1e10,1.2e15,1e16,2e16,2e17,4e17,9e17,2e18,4e18,1e19,2e19,5e19,1.3e20,2.4e20,6e20,1.4e20,3e21,1e22,1e23,1e24/) ! values of Na in cm-2 ! gamcool=(/8.79e15,8.79e15,8.79e15,1.87e16,1.94e17,3.79e17,8.04e17,1.62e18,2.92e18,6.22e18,9.72e18,1.45e19,2.21e19, 2.76e19,3.55e19,4.01e19,5.24e19,6.39e19,9.99e19,1.49e20/) ! Navalcool=(/1e10,1.12e16,1.12e16, 3.12e16, 2.31e17, 6.31e17,1.53e18, 3.53e18, 7.53e18, 1.75e19, 3.75e19,8.75e19, 2.175e20, 4.575e20, 1.057e21, 1.197e21,4.19e21, 1.419e22, 1.14e23,1.114e24/) ! gamcool=gam DO ig=1,ngrid eps23(ig,:)=0. eps33(ig,:)=0. eps76(ig,:)=0. phi3n4(ig,:)=0. phi2n4(ig,:)=0. phin3(ig,:)=0. phi4(ig,:)=0. nd(ig,:)=0. na(ig,:)=0. DO l=1,nlayer ! old !! nd(ig,l)=pplay(ig,l)/(cbol*temper) !! nd(ig,l)=pplay(ig,l)/(cbol*pt(ig,l)) ! number of mol density of the layer !! na(ig,l)=na(ig,l)+nd(ig,l) ! density : sum of each layer ! calculation of density ! rho(ig,l)=pplay(ig,l)/(r*temper) rho(ig,l)=pplay(ig,l)/(r*pt(ig,l)) ! calculation of column number density for ch4 (cm-2) na(ig,l)=pq(ig,l,igcm_ch4_gas)/(mmol(igcm_ch4_gas)*1e-3)* & pplev(ig,l)/g*avog*1.e-4 ! calculation of atmospheric number density (cm-3) for each layer nd(ig,l)=rho(ig,l)/(mmol(igcm_n2)*1.e-3)*avog*1.e-6 ! calculation of the phi : phin3(ig,l)=dphin3*nd(ig,l) phi3(ig,l)=dphi3*nd(ig,l) phi2n4(ig,l)=dphi2n4*nd(ig,l) phi3n4(ig,l)=dphi3n4*nd(ig,l) phi4(ig,l)=dphi4*nd(ig,l) ! prendre en compte Temp ! diff entre les gammaf4 et f3 ? ENDDO ! interpolation of na values in cm-2 from Strobel values : CALL interp_line(Naval2,gam2,12,na(ig,:),gammafi23,nlayer) CALL interp_line(Naval3,gam3,14,na(ig,:),gammafi33,nlayer) ! CALL interp_line(Naval,gam,20,na(ig,:),gammafi,nlayer) ! CALL interp_line(Naval,gamcool,20,na(ig,:),gammacool,nlayer) ! CALL interp_line(Navalcool,gam,20,na(ig,:),gammacool,nlayer) ! write(*,*) 'TBrad phi4 = ',phi4(1,:) ! write(*,*) 'TBrad gammafi = ',gammafi DO l=1,nlayer ! Heating : cf Strobel 1996 and Zalucha 2011 ! eps 2.3 eps23(ig,l)=phin3(ig,l)/(phin3(ig,l)+gammafi23(l))* & (beta23+phi3n4(ig,l)/(phi3n4(ig,l)+1)* & phi4(ig,l)/(phi4(ig,l)+gammafi23(l))* & (3.*phi2n4(ig,l))/(1.+phi2n4(ig,l))*(1-beta23)/3.) ! eps 3.3 eps33(ig,l)=phi3(ig,l)/(phi3(ig,l)+gammafi33(l))* & (beta33+phi2n4(ig,l)/(phi2n4(ig,l)+1)* & phi4(ig,l)/(phi4(ig,l)+gammafi33(l))* & (1-beta33)) ! Cooling from Zalucha et al 2013 eq10 eps76(ig,l)=1./(1.+gammafi23(l)*10./(2.*phi4(ig,l))) ENDDO ENDDO ! ngrid ELSE ! if (strobel) ! if strobel = false klw = 3.6E-4 ! fit strobel ! eps_sw1Pa = 0.6 ! fit strobel eps_sw1Pa = 0.90 ! to fit NH temperature en 1D [ch4]=0.5 ! eps_sw1Pa = 0.98 ! to fit NH temperature en 3D = fit 1D [ch4]=0.25 ksw = 8.5E-5*(0.9/(eps_sw1Pa -0.1) -1.) do l=1, nlayer do ig=1, ngrid rho(ig,l)=pplay(ig,l)/(r*pt(ig,l)) eps23(ig,l) = 0.1+ 0.9/(1+ ksw/rho(ig,l)) eps76(ig,l) = 1/(1+klw/rho(ig,l)) end do end do END IF end subroutine nlte_ch4