<p><b>laura@ucar.edu</b> 2012-08-23 12:22:03 -0600 (Thu, 23 Aug 2012)</p><p>Minor changes<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_atmos_physics/mpas_atmphys_initialize_real.F
===================================================================
--- branches/atmos_physics/src/core_atmos_physics/mpas_atmphys_initialize_real.F 2012-08-23 18:20:35 UTC (rev 2120)
+++ branches/atmos_physics/src/core_atmos_physics/mpas_atmphys_initialize_real.F 2012-08-23 18:22:03 UTC (rev 2121)
@@ -1,4 +1,4 @@
-!=============================================================================================
+!==================================================================================================
module mpas_atmphys_initialize_real
use mpas_kind_types
use mpas_configure, only: config_met_prefix, &
@@ -7,7 +7,9 @@
config_nsoillevels, &
config_start_time, &
config_sfc_prefix
+ use mpas_dmpar
use mpas_grid_types
+ use mpas_init_atm_surface
use init_atm_hinterp
use init_atm_llxy
use init_atm_read_met
@@ -21,140 +23,17 @@
contains
-!=============================================================================================
- subroutine physics_initialize_sst(mesh,fg)
-!=============================================================================================
-
+!==================================================================================================
+ subroutine physics_initialize_real(mesh,fg,dminfo)
+!==================================================================================================
!input arguments:
- type(mesh_type),intent(in) :: mesh
+ type(mesh_type),intent(in):: mesh
+ type(dm_info),intent(in) :: dminfo
!inout arguments:
type(fg_type),intent(inout):: fg
!local variables:
- character(len=StrKIND):: timeString
- integer:: i,j,iCell,istatus
- integer,dimension(5) :: interp_list
-
- type(met_data) :: field
- type(proj_info):: proj
-
- real(kind=RKIND):: lat, lon, x, y
- real(kind=RKIND),allocatable,dimension(:,:):: slab_r8
-
-!---------------------------------------------------------------------------------------------
-
- write(0,*) '--- enter subroutine physics_initialize_sst:'
-
-!set interpolation sequence to be used for SST/SEAICE field:
- interp_list(1) = FOUR_POINT
- interp_list(2) = SEARCH
- interp_list(3) = 0
-
-!open intermediate file:
- call read_met_init(trim(config_sfc_prefix),.false.,config_start_time(1:13),istatus)
- if(istatus /= 0) &
- write(0,*) 'Error reading ',trim(config_sfc_prefix)//':'//config_start_time(1:13)
- write(0,*) 'Processing ',trim(config_sfc_prefix)//':'//config_start_time(1:13)
-
-!scan through all the fields in the file:
- call read_next_met_field(field,istatus)
- do while (istatus == 0)
-
- !initialization of the sea-surface temperature (SST) and sea-ice fraction (XICE) arrays,
- !prior to reading the input data:
- fg % sst % array (1:mesh%nCells) = 0.0_RKIND
- fg % xice % array (1:mesh%nCells) = 0.0_RKIND
-
- if(index(field % field,'SKINTEMP') /= 0 .or. &
- index(field % field,'SST' ) /= 0 .or. &
- index(field % field,'SEAICE' ) /= 0 ) then
-
- !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 field to each MPAS grid cell:
- do iCell=1,mesh % nCells
- 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)
- if (y < 0.5) then
- y = 1.0
- else if (y >= real(field%ny)+0.5) then
- y = real(field % ny)
- endif
- 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
-
- if(index(field % field,'SST') /= 0) then
- 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)
- elseif(index(field % field,'SEAICE') /= 0) then
- 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)
- endif
- end do
-
- deallocate(slab_r8)
- deallocate(field % slab)
-! exit
- end if
- call read_next_met_field(field,istatus)
- enddo
-
- end subroutine physics_initialize_sst
-
-!=============================================================================================
- subroutine physics_initialize_real(mesh,fg)
-!=============================================================================================
-!input arguments:
- type(mesh_type),intent(in) :: mesh
-
-!inout arguments:
- type(fg_type),intent(inout):: fg
-
-!local variables:
character(len=StrKIND):: initial_date
integer:: iCell,nCellsSolve
@@ -171,7 +50,7 @@
real(kind=RKIND),dimension(:),pointer:: skintemp,sst
-!---------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
write(0,*)
write(0,*) '--- enter physics_initialize_real:'
@@ -200,34 +79,15 @@
!input file. calling this subroutine will overwrite the arrays sst and seaice already read
!in the file defined by config_input_name:
if(config_input_sst) then
- call physics_initialize_sst(mesh,fg)
+ write(0,*) '--- read sea-surface temperature from auxillary file:'
+ call interp_sfc_to_MPAS(config_start_time(1:13),mesh,fg,dminfo)
+ where(landmask .eq. 0) skintemp = sst
+ endif
- if(maxval(xice(1:nCellsSolve)) == 0._RKIND .and. minval(xice(1:nCellsSolve)) == 0._RKIND) then
- write(0,*)
- write(0,*) "The input file does not contain sea-ice data. We freeze the really cold ocean instead"
- do iCell = 1, nCellsSolve
- if(landmask(iCell).eq.0 .and. sst(iCell).lt.271._RKIND) xice(iCell) = 1._RKIND
- enddo
- endif
- write(0,*) 'max sst =',maxval(fg % sst % array(1:mesh%nCells))
- write(0,*) 'min sst =',minval(fg % sst % array(1:mesh%nCells))
- write(0,*) 'max xice =',maxval(fg % xice % array(1:mesh%nCells))
- write(0,*) 'min xice =',minval(fg % xice % array(1:mesh%nCells))
+!calculates the sea-ice flag:
+ seaice(:) = 0._RKIND
+ where(xice >= 0.5_RKIND) seaice = 1._RKIND
- do iCell = 1, nCellsSolve
- !recalculate the sea-ice flag:
- if(xice(iCell) .gt. 0._RKIND) then
- seaice(iCell) = 1._RKIND
- else
- seaice(iCell) = 0._RKIND
- endif
-
- !set the skin temperature to the sea-surface temperature over the oceans:
- if(landmask(iCell).eq.0 .and. sst(iCell).gt.170._RKIND .and. sst(iCell).lt.400._RKIND) &
- skintemp(iCell) = sst(iCell)
- enddo
- endif
-
!initialization of the surface background albedo: interpolation of the monthly values to the
!initial date:
initial_date = trim(config_start_time)
@@ -287,9 +147,9 @@
end subroutine physics_initialize_real
-!=============================================================================================
+!==================================================================================================
subroutine init_soil_layers(mesh,fg)
-!=============================================================================================
+!==================================================================================================
!input arguments:
type(mesh_type),intent(in):: mesh
@@ -297,7 +157,7 @@
!inout arguments:
type(fg_type),intent(inout):: fg
-!---------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
!adjust the annual mean deep soil temperature:
call adjust_input_soiltemps(mesh,fg)
@@ -310,9 +170,9 @@
end subroutine init_soil_layers
-!=============================================================================================
+!==================================================================================================
subroutine adjust_input_soiltemps(mesh,fg)
-!=============================================================================================
+!==================================================================================================
!input arguments:
type(mesh_type),intent(in) :: mesh
@@ -329,7 +189,7 @@
real(kind=RKIND),dimension(:),pointer :: skintemp,soiltemp,tmn
real(kind=RKIND),dimension(:,:),pointer:: st_fg
-!---------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
nCellsSolve = mesh % nCellsSolve
nFGSoilLevels = mesh % nFGSoilLevels
@@ -364,9 +224,9 @@
end subroutine adjust_input_soiltemps
-!=============================================================================================
+!==================================================================================================
subroutine init_soil_layers_depth(mesh,fg)
-!=============================================================================================
+!==================================================================================================
!input arguments:
type(mesh_type),intent(in) :: mesh
@@ -377,7 +237,7 @@
!local variables:
integer:: iCell,iSoil
-!---------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
write(0,*)
write(0,*) '--- enter subroutine init_soil_layers_depth:'
@@ -425,9 +285,9 @@
end subroutine init_soil_layers_depth
-!=============================================================================================
+!==================================================================================================
subroutine init_soil_layers_properties(mesh,fg)
-!=============================================================================================
+!==================================================================================================
!input arguments:
type(mesh_type),intent(in) :: mesh
@@ -446,7 +306,7 @@
real(kind=RKIND),dimension(:,:),pointer:: dzs,zs,tslb,smois,sh2o,smcrel
real(kind=RKIND),dimension(:,:),pointer:: sm_fg,st_fg,zs_fg
-!---------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
!write(0,*)
write(0,*) '--- enter subroutine init_soil_layers_properties:'
@@ -577,9 +437,9 @@
end subroutine init_soil_layers_properties
-!=============================================================================================
+!==================================================================================================
subroutine init_seaice_points(mesh,fg)
-!=============================================================================================
+!==================================================================================================
!input arguments:
type(mesh_type),intent(in) :: mesh
@@ -606,11 +466,10 @@
real(kind=RKIND),parameter:: xice_tsk_threshold = 271.
real(kind=RKIND),parameter:: total_depth = 3. ! 3-meter soil depth.
-!---------------------------------------------------------------------------------------------
+!--------------------------------------------------------------------------------------------------
write(0,*)
write(0,*) '--- enter init_seaice_points:'
- write(0,*) '--- config_frac_seaice :', config_frac_seaice
nCellsSolve = mesh % nCellsSolve
nSoilLevels = mesh % nSoilLevels
@@ -637,17 +496,19 @@
elseif(config_frac_seaice) then
xice_threshold = 0.02
endif
+ write(0,*) '--- config_frac_seaice :', config_frac_seaice
+ write(0,*) '--- xice_threshold :', xice_threshold
!make sure that all the cells flagged as sea-ice cells are defined as ocean cells:
num_seaice_changes = 0
do iCell = 1, nCellsSolve
- if((landmask(iCell).eq.1 .and. xice(iCell).gt.0.) .or. xice(iCell).gt.200.) then
+ if((landmask(iCell).eq.1 .and. xice(iCell).gt.0._RKIND) .or. xice(iCell).gt.200._RKIND) then
num_seaice_changes = num_seaice_changes + 1
- seaice(iCell) = 0.
- xice(iCell) = 0.
+ seaice(iCell) = 0._RKIND
+ xice(iCell) = 0._RKIND
endif
enddo
- write(mess,fmt='(A,i12)') 'number of seaice cells converted to land cells=', &
+ write(mess,fmt='(A,i12)') 'number of seaice cells converted to land cells 1 =', &
num_seaice_changes
call physics_message(mess)
@@ -655,38 +516,41 @@
do iCell =1 , nCellsSolve
if(xice(iCell) .ge. xice_threshold .or. &
- (landmask(iCell).eq.0 .and. skintemp(iCell).lt.xice_tsk_threshold)) then
+ (landmask(iCell).eq.0 .and. skintemp(iCell).lt.xice_tsk_threshold)) then
num_seaice_changes = num_seaice_changes + 1
!sea-ice points are converted to land points:
- if(.not. config_frac_seaice) xice(iCell) = 1.0
- if(landmask(iCell) .eq. 0) tmn(iCell) = 271.4
+ if(.not. config_frac_seaice) xice(iCell) = 1._RKIND
+ if(landmask(iCell) .eq. 0) tmn(iCell) = 271.4_RKIND
ivgtyp(iCell) = 24 ! (isice = 24)
isltyp(iCell) = 16
- vegfra(iCell) = 0.
- landmask(iCell) = 1.
+ vegfra(iCell) = 0._RKIND
+ landmask(iCell) = 1._RKIND
do iSoil = 1, nSoilLevels
mid_point_depth = total_depth/nSoilLevels/2. &
+ (iSoil-1)*(total_depth/nSoilLevels)
tslb(iSoil,iCell) = ((total_depth-mid_point_depth) * skintemp(iCell) &
+ mid_point_depth * tmn(iCell)) / total_depth
- smois(iSoil,iCell) = 1.0
- sh2o(iSoil,iCell) = 0.0
- smcrel(iSoil,iCell) = 0.0
+ smois(iSoil,iCell) = 1._RKIND
+ sh2o(iSoil,iCell) = 0._RKIND
+ smcrel(iSoil,iCell) = 0._RKIND
enddo
elseif(xice(iCell) .lt. xice_threshold) then
- xice(iCell) = 0.
+ xice(iCell) = 0._RKIND
endif
enddo
+ write(mess,fmt='(A,i12)') 'number of seaice cells converted to land cells 2 =', &
+ num_seaice_changes
+ call physics_message(mess)
end subroutine init_seaice_points
-!=============================================================================================
+!==================================================================================================
end module mpas_atmphys_initialize_real
-!=============================================================================================
+!==================================================================================================
</font>
</pre>