<p><b>duda</b> 2011-02-08 17:24:15 -0700 (Tue, 08 Feb 2011)</p><p>BRANCH COMMIT<br>
<br>
Add new interpolated fields: <br>
<br>
snoalb - maximum annual snow albedo<br>
greenfrac - monthly vegetation fraction<br>
vegfra - vegetation fraction interpolated in time to initial simulation date (eventually)<br>
shdmin - minimum vegetation fraction<br>
shdmax - maximum vegetation fraction<br>
albedo12m - monthly albedo<br>
sst - sea-surface temperature<br>
snowc - snow cover flag (0/1)<br>
<br>
Also, store soil moisture and soil temperature in 2d arrays dimensioned<br>
(nSoilLevels nCells), rather than as several 1d (nCells) arrays.<br>
<br>
<br>
M src/core_init_nhyd_atmos/module_test_cases.F<br>
M src/core_init_nhyd_atmos/Registry<br>
M namelist.input.nhyd_realdata<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/namelist.input.nhyd_realdata
===================================================================
--- branches/atmos_physics/namelist.input.nhyd_realdata 2011-02-08 23:13:39 UTC (rev 720)
+++ branches/atmos_physics/namelist.input.nhyd_realdata 2011-02-09 00:24:15 UTC (rev 721)
@@ -1,47 +1,29 @@
&nhyd_model
config_test_case = 7
- config_time_integration = 'SRK3'
- config_dt = 450
- config_ntimesteps = 1920
- config_output_interval = 192
- config_number_of_sub_steps = 6
- config_h_mom_eddy_visc2 = 0.0e+04
- config_h_mom_eddy_visc4 = 0.
- config_v_mom_eddy_visc2 = 00.0
- config_h_theta_eddy_visc2 = 0.0e+04
- config_h_theta_eddy_visc4 = 00.
- config_v_theta_eddy_visc2 = 00.0
- config_horiz_mixing = '2d_smagorinsky'
- config_len_disp = 60000.
- config_u_vadv_order = 3
- config_w_vadv_order = 3
- config_theta_vadv_order = 3
- config_scalar_vadv_order = 3
- config_theta_adv_order = 3
config_scalar_adv_order = 3
- config_scalar_advection = .false.
- config_positive_definite = .false.
- config_coef_3rd_order = 1.0
- config_monotonic = .false.
- config_epssm = 0.1
- config_smdiv = 0.1
/
&dimensions
- config_nvertlevels = 26
+ config_nvertlevels = 35
config_nfglevels = 27
+ config_nsoil_levels = 4
/
&data_sources
config_geog_data_path = '/data3/mp/wrfhelp/WPS_GEOG/'
- config_met_prefix = 'GFS'
- config_init_date = '2010-05-26_12'
+ config_met_prefix = 'FILE'
+ config_init_date = '2011-01-22_00'
/
+&preproc_stages
+ config_static_interp = .true.
+ config_vertical_grid = .true.
+ config_met_interp = .true.
+/
&io
config_input_name = 'grid.nc'
- config_output_name = 'output.nc'
+ config_output_name = 'init.nc'
config_restart_name = 'restart.nc'
/
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/Registry
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-02-08 23:13:39 UTC (rev 720)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/Registry 2011-02-09 00:24:15 UTC (rev 721)
@@ -32,6 +32,7 @@
namelist real nhyd_model config_smdiv 0.1
namelist integer dimensions config_nvertlevels 26
namelist integer dimensions config_nfglevels 27
+namelist integer dimensions config_nsoil_levels 4
namelist character data_sources config_geog_data_path /data3/mp/wrfhelp/WPS_GEOG/
namelist character data_sources config_met_prefix GFS
namelist character data_sources config_init_date 2010-05-26_12
@@ -61,8 +62,10 @@
dim TWENTYONE 21
dim R3 3
dim nVertLevels namelist:config_nvertlevels
+dim nSoilLevels namelist:config_nsoil_levels
dim nFGLevels namelist:config_nfglevels
dim nVertLevelsP1 nVertLevels+1
+dim nMonths 12
#
# var type name_in_file ( dims ) iro- name_in_code super-array array_class
@@ -126,11 +129,17 @@
var persistent real cf3 ( ) 0 io cf3 mesh - -
# static terrestrial fields
-var persistent real ter ( nCells ) 0 io ter mesh - -
-var persistent integer landmask ( nCells ) 0 io landmask mesh - -
-var persistent integer lu_index ( nCells ) 0 io lu_index mesh - -
+var persistent real ter ( nCells ) 0 io ter mesh - -
+var persistent integer landmask ( nCells ) 0 io landmask mesh - -
+var persistent integer lu_index ( nCells ) 0 io lu_index mesh - -
var persistent integer soilcat_top ( nCells ) 0 io soilcat_top mesh - -
var persistent integer soilcat_bot ( nCells ) 0 io soilcat_bot mesh - -
+var persistent real snoalb ( nCells ) 0 io snoalb mesh - -
+var persistent real greenfrac ( nMonths nCells ) 0 io greenfrac mesh - -
+var persistent real shdmin ( nCells ) 0 io shdmin mesh - -
+var persistent real shdmax ( nCells ) 0 io shdmax mesh - -
+var persistent real vegfra ( nCells ) 0 io vegfra mesh - -
+var persistent real albedo12m ( nMonths nCells ) 0 io albedo12m mesh - -
# description of the vertical grid structure
@@ -165,19 +174,16 @@
# Horizontally interpolated from first-guess data
# and should be read in by model
-var persistent real sm000010 ( nCells ) 0 o sm000010 fg - -
-var persistent real sm010040 ( nCells ) 0 o sm010040 fg - -
-var persistent real sm040100 ( nCells ) 0 o sm040100 fg - -
-var persistent real sm100200 ( nCells ) 0 o sm100200 fg - -
-var persistent real st000010 ( nCells ) 0 o st000010 fg - -
-var persistent real st010040 ( nCells ) 0 o st010040 fg - -
-var persistent real st040100 ( nCells ) 0 o st040100 fg - -
-var persistent real st100200 ( nCells ) 0 o st100200 fg - -
+var persistent real dzs ( nSoilLevels nCells ) 0 o dzs fg - -
+var persistent real sh2o ( nSoilLevels nCells ) 0 o sh2o fg - -
+var persistent real smois ( nSoilLevels nCells ) 0 o smois fg - -
+var persistent real tslb ( nSoilLevels nCells ) 0 o tslb fg - -
var persistent real skintemp ( nCells ) 0 o skintemp fg - -
+var persistent real sst ( nCells ) 0 o sst fg - -
var persistent real snow ( nCells ) 0 o snow fg - -
-var persistent real seaice ( nCells ) 0 o seaice fg - -
+var persistent real snowc ( nCells ) 0 o snowc fg - -
+var persistent real xice ( nCells ) 0 o xice fg - -
var persistent real gfs_z ( nVertLevels nCells ) 0 - gfs_z fg - -
-#var persistent real gfs_p ( nVertLevelsP1 nCells ) 0 - gfs_p fg - -
# Prognostic variables: read from input, saved in restart, and written to output
var persistent real u ( nVertLevels nEdges Time ) 2 io u state - -
Modified: branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-02-08 23:13:39 UTC (rev 720)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/module_test_cases.F 2011-02-09 00:24:15 UTC (rev 721)
@@ -2032,7 +2032,9 @@
integer, allocatable, dimension(:,:) :: ncat
real (kind=4) :: scalefactor ! NB: this should be a single-precision real
real (kind=RKIND) :: lat_pt, lon_pt, lon_pt_o
- real (kind=4), allocatable, dimension(:,:,:) :: rarray ! NB: this should be a single-precision real array
+ real (kind=4), allocatable, dimension(:,:,:) :: rarray ! NB: this should be a single-precision real array
+ real (kind=4), allocatable, dimension(:,:) :: maxsnowalb ! NB: this should be a single-precision real array
+ real (kind=4), allocatable, dimension(:,:,:) :: vegfra ! NB: this should be a single-precision real array
character (len=1024) :: fname
real (kind=RKIND) :: u, flux, fluxk, lat1, lat2, eta_v, r_pert, u_pert, lat_pert, lon_pert, r
@@ -2410,6 +2412,196 @@
deallocate(rarray)
deallocate(ncat)
+
+
+ !
+ ! Interpolate SNOALB
+ !
+ nx = 186
+ ny = 186
+ nzz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nzz))
+ allocate(maxsnowalb(360,180))
+ grid % snoalb % array(:) = 0.0
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = 1.0_4, &
+ loninc = 1.0_4, &
+ knowni = 1.0_4, &
+ knownj = 1.0_4, &
+ lat1 = -89.5_4, &
+ lon1 = -179.5_4)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'maxsnowalb/',1,'-',180,'.',1,'-',180
+write(0,*) trim(fname)
+ call read_geogrid(fname, len_trim(fname), &
+ rarray, &
+ nx, ny, nzz, &
+ isigned, endian, scalefactor, wordsize, istatus)
+write(0,*) istatus
+write(0,*) 'min/max = ',minval(rarray(:,:,:)),maxval(rarray(:,:,:))
+ maxsnowalb(1:180,1:180) = rarray(4:183,4:183,1)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'maxsnowalb/',181,'-',360,'.',1,'-',180
+write(0,*) trim(fname)
+ call read_geogrid(fname, len_trim(fname), &
+ rarray, &
+ nx, ny, nzz, &
+ isigned, endian, scalefactor, wordsize, istatus)
+write(0,*) istatus
+write(0,*) 'min/max = ',minval(rarray(:,:,:)),maxval(rarray(:,:,:))
+ maxsnowalb(181:360,1:180) = rarray(4:183,4:183,1)
+
+ do iCell=1,grid%nCells
+ lat = grid % latCell % array(iCell)*DEG_PER_RAD
+ lon = grid % lonCell % array(iCell)*DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ if (x < 0.0) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ grid % snoalb % array(iCell) = four_pt(360, 180, maxsnowalb, x, y)
+ end do
+
+ grid % snoalb % array(:) = grid % snoalb % array(:) / 100.0
+
+ deallocate(rarray)
+ deallocate(maxsnowalb)
+
+
+ !
+ ! Interpolate GREENFRAC
+ !
+ nx = 1256
+ ny = 1256
+ nzz = 12
+ isigned = 0
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nzz))
+ allocate(vegfra(2500,1250,12))
+ grid % vegfra % array(:) = 0.0
+ grid % greenfrac % array(:,:) = 0.0
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = 0.144_4, &
+ loninc = 0.144_4, &
+ knowni = 1.0_4, &
+ knownj = 1.0_4, &
+ lat1 = -89.928_4, &
+ lon1 = -179.928_4)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'greenfrac/',1,'-',1250,'.',1,'-',1250
+write(0,*) trim(fname)
+ call read_geogrid(fname, len_trim(fname), &
+ rarray, &
+ nx, ny, nzz, &
+ isigned, endian, scalefactor, wordsize, istatus)
+write(0,*) istatus
+write(0,*) 'min/max = ',minval(rarray(:,:,:)),maxval(rarray(:,:,:))
+ vegfra(1:1250,1:1250,1:12) = rarray(4:1253,4:1253,1:12)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'greenfrac/',1251,'-',2500,'.',1,'-',1250
+write(0,*) trim(fname)
+ call read_geogrid(fname, len_trim(fname), &
+ rarray, &
+ nx, ny, nzz, &
+ isigned, endian, scalefactor, wordsize, istatus)
+write(0,*) istatus
+write(0,*) 'min/max = ',minval(rarray(:,:,:)),maxval(rarray(:,:,:))
+ vegfra(1251:2500,1:1250,1:12) = rarray(4:1253,4:1253,1:12)
+
+ do iCell=1,grid%nCells
+ lat = grid % latCell % array(iCell)*DEG_PER_RAD
+ lon = grid % lonCell % array(iCell)*DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ if (x < 0.0) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ do k=1,12
+ grid % greenfrac % array(k,iCell) = four_pt(2500, 1250, vegfra(:,:,k), x, y)
+ end do
+ grid % shdmin % array(iCell) = minval(grid % greenfrac % array(:,iCell))
+ grid % shdmax % array(iCell) = maxval(grid % greenfrac % array(:,iCell))
+ end do
+
+ grid % greenfrac % array(:,:) = grid % greenfrac % array(:,:) / 100.0
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!TO-DO: The vegfra field should be interpolated in time from the greenfrac field
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ grid % vegfra % array(:) = grid % greenfrac % array(3,:)
+
+ deallocate(rarray)
+ deallocate(vegfra)
+
+
+ !
+ ! Interpolate ALBEDO12M
+ !
+ nx = 1256
+ ny = 1256
+ nzz = 12
+ isigned = 0
+ endian = 0
+ wordsize = 1
+ scalefactor = 1.0
+ allocate(rarray(nx,ny,nzz))
+ allocate(vegfra(2500,1250,12))
+ grid % albedo12m % array(:,:) = 0.0
+
+ call map_set(PROJ_LATLON, proj, &
+ latinc = 0.144_4, &
+ loninc = 0.144_4, &
+ knowni = 1.0_4, &
+ knownj = 1.0_4, &
+ lat1 = -89.928_4, &
+ lon1 = -179.928_4)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'albedo_ncep/',1,'-',1250,'.',1,'-',1250
+write(0,*) trim(fname)
+ call read_geogrid(fname, len_trim(fname), &
+ rarray, &
+ nx, ny, nzz, &
+ isigned, endian, scalefactor, wordsize, istatus)
+write(0,*) istatus
+write(0,*) 'min/max = ',minval(rarray(:,:,:)),maxval(rarray(:,:,:))
+ vegfra(1:1250,1:1250,1:12) = rarray(4:1253,4:1253,1:12)
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'albedo_ncep/',1251,'-',2500,'.',1,'-',1250
+write(0,*) trim(fname)
+ call read_geogrid(fname, len_trim(fname), &
+ rarray, &
+ nx, ny, nzz, &
+ isigned, endian, scalefactor, wordsize, istatus)
+write(0,*) istatus
+write(0,*) 'min/max = ',minval(rarray(:,:,:)),maxval(rarray(:,:,:))
+ vegfra(1251:2500,1:1250,1:12) = rarray(4:1253,4:1253,1:12)
+
+ do iCell=1,grid%nCells
+ lat = grid % latCell % array(iCell)*DEG_PER_RAD
+ lon = grid % lonCell % array(iCell)*DEG_PER_RAD
+ call latlon_to_ij(proj, lat, lon, x, y)
+ if (x < 0.0) then
+ lon = lon + 360.0
+ call latlon_to_ij(proj, lat, lon, x, y)
+ end if
+ do k=1,12
+ grid % albedo12m % array(k,iCell) = four_pt(2500, 1250, vegfra(:,:,k), x, y)
+ end do
+ end do
+
+ grid % albedo12m % array(:,:) = grid % albedo12m % array(:,:) / 100.0
+
+ deallocate(rarray)
+ deallocate(vegfra)
end if ! config_static_interp
@@ -2919,57 +3111,65 @@
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
- destField1d => fg % sm000010 % array
- ndims = 1
+ destField2d => fg % smois % array
+ k = 1
+ ndims = 2
else if (index(field % field, 'SM010040') /= 0) then
write(0,*) 'Interpolating SM010040'
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
- destField1d => fg % sm010040 % array
- ndims = 1
+ destField2d => fg % smois % array
+ k = 2
+ ndims = 2
else if (index(field % field, 'SM040100') /= 0) then
write(0,*) 'Interpolating SM040100'
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
- destField1d => fg % sm040100 % array
- ndims = 1
+ destField2d => fg % smois % array
+ k = 3
+ ndims = 2
else if (index(field % field, 'SM100200') /= 0) then
write(0,*) 'Interpolating SM100200'
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
- destField1d => fg % sm100200 % array
- ndims = 1
+ destField2d => fg % smois % array
+ k = 4
+ ndims = 2
else if (index(field % field, 'ST000010') /= 0) then
write(0,*) 'Interpolating ST000010'
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
- destField1d => fg % st000010 % array
- ndims = 1
+ destField2d => fg % tslb % array
+ k = 1
+ ndims = 2
else if (index(field % field, 'ST010040') /= 0) then
write(0,*) 'Interpolating ST010040'
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
- destField1d => fg % st010040 % array
- ndims = 1
+ destField2d => fg % tslb % array
+ k = 2
+ ndims = 2
else if (index(field % field, 'ST040100') /= 0) then
write(0,*) 'Interpolating ST040100'
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
- destField1d => fg % st040100 % array
- ndims = 1
+ destField2d => fg % tslb % array
+ k = 3
+ ndims = 2
else if (index(field % field, 'ST100200') /= 0) then
write(0,*) 'Interpolating ST100200'
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
- destField1d => fg % st100200 % array
- ndims = 1
+ destField2d => fg % tslb % array
+ k = 4
+ ndims = 2
else if (index(field % field, 'SNOW') /= 0) then
write(0,*) 'Interpolating SNOW'
nInterpPoints = grid % nCells
@@ -2982,7 +3182,7 @@
nInterpPoints = grid % nCells
latPoints => grid % latCell % array
lonPoints => grid % lonCell % array
- destField1d => fg % seaice % array
+ destField1d => fg % xice % array
ndims = 1
else if (index(field % field, 'SKINTEMP') /= 0) then
write(0,*) 'Interpolating SKINTEMP'
@@ -3007,6 +3207,10 @@
destField2d(k,i) = four_pt(field % nx, field % ny, field % slab, x, y)
end if
end do
+if (index(field % field, 'SEAICE') /= 0) then
+ write(0,*) minval(destfield1d(:)), maxval(destfield1d(:))
+end if
+
end if
deallocate(field % slab)
@@ -3027,7 +3231,17 @@
end do
end if
+ ! Set SST based on SKINTEMP field if it wasn't found in input data
+ if (minval(fg % sst % array) == 0.0 .and. maxval(fg % sst % array) == 0.0) then
+ write(0,*) 'Setting SST from SKINTEMP'
+ where (grid % landmask % array == 0) fg % sst % array = fg % skintemp % array
+ end if
+ ! Set SNOWC (snow-cover flag) based on SNOW
+ fg % snowc % array(:) = 0.0
+ where (fg % snow % array > 0.0) fg % snowc % array = 1.0
+
+
!
! Compute normal wind component and store in fg%u
!
</font>
</pre>