Software HRLDAS land surface data assimilation example single point

HRLDAS模型单点案例运行

Posted by CryoECNU on October 17, 2021

参考编译HRLDAS模型：【HRLDAS模型编译】。

编译完成后，有两个exe文件：HRLDAS_forcing/create_forcing.exe和run/hrldas.exe

这里先从最简单的单点案例开始。在HRLDAS_forcing/run/examples/single_point/目录下。

需要用NCL制作所需驱动数据文件，也就是运行这个目录下的两个NCL文件： create_setup.ncl 和 create_ldasin_files.ncl

从NLDAS的网站上下载一些数据：

1）NLDAS_NOAH0125_H.A20120601.0000.002.grb: create_setup.ncl会用这个文件来制作初始条件的输入。

2）一些气象驱动数据，类似于NLDAS_FORA0125_H.A20120601.0000.002.grb之类的文件。注意和上面文件的区别，这里是NLDAS_FORA0125_。

把这些数据都放在一个目录raw下面，另外创建目录LDASIN和LDASOUT。

第一步：修改create_setup.ncl文件。


inname = "raw/NLDAS_NOAH0125_H.A20120601.0000.002.grb"
setup_filename = "LDASIN/HRLDAS_setup_2012060100_d1"

需要给定grib的table表，在bash里设定环境变量NCL_GRIB_PTABLE_PATH：

> export NCL_GRIB_PTABLE_PATH=gribtab.nldas

然后，用NCL ncl create_setup.ncl，会在LDASIN下面生成HRLDAS_setup_2012060100_d1文件。

这个文件可以用ncdump查看一下：

netcdf HRLDAS_setup_2012060100_d1 {
dimensions:
	Time = UNLIMITED ; // (1 currently)
	south_north = 1 ;
	west_east = 1 ;
	soil_layers_stag = 4 ;
	DateStrLen = 19 ;
variables:
	char Times(Time, DateStrLen) ;
		Times:_FillValue = "" ;
	float XLAT(Time, south_north, west_east) ;
		XLAT:units = "degrees_north" ;
		XLAT:_FillValue = 9.96921e+36f ;
	float XLONG(Time, south_north, west_east) ;
		XLONG:units = "degrees_east" ;
		XLONG:_FillValue = 9.96921e+36f ;
	float TMN(Time, south_north, west_east) ;
		TMN:units = "K" ;
		TMN:_FillValue = 9.96921e+36f ;
	float HGT(Time, south_north, west_east) ;
		HGT:units = "m" ;
		HGT:_FillValue = 9.96921e+36f ;
	float SEAICE(Time, south_north, west_east) ;
		SEAICE:units = "" ;
		SEAICE:_FillValue = 9.96921e+36f ;
	float MAPFAC_MX(Time, south_north, west_east) ;
		MAPFAC_MX:units = "" ;
		MAPFAC_MX:_FillValue = 9.96921e+36f ;
	float MAPFAC_MY(Time, south_north, west_east) ;
		MAPFAC_MY:units = "" ;
		MAPFAC_MY:_FillValue = 9.96921e+36f ;
	float SHDMAX(Time, south_north, west_east) ;
		SHDMAX:units = "%" ;
		SHDMAX:_FillValue = 9.96921e+36f ;
	float SHDMIN(Time, south_north, west_east) ;
		SHDMIN:units = "%" ;
		SHDMIN:_FillValue = 9.96921e+36f ;
	float LAI(Time, south_north, west_east) ;
		LAI:units = "m^2/m^2" ;
		LAI:_FillValue = 9.96921e+36f ;
	int XLAND(Time, south_north, west_east) ;
		XLAND:units = "" ;
		XLAND:_FillValue = -2147483647 ;
	int IVGTYP(Time, south_north, west_east) ;
		IVGTYP:units = "" ;
		IVGTYP:_FillValue = -2147483647 ;
	int ISLTYP(Time, south_north, west_east) ;
		ISLTYP:units = "" ;
		ISLTYP:_FillValue = -2147483647 ;
	float SNOW(Time, south_north, west_east) ;
		SNOW:units = "kg/m^2" ;
		SNOW:_FillValue = 9.96921e+36f ;
	float CANWAT(Time, south_north, west_east) ;
		CANWAT:units = "kg/m^2" ;
		CANWAT:_FillValue = 9.96921e+36f ;
	float DZS(Time, soil_layers_stag) ;
		DZS:units = "m" ;
		DZS:_FillValue = 9.96921e+36f ;
	float ZS(Time, soil_layers_stag) ;
		ZS:units = "m" ;
		ZS:_FillValue = 9.96921e+36f ;
	float TSLB(Time, soil_layers_stag, south_north, west_east) ;
		TSLB:units = "K" ;
		TSLB:_FillValue = 9.96921e+36f ;
	float TSK(Time, south_north, west_east) ;
		TSK:units = "K" ;
		TSK:_FillValue = 9.96921e+36f ;
	float SMOIS(Time, soil_layers_stag, south_north, west_east) ;
		SMOIS:units = "m^3/m^3" ;
		SMOIS:_FillValue = 9.96921e+36f ;

// global attributes:
		:TITLE = "OUTPUT FROM VECTOR CREATION SCRIPTS: m.barlage v20150608" ;
		:DX = 1000.f ;
		:DY = 1000.f ;
		:TRUELAT1 = 45.f ;
		:TRUELAT2 = 45.f ;
		:STAND_LON = 45.f ;
		:MAP_PROJ = 1 ;
		:GRID_ID = 1 ;
		:ISWATER = 17 ;
		:ISURBAN = 13 ;
		:ISICE = 15 ;
		:MMINLU = "MODIFIED_IGBP_MODIS_NOAH" ;
data:

 Times =
  "2012-06-01_00:00:00" ;

 XLAT =
  40 ;

 XLONG =
  -105.25 ;

 TMN =
  10.9 ;

 HGT =
  1655 ;

 SEAICE =
  0 ;

 MAPFAC_MX =
  0 ;

 MAPFAC_MY =
  0 ;

 SHDMAX =
  100 ;

 SHDMIN =
  100 ;

 LAI =
  2 ;

 XLAND =
  1 ;

 IVGTYP =
  10 ;

 ISLTYP =
  9 ;

 SNOW =
  0 ;

 CANWAT =
  0 ;

 DZS =
  0.1, 0.3, 0.6, 1 ;

 ZS =
  0.05, 0.25, 0.7, 1.5 ;

 TSLB =
  295.14,
  283.53,
  280.13,
  273.59 ;

 TSK =
  295.14 ;

 SMOIS =
  0.2388,
  0.2729333,
  0.1840833,
  0.25646 ;
}

