<p><b>laura@ucar.edu</b> 2012-08-28 15:17:36 -0600 (Tue, 28 Aug 2012)</p><p>moved the call to subroutine init_seaice_points and renamed subroutine to physics_init_seaice. Changed argument list for subroutines physics_init_sst and physics_init_seaice so that they can be accessed from the model side too.<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-28 21:14:48 UTC (rev 2131)
+++ branches/atmos_physics/src/core_atmos_physics/mpas_atmphys_initialize_real.F 2012-08-28 21:17:36 UTC (rev 2132)
@@ -57,12 +57,12 @@
nCellsSolve = mesh % nCellsSolve
- landmask => mesh % landmask % array
- albedo12m => mesh % albedo12m % array
- greenfrac => mesh % greenfrac % array
- shdmin => mesh % shdmin % array
- shdmax => mesh % shdmax % array
- snoalb => mesh % snoalb % array
+ landmask => mesh % landmask % array
+ albedo12m => mesh % albedo12m % array
+ greenfrac => mesh % greenfrac % array
+ shdmin => mesh % shdmin % array
+ shdmax => mesh % shdmax % array
+ snoalb => mesh % snoalb % array
sfc_albbck => fg % sfc_albbck % array
vegfra => fg % vegfra % array
@@ -81,12 +81,11 @@
if(config_input_sst) then
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
+ call physics_init_sst(mesh,fg)
-!calculates the sea-ice flag:
- seaice(:) = 0._RKIND
- where(xice >= 0.5_RKIND) seaice = 1._RKIND
+!initialize seaice points:
+ call physics_init_seaice(mesh,fg)
!initialization of the surface background albedo: interpolation of the monthly values to the
!initial date:
@@ -129,11 +128,8 @@
enddo
!initialization of soil layers properties:
- call init_soil_layers(mesh,fg)
+ call init_soil_layers(mesh,fg,dminfo)
-!adjustment of all surface fields for seaice points:
- call init_seaice_points(mesh,fg)
-
!define xland over land and ocean:
do iCell = 1, nCellsSolve
if(landmask(iCell) .eq. 1 .or. (landmask(iCell).eq.0 .and. seaice(iCell).eq.1._RKIND)) then
@@ -148,11 +144,12 @@
end subroutine physics_initialize_real
!==================================================================================================
- subroutine init_soil_layers(mesh,fg)
+ subroutine init_soil_layers(mesh,fg,dminfo)
!==================================================================================================
!input arguments:
type(mesh_type),intent(in):: mesh
+ type(dm_info),intent(in) :: dminfo
!inout arguments:
type(fg_type),intent(inout):: fg
@@ -166,7 +163,7 @@
call init_soil_layers_depth(mesh,fg)
!initialize the temperature, moisture, and liquid water of the individual soil layers:
- call init_soil_layers_properties(mesh,fg)
+ call init_soil_layers_properties(mesh,fg,dminfo)
end subroutine init_soil_layers
@@ -207,12 +204,12 @@
if(landmask(iCell) .eq. 1) then
!adjust the annual deep mean soil temperature and skin temperatures over land:
- tmn(iCell) = soiltemp(iCell) - 0.0065 * ter(iCell)
- skintemp(iCell) = skintemp(iCell) - 0.0065 * (ter(iCell)-soilz(iCell))
+ tmn(iCell) = soiltemp(iCell) - 0.0065_RKIND * ter(iCell)
+ skintemp(iCell) = skintemp(iCell) - 0.0065_RKIND * (ter(iCell)-soilz(iCell))
!adjust the soil layer temperatures:
do ifgSoil = 1, nFGSoilLevels
- st_fg(ifgSoil,iCell) = st_fg(ifgSoil,iCell) - 0.0065 * (ter(iCell)-soilz(iCell))
+ st_fg(ifgSoil,iCell) = st_fg(ifgSoil,iCell) - 0.0065_RKIND * (ter(iCell)-soilz(iCell))
enddo
elseif(landmask(iCell) .eq. 0) then
@@ -229,7 +226,7 @@
!==================================================================================================
!input arguments:
- type(mesh_type),intent(in) :: mesh
+ type(mesh_type),intent(in):: mesh
!inout arguments:
type(fg_type),intent(inout):: fg
@@ -248,37 +245,28 @@
do iCell = 1, mesh % nCells
iSoil = 1
- fg % zs_fg % array(iSoil,iCell) = 0.5 * fg % dzs_fg % array(iSoil,iCell)
-! if(iCell .eq. 1) write(0,101) iSoil,fg % dzs_fg % array(iSoil,iCell), &
-! fg % zs_fg % array(iSoil,iCell)
+ fg % zs_fg % array(iSoil,iCell) = 0.5_RKIND * fg % dzs_fg % array(iSoil,iCell)
do iSoil = 2, mesh % nFGSoilLevels
fg % zs_fg % array(iSoil,iCell) = fg % zs_fg % array(iSoil-1,iCell) &
- + 0.5 * fg % dzs_fg % array(iSoil-1,iCell) &
- + 0.5 * fg % dzs_fg % array(iSoil,iCell)
-! if(iCell .eq. 1) write(0,101) iSoil,fg % dzs_fg % array(iSoil,iCell), &
-! fg % zs_fg % array(iSoil,iCell)
+ + 0.5_RKIND * fg % dzs_fg % array(iSoil-1,iCell) &
+ + 0.5_RKIND * fg % dzs_fg % array(iSoil,iCell)
enddo
enddo
101 format(i4,2(1x,e15.8))
do iCell = 1, mesh % nCells
- fg % dzs % array(1,iCell) = 0.10
- fg % dzs % array(2,iCell) = 0.30
- fg % dzs % array(3,iCell) = 0.60
- fg % dzs % array(4,iCell) = 1.00
+ fg % dzs % array(1,iCell) = 0.10_RKIND
+ fg % dzs % array(2,iCell) = 0.30_RKIND
+ fg % dzs % array(3,iCell) = 0.60_RKIND
+ fg % dzs % array(4,iCell) = 1.00_RKIND
iSoil = 1
- fg % zs % array(iSoil,iCell) = 0.5 * fg % dzs % array(iSoil,iCell)
-! if(iCell .eq. 1) write(0,101) iSoil,fg % dzs % array(iSoil,iCell), &
-! fg % zs % array(iSoil,iCell)
-
+ fg % zs % array(iSoil,iCell) = 0.5_RKIND * fg % dzs % array(iSoil,iCell)
do iSoil = 2, mesh % nSoilLevels
- fg % zs % array(iSoil,iCell) = fg % zs % array(iSoil-1,iCell) &
- + 0.5 * fg % dzs % array(iSoil-1,iCell) &
- + 0.5 * fg % dzs % array(iSoil,iCell)
-! if(iCell .eq. 1) write(0,101) iSoil,fg % dzs % array(iSoil,iCell), &
-! fg % zs % array(iSoil,iCell)
+ fg % zs % array(iSoil,iCell) = fg % zs % array(iSoil-1,iCell) &
+ + 0.5_RKIND * fg % dzs % array(iSoil-1,iCell) &
+ + 0.5_RKIND * fg % dzs % array(iSoil,iCell)
enddo
enddo
@@ -286,11 +274,12 @@
end subroutine init_soil_layers_depth
!==================================================================================================
- subroutine init_soil_layers_properties(mesh,fg)
+ subroutine init_soil_layers_properties(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
@@ -298,6 +287,7 @@
!local variables:
integer:: iCell,ifgSoil,iSoil,is
integer:: nCells,nFGSoilLevels,nSoilLevels
+ integer:: num_sm,num_st
integer,dimension(:),pointer:: landmask
real(kind=RKIND),dimension(:,:),allocatable:: zhave,sm_input,st_input
@@ -332,6 +322,24 @@
skintemp => fg % skintemp % array
tmn => fg % tmn % array
+!check that interpolation of the meteorological data to the MPAS grid did not create negative
+!values for the first-guess soil temperatures and soil moistures.
+ num_sm = 0
+ num_st = 0
+ do iCell = 1, nCells
+ do ifgSoil = 1, nFGSoilLevels
+ if(st_fg(ifgSoil,iCell) .le. 0._RKIND) num_st = num_st + 1
+ if(sm_fg(ifgSoil,iCell) .lt. 0._RKIND) num_sm = num_sm + 1
+ enddo
+ enddo
+ if(num_st .gt. 0) then
+ write(0,*) 'Error in interpolation of st_fg to MPAS grid: num_st =', num_st
+ call mpas_dmpar_abort(dminfo)
+ elseif(num_sm .gt. 0) then
+ write(0,*) 'Error in interpolation of sm_fg to MPAS grid: num_sm =', num_sm
+ call mpas_dmpar_abort(dminfo)
+ endif
+
if(config_nsoillevels .ne. 4) &
call physics_error_fatal('NOAH lsm uses 4 soil layers. Correct config_nsoillevels.')
@@ -342,17 +350,17 @@
do iCell = 1, nCells
ifgSoil = 1
- zhave(ifgSoil,iCell) = 0.
+ zhave(ifgSoil,iCell) = 0._RKIND
st_input(ifgSoil,iCell) = skintemp(iCell)
sm_input(ifgSoil,iCell) = sm_fg(ifgSoil+1,iCell)
do ifgSoil = 1, nFGSoilLevels
- zhave(ifgSoil+1,iCell) = zs_fg(ifgSoil,iCell) / 100.
+ zhave(ifgSoil+1,iCell) = zs_fg(ifgSoil,iCell) / 100._RKIND
st_input(ifgSoil+1,iCell) = st_fg(ifgSoil,iCell)
sm_input(ifgSoil+1,iCell) = sm_fg(ifgSoil,iCell)
enddo
- zhave(nFGSoilLevels+2,iCell) = 300./100.
+ zhave(nFGSoilLevels+2,iCell) = 300._RKIND/100._RKIND
st_input(nFGSoilLevels+2,iCell) = tmn(iCell)
sm_input(nFGSoilLevels+2,iCell) = sm_input(nFGSoilLevels,iCell)
@@ -391,8 +399,8 @@
+ sm_input(ifgSoil+1,iCell) * (zs(iSoil,iCell)-zhave(ifgSoil,iCell))) &
/ (zhave(ifgSoil+1,iCell)-zhave(ifgSoil,iCell))
- sh2o(iSoil,iCell) = 0.0
- smcrel(iSoil,iCell) = 0.0
+ sh2o(iSoil,iCell) = 0._RKIND
+ smcrel(iSoil,iCell) = 0._RKIND
exit input
endif
@@ -405,9 +413,9 @@
!fill the soil temperatures with the skin temperatures over oceans:
do iSoil = 1, nSoilLevels
tslb(iSoil,iCell) = skintemp(iCell)
- smois(iSoil,iCell) = 1.0
- sh2o(iSoil,iCell) = 1.0
- smcrel(iSoil,iCell) = 0.0
+ smois(iSoil,iCell) = 1._RKIND
+ sh2o(iSoil,iCell) = 1._RKIND
+ smcrel(iSoil,iCell) = 0._RKIND
enddo
endif
@@ -418,10 +426,10 @@
do iCell = 1, nCells
- if(landmask(iCell).eq. 1 .and. tslb(1,iCell).gt.170. .and. tslb(1,iCell).lt.400. .and. &
- smois(1,iCell).lt.0.005) then
+ if(landmask(iCell).eq. 1 .and. tslb(1,iCell).gt.170._RKIND .and. tslb(1,iCell).lt.400._RKIND &
+ .and. smois(1,iCell).lt.0.005_RKIND) then
do iSoil = 1, nSoilLevels
- smois(iSoil,iCell) = 0.005
+ smois(iSoil,iCell) = 0.005_RKIND
enddo
endif
@@ -438,18 +446,68 @@
end subroutine init_soil_layers_properties
!==================================================================================================
- subroutine init_seaice_points(mesh,fg)
+ subroutine physics_init_sst(mesh,input)
!==================================================================================================
+!input arguments:
+ type(mesh_type),intent(in):: mesh
+
+#if !defined(non_hydrostatic_core)
+!inout arguments: this subroutine is called from the MPAS initialization side.
+ type(fg_type),intent(inout):: input
+#else
+!inout arguments: this subroutine is called from the MPAS model side.
+ type(sfc_input_type),intent(inout):: input
+#endif
+
+!local variables:
+ integer:: iCell,nCells
+ integer,dimension(:),pointer:: landmask
+
+ real(kind=RKIND),dimension(:),pointer :: sst,tsk
+ real(kind=RKIND),dimension(:,:),pointer:: tslb
+
+!--------------------------------------------------------------------------------------------------
+ write(0,*)
+ write(0,*) '--- enter subroutine physics_update_sst:'
+
+!initialization:
+ nCells = mesh % nCells
+
+ landmask => mesh % landmask % array
+ sst => input % sst % array
+ tsk => input % skintemp % array
+ tslb => input % tslb % array
+
+!update the skin temperature and the soil temperature of the first soil layer with the updated
+!sea-surface temperatures:
+ do iCell = 1, nCells
+ if(landmask(iCell) == 0) then
+ tsk(iCell) = sst(iCell)
+ endif
+ enddo
+ write(0,*) '--- end subroutine physics_update_sst:'
+
+ end subroutine physics_init_sst
+
+!==================================================================================================
+ subroutine physics_init_seaice(mesh,input)
+!==================================================================================================
+
!input arguments:
type(mesh_type),intent(in) :: mesh
-!inout arguments:
- type(fg_type),intent(inout):: fg
+#if !defined(non_hydrostatic_core)
+!inout arguments: this subroutine is called from the MPAS initialization side.
+ type(fg_type),intent(inout):: input
+#else
+!inout arguments: this subroutine is called from the MPAS model side.
+ type(sfc_input_type),intent(inout):: input
+#endif
!local variables:
character(len=StrKIND):: mess
- integer:: iCell,iSoil,nCellsSolve,nSoilLevels
+ integer:: iCell,iSoil,nCells,nSoilLevels
integer:: num_seaice_changes
integer,dimension(:),pointer:: landmask,isltyp,ivgtyp
@@ -457,7 +515,7 @@
real(kind=RKIND):: mid_point_depth
real(kind=RKIND),dimension(:),pointer :: vegfra
real(kind=RKIND),dimension(:),pointer :: seaice,xice
- real(kind=RKIND),dimension(:),pointer :: skintemp,sst,tmn
+ real(kind=RKIND),dimension(:),pointer :: skintemp,tmn
real(kind=RKIND),dimension(:,:),pointer:: tslb,smois,sh2o,smcrel
!note that this threshold is also defined in module_physics_vars.F.It is defined here to avoid
@@ -467,41 +525,35 @@
real(kind=RKIND),parameter:: total_depth = 3. ! 3-meter soil depth.
!--------------------------------------------------------------------------------------------------
-
write(0,*)
- write(0,*) '--- enter init_seaice_points:'
+ write(0,*) '--- enter physics_init_seaice:'
- nCellsSolve = mesh % nCellsSolve
+ nCells = mesh % nCells
nSoilLevels = mesh % nSoilLevels
landmask => mesh % landmask % array
isltyp => mesh % soilcat_top % array
ivgtyp => mesh % lu_index % array
- seaice => fg % seaice % array
- xice => fg % xice % array
- vegfra => fg % vegfra % array
+ seaice => input % seaice % array
+ xice => input % xice % array
+ vegfra => input % vegfra % array
- skintemp => fg % skintemp % array
- sst => fg % sst % array
- tmn => fg % tmn % array
+ skintemp => input % skintemp % array
+ tmn => input % tmn % array
- tslb => fg % tslb % array
- smois => fg % smois % array
- sh2o => fg % sh2o % array
- smcrel => fg % smcrel % array
+ tslb => input % tslb % array
+ smois => input % smois % array
+ sh2o => input % sh2o % array
+ smcrel => input % smcrel % array
- if(.not. config_frac_seaice) then
- xice_threshold = 0.5
- elseif(config_frac_seaice) then
- xice_threshold = 0.02
- endif
- write(0,*) '--- config_frac_seaice :', config_frac_seaice
- write(0,*) '--- xice_threshold :', xice_threshold
+ do iCell = 1, nCells
+ seaice(iCell) = 0._RKIND
+ enddo
!make sure that all the cells flagged as sea-ice cells are defined as ocean cells:
num_seaice_changes = 0
- do iCell = 1, nCellsSolve
+ do iCell = 1, nCells
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._RKIND
@@ -512,8 +564,25 @@
num_seaice_changes
call physics_message(mess)
+!assign the threshold value for xice as a function of config_frac_seaice:
+ if(.not. config_frac_seaice) then
+ xice_threshold = 0.5_RKIND
+ do iCell = 1,nCells
+ if(xice(iCell) >= xice_threshold) then
+ xice(iCell) = 1._RKIND
+ else
+ xice(iCell) = 0._RKIND
+ endif
+ enddo
+ elseif(config_frac_seaice) then
+ xice_threshold = 0.02
+ endif
+ write(0,*) '--- config_frac_seaice :', config_frac_seaice
+ write(0,*) '--- xice_threshold :', xice_threshold
+
+!convert seaice points to land points:
num_seaice_changes = 0
- do iCell =1 , nCellsSolve
+ do iCell = 1, nCells
if(xice(iCell) .ge. xice_threshold .or. &
(landmask(iCell).eq.0 .and. skintemp(iCell).lt.xice_tsk_threshold)) then
@@ -548,8 +617,14 @@
num_seaice_changes
call physics_message(mess)
- end subroutine init_seaice_points
+!finally, update the sea-ice flag:
+ do iCell = 1, nCells
+ if(xice(iCell) > 0._RKIND) seaice(iCell) = 1._RKIND
+ enddo
+ write(0,*) '--- end physics_init_seaice:'
+ end subroutine physics_init_seaice
+
!==================================================================================================
end module mpas_atmphys_initialize_real
!==================================================================================================
</font>
</pre>