module gwf ! ! B. Foster: April, 2009 ! Read netcdf file containing Fx,Fy forcing from Sharon Vadas' ray ! trace program. This data is interpolated to model grid and time, ! and used in mgw.F. ! use params_module,only: | model_nlon => nlon, model_nlat => nlat, model_nlev => nlev implicit none #include save ! ! Spatial dimensions: integer,parameter :: | nlon=111, ! number of longitudes (nx) | nlat= 55, ! number of latitudes (ny) | nlev= 51, ! number of levels (nz) | ntime=54 ! number of times ! ! Spatial coordinate arrays of gwf data: real :: lon(nlon), lat(nlat), lev(nlev), time(ntime) ! ! Forcing fu,fv (nlon,nlat,nlev,ntime) real,save,dimension(nlon,nlat,nlev,ntime) :: | fu, ! Force in zonal direction (fx) | fv ! Force in meridional direction (fy) contains !----------------------------------------------------------------------- subroutine gwf_read(filename) implicit none ! ! Read netcdf gw force file. ! ! Args: character(len=*),intent(in) :: filename ! ! Local: integer :: istat,ncid,id,nlon_rd,nlat_rd,nlev_rd,ntime_rd ! ! Open file for reading: istat = nf_open(filename,NF_NOWRITE,ncid) if (istat /= NF_NOERR) | call handle_ncerr(istat,'Error from nf_open',1) write(6,"(/,72('-'))") write(6,"('gwf nc_read: Opened file ',a,' for reading.')") | trim(filename) istat = nf_inq_dimid(ncid,"lon",id) istat = nf_inq_dimlen(ncid,id,nlon_rd) if (nlon_rd /= nlon) then write(6,"('>>> gwf nc_read: nlon_rd=',i4,' but nlon=',i4)") | nlon_rd,nlon call shutdown('gwf nc_read') endif istat = nf_inq_varid(ncid,"lon",id) istat = nf_get_var_double(ncid,id,lon) istat = nf_inq_dimid(ncid,"lat",id) istat = nf_inq_dimlen(ncid,id,nlat_rd) if (nlat_rd /= nlat) then write(6,"('>>> gwf nc_read: nlat_rd=',i4,' but nlat=',i4)") | nlat_rd,nlat call shutdown('gwf nc_read') endif istat = nf_inq_varid(ncid,"lat",id) istat = nf_get_var_double(ncid,id,lat) istat = nf_inq_dimid(ncid,"lev",id) istat = nf_inq_dimlen(ncid,id,nlev_rd) if (nlev_rd /= nlev) then write(6,"('>>> gwf nc_read: nlev_rd=',i4,' but nlev=',i4)") | nlev_rd,nlev call shutdown('gwf nc_read') endif istat = nf_inq_varid(ncid,"lev",id) istat = nf_get_var_double(ncid,id,lev) istat = nf_inq_dimid(ncid,"time",id) istat = nf_inq_dimlen(ncid,id,ntime_rd) if (ntime_rd /= ntime) then write(6,"('>>> gwf nc_read: ntime_rd=',i4,' but ntime=',i4)") | ntime_rd,ntime call shutdown('gwf nc_read') endif istat = nf_inq_varid(ncid,"time",id) istat = nf_get_var_double(ncid,id,time) write(6,"('gwf nc_read: nlon =',i4,' lon min,max=',2f8.2)") | nlon,minval(lon),maxval(lon) write(6,"('gwf nc_read: nlat =',i4,' lat min,max=',2f8.2)") | nlat,minval(lat),maxval(lat) write(6,"('gwf nc_read: nlev =',i4,' lev min,max=',2f8.2)") | nlev,minval(lev),maxval(lev) write(6,"('gwf nc_read: ntime=',i4,' time min,max=',2f8.2)") | ntime,minval(time),maxval(time) istat = nf_inq_varid(ncid,"FU",id) istat = nf_get_var_double(ncid,id,fu) write(6,"('gwf nc_read: fu min,max=',2e12.4)") | minval(fu),maxval(fu) istat = nf_inq_varid(ncid,"FV",id) istat = nf_get_var_double(ncid,id,fv) write(6,"('gwf nc_read: fv min,max=',2e12.4)") | minval(fv),maxval(fv) istat = nf_close(ncid) write(6,"(72('-'),/)") end subroutine gwf_read !----------------------------------------------------------------------- subroutine gwf_interp(secs,lev0,lev1,lon0,lon1,lat0,lat1) use init_module,only: glon,glat,zpint use fields_module,only: itc,z ! full subdomain(k,i,j,2) implicit none ! ! Args: real,intent(in) :: secs integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1 ! ! Local: integer :: i,j,k,n,npts_in_gwf integer :: i0,j0,k0,it0 integer :: i1,j1,k1,it1 real :: | lons_in_gwf(model_nlon), | lats_in_gwf(model_nlat), | levs_in_gwf(model_nlev), | hts(lev0:lev1) real :: fracx,fracy,gwsfu_k0,gwsfu_k1 real,parameter :: spv=1.e36 real,parameter :: earthrad = 6.37122e8 ! earth radius (cm) ! ! Secondary forcing output ! (this will be global 3d data, used by mgw) real,dimension(lev0:lev1,lon0:lon1,lat0:lat1) :: | gwsfu,gwsfv ! ! Find 3d coords of global model grid points that fall within the ! gwf data box. ! lons_in_gwf = spv ! init lats_in_gwf = spv levs_in_gwf = spv npts_in_gwf = 0 ! ! Function bracket_1d returns j0, where glat(j) is between lat(j0) ! and lat(j0+1). If not bracketed, zero is returned. Bracket_1d ! is also used for lon and lev. ! do j=1,model_nlat j0 = bracket_1d(lat,nlat,glat(j)) if (j0==0) cycle do i=1,model_nlon i0 = bracket_1d(lon,nlon,glon(i)) if (i0==0) cycle ! ! If i,j is in my subdomain, check for vertical coord: if (i >= lon0 .and. i <= lon1 .and. | j >= lat0 .and. j <= lat1) then ! ! Convert geopotential (cm) to geometric ht (km) at this i,j location: hts(:) = (z(lev0:lev1,i,j,itc)* | (1.+z(lev0:lev1,i,j,itc)/earthrad))*1.e-5 ! write(6,"('gwf_interp: i=',i4,' j=',i4,' hts=',/, ! | (6e12.4))") i,j,hts do k=lev0,lev1 k0 = bracket_1d(lev,nlev,hts(k)) if (k0==0) cycle ! ! Save coords of this grid point, which is in my 3d subdomain: lats_in_gwf(j) = glat(j) lons_in_gwf(i) = glon(i) levs_in_gwf(k) = zpint(k) ! zp at interfaces npts_in_gwf = npts_in_gwf+1 enddo ! k=lev0,lev1 endif ! i,j is in my subdomain enddo ! i=1,model_nlon enddo ! j=1,model_nlat ! ! If no subdomain grid points are within the data, then no interpolation. if (npts_in_gwf==0) then ! write(6,"('gwf_interp: found no 3d subdomain grid ', ! | 'points inside gwf data coords')") return endif ! ! Report to stdout any subdomain points found within gwf data box: write(6,"('gwf_interp: found ',i4,' 3d subdomain grid points ', | ' within gwf data box:')") npts_in_gwf gwsfu = 0. ! init of 3d model subdomain interpolated output gwsfv = 0. n = 0 do j=lat0,lat1 do i=lon0,lon1 do k=lev0,lev1 if (lats_in_gwf(j) /= spv .and. lons_in_gwf(i) /= spv .and. | levs_in_gwf(k) /= spv) then n = n+1 j0 = bracket_1d(lat,nlat,glat(j)) i0 = bracket_1d(lon,nlon,glon(i)) k0 = bracket_1d(lev,nlev,hts(k)) j1 = j0+1 ; i1 = i0+1 ; k1 = k0+1 ! ! Bilinear interpolation to i,j (data are in i0,i1,j0,j1,k0,k1): if (abs(lon(i1)-lon(i0)) <= 1.e-20) then fracx = 1. else fracx = (glon(i)-lon(i0)) / abs(lon(i1)-lon(i0)) endif ! write(6,"('gwf_interp: i=',i4,' i0,1=',2i4,' glon(i)=', ! | f9.3,' lon(i0)=',f9.3,' lon(i1)=',f9.2,' fracx=', ! | e12.4)") i,i0,i1,glon(i),lon(i0),lon(i1),fracx if (abs(lat(j1)-lat(j0)) <= 1.e-20) then fracy = 1. else fracy = (glat(j)-lat(j0)) / abs(lat(j1)-lat(j0)) endif ! write(6,"('gwf_interp: j=',i4,' j0,1=',2i4,' glat(j)=', ! | f9.3,' lat(j0)=',f9.3,' lat(j1)=',f9.2,' fracy=', ! | e12.4)") j,j0,j1,glat(j),lat(j0),lat(j1),fracy ! ! Forcing data fu,fv (nlon,nlat,nlev,ntime) it0 = 10 ! temp for test gwsfu_k0 = | fracx*(fracy*fu(i1,j1,k0,it0)+ | (1.-fracy)*fracx*fu(i1,j0,k0,it0))+ | (1.-fracx)*(fracy*fu(i0,j1,k0,it0)+ | (1.-fracy)*(1.-fracx)*fu(i0,j0,k0,it0)) ! write(6,"('gwf_interp: n=',i5,' i0,1=',2i4,' j0,1=', ! | 2i4,' k0=',i4,' fu(i1,j1,k0)=',e12.4,' fu(i1,j0,k0)=', ! | e12.4,' fu(i0,j1,k0)=',e12.4,' fu(i0,j0,k0)=',e12.4)") ! | n,i0,i1,j0,j1,k0,fu(i1,j1,k0,it0),fu(i1,j0,k0,it0), ! | fu(i0,j1,k0,it0),fu(i0,j0,k0,it0) gwsfu_k1 = | fracx*(fracy*fu(i1,j1,k1,it0)+ | (1.-fracy)*fracx*fu(i1,j0,k1,it0))+ | (1.-fracx)*(fracy*fu(i0,j1,k1,it0)+ | (1.-fracy)*(1.-fracx)*fu(i0,j0,k1,it0)) ! ! Interpolate to height (output to gwsfu model subdomain) fracx = (hts(k)-lev(k0)) / (lev(k1)-lev(k0)) gwsfu(k,i,j) = fracx*gwsfu_k1 + (1.-fracx)*gwsfu_k0 write(6,"('gwf_interp: n=',i5,' model lon,lat=',2f9.3, | ' ht=',f9.3,' gwsfu_k0,k1=',2e12.4,' gwsfu(k,i,j)=', | e12.4)") n,glon(i),glat(j),hts(k),gwsfu_k0,gwsfu_k1, | gwsfu(k,i,j) endif ! model grid point is within 3d data box enddo ! k=lev0,lev1 enddo ! i=lon0,lon1 enddo ! j=lat0,lat1 end subroutine gwf_interp !----------------------------------------------------------------------- integer function bracket_1d(v,nv,x) ! ! If x lies between 2 adjacent points in v(i-1) and v(i), return i-1, ! otherwise return 0. ! integer,intent(in) :: nv real,intent(in) :: v(nv),x integer :: i ! bracket_1d = 0 if (nv < 2) then write(6,"('>>> bracket_1d: nv must be >= 2: nv=',i5)") nv return endif do i=2,nv if ((x >= v(i-1) .and. x <= v(i)).or. | (x <= v(i-1) .and. x >= v(i))) then bracket_1d = i-1 return endif enddo end function bracket_1d !----------------------------------------------------------------------- subroutine handle_ncerr(istat,msg,ifatal) implicit none ! ! Handle a netcdf lib error: ! integer,intent(in) :: istat,ifatal character(len=*),intent(in) :: msg ! write(6,"(/72('-'))") write(6,"('>>> Error from netcdf library:')") write(6,"(a)") trim(msg) write(6,"('istat=',i5)") istat write(6,"(a)") nf_strerror(istat) write(6,"(72('-')/)") if (ifatal > 0) stop 'handle_ncerr' return end subroutine handle_ncerr !----------------------------------------------------------------------- end module gwf