!------------------------------------------------------------------------- ! NASA GSFC Land Information Systems LIS 2.3 ! !------------------------------------------------------------------------- !BOP ! !ROUTINE: retgdas.F90 ! ! !DESCRIPTION: ! Defines forcing parameters, retrieves the fields using calls to ! getgb, and interpolates the fields to LDAS specifications ! ! !REVISION HISTORY: ! 14 Dec 2000: Urszula Jambor; Rewrote geteta.f to use GDAS forcing in GLDAS ! 15 Mar 2001: Jon Gottschalck; Added additional parameters and octets in ! which to search in GRIB files ! 01 Jun 2001: Urszula Jambor; Added option to get forcing from different ! files (F00 instantaneous and F06 six hour means) ! 29 Jan 2003: Urszula Jambor; Rewrote code, uses GETGB call to replace ! ungribgdas. Interpolation now occurs in interp_gdas. ! Using GETGB avoids problems with the Oct2002 GDAS ! grid update. ! 12 Nov 2003: Matt Rodell; Check to make sure input file exists before ! opening and thereby creating a new, empty file. ! 14 Nov 2003: Matt Rodell; Ensure lugb varies in call to baopen ! ! !INTERFACE: subroutine retgdas( order, ld, gindex, name, ferror,try ) ! !USES: use lis_module ! LDAS non-model-specific 1-D variables use lisdrv_module, only : lis use time_manager use baseforcing_module, only: glbdata1, glbdata2 use gdasdomain_module, only : gdasdrv implicit none ! !ARGUMENTS: type (lisdec) ld integer :: gindex(ld%d%lnc, ld%d%lnr) integer :: order ! 1 indicates lesser interp. bdry, 2 indicates greater character(len=80) :: name, nameF06 integer :: F00flag ! if 1, need for data from 2 files (name, nameF06) integer :: ferror ! set to zero if there's an error integer :: try !EOP !==== Local Variables======================= integer agrmetFlag, cmapFlag character(len=80) :: fname integer :: iv, c, r, ii, mx, my,i integer :: errorflag integer :: endloop, nforce integer :: j, lugb,iret,gbret,jret,jpds(200),jgds(200) integer :: lubi,kf,k,kpds(200),kgds(200) integer :: ngdas integer, dimension(gdasdrv%nmif) :: pds5, pds6, pds7, pds16 !logical*1, allocatable :: lb(:) logical*1 :: lb(gdasdrv%ncold*gdasdrv%nrold) logical :: file_exists !real, allocatable :: f(:) real :: f(gdasdrv%ncold*gdasdrv%nrold) !real, allocatable :: varfield(:,:) real :: varfield(ld%d%lnc,ld%d%lnr) real :: dswrf(ld%d%lnc,ld%d%lnr) real :: elev(ld%d%lnc,ld%d%lnr) real :: ism integer :: count ! IPOLATES VARIABLES, JESSE 20050225 integer :: ip, ipopt(20) integer :: kgdsi(25), kgdso(25) integer :: mi, mo integer, parameter :: km = 1 integer, parameter :: ibi = 1 integer :: no real, allocatable :: rlat(:) real, allocatable :: rlon(:) integer :: ibo logical*1, allocatable :: lo(:) !real, allocatable :: g(:) real :: g(gdasdrv%ncold*gdasdrv%nrold) !=== End Variable Definition ======================= !BOC print*,'MSG: Reading ',trim(lis%p%elevfile) open(30,file=lis%p%elevfile,form='unformatted',status='old') read(30) elev read(30) elev read(30) elev close(30) ngdas = (gdasdrv%ncold*gdasdrv%nrold) !-------------------------------------------------------------------------- ! Set the GRIB parameter specifiers ! FORCING() ARRAY: \\ ! 1. T 2m Temperature interpolated to 2 metres [$K$] \\ ! 2. q 2m Instantaneous specific humidity interpolated to 2 metres[$kg/kg$] \\ ! 3. radswg Downward shortwave flux at the ground [$W/m^2$] \\ ! 4. lwgdwn Downward longwave radiation at the ground [$W/m^2$] \\ ! 5. u 10m Instantaneous zonal wind interpolated to 10 metres [$m/s$] \\ ! 6. v 10m Instantaneous meridional wind interpolated to 10 metres[$m/s$] \\ ! 7. ps Instantaneous Surface Pressure [$Pa$] \\ ! 8. preacc Total precipitation [$mm/s$] \\ ! 9. gfrac Greenness fraction (0-1) ! 10. albedo Surface albedo (0-1) ! 11. zlvl Height of atmospheric forcing [m] ! 12. z0 Surface roughness [m] ! 13. ch Surface exchange coefficient [m/s] ! 14. t1 Skin Temperature [K] ! 15. sneqv SWE [m] ! 16. snowh Snow Depth [m] !-------------------------------------------------------------------------- !J if (get_nstep() .eq. 0) then !J pds5 = (/011,051,204,205,033,034,001,059,214,084,144,144, 011,011, 065/) !parameter !J pds7 = (/002,002,000,000,010,010,000,000,000,000,010,2760,010,2760,000/) !htlev2 !J nforce = gdasdrv%nmif !J else ! pds5 = (/ 011,051,204,205,033,034,001,059,087,084,007,083,208,011,065,066 /) !parameter ! pds6 = (/ 109,109,001,001,109,109,001,001,001,001,109,001,001,001,001,001 /) !level ! pds7 = (/ 001,001,000,000,001,001,000,000,000,000,001,000,000,000,000,000 /) !height ! pds16= (/ 010,010,010,010,010,010,010,003,010,003,010,010,010,010,010,010 /) !3-ave; 10-fcst !Jesse Meng 20191217 albedo=uswrf/dswrf pds5 = (/ 011,051,204,205,033,034,001,059,087,211,007,083,208,011,065,066/) !parameter pds6 = (/ 109,109,001,001,109,109,001,001,001,001,109,001,001,001,001,001/) !level pds7 = (/ 001,001,000,000,001,001,000,000,000,000,001,000,000,000,000,000/) !height pds16= (/ 010,010,010,010,010,010,010,003,010,010,010,010,010,010,010,010/) !3-ave; 10-fcst nforce = gdasdrv%nmif !J endif ferror = 1 !-------------------------------------------------------------------------- ! if there's a problem then ferror is set to zero !-------------------------------------------------------------------------- iv = 0 errorflag = 0 endloop = 0 fname = name read(fname(24:24),'(i1)') i read(fname(34:34),'(i1)') j read(fname(47:47),'(i1)') k inquire (file=fname, exist=file_exists) if (file_exists) then lugb = i+j+k*3 lubi = 0 !NOT INDEXED ! lubi = i+j+k*7 iret=0 call baopenr(lugb,fname,iret) print*,"J---retgdas()---baopenr()---lugb, iret= ",trim(fname),lugb,iret ! iret=0 ! call baopenr(lubi,trim(fname)//'.index',iret) ! print*,"J---retgdas()---baopenr()---lubi, iret= ",trim(fname)//'.index',lubi,iret else ferror = 0 return endif ! allocate(lb(gdasdrv%ncold*gdasdrv%nrold)) ! allocate( f(gdasdrv%ncold*gdasdrv%nrold)) mi = ngdas mo = ld%d%lnc*ld%d%lnr ! allocate(rlat(mo)) ! allocate(rlon(mo)) ! allocate( lo(mo)) ! allocate( g(mo)) ! allocate(varfield(ld%d%lnc,ld%d%lnr)) j=0 do iv = 1, nforce ! do ! if ( endloop == 1 ) exit ! iv = iv+1 ! fname = name ! inquire (file=fname, exist=file_exists) ! if (file_exists) then !!-------------------------------------------------------------------------- !! Set up to open file and retrieve specified field !!-------------------------------------------------------------------------- ! lugb = iv +try ! j = 0 jpds = -1 jpds(5) = pds5(iv) jpds(6) = pds6(iv) jpds(7) = pds7(iv) jpds(16)= pds16(iv) gbret=0 call getgb(lugb,lubi,ngdas,j,jpds,jgds,kf,k,kpds,kgds,lb,f,gbret) print*,"J---retgdas()---getgb()-----lugb,lubi,gbret= ", lugb,lubi,gbret ! call baopenr(lugb,fname,iret) ! print*,"J---retgdas()---baopenr()---iret= ", iret ! ! if(iret==0) then ! call getgb(lugb,lubi,ngdas,j,jpds,jgds,kf,k,kpds,kgds,lb,f,gbret) ! else ! gbret = 99 ! endif ! print*,"J---retgdas()---getgb()----gbret= ", gbret ! call baclose(lugb,jret) ! call getgb(lugb,lubi,ngdas,j,jpds,jgds,kf,k,kpds,kgds,lb,f,gbret) ! else ! ferror = 0 ! deallocate(f) ! deallocate(lb) ! deallocate(rlat) ! deallocate(rlon) ! deallocate(lo) ! deallocate(g) ! deallocate(varfield) ! return ! endif !-------------------------------------------------------------------------- ! If field successfully retrieved, interplate to LIS domain !-------------------------------------------------------------------------- ! JESSE 20050225 USE IPOLATES !-------------------------------------------------------------------------- if (gbret==0) then ip = 0 ipopt = 0 write(*,'(3I7)') iv, mi, mo write(*,'(11I7)') kpds(1:16) kgdsi(1:25) = kgds(1:25) write(*,'(11I7)') kgdsi(1:11) kgdso = 0 do i = 1, 10 if( i.EQ.4 .OR. i.EQ.7 .OR. i.EQ.8 .OR. i.EQ.9 ) then kgdso(i) = int(lis%d%gridDesc(i)*1000) else kgdso(i) = int(lis%d%gridDesc(i)) endif enddo kgdso(11) = 0 kgdso(20) = 255 write(*,'(11I7)') kgdso(1:11) !------------------------------------------------------ ! IPOLATE FORCING TO RUN DOMAIN !------------------------------------------------------ iret = 0 ! call ipolates (ip,ipopt,kgdsi,kgdso,mi,mo, & ! km,ibi,lb,f,no,rlat,rlon,ibo,lo,g,iret) g = f if(iret .NE. 0) then print*, "IPOLATES ERROR!! PROGRAM STOP!!" call exit(iret) end if !------------------------------------------------------ ! END JESSE 20050225 IPOLATES !------------------------------------------------------ count = 0 do r = ld%d%lnr, 1, -1 do c = 1, ld%d%lnc varfield(c,r) = g(c+count) if(iv==3) then dswrf(c,r)=varfield(c,r) endif if(iv==10) then !albedo=uswrfsfc/dswrfsfc*100. if(dswrf(c,r)<=0.) then varfield(c,r)=0. else varfield(c,r)=(varfield(c,r)/dswrf(c,r))*100. endif endif if(iv==11) then !covert hgt to hgt above sfc varfield(c,r)=varfield(c,r)-elev(c,r) endif end do count = count + ld%d%lnc end do else errorflag = 1 endif !endif (gbret==0) if ( errorflag == 1 ) then endloop = 1 ferror = 0 else count = 0 do r = 1, ld%d%lnr do c = 1, ld%d%lnc if (gindex(c,r).ne. -1) then if ( order == 1 ) then glbdata1(iv,gindex(c,r)) = varfield(c,r) else glbdata2(iv,gindex(c,r)) = varfield(c,r) end if endif end do count = count +ld%d%lnc end do end if if ( errorflag == 1 ) then print *, 'Could not find correct forcing parameter in file ',name end if if ( iv == nforce ) endloop = 1 end do jret=0 call baclose(lugb,jret) jret=0 call baclose(lubi,jret) !deallocate(lb) !deallocate(f) !deallocate(rlat) !deallocate(rlon) !deallocate(lo) !deallocate(g) !deallocate(varfield) return !EOC end subroutine retgdas !BOP ! !ROUTINE: interp_gdas ! ! !DESCRIPTION: ! This subroutine interpolates a given GDAS field ! to the LIS domain. Special treatment for some ! initialization fields. ! Code based on old ungribgdas.f ! ! !INTERFACE: subroutine interp_gdas(kpds,kgds,ngdas,f,lb,lis_gds,nc,nr, & varfield) ! !USES: ! use def_ipMod, only : w110,w120,w210,w220,n110,n120,n210,n220,& ! rlat0,rlon0, & ! w113,w123,w213,w223,n113,n123,n213,n223,rlat3,rlon3 use bilinear_interpMod use conserv_interpMod use gdasdomain_module, only : mi !EOP !w113,w123,w213,w223,n113,n123,n213,n223,rlat3, & ! rlon3 implicit none !=== Begin variable declarations integer :: nc, nr, ngdas, nglis integer :: kpds(200),kgds(200), lis_gds(200) integer :: ip, ipopt(20),ibi,km,iret integer :: no, ibo integer :: count,i,j,v real :: f(ngdas) real :: ism, udef real, dimension(nc,nr) :: varfield, geogtemp real, dimension(nc*nr) :: lis1d logical*1 :: geogmask(nc,nr) logical*1 :: lb(ngdas) logical*1 :: lo(nc*nr) !=== End variable declarations !BOC !----------------------------------------------------------------------- ! Setting interpolation options (ip=0,bilinear) ! (km=1, one parameter, ibi=1,use undefined bitmap ! (needed for soil moisture and temperature only) ! Use budget bilinear (ip=3) for precip forcing fields !----------------------------------------------------------------------- nglis = nc*nr if (kpds(5)==59 .or. kpds(5)==214) then ip=3 ipopt(1)=-1 ipopt(2)=-1 km=1 ibi=1 else ip=0 do i=1,20 ipopt(i)=0 enddo km=1 ibi=1 endif !----------------------------------------------------------------------- ! Initialize output bitmap. Important for soil moisture and temp. !----------------------------------------------------------------------- lo = .true. !----------------------------------------------------------------------- ! Interpolate to LIS grid !----------------------------------------------------------------------- if (kpds(5)==59 .or. kpds(5)==214) then call polates3(lis_gds,ibi,lb,f,ibo,lo,lis1d,mi, & rlat3,rlon3,w113,w123,w213,w223,n113,n123,n213,n223,iret) else call polates0 (lis_gds,ibi,lb,f,ibo,lo,lis1d,mi,& rlat0, rlon0,w110,w120,w210,w220,n110,n120,n210,n220,iret) endif !----------------------------------------------------------------------- ! Create 2D array for main program. Also define a "soil" mask ! due to different geography between GDAS & LDAS. For LDAS land ! points not included in GDAS geography dataset only. !----------------------------------------------------------------------- count = 0 do j = 1, nr do i = 1, nc varfield(i,j) = lis1d(i+count) geogmask(i,j) = lo(i+count) enddo count = count + nc enddo !----------------------------------------------------------------------- ! Save air tempertaure interpolated field for later use in ! initialization of soil temp where geography differs ! between GDAS and LDAS !----------------------------------------------------------------------- if (kpds(5) .eq. 11 .and. kpds(6) .eq. 105) then do i = 1, nc do j = 1, nr geogtemp(i,j) = varfield(i,j) enddo enddo endif !EOC end subroutine interp_gdas