source: trunk/LMDZ.PLUTO.old/libf/phypluto/setspv.F90

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