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

Last change on this file since 1974 was 1397, checked in by milmd, 10 years ago

In LMDZ.GENERIC replacement of all phystd .h files by module files.

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