source: trunk/LMDZ.GENERIC/libf/phystd/setspi.F90 @ 537

Last change on this file since 537 was 374, checked in by emillour, 13 years ago

Generic GCM: Upgrade: The location of the 'datagcm' directory can now be given in the callphys.def file ( datadir = /absolute/path/to/datagcm ). Changed "datafile.h" into a F90 module "datafile_mod.F90" and spread this change to all routines that used to use "datafile.h".
EM
File size: 6.8 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!=======================================================================
70!     Set up spectral bands - wavenumber [cm^(-1)]. Go from smaller to
71!     larger wavenumbers.
72
73      write(temp1,'(i2.2)') L_NSPECTI
74      !file_id='/corrk_data/' // corrkdir(1:LEN_TRIM(corrkdir)) // '/narrowbands_IR.in'
75      file_id='/corrk_data/'//trim(adjustl(banddir))//'/narrowbands_IR.in'
76      file_path=TRIM(datadir)//TRIM(file_id)
77
78      ! check that the file exists
79      inquire(FILE=file_path,EXIST=file_ok)
80      if(.not.file_ok) then
81         write(*,*)'The file ',TRIM(file_path)
82         write(*,*)'was not found by setspi.F90, exiting.'
83         write(*,*)'Check that your path to datagcm:',trim(datadir)
84         write(*,*)' is correct. You can change it in callphys.def with:'
85         write(*,*)' datadir = /absolute/path/to/datagcm'
86         write(*,*)'Also check that the corrkdir you chose in callphys.def exists.'
87         call abort
88      endif
89     
90      ! check that the file contains the right number of bands
91      call system('wc -l '//file_path//' > bandlen.txt')
92      open(131,file='bandlen.txt', form='formatted')
93      read(131,*) file_entries
94      close(131)
95      call system('rm -f bandlen.txt')
96
97      if(file_entries-1.ne.L_NSPECTI) then
98         write(*,*) 'setspi: L_NSPECTI = ',L_NSPECTI, 'in the model, but there are '
99         write(*,*) file_entries-1,'entries in ', &
100              TRIM(file_path),', exiting.'
101         call abort
102      endif
103
104      ! load and display the data
105      open(111,file=file_path,form='formatted')
106      read(111,*)
107      do M=1,L_NSPECTI-1
108         read(111,*) BWNI(M)
109      end do
110      read(111,*) lastband
111      close(111)
112      BWNI(L_NSPECTI)  =lastband(1)
113      BWNI(L_NSPECTI+1)=lastband(2)
114
115      print*,''
116      print*,'setspi: IR band limits:'
117      do M=1,L_NSPECTI+1
118         print*,m,'-->',BWNI(M),' cm^-1'
119      end do
120
121!     Set up mean wavenumbers and wavenumber deltas.  Units of
122!     wavenumbers is cm^(-1); units of wavelengths is microns.
123
124      do M=1,L_NSPECTI
125         WNOI(M)  = 0.5*(BWNI(M+1)+BWNI(M))
126         DWNI(M)  = BWNI(M+1)-BWNI(M)
127         WAVEI(M) = 1.0E+4/WNOI(M)
128         BLAMI(M) = 0.01/BWNI(M)         
129      end do
130      BLAMI(M) = 0.01/BWNI(M)
131!     note M=L_NSPECTI+1 after loop due to Fortran bizarreness
132
133!=======================================================================
134!     For each IR wavelength interval, compute the integral of B(T), the
135!     Planck function, divided by the wavelength interval, in cm-1.  The
136!     integration is in MKS units, the final answer is the same as the
137!     original planck.f; W m^-2 wavenumber^-1, where wavenumber is in CM^-1.
138
139      print*,''
140      print*,'setspi: Current Planck integration range:'
141      print*,'T = ',dble(NTstar)/1.0D+1, ' to ',dble(NTstop)/1.0D+1,' K.'
142
143      do NW=1,L_NSPECTI
144         a = 1.0D-2/BWNI(NW+1)
145         b = 1.0D-2/BWNI(NW)
146         bpa = (b+a)/2.0
147         bma = (b-a)/2.0
148         do nt=NTstar,NTstop
149            T   = dble(NT)/1.0D+1
150            ans = 0.0D0
151
152            do mm=1,12
153               y    = bma*x(mm)+bpa
154               ans  = ans + w(mm)*c1/(y**5*(exp(c2/(y*T))-1.0D0))
155            end do
156
157            planckir(NW,nt-NTstar+1) = ans*bma/(PI*DWNI(NW))
158         end do
159      end do
160         
161      ! force planck=sigma*eps*T^4 for each temperature in array
162      if(forceEC)then
163         print*,'setspi: Force F=sigma*eps*T^4 for all values of T!'
164         do nt=NTstar,NTstop
165            plancksum=0.0
166            T=dble(NT)/1.0D+1
167       
168            do NW=1,L_NSPECTI
169               plancksum=plancksum+  &
170                  planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
171            end do
172
173            do NW=1,L_NSPECTI
174               planckir(NW,nt-NTstar+1)=     &
175                  planckir(NW,nt-NTstar+1)*  &
176                          sigma*(dble(nt)/1.0D+1)**4/plancksum
177            end do
178         end do
179      endif
180
181      if(planckcheck)then
182         ! check energy conservation at lower temperature boundary
183         plancksum=0.0
184         nt=NTstar
185         do NW=1,L_NSPECTI
186            plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
187         end do
188         print*,'setspi: At lower limit:'
189         print*,'in model sig*T^4 = ',plancksum,' W m^-2'
190         print*,'actual sig*T^4   = ',sigma*(dble(nt)/1.0D+1)**4,' W m^-2'
191         
192         ! check energy conservation at upper temperature boundary
193         plancksum=0.0
194         nt=NTstop
195         do NW=1,L_NSPECTI
196            plancksum=plancksum+planckir(NW,nt-NTstar+1)*DWNI(NW)*pi
197         end do
198         print*,'setspi: At upper limit:'
199         print*,'in model sig*T^4 = ',plancksum,' W m^-2'
200         print*,'actual sig*T^4   = ',sigma*(dble(nt)/1.0D+1)**4,' W m^-2'
201         print*,''
202      endif
203
204      return
205    end subroutine setspi
Note: See TracBrowser for help on using the repository browser.