source: trunk/LMDZ.GENERIC/libf/phystd/setspv.F90 @ 253

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

Generic GCM

  • Massive update to version 0.7

EM+RW

File size: 4.0 KB
Line 
1      subroutine setspv
2
3!==================================================================
4!     
5!     Purpose
6!     -------
7!     Set up spectral intervals, stellar spectrum and Rayleigh
8!     opacity in the shortwave.
9!     
10!     Authors
11!     -------
12!     Adapted from setspv in the NASA Ames radiative code by
13!     Robin Wordsworth (2009).
14!
15!     Called by
16!     ---------
17!     callcorrk.F
18!     
19!     Calls
20!     -----
21!     ave_stelspec.F
22!     
23!==================================================================
24
25      use radinc_h,    only: L_NSPECTV, corrkdir, banddir
26      use radcommon_h, only: BWNV,BLAMV,WNOV,DWNV,WAVEV, &
27                             STELLARF,TAURAY,gastype
28
29      implicit none
30
31#include "comcstfi.h"
32#include "callkeys.h"
33#include "datafile.h"
34
35      logical file_ok
36
37      integer N, M, file_entries
38
39      character(len=30)  :: temp1
40      character(len=100) :: file_id
41      character(len=100) :: file_path
42
43      real*8 :: lastband(2)
44
45      real*8 STELLAR(L_NSPECTV)
46      real*8 sum
47
48!=======================================================================
49!     Set up spectral bands - wavenumber [cm^(-1)]. Go from smaller to
50!     larger wavenumbers, the same as in the IR.
51
52      write(temp1,'(i2.2)') L_NSPECTV
53      file_id='/corrk_data/'//trim(adjustl(banddir))//'/narrowbands_VI.in'
54      file_path=datafile(1:LEN_TRIM(datafile))//file_id(1:LEN_TRIM(file_id))
55
56      ! check that the file exists
57      inquire(FILE=file_path,EXIST=file_ok)
58      if(.not.file_ok) then
59         write(*,*)'The file ',file_path(1:LEN_TRIM(file_path))
60         write(*,*)'was not found by setspv.F90, exiting.'
61         call abort
62      endif
63     
64      ! check that the file contains the right number of bands
65      call system('wc -l '//file_path//' > bandlen.txt')
66      open(131,file='bandlen.txt', form='formatted')
67      read(131,*) file_entries
68      close(131)
69      call system('rm -f bandlen.txt')
70
71      if(file_entries-1.ne.L_NSPECTV) then
72         write(*,*) 'L_NSPECTV = ',L_NSPECTV, 'in the model, but there are '
73         write(*,*) file_entries-1,'entries in ', &
74              file_path(1:LEN_TRIM(file_path)),', exiting.'
75         call abort
76      endif
77
78      ! load and display the data
79      open(111,file=file_path,form='formatted')
80      read(111,*)
81       do M=1,L_NSPECTV-1
82         read(111,*) BWNV(M)
83      end do
84      read(111,*) lastband
85      close(111)
86      BWNV(L_NSPECTV)  =lastband(1)
87      BWNV(L_NSPECTV+1)=lastband(2)
88
89
90      print*,'VI band limits:'
91      do M=1,L_NSPECTV+1
92         print*,m,'-->',BWNV(M),' cm^-1'
93      end do
94      print*,' '
95
96!     Set up mean wavenumbers and wavenumber deltas.  Units of
97!     wavenumbers is cm^(-1); units of wavelengths is microns.
98
99      do M=1,L_NSPECTV
100         WNOV(M)  = 0.5*(BWNV(M+1)+BWNV(M))
101         DWNV(M)  = BWNV(M+1)-BWNV(M)
102         WAVEV(M) = 1.0E+4/WNOV(M)
103         BLAMV(M) = 0.01/BWNV(M)
104      end do
105      BLAMV(M) = 0.01/BWNV(M) ! wavelength in METERS for aerosol stuff
106!     note M=L_NSPECTV+1 after loop due to Fortran bizarreness
107
108!=======================================================================
109!     Set up stellar spectrum
110
111      write(*,*)'Interpolating stellar spectrum from the hires data...'
112      call ave_stelspec(STELLAR)
113
114!     Sum the stellar flux, and write out the result. 
115      sum = 0.0 
116      do N=1,L_NSPECTV
117         STELLARF(N) = STELLAR(N) * Fat1AU
118         sum         = sum+STELLARF(N)
119      end do
120      write(6,'(" Stellar flux at 1 AU = ",f7.2," W m-2")') sum
121      print*,' '
122
123
124!=======================================================================
125!     Set up the wavelength independent part of the Rayleigh scattering.
126!     The pressure dependent part will be computed elsewhere (OPTCV).
127!     WAVEV is in microns.  There is no Rayleigh scattering in the IR.
128
129      if(rayleigh) then
130         call calc_rayleigh
131      else
132         print*,'No Rayleigh scattering, check for NaN in output!'
133         do N=1,L_NSPECTV
134            TAURAY(N) = 1E-16
135         end do
136      endif
137
138      RETURN
139    END subroutine setspv
Note: See TracBrowser for help on using the repository browser.