<p><b>laura@ucar.edu</b> 2013-01-30 14:35:31 -0700 (Wed, 30 Jan 2013)</p><p>Moved the interpolation of the static fields from subroutine init_atm_test_case_gfs (in mpas_init_atm_test_cases.F) to subroutine init_atm_static. Added the interpolation of the static statistical fields needed for the parameterization of the gravity wave drag over orography (con,oa1,oa2,oa3,oa4,ol1,ol2,ol3,and var2d). The sourcecode is only available for uniform x1 meshes (10242,40962,163842).<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_init_nhyd_atmos/mpas_init_atm_static.F
===================================================================
--- branches/atmos_physics/src/core_init_nhyd_atmos/mpas_init_atm_static.F 2013-01-30 21:11:54 UTC (rev 2393)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/mpas_init_atm_static.F 2013-01-30 21:35:31 UTC (rev 2394)
@@ -7,6 +7,8 @@
use init_atm_hinterp
use init_atm_llxy
+ use mpas_atmphys_utilities
+
implicit none
private
public:: init_atm_static, &
@@ -649,7 +651,9 @@
type(proj_info):: proj
type(dm_info),pointer :: dminfo
+ character(len=StrKIND):: mess
character(len=StrKIND):: fname
+ character(len=StrKIND):: dir_gwdo,fname_g
integer:: nx,ny,nz
integer:: endian,isigned,istatus,wordsize
@@ -664,85 +668,181 @@
real(kind=RKIND):: lat,lon,x,y
real(kind=RKIND):: lat_pt,lon_pt
real(kind=RKIND):: dx,dy,known_lat,known_lon,known_x,known_y
+ real(kind=RKIND):: minMeshD,maxMeshD
+ real(kind=RKIND):: mindcEdge,maxdcEdge
real(kind=RKIND),dimension(:,:),allocatable:: xarray
!--------------------------------------------------------------------------------------------------
write(0,*)
write(0,*) '--- enter subroutine init_atm_static_orogwd:'
+!goto 100
!
! Interpolate VARSSO:
-!mesh%varsso%array(:) = 0.0_RKIND
-!nx = 600
-!ny = 600
-!nz = 1
-!isigned = 0
-!endian = 0
-!wordsize = 4
-!scalefactor = 1.0
+ mesh%varsso%array(:) = 0.0_RKIND
+ nx = 600
+ ny = 600
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 4
+ scalefactor = 1.0
-!dx = 0.00833333
-!dy = 0.00833333
-!known_x = 1.0
-!known_y = 1.0
-!known_lat = -59.99583
-!known_lon = -179.99583
+ dx = 0.00833333
+ dy = 0.00833333
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -59.99583
+ known_lon = -179.99583
-!allocate(rarray(nx,ny,nz))
-!allocate(nhs(mesh%nCells))
-!do jTileStart = 1,13801,ny
-! jTileEnd = jTileStart + ny - 1
+ allocate(rarray(nx,ny,nz))
+ allocate(nhs(mesh%nCells))
+ nhs(:) = 0
+ rarray(:,:,:) = 0._RKIND
+ do jTileStart = 1,13801,ny
+ jTileEnd = jTileStart + ny - 1
-! do iTileStart = 1,42601,nx
-! iTileEnd = iTileStart + nx -1
-! write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'varsso/', &
-! iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
-! write(0,*) trim(fname)
+ do iTileStart = 1,42601,nx
+ iTileEnd = iTileStart + nx -1
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'varsso/', &
+ iTileStart,'-',iTileEnd,'.',jTileStart,'-',jTileEnd
+ write(0,*) trim(fname)
-! call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
-! scalefactor,wordsize,istatus)
-! call init_atm_check_read_error(istatus,fname,dminfo)
+ call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,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 = known_lat + (jTileStart + j - 2) * dy
-! lon_pt = known_lon + (iTileStart + i - 2) * dx
-! lat_pt = lat_pt * PI / 180.0
-! lon_pt = lon_pt * PI / 180.0
+ iPoint = 1
+ do j = 1,ny
+ do i = 1,nx
+ lat_pt = known_lat + (jTileStart + j - 2) * dy
+ lon_pt = known_lon + (iTileStart + i - 2) * dx
+ lat_pt = lat_pt * PI / 180.0
+ lon_pt = lon_pt * PI / 180.0
-! iPoint = nearest_cell(lat_pt,lon_pt,iPoint,mesh%nCells,mesh%maxEdges, &
-! mesh%nEdgesOnCell%array,mesh%cellsOnCell%array, &
-! mesh%latCell%array,mesh%lonCell%array)
-! mesh%varsso%array(iPoint) = mesh%varsso%array(iPoint) + rarray(i,j,1)
-! nhs(iPoint) = nhs(iPoint) + 1
-! enddo
-! enddo
+ iPoint = nearest_cell(lat_pt,lon_pt,iPoint,mesh%nCells,mesh%maxEdges, &
+ mesh%nEdgesOnCell%array,mesh%cellsOnCell%array, &
+ mesh%latCell%array,mesh%lonCell%array)
+ mesh%varsso%array(iPoint) = mesh%varsso%array(iPoint) + rarray(i,j,1)
+ nhs(iPoint) = nhs(iPoint) + 1
+ enddo
+ enddo
-! enddo
-!enddo
+ enddo
+ enddo
-!do iCell = 1,mesh%nCells
-! mesh%varsso%array(iCell) = mesh%varsso%array(iCell) / real(nhs(iCell))
-!enddo
-!deallocate(rarray)
-!deallocate(nhs)
-!write(0,*) '--- end interpolate VARSSO'
+ do iCell = 1,mesh%nCells
+ if(nhs(iCell) .gt. 0) &
+ mesh%varsso%array(iCell) = mesh%varsso%array(iCell) / real(nhs(iCell))
+ enddo
+ deallocate(rarray)
+ deallocate(nhs)
+ write(0,*) '--- end interpolate VARSSO'
+ 100 continue
+!... statistic fields needed for the parameterization of gravity wavwe drag over orography. The
+!input directory depends on the mesh resolution, and the mesh must be a uniform mesh.
+ minMeshD = minval(mesh%meshDensity%array(1:mesh%nCells))
+ maxMeshD = maxval(mesh%meshDensity%array(1:mesh%nCells))
+ mindcEdge = minval(mesh%dcEdge%array(1:mesh%nEdges))
+ maxdcEdge = maxval(mesh%dcEdge%array(1:mesh%nEdges))
+
+ write(0,*)
+ write(0,*) 'BEGIN INTERPOLATION OF STATISTICAL FIELDS FOR GRAVITY WAVE DRAG OVER OROGRAPHY'
+ write(0,*) 'min MeshD =', minMeshD
+ write(0,*) 'max MeshD =', maxMeshD
+ write(0,*) 'min dcEdge =', mindcEdge
+ write(0,*) 'max dcEdge =', maxdcEdge
+
+ dir_gwdo = ' '
+ if(minMeshD == 1.0_RKIND .and. maxMeshD == 1.0_RKIND) then
+ !... uniform 10242 mesh:
+ if(mindcEdge .ge. 200000._RKIND .and. maxdcEdge .lt. 260000._RKIND) then
+ dir_gwdo = 'orogwd_2deg'
+ elseif(mindcEdge .ge. 90000._RKIND .and. maxdcEdge .lt. 150000_RKIND) then
+ dir_gwdo = 'orogwd_1deg'
+ elseif(mindcEdge .ge. 40000._RKIND .and. maxdcEdge .lt. 70000._RKIND) then
+ dir_gwdo = 'orogwd_30m'
+ else
+ write(0,*)
+ write(mess,*) 'GWDO: Interpolation not available. The initialization will abort'
+ call physics_error_fatal(mess)
+ endif
+ else
+ write(0,*)
+ write(mess,*) 'GWDO: The input mesh must be a uniform mesh. The initialization will abort'
+ call physics_error_fatal(mess)
+ endif
+ write(0,*) 'dir_gwdo = ', trim(dir_gwdo)
+ write(0,*)
+
!
! Interpolate CON:
!
mesh%con%array(:) = 0.0_RKIND
- nx = 2160
- ny = 1080
- nz = 1
- isigned = 0
- endian = 0
- wordsize = 2
- scalefactor = 0.025
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'orogwd_10m/con/', &
- 1,'-',2160,'.',1,'-',1080
+ con_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.025
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.025
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.025
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.025
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select con_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/con/',1,'-',nx,'.',1,'-',ny
write(0,*) trim(fname)
allocate(xarray(nx,ny))
@@ -750,14 +850,8 @@
call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
scalefactor,wordsize,istatus)
call init_atm_check_read_error(istatus,fname,dminfo)
- xarray(1:nx,1:ny) = rarray(1:nx,1:ny)
+ xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
- dx = 0.16666667
- dy = 0.16666667
- known_x = 1.0
- known_y = 1.0
- known_lat = -89.916667
- known_lon = 0.0833333
call map_set(PROJ_LATLON, proj, &
latinc = dy, &
loninc = dx, &
@@ -766,10 +860,10 @@
lat1 = known_lat, &
lon1 = known_lon)
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
interp_list(5) = 0
do iCell = 1,mesh%nCells
@@ -777,7 +871,7 @@
lat = mesh % latCell % array(iCell) * DEG_PER_RAD
lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- mesh % con % array(iCell) = interp_sequence(x,y,1,xarray,1,2160,1,1080,1,1, &
+ mesh % con % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
0.0_RKIND,interp_list,1)
endif
enddo
@@ -789,16 +883,69 @@
! Interpolate OA1:
!
mesh%oa1%array(:) = 0.0_RKIND
- nx = 2160
- ny = 1080
- nz = 1
- isigned = 1
- endian = 0
- wordsize = 2
- scalefactor = 0.0001
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'orogwd_10m/oa1/', &
- 1,'-',2160,'.',1,'-',1080
+ oa1_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select oa1_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/oa1/',1,'-',nx,'.',1,'-',ny
write(0,*) trim(fname)
allocate(xarray(nx,ny))
@@ -808,12 +955,6 @@
call init_atm_check_read_error(istatus,fname,dminfo)
xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
- dx = 0.16666667
- dy = 0.16666667
- known_x = 1.0
- known_y = 1.0
- known_lat = -89.916667
- known_lon = 0.0833333
call map_set(PROJ_LATLON, proj, &
latinc = dy, &
loninc = dx, &
@@ -822,10 +963,10 @@
lat1 = known_lat, &
lon1 = known_lon)
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
interp_list(5) = 0
do iCell = 1,mesh%nCells
@@ -833,7 +974,7 @@
lat = mesh % latCell % array(iCell) * DEG_PER_RAD
lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- mesh % oa1 % array(iCell) = interp_sequence(x,y,1,xarray,1,2160,1,1080,1,1, &
+ mesh % oa1 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
0.0_RKIND,interp_list,1)
endif
enddo
@@ -844,16 +985,69 @@
!
! Interpolate OA2:
mesh%oa2%array(:) = 0.0_RKIND
- nx = 2160
- ny = 1080
- nz = 1
- isigned = 1
- endian = 0
- wordsize = 2
- scalefactor = 0.0001
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'orogwd_10m/oa2/', &
- 1,'-',2160,'.',1,'-',1080
+ oa2_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select oa2_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/oa2/',1,'-',nx,'.',1,'-',ny
write(0,*) trim(fname)
allocate(xarray(nx,ny))
@@ -863,12 +1057,6 @@
call init_atm_check_read_error(istatus,fname,dminfo)
xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
- dx = 0.16666667
- dy = 0.16666667
- known_x = 1.0
- known_y = 1.0
- known_lat = -89.916667
- known_lon = 0.0833333
call map_set(PROJ_LATLON, proj, &
latinc = dy, &
loninc = dx, &
@@ -877,10 +1065,10 @@
lat1 = known_lat, &
lon1 = known_lon)
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
interp_list(5) = 0
do iCell = 1,mesh%nCells
@@ -888,7 +1076,7 @@
lat = mesh % latCell % array(iCell) * DEG_PER_RAD
lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- mesh % oa2 % array(iCell) = interp_sequence(x,y,1,xarray,1,2160,1,1080,1,1, &
+ mesh % oa2 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
0.0_RKIND,interp_list,1)
endif
enddo
@@ -900,16 +1088,69 @@
! Interpolate OA3:
!
mesh%oa3%array(:) = 0.0_RKIND
- nx = 2160
- ny = 1080
- nz = 1
- isigned = 1
- endian = 0
- wordsize = 2
- scalefactor = 0.0001
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'orogwd_10m/oa3/', &
- 1,'-',2160,'.',1,'-',1080
+ oa3_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select oa3_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/oa3/',1,'-',nx,'.',1,'-',ny
write(0,*) trim(fname)
allocate(xarray(nx,ny))
@@ -919,12 +1160,6 @@
call init_atm_check_read_error(istatus,fname,dminfo)
xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
- dx = 0.16666667
- dy = 0.16666667
- known_x = 1.0
- known_y = 1.0
- known_lat = -89.916667
- known_lon = 0.0833333
call map_set(PROJ_LATLON, proj, &
latinc = dy, &
loninc = dx, &
@@ -933,10 +1168,10 @@
lat1 = known_lat, &
lon1 = known_lon)
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
interp_list(5) = 0
do iCell = 1,mesh%nCells
@@ -944,7 +1179,7 @@
lat = mesh % latCell % array(iCell) * DEG_PER_RAD
lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- mesh % oa3 % array(iCell) = interp_sequence(x,y,1,xarray,1,2160,1,1080,1,1, &
+ mesh % oa3 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
0.0_RKIND,interp_list,1)
endif
enddo
@@ -956,16 +1191,69 @@
! Interpolate OA4:
!
mesh%oa4%array(:) = 0.0_RKIND
- nx = 2160
- ny = 1080
- nz = 1
- isigned = 1
- endian = 0
- wordsize = 2
- scalefactor = 0.0001
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'orogwd_10m/oa4/', &
- 1,'-',2160,'.',1,'-',1080
+ oa4_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 1
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select oa4_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/oa4/',1,'-',nx,'.',1,'-',ny
write(0,*) trim(fname)
allocate(xarray(nx,ny))
@@ -975,12 +1263,6 @@
call init_atm_check_read_error(istatus,fname,dminfo)
xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
- dx = 0.16666667
- dy = 0.16666667
- known_x = 1.0
- known_y = 1.0
- known_lat = -89.916667
- known_lon = 0.0833333
call map_set(PROJ_LATLON, proj, &
latinc = dy, &
loninc = dx, &
@@ -989,10 +1271,10 @@
lat1 = known_lat, &
lon1 = known_lon)
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
interp_list(5) = 0
do iCell = 1,mesh%nCells
@@ -1000,7 +1282,7 @@
lat = mesh % latCell % array(iCell) * DEG_PER_RAD
lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- mesh % oa4 % array(iCell) = interp_sequence(x,y,1,xarray,1,2160,1,1080,1,1, &
+ mesh % oa4 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
0.0_RKIND,interp_list,1)
endif
enddo
@@ -1012,16 +1294,69 @@
! Interpolate OL1:
!
mesh%ol1%array(:) = 0.0_RKIND
- nx = 2160
- ny = 1080
- nz = 1
- isigned = 0
- endian = 0
- wordsize = 2
- scalefactor = 0.0001
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'orogwd_10m/ol1/', &
- 1,'-',2160,'.',1,'-',1080
+ ol1_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select ol1_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/ol1/',1,'-',nx,'.',1,'-',ny
write(0,*) trim(fname)
allocate(xarray(nx,ny))
@@ -1031,12 +1366,6 @@
call init_atm_check_read_error(istatus,fname,dminfo)
xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
- dx = 0.16666667
- dy = 0.16666667
- known_x = 1.0
- known_y = 1.0
- known_lat = -89.916667
- known_lon = 0.0833333
call map_set(PROJ_LATLON, proj, &
latinc = dy, &
loninc = dx, &
@@ -1045,10 +1374,10 @@
lat1 = known_lat, &
lon1 = known_lon)
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
interp_list(5) = 0
do iCell = 1,mesh%nCells
@@ -1056,7 +1385,7 @@
lat = mesh % latCell % array(iCell) * DEG_PER_RAD
lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- mesh % ol1 % array(iCell) = interp_sequence(x,y,1,xarray,1,2160,1,1080,1,1, &
+ mesh % ol1 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
0.0_RKIND,interp_list,1)
endif
enddo
@@ -1068,16 +1397,69 @@
! Interpolate OL2:
!
mesh%ol2%array(:) = 0.0_RKIND
- nx = 2160
- ny = 1080
- nz = 1
- isigned = 0
- endian = 0
- wordsize = 2
- scalefactor = 0.0001
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'orogwd_10m/ol2/', &
- 1,'-',2160,'.',1,'-',1080
+ ol2_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select ol2_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/ol2/',1,'-',nx,'.',1,'-',ny
write(0,*) trim(fname)
allocate(xarray(nx,ny))
@@ -1087,12 +1469,6 @@
call init_atm_check_read_error(istatus,fname,dminfo)
xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
- dx = 0.16666667
- dy = 0.16666667
- known_x = 1.0
- known_y = 1.0
- known_lat = -89.916667
- known_lon = 0.0833333
call map_set(PROJ_LATLON, proj, &
latinc = dy, &
loninc = dx, &
@@ -1101,10 +1477,10 @@
lat1 = known_lat, &
lon1 = known_lon)
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
interp_list(5) = 0
do iCell = 1,mesh%nCells
@@ -1112,7 +1488,7 @@
lat = mesh % latCell % array(iCell) * DEG_PER_RAD
lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- mesh % ol2 % array(iCell) = interp_sequence(x,y,1,xarray,1,2160,1,1080,1,1, &
+ mesh % ol2 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
0.0_RKIND,interp_list,1)
endif
enddo
@@ -1124,16 +1500,69 @@
! Interpolate OL3:
!
mesh%ol3%array(:) = 0.0_RKIND
- nx = 2160
- ny = 1080
- nz = 1
- isigned = 0
- endian = 0
- wordsize = 2
- scalefactor = 0.0001
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'orogwd_10m/ol3/', &
- 1,'-',2160,'.',1,'-',1080
+ ol3_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select ol3_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/ol3/',1,'-',nx,'.',1,'-',ny
write(0,*) trim(fname)
allocate(xarray(nx,ny))
@@ -1143,12 +1572,6 @@
call init_atm_check_read_error(istatus,fname,dminfo)
xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
- dx = 0.16666667
- dy = 0.16666667
- known_x = 1.0
- known_y = 1.0
- known_lat = -89.916667
- known_lon = 0.0833333
call map_set(PROJ_LATLON, proj, &
latinc = dy, &
loninc = dx, &
@@ -1157,10 +1580,10 @@
lat1 = known_lat, &
lon1 = known_lon)
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
interp_list(5) = 0
do iCell = 1,mesh%nCells
@@ -1168,7 +1591,7 @@
lat = mesh % latCell % array(iCell) * DEG_PER_RAD
lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- mesh % ol3 % array(iCell) = interp_sequence(x,y,1,xarray,1,2160,1,1080,1,1, &
+ mesh % ol3 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
0.0_RKIND,interp_list,1)
endif
enddo
@@ -1180,16 +1603,69 @@
! Interpolate OL4:
!
mesh%ol4%array(:) = 0.0_RKIND
- nx = 2160
- ny = 1080
- nz = 1
- isigned = 0
- endian = 0
- wordsize = 2
- scalefactor = 0.0001
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'orogwd_10m/ol4/', &
- 1,'-',2160,'.',1,'-',1080
+ ol4_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.0001
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select ol4_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/ol4/',1,'-',nx,'.',1,'-',ny
write(0,*) trim(fname)
allocate(xarray(nx,ny))
@@ -1199,12 +1675,6 @@
call init_atm_check_read_error(istatus,fname,dminfo)
xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
- dx = 0.16666667
- dy = 0.16666667
- known_x = 1.0
- known_y = 1.0
- known_lat = -89.916667
- known_lon = 0.0833333
call map_set(PROJ_LATLON, proj, &
latinc = dy, &
loninc = dx, &
@@ -1213,10 +1683,10 @@
lat1 = known_lat, &
lon1 = known_lon)
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
interp_list(5) = 0
do iCell = 1,mesh%nCells
@@ -1224,7 +1694,7 @@
lat = mesh % latCell % array(iCell) * DEG_PER_RAD
lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- mesh % ol4 % array(iCell) = interp_sequence(x,y,1,xarray,1,2160,1,1080,1,1, &
+ mesh % ol4 % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
0.0_RKIND,interp_list,1)
endif
enddo
@@ -1236,18 +1706,72 @@
! Interpolate VAR2D:
!
mesh%var2d%array(:) = 0.0_RKIND
- nx = 2160
- ny = 1080
- nz = 1
- isigned = 0
- endian = 0
- wordsize = 2
- scalefactor = 0.02
- write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') trim(config_geog_data_path)//'orogwd_10m/var/', &
- 1,'-',2160,'.',1,'-',1080
+ var2d_select: select case(dir_gwdo)
+ case("orogwd_2deg")
+ nx = 180
+ ny = 90
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 4
+ scalefactor = 0.02
+ dx = 2.0
+ dy = 2.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.0
+ known_lon = 1.0
+ case("orogwd_1deg")
+ nx = 360
+ ny = 180
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 4
+ scalefactor = 0.02
+ dx = 1.0
+ dy = 1.0
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.5
+ known_lon = 0.5
+ case("orogwd_30m")
+ nx = 720
+ ny = 360
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 4
+ scalefactor = 0.02
+ dx = 0.5
+ dy = 0.5
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.75
+ known_lon = 0.25
+ case("orogwd_10m")
+ nx = 2160
+ ny = 1080
+ nz = 1
+ isigned = 0
+ endian = 0
+ wordsize = 2
+ scalefactor = 0.02
+ dx = 0.16666667
+ dy = 0.16666667
+ known_x = 1.0
+ known_y = 1.0
+ known_lat = -89.916667
+ known_lon = 0.0833333
+ case default
+ end select var2d_select
+
+ write(fname,'(a,i5.5,a1,i5.5,a1,i5.5,a1,i5.5)') &
+ trim(config_geog_data_path)//trim(dir_gwdo)//'/var/',1,'-',nx,'.',1,'-',ny
write(0,*) trim(fname)
+
allocate(xarray(nx,ny))
allocate(rarray(nx,ny,nz))
call read_geogrid(fname,len_trim(fname),rarray,nx,ny,nz,isigned,endian, &
@@ -1255,12 +1779,6 @@
call init_atm_check_read_error(istatus,fname,dminfo)
xarray(1:nx,1:ny) = rarray(1:nx,1:ny,1)
- dx = 0.16666667
- dy = 0.16666667
- known_x = 1.0
- known_y = 1.0
- known_lat = -89.916667
- known_lon = 0.0833333
call map_set(PROJ_LATLON, proj, &
latinc = dy, &
loninc = dx, &
@@ -1269,10 +1787,10 @@
lat1 = known_lat, &
lon1 = known_lon)
- interp_list(1) = FOUR_POINT
- interp_list(2) = W_AVERAGE4
- interp_list(3) = W_AVERAGE16
- interp_list(4) = SEARCH
+ interp_list(1) = AVERAGE4
+ interp_list(2) = AVERAGE4
+ interp_list(3) = AVERAGE4
+ interp_list(4) = AVERAGE4
interp_list(5) = 0
do iCell = 1,mesh%nCells
@@ -1280,7 +1798,7 @@
lat = mesh % latCell % array(iCell) * DEG_PER_RAD
lon = mesh % lonCell % array(iCell) * DEG_PER_RAD
call latlon_to_ij(proj, lat, lon, x, y)
- mesh % var2d % array(iCell) = interp_sequence(x,y,1,xarray,1,2160,1,1080,1,1, &
+ mesh % var2d % array(iCell) = interp_sequence(x,y,1,xarray,1,nx,1,ny,1,1, &
0.0_RKIND,interp_list,1)
endif
enddo
</font>
</pre>