PROGRAM wrfstatic,48
!
!##################################################################
!##################################################################
!###### ######
!###### 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 :: nmax_domains = 100
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 :: nhisfile,lengbf
CHARACTER(LEN=256), DIMENSION(nmax_domains) :: adashisfn
CHARACTER(LEN=256), DIMENSION(nmax_domains) :: adasbasfn
INTEGER, DIMENSION(nmax_domains) :: hinfmt
LOGICAL, DIMENSION(nmax_domains) :: input_from_file
INTEGER :: max_dom
INTEGER :: finexist(nhisfile_max)
CHARACTER(LEN=80) :: execname
!
!-----------------------------------------------------------------------
!
! ARPS include files:
!
!-----------------------------------------------------------------------
!
INCLUDE 'globcst.inc'
INCLUDE 'mp.inc'
!
!-----------------------------------------------------------------------
!
! Variables for split files
!
!-----------------------------------------------------------------------
INTEGER :: nprocx_in, nprocy_in
!-----------------------------------------------------------------------
!
! ARPS variables to be read in:
!
!-----------------------------------------------------------------------
!
INTEGER, DIMENSION(nmax_domains) :: nxnm,nynm,nznm
INTEGER, DIMENSION(nmax_domains) :: nzsoilnm
INTEGER, DIMENSION(nmax_domains) :: nstypsnm
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, DIMENSION(nmax_domains) :: nx_wrfnm
INTEGER, DIMENSION(nmax_domains) :: ny_wrfnm
INTEGER, DIMENSION(nmax_domains) :: i_parent_start
INTEGER, DIMENSION(nmax_domains) :: j_parent_start
INTEGER, DIMENSION(nmax_domains) :: parent_id
INTEGER :: i_parent_end ! they are use as temporary variable
INTEGER :: j_parent_end ! So do not have to be saved
INTEGER :: parent_grid_ratio
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), DIMENSION(nmax_domains) :: sfcinitopt ! either "ARPS" or "WRFSI"
CHARACTER(LEN=256) :: geogdir
CHARACTER(LEN=19), DIMENSION(nmax_domains) :: start_date
CHARACTER(LEN=256),DIMENSION(nmax_domains) :: sfcdtfn
REAL :: wvln, silwt
INTEGER, DIMENSION(nmax_domains) :: 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
REAL, DIMENSION(nmax_domains) :: dx_wrfnm
REAL, DIMENSION(nmax_domains) :: dy_wrfnm
INTEGER :: spec_bdy_width
INTEGER, DIMENSION(nmax_domains) :: mp_physics
INTEGER :: io_form
INTEGER, DIMENSION(nmax_domains) :: idumar
REAL, DIMENSION(nmax_domains) :: rdumar
INTEGER :: idummy
REAL :: rdummy, rdumary(max_vertical_levels)
LOGICAL :: ldummy
CHARACTER(LEN=256) :: cdummy
INTEGER :: mapproj_moad
REAL :: ctrlat_moad
REAL :: trulat1_moad, trulat2_moad, trulon_moad
!-----------------------------------------------------------------------
!
! 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, DIMENSION(nmax_domains) :: 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
INTEGER :: domid
LOGICAL :: fexist,dem_data
REAL :: latnot(2),ctrx, ctry, swx, swy
CHARACTER(LEN=256) :: 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=256) :: 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,max_dom,input_from_file, &
hinfmt,adashisfn,adasbasfn,bdybasfn,hisfile,nhisfile,finexist, &
nxnm,nynm,idumar,idumar,idumar,nprocx_in,nprocy_in,idummy,idummy, &
use_arps_grid,nx_wrfnm,ny_wrfnm,idummy,rdumary,rdumar, &
i_parent_start,j_parent_start,parent_id, &
mapproj_wrf,sclfct_wrf,lattru_wrf,lontru_wrf, &
ctrlat_wrf,ctrlon_wrf,dx_wrfnm,dy_wrfnm,rdummy, &
rdummy,rdummy,rdummy, &
sfcinitopt,idumar,sfcdtfn,geogdir,start_date,silwt,wvln, &
idummy,idummy,idummy,idummy,spec_bdy_width, &
idummy,idummy,rdumar,rdumar,mp_physics,idumar, &
idumar,idumar,idumar,idumar, &
idumar,idummy,idummy,idumar, &
idummy,rdumar,rdumar,idummy,idummy,idumar, &
idummy,idummy,io_form,idummy,idummy,cdummy, &
idummy,idummy,idumar,cdummy,cdummy,cdummy,ldummy,istatus)
trulat1_wrf = lattru_wrf(1)
trulat2_wrf = lattru_wrf(2)
trulon_wrf = lontru_wrf
nz_wrf = 1
nzsoil_wrf = 1
io_form = 7
nproc_x = 1
nproc_y = 1
loc_x = 1
loc_y = 1
myproc = 0
readsplit = 0
execname = 'WRFSTATIC'
!--------------------------------------------------------------------
!
! Begin to set WRF horizontal grid
!
!--------------------------------------------------------------------
DO domid = 1, max_dom
WRITE(6,'(1x,a,I2,3a)') '========= Begin of domain: ',domid, &
', for file - ',TRIM(sfcdtfn(domid)),' ==========='
WRITE(grdbasfn,'(a)') adasbasfn(domid)
IF (use_arps_grid(domid) == 1) THEN ! reset nx_wrf,ny_wrf
lengbf = LEN_TRIM(grdbasfn)
WRITE(6,'(1x,/2a/)') ' Getting dimensions from file: ',grdbasfn(1:lengbf)
CALL get_arps_dims_atts(nchr,hinfmt(domid),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('get_arps_dims_atts errors.',1)
nx_wrf = nx - 2
ny_wrf = ny - 2
dx_wrf = dx
dy_wrf = dy
mapproj_wrf = mapproj
sclfct_wrf = sclfct
trulat1_wrf = trulat1
trulat2_wrf = trulat2
trulon_wrf = trulon
ctrlat_wrf = ctrlat
ctrlon_wrf = ctrlon
nx_wrfnm(domid) = nx_wrf
ny_wrfnm(domid) = ny_wrf
dx_wrfnm(domid) = dx_wrf
dy_wrfnm(domid) = dy_wrf
IF (domid == 1) THEN
mapproj_moad = mapproj
trulat1_moad = trulat1
trulat2_moad = trulat2
trulon_moad = trulon
ctrlat_moad = ctrlat_wrf
ELSE
IF ( ABS(mapproj_wrf-mapproj_moad) > a_small_real_number .OR. &
ABS(trulat1_wrf-trulat1_moad) > a_small_real_number .OR. &
ABS(trulat2_wrf-trulat2_moad) > a_small_real_number .OR. &
ABS(trulon_wrf -trulon_moad) > a_small_real_number ) THEN
WRITE(6,'(1x,a,I2,a,/)') &
'Inconsistent map project for subdomian ',domid,' with the mother domain.'
WRITE(6,'(a38)') ' Mother grid subdomain grid'
WRITE(6,'(a38)') ' ========== =========='
WRITE(6,'(a,I10,4x,I10)') ' mapproj = ',mapproj_moad,mapproj_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
CALL arpsstop('Inconsistent map projection',1)
END IF
END IF
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(domid),'(I4.4,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2)') &
year,'-',month,'-',day,'_',hour,':',minute,':',second
ELSE ! use parameters got from namelist above.
nx_wrf = nx_wrfnm(domid)
ny_wrf = ny_wrfnm(domid)
dx_wrf = dx_wrfnm(domid)
dy_wrf = dy_wrfnm(domid)
WRITE(6,'(/a/)') ' The WRF grid settings are:'
WRITE(6,'(a)') ' WRF grid'
WRITE(6,'(a)') ' =========='
WRITE(6,'(a,I10)') ' nx_wrf = ',nx_wrf
WRITE(6,'(a,I10)') ' ny_wrf = ',ny_wrf
WRITE(6,'(a,I10)') ' mapproj = ',mapproj_wrf
WRITE(6,'(a,F10.2)') ' sclfct = ',sclfct_wrf
WRITE(6,'(a,F10.2)') ' trulat1 = ',trulat1_wrf
WRITE(6,'(a,F10.2)') ' trulat2 = ',trulat2_wrf
WRITE(6,'(a,F10.2)') ' trulon = ',trulon_wrf
WRITE(6,'(a,F10.2)') ' ctrlat = ',ctrlat_wrf
WRITE(6,'(a,F10.2)') ' ctrlon = ',ctrlon_wrf
WRITE(6,'(a,F10.0)') ' dx_wrf = ',dx_wrf
WRITE(6,'(a,F10.0)') ' dy_wrf = ',dy_wrf
END IF
CALL get_check_grid_ratio(domid,nx_wrf,ny_wrf,dx_wrf,dy_wrf, &
dx_wrfnm(parent_id(domid)),dy_wrfnm(parent_id(domid)), &
a_small_real_number,parent_grid_ratio,istatus)
IF (istatus /= 0) CALL arpsstop('Unacceptable parent_grid_ratio or nesting grid size.',1)
!---------------------------------------------------------------------------
!
! Allocate WRF arrays
!
!---------------------------------------------------------------------------
IF (ALLOCATED(x_wrf)) DEALLOCATE(x_wrf)
IF (ALLOCATED(y_wrf)) DEALLOCATE(y_wrf)
IF (ALLOCATED(xs_wrf)) DEALLOCATE(xs_wrf)
IF (ALLOCATED(ys_wrf)) DEALLOCATE(ys_wrf)
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)
IF (ALLOCATED(lat_wrf)) DEALLOCATE(lat_wrf)
IF (ALLOCATED(lon_wrf)) DEALLOCATE(lon_wrf)
IF (ALLOCATED(msft_wrf)) DEALLOCATE(msft_wrf)
IF (ALLOCATED(msfu_wrf)) DEALLOCATE(msfu_wrf)
IF (ALLOCATED(msfv_wrf)) DEALLOCATE(msfv_wrf)
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)
IF (ALLOCATED(f)) DEALLOCATE(f)
IF (ALLOCATED(e)) DEALLOCATE(e)
IF (ALLOCATED(cosalpha)) DEALLOCATE(cosalpha)
IF (ALLOCATED(sinalpha)) DEALLOCATE(sinalpha)
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)
IF (ALLOCATED(hgt)) DEALLOCATE(hgt)
IF (ALLOCATED(topostdv)) DEALLOCATE(topostdv)
IF (ALLOCATED(toposlpx)) DEALLOCATE(toposlpx)
IF (ALLOCATED(toposlpy)) DEALLOCATE(toposlpy)
IF (ALLOCATED(landuse)) DEALLOCATE(landuse)
IF (ALLOCATED(landmask)) DEALLOCATE(landmask)
IF (ALLOCATED(slopecat)) DEALLOCATE(slopecat)
IF (ALLOCATED(tmn)) DEALLOCATE(tmn)
IF (ALLOCATED(snoalb)) DEALLOCATE(snoalb)
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)
IF (ALLOCATED(shdmax)) DEALLOCATE(shdmax)
IF (ALLOCATED(shdmin)) DEALLOCATE(shdmin)
IF (ALLOCATED(landusef)) DEALLOCATE(landusef)
IF (ALLOCATED(soilctop)) DEALLOCATE(soilctop)
IF (ALLOCATED(soilcbot)) DEALLOCATE(soilcbot)
IF (ALLOCATED(green12m)) DEALLOCATE(green12m)
IF (ALLOCATED(albdo12m)) DEALLOCATE(albdo12m)
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)
IF (ALLOCATED(tem2d1)) DEALLOCATE(tem2d1)
IF (ALLOCATED(tem2d2)) DEALLOCATE(tem2d2)
IF (ALLOCATED(tem2d3)) DEALLOCATE(tem2d3)
IF (ALLOCATED(tem3d)) DEALLOCATE(tem3d)
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
IF (domid == 1) THEN
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(domid) = 0
ysub0(domid) = 0
ELSE
IF (use_arps_grid(domid) == 1) THEN
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
xsub0(domid) = swx
ysub0(domid) = swy
ctrx = ( xsub0(domid)-xsub0(parent_id(domid)) ) / dx_wrfnm(parent_id(domid)) + 1
ctry = ( ysub0(domid)-ysub0(parent_id(domid)) ) / dy_wrfnm(parent_id(domid)) + 1
i_parent_start(domid) = NINT(ctrx)
j_parent_start(domid) = NINT(ctry)
IF ( ABS(ctrx - i_parent_start(domid)) > a_small_real_number*i_parent_start(domid) .OR. &
ABS(ctry - j_parent_start(domid)) > a_small_real_number*j_parent_start(domid) ) THEN
WRITE(6,'(1x,a,I2,2(a,F12.2),/,1x,2(a,F12.2),a,I2,/,2(1x,2(a,F12.2),a,/))') &
'Cannot align subdomain ',domid,' with central lat/lon = ',ctrlat_wrf,', ',ctrlon_wrf, &
'and low-left corner = ',swx,', ',swy,', inside the parent grid ',parent_id(domid), &
'with grid spacing = ',dx_wrfnm(parent_id(domid)),', ',dy_wrfnm(parent_id(domid)),'.',&
'The parent grid is starting from ',xsub0(parent_id(domid)),', ',ysub0(parent_id(domid)),'.'
CALL arpsstop('Inconsistent nesting grid',1)
END IF
WRITE(6,'(14x,2(a,i6),/)') 'Starting index are: ',i_parent_start(domid), &
', ',j_parent_start(domid)
ELSE
xsub0(domid) = xsub0(parent_id(domid)) + (i_parent_start(domid)-1)*dx_wrfnm(parent_id(domid))
ysub0(domid) = ysub0(parent_id(domid)) + (j_parent_start(domid)-1)*dy_wrfnm(parent_id(domid))
ctrx = xsub0(domid) + 0.5*(nx_wrf-1)*dx_wrfscl
ctry = ysub0(domid) + 0.5*(ny_wrf-1)*dy_wrfscl
CALL xytoll( 1,1, ctrx, ctry, ctrlat_wrf, ctrlon_wrf)
WRITE(6,'(14x,2(a,F12.7),/)') 'domain center lat/lon is: ',ctrlat_wrf, ', ', ctrlon_wrf
END IF
END IF
i_parent_end = i_parent_start(domid) + (nx_wrf-1)/parent_grid_ratio
j_parent_end = j_parent_start(domid) + (ny_wrf-1)/parent_grid_ratio
IF (i_parent_end > nx_wrfnm(parent_id(domid)) .OR. &
j_parent_end > ny_wrfnm(parent_id(domid)) ) THEN
WRITE(6,'(1x,2(a,I2),a)') 'Subdomain ',domid,' exceed its parent domain ', &
parent_id(domid),' in size.'
WRITE(6,'(1x,4(a,I8),a)') 'i_parent_start-end, j_parent_start-end = ', &
i_parent_start(domid),' - ',i_parent_end,', ', &
j_parent_start(domid),' - ',j_parent_end,'.'
WRITE(6,'(1x,a,I2,a,2(I8,a),/)') 'Parent domain ',parent_id(domid),' size: ', &
nx_wrfnm(parent_id(domid)),', ', &
ny_wrfnm(parent_id(domid)),'.'
CALL arpsstop('Subdomain is too large.',1)
END IF
DO i=1,nx_wrf
x_wrf(i)= sclf_wrf*xsub0(domid) + (i-1)*dx_wrfscl
END DO
DO j=1,ny_wrf
y_wrf(j)= sclf_wrf*ysub0(domid) + (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
times_str = start_date(domid)
WRITE(6,'(/2a/)') ' Preparing static data set starting at ',times_str
!---------------------- Terrain USGS 30 sec --------------------------
! type = U
! symbol = TOPO
lenstr = LEN_TRIM(geogdir)
WRITE(tmpstr,'(3a)') geogdir(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)') geogdir(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)') geogdir(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)') geogdir(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)') geogdir(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)') geogdir(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)') geogdir(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)') geogdir(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)') geogdir(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 = sfcdtfn(domid)
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,ncid)
!-----------------------------------------------------------------------
!
! Initialize and write global attributes
!
!-----------------------------------------------------------------------
CALL write_static_attribute(ncid,io_form,start_date(domid),domid, &
parent_id(domid),i_parent_start(domid),j_parent_start(domid), &
i_parent_end,j_parent_end,parent_grid_ratio, &
nx_wrf,ny_wrf,dx_wrf,dy_wrf, &
mapproj_wrf,trulat1_wrf,trulat2_wrf,trulon_wrf,&
ctrlat_moad,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(sfcdtfn(domid)),"_ready"
WRITE(filename,'(a)') TRIM(dirname) // TRIM(filename)
OPEN (UNIT=nchout,FILE=trim(filename))
WRITE (nchout,'(a)') trim(sfcdtfn(domid))
CLOSE (nchout)
CALL retunit (nchout )
END IF
END DO ! max_dom
CALL arpsstop(' ==== WRFSTATIC terminated normally. ====',0)
END PROGRAM wrfstatic