module read_ncfile ! ! Read netcdf input file. ! use netcdf use shr_kind_mod ,only: r8 => shr_kind_r8 ! 8-byte reals use params ,only: iulog,finit,fillvalue,pi,dtr,kbotdyn,zbotdyn,kboltz ! use my_mpi ,only: lev0,lev1,lon0,lon1,lat0,lat1 use namelist ,only: ctpoten_user => ctpoten,inLat,inLon use geogrid ,only: & ntime,nlon,nlat,nlev,nilev,time,glat,glon,ylatg,ylong,& dlatg,dlong,zlev,zilev,lev_sequence,lon_sequence, & jspole,jnpole,nlonp1,nlonp2,nlatp1,nmlon,nmlat,pMid,pInt,p0,pver,pverp ! use maggrid ,only: zmlev,nmlev implicit none save integer :: ncid ! netcdf file id integer,parameter :: nf=19 ! number of input fields character(len=8) :: fnames(nf) ! field names as on input file integer,parameter :: nf2D=9 ! number of input fields character(len=8) :: fnames2D(nf2D) ! 2D field names as on input file ! ! Inputs to edynamo on 3d distributed geographic subdomains ! (nlev,nlon,nlat) ! These are read from the input history file for each time ! ! real(r8),allocatable,dimension(:,:,:),save :: & ! tn, & ! neutral temperature (deg K) ! un, & ! neutral zonal wind (cm/s) ! vn, & ! neutral meridional wind (cm/s) ! omega, & ! vertical motion (/sec) ! zpot, & ! geopotential height (cm) ! barm, & ! mean molecular weight ! ped, & ! pedersen conductivity ! hall ! hall conductivity ! ! Inputs to WACCM-X routine on 3d distributed geographic subdomains ! These are read from the input history file for each time ! real(r8),allocatable,dimension(:,:,:),save :: & tNIn, & ! neutral temperature (deg K) uNIn, & ! neutral zonal wind (cm/s) vNIn, & ! neutral meridional wind (cm/s) o1, & ! O number density (cm-3) o2, & ! O2 number density (cm-3) nO, & ! NO number density (cm-3) h, & ! H number density (cm-3) oP, & ! O+ number density (cm-3) o2P, & ! O2+ number density (cm-3) nOP, & ! NO+ number density (cm-3) nE, & ! Electon number density (cm-3) omegaIn, & ! vertical motion (/sec) tEIn, & ! electron temperature (deg K) tIIn, & ! ion temperature (deg K) zGm, & ! geometric altitude zGp, & ! geopotential altitude eX, & ! eastward component of electric field (V) eY, & ! northward component of electric field (V) eZ ! upward component of electric field (V) real(r8),allocatable,dimension(:,:),save :: & ed1, & ! zonal electric field ed2, & ! meridional electric field zenAng, & ! zenith angle bEastIn, & ! eastward magnetic field (nT) bNorthIn,& ! northward magnetic field (nT) bDownIn, & ! upward magnetic field (nT) bMagIn, & ! magnitude of magnetic field (nT) gLatM, & ! geomagnetic latitude on geographic grid gLonM ! geomagnetic longitude on geographic grid ! ! 2D variables used in WACCM-X routines ! real(r8),allocatable,dimension(:,:),save :: & tN, & ! neutral temperature (deg K) uN, & ! neutral temperature (deg K) vN, & ! neutral temperature (deg K) mmrO1, & ! atomic oxygen mass mixing ratio mmrO2, & ! molecular oxygen mass mixing ratio mmrNO, & ! nitric oxide mass mixing ratio ndensOp, & ! O+ ion number density ndenso2p,& ! O2+ ion number density ndensnop,& ! NO+ ion number density ndensE, & ! electron number density omega, & ! vertical motion (sec-1) tE, & ! electron temperature (deg K) tI ! ion temperature (deg K) ! ! 1D variables used in WACCM-X routines ! real(r8),allocatable,dimension(:),save :: & ed1_geo, & ! eastward component of electric field on geographic coords (V) ed2_geo ! northward component of electric field on geographic coords (V) ! ! Scalar variables used in WACCM-X routines ! real(r8),save :: & zenAngR, & ! zenith angle(radians) geoMagLat, & ! geomagnetic latitude (radians) geoMagLon, & ! geomagnetic longitude (radians) bMag ! magnitude of magnetic field (nT) real(r8) :: year,ut,ctpoten,f107,bEast,bNorth,bDown ! integer :: iday ! ! TIMEGCM variables needed in output module: ! integer :: mtime(3) ! TIMEGCM model time (day,hour,min) real(r8),allocatable,dimension(:,:) :: & ! (nlon,nlat) tlbc, & ! TN lower boundary ulbc, & ! UN lower boundary vlbc ! VN lower boundary real(r8) :: lbc ! bottom interface level (scalar) ! ! WACCM variables needed on output file: ! integer :: & date ,& ! current date (YYYYMMDD) datesec ,& ! current seconds of current date nsteph ,& ! current timestep ndcur ,& ! current day (from base day) nscur ! current seconds of current day real(r8),parameter :: eps = 1.e-6_r8 private eps contains !----------------------------------------------------------------------- subroutine get_geogrid(ncfile,model) ! ! Open netcdf input file, and get dimensions and coordinate variables ! of 3d global geographic grid, and time. This is called once per run, ! and works for either TGCM or WACCM files. ! ! Args: character(len=*),intent(in) :: & ncfile, & ! file path model ! model name ! ! Local: integer :: istat,idunlim,id, ilev integer :: idv_lon,idv_lat,idv_lev,idv_ilev,idv_time character(len=1024) :: msg ! message for error handler character(len=NF90_MAX_NAME) :: varname ! ! Open netcdf history file. ! istat = nf90_open(ncfile,NF90_NOWRITE,ncid) if (istat /= NF90_NOERR) then write(msg,"('Error opening file ',a)") trim(ncfile) call handle_ncerr(istat,trim(msg),1) else write(iulog,"(/,'Opened ',a,' file ',a)") model,trim(ncfile) ! write(iulog,"(/,'Opened file ',a)") trim(ncfile) endif ! ! Get number of times on the file (length of unlimited variable): ! istat = nf90_inq_dimid(ncid,'time',idunlim) istat = nf90_inquire_dimension(ncid,idunlim,varname,ntime) ! ! Get grid dimensions from the file: ! istat = nf90_inq_dimid(ncid,'lon',id) istat = nf90_inquire_dimension(ncid,id,varname,nlon) istat = nf90_inq_dimid(ncid,'lat',id) istat = nf90_inquire_dimension(ncid,id,varname,nlat) istat = nf90_inq_dimid(ncid,'lev',id) istat = nf90_inquire_dimension(ncid,id,varname,nlev) istat = nf90_inq_dimid(ncid,'ilev',id) istat = nf90_inquire_dimension(ncid,id,varname,nilev) istat = nf90_inq_dimid(ncid,'mlon',id) istat = nf90_inquire_dimension(ncid,id,varname,nmlon) istat = nf90_inq_dimid(ncid,'mlat',id) istat = nf90_inquire_dimension(ncid,id,varname,nmlat) ! ! Set vertical grid names for WACCM-X ! pver = nlev pverp = nlev + 1 ! ! Allocate coordinate variables: ! allocate(glon(nlon)) allocate(glat(nlat)) allocate(zlev(nlev)) allocate(zilev(nilev)) allocate(time(ntime)) allocate(pMid(1,nlev)) allocate(pInt(1,nilev+1)) ! ! Read coordinate vars: ! istat = nf90_inq_varid(ncid,'lon',idv_lon) istat = nf90_get_var(ncid,idv_lon,glon,(/1/),(/nlon/)) istat = nf90_inq_varid(ncid,'lat',idv_lat) istat = nf90_get_var(ncid,idv_lat,glat,(/1/),(/nlat/)) istat = nf90_inq_varid(ncid,'lev',idv_lev) istat = nf90_get_var(ncid,idv_lev,zlev,(/1/),(/nlev/)) istat = nf90_inq_varid(ncid,'ilev',idv_ilev) istat = nf90_get_var(ncid,idv_ilev,zilev,(/1/),(/nilev/)) istat = nf90_inq_varid(ncid,'time',idv_time) istat = nf90_get_var(ncid,idv_time,time,(/1/),(/ntime/)) istat = nf90_inq_varid(ncid,'p0',id) istat = nf90_get_var(ncid,id,p0) ! ! Get pressure from TIME-GCM ln(P0/P) vertical grid ! ionBot and oPBot are bottom vertical levels to calculate Te/Ti and O+/e respectively. ! ionBot should correspond to about 60 km or 0.5 hPa. ! oPBot should correspond to about 130 km or 0.00001 hPa. ! write(*,*) 'Converting pressure from ln(p0/p) to Pa' do ilev = 1, nlev pMid(1,ilev) = p0/1000._r8 * exp(-zlev(ilev)) * 100._r8 enddo do ilev = 1, nilev pInt(1,ilev) = p0/1000._r8 * exp(-zilev(ilev)) * 100._r8 enddo ! write(iulog,"('get_geogrid: pMid before reverse ',2e12.4)") pMid(1,1), pMid(1,nlev) ! write(iulog,"('get_geogrid: pInt before reverse ',2e12.4)") pInt(1,1), pInt(1,nilev+1) ! ! Reverse order of vertical dimension: ! call reverse_vec(pMid(1,:),nlev) call reverse_vec(pInt(1,:),nilev+1) pInt(1,1) = pint(1,2) - (pInt(1,3) - pInt(1,2)) ! write(iulog,"('get_geogrid: pMid after reverse ',2e12.4)") pMid(1,1), pMid(1,nlev) ! write(iulog,"('get_geogrid: pInt after reverse ',2e12.4)") pInt(1,1), pInt(1,nilev+1) ! ! write(iulog,*) 'get_geogrid: pInt after reverse all ', pInt(1,:) ! write(iulog,"('get_geogrid: pMid min,max=',2e12.4)") minval(pMid), maxval(pMid) write(iulog,"('get_geogrid: pInt min,max=',2e12.4)") minval(pInt), maxval(pInt) ! ! Model-independent sizes/dimensions: nlonp1 = nlon+1 nlonp2 = nlon+2 nlatp1 = nlat+1 ! write(iulog,"('Read ',a,' input file: nlon=' ,i4,' glon=' ,/,(6f13.4))") model,nlon,glon ! write(iulog,"('Read ',a,' input file: nlat=' ,i4,' glat=' ,/,(6f13.4))") model,nlat,glat ! write(iulog,"('Read ',a,' input file: nlev=' ,i4,' zlev=' ,/,(6e13.4))") model,nlev,zlev ! write(iulog,"('Read ',a,' input file: nilev=',i4,' zilev=',/,(6e13.4))") model,nilev,zilev ! write(iulog,"('Read ',a,' input file: ntime=',i4,' time =',/,(6e13.4))") model,ntime,time write(iulog,"('get_geogrid for ',a,' ntime=',i4,' nlon=',i4,' nlat=',i4,' nlev=',i4,' nilev=',i4)") & trim(model),ntime,nlon,nlat,nlev,nilev end subroutine get_geogrid !----------------------------------------------------------------------- subroutine read_tgcm_wx(ncfile,itime) ! ! Args: character(len=*),intent(in) :: ncfile integer,intent(in) :: itime ! ! Local: integer :: i,ii,j,k,istat,id,n,start(4),start2D(3),count(4),count2D(3),lonbeg,lonend,looplat,looplon integer :: iyear,iyear1(1),iday1(1),indLat,indLon,iLt,iLn real(r8),allocatable,save :: f3d(:,:,:) ! global field (lev,lat,lon) real(r8),allocatable,save :: f2D(:,:) ! 2D field (lat,lon) real(r8),allocatable,save :: f2DMag(:,:) ! 2D mag field (mlat,mlon) real(r8) :: ut1(1),ctpoten1(1),f1071(1) real(r8) :: real8,fmin,fmax logical,save :: first = .true. fnames = (/'TN','UN','VN','O1','O2','NO','H','OP','O2P','NOP','NE',& 'OMEGA','TE','TI','ZG','Z','EX','EY','EZ'/) fnames2D = (/'ED1','ED2','CHI','bx','by','bz',& 'bmag','glatm','glonm'/) istat = nf90_inq_varid(ncid,'mtime',id) istat = nf90_get_var(ncid,id,mtime,(/1,itime/),(/3,1/)) istat = nf90_inq_varid(ncid,'day',id) istat = nf90_get_var(ncid,id,iday1,(/itime/),(/1/)) iday = iday1(1) istat = nf90_inq_varid(ncid,'ut',id) istat = nf90_get_var(ncid,id,ut1,(/itime/),(/1/)) ut = ut1(1) istat = nf90_inq_varid(ncid,'year',id) istat = nf90_get_var(ncid,id,iyear1,(/itime/),(/1/)) iyear = iyear1(1) year = real(iyear) istat = nf90_inq_varid(ncid,'f107d',id) istat = nf90_get_var(ncid,id,f1071,(/itime/),(/1/)) f107 = f1071(1) ! ! If ctpoten was not provided by user via namelist, ! then read from input file: ! ! if (ctpoten_user == fillvalue) then ! read from file istat = nf90_inq_varid(ncid,'ctpoten',id) istat = nf90_get_var(ncid,id,ctpoten1,(/itime/),(/1/)) ctpoten = ctpoten1(1) ! else ! use from namelist ! ctpoten = ctpoten_user ! endif ! write(iulog,"('read_tgcm: itime=',i4,' mtime=',3i4,' ut=',f8.3,' iday=',i4,' year=',f7.0)") & ! itime,mtime,ut,iday,year ! ! Allocate a 3d/2D global fields for reading: if (first) then call alloc_inputs ! allocate(tlbc(nlon,nlat)) ! TN lower boundary ! allocate(ulbc(nlon,nlat)) ! UN lower boundary ! allocate(vlbc(nlon,nlat)) ! VN lower boundary allocate(f3d(nlon,nlat,nlev)) ! 3d global read array f3d = finit allocate(f2D(nlon,nlat)) ! 2D read array f2D = finit allocate(f2DMag(nmlon,nmlat)) ! 2D mag coords read array f2DMag = finit ! ! Set horizontal geographic grid in radians (for apex code): allocate(ylatg(0:nlat+1)) ! include poles for TIMEGCM allocate(ylong(nlon+1)) ! single periodic point real8 = real(nlat) ; dlatg = pi/real8 real8 = real(nlon) ; dlong = 2._r8*pi/real8 do j=1,nlat real8 = real(j-1) ylatg(j) = -0.5_r8*(pi-dlatg)+real8*dlatg enddo ylatg(0) = -pi/2._r8+eps ylatg(nlat+1) = pi/2._r8-eps do i=1,nlon+1 real8 = real(i-1) ylong(i) = -pi+real8*dlong enddo ! ! j-index to geographic poles (timegcm): jspole = 0 jnpole = nlat+1 ! write(iulog,"('read_tgcm: nlat=',i4,' ylatg(0:nlat+1)=',/,(6e12.4))") nlat,ylatg ! write(iulog,"('read_tgcm: nlon=',i4,' ylong(1:nlon+1)=',/,(6e12.4))") nlon,ylong endif ! first call ! ! Read global input fields at current time: ! do n=1,nf2D istat = nf90_inq_varid(ncid,fnames2D(n),id) if (istat /= NF90_NOERR) then write(iulog,"('>>> Error getting var id for field ',a)") & fnames2D(n) cycle endif ! ! Read variable into 2D global array: ! start2D(1:2) = 1 start2D(3) = itime if (fnames2D(n) == 'ED1' .or. fnames2D(n) == 'ED2') then count2D(1) = nmlon count2D(2) = nmlat count2D(3) = 1 looplat = nmlat looplon = nmlon write(*,*) 'nmlat, nmlon ', nmlat, nmlon istat = nf90_get_var(ncid,id,f2DMag,start=start2D,count=count2D) if (istat /= NF90_NOERR) then write(iulog,"('>>> Error reading var ')") fnames2D(n) call handle_ncerr(istat,'Error from nf90_get_var',1) endif else if (fnames2D(n) == 'CHI') then count2D(1) = nlon count2D(2) = nlat count2D(3) = 1 looplat = nlat looplon = nlon istat = nf90_get_var(ncid,id,f2D,start=start2D,count=count2D) if (istat /= NF90_NOERR) then write(iulog,"('>>> Error reading var ')") fnames2D(n) call handle_ncerr(istat,'Error from nf90_get_var',1) endif else looplat = nlat looplon = nlon istat = nf90_get_var(ncid,id,f2D) if (istat /= NF90_NOERR) then write(iulog,"('>>> Error reading var ')") fnames2D(n) call handle_ncerr(istat,'Error from nf90_get_var',1) endif endif if (first) write(iulog,"('Read TIMEGCM field ',a,' 2D global min,max=',2e12.4)") & fnames2D(n),minval(f2D),maxval(f2D) ! ! Loop and handle periodic points and store data ! do j=1,looplat do i=1,looplon ii = i ! don't need to worry about periodic points for just one column ! ii = i-2 ! body of longitude domain ! if (i <= 2) ii = looplon-2+i ! i=1,2 -> ii= ! west end periodic points ! if (i >= looplon+3) ii = i-looplon-2 ! i= -> ii=1,2 ! east end periodic points if (fnames2D(n)=='ED1 ') ed1(i,j) = f2DMag(ii,j) if (fnames2D(n)=='ED2 ') ed2(i,j) = f2DMag(ii,j) if (fnames2D(n)=='CHI ') zenAng(i,j) = f2D(ii,j) if (fnames2D(n)=='bx ') bEastIn(i,j) = f2D(ii,j) if (fnames2D(n)=='by ') bNorthIn(i,j) = f2D(ii,j) if (fnames2D(n)=='bz ') bDownIn(i,j) = -f2D(ii,j) if (fnames2D(n)=='bmag ') bMagIn(i,j) = f2D(ii,j) if (fnames2D(n)=='glatm ') gLatM(i,j) = f2D(ii,j) if (fnames2D(n)=='glonm ') gLonM(i,j) = f2D(ii,j) enddo enddo ! lon loop enddo ! nf2D lat loop if (first) then write(iulog,"(/,'read_tgcm_wx (first call): nlev=',i4,' nlon=',i4,' nlat=',i4)") & nlev,nlon,nlat write(iulog,"('read_tgcm_wx: ED1 subdomain min,max=',2e12.4)") minval(ed1), maxval(ed1) write(iulog,"('read_tgcm_wx: ED2 subdomain min,max=',2e12.4)") minval(ed2), maxval(ed2) write(iulog,"('read_tgcm_wx: CHI subdomain min,max=',2e12.4)") minval(zenAng), maxval(zenAng) write(iulog,"('read_tgcm_wx: bx subdomain min,max=',2e12.4)") minval(bEastIn),maxval(bEastIn) write(iulog,"('read_tgcm_wx: by subdomain min,max=',2e12.4)") minval(bNorthIn), maxval(bNorthIn) write(iulog,"('read_tgcm_wx: bz subdomain min,max=',2e12.4)") minval(bDownIn), maxval(bDownIn) write(iulog,"('read_tgcm_wx: bmag subdomain min,max=',2e12.4)") minval(bMagIn), maxval(bMagIn) write(iulog,"('read_tgcm_wx: glatm subdomain min,max=',2e12.4)") minval(gLatM), maxval(gLatM) write(iulog,"('read_tgcm_wx: glonm subdomain min,max=',2e12.4)") minval(gLonM), maxval(gLatM) endif ! ! Read global input fields at current time: ! do n=1,nf istat = nf90_inq_varid(ncid,fnames(n),id) if (istat /= NF90_NOERR) then write(iulog,"('>>> Error getting var id for field ',a)") & fnames(n) cycle endif ! ! Read variable into 3d global array: ! start(1:3) = 1 start(4) = itime count(1) = nlon count(2) = nlat count(3) = nlev count(4) = 1 istat = nf90_get_var(ncid,id,f3d,start=start,count=count) if (istat /= NF90_NOERR) then write(iulog,"('>>> Error reading var ')") fnames(n) call handle_ncerr(istat,'Error from nf90_get_var',1) endif if (first) write(iulog,"('Read TIMEGCM field ',a,' 3d global min,max=',2e12.4)") & fnames(n),minval(f3d),maxval(f3d) ! ! Transfer to geographic subdomains: ! ! lonbeg = lon0 ! if (lon0==1) lonbeg = 3 ! lonend = lon1 ! if (lon1==nlon+4) lonend = lon1-2 ! ! do j=lat0,lat1 ! do i=lonbeg,lonend ! ! Reverse order of vertical dimension (as in zlev above): ! do j=1,nlat do i=1,nlon call reverse_vec(f3d(i,j,:),nlev) enddo enddo do j=1,nlat do i=1,nlon ii = i ! don't need to worry about periodic points for just one column ! ii = i-2 ! body of longitude domain ! if (i <= 2) ii = nlon-2+i ! i=1,2 -> ii=71,72 ! west end periodic points ! if (i >= nlon+3) ii = i-nlon-2 ! i=75,76 -> ii=1,2 ! east end periodic points ! Subscript #3 of the array F3D has value 50 which is greater than the upper bound of 49 ! if (fnames(n)=='TN_DYN ') tn(:,i,j) = f3d(ii,j,:) ! if (fnames(n)=='UN_DYN ') un(:,i,j) = f3d(ii,j,:) ! if (fnames(n)=='VN_DYN ') vn(:,i,j) = f3d(ii,j,:) ! if (fnames(n)=='W_DYN ') omegaIn(:,i,j) = f3d(ii,j,:) ! if (fnames(n)=='Z_DYN ') zpot(:,i,j) = f3d(ii,j,:) ! if (fnames(n)=='BARM_DYN') barm(:,i,j) = f3d(ii,j,:) ! if (fnames(n)=='PED_DYN ') ped(:,i,j) = f3d(ii,j,:) ! if (fnames(n)=='HAL_DYN ') hall(:,i,j) = f3d(ii,j,:) if (fnames(n)=='TN ') tNIn(:,i,j) = f3d(ii,j,:) if (fnames(n)=='UN ') uNIn(:,i,j) = f3d(ii,j,:) if (fnames(n)=='VN ') vNIn(:,i,j) = f3d(ii,j,:) if (fnames(n)=='O1 ') o1(:,i,j) = f3d(ii,j,:) if (fnames(n)=='O2 ') o2(:,i,j) = f3d(ii,j,:) if (fnames(n)=='NO ') nO(:,i,j) = f3d(ii,j,:) if (fnames(n)=='H ') h(:,i,j) = f3d(ii,j,:) if (fnames(n)=='OP ') oP(:,i,j) = f3d(ii,j,:) if (fnames(n)=='O2P ') o2P(:,i,j) = f3d(ii,j,:) if (fnames(n)=='NOP ') nOP(:,i,j) = f3d(ii,j,:) if (fnames(n)=='NE ') nE(:,i,j) = f3d(ii,j,:) if (fnames(n)=='OMEGA ') omegaIn(:,i,j) = f3d(ii,j,:) if (fnames(n)=='TE ') tEIn(:,i,j) = f3d(ii,j,:) if (fnames(n)=='TI ') tIIn(:,i,j) = f3d(ii,j,:) if (fnames(n)=='ZG ') zGm(:,i,j) = f3d(ii,j,:) if (fnames(n)=='Z ') zGp(:,i,j) = f3d(ii,j,:) if (fnames(n)=='EX ') eX(:,i,j) = f3d(ii,j,:) if (fnames(n)=='EY ') eY(:,i,j) = f3d(ii,j,:) if (fnames(n)=='EZ ') eZ(:,i,j) = f3d(ii,j,:) enddo enddo ! j=lat0,lat1 ! ! Set kbotdyn (params.F90): ! ! if (fnames(n)=='ZG '.and.first) then ! kbotdyn = 0 ! do k=1,nlev ! fmin = minval(f3d(:,:,k)*1.e-5) ! minimum km ! fmax = maxval(f3d(:,:,k)*1.e-5) ! maximum km ! if (fmin < zbotdyn.and.fmax > zbotdyn) then ! kbotdyn = k ! exit ! endif ! if (fmin > zbotdyn.and.fmax > zbotdyn.and.k > 1) then ! kbotdyn = k-1 ! exit ! endif ! enddo ! k=1,nlev ! if (kbotdyn > 0) then ! write(iulog,"('Set kbotdyn=',i4,' for TIMEGCM: zlev(kbotdyn)=',e12.4,' zbotdyn=',f10.2,' Z min,max at kbotdyn=',2f10.2)") & ! kbotdyn,zlev(kbotdyn),zbotdyn,fmin,fmax ! else ! call shutdown('Could not set kbotdyn for TIMEGCM file') ! endif ! endif ! Z3 and first call enddo ! n=1,nf ! ! Report min,max to stdout if first timestep: ! if (first) then write(iulog,"(/,'read_tgcm (first call): nlev=',i4,' nlon=',i4,' nlat=',i4)") & nlev,nlon,nlat write(iulog,"('read_tgcm_wx: TN subdomain min,max=',2e12.4)") minval(tnIn), maxval(tnIn) write(iulog,"('read_tgcm_wx: UN subdomain min,max=',2e12.4)") minval(unIn), maxval(unIn) write(iulog,"('read_tgcm_wx: VN subdomain min,max=',2e12.4)") minval(vnIn), maxval(vnIn) write(iulog,"('read_tgcm_wx: OMEGA subdomain min,max=',2e12.4)") minval(omegaIn),maxval(omegaIn) write(iulog,"('read_tgcm_wx: O1 subdomain min,max=',2e12.4)") minval(o1), maxval(o1) write(iulog,"('read_tgcm_wx: O2 subdomain min,max=',2e12.4)") minval(o2), maxval(o2) write(iulog,"('read_tgcm_wx: NO subdomain min,max=',2e12.4)") minval(nO), maxval(nO) write(iulog,"('read_tgcm_wx: Op subdomain min,max=',2e12.4)") minval(oP), maxval(oP) write(iulog,"('read_tgcm_wx: TE subdomain min,max=',2e12.4)") minval(tEIn), maxval(tEIn) write(iulog,"('read_tgcm_wx: TI subdomain min,max=',2e12.4)") minval(tIIn), maxval(tIIn) ! write(iulog,"('read_tgcm: TLBC global min,max=',2e12.4)") minval(tlbc), maxval(tlbc) ! write(iulog,"('read_tgcm: ULBC global min,max=',2e12.4)") minval(ulbc), maxval(ulbc) ! write(iulog,"('read_tgcm: VLBC global min,max=',2e12.4)") minval(vlbc), maxval(vlbc) endif first = .false. ! ! Find latitude and longitude indices of namelist requested latitude/longitude ! do iLt = 1, nlat-1 if (glat(iLt+1) >= inLat .and. glat(iLt) <= inLat) indLat = iLt enddo do iLn = 1, nlon-1 if (glon(iLn+1) >= inLon .and. glon(iLn) <= inLon) indLon = iLn enddo write(iulog,*) 'inLat, indLat, glat(indLat+1), glat(indLat) ', inLat, indLat, glat(indLat+1), glat(indLat) write(iulog,*) 'inLon, indLon, glon(indLon-1), glon(indLon) ', inLon, indLon, glon(indLon+1), glon(indLon) ! ! Convert variables to single vertical column for input latitude/longitude ! zenAngR = zenAng(indLon,indLat) bEast = bEastIn(indLon,indLat) bNorth = bNorthIn(indLon,indLat) bDown = bDownIn(indLon,indLat) bMag = bMagIn(indLon,indLat) tN(1,:) = tNIn(:,indLon,indLat) uN(1,:) = uNIn(:,indLon,indLat) / 100._r8 vN(1,:) = vNIn(:,indLon,indLat) / 100._r8 mmrO1(1,:) = o1(:,indLon,indLat) mmrO2(1,:) = o2(:,indLon,indLat) mmrNO(1,:) = nO(:,indLon,indLat) ndensOp(1,:) = oP(:,indLon,indLat) ndensO2p(1,:) = o2P(:,indLon,indLat) ndensNOp(1,:) = nOP(:,indLon,indLat) ndensE(1,:) = NE(:,indLon,indLat) omega(1,:) = omegaIn(:,indLon,indLat) tE(1,:) = tEIn(:,indLon,indLat) tI(1,:) = tIIn(:,indLon,indLat) !! The upper BC for tE and tI Op and NOp do not look right from TIME-GCM tI(:,1) = 1.5*tI(:,2)-.5*tI(:,3) tE(:,1) = 1.5*tE(:,2)-.5*tE(:,3) ndensOp(:,1) = 1.5*ndensOp(:,2)-.5*ndensOp(:,3) write(iulog,*)'ndensOp top 3',ndensOp(1,1:3) ndensNOp(:,1) = 1.5*ndensNOp(:,2)-.5*ndensOp(:,3) write(iulog,*)'ndensNOp top 3',ndensNOp(1,1:3) ed1_geo(1) = eX(nlev/2,indLon,indLat) ed2_geo(1) = eY(nlev/2,indLon,indLat) geoMagLat = gLatM(indLon,indLat) geoMagLon = gLonM(indLon,indLat) write(iulog,*) zenAngR,bEast,bNorth,bDown,bMag write(iulog,"('read_tgcm_wx: TN converted min,max=',2e12.4)") minval(tn), maxval(tn) write(iulog,"('read_tgcm_wx: UN converted min,max=',2e12.4)") minval(un), maxval(un) write(iulog,"('read_tgcm_wx: VN converted min,max=',2e12.4)") minval(vn), maxval(vn) write(iulog,"('read_tgcm_wx: OMEGA converted min,max=',2e12.4)") minval(omega),maxval(omega) write(iulog,"('read_tgcm_wx: mmrO1 converted min,max=',2e12.4)") minval(mmrO1), maxval(mmrO1) write(iulog,"('read_tgcm_wx: mmrO2 converted min,max=',2e12.4)") minval(mmrO2), maxval(mmrO2) write(iulog,"('read_tgcm_wx: mmrNO converted min,max=',2e12.4)") minval(mmrNO), maxval(mmrNO) write(iulog,"('read_tgcm_wx: ndensOp converted min,max=',2e12.4)") minval(ndensOp), maxval(ndensOp) write(iulog,"('read_tgcm_wx: TE converted min,max=',2e12.4)") minval(tE), maxval(tE) write(iulog,"('read_tgcm_wx: TI converted min,max=',2e12.4)") minval(tI), maxval(tI) write(iulog,*) ed1_geo(1),ed2_geo(1),geoMagLat,geoMagLon end subroutine read_tgcm_wx !----------------------------------------------------------------------- !----------------------------------------------------------------------- !----------------------------------------------------------------------- subroutine alloc_inputs ! ! Allocate input fields on geographic subdomain, and init: ! ! write(iulog,"('allocate input fields on subdomains: lev0,1=',2i4,' lon0,1=',2i4,' lat0,1=',2i4)") & ! lev0,lev1,lon0,lon1,lat0,lat1 ! allocate(tn (nlev,nlon,nlat)) ! allocate(un (nlev,nlon,nlat)) ! allocate(vn (nlev,nlon,nlat)) ! allocate(omega(nlev,nlon,nlat)) ! allocate(zpot (nlev,nlon,nlat)) ! allocate(barm (nlev,nlon,nlat)) ! allocate(ped (nlev,nlon,nlat)) ! allocate(hall (nlev,nlon,nlat)) ! tn =finit ! un =finit ! vn =finit ! omega=finit ! zpot =finit ! barm =finit ! ped =finit ! hall =finit allocate(tNIn (nlev,nlon,nlat)) allocate(uNIn (nlev,nlon,nlat)) allocate(vNIn (nlev,nlon,nlat)) allocate(o1 (nlev,nlon,nlat)) allocate(o2 (nlev,nlon,nlat)) allocate(nO (nlev,nlon,nlat)) allocate(h (nlev,nlon,nlat)) allocate(oP (nlev,nlon,nlat)) allocate(o2P (nlev,nlon,nlat)) allocate(nOP (nlev,nlon,nlat)) allocate(nE (nlev,nlon,nlat)) allocate(omegaIn (nlev,nlon,nlat)) allocate(tEIn (nlev,nlon,nlat)) allocate(tIIn (nlev,nlon,nlat)) allocate(zGm (nlev,nlon,nlat)) allocate(zGp (nlev,nlon,nlat)) allocate(eX (nlev,nlon,nlat)) allocate(eY (nlev,nlon,nlat)) allocate(eZ (nlev,nlon,nlat)) allocate(ed1 (nmlon,nmlat)) allocate(ed2 (nmlon,nmlat)) allocate(zenAng (nlon,nlat)) allocate(bEastIn (nlon,nlat)) allocate(bNorthIn (nlon,nlat)) allocate(bDownIn (nlon,nlat)) allocate(bMagIn (nlon,nlat)) allocate(gLatM (nlon,nlat)) allocate(gLonM (nlon,nlat)) allocate(tN (1,nlev)) allocate(uN (1,nlev)) allocate(vN (1,nlev)) allocate(mmro1 (1,nlev)) allocate(mmro2 (1,nlev)) allocate(mmrNO (1,nlev)) allocate(ndensOp (1,nlev)) allocate(ndensO2p (1,nlev)) allocate(ndensNOp (1,nlev)) allocate(ndensE (1,nlev)) allocate(omega (1,nlev)) allocate(tE (1,nlev)) allocate(tI (1,nlev)) allocate(ed1_geo(1)) allocate(ed2_geo(1)) tNIn =finit uNIn =finit vNIn =finit o1 =finit o2 =finit nO =finit h =finit oP =finit o2P =finit nOP =finit nE =finit omegaIn =finit tEIn =finit tIIn =finit zenAng =finit zGm =finit zGp =finit ed1 =finit ed2 =finit zenAng =finit bEastIn =finit bNorthIn =finit bDownIn =finit bMagIn =finit gLatM =finit gLonM =finit tN =finit uN =finit vN =finit mmro1 =finit mmro2 =finit ndensoP =finit omega =finit tE =finit tI =finit end subroutine alloc_inputs !----------------------------------------------------------------------- end module read_ncfile