diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml
index 86e37f24..6626dcc1 100644
--- a/src/core_atmosphere/Registry.xml
+++ b/src/core_atmosphere/Registry.xml
@@ -591,6 +591,8 @@
+
+
@@ -883,6 +885,8 @@
+
+
@@ -1240,6 +1244,8 @@
+
+
@@ -3898,6 +3904,12 @@
+
+
+
+
diff --git a/src/core_atmosphere/mpas_atm_core.F b/src/core_atmosphere/mpas_atm_core.F
index 7f9be119..5f6f1d48 100644
--- a/src/core_atmosphere/mpas_atm_core.F
+++ b/src/core_atmosphere/mpas_atm_core.F
@@ -572,7 +572,7 @@ subroutine atm_mpas_init_block(dminfo, stream_manager, block, mesh, dt)
call mpas_pool_get_subpool(block % structs, 'atm_input', atm_input)
call mpas_pool_get_subpool(block % structs, 'output_noahmp', output_noahmp)
call physics_tables_init(dminfo, block % configs)
- call physics_registry_init(mesh, block % configs, sfc_input)
+ call physics_registry_init(mesh, block % dimensions, block % configs, sfc_input)
call physics_run_init(block % configs, mesh, state, clock, stream_manager)
!initialization of all physics:
diff --git a/src/core_atmosphere/physics/Registry_noahmp.xml b/src/core_atmosphere/physics/Registry_noahmp.xml
index 89d980f7..ec4c7b03 100644
--- a/src/core_atmosphere/physics/Registry_noahmp.xml
+++ b/src/core_atmosphere/physics/Registry_noahmp.xml
@@ -16,7 +16,7 @@
-
diff --git a/src/core_atmosphere/physics/mpas_atmphys_control.F b/src/core_atmosphere/physics/mpas_atmphys_control.F
index bcdacad7..e898cb11 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_control.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_control.F
@@ -358,11 +358,12 @@ subroutine physics_namelist_check(configs)
end subroutine physics_namelist_check
!=================================================================================================================
- subroutine physics_registry_init(mesh,configs,sfc_input)
+ subroutine physics_registry_init(mesh,dims,configs,sfc_input)
!=================================================================================================================
!input arguments:
type(mpas_pool_type),intent(in):: mesh
+ type(mpas_pool_type),intent(in):: dims
type(mpas_pool_type),intent(in):: configs
!inout arguments:
@@ -372,40 +373,52 @@ subroutine physics_registry_init(mesh,configs,sfc_input)
logical,pointer:: config_do_restart
character(len=StrKIND),pointer:: config_lsm_scheme
integer,pointer:: nCells
+ integer, pointer :: nSoilLevels
integer,dimension(:),pointer:: landmask
- real(kind=RKIND),dimension(:,:),pointer:: dzs
+ real(kind=RKIND),dimension(:,:), pointer :: dzs
+ real(kind=RKIND),dimension(:) , pointer :: dzstop, kzs
!local variables:
integer:: iCell
+ integer:: iSoil
!-----------------------------------------------------------------------------------------------------------------
call mpas_pool_get_config(configs,'config_do_restart',config_do_restart)
call mpas_pool_get_config(configs,'config_lsm_scheme',config_lsm_scheme)
+ call mpas_pool_get_dimension(dims, 'nSoilLevels', nSoilLevels)
call mpas_pool_get_dimension(mesh,'nCells',nCells)
call mpas_pool_get_array(sfc_input,'landmask',landmask)
call mpas_pool_get_array(sfc_input,'dzs' , dzs )
+ call mpas_pool_get_array(sfc_input,'dzstop' , dzstop )
+ call mpas_pool_get_array(sfc_input,'kzs' , kzs )
-!initialization of input variables, if needed:
+
+
+!initialization of input variables, if needed:
if(.not. config_do_restart) then
lsm_select: select case(trim(config_lsm_scheme))
+ case("sf_noah","sf_noahmp")
+ ! Report soil layers as defined during the initialisation step.
+ call mpas_log_write(' ')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' Soil depth layers (physics_registry_init)')
+ call mpas_log_write('---~---')
+ iCell = min(2,nCells)
+ do iSoil=1, nSoilLevels
+ call mpas_log_write(' Layer $i - Thickness = $r' &
+ ,intArgs=(/iSoil/),realArgs=(/dzs(iSoil,iCell)/))
+ end do
+ call mpas_log_write('---~---')
+ call mpas_log_write(' ')
+
+ case default
- case("sf_noah","sf_noahmp")
- !initialize the thickness of the soil layers for the Noah scheme:
- do iCell = 1, nCells
- dzs(1,iCell) = 0.10_RKIND
- dzs(2,iCell) = 0.30_RKIND
- dzs(3,iCell) = 0.60_RKIND
- dzs(4,iCell) = 1.00_RKIND
- enddo
-
- case default
-
end select lsm_select
endif
diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F
index 7b93e7cf..e2eb8a8c 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F
@@ -842,14 +842,10 @@ subroutine lsm_noahmp_toMPAS(diag_physics,diag_physics_noahmp,output_noahmp,sfc_
snicexy(ns,i) = mpas_noahmp%snicexy(i,n)
snliqxy(ns,i) = mpas_noahmp%snliqxy(i,n)
enddo
- do ns = 1,nsnow
+ do ns = 1,nzsnow
n = ns - nsnow
zsnsoxy(ns,i) = mpas_noahmp%zsnsoxy(i,n)
enddo
- do ns = nsnow+1,nzsnow
- n = ns - nsoil + 1
- zsnsoxy(ns,i) = mpas_noahmp%zsnsoxy(i,n)
- enddo
enddo
diff --git a/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F b/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F
index d6dc1bc0..0cb2b31c 100644
--- a/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F
+++ b/src/core_atmosphere/physics/mpas_atmphys_initialize_real.F
@@ -285,9 +285,13 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs)
integer :: iCell,iSoil
integer, pointer :: nCellsSolve,nSoilLevels,nFGSoilLevels
integer, pointer :: config_nsoillevels
+ real(kind=RKIND), pointer :: config_dzstop
+ real(kind=RKIND), pointer :: config_kzs
real(kind=RKIND),dimension(:,:),pointer:: dzs_fg,zs_fg
real(kind=RKIND),dimension(:,:),pointer:: dzs,zs
+ real(kind=RKIND),dimension(:),pointer:: dzstop,kzs
+ logical :: config_fine
!-----------------------------------------------------------------------------------------------------------------
!call mpas_log_write('')
@@ -301,11 +305,32 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs)
call mpas_pool_get_array(fg, 'dzs_fg', dzs_fg)
call mpas_pool_get_array(fg, 'zs', zs)
call mpas_pool_get_array(fg, 'dzs', dzs)
+ call mpas_pool_get_array(fg, 'dzstop', dzstop)
+ call mpas_pool_get_array(fg, 'kzs', kzs)
call mpas_pool_get_config(configs, 'config_nsoillevels', config_nsoillevels)
- if(config_nsoillevels .ne. 4) &
- call physics_error_fatal('NOAH lsm uses 4 soil layers. Correct config_nsoillevels.')
+
+ call mpas_pool_get_config(configs, 'config_dzstop', config_dzstop)
+ call mpas_pool_get_config(configs, 'config_kzs' , config_kzs)
+
+ config_fine = ( config_nsoillevels >= 4 ) .and. ( config_dzstop >= 0.001 ) &
+ .and. ( config_kzs > 0.0 )
+
+ if (.not. config_fine) then
+ call mpas_log_write( ' ' )
+ call mpas_log_write( '---~---' )
+ call mpas_log_write( ' Invalid soil depth settings.' )
+ call mpas_log_write( '---~---' )
+ call mpas_log_write( ' config_nsoillevels = $i (It must be at least 4).' &
+ , intArgs=(/config_nsoillevels/))
+ call mpas_log_write( ' config_dztop = $r (It must be >= 0.001).' &
+ , realArgs=(/config_dzstop/))
+ call mpas_log_write( ' config_kzs = $r (It must be >= 0.0).' &
+ , realArgs=(/ config_kzs /))
+ call mpas_log_write(' ')
+ call physics_error_fatal(' At least one of the settings for soil depth is incorrect.')
+ end if
do iCell = 1, nCellsSolve
iSoil = 1
@@ -318,10 +343,11 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs)
enddo
do iCell = 1, nCellsSolve
- dzs(1,iCell) = 0.10_RKIND
- dzs(2,iCell) = 0.30_RKIND
- dzs(3,iCell) = 0.60_RKIND
- dzs(4,iCell) = 1.00_RKIND
+ dzstop(iCell) = config_dzstop
+ kzs(iCell) = config_kzs
+ do iSoil = 1, nSoilLevels
+ dzs(iSoil,iCell) = dzstop(iCell) * exp( kzs(iCell) * real(iSoil-1,kind=RKIND) )
+ end do
iSoil = 1
zs(iSoil,iCell) = 0.5_RKIND * dzs(iSoil,iCell)
@@ -333,6 +359,18 @@ subroutine init_soil_layers_depth(mesh, fg, dims, configs)
enddo
+ call mpas_log_write(' ')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' Soil depth layers (init_soil_layers_depth)')
+ call mpas_log_write('---~---')
+ iCell = min(2,nCellsSolve)
+ do iSoil=1, nSoilLevels
+ call mpas_log_write(' Layer $i -- Thickness = $r -- MidDepth = $r ' &
+ ,intArgs=(/iSoil/),realArgs=(/dzs(iSoil,iCell),zs(iSoil,iCell)/))
+ end do
+ call mpas_log_write('---~---')
+ call mpas_log_write(' ')
+
end subroutine init_soil_layers_depth
!=================================================================================================================
@@ -351,7 +389,7 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
type (mpas_pool_type), intent(inout) :: fg
!local variables:
- integer:: iCell,ifgSoil,iSoil
+ integer:: iCell,ifgSoil,iSoil,iCellShow
integer, pointer:: nCellsSolve,nFGSoilLevels,nSoilLevels
integer:: num_sm,num_st
integer,dimension(:),pointer:: landmask
@@ -409,8 +447,8 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
call mpas_log_write('Error in interpolation of sm_fg to MPAS grid: num_sm = $i', messageType=MPAS_LOG_CRIT, intArgs=(/num_sm/))
endif
- if(config_nsoillevels .ne. 4) &
- call physics_error_fatal('NOAH lsm uses 4 soil layers. Correct config_nsoillevels.')
+ if(config_nsoillevels < 4) &
+ call physics_error_fatal('NOAH and Noah-MP require at least 4 soil layers. Correct config_nsoillevels.')
if(.not.allocated(zhave) ) allocate(zhave(nFGSoilLevels+2,nCellsSolve) )
if(.not.allocated(st_input)) allocate(st_input(nFGSoilLevels+2,nCellsSolve))
@@ -429,13 +467,13 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
sm_input(ifgSoil+1,iCell) = sm_fg(ifgSoil,iCell)
enddo
- zhave(nFGSoilLevels+2,iCell) = 300._RKIND/100._RKIND
+ zhave(nFGSoilLevels+2,iCell) = max( zs(nSoilLevels,iCell) + dzs(nSoilLevels,iCell), 300._RKIND/100._RKIND)
st_input(nFGSoilLevels+2,iCell) = tmn(iCell)
sm_input(nFGSoilLevels+2,iCell) = sm_input(nFGSoilLevels,iCell)
if(iCell .eq. 1) then
do ifgSoil = 1,nFGSoilLevels+2
- call mpas_log_write('$i $r', intArgs=(/ifgSoil/), realArgs=(/zhave(ifgSoil,iCell)/))
+ call mpas_log_write(' ifgSoil = $i -- zhave = $r', intArgs=(/ifgSoil/), realArgs=(/zhave(ifgSoil,iCell)/))
enddo
endif
@@ -443,14 +481,15 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
!... interpolate the soil temperature, soil moisture, and soil liquid temperature to the four
! layers used in the NOAH land surface scheme:
-
+ iCellShow = 0
do iCell = 1, nCellsSolve
if(landmask(iCell) .eq. 1) then
+ if (iCellShow == 0) iCellShow = iCell
noah: do iSoil = 1 , nSoilLevels
input: do ifgSoil = 1 , nFGSoilLevels+2-1
- if(iCell .eq. 1) call mpas_log_write('$i $i $r $r $r', &
+ if(iCell .eq. 1) call mpas_log_write(' iSoil = $i -- ifgSoil = $i -- zs = $r -- zhave = $r zhaveBelow = $r', &
intArgs=(/iSoil,ifgSoil/), &
realArgs=(/zs(iSoil,iCell), zhave(ifgSoil,iCell),zhave(ifgSoil+1,iCell)/))
@@ -461,7 +500,7 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
(st_input(ifgSoil,iCell) * (zhave(ifgSoil+1,iCell)-zs(iSoil,iCell)) &
+ st_input(ifgSoil+1,iCell) * (zs(iSoil,iCell)-zhave(ifgSoil,iCell))) &
/ (zhave(ifgSoil+1,iCell)-zhave(ifgSoil,iCell))
- if(iCell .eq. 1) call mpas_log_write('$i $i $r $r $r', &
+ if(iCell .eq. 1) call mpas_log_write(' iSoil = $i -- ifgSoil = $i -- zs = $r -- zhave = $r zhaveBelow = $r', &
intArgs=(/iSoil,ifgSoil/), &
realArgs=(/zs(iSoil,iCell), zhave(ifgSoil,iCell),zhave(ifgSoil+1,iCell)/))
@@ -506,6 +545,18 @@ subroutine init_soil_layers_properties(mesh, fg, dminfo, dims, configs)
enddo
+ call mpas_log_write(' ')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' Soil initial conditions (init_soil_layers_properties)')
+ call mpas_log_write('---~---')
+ iCell = max(1,iCellShow)
+ do iSoil=1, nSoilLevels
+ call mpas_log_write(' Layer $i -- Thickness = $r -- midDepth = $r -- Temperature = $r -- Moisture = $r' &
+ ,intArgs=(/iSoil/),realArgs=(/dzs(iSoil,iCell),zs(iSoil,iCell),tslb(iSoil,iCell),smois(iSoil,iCell)/))
+ end do
+ call mpas_log_write('---~---')
+ call mpas_log_write(' ')
+
if(allocated(zhave) ) deallocate(zhave )
if(allocated(st_input)) deallocate(st_input)
if(allocated(sm_input)) deallocate(sm_input)
diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90
index 2de35ed9..e0cea648 100644
--- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90
@@ -114,10 +114,10 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO)
if ( noahmp%config%nmlist%OptSoilProperty == 1 ) then
noahmp%config%domain%SoilType(1:NumSoilLayer) = NoahmpIO%ISLTYP(I) ! soil type same in all layers
elseif ( noahmp%config%nmlist%OptSoilProperty == 2 ) then
- noahmp%config%domain%SoilType(1) = nint(NoahmpIO%SOILCL1(I)) ! soil type in layer1
- noahmp%config%domain%SoilType(2) = nint(NoahmpIO%SOILCL2(I)) ! soil type in layer2
- noahmp%config%domain%SoilType(3) = nint(NoahmpIO%SOILCL3(I)) ! soil type in layer3
- noahmp%config%domain%SoilType(4) = nint(NoahmpIO%SOILCL4(I)) ! soil type in layer4
+ noahmp%config%domain%SoilType(1) = nint(NoahmpIO%SOILCL1(I)) ! soil type in layer1
+ noahmp%config%domain%SoilType(2) = nint(NoahmpIO%SOILCL2(I)) ! soil type in layer2
+ noahmp%config%domain%SoilType(3) = nint(NoahmpIO%SOILCL3(I)) ! soil type in layer3
+ noahmp%config%domain%SoilType(4:NumSoilLayer) = nint(NoahmpIO%SOILCL4(I)) ! soil type in layers 4 and below
elseif ( noahmp%config%nmlist%OptSoilProperty == 3 ) then
noahmp%config%domain%SoilType(1:NumSoilLayer) = NoahmpIO%ISLTYP(I) ! to initialize with default
endif
diff --git a/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckGlacierMod.F90 b/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckGlacierMod.F90
index 7b5e8391..a082d983 100644
--- a/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckGlacierMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckGlacierMod.F90
@@ -5,6 +5,8 @@ module BalanceErrorCheckGlacierMod
use Machine
use NoahmpVarType
use ConstantDefineMod
+ use NoahmpFatalErrorMod,only: Noahmp_error_fatal
+ use mpas_log
implicit none
@@ -77,19 +79,22 @@ subroutine BalanceWaterCheckGlacier(noahmp)
#ifndef WRF_HYDRO
if ( abs(WaterBalanceError) > 0.1 ) then
- if ( WaterBalanceError > 0) then
- write(*,*) "The model is gaining water (WaterBalanceError is positive)"
- else
- write(*,*) "The model is losing water (WaterBalanceError is negative)"
- endif
- write(*,*) "WaterBalanceError = ",WaterBalanceError, "kg m{-2} timestep{-1}"
- write(*, &
- '(" GridIndexI GridIndexJ WaterStorageTotEnd WaterStorageTotBeg PrecipTotRefHeight &
- EvapGroundNet RunoffSurface RunoffSubsurface")')
- write(*,'(i6,1x,i6,1x,2f15.3,9f11.5)') GridIndexI, GridIndexJ, WaterStorageTotEnd, WaterStorageTotBeg, &
- PrecipTotRefHeight*MainTimeStep, EvapGroundNet*MainTimeStep, &
- RunoffSurface*MainTimeStep, RunoffSubsurface*MainTimeStep
- stop "Error: Water budget problem in NoahMP LSM"
+ call mpas_log_write(' ')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' Noah-MP water budget conservation error (glacier):')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) )
+ call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) )
+ call mpas_log_write(' WaterStorageTotBeg = $r ', realArgs = (/ WaterStorageTotBeg /) )
+ call mpas_log_write(' WaterStorageTotEnd = $r ', realArgs = (/ WaterStorageTotEnd /) )
+ call mpas_log_write(' WaterBalanceError = $r ', realArgs = (/ WaterBalanceError /) )
+ call mpas_log_write(' (positive value above means water gain).' )
+ call mpas_log_write(' PrecipTotRefHeight = $r ', realArgs = (/ PrecipTotRefHeight*MainTimeStep /) )
+ call mpas_log_write(' EvapGroundNet = $r ', realArgs = (/ EvapGroundNet*MainTimeStep /) )
+ call mpas_log_write(' RunoffSurface = $r ', realArgs = (/ RunoffSurface*MainTimeStep /) )
+ call mpas_log_write(' RunoffSubsurface = $r ', realArgs = (/ RunoffSubsurface*MainTimeStep /) )
+ call mpas_log_write(' ')
+ call Noahmp_error_fatal("Error: Water budget problem in NoahMP LSM (glacier)")
endif
#endif
@@ -133,27 +138,45 @@ subroutine BalanceEnergyCheckGlacier(noahmp)
RadSwBalanceError = RadSwDownRefHeight - (RadSwAbsSfc + RadSwReflSfc)
! print out diagnostics when error is large
if ( abs(RadSwBalanceError) > 0.01 ) then
- write(*,*) "GridIndexI, GridIndexJ = ", GridIndexI, GridIndexJ
- write(*,*) "RadSwBalanceError = ", RadSwBalanceError
- write(*,*) "RadSwDownRefHeight = ", RadSwDownRefHeight
- write(*,*) "RadSwReflSfc = ", RadSwReflSfc
- write(*,*) "RadSwAbsGrd = ", RadSwAbsGrd
- write(*,*) "RadSwAbsSfc = ", RadSwAbsSfc
- stop "Error: Solar radiation budget problem in NoahMP LSM"
+ call mpas_log_write(' ')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' Noah-MP solar radiation budget conservation error (glacier):')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) )
+ call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) )
+ call mpas_log_write(' RadSwBalanceError = $r ', realArgs = (/ RadSwBalanceError /) )
+ call mpas_log_write(' (positive value above means energy gain).' )
+ call mpas_log_write(' RadSwBalanceError = $r ', realArgs = (/ RadSwBalanceError /) )
+ call mpas_log_write(' RadSwDownRefHeight = $r ', realArgs = (/ RadSwDownRefHeight /) )
+ call mpas_log_write(' RadSwReflSfc = $r ', realArgs = (/ RadSwReflSfc /) )
+ call mpas_log_write(' RadSwAbsGrd = $r ', realArgs = (/ RadSwAbsGrd /) )
+ call mpas_log_write(' RadSwAbsSfc = $r ', realArgs = (/ RadSwAbsSfc /) )
+ call mpas_log_write('---~---')
+ call mpas_log_write(' ')
+ call Noahmp_error_fatal("Error: Solar radiation budget problem in NoahMP LSM (glacier)")
endif
! error in surface energy balance should be <0.01 W/m2
EnergyBalanceError = RadSwAbsGrd + HeatPrecipAdvSfc - (RadLwNetSfc + HeatSensibleSfc + HeatLatentGrd + HeatGroundTot)
! print out diagnostics when error is large
if ( abs(EnergyBalanceError) > 0.01 ) then
- write(*,*) 'EnergyBalanceError = ', EnergyBalanceError, ' at GridIndexI,GridIndexJ: ', GridIndexI, GridIndexJ
- write(*,'(a17,F10.4)' ) "Net longwave: ", RadLwNetSfc
- write(*,'(a17,F10.4)' ) "Total sensible: ", HeatSensibleSfc
- write(*,'(a17,F10.4)' ) "Ground evap: ", HeatLatentGrd
- write(*,'(a17,F10.4)' ) "Total ground: ", HeatGroundTot
- write(*,'(a17,4F10.4)') "Precip advected: ", HeatPrecipAdvSfc
- write(*,'(a17,F10.4)' ) "absorbed shortwave: ", RadSwAbsGrd
- stop "Error: Surface energy budget problem in NoahMP LSM"
+ call mpas_log_write(' ')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' Noah-MP energy budget conservation error (glacier):')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) )
+ call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) )
+ call mpas_log_write(' EnergyBalanceError = $r ', realArgs = (/ EnergyBalanceError /) )
+ call mpas_log_write(' (positive value above energy gain).' )
+ call mpas_log_write(' Net longwave = $r ', realArgs = (/ RadLwNetSfc /) )
+ call mpas_log_write(' Total sensible = $r ', realArgs = (/ HeatSensibleSfc /) )
+ call mpas_log_write(' Ground evap = $r ', realArgs = (/ HeatLatentGrd /) )
+ call mpas_log_write(' Total ground = $r ', realArgs = (/ HeatGroundTot /) )
+ call mpas_log_write(' Precip advected = $r ', realArgs = (/ HeatPrecipAdvSfc /) )
+ call mpas_log_write(' Absorbed shortwave = $r ', realArgs = (/ RadSwAbsGrd /) )
+ call mpas_log_write('---~---')
+ call mpas_log_write(' ')
+ call Noahmp_error_fatal("Error: Energy budget problem in NoahMP LSM (glacier)")
endif
end associate
diff --git a/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckMod.F90 b/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckMod.F90
index f076e2a5..ce58f61d 100644
--- a/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckMod.F90
+++ b/src/core_atmosphere/physics/physics_noahmp/src/BalanceErrorCheckMod.F90
@@ -5,6 +5,8 @@ module BalanceErrorCheckMod
use Machine
use NoahmpVarType
use ConstantDefineMod
+ use NoahmpFatalErrorMod,only: Noahmp_error_fatal
+ use mpas_log
implicit none
@@ -68,6 +70,7 @@ subroutine BalanceWaterCheck(noahmp)
! local variable
integer :: LoopInd ! loop index
+ real(kind=kind_noahmp) :: SoilWaterStorage
! --------------------------------------------------------------------
associate( &
@@ -116,10 +119,12 @@ subroutine BalanceWaterCheck(noahmp)
! only water balance check for every soil timestep
! Error in water balance should be < 0.1 mm
if ( SurfaceType == 1 ) then ! soil
- WaterStorageTotEnd = CanopyLiqWater + CanopyIce + SnowWaterEquiv + WaterStorageAquifer
+ SoilWaterStorage = 0.0
do LoopInd = 1, NumSoilLayer
- WaterStorageTotEnd = WaterStorageTotEnd + SoilMoisture(LoopInd) * ThicknessSnowSoilLayer(LoopInd) * 1000.0
+ SoilWaterStorage = SoilWaterStorage + SoilMoisture(LoopInd) * ThicknessSnowSoilLayer(LoopInd) * 1000.0
enddo
+ WaterStorageTotEnd = SoilWaterStorage + CanopyLiqWater + CanopyIce + SnowWaterEquiv + WaterStorageAquifer
+
! accumualted water change (only for canopy and snow during non-soil timestep)
SfcWaterTotChgAcc = SfcWaterTotChgAcc + (WaterStorageTotEnd - WaterStorageTotBeg) ! snow, canopy, and soil water change
PrecipTotAcc = PrecipTotAcc + PrecipTotRefHeight * MainTimeStep ! accumulated precip
@@ -134,21 +139,44 @@ subroutine BalanceWaterCheck(noahmp)
TileDrain)
#ifndef WRF_HYDRO
if ( abs(WaterBalanceError) > 0.1 ) then
- if ( WaterBalanceError > 0 ) then
- write(*,*) "The model is gaining water (WaterBalanceError is positive)"
- else
- write(*,*) "The model is losing water (WaterBalanceError is negative)"
- endif
- write(*,*) 'WaterBalanceError = ',WaterBalanceError, "kg m{-2} timestep{-1}"
- write(*, &
- '(" GridIndexI GridIndexJ SfcWaterTotChgAcc PrecipTotRefHeightAcc IrrigationRateMicro &
- IrrigationRateFlood EvapCanopyNetAcc EvapGroundNetAcc TranspirationAcc RunoffSurface &
- RunoffSubsurface WaterTableDepth TileDrain")')
- write(*,'(i6,i6,f10.3,10f10.5)') GridIndexI, GridIndexJ, SfcWaterTotChgAcc, PrecipTotAcc, &
- IrrigationRateMicro*1000.0, IrrigationRateFlood*1000.0, &
- EvapCanopyNetAcc, EvapGroundNetAcc, TranspirationAcc, RunoffSurface, &
- RunoffSubsurface, WaterTableDepth, TileDrain
- stop "Error: Water budget problem in NoahMP LSM"
+ call mpas_log_write(' ')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' Noah-MP water budget conservation error (land):')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) )
+ call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) )
+ call mpas_log_write(' WaterStorageTotBeg = $r ', realArgs = (/ WaterStorageTotBeg /) )
+ call mpas_log_write(' WaterStorageTotEnd = $r ', realArgs = (/ WaterStorageTotEnd /) )
+ call mpas_log_write(' WaterBalanceError = $r ', realArgs = (/ WaterBalanceError /) )
+ call mpas_log_write(' (positive value above means water gain).' )
+ call mpas_log_write('--- State:')
+ call mpas_log_write(' SoilWaterStorage = $r ', realArgs = (/ SoilWaterStorage /) )
+ call mpas_log_write(' CanopyLiqWater = $r ', realArgs = (/ CanopyLiqWater /) )
+ call mpas_log_write(' CanopyIce = $r ', realArgs = (/ CanopyIce /) )
+ call mpas_log_write(' SnowWaterEquiv = $r ', realArgs = (/ SnowWaterEquiv /) )
+ call mpas_log_write(' WaterStorageAquifer = $r ', realArgs = (/ WaterStorageAquifer /) )
+ call mpas_log_write('--- Fluxes:')
+ call mpas_log_write(' SfcWaterTotChgAcc = $r ', realArgs = (/ SfcWaterTotChgAcc /) )
+ call mpas_log_write(' PrecipTotAcc = $r ', realArgs = (/ PrecipTotAcc /) )
+ call mpas_log_write(' IrrigationRateMicro*1000.0 = $r ', realArgs = (/ IrrigationRateMicro*1000.0 /) )
+ call mpas_log_write(' IrrigationRateFlood*1000.0 = $r ', realArgs = (/ IrrigationRateFlood*1000.0 /) )
+ call mpas_log_write(' EvapCanopyNetAcc = $r ', realArgs = (/ EvapCanopyNetAcc /) )
+ call mpas_log_write(' EvapGroundNetAcc = $r ', realArgs = (/ EvapGroundNetAcc /) )
+ call mpas_log_write(' TranspirationAcc = $r ', realArgs = (/ TranspirationAcc /) )
+ call mpas_log_write(' RunoffSurface = $r ', realArgs = (/ RunoffSurface /) )
+ call mpas_log_write(' RunoffSubsurface = $r ', realArgs = (/ RunoffSubsurface /) )
+ call mpas_log_write(' WaterTableDepth = $r ', realArgs = (/ WaterTableDepth /) )
+ call mpas_log_write(' TileDrain = $r ', realArgs = (/ TileDrain /) )
+ call mpas_log_write('--- Soil state:')
+ call mpas_log_write(' NumSoilLayer = $i ',intArgs=(/ NumSoilLayer /))
+ do LoopInd=1,NumSoilLayer
+ call mpas_log_write(' Layer $i - Thickness = $r ; SoilMoisture = $r' &
+ , intArgs = (/ LoopInd /) &
+ , realArgs = (/ ThicknessSnowSoilLayer(LoopInd), SoilMoisture(LoopInd) /) )
+ end do
+ call mpas_log_write('---~---')
+ call mpas_log_write(' ')
+ call Noahmp_error_fatal("Error: Water budget problem in NoahMP LSM (land)")
endif
#endif
endif ! FlagSoilProcess
@@ -208,23 +236,31 @@ subroutine BalanceEnergyCheck(noahmp)
RadSwBalanceError = RadSwDownRefHeight - (RadSwAbsSfc + RadSwReflSfc)
! print out diagnostics when error is large
if ( abs(RadSwBalanceError) > 0.01 ) then
- write(*,*) "GridIndexI, GridIndexJ = ", GridIndexI, GridIndexJ
- write(*,*) "RadSwBalanceError = ", RadSwBalanceError
- write(*,*) "VEGETATION ---------"
- write(*,*) "RadSwDownRefHeight * VegFrac = ", RadSwDownRefHeight*VegFrac
- write(*,*) "VegFrac*RadSwAbsVeg + RadSwAbsGrd = ", VegFrac*RadSwAbsVeg+RadSwAbsGrd
- write(*,*) "VegFrac*RadSwReflVeg + RadSwReflGrd = ", VegFrac*RadSwReflVeg+RadSwReflGrd
- write(*,*) "GROUND -------"
- write(*,*) "(1 - VegFrac) * RadSwDownRefHeight = ", (1.0-VegFrac)*RadSwDownRefHeight
- write(*,*) "(1 - VegFrac) * RadSwAbsGrd = ", (1.0-VegFrac)*RadSwAbsGrd
- write(*,*) "(1 - VegFrac) * RadSwReflGrd = ", (1.0-VegFrac)*RadSwReflGrd
- write(*,*) "RadSwReflVeg = ", RadSwReflVeg
- write(*,*) "RadSwReflGrd = ", RadSwReflGrd
- write(*,*) "RadSwReflSfc = ", RadSwReflSfc
- write(*,*) "RadSwAbsVeg = ", RadSwAbsVeg
- write(*,*) "RadSwAbsGrd = ", RadSwAbsGrd
- write(*,*) "RadSwAbsSfc = ", RadSwAbsSfc
- stop "Error: Solar radiation budget problem in NoahMP LSM"
+ call mpas_log_write(' ')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' Noah-MP solar radiation budget conservation error (land):')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) )
+ call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) )
+ call mpas_log_write(' RadSwBalanceError = $r ', realArgs = (/ RadSwBalanceError /) )
+ call mpas_log_write(' (positive value above means energy gain).' )
+ call mpas_log_write('--- Vegetation ')
+ call mpas_log_write(' RadSwDownRefHeight * VegFrac = $r ', realArgs = (/ RadSwDownRefHeight*VegFrac /) )
+ call mpas_log_write(' VegFrac*RadSwAbsVeg + RadSwAbsGrd = $r ', realArgs = (/ VegFrac*RadSwAbsVeg+RadSwAbsGrd /) )
+ call mpas_log_write(' VegFrac*RadSwReflVeg + RadSwReflGrd = $r ', realArgs = (/ VegFrac*RadSwReflVeg+RadSwReflGrd /) )
+ call mpas_log_write('--- Ground ')
+ call mpas_log_write(' (1 - VegFrac) * RadSwDownRefHeight = $r ', realArgs = (/ (1.0-VegFrac)*RadSwDownRefHeight /) )
+ call mpas_log_write(' (1 - VegFrac) * RadSwAbsGrd = $r ', realArgs = (/ (1.0-VegFrac)*RadSwAbsGrd /) )
+ call mpas_log_write(' (1 - VegFrac) * RadSwReflGrd = $r ', realArgs = (/ (1.0-VegFrac)*RadSwReflGrd /) )
+ call mpas_log_write(' RadSwReflVeg = $r ', realArgs = (/ RadSwReflVeg /) )
+ call mpas_log_write(' RadSwReflGrd = $r ', realArgs = (/ RadSwReflGrd /) )
+ call mpas_log_write(' RadSwReflSfc = $r ', realArgs = (/ RadSwReflSfc /) )
+ call mpas_log_write(' RadSwAbsVeg = $r ', realArgs = (/ RadSwAbsVeg /) )
+ call mpas_log_write(' RadSwAbsGrd = $r ', realArgs = (/ RadSwAbsGrd /) )
+ call mpas_log_write(' RadSwAbsSfc = $r ', realArgs = (/ RadSwAbsSfc /) )
+ call mpas_log_write('---~---')
+ call mpas_log_write(' ')
+ call Noahmp_error_fatal("Error: Solar radiation budget problem in NoahMP LSM (land)")
endif
! error in surface energy balance should be <0.01 W/m2
@@ -233,19 +269,32 @@ subroutine BalanceEnergyCheck(noahmp)
HeatLatentTransp + HeatGroundTot + HeatLatentIrriEvap + HeatCanStorageChg)
! print out diagnostics when error is large
if ( abs(EnergyBalanceError) > 0.01 ) then
- write(*,*) 'EnergyBalanceError = ', EnergyBalanceError, ' at GridIndexI,GridIndexJ: ', GridIndexI, GridIndexJ
- write(*,'(a17,F10.4)' ) "Net solar: ", RadSwAbsSfc
- write(*,'(a17,F10.4)' ) "Net longwave: ", RadLwNetSfc
- write(*,'(a17,F10.4)' ) "Total sensible: ", HeatSensibleSfc
- write(*,'(a17,F10.4)' ) "Canopy evap: ", HeatLatentCanopy
- write(*,'(a17,F10.4)' ) "Ground evap: ", HeatLatentGrd
- write(*,'(a17,F10.4)' ) "Transpiration: ", HeatLatentTransp
- write(*,'(a17,F10.4)' ) "Total ground: ", HeatGroundTot
- write(*,'(a17,F10.4)' ) "Sprinkler: ", HeatLatentIrriEvap
- write(*,'(a17,F10.4)' ) "Canopy heat storage change: ", HeatCanStorageChg
- write(*,'(a17,4F10.4)') "Precip advected: ", HeatPrecipAdvSfc,HeatPrecipAdvCanopy,HeatPrecipAdvVegGrd,HeatPrecipAdvBareGrd
- write(*,'(a17,F10.4)' ) "Veg fraction: ", VegFrac
- stop "Error: Energy budget problem in NoahMP LSM"
+ call mpas_log_write(' ')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' Noah-MP energy budget conservation error (land):')
+ call mpas_log_write('---~---')
+ call mpas_log_write(' GridIndexI = $i ', intArgs = (/ GridIndexI /) )
+ call mpas_log_write(' GridIndexJ = $i ', intArgs = (/ GridIndexJ /) )
+ call mpas_log_write(' EnergyBalanceError = $r ', realArgs = (/ EnergyBalanceError /) )
+ call mpas_log_write(' (positive value above energy gain).' )
+ call mpas_log_write(' Net solar = $r ', realArgs = (/ RadSwAbsSfc /) )
+ call mpas_log_write(' Net longwave = $r ', realArgs = (/ RadLwNetSfc /) )
+ call mpas_log_write(' Total sensible = $r ', realArgs = (/ HeatSensibleSfc /) )
+ call mpas_log_write(' Canopy evap = $r ', realArgs = (/ HeatLatentCanopy /) )
+ call mpas_log_write(' Ground evap = $r ', realArgs = (/ HeatLatentGrd /) )
+ call mpas_log_write(' Transpiration = $r ', realArgs = (/ HeatLatentTransp /) )
+ call mpas_log_write(' Total ground = $r ', realArgs = (/ HeatGroundTot /) )
+ call mpas_log_write(' Sprinkler = $r ', realArgs = (/ HeatLatentIrriEvap /) )
+ call mpas_log_write(' Canopy heat storage change = $r ', realArgs = (/ HeatCanStorageChg /) )
+
+ call mpas_log_write(' Precip advected (surface) = $r ', realArgs = (/ HeatPrecipAdvSfc /) )
+ call mpas_log_write(' Precip advected (canopy) = $r ', realArgs = (/ HeatPrecipAdvCanopy /) )
+ call mpas_log_write(' Precip advected (veg. grnd) = $r ', realArgs = (/ HeatPrecipAdvVegGrd /) )
+ call mpas_log_write(' Precip advected (bare grnd) = $r ', realArgs = (/ HeatPrecipAdvBareGrd /) )
+ call mpas_log_write(' Veg fraction = $r ', realArgs = (/ VegFrac /) )
+ call mpas_log_write('---~---')
+ call mpas_log_write(' ')
+ call Noahmp_error_fatal("Error: Energy budget problem in NoahMP LSM (land)")
endif
end associate
diff --git a/src/core_atmosphere/physics/physics_noahmp/src/Makefile b/src/core_atmosphere/physics/physics_noahmp/src/Makefile
index 675bdf9d..6cad5478 100644
--- a/src/core_atmosphere/physics/physics_noahmp/src/Makefile
+++ b/src/core_atmosphere/physics/physics_noahmp/src/Makefile
@@ -31,6 +31,7 @@ OBJS = ConstantDefineMod.o \
IrrigationInfilPhilipMod.o \
IrrigationMicroMod.o \
MatrixSolverTriDiagonalMod.o \
+ NoahmpFatalErrorMod.o \
RunoffSubSurfaceDrainageMod.o \
RunoffSubSurfaceEquiWaterTableMod.o \
RunoffSubSurfaceGroundWaterMod.o \
@@ -168,6 +169,7 @@ IrrigationInfilPhilipMod.o: ../utility/Machine.o NoahmpVarType.o Const
IrrigationMicroMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o \
IrrigationInfilPhilipMod.o
MatrixSolverTriDiagonalMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o
+NoahmpFatalErrorMod.o: ../utility/Machine.o
RunoffSubSurfaceDrainageMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o
RunoffSubSurfaceEquiWaterTableMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o \
WaterTableEquilibriumMod.o
@@ -301,10 +303,10 @@ BiochemNatureVegMainMod.o: ../utility/Machine.o NoahmpVarType.o Co
BiochemCropMainMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o CarbonFluxCropMod.o \
CropGrowDegreeDayMod.o CropPhotosynthesisMod.o
IrrigationPrepareMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o IrrigationTriggerMod.o
-BalanceErrorCheckMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o
+BalanceErrorCheckMod.o: ../utility/Machine.o NoahmpFatalErrorMod.o NoahmpVarType.o ConstantDefineMod.o
GeneralInitMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o
GroundWaterMmfMod.o: ../utility/Machine.o NoahmpVarType.o ../drivers/hrldas/NoahmpIOVarType.o
-BalanceErrorCheckGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o
+BalanceErrorCheckGlacierMod.o: ../utility/Machine.o NoahmpFatalErrorMod.o NoahmpVarType.o ConstantDefineMod.o
EnergyMainGlacierMod.o: ../utility/Machine.o NoahmpVarType.o ConstantDefineMod.o SnowCoverGlacierMod.o \
GroundRoughnessPropertyGlacierMod.o GroundThermalPropertyGlacierMod.o \
SurfaceAlbedoGlacierMod.o SurfaceRadiationGlacierMod.o SurfaceEmissivityGlacierMod.o \
diff --git a/src/core_atmosphere/physics/physics_noahmp/src/NoahmpFatalErrorMod.F90 b/src/core_atmosphere/physics/physics_noahmp/src/NoahmpFatalErrorMod.F90
new file mode 100644
index 00000000..b42c4be3
--- /dev/null
+++ b/src/core_atmosphere/physics/physics_noahmp/src/NoahmpFatalErrorMod.F90
@@ -0,0 +1,39 @@
+module NoahmpFatalErrorMod
+ ! This module defines an error message for Noah-MP that allows for gracefully failing.
+ implicit none
+ private
+ public :: NoahMP_error_fatal
+
+ ! I would rather call mpas_derived_types from within this module, but I tried so many
+ ! things and I still cannot figure out a way of loading it (somehow mpas_log seems to
+ ! be accessible, but not mpas_derived_types). For the time being I will duplicate the
+ ! values, fully aware that this is not the correct thing to do, but please, if you
+ ! know how to make this neat, let me know...
+
+ integer, parameter :: NOAHMP_LOG_OUT = 1 !< code for message type "output"
+ integer, parameter :: NOAHMP_LOG_WARN = 2 !< code for message type "warning"
+ integer, parameter :: NOAHMP_LOG_ERR = 3 !< code for message type "error"
+ integer, parameter :: NOAHMP_LOG_CRIT = 4 !< code for message type "critical error"
+
+ contains
+
+ !---~---
+ ! Print error message then gracefully exit.
+ !---~---
+ subroutine Noahmp_error_fatal(str)
+
+ use mpas_log , only : mpas_log_write
+
+ ! input arguments:
+ character(len=*),intent(in):: str
+
+
+ call mpas_log_write(' ' , messageType=NOAHMP_LOG_ERR)
+ call mpas_log_write('---~---' , messageType=NOAHMP_LOG_ERR)
+ call mpas_log_write(' Noah-MP FATAL ERROR', messageType=NOAHMP_LOG_ERR)
+ call mpas_log_write('---~---' , messageType=NOAHMP_LOG_ERR)
+ call mpas_log_write(trim(str) , messageType=NOAHMP_LOG_ERR)
+ call mpas_log_write('Noah-MP abort' , messageType=NOAHMP_LOG_CRIT)
+ end subroutine Noahmp_error_fatal
+ !---~---
+end module NoahmpFatalErrorMod
diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml
index 7b8905ff..efb732fd 100644
--- a/src/core_init_atmosphere/Registry.xml
+++ b/src/core_init_atmosphere/Registry.xml
@@ -262,6 +262,19 @@
+
+
+
+
+
+
+
+
+
@@ -1218,6 +1233,14 @@
+
+
+
+