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

Last change on this file since 1384 was 1384, checked in by emillour, 10 years ago

Generic GCM:

  • Some code cleanup: turning comcstfi.h into module comcstfi_mod.F90

EM

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