| 1 | program generate_kmatrix |
|---|
| 2 | implicit none |
|---|
| 3 | |
|---|
| 4 | ! ------------------------------------------------------------- |
|---|
| 5 | ! Purpose: to read hires spectra and produce kdistribution data |
|---|
| 6 | ! Authour: Robin Wordsworth (2010) |
|---|
| 7 | ! ------------------------------------------------------------- |
|---|
| 8 | |
|---|
| 9 | integer nmolec,mol |
|---|
| 10 | integer i,j,k,im |
|---|
| 11 | integer Nb |
|---|
| 12 | |
|---|
| 13 | character*4 fname |
|---|
| 14 | character*32 fnum |
|---|
| 15 | |
|---|
| 16 | integer iT, iP, iQ, Nlevs |
|---|
| 17 | |
|---|
| 18 | double precision nuT, sigT, kT |
|---|
| 19 | |
|---|
| 20 | integer Nq |
|---|
| 21 | parameter (Nq=16) ! we choose this here |
|---|
| 22 | |
|---|
| 23 | integer Nmax |
|---|
| 24 | parameter (Nmax=1000) ! maximum number of levels |
|---|
| 25 | |
|---|
| 26 | integer Nmol_max |
|---|
| 27 | parameter (Nmol_max=10) ! maximum number of radiatively active species |
|---|
| 28 | |
|---|
| 29 | integer Nkmax |
|---|
| 30 | ! large array = slow calculation... this makes a huge difference |
|---|
| 31 | !parameter (Nkmax=50000) ! works for IR, 27 bands |
|---|
| 32 | !parameter (Nkmax=100000) ! works for IR, 27 bands |
|---|
| 33 | !parameter (Nkmax=140000) ! works for IR, 32 bands |
|---|
| 34 | !parameter (Nkmax=250000) ! works for VI, 36 bands |
|---|
| 35 | parameter (Nkmax=400000) ! works for IR, 8 bands |
|---|
| 36 | !parameter (Nkmax=400000) ! works for VI, 12 bands |
|---|
| 37 | !parameter (Nkmax=200000) ! works for IR, 30 bands |
|---|
| 38 | |
|---|
| 39 | integer NbMAX |
|---|
| 40 | parameter (NbMAX=100) |
|---|
| 41 | |
|---|
| 42 | integer Ngres |
|---|
| 43 | parameter (Ngres=200) |
|---|
| 44 | |
|---|
| 45 | character*10 molec_names(1:Nmol_max) |
|---|
| 46 | double precision T(1:Nmax), P(1:Nmax), Q(1:Nmax) |
|---|
| 47 | |
|---|
| 48 | |
|---|
| 49 | double precision pi |
|---|
| 50 | parameter (pi=3.14159265) |
|---|
| 51 | double precision x |
|---|
| 52 | double precision y(1:Nq+2) |
|---|
| 53 | |
|---|
| 54 | double precision nu_min(1:NbMAX) ! make allocatable? |
|---|
| 55 | double precision nu_max(1:NbMAX) |
|---|
| 56 | |
|---|
| 57 | integer quad, band |
|---|
| 58 | |
|---|
| 59 | double precision ktemp(1:Nq) |
|---|
| 60 | double precision w(1:Nq), gw(1:Nq) |
|---|
| 61 | |
|---|
| 62 | double precision, allocatable :: kdist_temp(:) |
|---|
| 63 | double precision, allocatable :: kdist(:,:,:,:,:) |
|---|
| 64 | |
|---|
| 65 | double precision nuband(1:Nkmax), kband(1:Nkmax), nuwband(1:Nkmax) ! better static I think |
|---|
| 66 | double precision nuwbtot |
|---|
| 67 | |
|---|
| 68 | double precision gfine(1:Ngres) |
|---|
| 69 | double precision klevs(1:Ngres) |
|---|
| 70 | double precision addline(1:Nkmax) |
|---|
| 71 | |
|---|
| 72 | double precision kmax, kmin, lkmax, lkmin, dlogK, klev |
|---|
| 73 | integer ic, ib, ig, ik1, ik2 |
|---|
| 74 | |
|---|
| 75 | integer filerr |
|---|
| 76 | logical filexi |
|---|
| 77 | |
|---|
| 78 | logical addCIA |
|---|
| 79 | integer iCO2 |
|---|
| 80 | double precision k1,k2,kCIA,amg,Pp |
|---|
| 81 | |
|---|
| 82 | real truemean, distmean, etot |
|---|
| 83 | |
|---|
| 84 | |
|---|
| 85 | addCIA=.false. |
|---|
| 86 | |
|---|
| 87 | print*,'Creating the GCM correlated-k data...' |
|---|
| 88 | print*,' ' |
|---|
| 89 | |
|---|
| 90 | !-------------------------------------------------------- |
|---|
| 91 | ! load temperature, pressure, mixing ratio data |
|---|
| 92 | open(9,file='T.dat') |
|---|
| 93 | read(9,*) iT |
|---|
| 94 | do i=1,iT |
|---|
| 95 | read(9,*) T(i) |
|---|
| 96 | enddo |
|---|
| 97 | close(9) |
|---|
| 98 | |
|---|
| 99 | open(10,file='p.dat') |
|---|
| 100 | read(10,*) iP |
|---|
| 101 | do i=1,iP |
|---|
| 102 | read(10,*) P(i) |
|---|
| 103 | P(i)=10.**P(i) ! note we keep P in millibars for the conversion later |
|---|
| 104 | enddo |
|---|
| 105 | close(10) |
|---|
| 106 | |
|---|
| 107 | open(11,file='Q.dat') |
|---|
| 108 | read(11,*) nmolec |
|---|
| 109 | do mol=1,nmolec |
|---|
| 110 | read(11,*) molec_names(mol) |
|---|
| 111 | enddo |
|---|
| 112 | read(11,*) iQ |
|---|
| 113 | do i=1,iQ |
|---|
| 114 | read(11,*) Q(i) |
|---|
| 115 | enddo |
|---|
| 116 | close(11) |
|---|
| 117 | |
|---|
| 118 | Nlevs=iT*iP*iQ ! total number of layers |
|---|
| 119 | |
|---|
| 120 | print*,'Temperature layers: ',iT |
|---|
| 121 | print*,'Pressure layers: ',iP |
|---|
| 122 | print*,'Mixing ratio layers: ',iQ |
|---|
| 123 | print*,'Total: ',Nlevs |
|---|
| 124 | |
|---|
| 125 | |
|---|
| 126 | if(addCIA)then |
|---|
| 127 | |
|---|
| 128 | if(nmolec.gt.2)then |
|---|
| 129 | print*,'I cant deal with this situation yet' |
|---|
| 130 | call abort |
|---|
| 131 | endif |
|---|
| 132 | |
|---|
| 133 | iCO2=0 |
|---|
| 134 | do mol=1,nmolec |
|---|
| 135 | if(molec_names(mol).eq.'CO2')then |
|---|
| 136 | iCO2=mol |
|---|
| 137 | endif |
|---|
| 138 | enddo |
|---|
| 139 | |
|---|
| 140 | !if(iCO2.eq.1)then |
|---|
| 141 | print*,' ' |
|---|
| 142 | print*,'--- Adding CO2 CIA ---' |
|---|
| 143 | !else |
|---|
| 144 | ! print*,'Warning, not CO2+H2O mixture, need CO2 mixing ratio' |
|---|
| 145 | ! call abort |
|---|
| 146 | !endif |
|---|
| 147 | endif |
|---|
| 148 | |
|---|
| 149 | !-------------------------------------------------------- |
|---|
| 150 | ! load spectral data |
|---|
| 151 | open(12,file='narrowbands.in') |
|---|
| 152 | read(12,*) Nb |
|---|
| 153 | |
|---|
| 154 | do band=1,Nb |
|---|
| 155 | read(12,*) nu_min(band),nu_max(band) |
|---|
| 156 | enddo |
|---|
| 157 | close(12) |
|---|
| 158 | |
|---|
| 159 | print*,' ' |
|---|
| 160 | print*,'Number of spectral bands: ',Nb |
|---|
| 161 | |
|---|
| 162 | print*,'Band limits:' |
|---|
| 163 | do band=1,Nb |
|---|
| 164 | print*,band,'-->',nu_min(band),nu_max(band),' cm^-1' |
|---|
| 165 | end do |
|---|
| 166 | |
|---|
| 167 | !-------------------------------------------------------- |
|---|
| 168 | ! create g-space data |
|---|
| 169 | do quad=0,Nq+1 |
|---|
| 170 | x=dble(quad)/dble(Nq+2) |
|---|
| 171 | !y(quad+1)=sin(pi*x)*exp(-5*x) |
|---|
| 172 | !y(quad+1)=sin(pi*x)*exp(-9*x) ! CO2_H2Ovar 38x36 |
|---|
| 173 | y(quad+1)=sin(pi*x)*exp(-12*x) |
|---|
| 174 | !y(quad+1)=sin(pi*x)*exp(-15*x) ! CO2_CH4 |
|---|
| 175 | enddo |
|---|
| 176 | |
|---|
| 177 | do quad=1,Nq |
|---|
| 178 | w(quad)=y(quad+1) |
|---|
| 179 | enddo |
|---|
| 180 | w=w/sum(w) |
|---|
| 181 | |
|---|
| 182 | print*,' ' |
|---|
| 183 | print*,'Number of g-space intervals: ',Nq |
|---|
| 184 | |
|---|
| 185 | print*,'g-space weighting:' |
|---|
| 186 | do quad=1,Nq |
|---|
| 187 | print*,quad,'-->',w(quad) |
|---|
| 188 | end do |
|---|
| 189 | |
|---|
| 190 | ! write g-space data |
|---|
| 191 | open(14,file='g.dat') |
|---|
| 192 | write(14,*) Nq+1 |
|---|
| 193 | |
|---|
| 194 | do quad=1,Nq |
|---|
| 195 | write(14,*) w(quad) |
|---|
| 196 | enddo |
|---|
| 197 | write(14,*) 0.0 |
|---|
| 198 | close(14) |
|---|
| 199 | |
|---|
| 200 | ! cumulative sum for g-vector |
|---|
| 201 | gw(1)=w(1) |
|---|
| 202 | do quad=2,Nq |
|---|
| 203 | gw(quad)=gw(quad-1)+w(quad) |
|---|
| 204 | enddo |
|---|
| 205 | |
|---|
| 206 | ! allocate 5D matrices |
|---|
| 207 | allocate(kdist_temp(Nq)) |
|---|
| 208 | allocate(kdist(iT,iP,iQ,Nb,Nq+1)) |
|---|
| 209 | |
|---|
| 210 | |
|---|
| 211 | nuT=0.0 |
|---|
| 212 | do i=1,iQ |
|---|
| 213 | do j=1,iP |
|---|
| 214 | |
|---|
| 215 | |
|---|
| 216 | if(addCIA.and.(iCO2.eq.1))then |
|---|
| 217 | Pp=P(j)*(1.0-Q(i)) |
|---|
| 218 | else |
|---|
| 219 | Pp=P(j)*Q(i) |
|---|
| 220 | endif |
|---|
| 221 | !if(addCIA) Pp=P(j)*(1.0-Q(i)) |
|---|
| 222 | |
|---|
| 223 | |
|---|
| 224 | do k=1,iT |
|---|
| 225 | im = (i-1)*iP*iT + (j-1)*iT + k |
|---|
| 226 | |
|---|
| 227 | !----------------------------------------------- |
|---|
| 228 | ! open kspectrum file |
|---|
| 229 | write(fnum,*) im |
|---|
| 230 | if(im<10)then |
|---|
| 231 | fname='k00'//trim(adjustl(fnum)) |
|---|
| 232 | elseif(im<100)then |
|---|
| 233 | fname='k0'//trim(adjustl(fnum)) |
|---|
| 234 | else |
|---|
| 235 | fname='k'//trim(adjustl(fnum)) |
|---|
| 236 | endif |
|---|
| 237 | |
|---|
| 238 | ! check that the file exists |
|---|
| 239 | inquire(FILE='hires_spectrum/'//fname,EXIST=filexi) |
|---|
| 240 | if(.not.filexi) then |
|---|
| 241 | write(*,*)'The file ',fname(1:LEN_TRIM(fname)) |
|---|
| 242 | write(*,*)'was not found in the folder hires_spectrum, exiting.' |
|---|
| 243 | call abort |
|---|
| 244 | endif |
|---|
| 245 | |
|---|
| 246 | print*,'Opening file ',fname(1:LEN_TRIM(fname)) |
|---|
| 247 | open(17,file='hires_spectrum/'//fname) |
|---|
| 248 | |
|---|
| 249 | nuT=0.0 |
|---|
| 250 | do band=1,Nb |
|---|
| 251 | print*,'Band ---> ',band |
|---|
| 252 | |
|---|
| 253 | !-------------------------------------------- |
|---|
| 254 | ! write kdata to temporary band arrays |
|---|
| 255 | ! and find max / min k-values in band |
|---|
| 256 | do while (nuT.lt.nu_min(band)) |
|---|
| 257 | read(17,*,IOSTAT=filerr) nuT, sigT, kT |
|---|
| 258 | ! ideally need to catch all errors here (ask Ehouarn) |
|---|
| 259 | if(filerr.lt.0)then |
|---|
| 260 | print*, & |
|---|
| 261 | 'nu_min greater than highest value in file' |
|---|
| 262 | print*,'IOSTAT=',filerr |
|---|
| 263 | print*,'nu_min=',nu_min(band) |
|---|
| 264 | print*,'nu_endoffile=',nuT |
|---|
| 265 | call abort |
|---|
| 266 | endif |
|---|
| 267 | enddo |
|---|
| 268 | |
|---|
| 269 | ic=1 |
|---|
| 270 | kmin=kT |
|---|
| 271 | kmax=kT |
|---|
| 272 | nuband(:)=0 |
|---|
| 273 | kband(:)=0 |
|---|
| 274 | |
|---|
| 275 | do while (nuT.lt.nu_max(band)) |
|---|
| 276 | |
|---|
| 277 | nuband(ic)=nuT |
|---|
| 278 | kband(ic)=kT |
|---|
| 279 | |
|---|
| 280 | ! insert CIA option here |
|---|
| 281 | if(addCIA)then |
|---|
| 282 | |
|---|
| 283 | k1=0.0 |
|---|
| 284 | k2=0.0 |
|---|
| 285 | if(nuT.lt.500.0) call getspc(T(k),nuT,k1) |
|---|
| 286 | if((nuT.gt.1000.0).and.(nuT.lt.1800.0)) call baranov(T(k),nuT,k2) |
|---|
| 287 | ! cm^-1.amagat^2 |
|---|
| 288 | |
|---|
| 289 | amg=273.15D+0/T(k)*(Pp/1013.25) ! amagats of CO2 |
|---|
| 290 | |
|---|
| 291 | kCIA=(k1+k2)*1.0D+2*amg**2.0D+0 ! m^-1 |
|---|
| 292 | kband(ic)=kband(ic)+kCIA |
|---|
| 293 | kT=kband(ic) |
|---|
| 294 | |
|---|
| 295 | !kband(ic)=kCIA |
|---|
| 296 | !kT=kCIA |
|---|
| 297 | !if(nuband(ic).gt.1250.0)then |
|---|
| 298 | ! print*,'P=',P(i) |
|---|
| 299 | ! print*,'T=',T(k) |
|---|
| 300 | ! print*,'amg=',amg |
|---|
| 301 | ! print*,'nu=',nuband(ic) |
|---|
| 302 | ! print*,'k1=',k1 |
|---|
| 303 | ! print*,'k2=',k2 |
|---|
| 304 | ! print*,'kCIA=',kCIA |
|---|
| 305 | ! call abort |
|---|
| 306 | !endif |
|---|
| 307 | |
|---|
| 308 | endif |
|---|
| 309 | |
|---|
| 310 | if(kT.lt.kmin)then |
|---|
| 311 | kmin=kT |
|---|
| 312 | endif |
|---|
| 313 | if(kT.gt.kmax)then ! bug corrected! |
|---|
| 314 | kmax=kT |
|---|
| 315 | endif |
|---|
| 316 | |
|---|
| 317 | read(17,*,IOSTAT=filerr) nuT, sigT, kT |
|---|
| 318 | if(filerr.lt.0)then |
|---|
| 319 | print*, & |
|---|
| 320 | 'nu_max greater than highest value in file' |
|---|
| 321 | print*,'IOSTAT=',filerr |
|---|
| 322 | print*,'nu_max=',nu_max(band) |
|---|
| 323 | print*,'nu_endoffile=',nuT |
|---|
| 324 | call abort |
|---|
| 325 | endif |
|---|
| 326 | |
|---|
| 327 | |
|---|
| 328 | ic=ic+1 |
|---|
| 329 | |
|---|
| 330 | if(ic.gt.Nkmax)then |
|---|
| 331 | print*,'Spectral resolution too low for these bands, exiting!' |
|---|
| 332 | print*,'Increase Nkmax in generate_kmatrix.F90.' |
|---|
| 333 | call abort |
|---|
| 334 | endif |
|---|
| 335 | |
|---|
| 336 | enddo |
|---|
| 337 | nuband(ic)=nuT |
|---|
| 338 | kband(ic)=kT |
|---|
| 339 | |
|---|
| 340 | if(kmax.le.0.0)then |
|---|
| 341 | !if(kmax.le.1.0e-15)then |
|---|
| 342 | print*,'kmax=',kmax |
|---|
| 343 | print*,'Setting values to zero in this band.' |
|---|
| 344 | kdist(k,j,i,band,1:(Nq+1))=0.0 |
|---|
| 345 | else |
|---|
| 346 | |
|---|
| 347 | if(kmin.lt.1e-30)then |
|---|
| 348 | kmin=1.0e-30 |
|---|
| 349 | endif |
|---|
| 350 | |
|---|
| 351 | !-------------------------------------------- |
|---|
| 352 | ! calculate delta nu values |
|---|
| 353 | nuwband(:)=0.0 |
|---|
| 354 | do ib=2,(ic-1) |
|---|
| 355 | nuwband(ib)=(nuband(ib+1)-nuband(ib-1))/2 |
|---|
| 356 | enddo |
|---|
| 357 | nuwband(1)=nuwband(2) |
|---|
| 358 | nuwband(ic)=nuwband(ic-1) ! causes tiny boundary error; whatever |
|---|
| 359 | nuwbtot=sum(nuwband) |
|---|
| 360 | |
|---|
| 361 | lkmax=log10(kmax) |
|---|
| 362 | lkmin=log10(kmin) |
|---|
| 363 | |
|---|
| 364 | dlogK=(lkmax-lkmin)/(Ngres-1) |
|---|
| 365 | |
|---|
| 366 | !-------------------------------------------- |
|---|
| 367 | ! calculate the g function |
|---|
| 368 | gfine(:)=0.0 |
|---|
| 369 | do ig=1,Ngres |
|---|
| 370 | klev=lkmin + dble(ig)*dlogK |
|---|
| 371 | klev=10**klev |
|---|
| 372 | klevs(ig)=klev |
|---|
| 373 | |
|---|
| 374 | addline(:)=0 |
|---|
| 375 | do ib=1,ic |
|---|
| 376 | if(kband(ib).lt.klev)then |
|---|
| 377 | addline(ib)=1.0 |
|---|
| 378 | endif |
|---|
| 379 | enddo |
|---|
| 380 | gfine(ig)=sum(addline(:)*nuwband(:))/nuwbtot |
|---|
| 381 | |
|---|
| 382 | enddo |
|---|
| 383 | print*,'gnorm = ',gfine(Ngres) |
|---|
| 384 | |
|---|
| 385 | !-------------------------------------------- |
|---|
| 386 | ! invert the g function |
|---|
| 387 | ik1=1 |
|---|
| 388 | ik2=1 |
|---|
| 389 | kdist_temp(:)=0.0 |
|---|
| 390 | do quad=1,Nq |
|---|
| 391 | |
|---|
| 392 | ik2=ik1 |
|---|
| 393 | |
|---|
| 394 | gfind: do while (gfine(ik2).lt.gw(quad)) |
|---|
| 395 | ik2=ik2+1 |
|---|
| 396 | if(ik2.gt.Ngres)then |
|---|
| 397 | print*,'g-inversion overflow...' |
|---|
| 398 | ik2=Ngres |
|---|
| 399 | exit gfind |
|---|
| 400 | endif |
|---|
| 401 | enddo gfind |
|---|
| 402 | |
|---|
| 403 | !if(ik2.gt.Ngres)then ! avoid out-of-bounds on final step |
|---|
| 404 | ! ik2=Ngres |
|---|
| 405 | !endif |
|---|
| 406 | |
|---|
| 407 | |
|---|
| 408 | if(.not.(ik1.ge.ik2))then ! this could probably be improved |
|---|
| 409 | !kdist_temp(quad)=sum(klevs(ik1:ik2))/(ik2-ik1) ! a bug was here! |
|---|
| 410 | kdist_temp(quad)=sum(klevs(ik1:ik2))/(ik2+1-ik1) |
|---|
| 411 | |
|---|
| 412 | |
|---|
| 413 | elseif(quad.ne.1)then |
|---|
| 414 | kdist_temp(quad)=kdist_temp(quad-1) |
|---|
| 415 | else |
|---|
| 416 | kdist_temp(quad)=kmin |
|---|
| 417 | endif |
|---|
| 418 | |
|---|
| 419 | !----------------------------------------- |
|---|
| 420 | ! assign to 5D matrix and convert to cm^2 / molecule |
|---|
| 421 | kdist(k,j,i,band,quad)=1.3807e-21 * & |
|---|
| 422 | kdist_temp(quad) * T(k)/P(j) |
|---|
| 423 | |
|---|
| 424 | ik1=ik2 |
|---|
| 425 | enddo ! g-space |
|---|
| 426 | kdist(k,j,i,band,Nq+1)=0.0 |
|---|
| 427 | |
|---|
| 428 | truemean = sum(kband(1:ic)*nuwband(1:ic))/sum(nuwband(1:ic)) |
|---|
| 429 | distmean = sum(kdist_temp(:)*w(:))/sum(w) |
|---|
| 430 | etot = 100*abs(truemean-distmean)/(distmean+truemean) |
|---|
| 431 | print*,'truemean = ',truemean |
|---|
| 432 | print*,'distmean = ',distmean |
|---|
| 433 | !print*,'kdisttemp = ',kdist_temp(Nq) |
|---|
| 434 | !print*,'kdistmax = ',kmax |
|---|
| 435 | print*,'etot (%) = ',etot |
|---|
| 436 | !print*,'w = ',w |
|---|
| 437 | |
|---|
| 438 | endif ! if kmax .le. 0.0 |
|---|
| 439 | |
|---|
| 440 | enddo ! bands |
|---|
| 441 | |
|---|
| 442 | enddo |
|---|
| 443 | enddo |
|---|
| 444 | enddo |
|---|
| 445 | |
|---|
| 446 | !-------------------------------------------------------- |
|---|
| 447 | ! save data and clean up |
|---|
| 448 | close(17) |
|---|
| 449 | open(14,file='corrk_gcm.dat',form='formatted') |
|---|
| 450 | write(14,*) kdist |
|---|
| 451 | close(14) |
|---|
| 452 | print*,' ' |
|---|
| 453 | print*,'Correlated-k data saved in corrk_gcm.dat' |
|---|
| 454 | print*,' ' |
|---|
| 455 | print*,'Now you must change the values of L_NTREF etc. in radinc_h.F90 (if necessary)' |
|---|
| 456 | print*,'and the variable "corrkdata" in callphys.def' |
|---|
| 457 | |
|---|
| 458 | deallocate(kdist_temp) |
|---|
| 459 | deallocate(kdist) |
|---|
| 460 | |
|---|
| 461 | end program generate_kmatrix |
|---|
