<p><b>laura@ucar.edu</b> 2012-08-23 12:15:16 -0600 (Thu, 23 Aug 2012)</p><p>changed calls to physics_initialize_real and init_atm_test_case_sfc. In subroutine init_atm_test_case_gfs, changed allocation of maskslab and rslab to -2:field % nx+3, and changed call to function interp_sequence accordingly. In subroutine init_atm_test_case_gfs, added limit to xice between 0. and 1. Removed subroutine init_atm_test_case_sfc (moved to mpas_init_atm_surface.F).<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 2012-08-23 17:29:19 UTC (rev 2117)
+++ branches/atmos_physics/src/core_init_nhyd_atmos/mpas_init_atm_test_cases.F 2012-08-23 18:15:16 UTC (rev 2118)
@@ -10,6 +10,7 @@
use mpas_RBF_interpolation
use mpas_vector_reconstruction
use mpas_timer
+ use mpas_init_atm_surface
! Added only clause to keep xlf90 from getting confused from the overloaded abs intrinsic in mpas_timekeeping
use mpas_timekeeping !, only: MPAS_Time_type, MPAS_TimeInterval_type, MPAS_Clock_type, &
@@ -107,7 +108,7 @@
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)
+ if (config_met_interp) call physics_initialize_real(block_ptr % mesh, block_ptr % fg, domain % dminfo)
block_ptr => block_ptr % next
end do
@@ -116,8 +117,7 @@
write(0,*) 'real-data surface (SST) update test case '
block_ptr => domain % blocklist
do while (associated(block_ptr))
- call init_atm_test_case_sfc(domain, domain % dminfo, block_ptr % mesh, block_ptr % fg, block_ptr % state % time_levs(1) % state, &
- block_ptr % diag, config_test_case, block_ptr % parinfo)
+ call init_atm_test_case_sfc(domain, domain % dminfo, block_ptr % mesh,block_ptr % fg, block_ptr % state % time_levs(1) % state)
block_ptr => block_ptr % next
end do
@@ -3379,7 +3379,7 @@
do while (istatus == 0)
if (index(field % field, 'LANDSEA') /= 0) then
- allocate(maskslab(-3:field % nx+3, field % ny))
+ allocate(maskslab(-2:field % nx+3, field % ny))
maskslab(1:field % nx, 1:field % ny) = field % slab(1:field % nx, 1:field % ny)
maskslab(0, 1:field % ny) = field % slab(field % nx, 1:field % ny)
maskslab(-1, 1:field % ny) = field % slab(field % nx-1, 1:field % ny)
@@ -3824,7 +3824,7 @@
ndims = 1
end if
- allocate(rslab(-3:field % nx+3, field % ny))
+ allocate(rslab(-2:field % nx+3, field % ny))
rslab(1:field % nx, 1:field % ny) = field % slab(1:field % nx, 1:field % ny)
rslab(0, 1:field % ny) = field % slab(field % nx, 1:field % ny)
rslab(-1, 1:field % ny) = field % slab(field % nx-1, 1:field % ny)
@@ -3843,9 +3843,9 @@
call latlon_to_ij(proj, lat, lon, x, y)
end if
if (ndims == 1) then
- destField1d(i) = interp_sequence(x, y, 1, rslab, -3, field % nx + 3, 1, field % ny, 1, 1, msgval, interp_list, 1, maskval=maskval, mask_array=maskslab)
+ destField1d(i) = interp_sequence(x, y, 1, rslab, -2, field % nx + 3, 1, field % ny, 1, 1, msgval, interp_list, 1, maskval=maskval, mask_array=maskslab)
else if (ndims == 2) then
- destField2d(k,i) = interp_sequence(x, y, 1, rslab, -3, field % nx + 3, 1, field % ny, 1, 1, msgval, interp_list, 1, maskval=maskval, mask_array=maskslab)
+ destField2d(k,i) = interp_sequence(x, y, 1, rslab, -2, field % nx + 3, 1, field % ny, 1, 1, msgval, interp_list, 1, maskval=maskval, mask_array=maskslab)
end if
else
if (ndims == 1) then
@@ -4001,6 +4001,12 @@
! Freeze really cold ocean
where (fg % sst % array < 271.0 .and. grid % landmask % array == 0) fg % xice % array = 1.0
+ ! Limit XICE to values between 0 and 1. Although the input meteorological field is between 0.
+ ! and 1., interpolation to the MPAS grid can yield values of XiCE less than 0. and greater
+ ! than 1.:
+ where (fg % xice % array < 0._RKIND) fg % xice % array = 0._RKIND
+ where (fg % xice % array > 1._RKIND) fg % xice % array = 1._RKIND
+
! Set SEAICE (0/1 flag) based on XICE (fractional ice coverage)
fg % seaice % array(:) = 0.0
where (fg % xice % array >= 0.5) fg % seaice % array = 1.0
@@ -4337,242 +4343,6 @@
end subroutine init_atm_test_case_gfs
- subroutine init_atm_test_case_sfc(domain, dminfo, grid, fg, state, diag, test_case, parinfo)
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! Real-data test case using SST data
- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
- use mpas_dmpar
- use mpas_io_output
- use init_atm_read_met
- use init_atm_llxy
- use init_atm_hinterp
-
- implicit none
-
- type (domain_type), intent(inout) :: domain
- type (dm_info), intent(in) :: dminfo
- type (mesh_type), intent(inout) :: grid
- type (fg_type), intent(inout) :: fg
- type (state_type), intent(inout) :: state
- type (diag_type), intent(inout) :: diag
- integer, intent(in) :: test_case
- type (parallel_info), pointer :: parinfo
-
- integer :: istatus
- integer :: iCell, i, j
- type (met_data) :: field
- type (proj_info) :: proj
- real (kind=RKIND) :: lat, lon, x, y
- integer, dimension(5) :: interp_list
- real (kind=RKIND), allocatable, dimension(:,:) :: slab_r8
- type (io_output_object) :: sfc_update_obj
- type (MPAS_Clock_type) :: fg_clock
- type (MPAS_Time_type) :: start_time, stop_time, curr_time
- type (MPAS_TimeInterval_type) :: fg_interval
- character (len=StrKIND) :: timeString
-
-
- ! Set interpolation sequence to be used for SST/SKINTEMP field
- interp_list(1) = FOUR_POINT
- interp_list(2) = SEARCH
- interp_list(3) = 0
-
-
- ! Set up clock to step through all intermediate file dates to be processed
- call mpas_set_time(start_time, dateTimeString=trim(config_start_time))
- call mpas_set_time(stop_time, dateTimeString=trim(config_stop_time))
- call mpas_set_timeInterval(fg_interval, S=config_fg_interval)
-
- call mpas_create_clock(fg_clock, start_time, fg_interval, stopTime=stop_time)
-
-
- ! Initialize the output file
- sfc_update_obj % time = 1
- sfc_update_obj % filename = trim(config_sfc_update_name)
-
- call mpas_output_state_init(sfc_update_obj, domain, "SFC")
-
- ! Loop over all times, interpolating the SST/SKINTEMP field from each intermediate file
- curr_time = mpas_get_clock_time(fg_clock, MPAS_NOW)
- do while (curr_time <= stop_time)
- call mpas_get_time(curr_time, dateTimeString=timeString)
- write(0,*) 'Processing ',trim(config_sfc_prefix)//':'//timeString(1:13)
-
- ! Open intermediate file
- call read_met_init(trim(config_sfc_prefix), .false., timeString(1:13), istatus)
- if (istatus /= 0) then
- write(0,*) 'Error reading ',trim(config_sfc_prefix)//':'//timeString(1:13)
- exit
- end if
-
- ! Scan through all fields in the file, looking for the SST or SKINTEMP field
- call read_next_met_field(field, istatus)
- do while (istatus == 0)
-
- !initialization of sea-surface temperature (SST) and sea-ice fraction (XICE) arrays,
- !prior to reading the input data:
-
- if (index(field % field, 'SKINTEMP') /= 0 .or. index(field % field, 'SST') /= 0) then
- fg % sst % array (1:grid%nCells) = 0.0_RKIND
-
- ! Interpolation routines use real(kind=RKIND), so copy from default real array
- allocate(slab_r8(field % nx, field % ny))
- do j=1,field % ny
- do i=1,field % nx
- slab_r8(i,j) = field % slab(i,j)
- end do
- end do
-
- !
- ! Set up map projection
- !
- call map_init(proj)
-
- if (field % iproj == PROJ_LATLON) then
- call map_set(PROJ_LATLON, proj, &
- latinc = real(field % deltalat,RKIND), &
- loninc = real(field % deltalon,RKIND), &
- knowni = 1.0_RKIND, &
- knownj = 1.0_RKIND, &
- lat1 = real(field % startlat,RKIND), &
- lon1 = real(field % startlon,RKIND))
- else if (field % iproj == PROJ_GAUSS) then
- call map_set(PROJ_GAUSS, proj, &
- nlat = nint(field % deltalat), &
- loninc = real(field % deltalon,RKIND), &
- lat1 = real(field % startlat,RKIND), &
- lon1 = real(field % startlon,RKIND))
-! nxmax = nint(360.0 / field % deltalon), &
- else if (field % iproj == PROJ_PS) then
- call map_set(PROJ_PS, proj, &
- dx = real(field % dx,RKIND), &
- truelat1 = real(field % truelat1,RKIND), &
- stdlon = real(field % xlonc,RKIND), &
- knowni = real(field % nx / 2.0,RKIND), &
- knownj = real(field % ny / 2.0,RKIND), &
- lat1 = real(field % startlat,RKIND), &
- lon1 = real(field % startlon,RKIND))
- end if
-
- ! Interpolate SST/SKINTEMP field to each MPAS grid cell
- 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 (y < 0.5) then
- y = 1.0
- else if (y >= real(field%ny)+0.5) then
- y = real(field % ny)
- end if
- if (x < 0.5) then
- lon = lon + 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- else if (x >= real(field%nx)+0.5) then
- lon = lon - 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- end if
- fg % sst % array(iCell) = interp_sequence(x, y, 1, slab_r8, 1, field % nx, 1, field % ny, 1, 1, -1.e30_RKIND, interp_list, 1)
- end do
-
- deallocate(slab_r8)
- deallocate(field % slab)
-
- else if (index(field % field, 'SEAICE') /= 0) then
- fg % xice % array (1:grid%nCells) = 0.0_RKIND
-
- ! Interpolation routines use real(kind=RKIND), so copy from default real array
- allocate(slab_r8(field % nx, field % ny))
- do j=1,field % ny
- do i=1,field % nx
- slab_r8(i,j) = field % slab(i,j)
- end do
- end do
-
- !
- ! Set up map projection
- !
- call map_init(proj)
-
- if (field % iproj == PROJ_LATLON) then
- call map_set(PROJ_LATLON, proj, &
- latinc = real(field % deltalat,RKIND), &
- loninc = real(field % deltalon,RKIND), &
- knowni = 1.0_RKIND, &
- knownj = 1.0_RKIND, &
- lat1 = real(field % startlat,RKIND), &
- lon1 = real(field % startlon,RKIND))
- else if (field % iproj == PROJ_GAUSS) then
- call map_set(PROJ_GAUSS, proj, &
- nlat = nint(field % deltalat), &
- loninc = real(field % deltalon,RKIND), &
- lat1 = real(field % startlat,RKIND), &
- lon1 = real(field % startlon,RKIND))
-! nxmax = nint(360.0 / field % deltalon), &
- else if (field % iproj == PROJ_PS) then
- call map_set(PROJ_PS, proj, &
- dx = real(field % dx,RKIND), &
- truelat1 = real(field % truelat1,RKIND), &
- stdlon = real(field % xlonc,RKIND), &
- knowni = real(field % nx / 2.0,RKIND), &
- knownj = real(field % ny / 2.0,RKIND), &
- lat1 = real(field % startlat,RKIND), &
- lon1 = real(field % startlon,RKIND))
- end if
-
- ! Interpolate SEAICE/SKINTEMP field to each MPAS grid cell
- 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 (y < 0.5) then
- y = 1.0
- else if (y >= real(field%ny)+0.5) then
- y = real(field % ny)
- end if
- if (x < 0.5) then
- lon = lon + 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- else if (x >= real(field%nx)+0.5) then
- lon = lon - 360.0
- call latlon_to_ij(proj, lat, lon, x, y)
- end if
- fg % xice % array(iCell) = interp_sequence(x, y, 1, slab_r8, 1, field % nx, 1, field % ny, 1, 1, -1.e30_RKIND, interp_list, 1)
- if (fg % xice % array(iCell) == -1.e30_RKIND) fg % xice % array(iCell) = 0.0_RKIND
-
- end do
-
- deallocate(slab_r8)
- deallocate(field % slab)
-
- else
-
- deallocate(field % slab)
- end if
-
- call read_next_met_field(field, istatus)
- end do
-
- ! Close intermediate file
- call read_met_close()
-
- ! Write the interpolated SST/SKINTEMP field as a new time slice in the MPAS output file
- call mpas_output_state_for_domain(sfc_update_obj, domain, sfc_update_obj % time)
- sfc_update_obj % time = sfc_update_obj % time + 1
-
- call mpas_advance_clock(fg_clock)
- curr_time = mpas_get_clock_time(fg_clock, MPAS_NOW)
-
- call mpas_get_time(curr_time, dateTimeString=timeString)
- state % xtime % scalar = timeString
-
- end do
-
- call mpas_output_state_finalize(sfc_update_obj, dminfo)
-
- end subroutine init_atm_test_case_sfc
-
-
integer function nearest_cell(target_lat, target_lon, &
start_cell, &
nCells, maxEdges, nEdgesOnCell, cellsOnCell, latCell, lonCell)
</font>
</pre>