subroutine setspi !================================================================== ! ! Purpose ! ------- ! Set up spectral intervals and Planck function in the longwave. ! ! Authors ! ------- ! Adapted from setspi in the NASA Ames radiative code by ! Robin Wordsworth (2009). ! ! Called by ! --------- ! callcorrk.F ! ! Calls ! ----- ! none ! !================================================================== use radinc_h, only: L_NSPECTI,corrkdir,banddir,NTstar,NTstop use radcommon_h, only: BWNI,BLAMI,WNOI,DWNI,WAVEI,planckir,sigma use datafile_mod, only: datadir implicit none #include "callkeys.h" #include "comcstfi.h" logical file_ok integer nw, nt, m, mm, file_entries real*8 a, b, ans, y, bpa, bma, T character(len=30) :: temp1 character(len=100) :: file_id character(len=100) :: file_path ! C1 and C2 values from Goody and Yung (2nd edition) MKS units ! These values lead to a "sigma" (sigma*T^4) of 5.67032E-8 W m^-2 K^-4 real*8 :: c1 = 3.741832D-16 ! W m^-2 real*8 :: c2 = 1.438786D-2 ! m K real*8 :: lastband(2), plancksum logical forceEC, planckcheck real*8 :: x(12) = [ -0.981560634246719D0, -0.904117256370475D0, & -0.769902674194305D0, -0.587317954286617D0, & -0.367831498998180D0, -0.125233408511469D0, & 0.125233408511469D0, 0.367831498998180D0, & 0.587317954286617D0, 0.769902674194305D0, & 0.904117256370475D0, 0.981560634246719D0 ] real*8 :: w(12) = [ 0.047175336386512D0, 0.106939325995318D0, & 0.160078328543346D0, 0.203167426723066D0, & 0.233492536538355D0, 0.249147045813403D0, & 0.249147045813403D0, 0.233492536538355D0, & 0.203167426723066D0, 0.160078328543346D0, & 0.106939325995318D0, 0.047175336386512D0 ] mm=0 forceEC=.false. planckcheck=.true. !======================================================================= ! Set up spectral bands - wavenumber [cm^(-1)]. Go from smaller to ! larger wavenumbers. write(temp1,'(i2.2)') L_NSPECTI !file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/narrowbands_IR.in' file_id='/corrk_data/'//trim(adjustl(banddir))//'/narrowbands_IR.in' file_path=TRIM(datadir)//TRIM(file_id) ! check that the file exists inquire(FILE=file_path,EXIST=file_ok) if(.not.file_ok) then write(*,*)'The file ',TRIM(file_path) write(*,*)'was not found by setspi.F90, exiting.' write(*,*)'Check that your path to datagcm:',trim(datadir) write(*,*)' is correct. You can change it in callphys.def with:' write(*,*)' datadir = /absolute/path/to/datagcm' write(*,*)'Also check that the corrkdir you chose in callphys.def exists.' call abort endif ! check that the file contains the right number of bands call system('wc -l '//file_path//' > bandlen.txt') open(131,file='bandlen.txt', form='formatted') read(131,*) file_entries close(131) call system('rm -f bandlen.txt') if(file_entries-1.ne.L_NSPECTI) then write(*,*) 'setspi: L_NSPECTI = ',L_NSPECTI, 'in the model, but there are ' write(*,*) file_entries-1,'entries in ', & TRIM(file_path),', exiting.' call abort endif ! load and display the data open(111,file=file_path,form='formatted') read(111,*) do M=1,L_NSPECTI-1 read(111,*) BWNI(M) end do read(111,*) lastband close(111) BWNI(L_NSPECTI) =lastband(1) BWNI(L_NSPECTI+1)=lastband(2) print*,'' print*,'setspi: IR band limits:' do M=1,L_NSPECTI+1 print*,m,'-->',BWNI(M),' cm^-1' end do ! Set up mean wavenumbers and wavenumber deltas. Units of ! wavenumbers is cm^(-1); units of wavelengths is microns. do M=1,L_NSPECTI WNOI(M) = 0.5*(BWNI(M+1)+BWNI(M)) DWNI(M) = BWNI(M+1)-BWNI(M) WAVEI(M) = 1.0E+4/WNOI(M) BLAMI(M) = 0.01/BWNI(M) end do BLAMI(M) = 0.01/BWNI(M) ! note M=L_NSPECTI+1 after loop due to Fortran bizarreness !======================================================================= ! For each IR wavelength interval, compute the integral of B(T), the ! Planck function, divided by the wavelength interval, in cm-1. The ! integration is in MKS units, the final answer is the same as the ! original planck.f; W m^-2 wavenumber^-1, where wavenumber is in CM^-1. print*,'' print*,'setspi: Current Planck integration range:' print*,'T = ',dble(NTstar)/1.0D+1, ' to ',dble(NTstop)/1.0D+1,' K.' do NW=1,L_NSPECTI a = 1.0D-2/BWNI(NW+1) b = 1.0D-2/BWNI(NW) bpa = (b+a)/2.0 bma = (b-a)/2.0 do nt=NTstar,NTstop T = dble(NT)/1.0D+1 ans = 0.0D0 do mm=1,12 y = bma*x(mm)+bpa ans = ans + w(mm)*c1/(y**5*(exp(c2/(y*T))-1.0D0)) end do planckir(NW,nt-NTstar+1) = ans*bma/(PI*DWNI(NW)) end do end do ! force planck=sigma*eps*T^4 for each temperature in array if(forceEC)then print*,'setspi: Force F=sigma*eps*T^4 for all values of T!' do nt=NTstar,NTstop plancksum=0.0 T=dble(NT)/1.0D+1 do NW=1,L_NSPECTI plancksum=plancksum+ & planckir(NW,nt-NTstar+1)*DWNI(NW)*pi end do do NW=1,L_NSPECTI planckir(NW,nt-NTstar+1)= & planckir(NW,nt-NTstar+1)* & sigma*(dble(nt)/1.0D+1)**4/plancksum end do end do endif if(planckcheck)then ! check energy conservation at lower temperature boundary plancksum=0.0 nt=NTstar do NW=1,L_NSPECTI plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi end do print*,'setspi: At lower limit:' print*,'in model sig*T^4 = ',plancksum,' W m^-2' print*,'actual sig*T^4 = ',sigma*(dble(nt)/1.0D+1)**4,' W m^-2' ! check energy conservation at upper temperature boundary plancksum=0.0 nt=NTstop do NW=1,L_NSPECTI plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi end do print*,'setspi: At upper limit:' print*,'in model sig*T^4 = ',plancksum,' W m^-2' print*,'actual sig*T^4 = ',sigma*(dble(nt)/1.0D+1)**4,' W m^-2' print*,'' endif return end subroutine setspi