From 8fcda2b8ec12604f684428fa4f66fcd07b7bcf00 Mon Sep 17 00:00:00 2001 From: Marcos Longo Date: Tue, 26 May 2026 07:55:12 -0300 Subject: [PATCH 1/2] This fixes a cryptic bug that only occurs when the code is modified to run with more than four soil layers. The logic in mpas_atmphys_driver_lsm_noahmp.F can be simplified to a single loop that will always work. Likewise, Registry_noahmp.xml does not need to have the combined number of soil plus snow levels (nzSoilLevels) hardcoded, it can be written as the sum of the levels so it always works. --- src/core_atmosphere/physics/Registry_noahmp.xml | 2 +- .../physics/mpas_atmphys_driver_lsm_noahmp.F | 6 +----- 2 files changed, 2 insertions(+), 6 deletions(-) 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_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 From b592e0cd5f4b8dc93bc36291e21453fbc71fd657 Mon Sep 17 00:00:00 2001 From: Marcos Longo Date: Tue, 26 May 2026 09:19:58 -0300 Subject: [PATCH 2/2] First commit implementing a flexible soil grid. --- src/core_atmosphere/Registry.xml | 12 ++ src/core_atmosphere/mpas_atm_core.F | 2 +- .../physics/mpas_atmphys_control.F | 41 +++-- .../physics/mpas_atmphys_initialize_real.F | 79 ++++++++-- .../drivers/mpas/ConfigVarInTransferMod.F90 | 8 +- .../src/BalanceErrorCheckGlacierMod.F90 | 79 ++++++---- .../src/BalanceErrorCheckMod.F90 | 143 ++++++++++++------ .../physics/physics_noahmp/src/Makefile | 6 +- .../src/NoahmpFatalErrorMod.F90 | 39 +++++ src/core_init_atmosphere/Registry.xml | 23 +++ 10 files changed, 322 insertions(+), 110 deletions(-) create mode 100644 src/core_atmosphere/physics/physics_noahmp/src/NoahmpFatalErrorMod.F90 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/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_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 @@ + + + +