第二步，修改create_ldasin_files.ncl文件：

indir = "raw/"
outdir = "LDASIN/"

然后，用NCL ncl create_ldasin_files.ncl，会在LDASIN下面生成类似于2012060100.LDASIN_DOMAIN1的一系列文件。

第三步，创建namelist.hrldas文件:

&NOAHLSM_OFFLINE

 HRLDAS_SETUP_FILE = "LDASIN/HRLDAS_setup_2012060100_d1"
 ! MMF_RUNOFF_FILE = "/blah/WRF/wrfleafxyz_d01.nc"
 INDIR = "LDASIN"
 OUTDIR = "LDASOUT/"

 START_YEAR  = 2012
 START_MONTH = 06
 START_DAY   = 01
 START_HOUR  = 00
 START_MIN   = 00
 
! RESTART_FILENAME_REQUESTED = "/d1/barlage/data/mptest/vector/LDASOUT/RESTART.2012020200_DOMAIN1"

 KDAY = 2

 DYNAMIC_VEG_OPTION                = 4
 CANOPY_STOMATAL_RESISTANCE_OPTION = 1
 BTR_OPTION                        = 1
 RUNOFF_OPTION                     = 3
 SURFACE_DRAG_OPTION               = 1
 FROZEN_SOIL_OPTION                = 1
 SUPERCOOLED_WATER_OPTION          = 1
 RADIATIVE_TRANSFER_OPTION         = 3
 SNOW_ALBEDO_OPTION                = 2
 PCP_PARTITION_OPTION              = 1
 TBOT_OPTION                       = 2
 TEMP_TIME_SCHEME_OPTION           = 1
 GLACIER_OPTION                    = 1
 SURFACE_RESISTANCE_OPTION         = 1

 FORCING_TIMESTEP = 3600
 NOAH_TIMESTEP    = 3600
 OUTPUT_TIMESTEP  = 3600

 SPLIT_OUTPUT_COUNT = 0
 RESTART_FREQUENCY_HOURS = 0

 ! XSTART = 20
 ! XEND   = 20
 ! YSTART = 203
 ! YEND   = 203

 NSOIL=4
 soil_thick_input(1) = 0.10
 soil_thick_input(2) = 0.30
 soil_thick_input(3) = 0.60
 soil_thick_input(4) = 1.00

 ZLVL = 10.0

/

&URBAN_OFFLINE
 UCMCALL = 0
 ZLVL_URBAN = 15.0
/

第四步，运行hrldas.exe

可以把编译出来的hrldas.exe