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) ! ! 4-d interpolation of gwf data to model subdomains. ! use init_module,only: glon,glat use fields_module,only: itc,z ! full subdomain(k,i,j,2) use addfld_module,only: addfld implicit none ! ! Args: real,intent(in) :: secs integer,intent(in) :: lev0,lev1,lon0,lon1,lat0,lat1 ! ! Local: integer :: i,j,k,npts_in_gwf integer :: i0,j0,k0,it0 integer :: i1,j1,k1,it1 real :: hts(lev0:lev1) ! geometric heights (km) real :: f0,gwsfu_k0,gwsfu_k1,gwsfv_k0,gwsfv_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 ! ! Init 3d model subdomain interpolated output: gwsfu = 0. ! whole array gwsfv = 0. it0 = 10 ! temp for test ! ! Loop over full global model domain, and interpolate data only ! when 3d model coords fall within the data block. ! npts_in_gwf = 0 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 ! ! Continue if i,j is in my subdomain: 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 ! ! Loop over model vertical dimension: do k=lev0,lev1 k0 = bracket_1d(lev,nlev,hts(k)) if (k0==0) cycle ! ! Current grid point is within gwf data, and in my 3d subdomain. npts_in_gwf = npts_in_gwf+1 j1 = j0+1 ; i1 = i0+1 ; k1 = k0+1 ! ! Bilinear interpolation of data at k0,k1 to i,j ! (data are in i0,i1,j0,j1,k0,k1): ! gwsfu_k0 = bilin(fu(:,:,k0,it0),lon,lat,nlon,nlat, | i0,i1,j0,j1,glon(i),glat(j)) gwsfu_k1 = bilin(fu(:,:,k1,it0),lon,lat,nlon,nlat, | i0,i1,j0,j1,glon(i),glat(j)) gwsfv_k0 = bilin(fv(:,:,k0,it0),lon,lat,nlon,nlat, | i0,i1,j0,j1,glon(i),glat(j)) gwsfv_k1 = bilin(fv(:,:,k1,it0),lon,lat,nlon,nlat, | i0,i1,j0,j1,glon(i),glat(j)) ! ! Interpolate to height (output to gwsfu,v in model subdomain) f0 = (hts(k)-lev(k0)) / (lev(k1)-lev(k0)) gwsfu(k,i,j) = f0*gwsfu_k1 + (1.-f0)*gwsfu_k0 gwsfv(k,i,j) = f0*gwsfv_k1 + (1.-f0)*gwsfv_k0 write(6,"('gwf_interp: n=',i5,' model lon,lat=',2f9.3, | ' ht=',f9.3,' gwsfu=',e12.4,' gwsfv=',e12.4)") | npts_in_gwf,glon(i),glat(j),hts(k),gwsfu(k,i,j), | gwsfv(k,i,j) enddo ! k=lev0,lev1 endif ! i,j is in my subdomain enddo ! i=1,model_nlon enddo ! j=1,model_nlat ! ! Save to secondary history: do j=lat0,lat1 call addfld('GWSFU',' ',' ',gwsfu(lev0:lev1,lon0:lon1,j), | 'lev',lev0,lev1,'lon',lon0,lon1,j) call addfld('GWSFV',' ',' ',gwsfv(lev0:lev1,lon0:lon1,j), | 'lev',lev0,lev1,'lon',lon0,lon1,j) enddo end subroutine gwf_interp !----------------------------------------------------------------------- real function bilin(f,x,y,nx,ny,ix0,ix1,iy0,iy1,xx,yy) ! ! Bilinear interpolation of f(nx,ny) to xx,yy. ! ! Args: integer,intent(in) :: nx,ny real,intent(in) :: f(nx,ny) real,intent(in) :: x(nx),y(ny),xx,yy integer,intent(in) :: ix0,ix1,iy0,iy1 ! ! Local: real :: fracx,fracy ! bilin = 0. if (abs(x(ix1)-x(ix0)) <= 1.e-20) then fracx = 1. else fracx = (xx-x(ix0)) / abs(x(ix1)-x(ix0)) endif if (abs(y(iy1)-y(iy0)) <= 1.e-20) then fracy = 1. else fracy = (yy-y(iy0)) / abs(y(iy1)-y(iy0)) endif ! ! Bilinear interpolation: bilin = fracx*(fracy*f(ix1,iy1)+ | (1.-fracy)*fracx*f(ix1,iy0))+ | (1.-fracx)*(fracy*f(ix0,iy1)+ | (1.-fracy)*(1.-fracx)*f(ix0,iy0)) end function bilin !----------------------------------------------------------------------- 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