<p><b>laura@ucar.edu</b> 2013-01-30 14:39:07 -0700 (Wed, 30 Jan 2013)</p><p>Removed the interpolation of the static fields from subroutine init_atm_test_case_gfs to subroutine init_atm_static (in mpas_init_atm_static.F). Added separate calls to the subroutines init_atm_static and static_atm_static_orogwd when config_test_case equals 7 and config_static_interp is true.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2013-01-30 21:35:31 UTC (rev 2394)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2013-01-30 21:39:07 UTC (rev 2395)
@@ -10,6 +10,7 @@
use mpas_RBF_interpolation
use mpas_vector_reconstruction
use mpas_timer
+ use mpas_init_atm_static
use mpas_init_atm_surface
! Added only clause to keep xlf90 from getting confused from the overloaded abs intrinsic in mpas_timekeeping
@@ -105,10 +106,15 @@
write(0,*) ' real-data GFS test case '
block_ptr => domain % blocklist
do while (associated(block_ptr))
+ if (config_static_interp) then
+ call init_atm_static(block_ptr % mesh)
+ call init_atm_static_orogwd(block_ptr % mesh)
+ endif
call init_atm_test_case_gfs(block_ptr % mesh, block_ptr % fg, &
block_ptr % state % time_levs(1) % state, block_ptr % diag, &
config_test_case)
if (config_met_interp) call physics_initialize_real(block_ptr % mesh, block_ptr % fg, domain % dminfo)
+
block_ptr => block_ptr % next
end do
@@ -2292,10 +2298,9 @@
!This is temporary variable here. It just need when calculate tangential velocity v.
integer :: eoe, j
- integer, dimension(:), pointer :: nEdgesOnEdge, nEdgesOnCell
+ integer, dimension(:), pointer :: nEdgesOnCell
integer, dimension(:,:), pointer :: edgesOnEdge, cellsOnEdge, edgesOnCell, cellsOnCell
real (kind=RKIND), dimension(:), pointer :: dvEdge, dcEdge, AreaCell
- real (kind=RKIND), dimension(:,:), pointer :: weightsOnEdge
real (kind=RKIND), dimension(:,:), pointer :: v
real (kind=RKIND), dimension(:,:), pointer :: sorted_arr
@@ -2311,19 +2316,10 @@
real (kind=RKIND) :: p_check
! For interpolating terrain and land use
- integer :: nx, ny, nzz, iPoint, subx, suby
- integer :: isigned, endian, wordsize, istatus
- integer :: iTileStart, iTileEnd
- integer :: jTileStart, jTileEnd
- integer, allocatable, dimension(:) :: nhs
- 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
+ integer :: nx, ny
+ integer :: istatus
+
real (kind=RKIND), allocatable, dimension(:,:) :: rslab, maskslab
- real (kind=RKIND), allocatable, dimension(:,:) :: maxsnowalb
- real (kind=RKIND), allocatable, dimension(:,:) :: soiltemp_1deg
- real (kind=RKIND), allocatable, dimension(:,:,:) :: vegfra
integer, dimension(:), pointer :: mask_array
integer, dimension(grid % nEdges), target :: edge_mask
character (len=StrKIND) :: fname
@@ -2365,8 +2361,6 @@
parinfo => block % parinfo
dminfo => block % domain % dminfo
- weightsOnEdge => grid % weightsOnEdge % array
- nEdgesOnEdge => grid % nEdgesOnEdge % array
nEdgesOnCell => grid % nEdgesOnCell % array
edgesOnEdge => grid % edgesOnEdge % array
edgesOnCell => grid % edgesOnCell % array
@@ -2422,622 +2416,9 @@
omega_e = omega
p0 = 1.e+05
- interp_list(1) = FOUR_POINT
- interp_list(2) = SEARCH
- interp_list(3) = 0
-
-
- !
- ! Scale all distances and areas from a unit sphere to one with radius sphere_radius
- !
-
- if (config_static_interp) then
-
- grid % xCell % array = grid % xCell % array * r_earth
- grid % yCell % array = grid % yCell % array * r_earth
- grid % zCell % array = grid % zCell % array * r_earth
- grid % xVertex % array = grid % xVertex % array * r_earth
- grid % yVertex % array = grid % yVertex % array * r_earth
- grid % zVertex % array = grid % zVertex % array * r_earth
- grid % xEdge % array = grid % xEdge % array * r_earth
- grid % yEdge % array = grid % yEdge % array * r_earth
- grid % zEdge % array = grid % zEdge % array * r_earth
- grid % dvEdge % array = grid % dvEdge % array * r_earth
- grid % dcEdge % array = grid % dcEdge % array * r_earth
- grid % areaCell % array = grid % areaCell % array * r_earth**2.0
- grid % areaTriangle % array = grid % areaTriangle % array * r_earth**2.0
- grid % kiteAreasOnVertex % array = grid % kiteAreasOnVertex % array * r_earth**2.0
-
scalars(:,:,:) = 0.
- call atm_initialize_advection_rk(grid)
- call atm_initialize_deformation_weights(grid)
-
-
- !
- ! Interpolate HGT
- !
-! nx = 126
-! ny = 126
- nx = 1206
- ny = 1206
- nzz = 1
- isigned = 1
- endian = 0
- wordsize = 2
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(nhs(grid % nCells))
- nhs(:) = 0
- ter(:) = 0.0
-
- do jTileStart=1,20401,ny-6
-! do jTileStart=1,961,ny-6
- jTileEnd = jTileStart + ny - 1 - 6
- do iTileStart=1,42001,nx-6
-! do iTileStart=1,2041,nx-6
- iTileEnd = iTileStart + nx - 1 - 6
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'topo_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-! write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'topo_10m/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-write(0,*) trim(fname)
-
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- iPoint = 1
- do j=4,ny-3
- do i=4,nx-3
- lat_pt = (-90.0 + 1.0/240.0) + (jTileStart + j - 5) / 120.0
- lon_pt = (-180.0 + 1.0/240.0) + (iTileStart + i - 5) / 120.0
-! lat_pt = -89.91667 + (jTileStart + j - 5) * 0.166667
-! lon_pt = -179.91667 + (iTileStart + i - 5) * 0.166667
- lat_pt = lat_pt * pii / 180.0
- lon_pt = lon_pt * pii / 180.0
-
- iPoint = nearest_cell(lat_pt, lon_pt, &
- iPoint, &
- grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
- grid % latCell % array, grid % lonCell % array)
-
- ter(iPoint) = ter(iPoint) + rarray(i,j,1)
- nhs(iPoint) = nhs(iPoint) + 1
-
- end do
- end do
-
- end do
- end do
-
- do iCell=1, grid % nCells
- ter(iCell) = ter(iCell) / real(nhs(iCell))
- end do
-
- deallocate(rarray)
- deallocate(nhs)
-
-
- !
- ! Interpolate LU_INDEX
- !
- nx = 1200
- ny = 1200
- nzz = 1
- isigned = 1
- endian = 0
- wordsize = 1
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(ncat(24,grid % nCells))
- ncat(:,:) = 0
- grid % lu_index % array(:) = 0.0
-
- do jTileStart=1,20401,ny
- jTileEnd = jTileStart + ny - 1
- do iTileStart=1,42001,nx
- iTileEnd = iTileStart + nx - 1
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'/landuse_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-write(0,*) trim(fname)
-
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- iPoint = 1
- do j=1,ny
- do i=1,nx
- lat_pt = (-90.0 + 1.0/240.0) + (jTileStart + j - 2) / 120.0
- lon_pt = (-180.0 + 1.0/240.0) + (iTileStart + i - 2) / 120.0
- lat_pt = lat_pt * pii / 180.0
- lon_pt = lon_pt * pii / 180.0
-
- iPoint = nearest_cell(lat_pt, lon_pt, &
- iPoint, &
- grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
- grid % latCell % array, grid % lonCell % array)
-
- ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
-
- end do
- end do
-
- end do
- end do
-
- do iCell=1, grid % nCells
- grid % lu_index % array(iCell) = 1
- do i=2,24
- if (ncat(i,iCell) > ncat(grid % lu_index % array(iCell),iCell)) then
- grid % lu_index % array(iCell) = i
- end if
- end do
- end do
-
- deallocate(rarray)
- deallocate(ncat)
-
-
- !
- ! Interpolate SOILCAT_TOP
- !
- nx = 1200
- ny = 1200
- nzz = 1
- isigned = 1
- endian = 0
- wordsize = 1
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(ncat(16,grid % nCells))
- ncat(:,:) = 0
- grid % soilcat_top % array(:) = 0.0
-
- do jTileStart=1,20401,ny
- jTileEnd = jTileStart + ny - 1
- do iTileStart=1,42001,nx
- iTileEnd = iTileStart + nx - 1
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'/soiltype_top_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-write(0,*) trim(fname)
-
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- iPoint = 1
- do j=1,ny
- do i=1,nx
- lat_pt = (-90.0 + 1.0/240.0) + (jTileStart + j - 2) / 120.0
- lon_pt = (-180.0 + 1.0/240.0) + (iTileStart + i - 2) / 120.0
- lat_pt = lat_pt * pii / 180.0
- lon_pt = lon_pt * pii / 180.0
-
- iPoint = nearest_cell(lat_pt, lon_pt, &
- iPoint, &
- grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
- grid % latCell % array, grid % lonCell % array)
-
- ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
-
- end do
- end do
-
- end do
- end do
-
- do iCell=1, grid % nCells
- grid % soilcat_top % array(iCell) = 1
- do i=2,16
- if (ncat(i,iCell) > ncat(grid % soilcat_top % array(iCell),iCell)) then
- grid % soilcat_top % array(iCell) = i
- end if
- end do
- end do
-
- deallocate(rarray)
- deallocate(ncat)
-
-
- !
- ! Interpolate SOILCAT_BOT
- !
- nx = 1200
- ny = 1200
- nzz = 1
- isigned = 1
- endian = 0
- wordsize = 1
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(ncat(16,grid % nCells))
- ncat(:,:) = 0
- grid % soilcat_bot % array(:) = 0.0
-
- do jTileStart=1,20401,ny
- jTileEnd = jTileStart + ny - 1
- do iTileStart=1,42001,nx
- iTileEnd = iTileStart + nx - 1
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'/soiltype_bot_30s/',iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-write(0,*) trim(fname)
-
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- iPoint = 1
- do j=1,ny
- do i=1,nx
- lat_pt = (-90.0 + 1.0/240.0) + (jTileStart + j - 2) / 120.0
- lon_pt = (-180.0 + 1.0/240.0) + (iTileStart + i - 2) / 120.0
- lat_pt = lat_pt * pii / 180.0
- lon_pt = lon_pt * pii / 180.0
-
- iPoint = nearest_cell(lat_pt, lon_pt, &
- iPoint, &
- grid % nCells, grid % maxEdges, grid % nEdgesOnCell % array, grid % cellsOnCell % array, &
- grid % latCell % array, grid % lonCell % array)
-
- ncat(int(rarray(i,j,1)),iPoint) = ncat(int(rarray(i,j,1)),iPoint) + 1
-
- end do
- end do
-
- end do
- end do
-
- do iCell=1, grid % nCells
- grid % soilcat_bot % array(iCell) = 1
- do i=2,16
- if (ncat(i,iCell) > ncat(grid % soilcat_bot % array(iCell),iCell)) then
- grid % soilcat_bot % array(iCell) = i
- end if
- end do
- end do
-
- deallocate(rarray)
- deallocate(ncat)
-
-
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! KLUDGE TO FIX SOIL TYPE OVER ANTARCTICA
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- where (grid % lu_index % array == 24) grid % soilcat_top % array = 16
- where (grid % lu_index % array == 24) grid % soilcat_bot % array = 16
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! CORRECT INCONSISTENT SOIL AND LAND USE DATA
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- do iCell = 1,grid % nCells
- if (grid % lu_index % array(iCell) == 16 .or. &
- grid % soilcat_top % array(iCell) == 14 .or. &
- grid % soilcat_bot % array(iCell) == 14) then
- if (grid % lu_index % array(iCell) /= 16) then
- write(0,*) 'Turning lu_index into water at ', iCell
- grid % lu_index % array(iCell) = 16
- end if
- if (grid % soilcat_top % array(iCell) /= 14) then
- write(0,*) 'Turning soilcat_top into water at ', iCell
- grid % soilcat_top % array(iCell) = 14
- end if
- if (grid % soilcat_bot % array(iCell) /= 14) then
- write(0,*) 'Turning soilcat_bot into water at ', iCell
- grid % soilcat_bot % array(iCell) = 14
- end if
- end if
- end do
-
-
- !
- ! Derive LANDMASK
- !
- grid % landmask % array(:) = 0
- do iCell=1, grid % nCells
- if (grid % lu_index % array(iCell) /= 16) grid % landmask % array(iCell) = 1
- end do
-
-
- !
- ! Interpolate SOILTEMP:
- !
- nx = 186
- ny = 186
- nzz = 1
- isigned = 0
- endian = 0
- wordsize = 2
- scalefactor = 0.01
- allocate(rarray(nx,ny,nzz))
- allocate(soiltemp_1deg(-2:363,-2:183))
- grid % soiltemp % array(:) = 0.0
-
- call map_set(PROJ_LATLON, proj, &
- latinc = 1.0_RKIND, &
- loninc = 1.0_RKIND, &
- knowni = 1.0_RKIND, &
- knownj = 1.0_RKIND, &
- lat1 = -89.5_RKIND, &
- lon1 = -179.5_RKIND)
-
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'soiltemp_1deg/',1,'-',180,'.',1,'-',180
-write(0,*) trim(fname)
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- soiltemp_1deg(-2:180,-2:183) = rarray(1:183,1:186,1)
-
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'soiltemp_1deg/',181,'-',360,'.',1,'-',180
-write(0,*) trim(fname)
- call read_geogrid(fname, len_trim(fname), &
- rarray, &
- nx, ny, nzz, &
- isigned, endian, scalefactor, wordsize, istatus)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- soiltemp_1deg(181:363,-2:183) = rarray(4:186,1:186,1)
-
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
- interp_list(5) = 0
-
- do iCell=1,grid%nCells
-
- if (grid % landmask % array(iCell) == 1) then
- 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.5) then
- lon = lon + 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- else if (x >= 360.5) then
- lon = lon - 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- end if
-if (y < 1.0) y = 1.0
-if (y > 179.0) y = 179.0
-! grid % soiltemp % array(iCell) = interp_sequence(x, y, 1, soiltemp_1deg, 1, 360, 1, 180, 1, 1, -1.e30_RKIND, interp_list, 1)
- grid % soiltemp % array(iCell) = interp_sequence(x, y, 1, soiltemp_1deg, -2, 363, -2, 183, 1, 1, 0.0_RKIND, interp_list, 1)
- else
- grid % soiltemp % array(iCell) = 0.0
- end if
-
- end do
-
- deallocate(rarray)
- deallocate(soiltemp_1deg)
-
-
- !
- ! Interpolate SNOALB
- !
- nx = 186
- ny = 186
- nzz = 1
- isigned = 0
- endian = 0
- wordsize = 1
- scalefactor = 1.0
- allocate(rarray(nx,ny,nzz))
- allocate(maxsnowalb(-2:363,-2:183))
- grid % snoalb % array(:) = 0.0
-
- call map_set(PROJ_LATLON, proj, &
- latinc = 1.0_RKIND, &
- loninc = 1.0_RKIND, &
- knowni = 1.0_RKIND, &
- knownj = 1.0_RKIND, &
- lat1 = -89.5_RKIND, &
- lon1 = -179.5_RKIND)
-
- 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)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- maxsnowalb(-2:180,-2:183) = rarray(1:183,1:186,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)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- maxsnowalb(181:363,-2:183) = rarray(4:186,1:186,1)
-
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
- interp_list(5) = 0
-
- do iCell=1,grid%nCells
-
- if (grid % landmask % array(iCell) == 1) then
- 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.5) then
- lon = lon + 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- else if (x >= 360.5) then
- lon = lon - 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- end if
-if (y < 1.0) y = 1.0
-if (y > 179.0) y = 179.0
-! grid % snoalb % array(iCell) = interp_sequence(x, y, 1, maxsnowalb, 1, 360, 1, 180, 1, 1, -1.e30_RKIND, interp_list, 1)
- grid % snoalb % array(iCell) = interp_sequence(x, y, 1, maxsnowalb, -2, 363, -2, 183, 1, 1, 0.0_RKIND, interp_list, 1)
- else
- grid % snoalb % array(iCell) = 0.0
- end if
-
- 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(-2:2503,-2:1253,12))
-! grid % vegfra % array(:) = 0.0
- grid % greenfrac % array(:,:) = 0.0
-
- call map_set(PROJ_LATLON, proj, &
- latinc = 0.144_RKIND, &
- loninc = 0.144_RKIND, &
- knowni = 1.0_RKIND, &
- knownj = 1.0_RKIND, &
- lat1 = -89.928_RKIND, &
- lon1 = -179.928_RKIND)
-
- 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)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- vegfra(-2:1250,-2:1253,1:12) = rarray(1:1253,1:1256,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)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- vegfra(1251:2503,-2:1253,1:12) = rarray(4:1256,1:1256,1:12)
-
- do iCell=1,grid%nCells
- if (grid % landmask % array(iCell) == 1) then
- 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.5) then
- lon = lon + 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- else if (x >= 2500.5) then
- lon = lon - 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- end if
-if (y < 1.0) y = 1.0
-if (y > 1249.0) y = 1249.0
- do k=1,12
- grid % greenfrac % array(k,iCell) = interp_sequence(x, y, k, vegfra, -2, 2503, -2, 1253, 1, 12, -1.e30_RKIND, interp_list, 1)
- end do
- else
- grid % greenfrac % array(:,iCell) = 0.0
- end if
- grid % shdmin % array(iCell) = minval(grid % greenfrac % array(:,iCell))
- grid % shdmax % array(iCell) = maxval(grid % greenfrac % array(:,iCell))
-
- end do
-
- 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(-2:2503,-2:1253,12))
- grid % albedo12m % array(:,:) = 0.0
-
- call map_set(PROJ_LATLON, proj, &
- latinc = 0.144_RKIND, &
- loninc = 0.144_RKIND, &
- knowni = 1.0_RKIND, &
- knownj = 1.0_RKIND, &
- lat1 = -89.928_RKIND, &
- lon1 = -179.928_RKIND)
-
- 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)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- vegfra(-2:1250,-2:1253,1:12) = rarray(1:1253,1:1256,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)
- call init_atm_check_read_error(istatus, fname, dminfo)
-
- vegfra(1251:2503,-2:1253,1:12) = rarray(4:1256,1:1256,1:12)
-
- do iCell=1,grid%nCells
- if (grid % landmask % array(iCell) == 1) then
- 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.5) then
- lon = lon + 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- else if (x >= 2500.5) then
- lon = lon - 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- end if
-if (y < 1.0) y = 1.0
-if (y > 1249.0) y = 1249.0
- do k=1,12
- grid % albedo12m % array(k,iCell) = interp_sequence(x, y, k, vegfra, -2, 2503, -2, 1253, 1, 12, 0.0_RKIND, interp_list, 1)
- end do
- else
- grid % albedo12m % array(:,iCell) = 8.0
- end if
- end do
-
- deallocate(rarray)
- deallocate(vegfra)
-
-
- end if ! config_static_interp
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN ADOPT GFS TERRAIN HEIGHT
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -5977,48 +5358,6 @@
end subroutine init_atm_test_case_resting_atmosphere
- integer function nearest_cell(target_lat, target_lon, &
- start_cell, &
- nCells, maxEdges, nEdgesOnCell, cellsOnCell, latCell, lonCell)
-
- implicit none
-
- real (kind=RKIND), intent(in) :: target_lat, target_lon
- integer, intent(in) :: start_cell
- integer, intent(in) :: nCells, maxEdges
- integer, dimension(nCells), intent(in) :: nEdgesOnCell
- integer, dimension(maxEdges,nCells), intent(in) :: cellsOnCell
- real (kind=RKIND), dimension(nCells), intent(in) :: latCell, lonCell
-
- integer :: i
- integer :: iCell
- integer :: current_cell
- real (kind=RKIND) :: current_distance, d
- real (kind=RKIND) :: nearest_distance
-
- nearest_cell = start_cell
- current_cell = -1
-
- do while (nearest_cell /= current_cell)
- current_cell = nearest_cell
- current_distance = sphere_distance(latCell(current_cell), lonCell(current_cell), target_lat, target_lon, 1.0_RKIND)
- nearest_cell = current_cell
- nearest_distance = current_distance
- do i = 1, nEdgesOnCell(current_cell)
- iCell = cellsOnCell(i,current_cell)
- if (iCell <= nCells) then
- d = sphere_distance(latCell(iCell), lonCell(iCell), target_lat, target_lon, 1.0_RKIND)
- if (d < nearest_distance) then
- nearest_cell = iCell
- nearest_distance = d
- end if
- end if
- end do
- end do
-
- end function nearest_cell
-
-
integer function nearest_edge(target_lat, target_lon, &
start_edge, &
nCells, nEdges, maxEdges, nEdgesOnCell, edgesOnCell, cellsOnEdge, latCell, lonCell, latEdge, lonEdge)
@@ -6160,44 +5499,8 @@
end function vertical_interp
- subroutine init_atm_check_read_error(istatus, fname, dminfo)
-
- implicit none
-
- integer, intent(in) :: istatus
- character (len=*), intent(in) :: fname
- type (dm_info), intent(in) :: dminfo
-
- if (istatus /= 0) then
- write(0,*) 'ERROR: Could not read file '//trim(fname)
- call mpas_dmpar_abort(dminfo)
- end if
-
- end subroutine init_atm_check_read_error
-
-
!----------------------------------------------------------------------------------------------------------
- real (kind=RKIND) function sphere_distance(lat1, lon1, lat2, lon2, radius)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Compute the great-circle distance between (lat1, lon1) and (lat2, lon2) on a
- ! sphere with given radius.
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- implicit none
-
- real (kind=RKIND), intent(in) :: lat1, lon1, lat2, lon2, radius
-
- real (kind=RKIND) :: arg1
-
- arg1 = sqrt( sin(0.5*(lat2-lat1))**2 + &
- cos(lat1)*cos(lat2)*sin(0.5*(lon2-lon1))**2 )
- sphere_distance = 2.*radius*asin(arg1)
-
- end function sphere_distance
-
-!--------------------------------------------------------------------
-
real (kind=RKIND) function env_qv( z, temperature, pressure, rh_max )
implicit none
</font>
</pre>