PROGRAM wrfstatic,41
!
!##################################################################
!##################################################################
!###### ######
!###### PROGRAM WRFSTATIC ######
!###### ######
!###### Developed by ######
!###### Center for Analysis and Prediction of Storms ######
!###### University of Oklahoma ######
!###### ######
!##################################################################
!##################################################################
!
!
!-----------------------------------------------------------------------
!
!
! PURPOSE:
!
! Read static file sets and produce static file for ARPS2WRF & WRF.
!
! NOTE:
! o It will replace both grid_gen.exe and staticpost in WRFSI.
! o This program shares the same namelist input file as ARPS2WRF
! for convenience with an extra namelist variable for static
! dataset directory in input/arps2wrf.input.
!
! Required Static datasets (the same as WRFSI requirements):
!
! 1. albedo_ncep
! 2. greenfrac
! 3. islope
! 4. landuse_30s
! 5. maxsnowalb
! 6. soiltemp_1deg
! 7. soiltype_bot_30s
! 8. soiltype_top_30s
! *. topo_30s (not required for ARPS2WRF, because terrain will
! be provided by arpstrn/arpstern.)
!
!-----------------------------------------------------------------------
!
! AUTHOR: Yunheng Wang (09/28/2005)
! Created initially based on WRFSIV2.1.
!
! MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
USE wrf_metadata
! WRF constants and metadata
!-----------------------------------------------------------------------
!
! Variable Declarations:
!
!-----------------------------------------------------------------------
!
IMPLICIT NONE
!----------------------------------------------------------------------
!
! ARPS grid variables
!
!---------------------------------------------------------------------
INTEGER, PARAMETER :: nhisfile_max = 100
INTEGER, PARAMETER :: max_vertical_levels = 100
CHARACTER(LEN=256) :: grdbasfn
CHARACTER(LEN=256) :: bdybasfn(nhisfile_max)
CHARACTER(LEN=256) :: hisfile(nhisfile_max)
INTEGER :: hinfmt,nhisfile,lengbf
CHARACTER(LEN=80) :: execname
!
!-----------------------------------------------------------------------
!
! ARPS include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc'
!
!-----------------------------------------------------------------------
!
! Variables for split files
!
!-----------------------------------------------------------------------
INTEGER :: nprocx_in, nprocy_in
INTEGER :: nxsm, nysm
INTEGER, PARAMETER :: fzone_arps = 3
INTEGER, PARAMETER :: fzone_wrf = 1
!-----------------------------------------------------------------------
!
! ARPS variables to be read in:
!
!-----------------------------------------------------------------------
!
INTEGER :: nx,ny,nz ! Grid dimensions for ARPS.
INTEGER :: nzsoil ! Soil levels
INTEGER :: nstyps ! Maximum number of soil types.
INTEGER :: mapproj
REAL :: sclfct
REAL :: trulat1, trulat2, trulon
REAL :: ctrlat, ctrlon
REAL :: dx, dy
REAL :: time
!-----------------------------------------------------------------------
!
! Namelist definitions for ARPS2WRF.input
!
! sfcdt Specify surface characteristics
! bdyspc Obtain boundary input files (ARPS format)
! wrf_grid Define WRF horizontal and vertical grid
! interp_options Choose interpolation scheme
! wrf_opts WRF options from namelist.input
! output Ouput options
!
!-----------------------------------------------------------------------
!
INTEGER :: nx_wrf ! = nx-2 if the same grid as ARPS
INTEGER :: ny_wrf ! = ny-2 if the same grid as ARPS
INTEGER :: nz_wrf
INTEGER :: nzsoil_wrf
REAL :: dx_wrfscl ! grid length times map scale
REAL :: dy_wrfscl
REAL :: lattru_wrf(2) ! array of true latitude of WRF map projection
REAL :: lontru_wrf ! true longitude of WRF map projection
! = trulon_wrf
REAL :: sclf_wrf ! map scale factor (usually = 1.0)
! Namelist variable declaration
CHARACTER(LEN=5) :: sfcinitopt ! either "ARPS" or "WRFSI"
CHARACTER(LEN=256):: static_dir
CHARACTER(LEN=19) :: start_date
REAL :: wvln, silwt
INTEGER :: use_arps_grid ! Use ARPS horizontal grid as WRF grid
INTEGER :: mapproj_wrf ! Type of map projection in WRF model grid
! modproj = 1 Polar Stereographic
! projection.
! = 2 Mercator projection.
! = 3 Lambert projection.
REAL :: sclfct_wrf ! Map scale factor.
! Distance on map, between two latitudes
! trulat1 and trulat2,
! is = (Distance on earth)*sclfct.
! For ARPS model runs,
! generally this is 1.0
REAL :: trulat1_wrf, trulat2_wrf, trulon_wrf
! 1st, 2nd real true latitude and true longitude
! of WRF map projection
REAL :: ctrlat_wrf ! Center latitude of WRF model domain (deg. N)
REAL :: ctrlon_wrf ! Center longitude of WRF model domain (deg. E)
REAL :: dx_wrf ! WRF Grid spacing in x-direction
REAL :: dy_wrf ! WRF Grid spacing in y-direction
INTEGER :: spec_bdy_width
INTEGER :: mp_physics
INTEGER :: io_form
INTEGER :: idummy
REAL :: rdummy
REAL :: rdumary(max_vertical_levels)
CHARACTER(LEN=256) :: cdummy
!-----------------------------------------------------------------------
!
! WRF grid related variables
!
!-----------------------------------------------------------------------
REAL, ALLOCATABLE :: x_wrf(:) ! X coordinate of WRF U points
REAL, ALLOCATABLE :: y_wrf(:) ! Y coordinate of WRF V points
REAL, ALLOCATABLE :: xs_wrf(:) ! X coordinate of WRF mass points
REAL, ALLOCATABLE :: ys_wrf(:) ! Y coordinate of WRF mass points
REAL, ALLOCATABLE :: lat_wrf(:,:,:) ! WRF grid point latitudes
REAL, ALLOCATABLE :: lon_wrf(:,:,:) ! WRF grid point lontitudes
REAL, ALLOCATABLE :: msft_wrf(:,:) ! Map scale factor on mass grid
REAL, ALLOCATABLE :: msfu_wrf(:,:) ! Map scale factor on u grid
REAL, ALLOCATABLE :: msfv_wrf(:,:) ! Map scale factor on v grid
REAL, ALLOCATABLE :: cosalpha(:,:)
REAL, ALLOCATABLE :: sinalpha(:,:)
REAL, ALLOCATABLE :: f(:,:)
REAL, ALLOCATABLE :: e(:,:)
REAL, ALLOCATABLE :: hgt(:,:)
REAL, ALLOCATABLE :: topostdv(:,:)
REAL, ALLOCATABLE :: toposlpx(:,:)
REAL, ALLOCATABLE :: toposlpy(:,:)
REAL, ALLOCATABLE :: landuse(:,:)
REAL, ALLOCATABLE :: landmask(:,:)
REAL, ALLOCATABLE :: landusef(:,:,:)
REAL, ALLOCATABLE :: soilctop(:,:,:)
REAL, ALLOCATABLE :: soilcbot(:,:,:)
REAL, ALLOCATABLE :: green12m(:,:,:)
REAL, ALLOCATABLE :: albdo12m(:,:,:)
REAL, ALLOCATABLE :: shdmax(:,:)
REAL, ALLOCATABLE :: shdmin(:,:)
REAL, ALLOCATABLE :: tmn(:,:)
REAL, ALLOCATABLE :: slopecat(:,:)
REAL, ALLOCATABLE :: snoalb(:,:)
REAL :: xsub0, ysub0
REAL :: lat_ll(4), lat_lr(4), lat_ul(4), lat_ur(4)
REAL :: lon_ll(4), lon_lr(4), lon_ul(4), lon_ur(4)
!-----------------------------------------------------------------------
!
! Temporary working arrays
!
!-----------------------------------------------------------------------
REAL, ALLOCATABLE :: tem2d1(:,:), tem2d2(:,:), tem2d3(:,:)
REAL, ALLOCATABLE :: tem3d(:,:,:)
!-----------------------------------------------------------------------
!
! Misc local variables
!
!-----------------------------------------------------------------------
!
INTEGER :: i,j,k,n
INTEGER :: istatus, lenstr, ireturn
INTEGER :: abstime
INTEGER :: iter
LOGICAL :: fexist,dem_data
REAL :: latnot(2),ctrx, ctry, swx, swy
CHARACTER(LEN=200) :: filename, tmpstr
CHARACTER(LEN=24) :: times_str
INTEGER :: nchr, ncid, nchout
REAL, PARAMETER :: eps = 1.0E-5
CHARACTER(LEN=16), PARAMETER :: setnames(9) = (/ 'topo_30s ', &
'landuse_30s ', 'islope ', 'greenfrac ', &
'maxsnowalb ', 'soiltemp_1deg ', 'soiltype_bot_30s', &
'soiltype_top_30s', 'albedo_ncep ' /)
CHARACTER(LEN=2), PARAMETER :: setindx(9) = (/ '/U', &
'/V', '/I', '/G', &
'/M', '/T', '/O', &
'/O', '/A' /)
INTEGER, PARAMETER :: TOPO = 1, LUSE = 2, ISLOPE = 3, &
GFRAC= 4, SNOWALB = 5, SOILTMP= 6, &
SOILTB=7, SOILTT = 8, ALBEDO = 9
CHARACTER(LEN=200) :: path_to_soiltemp_1deg
COMMON /lapscmn/ path_to_soiltemp_1deg
REAL :: projrot_latlon
REAL :: projrot_deg
REAL :: bilinear_interp
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
! Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
WRITE(6,'(10(/5x,a),/)') &
'###############################################################',&
'###############################################################',&
'#### ####',&
'#### Welcome to WRFSTATIC ####',&
'#### ####',&
'#### Program that reads in static datasets and generates ####',&
'#### WRF static file for ARPS2WRF or WRF. ####',&
'#### ####',&
'###############################################################',&
'###############################################################'
!-----------------------------------------------------------------------
!
! First, read namelist input
!
!-----------------------------------------------------------------------
CALL readnamelist(1,hinfmt,bdybasfn,hisfile,nhisfile, &
nx,ny,idummy,idummy,idummy,nprocx_in,nprocy_in,idummy,idummy, &
use_arps_grid,nx_wrf,ny_wrf,idummy,rdumary,rdummy, &
mapproj_wrf,sclfct_wrf,lattru_wrf,lontru_wrf, &
ctrlat_wrf,ctrlon_wrf,dx_wrf,dy_wrf,rdummy,rdummy, &
sfcinitopt,static_dir,start_date,silwt,wvln, &
idummy,idummy,rdummy,spec_bdy_width, &
idummy,rdummy,rdummy,rdummy,mp_physics,idummy, &
idummy,idummy,idummy,idummy, &
idummy,idummy,idummy,idummy,idummy,io_form, &
idummy,cdummy,istatus)
trulat1_wrf = lattru_wrf(1)
trulat2_wrf = lattru_wrf(2)
trulon_wrf = lontru_wrf
WRITE(grdbasfn,'(a)') bdybasfn(1)
nz_wrf = 1
nzsoil_wrf = 1
nproc_x = 1
nproc_y = 1
loc_x = 1
loc_y = 1
myproc = 0
readsplit = 0
execname = 'WRFSTATIC'
!--------------------------------------------------------------------
!
! Begin to set WRF horizontal grid
!
!--------------------------------------------------------------------
IF (use_arps_grid == 1) THEN
lengbf = LEN_TRIM(grdbasfn)
IF (nprocx_in > 1 .OR. nprocy_in > 1) THEN
WRITE(grdbasfn(lengbf+1:lengbf+5),'(a,i2.2,i2.2)') '_',1,1
lengbf = lengbf + 5
END IF
WRITE(6,'(1x,/2a)') ' Reading file: ',grdbasfn(1:lengbf)
CALL get_arps_dims_atts(nchr,hinfmt,grdbasfn,lengbf, &
nx,ny,nz,nzsoil,nstyps, &
year,month,day,hour,minute,second,time, &
mapproj, sclfct,trulat1,trulat2,trulon, &
ctrlat,ctrlon,dx,dy,ireturn)
IF( ireturn /= 0 ) CALL arpsstop('dtaread errors.',1)
nxsm = nx
nysm = ny
IF (nprocx_in > 1 .OR. nprocy_in > 1) THEN
nx = (nxsm-fzone_arps)*nprocx_in + fzone_arps
ny = (nysm-fzone_arps)*nprocy_in + fzone_arps
END IF
nx_wrf = nx - 2
ny_wrf = ny - 2
mapproj_wrf = mapproj
sclfct_wrf = sclfct
trulat1_wrf = trulat1
trulat2_wrf = trulat2
trulon_wrf = trulon
ctrlat_wrf = ctrlat
ctrlon_wrf = ctrlon
dx_wrf = dx
dy_wrf = dy
WRITE(6,'(a38)') ' ARPS grid WRF grid'
WRITE(6,'(a38)') ' ========== =========='
WRITE(6,'(a,I10,4x,I10)') ' nx = ',nx,nx_wrf
WRITE(6,'(a,I10,4x,I10)') ' ny = ',ny,ny_wrf
WRITE(6,'(a,I10,4x,I10)') ' mapproj = ',mapproj,mapproj_wrf
WRITE(6,'(a,F10.2,4x,F10.2)') ' sclfct = ',sclfct,sclfct_wrf
WRITE(6,'(a,F10.2,4x,F10.2)') ' trulat1 = ',trulat1,trulat1_wrf
WRITE(6,'(a,F10.2,4x,F10.2)') ' trulat2 = ',trulat2,trulat2_wrf
WRITE(6,'(a,F10.2,4x,F10.2)') ' trulon = ',trulon,trulon_wrf
WRITE(6,'(a,F10.2,4x,F10.2)') ' ctrlat = ',ctrlat,ctrlat_wrf
WRITE(6,'(a,F10.2,4x,F10.2)') ' ctrlon = ',ctrlon,ctrlon_wrf
WRITE(6,'(a,F10.0,4x,F10.0)') ' dx = ',dx,dx_wrf
WRITE(6,'(a,F10.0,4x,F10.0)') ' dy = ',dy,dy_wrf
CALL ctim2abss( year,month,day,hour,minute,second, abstime )
abstime = abstime + INT(time)
CALL abss2ctim( abstime, year, month, day, hour, minute, second )
WRITE(start_date,'(I4.4,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2)') &
year,'-',month,'-',day,'_',hour,':',minute,':',second
END IF
times_str = start_date
!---------------------------------------------------------------------------
!
! Allocate WRF arrays
!
!---------------------------------------------------------------------------
ALLOCATE(x_wrf (nx_wrf), STAT = istatus)
ALLOCATE(y_wrf (ny_wrf), STAT = istatus)
ALLOCATE(xs_wrf(nx_wrf), STAT = istatus)
ALLOCATE(ys_wrf(ny_wrf), STAT = istatus)
ALLOCATE(lat_wrf(nx_wrf,ny_wrf,3), STAT = istatus) ! mass grid
ALLOCATE(lon_wrf(nx_wrf,ny_wrf,3), STAT = istatus)
ALLOCATE(msft_wrf(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(msfu_wrf(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(msfv_wrf(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(f(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(e(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(cosalpha(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(sinalpha(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(hgt(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(topostdv(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(toposlpx(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(toposlpy(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(landuse (nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(landmask(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(slopecat(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(tmn(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(snoalb(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(shdmax(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(shdmin(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(landusef(nx_wrf,ny_wrf,LanduseCategories), STAT = istatus)
ALLOCATE(soilctop(nx_wrf,ny_wrf,SoilCategories), STAT = istatus)
ALLOCATE(soilcbot(nx_wrf,ny_wrf,SoilCategories), STAT = istatus)
ALLOCATE(green12m(nx_wrf,ny_wrf,12), STAT = istatus)
ALLOCATE(albdo12m(nx_wrf,ny_wrf,12), STAT = istatus)
ALLOCATE(tem2d1(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(tem2d2(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(tem2d3(nx_wrf,ny_wrf), STAT = istatus)
ALLOCATE(tem3d (nx_wrf,ny_wrf,LanduseCategories), STAT = istatus)
msft_wrf = 0.0
msfu_wrf = 0.0
msfv_wrf = 0.0
!
!-----------------------------------------------------------------------
!
! Establish WRF map projection
!
!-----------------------------------------------------------------------
!
lattru_wrf(1) = trulat1_wrf
lattru_wrf(2) = trulat2_wrf
lontru_wrf = trulon_wrf
sclf_wrf = 1.0/sclfct_wrf
dx_wrfscl = dx_wrf*sclf_wrf
dy_wrfscl = dy_wrf*sclf_wrf
CALL setmapr(mapproj_wrf,sclf_wrf,lattru_wrf,lontru_wrf)
CALL lltoxy( 1,1, ctrlat_wrf,ctrlon_wrf, ctrx, ctry )
swx = ctrx - 0.5*(nx_wrf-1) * dx_wrfscl
swy = ctry - 0.5*(ny_wrf-1) * dy_wrfscl
CALL setorig( 1, swx, swy)
xsub0 = dx_wrf * (nx_wrf-fzone_wrf) * (loc_x-1)
ysub0 = dy_wrf * (ny_wrf-fzone_wrf) * (loc_y-1)
DO i=1,nx_wrf
x_wrf(i)= sclf_wrf*xsub0 + (i-1)*dx_wrfscl
END DO
DO j=1,ny_wrf
y_wrf(j)= sclf_wrf*ysub0 + (j-1)*dy_wrfscl
END DO
DO i=1,nx_wrf-1
xs_wrf(i)=0.5*(x_wrf(i)+x_wrf(i+1))
END DO
xs_wrf(nx_wrf)=2.*xs_wrf(nx_wrf-1)-xs_wrf(nx_wrf-2)
DO j=1,ny_wrf-1
ys_wrf(j)=0.5*(y_wrf(j)+y_wrf(j+1))
END DO
ys_wrf(ny_wrf)=2.*ys_wrf(ny_wrf-1)-ys_wrf(ny_wrf-2)
!-----------------------------------------------------------------------
!
! Find latitude and longitude of WRF grid.
!
!-----------------------------------------------------------------------
CALL xytoll(nx_wrf,ny_wrf,xs_wrf,ys_wrf, &
lat_wrf(:,:,1),lon_wrf(:,:,1)) ! T points
CALL xytoll(nx_wrf,ny_wrf,x_wrf,ys_wrf, &
lat_wrf(:,:,2),lon_wrf(:,:,2)) ! U points
CALL xytoll(nx_wrf,ny_wrf,xs_wrf,y_wrf, &
lat_wrf(:,:,3),lon_wrf(:,:,3)) ! V points
!-----------------------------------------------------------------------
!
! Find Map scale factor
!
!-----------------------------------------------------------------------
CALL lattomf(nx_wrf,ny_wrf,lat_wrf(:,:,1),msft_wrf) !mass points
CALL lattomf(nx_wrf,ny_wrf,lat_wrf(:,:,2),msfu_wrf) !U points
CALL lattomf(nx_wrf,ny_wrf,lat_wrf(:,:,3),msfv_wrf) !V points
!-----------------------------------------------------------------------
!
! Latitude and Longitude variables at corns
!
!-----------------------------------------------------------------------
!
! Latitude, 1 -- T point, 2 -- U point, 3 -- V point, 4 -- massless point
!
lat_ll(1) = lat_wrf( 1, 1,1)
lat_ul(1) = lat_wrf( 1,ny_wrf-1,1)
lat_ur(1) = lat_wrf(nx_wrf-1,ny_wrf-1,1)
lat_lr(1) = lat_wrf(nx_wrf-1, 1,1)
lon_ll(1) = lon_wrf( 1, 1,1)
lon_ul(1) = lon_wrf( 1,ny_wrf-1,1)
lon_ur(1) = lon_wrf(nx_wrf-1,ny_wrf-1,1)
lon_lr(1) = lon_wrf(nx_wrf-1, 1,1)
lat_ll(2) = lat_wrf( 1, 1,2)
lat_ul(2) = lat_wrf( 1,ny_wrf-1,2)
lat_ur(2) = lat_wrf(nx_wrf, ny_wrf-1,2)
lat_lr(2) = lat_wrf(nx_wrf, 1,2)
lon_ll(2) = lon_wrf( 1, 1,2)
lon_ul(2) = lon_wrf( 1,ny_wrf-1,2)
lon_ur(2) = lon_wrf(nx_wrf, ny_wrf-1,2)
lon_lr(2) = lon_wrf(nx_wrf, 1,2)
lat_ll(3) = lat_wrf( 1, 1,3)
lat_ul(3) = lat_wrf( 1, ny_wrf,3)
lat_ur(3) = lat_wrf(nx_wrf-1, ny_wrf,3)
lat_lr(3) = lat_wrf(nx_wrf-1, 1,3)
lon_ll(3) = lon_wrf( 1, 1,3)
lon_ul(3) = lon_wrf( 1, ny_wrf,3)
lon_ur(3) = lon_wrf(nx_wrf-1, ny_wrf,3)
lon_lr(3) = lon_wrf(nx_wrf-1, 1,3)
CALL xytoll(1,1,x_wrf(1), y_wrf(1), lat_ll(4),lon_ll(4))
CALL xytoll(1,1,x_wrf(1), y_wrf(ny_wrf),lat_ul(4),lon_ul(4))
CALL xytoll(1,1,x_wrf(nx_wrf),y_wrf(ny_wrf),lat_ur(4),lon_ur(4))
CALL xytoll(1,1,x_wrf(nx_wrf),y_wrf(1), lat_lr(4),lon_lr(4))
DO j=1,ny_wrf-1
DO i=1,nx_wrf-1
f(i,j) =2*omega_ear*SIN(lat_wrf(i,j,1)*d2rfactor)
e(i,j) =2*omega_ear*COS(lat_wrf(i,j,1)*d2rfactor)
END DO
END DO
DO j=1,ny_wrf-1
DO i=1,nx_wrf-1
projrot_deg = projrot_latlon(mapproj_wrf,trulat1_wrf, &
trulat2_wrf,trulon_wrf,ctrlat_wrf,ctrlon_wrf,&
lat_wrf(i,j,1),lon_wrf(i,j,1),istatus)
IF(istatus == 1)THEN
sinalpha(i,j) = SIN(d2rfactor*projrot_deg)
cosalpha(i,j) = COS(d2rfactor*projrot_deg)
ELSE
CALL wrf_debug(1,'WARNING: Problem in projrot_latlon.')
END IF
END DO
END DO
!---------------------- Terrain USGS 30 sec --------------------------
! type = U
! symbol = TOPO
lenstr = LEN_TRIM(static_dir)
WRITE(tmpstr,'(3a)') static_dir(1:lenstr), &
TRIM(setnames(TOPO)),setindx(TOPO)
CALL rdgeodat(nx_wrf,ny_wrf,x_wrf,y_wrf,dx_wrf, &
dy_wrf,TRIM(tmpstr),wvln,silwt,1, &
tem2d1,tem3d,tem2d2,tem2d3,dem_data,istatus)
hgt(:,:)=0.0
DO j = 2,ny_wrf
DO i = 2,nx_wrf
hgt(i-1,j-1) = bilinear_interp(i,j,nx_wrf,ny_wrf,tem2d1)
END DO
END DO
WHERE( (hgt(:,:)< 0.01) .AND. (hgt(:,:)> -0.01) ) hgt(:,:) = 0.0
topostdv(:,:) = tem3d (:,:,1) ! Actually, they are over non-stagger points, why?
toposlpx(:,:) = tem2d2(:,:)
toposlpy(:,:) = tem2d3(:,:)
! ------------------- 30s Landuse -------------------------
! type = V
! symbol = LUSE
wvln = 2.0 ! they will be fixed from now on
silwt = 0.0
WRITE(tmpstr,'(3a)') static_dir(1:lenstr), &
TRIM(setnames(LUSE)),setindx(LUSE)
CALL rdgeodat(nx_wrf,ny_wrf,xs_wrf,ys_wrf,dx_wrf, &
dy_wrf,TRIM(tmpstr),wvln,silwt,LanduseCategories, &
tem2d1,tem3d,tem2d2,tem2d3,dem_data,istatus)
landuse(:,:) = tem2d1(:,:)
landmask(:,:) = 1.
where(landuse(:,:) == ISWATER) landmask(:,:) = 0.
DO n = 1,LanduseCategories
landusef(:,:,n) = tem3d(:,:,n)
END DO
! -------------------------------------------------------------------
! potential fix of fictitous islands for certain resolution domains.
! Story here is that west coast terrain can be "funky" possibly due
! to steepness and method to compute avg terrain.
tem2d1(:,:) = 1.- tem3d(:,:,ISWATER)
! Filter land fraction with 2dx,iter
iter = MIN(10,INT(8000./dx_wrf))
DO n = 1,iter
CALL filter_2dx(tem2d1,nx_wrf,ny_wrf,1, 0.5)
CALL filter_2dx(tem2d1,nx_wrf,ny_wrf,1,-0.5)
END DO
WHERE(tem2d1(:,:) <= 0.1 .AND. hgt(:,:) < 5.0) hgt(:,:) = 0.0
! -------------- Top Layer Soil Type -------------------------
! type = O
!
WRITE(tmpstr,'(3a)') static_dir(1:lenstr), &
TRIM(setnames(SOILTT)),setindx(SOILTT)
CALL rdgeodat(nx_wrf,ny_wrf,xs_wrf,ys_wrf,dx_wrf, &
dy_wrf,TRIM(tmpstr),wvln,silwt,SoilCategories, &
tem2d1,tem3d,tem2d2,tem2d3,dem_data,ireturn)
! make water points = 0 for adjust_geog. later well put it back to original
WHERE(tem2d1 == 14) tem2d1 = 0.0
CALL adjust_geog(nx_wrf-1,ny_wrf-1,1,'soiltype',ireturn, &
lat_wrf (1:nx_wrf-1,1:ny_wrf-1,1), hgt(1:nx_wrf-1,1:ny_wrf-1), &
landmask(1:nx_wrf-1,1:ny_wrf-1),tem2d1(1:nx_wrf-1,1:ny_wrf-1), &
istatus)
WHERE(tem2d1 == 0.0) tem2d1 = 14.
DO n = 1,SoilCategories
soilctop(:,:,n) = tem3d(:,:,n)
END DO
! -------------- Bottom Layer Soil Type -------------------------
! type = O
!
WRITE(tmpstr,'(3a)') static_dir(1:lenstr), &
TRIM(setnames(SOILTB)),setindx(SOILTB)
CALL rdgeodat(nx_wrf,ny_wrf,xs_wrf,ys_wrf,dx_wrf, &
dy_wrf,TRIM(tmpstr),wvln,silwt,SoilCategories, &
tem2d1,tem3d,tem2d2,tem2d3,dem_data,ireturn)
! make water points = 0 for adjust_geog. later well put it back to original
WHERE(tem2d1 == 14) tem2d1 = 0.0
CALL adjust_geog(nx_wrf-1,ny_wrf-1,1,'soiltype',ireturn, &
lat_wrf (1:nx_wrf-1,1:ny_wrf-1,1), hgt(1:nx_wrf-1,1:ny_wrf-1), &
landmask(1:nx_wrf-1,1:ny_wrf-1),tem2d1(1:nx_wrf-1,1:ny_wrf-1), &
istatus)
WHERE(tem2d1 == 0.0) tem2d1 = 14.
DO n = 1,SoilCategories
soilcbot(:,:,n) = tem3d(:,:,n)
END DO
! ----------------- Greenness Fraction --------------------
! type = G
!
WRITE(tmpstr,'(3a)') static_dir(1:lenstr), &
TRIM(setnames(GFRAC)),setindx(GFRAC)
tem3d(:,:,:) = 0.0
CALL proc_geodat(nx_wrf,ny_wrf,12,TRIM(tmpstr), &
lat_wrf(:,:,1),lon_wrf(:,:,1),landmask,tem3d,istatus)
DO n = 1,12
green12m(:,:,n) = tem3d(:,:,n)
END DO
! annual max/min greenness fraction in domain
! --------------------------------------------
print*,' compute max/min greenness frac at grid points'
DO j=1,ny_wrf
DO i=1,nx_wrf
shdmax(i,j) = MAXVAL(green12m(i,j,:))
shdmin(i,j) = MINVAL(green12m(i,j,:))
END DO
END DO
! --------------- Deep Soil Temperature -------------------
! type = T
!
WRITE(tmpstr,'(3a)') static_dir(1:lenstr), &
TRIM(setnames(SOILTMP)),setindx(SOILTMP)
path_to_soiltemp_1deg = tmpstr
CALL proc_geodat(nx_wrf,ny_wrf,1,TRIM(tmpstr), &
lat_wrf(:,:,1),lon_wrf(:,:,1),landmask, &
tmn,ireturn)
CALL adjust_geog(nx_wrf-1,ny_wrf-1,1,'soiltemp',ireturn, &
lat_wrf(1:nx_wrf-1,1:ny_wrf-1,1),hgt(1:nx_wrf-1,1:ny_wrf-1), &
landmask(1:nx_wrf-1,1:ny_wrf-1),tmn(1:nx_wrf-1,1:ny_wrf-1), &
istatus)
! -------------- Terrain slope index categories -------------------------
! type = I
!
WRITE(tmpstr,'(3a)') static_dir(1:lenstr), &
TRIM(setnames(ISLOPE)),setindx(ISLOPE)
CALL rdgeodat(nx_wrf,ny_wrf,xs_wrf,ys_wrf,dx_wrf, &
dy_wrf,TRIM(tmpstr),wvln,silwt,9, &
tem2d1,tem3d,tem2d2,tem2d3,dem_data,ireturn)
slopecat(:,:) = tem2d1(:,:)
! put the categories back to the original raw data. if it is a land
! point but islope indicates water, force islope = 1.
WHERE(slopecat .eq. 8) slopecat = 13.0
WHERE(slopecat .eq. 9) slopecat = 0.0
CALL adjust_geog(nx_wrf-1,ny_wrf-1,1,'islope',ireturn, &
lat_wrf(1:nx_wrf-1,1:ny_wrf-1,1),hgt(1:nx_wrf-1,1:ny_wrf-1), &
landmask(1:nx_wrf-1,1:ny_wrf-1),slopecat(1:nx_wrf-1,1:ny_wrf-1), &
istatus)
! force land points to have the correct (default) islope value.
WHERE(slopecat(:,:) == 0.0 .AND. landmask(:,:) == 1.0) &
slopecat(:,:) = 1.0
! force water points to have the correct category for islope
WHERE(landmask(:,:) == 0.0) slopecat(:,:)=0.0
! ---------------Monthly Albedo Climatology--------------------------
! type = A
!
WRITE(tmpstr,'(3a)') static_dir(1:lenstr), &
TRIM(setnames(ALBEDO)),setindx(ALBEDO)
CALL proc_geodat(nx_wrf,ny_wrf,12,TRIM(tmpstr), &
lat_wrf(:,:,1),lon_wrf(:,:,1),landmask(:,:),tem3d,istatus)
DO n=1,12
albdo12m(:,:,n) = tem3d(:,:,n)
WHERE(landmask(:,:) == 0.0) albdo12m(:,:,n) = 0.08
END DO
! ---------------- Max Snow Albedo ------------------
! type = M
!
WRITE(tmpstr,'(3a)') static_dir(1:lenstr), &
TRIM(setnames(SNOWALB)),setindx(SNOWALB)
CALL proc_geodat(nx_wrf,ny_wrf,1,TRIM(tmpstr), &
lat_wrf(:,:,1),lon_wrf(:,:,1),landmask(:,:), &
tem2d1,istatus)
snoalb(:,:) = tem2d1
!
! force max albedo = 0.08 over water. force max albedo = 0.7 over ice
!
WHERE(landmask(:,:) == 0.0) snoalb = 0.08
WHERE(landuse(:,:) == ISICE)snoalb = 0.7
! ---------------------------------------------------------------------------------
! Lets compare the grids to landmask to assess their consistency (or lack thereof)
! ---------------------------------------------------------------------------------
! CALL gridcompare(nx_wrf,ny_wrf,i,xxx,landmask)
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!
! OUTPUT of WRF static begin
!
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filename = sfcdtfl
lenstr = LEN_TRIM(dirname)
IF(lenstr > 0) filename = dirname(1:lenstr) // TRIM(filename)
PRINT *, ' Output file name is ',TRIM(filename)
CALL open_output_file(filename,'STATIC',io_form,nx_wrf,ny_wrf, &
nz_wrf, nzsoil_wrf,spec_bdy_width,mp_physics,ncid)
!-----------------------------------------------------------------------
!
! Initialize and write global attributes
!
!-----------------------------------------------------------------------
CALL write_static_attribute(ncid,io_form,start_date, &
nx_wrf,ny_wrf,dx_wrf,dy_wrf, &
mapproj_wrf,trulat1_wrf,trulat2_wrf,trulon_wrf,&
ctrlat_wrf,ctrlon_wrf, &
lat_ll,lat_ul,lat_ur,lat_lr, &
lon_ll,lon_ul,lon_ur,lon_lr, istatus)
CALL write_times_str(ncid,io_form,'Times', &
times_str(1:19),times_str(1:19),1)
CALL write_wrf_static(ncid,io_form,times_str(1:19), &
nx_wrf,ny_wrf,LanduseCategories,SoilCategories, &
lat_wrf(:,:,1),lon_wrf(:,:,1),f,e, &
msft_wrf,msfu_wrf,msfv_wrf,cosalpha,sinalpha, &
hgt,topostdv,toposlpx,toposlpy,landmask,landusef, &
soilctop,soilcbot,tmn,slopecat,snoalb, &
shdmax,shdmin,green12m,albdo12m, istatus)
CALL close_output_file(ncid,io_form) ! close input file
CALL io_shutdown(io_form)
!-----------------------------------------------------------------------
!
! Generate a ready file if needed
!
!-----------------------------------------------------------------------
IF ( readyfl == 1 ) THEN
CALL getunit( nchout )
WRITE (filename,'(2a)') trim(sfcdtfl),"_ready"
WRITE(filename,'(a)') TRIM(dirname) // TRIM(filename)
OPEN (UNIT=nchout,FILE=trim(filename))
WRITE (nchout,'(a)') trim(sfcdtfl)
CLOSE (nchout)
CALL retunit (nchout )
END IF
CALL arpsstop(' ==== WRFSTATIC terminated normally. ====',0)
END PROGRAM wrfstatic