source: trunk/LMDZ.PLUTO.old/libf/phypluto/setspi.F90 @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM
File size: 6.6 KB
Line 
1      subroutine setspi
2
3!==================================================================
4!     
5!     Purpose
6!     -------
7!     Set up spectral intervals and Planck function in the longwave.
8!     
9!     Authors
10!     -------
11!     Adapted from setspi in the NASA Ames radiative code by
12!     Robin Wordsworth (2009).
13!     
14!     Called by
15!     ---------
16!     callcorrk.F
17!     
18!     Calls
19!     -----
20!     none
21!     
22!==================================================================
23
24      use radinc_h,    only: L_NSPECTI,corrkdir,banddir,NTstar,NTstop
25      use radcommon_h, only: BWNI,BLAMI,WNOI,DWNI,WAVEI,planckir,sigma
26      use datafile_mod, only: datadir
27
28      implicit none
29
30#include "callkeys.h"
31#include "comcstfi.h"
32
33      logical file_ok
34      integer nw, nt, m, mm, file_entries
35      real*8 a, b, ans, y, bpa, bma, T
36
37      character(len=30)  :: temp1
38      character(len=100) :: file_id
39      character(len=100) :: file_path
40
41!     C1 and C2 values from Goody and Yung (2nd edition)  MKS units
42!     These values lead to a "sigma" (sigma*T^4) of 5.67032E-8 W m^-2 K^-4
43
44      real*8 :: c1 = 3.741832D-16 ! W m^-2
45      real*8 :: c2 = 1.438786D-2  ! m K
46     
47      real*8 :: lastband(2), plancksum
48
49      logical forceEC, planckcheck
50
51      real*8 :: x(12) = [ -0.981560634246719D0,  -0.904117256370475D0, &
52      -0.769902674194305D0,  -0.587317954286617D0,                     &
53      -0.367831498998180D0,  -0.125233408511469D0,                     &
54       0.125233408511469D0,   0.367831498998180D0,                     &
55       0.587317954286617D0,   0.769902674194305D0,                     &
56       0.904117256370475D0,   0.981560634246719D0  ]
57
58      real*8 :: w(12) = [  0.047175336386512D0,   0.106939325995318D0, &
59           0.160078328543346D0,   0.203167426723066D0,                 &
60           0.233492536538355D0,   0.249147045813403D0,                 &
61           0.249147045813403D0,   0.233492536538355D0,                 &
62           0.203167426723066D0,   0.160078328543346D0,                 &
63           0.106939325995318D0,   0.047175336386512D0  ]
64      mm=0
65
66      forceEC=.false.
67      planckcheck=.true.
68
69      planckir(1:L_NSPECTI,1:NTstop-NTstar+1) = 0.    !TB16 bug if not set when changing layers
70!=======================================================================
71!     Set up spectral bands - wavenumber [cm^(-1)]. Go from smaller to
72!     larger wavenumbers.
73
74      write(temp1,'(i2.2)') L_NSPECTI
75      !file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/narrowbands_IR.in'
76      file_id='/corrk_data/'//trim(adjustl(banddir))//'/narrowbands_IR.in'
77      file_path=TRIM(datadir)//TRIM(file_id)
78     
79      ! check that the file exists
80      inquire(FILE=file_path,EXIST=file_ok)
81      if(.not.file_ok) then
82         write(*,*)'The file ',file_path(1:LEN_TRIM(file_path))
83         write(*,*)'was not found by setspi.F90, exiting.'
84         call abort
85      endif
86     
87      ! check that the file contains the right number of bands
88      call system('wc -l '//file_path//' > bandlen.txt')
89      open(131,file='bandlen.txt', form='formatted')
90      read(131,*) file_entries
91      close(131)
92      call system('rm -f bandlen.txt')
93
94      if(file_entries-1.ne.L_NSPECTI) then
95         write(*,*) 'L_NSPECTI = ',L_NSPECTI, 'in the model, but there are '
96         write(*,*) file_entries-1,'entries in ', &
97              file_path(1:LEN_TRIM(file_path)),', exiting.'
98         call abort
99      endif
100
101      ! load and display the data
102      open(111,file=file_path,form='formatted')
103      read(111,*)
104      do M=1,L_NSPECTI-1
105         read(111,*) BWNI(M)
106      end do
107      read(111,*) lastband
108      close(111)
109      BWNI(L_NSPECTI)  =lastband(1)
110      BWNI(L_NSPECTI+1)=lastband(2)
111
112      print*,''
113      print*,'IR band limits:'
114      do M=1,L_NSPECTI+1
115         print*,m,'-->',BWNI(M),' cm^-1'
116      end do
117
118!     Set up mean wavenumbers and wavenumber deltas.  Units of
119!     wavenumbers is cm^(-1); units of wavelengths is microns.
120
121      do M=1,L_NSPECTI
122         WNOI(M)  = 0.5*(BWNI(M+1)+BWNI(M))
123         DWNI(M)  = BWNI(M+1)-BWNI(M)
124         WAVEI(M) = 1.0E+4/WNOI(M)
125         BLAMI(M) = 0.01/BWNI(M)         
126      end do
127      BLAMI(M) = 0.01/BWNI(M)
128!     note M=L_NSPECTI+1 after loop due to Fortran bizarreness
129
130!=======================================================================
131!     For each IR wavelength interval, compute the integral of B(T), the
132!     Planck function, divided by the wavelength interval, in cm-1.  The
133!     integration is in MKS units, the final answer is the same as the
134!     original planck.f; W m^-2 wavenumber^-1, where wavenumber is in CM^-1.
135
136      print*,''
137      print*,'Current Planck integration range:'
138      print*,'T = ',dble(NTstar)/1.0D+1, ' to ',dble(NTstop)/1.0D+1,' K.'
139
140      do NW=1,L_NSPECTI
141         a = 1.0D-2/BWNI(NW+1)
142         b = 1.0D-2/BWNI(NW)
143         bpa = (b+a)/2.0
144         bma = (b-a)/2.0
145         do nt=NTstar,NTstop
146            T   = dble(NT)/1.0D+1
147            ans = 0.0D0
148
149            do mm=1,12
150               y    = bma*x(mm)+bpa
151               ans  = ans + w(mm)*c1/(y**5*(exp(c2/(y*T))-1.0D0))
152            end do
153
154            planckir(NW,nt-NTstar+1) = ans*bma/(PI*DWNI(NW))
155         end do
156      end do
157         
158      ! force planck=sigma*eps*T^4 for each temperature in array
159      if(forceEC)then
160         do nt=NTstar,NTstop
161            plancksum=0.0
162            T=dble(NT)/1.0D+1
163       
164            do NW=1,L_NSPECTI
165               plancksum=plancksum+  &
166                  planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
167            end do
168
169            do NW=1,L_NSPECTI
170               planckir(NW,nt-NTstar+1)=     &
171                  planckir(NW,nt-NTstar+1)*  &
172                          sigma*(dble(nt)/1.0D+1)**4/plancksum
173            end do
174         end do
175      endif
176
177      if(planckcheck)then
178         ! check energy conservation at lower temperature boundary
179         plancksum=0.0
180         nt=NTstar
181         do NW=1,L_NSPECTI
182            plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
183         end do
184         print*,'At lower limit:'
185         print*,'in model sig*T^4 = ',plancksum,' W m^-2'
186         print*,'actual sig*T^4   = ',sigma*(dble(nt)/1.0D+1)**4,' W m^-2'
187         
188         ! check energy conservation at upper temperature boundary
189         plancksum=0.0
190         nt=NTstop
191         do NW=1,L_NSPECTI
192            plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
193         end do
194         print*,'At upper limit:'
195         print*,'in model sig*T^4 = ',plancksum,' W m^-2'
196         print*,'actual sig*T^4   = ',sigma*(dble(nt)/1.0D+1)**4,' W m^-2'
197         print*,''
198      endif
199
200      return
201    end subroutine setspi
Note: See TracBrowser for help on using the repository browser.