! module guydat implicit none ! ! Data from Guy Brassier (formerly in guydata.h): ! This data is read from netcdf file TGCM.glbmean.guyprof.nc ! (mss:/TGCM/glbmean/guyprof.nc) by subroutine rdguydat (see below). ! ! character(len=32) :: guydat_ncfile = 'TGCM.glbmean.guyprof.nc' ! ! Order of 24 temp-independent species (fsig(nsigbin,nfsig)): ! O2 H2O CO2 CH4 HO2 HCl ClO HOCl ! HO2NO2 CHF2Cl CHFCl2 NO3 HONO CH3ONO2 CH3CHO PAN ! CH3OOH OClO Cl2O2 Cl2 ClNO2 BrONO2 BrCl HOBr ! ! Order of the 20 temp-dependent species (fsigt(nsigbin,nfsigt)): ! O3 N2O H2O2 NO2 CH2O HNO3 CF2Cl2 CFCl3 ! CCl4 CH3CCl3 CH3Cl ClONO2 N2O5 C2F3Cl3 O3-1D C2HF3Cl2 ! C2H3FCl2CH2O-2H CH2O-H2 CH3Br ! ! Temperatures at which temp-dependent species are available (sigtn): ! (routine interp_sigt will interpolate species to required tn within ! this range) 180.00 200.00 220.00 240.00 260.00 280.00 (deg K) ! integer,parameter :: | nfsig=24, ! number of temp-independent species | nfsigt=20, ! number of temp-dependent species | nsigt=6, ! number of available temperatures | nsigbin=170 ! number of wavelength bins ! real :: | fsig(nsigbin,nfsig), ! temp-independent fields | fsigt(nsigbin,nfsigt,nsigt), ! temp-dependent fields | sigtn(nsigt) ! tn (deg K) for dependent fields ! character(len=8) :: | fsig_names(nfsig), ! names of tn-independent fields | fsigt_names(nfsigt) ! names of tn-dependent fields contains !----------------------------------------------------------------------- subroutine rdguydat use input,only: guydat_ncfile #include ! ! Read Guy Brassier data from netcdf file (called from main glbmean). ! ! Local: integer :: i,it,ip,istat,isig,isigt real :: fsigtrd(nsigbin,nsigt) ! local for reads character(len=8) :: name character(len=120) :: char120 integer :: ncid,ndims,ndimsvar,nvars,ngatts,itype,idunlim,natts, | idimids(NF_MAX_VAR_DIMS),idimlens(NF_MAX_VAR_DIMS),start1(1), | count1(1),start2(2),count2(2) ! ! Open netcdf file: write(6,"(/,72('-'))") write(6,"('rdguydat (guydat.F): guydat_ncfile=',a)") | trim(guydat_ncfile) call expand_path(guydat_ncfile) ! expand any env vars write(6,"('guydat_ncfile after expand_path = ',a)") | trim(guydat_ncfile) istat = nf_open(guydat_ncfile,NF_NOWRITE,ncid) if (istat /= NF_NOERR) then write(char120,"('>>> Error return from nf_open for netcdf ', | 'file ',a,' istat=',a)") trim(guydat_ncfile),istat call handle_ncerr(istat,char120) return endif write(6,"('rdguydat: Opened netcdf file ',a,' ncid=',i3)") | trim(guydat_ncfile),ncid ! ! Get number of dims, vars, and atts: istat = nf_inq(ncid,ndims,nvars,ngatts,idunlim) isig = 0 isigt = 0 ! ! Loop over all vars in the file: field_loop: do ip=1,nvars idimids = -1 ; idimlens = -1 ! ! Get info about this variable: istat = nf_inq_var(ncid,ip,name,itype,ndimsvar,idimids,natts) ! write(6,"('rdguydat: ip=',i3,' Field ',a,' ndims=',i3)") ! | ip,name,ndimsvar if (ndimsvar > 0) then do i=1,ndimsvar istat = nf_inq_dimlen(ncid,idimids(i),idimlens(i)) ! write(6,"(' dimension ',i2,': len=',i3)") i,idimlens(i) enddo endif if (ndimsvar<=0) cycle field_loop ! ! Read temperature-independent field: if (ndimsvar==1.and.idimlens(1)==nsigbin) then isig = isig+1 if (isig <= nfsig) then start1(1) = 1 count1(1) = nsigbin istat = nf_get_vara_double(ncid,ip,start1,count1, | fsig(1,isig)) if (istat /= NF_NOERR) then write(char120,"('>>> rdvar: Error return from ', | 'nf_get_vara_double for temp-indep field ',a)") name call handle_ncerr(istat,char120) return endif write(6,"('Read temp-independent field ',a,' isig=',i3, | ' min,max=',2e12.4)") name,isig,minval(fsig(:,isig)), | maxval(fsig(:,isig)) else write(6,"('>>> rdguydat: too many temperature-independent', | ' files: isig=',i3)") isig cycle field_loop endif ! ! Read available temperatures: elseif (ndimsvar==1.and.idimlens(1)==nsigt) then start1(1) = 1 count1(1) = nsigt istat = nf_get_vara_double(ncid,ip,start1,count1,sigtn) if (istat /= NF_NOERR) then write(char120,"('>>> rdvar: Error return from ', | 'nf_get_vara_double for field ',a)") name call handle_ncerr(istat,char120) return endif write(6,"(/,'Available temperatures ',a,'=',(6f8.2))") | name,sigtn ! ! Read temperature-dependent field: elseif (ndimsvar==2.and. | idimlens(1)==nsigbin.and.idimlens(2)==nsigt) then isigt = isigt+1 if (isigt <= nfsigt) then start2(:) = 1 count2(1) = nsigbin count2(2) = nsigt istat = nf_get_vara_double(ncid,ip,start2,count2,fsigtrd) if (istat /= NF_NOERR) then write(char120,"('>>> rdvar: Error return from ', | 'nf_get_vara_double for temp-dep field ',a)") name call handle_ncerr(istat,char120) return endif do it=1,nsigt fsigt(:,isigt,it) = fsigtrd(:,it) enddo write(6,"('Read temp-dependent field ',a,' isigt=',i3, | ' min,max=',2e12.4)") name,isigt, | minval(fsigt(:,isigt,:)),maxval(fsigt(:,isigt,:)) else write(6,"('>>> rdguydat: too many temperature-dependent', | ' files: isigt=',i3)") isigt cycle field_loop endif else write(6,"('rdguydat: not reading variable ',a)") name endif enddo field_loop ! ip=1,nvars write(6,"(72('-')/)") end subroutine rdguydat !------------------------------------------------------------------- c subroutine interp_sigt(f,ifget,tn,loginterp,spval) c c Get Guy Brassier temperature dependent species, and interpolate c to given temperature tn. (Species were read by rdguydat.f, and c are held in common in guydata.h). c Return interpolated fields in f(*,ip) when ifget(ip) > 0. c (if ifget(ip) le 0 set f(*,ip) = spval) c Do log (exponential) interpolation if loginterp > 0, otherwise linear c If tn is <= sigt(1), return data at sigt(1), c likewise if tn >= sigt(nsigt) then return data at sigt(nsigt) c If other error occurs in loops, then return f(i,ip) = spval c ! Args: integer,intent(in) :: ifget(nfsigt) ! input get flags integer,intent(in) :: loginterp real,intent(in) :: tn,spval real,intent(out) :: f(nsigbin,nfsigt) ! output ! ! Local: integer :: k0,k1,it,ip,i real :: f1,exparg c c Initialize output: c f(:,:) = spval c c Check validity of tn input: c c if (tn.lt.sigtn(1).or.tn.gt.sigtn(nsigt)) then c write(6,"(/,'>>>>> WARNING interp_sigt: requested temperature ', c + f10.3,' is out of range.',/,' tn must be >= ',f7.2, c + ' and <= ',f7.2)") tn,sigtn(1),sigtn(nsigt) c write(6,"(' Available temperatures from which to ', c + 'interpolate to tn are: ',/(8f9.2))") sigtn c write(6,"(' Returning all fields = ',e12.4)") spval c write(6,"('>>>>> END WARNING from interp_sigt'/)") c return c endif c c If tn is out of range, return data at lowest or highest available c temperature, as appropriate: c if (tn.le.sigtn(1)) then f(:,:) = fsigt(:,:,1) return endif if (tn.ge.sigtn(nsigt)) then f(:,:) = fsigt(:,:,nsigt) return endif c c Bracket tn (if tn == an available sigtn, then return data at that tn): c k0 = 0 k1 = 0 do it=1,nsigt-1 if (tn.ge.sigtn(it).and.tn.le.sigtn(it+1)) then if (tn.eq.sigtn(it)) then ! interpolation not necessary f(:,:) = fsigt(:,:,it) return endif if (tn.eq.sigtn(it+1)) then ! interpolation not necessary f(:,:) = fsigt(:,:,it+1) return endif k0 = it k1 = it+1 endif enddo c c Do linear interpolation: c if (loginterp.le.0) then do ip=1,nfsigt ! number of temp-dep species if (ifget(ip).gt.0) then ! is a requested species do i=1,nsigbin ! number of wavelength bins if (fsigt(i,ip,k0).eq.0..and.fsigt(i,ip,k1).eq.0.) then f(i,ip) = 0. else f1 = (tn-sigtn(k0)) / (sigtn(k1)-sigtn(k0)) f(i,ip) = f1*fsigt(i,ip,k1) + (1.-f1)*fsigt(i,ip,k0) endif enddo endif enddo c c Do log interpolation: c (check for low densities to avoid divide by small number) c else do ip=1,nfsigt ! number of temp-dep species if (ifget(ip).gt.0) then ! is a requested species do i=1,nsigbin ! number of wavelength bins if (fsigt(i,ip,k0).eq.0..and.fsigt(i,ip,k1).eq.0.) then f(i,ip) = 0. elseif (fsigt(i,ip,k0).ge.0..and.fsigt(i,ip,k0).le.1.e-20) + then ! density too low for log interp c write(6,"('>>> WARNING interp_sigt: ',a, c + ' density ',e12.4,' too low for log interp:')") c + fsigt_names(ip),fsigt(i,ip,k0) c write(6,"(' wavelength bin:',f8.2,' to ',f8.2, c + ' tn=',f9.2)") wv_low(i),wv_high(i),tn f(i,ip) = spval else ! exponential interp exparg = (alog(fsigt(i,ip,k1) / fsigt(i,ip,k0)) / + (sigtn(k1)-sigtn(k0))) * (tn-sigtn(k0)) f(i,ip) = fsigt(i,ip,k0)*exp(exparg) endif enddo endif enddo endif end subroutine interp_sigt !------------------------------------------------------------------- end module guydat