c kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum c Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet c c This file must be used under the terms of the CeCILL license. c This source file is licensed as described in the file COPYING, which c you should have received as part of this distribution. The terms c are also available at c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt c subroutine write_bandk(Nb,band, & nu_low,nu_hi,m,i,nk,nu, & Nmol,P,T,x,mol_name, & band_sigma,band_k, & sens_P,sens_T,sens_x, & band_dsigma_dP,band_dsigma_dT,band_dsigma_dx, & band_dk_dP,band_dk_dT,band_dk_dx) implicit none c c Purpose: to write the array of k coefficients into the appropriate file c c Inputs: c + Nb: number of narrowbands c + band: index of the narrowband to discretize c + nu_low: lower limit of each narrowband [cm¯¹] c + nu_hi: upper limit of each narrowband [cm¯¹] c + m : number of atmospheric levels c + i: index of the atmospheric level c + nk: number of wavenumber values in the band c + nu: array of wavenumber values [cm¯¹] c + Nmol: number of molecular species c + P: array of pressure levels c + T: array of temperature levels c + x: array of molecular species concentrations c + mol_name: array of molecular species names c + band_sigma: array of cross-section values [cm²/molecule] c + band_k: array of k values [m¯¹] c + sens_P: option for computation of sensitivities to total pressure c + sens_T: option for computation of sensitivities to temperature c + sens_x: option for computation of sensitivities to concentrations c + band_dsigma_dP: sensitivity of sigma to total pressure c + band_dsigma_dT: sensitivity of sigma to temperature c + band_dsigma_dx: sensitivity of sigma to species concentrations c + band_dk_dP: sensitivity of k to total pressure c + band_dk_dT: sensitivity of k to temperature c + band_dk_dx: sensitivity of k to species concentrations c include 'max.inc' include 'formats.inc' integer nl integer pass,Nb,band,m,i,kmin double precision nu_low(1:Nbmx) double precision nu_hi(1:Nbmx) integer nk,ik double precision nu(1:Nkmx) integer Nmol double precision P(1:Nmax) double precision T(1:Nmax) double precision x(1:Nmax,1:Nmol_max) character*10 mol_name(1:Nmol_max) double precision band_sigma(1:Nkmx) double precision band_k(1:Nkmx) logical sens_P,sens_T,sens_x double precision band_dsigma_dP(1:Nkmx) double precision band_dsigma_dT(1:Nkmx) double precision band_dsigma_dx(1:Nkmx,1:Nmol_max) double precision band_dk_dP(1:Nkmx) double precision band_dk_dT(1:Nkmx) double precision band_dk_dx(1:Nkmx,1:Nmol_max) character*(Nchar_mx) kfile,kfile_dP,kfile_dT,kfile_dx,nlines_file double precision sigma_min,k_min,dsigma_min,dk_min parameter(sigma_min=1.0D-99) parameter(k_min=1.0D-99) parameter(dsigma_min=1.0D-99) parameter(dk_min=1.0D-99) integer mol character*(Nchar_mx) lab_line character*(Nchar_mx) lab_lineP character*(Nchar_mx) lab_lineT character*(Nchar_mx) lab_linex integer wh logical ex integer strlen character*(Nchar_mx) label label='subroutine write_bandk' lab_line='nu (cm^{-1}) / sigma (cm^{2}/molecule) / k (m^{-1})' lab_lineP='nu (cm^{-1}) / d(sigma)/d(P) (cm^{2}/molecule/atm)' & //' / d(k)/d(P) (m^{-1}/atm)' lab_lineT='nu (cm^{-1}) / d(sigma)/d(T) (cm^{2}/molecule/K)' & //' / d(k)/d(T) (m^{-1}/K)' lab_linex='nu (cm^{-1}) / d(sigma)/d(x) (cm^{2}/molecule)' & //' / d(k)/d(x) (m^{-1})' call result_file_name(i,kfile,kfile_dP,kfile_dT,kfile_dx) c Debug c write(*,*) 'In:',label(1:strlen(label)) c write(*,*) 'band_dk_dT=',band_dk_dT(1) c write(*,*) 'band_dk_dP=',band_dk_dP(1) c Debug inquire(file=kfile(1:strlen(kfile)),exist=ex) if (ex) then wh=0 nlines_file='./nlines' call get_nlines(kfile,nlines_file,nl) else wh=1 nl=0 endif open(11,file=kfile(1:strlen(kfile)) & ,access='append') c write header in the case the file has just been created if (wh.eq.1) then write(11,10) 'Pressure (atm):' write(11,50) P(i) write(11,10) 'Temperature (K):' write(11,50) T(i) write(11,10) 'Number of molecules:' write(11,*) Nmol write(11,*) do mol=1,Nmol write(11,23) 'Molecule index:',mol write(11,10) 'Name:' write(11,10) mol_name(mol)(1:strlen(mol_name(mol))) write(11,10) 'Concentration:' write(11,50) x(i,mol) enddo ! mol write(11,*) write(11,10) lab_line(1:strlen(lab_line)) endif c do ik=1,nl c read(11,*) c enddo if (nl.eq.0) then kmin=1 else kmin=2 endif do ik=kmin,nk if (band_sigma(ik).lt.sigma_min) then band_sigma(ik)=0.0D+0 endif if (band_k(ik).lt.k_min) then c Debug c write(*,*) 'band_k(',ik,')=',band_k(ik),' -> 0' c Debug band_k(ik)=0.0D+0 endif write(11,52) nu(ik),band_sigma(ik),band_k(ik) enddo close(11) if (nl.eq.0) then write(*,*) kfile(1:strlen(kfile)) & ,' has been successfully generated' else write(*,*) kfile(1:strlen(kfile)) & ,' has been successfully updated' endif c results: dsigma_dP and dk_dP if (sens_P) then inquire(file=kfile(1:strlen(kfile_dP)),exist=ex) if (ex) then wh=0 nlines_file='./nlines' call get_nlines(kfile_dP,nlines_file,nl) else wh=1 nl=0 endif open(12,file=kfile_dP(1:strlen(kfile_dP)) & ,access='append') c write header in the case the file has just been created if (wh.eq.1) then write(12,10) 'Pressure (atm):' write(12,50) P(i) write(12,10) 'Temperature (K):' write(12,50) T(i) write(12,10) 'Number of molecules:' write(12,*) Nmol write(12,*) do mol=1,Nmol write(12,23) 'Molecule index:',mol write(12,10) 'Name:' write(12,10) mol_name(mol)(1:strlen(mol_name(mol))) write(12,10) 'Concentration:' write(12,50) x(i,mol) enddo ! mol write(12,*) write(12,10) lab_lineP(1:strlen(lab_lineP)) endif c do ik=1,nl c read(12,*) c enddo if (nl.eq.0) then kmin=1 else kmin=2 endif do ik=kmin,nk if (dabs(band_dsigma_dP(ik)).lt.dsigma_min) then band_dsigma_dP(ik)=0.0D+0 endif if (dabs(band_dk_dP(ik)).lt.dk_min) then band_dk_dP(ik)=0.0D+0 endif write(12,52) nu(ik),band_dsigma_dP(ik),band_dk_dP(ik) enddo close(12) if (nl.eq.0) then write(*,*) kfile_dP(1:strlen(kfile_dP)) & ,' has been successfully generated' else write(*,*) kfile_dP(1:strlen(kfile_dP)) & ,' has been successfully updated' endif endif c results: dsigma_dT and dk_dT if (sens_T) then inquire(file=kfile(1:strlen(kfile_dT)),exist=ex) if (ex) then wh=0 nlines_file='./nlines' call get_nlines(kfile_dT,nlines_file,nl) else wh=1 nl=0 endif open(13,file=kfile_dT(1:strlen(kfile_dT)) & ,access='append') c write header in the case the file has just been created if (wh.eq.1) then write(13,10) 'Pressure (atm):' write(13,50) P(i) write(13,10) 'Temperature (K):' write(13,50) T(i) write(13,10) 'Number of molecules:' write(13,*) Nmol write(13,*) do mol=1,Nmol write(13,23) 'Molecule index:',mol write(13,10) 'Name:' write(13,10) mol_name(mol)(1:strlen(mol_name(mol))) write(13,10) 'Concentration:' write(13,50) x(i,mol) enddo ! mol write(13,*) write(13,10) lab_lineT(1:strlen(lab_lineT)) endif c do ik=1,nl c read(13,*) c enddo if (nl.eq.0) then kmin=1 else kmin=2 endif do ik=kmin,nk if (dabs(band_dsigma_dT(ik)).lt.dsigma_min) then band_dsigma_dT(ik)=0.0D+0 endif if (dabs(band_dk_dT(ik)).lt.dk_min) then band_dk_dT(ik)=0.0D+0 endif write(13,52) nu(ik),band_dsigma_dT(ik),band_dk_dT(ik) enddo close(13) if (nl.eq.0) then write(*,*) kfile_dT(1:strlen(kfile_dT)) & ,' has been successfully generated' else write(*,*) kfile_dT(1:strlen(kfile_dT)) & ,' has been successfully updated' endif endif c results: dsigma_dx and dk_dx if (sens_x) then inquire(file=kfile(1:strlen(kfile_dx)),exist=ex) if (ex) then wh=0 nlines_file='./nlines' call get_nlines(kfile_dx,nlines_file,nl) else wh=1 nl=0 endif open(14,file=kfile_dx(1:strlen(kfile_dx)) & ,access='append') c write header in the case the file has just been created if (wh.eq.1) then write(14,10) 'Pressure (atm):' write(14,50) P(i) write(14,10) 'Temperature (K):' write(14,50) T(i) write(14,10) 'Number of molecules:' write(14,*) Nmol write(14,*) do mol=1,Nmol write(14,23) 'Molecule index:',mol write(14,10) 'Name:' write(14,10) mol_name(mol)(1:strlen(mol_name(mol))) write(14,10) 'Concentration:' write(14,50) x(i,mol) enddo ! mol write(14,*) write(14,10) lab_linex(1:strlen(lab_linex)) endif c do ik=1,nl c read(14,*) c enddo if (nl.eq.0) then kmin=1 c write(14,*) Nmol else kmin=2 endif do ik=kmin,nk do mol=1,Nmol if (dabs(band_dsigma_dx(ik,mol)).lt.dsigma_min) then band_dsigma_dx(ik,mol)=0.0D+0 endif if (dabs(band_dk_dx(ik,mol)).lt.dk_min) then band_dk_dx(ik,mol)=0.0D+0 endif enddo ! mol write(14,52) nu(ik), & (band_dsigma_dx(ik,mol),mol=1,Nmol), & (band_dk_dx(ik,mol),mol=1,Nmol) enddo close(14) if (nl.eq.0) then write(*,*) kfile_dx(1:strlen(kfile_dx)) & ,' has been successfully generated' else write(*,*) kfile_dx(1:strlen(kfile_dx)) & ,' has been successfully updated' endif endif return end subroutine result_file_name(i,kfile,kfile_dP,kfile_dT,kfile_dx) implicit none include 'max.inc' include 'formats.inc' integer strlen,i character*3 kch character*(Nchar_mx) kf,kf_dP,kf_dT,kf_dx character*(Nchar_mx) base,kfile,kfile_dP,kfile_dT,kfile_dx base='./results/' call num2str3(i,kch) kf='k'//kch(1:strlen(kch)) kf_dP='dk_dP'//kch(1:strlen(kch)) kf_dT='dk_dT'//kch(1:strlen(kch)) kf_dx='dk_dx'//kch(1:strlen(kch)) kfile=base(1:strlen(base))//kf(1:strlen(kf)) kfile_dP=base(1:strlen(base))//kf_dP(1:strlen(kf_dP)) kfile_dT=base(1:strlen(base))//kf_dT(1:strlen(kf_dT)) kfile_dx=base(1:strlen(base))//kf_dx(1:strlen(kf_dx)) return end