! module radcool_module ! ! Radiative cooling: ! O3 at 9.6 microns ! CO2 at 15 microns (s.a., cco2gr.F) ! NO at 5.3 microns (see sub nocool in this module) ! O(3P) at 63 microns (see sub o3pcool in this module) ! use addfld_module,only: addfld implicit none real,parameter :: | dzcool = 0.25, ! delta zp resolution for cooling (includes tndown) | zpsrf = -21.5, ! zp at ground surface | zpm9 = -9., ! zp -9.0 level | zpm5 = -5., ! zp -5.0 level | zpm1925 = -19.25 ! zp -19.25 level integer,parameter :: | lmax = 59, ! number of wavelengths (for hco2,ho3, see vicool) | zpmax = 200 ! dimension for vars using nzpsrf_zpm5 real,parameter :: | a10=1.5988, | const = 2.63187e11 ! ! SGI compiler does not allow expressions in parameter initialization, ! so cannot dimension arrays with nzpsrf_zpm5 or nzpm9_zpm5. Therefore, ! the arrays are dimensioned zpmax (nzpxxxx < zpmax). integer :: | nzpsrf_zpm5, ! = int((zpm5-zpsrf)/dzcool)+1 ! 67 for tv,o3 | nzpm9_zpm5 ! = int((zpm5-zpm9 )/dzcool)+1 ! 17 for sn2,so2,etc ! contains !----------------------------------------------------------------------- subroutine radcool(tn,o2,ox,o1,n2,o3,no,co2,barm,xnmbar,cp, | cool_implicit,cool_explicit,cool_no, | lev0,lev1,lon0,lon1,lat) use params_module,only: dz,zibot,nlevp1 use cons_module,only: avo,p0,expz,expzmid,rmass_co2,grav, | rmassinv_o2,rmassinv_o1,rmassinv_co2,rmassinv_o3,rmassinv_n2 use solgar_module,only: | ntndown, | tndown ! (nlat,ntndown) ! ! Input args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: | tn, ! neutral temperature (f(:,nj+nt )) | o2, ! molecular oxygen (f(:,nj+nps )) | ox, ! ox (mmr) (f(:,nj+nps2)) | o1, ! o1 (mmr) (f(:,nj+npo1)) | n2, ! n2 (mmr) (f(:,nj+npn2)) | o3, ! o3 (mmr) (f(:,nj+npo3)) | no, ! no (mmr) (f(:,nj+npno)) | co2, ! co2(mmr) (f(:,nj+npco2)) | barm, ! mean mol weight (f(:,nj+nms)) | xnmbar, ! p0e(-z)*barm/(boltz*T) s3 | cp ! specific heat (f(:,nj+ncp)) ! ! Output args: real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) :: | cool_implicit, ! implicit cooling (NO) | cool_explicit, ! explicit cooling (NO) | cool_no ! NO cooling ! ! Local: integer :: i,k,kk,i0,i1 integer :: kzpm9,kzpm5,kzlb,kstep,kzpm1925 real,dimension(lev0:lev1,lon0:lon1) :: | vo2, | vo , | vn2, | vo3, | vco2, | den, | alam, | colco2, ! co2 column number density | cool_total,cool_total_deg, ! total cooling (diagnostic) | cool_co2,cool_o3,cool_o3p, ! co2, o3, o3p cooling | o3p_implicit, ! output by o3pcool | o3p_explicit ! output by o3pcool ! ! Cooling in degrees per day for secondary history diagnostics: real,dimension(lev0:lev1,lon0:lon1) :: | cooltot_deg, ! total cooling (deg/day) | coolco2_deg, ! 15 micron co2 | coolno_deg, ! 5.3 micron no | coolo3_deg, ! 9.6 micron o3 | coolo3p_deg ! 63 micron o3+ ! real,dimension(zpmax,lon0:lon1) :: | ztv , ! surface to zp -5 | zo3 , ! surface to zp -5 | zsn2 , ! zp -9 to zp -5 | zso2 , ! zp -9 to zp -5 | zo , ! zp -9 to zp -5 | zco2 , ! zp -9 to zp -5 | zuco2, ! zp -9 to zp -5 | al , ! zp -9 to zp -5 | am , ! zp -9 to zp -5 | zden , ! zp -9 to zp -5 | zhco2(lmax,lon0:lon1), ! set by vicool | zho3 (lmax,lon0:lon1), ! set by vicool | zflux(lon0:lon1) ! set by vicool real,dimension(lon0:lon1) :: y,zn2,zo2,zz,cppf real :: rko ! temperature dependent O-CO2 rate real :: factor,o3lb(ntndown+1) real,parameter :: ur=8.3144e+7 ! write(6,"('Enter radcool: lat=',i3)") lat ! o3lb=( | /0.028,0.034,0.061,0.106,0.166,0.239,0.321,0.428,0.523,0.661, | 1.051,1.606,2.335,3.223,4.201,5.190,6.113,6.925,7.646/) ! nzpsrf_zpm5 = int((zpm5-zpsrf)/dzcool)+1 ! 67 for tv,o3 nzpm9_zpm5 = int((zpm5-zpm9 )/dzcool)+1 ! 17 for sn2,so2,etc kzpm9 = int((zpm9-zibot)/dz)+1 ! col index at zp -9 from model lb kzpm5 = int((zpm5-zibot)/dz)+1 ! col index at zp -5 from model lb kzlb = int((zibot-zpsrf)/dzcool)+1 ! col index to model lb from surface kstep = int((dz+.01)/dzcool) ! 2 for dz 0.5, 1 for dz 0.25 kzpm1925 = int((zpm1925-zpsrf)/dzcool)+1 ! k at zp -19.25 from surface i0=lon0 ; i1=lon1 ! call addfld('XNMBAR',' ',' ',xnmbar(:,i0:i1), ! | 'lev',lev0,lev1,'lon',i0,i1,lat) ! call addfld('RC_CO2',' ',' ',co2,'lev',lev0,lev1,'lon',i0,i1,lat) do i=lon0,lon1 do k=lev0,lev1 vo2 (k,i) = o2 (k,i)*barm(k,i)*rmassinv_o2 vo (k,i) = o1 (k,i)*barm(k,i)*rmassinv_o1 vo3 (k,i) = o3 (k,i)*barm(k,i)*rmassinv_o3 vco2(k,i) = co2(k,i)*barm(k,i)*rmassinv_co2 vn2 (k,i) = n2 (k,i)*barm(k,i)*rmassinv_n2 den (k,i) = xnmbar(k,i)/barm(k,i) enddo ! k=lev0,lev1 enddo ! i=lon0,lon1 ! call addfld('VO2' ,' ',' ',vo2 ,'lev',lev0,lev1,'lon',i0,i1,lat) ! call addfld('VO' ,' ',' ',vo ,'lev',lev0,lev1,'lon',i0,i1,lat) ! call addfld('VO3' ,' ',' ',vo3 ,'lev',lev0,lev1,'lon',i0,i1,lat) ! call addfld('VCO2',' ',' ',vco2,'lev',lev0,lev1,'lon',i0,i1,lat) ! call addfld('VN2' ,' ',' ',vn2 ,'lev',lev0,lev1,'lon',i0,i1,lat) ! call addfld('DEN' ,' ',' ',den ,'lev',lev0,lev1,'lon',i0,i1,lat) ! ! Calculate co2 vertical column number density for reccurance relations: ! ! call addfld('CO2MMR' ,' ',' ',co2(:,i0:i1), ! | 'lev',lev0,lev1,'lon',i0,i1,lat) ! ! Upper boundary: factor = avo*p0*expz(lev1-1)*expzmid/(rmass_co2**2*grav) do i=lon0,lon1 colco2(lev1,i) = factor*.5*(co2(lev1-1,i)+co2(lev1,i))* | barm(lev1,i) enddo ! ! Integrate down to level 1: ! (loop order must be reversed from normal, i.e., i loop on inside): ! do k=lev1-1,lev0,-1 factor = avo*p0*expz(k)/(rmass_co2*grav)*dz do i=lon0,lon1 colco2(k,i) = colco2(k+1,i)+factor*co2(k,i) enddo ! i=lon0,lon1 ! These appear to be ok: ! write(6,"('radcool: lat=',i3,' k=',i3,' colco2(k,:)=',/, ! | (6e12.4))") lat,k,colco2(k,:) enddo ! k=lev0,lev1 ! call addfld('COLCO2' ,' ',' ',colco2, ! | 'lev',lev0,lev1,'lon',i0,i1,lat) do i=lon0,lon1 do kk=1,ntndown ! zp -21.5 to -17.25 by 0.25 zo3(kk,i) = o3lb(kk) ztv(kk,i) = tndown(lat,1) ! constant in vertical enddo ! kk=1,ntndown do kk=kzpm1925,ntndown ! zp -19.25 to -17.25 ztv(kk,i) = tndown(lat,kk) ! rewrite tv with tndown from zp -19.25 enddo ! kk=kzpm1925,ntndown do k=1,kzpm5 ! zp -17 (model lb) to zp -5 kk = kzlb+(k-1)*kstep ztv(kk,i) = tn(k,i) zo3(kk,i) = vo3(k,i) ! ! Interpolate to midpoints for dz=0.50 model: if (kstep > 1 .and. k > 1) then ! true for dz0.5 (kstep==2) ztv(kk-1,i) = ( tn(k,i) + tn(k-1,i))*0.5 zo3(kk-1,i) = (vo3(k,i) + vo3(k-1,i))*0.5 endif enddo ! k=1,kzpm5 ! do k=kzpm9,kzpm5 ! 17-25 for dz=0.25 model, 33-49 for dz=0.5 model kk = k-kzpm9+1 ! 1- 9 for dz=0.25 model, 1 -17 for dz=0.5 model if (kstep > 1) kk = kk*kstep-1 ! 1,3,5...17 for dz=0.5 model zsn2 (kk,i) = vn2 (k,i) zso2 (kk,i) = vo2 (k,i) zo (kk,i) = vo (k,i) zco2 (kk,i) = vco2 (k,i) zuco2(kk,i) = colco2(k,i) am (kk,i) = barm (k,i) zden (kk,i) = xnmbar(k,i)/barm(k,i) ! ! Interpolate to midpoints for dz=0.50 model: if (kstep > 1 .and. k > kzpm9) then ! true for dz=0.5 model zsn2 (kk-1,i) = (vn2 (k,i) + vn2 (k-1,i)) *0.5 zso2 (kk-1,i) = (vo2 (k,i) + vo2 (k-1,i)) *0.5 zo (kk-1,i) = (vo (k,i) + vo (k-1,i)) *0.5 zco2 (kk-1,i) = (vco2 (k,i) + vco2 (k-1,i)) *0.5 zuco2(kk-1,i) = (colco2(k,i) + colco2(k-1,i)) *0.5 am (kk-1,i) = (barm (k,i) + barm (k-1,i)) *0.5 zden (kk-1,i) = ((xnmbar(k,i)/barm(k,i))+ | (xnmbar(k-1,i)/barm(k-1,i)))*0.5 endif enddo ! k=kzpm9,kzpm5 enddo ! i=lon0,lon1 ! do k=1,nzpm9_zpm5 ! write(6,"('radcool: lat=',i3,' k=',i3,' zsn2=',/,(6e12.4))") ! | lat,k,zsn2(k,:) ! write(6,"('radcool: lat=',i3,' k=',i3,' zso2=',/,(6e12.4))") ! | lat,k,zso2(k,:) ! write(6,"('radcool: lat=',i3,' k=',i3,' zo=',/,(6e12.4))") ! | lat,k,zo(k,:) ! write(6,"('radcool: lat=',i3,' k=',i3,' zden=',/,(6e12.4))") ! | lat,k,zden(k,:) ! write(6,"('radcool: lat=',i3,' k=',i3,' zuco2=',/,(6e12.4))") ! | lat,k,zuco2(k,:) enddo ! subroutine recur(al,zuco2,nzpm9_zpm5,lev0,lev1,lon0,lon1,lat) ! integer,intent(in) :: nzpm9_zpm5,lev0,lev1,lon0,lon1,lat ! real,intent(in) :: ! | zuco2(nzpm9_zpm5,lon0:lon1) ! zp -9 to zp -5 ! real,intent(out) :: ! | al (nzpm9_zpm5,lon0:lon1) ! zp -9 to zp -5 ! do k=1,nzpm9_zpm5 ! write(6,"('Before recur: lat=',i3,' k=',i3, ! | ' zuco2(k,:) min,max=',2e12.4)") lat,k,minval(zuco2(k,:)), ! | maxval(zuco2(k,:)) ! enddo call recur(al(1:nzpm9_zpm5,:),zuco2(1:nzpm9_zpm5,:), | nzpm9_zpm5,lev0,lev1,lon0,lon1,lat) ! do k=1,nzpm9_zpm5 ! write(6,"('After recur: lat=',i3,' k=',i3, ! | ' al(k,:) min,max=',2e12.4)") lat,k,minval(al(k,:)), ! | maxval(al(k,:)) ! enddo call vicool( | ztv,zo3,zsn2,zso2,zo,zco2,al,am,zden, ! input to vicool | zflux, zhco2, zho3, ! output from vicool | lev0,lev1,lon0,lon1,lat) ! write(6,"('radcool: lat=',i3,' zflux mnmx=',2e12.4,' zhco2 mnmx=', ! | 2e12.4,' zho3 mnmx=',2e12.4)") lat,minval(zflux),maxval(zflux), ! | minval(zhco2),maxval(zhco2),minval(zho3),maxval(zho3) do k = 1,kzpm5 ! zp -17 to -5: 1->25 for dz=0.5, 1->49 for dz=0.25 kk = kzpm1925+k ! 11,12,13...59 for dz=0.25 model if (kstep > 1) | kk = kzpm1925+kstep*k-1 ! 11,13,15...59 for dz=0.5 model do i=lon0,lon1 cool_co2(k,i) = -zhco2(kk,i) ! s12 cool_o3 (k,i) = -zho3 (kk,i) ! s13 enddo ! i=lon0,lon1 enddo ! k=1,kzpm5 ! ! co2 cool to space approximation: ! alam = 0. ! whole array do k=kzpm5+1,lev1-1 do i=lon0,lon1 y(i) = tn(k,i)**(-1./3.) zn2(i) = 5.5e-17*sqrt(tn(k,i))+6.7e-10*exp(-83.8*y(i)) zo2(i) = 1.e-15*exp(23.37-230.9*y(i)+564.*y(i)*y(i)) rko = setrko(tn(k,i)) ! function setrko is in this module zz(i) = (vn2(k,i)*zn2(i)+vo2(k,i)*zo2(i)+vo(k,i)*rko)* | den(k,i) alam(k,i) = a10/(a10+zz(i)) cool_co2(k,i) = -const/barm(k,i)*vco2(k,i)*(1.-alam(k,i))* | (zflux(i)-exp(-960.217/tn(k,i))) cool_o3(k,i) = 0. ! s13 enddo ! i=lon0,lon1 enddo ! k=kzpm5+1,lev1-1 ! ! NO cooling: ! Sub nocooling returns cool_implicit, cool_explicit, and cool_no. ! call nocooling(no,tn,o1,o2,barm,cp,cool_co2,cool_o3, | cool_implicit,cool_explicit,cool_no, | lev0,lev1,lon0,lon1,lat) ! ! Add o(3p) cooling to total implicit and explicit terms. ! Sub o3pcool also returns o3+ implicit and explicit terms: ! call o3pcool(tn,ox,cp,cool_implicit,cool_explicit, | o3p_implicit,o3p_explicit,lev0,lev1,lon0,lon1,lat) ! ! Combine implicit and explicit terms for total and o3+ cooling: do i=lon0,lon1 do k=lev0,lev1-1 cool_total(k,i) = cool_explicit(k,i) / expz(k) + | cool_implicit(k,i)*cp(k,i)*tn(k,i) cool_o3p(k,i) = o3p_explicit(k,i) / expz(k) + | o3p_implicit(k,i)*cp(k,i)*tn(k,i) enddo ! k=lev0,lev1-1 enddo ! i=lon0,lon1 ! ! Save cooling terms on secondary history (ergs): ! call addfld('COOL_TOT','TOTAL COOLING','ERGS/GM/SEC', ! | cool_total(lev0:lev1-1,:),'lev',lev0,lev1-1,'lon',i0,i1,lat) ! call addfld('COOL_CO2','15 MICRON CO2 COOLING','ERGS/GM/SEC', ! | cool_co2(lev0:lev1-1,:),'lev',lev0,lev1-1,'lon',i0,i1,lat) ! call addfld('COOL_NO','5.3 MICRON NO COOLING','ERGS/GM/SEC', ! | cool_no(lev0:lev1-1,i0:i1),'lev',lev0,lev1-1,'lon',i0,i1,lat) ! call addfld('COOL_O3','9.6 MICRON O3 COOLING','ERGS/GM/SEC', ! | cool_o3(lev0:lev1-1,:),'lev',lev0,lev1-1,'lon',i0,i1,lat) ! call addfld('COOL_O3P','63 MICRON O3+ COOLING ','ERGS/GM/SEC', ! | cool_o3p(lev0:lev1-1,:),'lev',lev0,lev1-1,'lon',i0,i1,lat) ! ! Save cooling terms on secondary history (deg/day): ! cppf: conversion factor from ergs to deg/day ! do i=lon0,lon1 do k=lev0,lev1-1 cppf(i) = 86400./(0.5*ur*(7.*rmassinv_o2*vo2(k,i)+ | 7.*rmassinv_n2*vn2(k,i)+ | 5.*rmassinv_o1*vo (k,i))) cooltot_deg(k,i) = cool_total(k,i)*cppf(i) coolco2_deg(k,i) = cool_co2(k,i)*cppf(i) coolno_deg (k,i) = cool_no (k,i)*cppf(i) coolo3_deg (k,i) = cool_o3 (k,i)*cppf(i) coolo3p_deg(k,i) = cool_o3p(k,i)*cppf(i) enddo ! k=lev0,lev1-1 enddo ! i=lon0,lon1 ! ! call addfld('COOL_DEG','TOTAL COOLING','DEG/DAY', ! | cooltot_deg(lev0:lev1-1,:),'lev',lev0,lev1-1,'lon',i0,i1,lat) ! call addfld('CO2_DEG','15 MICRON CO2 COOLING','DEG/DAY', ! | coolco2_deg(lev0:lev1-1,:),'lev',lev0,lev1-1,'lon',i0,i1,lat) ! call addfld('NO_DEG','5.3 MICRON NO COOLING','DEG/DAY', ! | coolno_deg(lev0:lev1-1,:),'lev',lev0,lev1-1,'lon',i0,i1,lat) ! call addfld('O3_DEG','9.6 MICRON O3 COOLING','DEG/DAY', ! | coolo3_deg(lev0:lev1-1,:),'lev',lev0,lev1-1,'lon',i0,i1,lat) ! call addfld('O3P_DEG','63 MICRON O3+ COOLING ','DEG/DAY', ! | coolo3p_deg(lev0:lev1-1,:),'lev',lev0,lev1-1,'lon',i0,i1,lat) ! end subroutine radcool !----------------------------------------------------------------------- real function setrko(tn) ! ! Set temperature dependent O-CO2 rate (as per Vladimir Ogibalov, 11/8/02). ! V-T constant for co2(0110) quenching by O ! (measurements by Khvorostovskaya et al. 2002) ! This should be inlined. ! real,intent(in) :: tn ! if (tn < 260.) then setrko = 1.56e-12 elseif (tn >= 260..and.tn <= 300.) then setrko = (2.6-0.004*tn)*1.0e-12 else setrko = 1.4e-12 endif end function setrko !----------------------------------------------------------------------- subroutine recur(al,zuco2,nzpm9_zpm5,lev0,lev1,lon0,lon1,lat) use cco2gr_module,only: | cor150, cor360, cor540, cor720, uco2co, uco2ro, alo ! ! Args: integer,intent(in) :: nzpm9_zpm5,lev0,lev1,lon0,lon1,lat real,intent(in) :: | zuco2(nzpm9_zpm5,lon0:lon1) ! zp -9 to zp -5 real,intent(out) :: | al (nzpm9_zpm5,lon0:lon1) ! zp -9 to zp -5 ! ! Local: integer :: k,i,ii,nlons real :: co2int(4),uref(4) real,dimension(lon1-lon0+1) :: a, cor, uc ! do k=1,nzpm9_zpm5 ! write(6,"('recur: lat=',i3,' k=',i3,' zuco2=',/,(6e12.4))") ! | lat,k,zuco2(k,:) ! enddo ! k=1,nzpm9_zpm5 ! nlons = lon1-lon0+1 do k=1,6 co2int(1) = cor150(k) co2int(2) = cor360(k) co2int(3) = cor540(k) co2int(4) = cor720(k) ! uref(1) = uco2co(k)*150./360. uref(2) = uco2co(k) uref(3) = uco2co(k)*540./360. uref(4) = uco2co(k)*720./360. ! do i=lon0,lon1 ii = i-lon0+1 uc(ii) = zuco2(k,i) enddo ! i=lon0,lon1 ! SUBROUTINE A18LINV(X,Y,XN,YN,N,IMAX) ! integer,intent(in) :: n,imax ! real,intent(out) :: Y(IMAX) ! real,intent(in) :: X(IMAX), XN(N), YN(N) ! write(6,"('recur before a18linv: lat=',i3,' k=',i3,' uc=',/, ! | (6e12.4))") lat,k,uc call a18linv(uc,a,uco2ro,alo,51,nlons) ! a(nlons) is output ! write(6,"('recur: lat=',i3,' k=',i3,' a=',/,(6e12.4))") ! | lat,k,a call a18linv(uc,cor,uref,co2int,4,nlons) ! cor(nlons) is output ! write(6,"('recur: lat=',i3,' k=',i3,' cor=',/,(6e12.4))") ! | lat,k,cor do i=lon0,lon1 ii = i-lon0+1 al(k,i) = exp(cor(ii)+a(ii)) enddo ! i=lon0,lon1 enddo ! k=1,6 do k=7,nzpm9_zpm5 do i=lon0,lon1 ii = i-lon0+1 uc(ii) = zuco2(k,i) enddo call a18linv(uc,a,uco2ro,alo,51,nlons) ! a(nlons) is output ! write(6,"('recur: lat=',i3,' k=',i3,' a=',/,(6e12.4))") ! | lat,k,a do i=lon0,lon1 ii = i-lon0+1 al(k,i) = exp(a(ii)) enddo enddo ! do k=7,nzpm9_zpm5 ! do k=1,nzpm9_zpm5 ! write(6,"('recur: lat=',i3,' k=',i3,' zuco2(k,:)=',/, ! | (6e12.4))") lat,k,zuco2(k,:) ! write(6,"('recur: lat=',i3,' lon0,1=',2i3,' k=',i3,' al(k,:)=', ! | /,(6e12.4))") lat,lon0,lon1,k,al(k,:) ! enddo ! k=1,nzpm9_zpm5 end subroutine recur !----------------------------------------------------------------------- subroutine vicool(ztv,zo3,zsn2,zso2,zo,zco2,al,am,den, | zflux,zhco2,zho3,lev0,lev1,lon0,lon1,lat) use cco2gr_module,only: amat,bmat,ao3,ig use params_module,only: nlevp1 ! ! Args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat real,dimension(zpmax,lon0:lon1),intent(in) :: | ztv , ! surface to zp -5 | zo3 , ! surface to zp -5 | zsn2, ! zp -9 to zp -5 | zso2, ! zp -9 to zp -5 | zo , ! zp -9 to zp -5 | zco2, ! zp -9 to zp -5 | al , ! zp -9 to zp -5 | am , ! zp -9 to zp -5 | den ! zp -9 to zp -5 real,intent(out) :: | zflux(lon0:lon1), ! | zhco2(lmax,lon0:lon1), ! | zho3 (lmax,lon0:lon1) ! ! ! Local: integer :: i,k,j,js,ks,km real,dimension(zpmax,lon0:lon1) :: fu,fo3 real,dimension(nlevp1,lon0:lon1) :: alam real,dimension(lon0:lon1) :: h1,h2,h3,tt,zo2,zn2,y,zz, | aa1,aa2,d1,d2 real,parameter :: constb = 9.08795e9 real :: rko ! write(6,"('Enter vicool: lat=',i3,' lon0,1=',2i3,' nzpsrf_zpm5=', ! | i3,' nlevp1=',i3)") lat,lon0,lon1,nzpsrf_zpm5,nlevp1 ! ! Grid levels for height integration (p.s.h. distance = 0.25*IG) ! FU(nzpsrf_zpm5) - the LTE source functions for 15 um CO2 band at ! p.s.h.=0-16.5, ! FO3(45)- the same for 9.6 um O3 band at p.s.h. = 0-11, 0.25 ! AL(nzpm9_zpm5) - the quantum survival probability for ! p.s.h. = 12.5-16.5, 0.25 ! (for 15 um CO2 band only) do i=lon0,lon1 do k=1,nzpsrf_zpm5 fu (k,i) = exp(-960.217/ztv(k,i)) fo3(k,i) = exp(-1500. /ztv(k,i)) enddo ! k=1,nzpsrf_zpm5 enddo ! i=lon0,lon1 ! do k=1,nzpsrf_zpm5 ! write(6,"('vicool: lat=',i3,' k=',i3,' fu(k,:)=',/,(6e12.4))") ! | lat,k,fu(k,:) ! write(6,"('vicool: lat=',i3,' k=',i3,' fo3(k,:)=',/,(6e12.4))") ! | lat,k,fo3(k,:) ! enddo ! k=1,nzpsrf_zpm5 ! ! Calculate the heating rates for layer below s.h.p. = 12.5 ! 15 um CO2 + 9.6 um O3: ! Cooling rate in both o3 and co2 bands (x=2-10.5) matrix approach. ! ! cco2gr_module: real :: amat(43,9), bmat(43,9), ao3(35,9) ! do k=1,5 ks = k+8 do i=lon0,lon1 h2(i) = (amat(k,1)+bmat(k,1)*fu(ks,i))*fu(1,i) h3(i) = ao3(k,1)*fo3(1,i) enddo ! i=lon0,lon1 ! write(6,"('vicool: lat=',i3,' k=',i3,' h2=',/,(6e12.4))") ! | lat,k,h2 ! write(6,"('vicool: lat=',i3,' k=',i3,' h3=',/,(6e12.4))") ! | lat,k,h3 do j=3,9 js = ks+ig(j) do i=lon0,lon1 h2(i) = h2(i)+(amat(k,j)+bmat(k,j)*fu(ks,i))*fu(js,i) h3(i) = h3(i)+ao3(k,j)*fo3(js,i) enddo ! i=lon0,lon1 enddo ! j=3,9 do i=lon0,lon1 zhco2(k,i) = h2(i) zho3(k,i) = h3(i)*zo3(ks,i) enddo ! i=lon0,lon1 ! write(6,"('vicool: lat=',i3,' k=',i3,' zhco2=',/,(6e12.4))") ! | lat,k,zhco2(k,:) ! write(6,"('vicool: zho3=',/,(6e12.4))") zho3(k,:) enddo ! k=1,5 ! do k=6,18 ks = k+8 do i=lon0,lon1 h2(i) = (amat(k,1)+bmat(k,1)*fu(ks,i))*fu(1,i) h3(i) = ao3(k,1)*fo3(1,i) enddo ! i=lon0,lon1 do j=2,9 js = ks+ig(j) do i=lon0,lon1 h2(i) = h2(i)+(amat(k,j)+bmat(k,j)*fu(ks,i))*fu(js,i) h3(i) = h3(i)+ao3(k,j)*fo3(js,i) enddo ! i=lon0,lon1 enddo ! j=3,9 do i=lon0,lon1 zhco2(k,i) = h2(i) zho3(k,i) = h3(i)*zo3(ks,i) enddo ! i=lon0,lon1 ! write(6,"('vicool: lat=',i3,' k=',i3,' zhco2=',/,(6e12.4))") ! | lat,k,zhco2(k,:) ! write(6,"('vicool: lat=',i3,' k=',i3,' zho3(ks,:)=',/, ! | (6e12.4))") lat,k,zho3(k,:) enddo ! k=6,18 ! do k=19,35 ks = k+8 do i=lon0,lon1 h2(i) = 0. h3(i) = 0. enddo ! i=lon0,lon1 do j=1,9 js = ks+ig(j) do i=lon0,lon1 h2(i) = h2(i)+(amat(k,j)+bmat(k,j)*fu(ks,i))*fu(js,i) h3(i) = h3(i)+ao3(k,j)*fo3(js,i) enddo enddo ! j=1,9 do i=lon0,lon1 zhco2(k,i) = h2(i) zho3(k,i) = h3(i)*zo3(ks,i) enddo ! i=lon0,lon1 enddo ! k=19,35 ! ! Cooling rate in co2 bands (x=10.75-12.5, matrix approach) ! do k=36,43 ks = k+8 do i=lon0,lon1 h2(i) = 0. enddo ! i=lon0,lon1 do j=1,9 js = ks+ig(j) do i=lon0,lon1 h2(i) = h2(i)+(amat(k,j)+bmat(k,j)*fu(ks,i))*fu(js,i) enddo ! i=lon0,lon1 enddo ! j=1,9 do i=lon0,lon1 zhco2(k,i) = h2(i) zho3(k,i) = 0. enddo ! i=lon0,lon1 enddo ! k=36,43 ! ! Determine lambda(17) ! do k=1,nzpm9_zpm5 do i=lon0,lon1 ! ! co2-o2 and co2-n2 v-t constants. ! tt(i) = ztv(k+50,i) ! rko = setrko(tt(i)) ! function setrko is in this module ! y(i) = tt(i)**(-1./3.) zn2(i) = 5.5e-17*sqrt(tt(i))+6.7e-10*exp(-83.8*y(i)) zo2(i) = 1.e-15*exp(23.37-230.9*y(i)+564.*y(i)*y(i)) zz(i) = (zsn2(k,i)*zn2(i)+zso2(k,i)*zo2(i)+zo(k,i)*rko)* | den(k,i) alam(k,i) = a10/(a10+zz(i)) enddo ! i=lon0,lon1 enddo ! k=1,nzpm9_zpm5 ! do i=lon0,lon1 h1(i) = zhco2(43,i)/(zco2(1,i)*(1.-alam(1,i))*constb) enddo ! i=lon0,lon1 do k=2,nzpm9_zpm5 km = k-1 do i=lon0,lon1 aa1(i)= 1.-alam(km,i)*(1.-.25*al(k,i)-.75*al(km,i)) aa2(i)= 1.-alam(k ,i)*(1.-.75*al(k,i)-.25*al(km,i)) d1(i) = -.25*( al(k,i)+3.*al(km,i)) d2(i) = .25*(3.*al(k,i)+ al(km,i)) h2(i) = (aa1(i)*h1(i)-d1(i)*fu(km+50,i)-d2(i)*fu(k+50,i))/ | aa2(i) zhco2(k+42,i) = h2(i)*zco2(k,i)*(1.-alam(k,i))/am(k,i)*const zho3 (k+42,i) = 0. h1(i) = h2(i) enddo ! i=lon0,lon1 ! write(6,"('vicool: lat=',i3,' k=',i3,' km=',i3,' fu(km+50,:)=', ! | /,(6e12.4))") lat,k,km,fu(km+50,:) ! write(6,"('vicool: lat=',i3,' k=',i3,' km=',i3,' fu(k+50,:)=', ! | /,(6e12.4))") lat,k,km,fu(k+50,:) ! write(6,"('vicool: lat=',i3,' k=',i3,' aa1=',/,(6e12.4))") ! | lat,k,aa1 ! write(6,"('vicool: lat=',i3,' k=',i3,' aa2=',/,(6e12.4))") ! | lat,k,aa2 ! ! write(6,"('vicool: lat=',i3,' k=',i3,' km=',i3,' al(k,:)=',/, ! | (6e12.4))") lat,k,km,al(k,:) ! write(6,"('vicool: lat=',i3,' k=',i3,' km=',i3,' al(km,:)=',/, ! | (6e12.4))") lat,k,km,al(km,:) ! ! write(6,"('vicool: lat=',i3,' k=',i3,' d1=',/,(6e12.4))") ! | lat,k,d1 ! write(6,"('vicool: lat=',i3,' k=',i3,' d2=',/,(6e12.4))") ! | lat,k,d2 ! ! write(6,"('vicool: lat=',i3,' k=',i3,' h2=',/,(6e12.4))") ! wrong ! | lat,k,h2 ! write(6,"('vicool: lat=',i3,' k=',i3,' zco2=',/,(6e12.4))") ! | lat,k,zco2(k,:) ! write(6,"('vicool: lat=',i3,' k=',i3,' alam=',/,(6e12.4))") ! | lat,k,alam(k,:) ! write(6,"('vicool: lat=',i3,' k=',i3,' am=',/,(6e12.4))") ! | lat,k,am(k,:) ! write(6,"('vicool: lat=',i3,' k=',i3,' zhco2(k+42,:)=',/, ! | (6e12.4))") lat,k,zhco2(k+42,:) enddo ! k=2,nzpm9_zpm5 do i=lon0,lon1 zflux(i) = h2(i)+fu(nzpsrf_zpm5,i) enddo ! i=lon0,lon1 ! write(6,"('vicool: lat=',i3,' zflux=',/,(6e12.4))") lat,zflux ! do k=1,lmax ! write(6,"('vicool: lat=',i3,' k=',i3,' lmax=',i3,' zhco2=',/, ! | (6e12.4))") lat,k,lmax,zhco2(k,:) ! enddo ! k=1,lmax ! do k=36,lmax ! write(6,"('vicool: lat=',i3,' k=',i3,' lam0=36 lmax=',i3, ! | ' zhco2(k,:)=',/,(6e12.4))") lat,k,lmax,zhco2(k,:) ! enddo ! k=1,lmax end subroutine vicool !----------------------------------------------------------------------- subroutine nocooling(no,tn,o1,o2,barm,cp,cool_co2,cool_o3, | cool_implicit,cool_explicit,cool_no,lev0,lev1,lon0,lon1,lat) ! ! Calculate NO cooling on 3-d subdomain grids. ! See Kockarts,G, Nitric Oxide Cooling in the Terrestrial Thermosphere, ! GRL Vol. 7, No. 2, pages 137-140, February 1980 ! use cons_module,only: | avo, ! avogadro number 6.023e23 | p0, ! standard pressure 5.e-4 | expz, ! exp(-z) at midpoints (dimensioned nlev+1) | gask, ! gas constant 8.314e7 | rmassinv_o1, ! inverse of o1 mass (1./16.) | rmassinv_o2, ! inverse of o2 mass (1./32.) | rmassinv_no ! inverse of no mass (1./30.) ! ! Input args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: | no, ! 3d nitric oxide (mmr) | tn, ! 3d neutral temperature (deg K) | o1, ! 3d atomic oxygen (mmr) | o2, ! 3d molecular oxygen (mmr) | barm, ! mean molecular weight | cp ! specific heat real,dimension(lev0:lev1,lon0:lon1),intent(in) :: | cool_co2, cool_o3 ! ! Output args: real,dimension(lev0:lev1,lon0-2:lon1+2),intent(out) :: | cool_implicit, ! implicit cooling (NO) | cool_explicit, ! explicit cooling (NO) | cool_no ! NO cooling output (deg/day) ! ! Local: integer :: i,j,k real,dimension(lev0:lev1,lon0:lon1) :: | nocoolrc, ! NO cooling recombined from implicit and explicit terms | nocooldd ! NO cooling (deg/day) ! ! Longitude subdomain: do i=lon0,lon1 ! ! Pressure (full column, cool_no defined to lev1-1): do k=lev0,lev1-1 ! ! ergs/cm3/sec: ! (13.3 is the Einstein radiative transition probability) ! cool_no(k,i) = 0.5*(barm(k,i)+barm(k+1,i))*avo*p0* ! s11 ! | expz(k)/(gask*tn(k,i))*(2.7e-11*o1(k,i)*rmassinv_o1+ | expz(k)/(gask*tn(k,i))*(4.2e-11*o1(k,i)*rmassinv_o1+ ! | expz(k)/(gask*tn(k,i))*(1.0e-13*o1(k,i)*rmassinv_o1+ | 2.4e-14*o2(k,i)*rmassinv_o2) cool_no(k,i) = 4.956e-12*(avo*no(k,i)*rmassinv_no)* ! s11 | (cool_no(k,i)/(cool_no(k,i)+13.3))*exp(-2700./tn(k,i)) ! ! Implicit and explicit terms: ! cool_implicit(k,i) = (2700.*cool_no(k,i)+960.*cool_co2(k,i)+ ! f(:,ntik) | 1500.*cool_o3(k,i))/ | (.5*(cp(k,i)+cp(k+1,i))*tn(k,i)**2.) cool_explicit(k,i)= ((1.-2700./tn(k,i))*cool_no (k,i)+ ! f(:,ntek) | (1.-960. /tn(k,i))*cool_co2(k,i)+ | (1.-1500./tn(k,i))*cool_o3 (k,i))*expz(k) ! ! The following was provided to F. Sassi for waccm in July, 2003: ! noimp(k,i) = 2700.*cool_no(k,i)/ ! | (.5*(cp(k,i)+cp(k+1,i))*tn(k,i)**2.) ! noexp(k,i) = (1.-2700./tn(k,i))*cool_no(k,i)*expz(k) ! ! Recombine implicit and explicit terms (ergs/gm/sec): nocoolrc(k,i) = cool_explicit(k,i) / expz(k) + | cool_implicit(k,i) * cp(k,i) * tn(k,i) ! ! Convert from ergs to degrees/day: ! nocooldd(k,i) = cool_no(k,i)*86400./cp(k,i) enddo ! k=lev0,lev1 enddo ! i=lon0,lon1 ! end subroutine nocooling !----------------------------------------------------------------------- subroutine o3pcool(tn,ox,cp,cool_implicit,cool_explicit, | o3p_implicit,o3p_explicit,lev0,lev1,lon0,lon1,lat) use cons_module,only: avo,rmassinv_o1,expz ! ! Calculate implicit and explicit cooling terms for O(3P), and ! add to total cool_implicit and cool_explicit. ! (This was formerly subroutine newto3p, newto3p.F (tgcm24)) ! ! E.C.RIDLEY comment from newto3p.F (7/9/87): ! CHANGE O(3P) COOLING FACTOR FROM 0.5 TO 1.0 OF THE BATES EXPRESSION ! ! Input args: integer,intent(in) :: lev0,lev1,lon0,lon1,lat real,dimension(lev0:lev1,lon0-2:lon1+2),intent(in) :: | tn, ! neutral temperature (f(:,nj+nt )) | ox, ! ox (mmr) (f(:,nj+nps2)) | cp ! specific heat (f(:,ncp)) ! ! Inout args: real,dimension(lev0:lev1,lon0-2:lon1+2),intent(inout) :: | cool_implicit, ! output with o3p implicit cooling added | cool_explicit ! output with o3p explicit cooling added ! ! Output args: real,dimension(lev0:lev1,lon0:lon1),intent(out) :: | o3p_implicit, ! o3+ implicit cooling | o3p_explicit ! o3+ explicit cooling ! ! Local: integer :: k,i real,dimension(lev0:lev1,lon0:lon1) :: | work1, ! s11 ! work2 ! s12 real,parameter :: | an(3) = (/0.835E-18, 0.6, 0.2/), | bn(3) = (/228.,228.,325./) real :: xfac(lev0:lev1) ! ! This should be generalized for dz=0.5 and dz=0.25. xfac(lev0:11)=.01 xfac(12:18)=(/.05,.1,.2,.4,.55,.7,.75/) xfac(19:lev1) = .8 ! do i=lon0,lon1 do k=lev0,lev1-1 work1(k,i) = an(1)*.5*(xfac(k)+xfac(k+1))*avo*rmassinv_o1 ! s11 enddo ! k=lev0,lev1 enddo ! i=lon0,lon1 ! do i=lon0,lon1 do k=lev0,lev1-1 work1(k,i) = work1(k,i)*ox(k,i)*exp(-bn(1)/tn(k,i)) ! s11 work2(k,i) = 1.+an(2)*exp(-bn(2)/tn(k,i))+ ! s12 | an(3)*exp(-bn(3)/tn(k,i)) work1(k,i) = work1(k,i)/work2(k,i) work2(k,i) = work1(k,i)/work2(k,i)/tn(k,i)**2* | (bn(1)+(bn(1)-bn(3))*an(3)*exp(-bn(3)/tn(k,i))) ! ! Add implicit contributions to cool_implicit: ! cool_implicit(k,i) = cool_implicit(k,i)+work2(k,i)/ | (0.5*(cp(k,i)+cp(k+1,i))) o3p_implicit(k,i) = work2(k,i)/(0.5*(cp(k,i)+cp(k+1,i))) ! ! Add explicit contributions to cool_explicit: ! work1(k,i) = work1(k,i)-tn(k,i)*work2(k,i) cool_explicit(k,i) = cool_explicit(k,i)+work1(k,i)*expz(k) o3p_explicit(k,i) = work1(k,i)*expz(k) enddo ! k=lev0,lev1 enddo ! i=lon0,lon1 ! end subroutine o3pcool !----------------------------------------------------------------------- SUBROUTINE A18LINV(X,Y,XN,YN,N,IMAX) implicit none C **** C **** This procedure performs linear interpolation within the C **** table defined by the N points (XN(NN),Y(NN)). C **** Where: C **** C **** NN = 1,N,1 C **** C **** XN(NN) < XN(NN+1) for NN = 1,N-1 C **** C **** Parameters: C **** C **** X(IMAX) = array of IMAX x-values at which linear C **** interpolation is required C **** C **** XN(N) = array of N abscissae at which function values C **** are given C **** C **** YN(N) = function values corresponding to abscissae, C **** XN(N) C **** C **** Output: C **** C **** Y(IMAX) The IMAX interpolated values are C **** returned in this array C **** ! ! Args: integer,intent(in) :: n,imax real,intent(out) :: Y(IMAX) real,intent(in) :: X(IMAX), XN(N), YN(N) ! ! Local: integer :: KK(IMAX),i,nn C **** C **** Where: C **** Y(IMAX) is vector output C **** C **** KK is work space C **** C **** C **** Initialize array KK C **** DO I = 1,IMAX KK(I) = 0 ENDDO C **** C **** Locate interval in (XN,YN) in table containing X(I) C **** DO NN = 1,N-1 DO I = 1,IMAX ! KK(I) = merge(NN+1,KK(I),(XN(NN+1)-X(I))*(X(I)-XN(NN))>=0.) if ((xn(nn+1)-x(i))*(x(i)-xn(nn)) >= 0.) kk(i) = nn+1 ENDDO ENDDO C **** C **** Check for C **** C **** X(I) < XN(1), X(I) > X(N) C **** C **** and use linear extrapolation if necessary C **** DO I = 1,IMAX ! KK(I) = merge(2,KK(I),XN(1)-X(I)>=0.) ! KK(I) = merge(N,KK(I),X(I)-XN(N)>=0.) if (xn(1)-x(i) >= 0.) kk(i) = 2 if (x(i) >= xn(n)) kk(i) = n ENDDO C **** C **** Perform interpolation prescribed above C **** DO I = 1,IMAX ! 24-31:forrtl: severe (408): fort: (3): Subscript #1 of the array YN has value -1 which is less than the lower bound of 1 Y(I) = (YN(KK(I)-1)*(XN(KK(I))-X(I)) + YN(KK(I))* 1 (X(I)-XN(KK(I)-1)))/(XN(KK(I))-XN(KK(I)-1)) ENDDO END subroutine a18linv !----------------------------------------------------------------------- end module radcool_module