diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 1327acb..c7acc0c 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -503,6 +503,7 @@ #ifdef DO_PHYSICS + @@ -4295,7 +4296,10 @@ + description="dominant soil texture category"/> + + diff --git a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F index 0116dcf..ba5faf2 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm.F @@ -139,6 +139,7 @@ subroutine allocate_lsm if(.not.allocated(cpm_p) ) allocate(cpm_p(ims:ime,jms:jme) ) if(.not.allocated(cqs2_p) ) allocate(cqs2_p(ims:ime,jms:jme) ) if(.not.allocated(isltyp_p) ) allocate(isltyp_p(ims:ime,jms:jme) ) + if(.not.allocated(isctyp_p) ) allocate(isctyp_p(ims:ime,jms:jme) ) if(.not.allocated(ivgtyp_p) ) allocate(ivgtyp_p(ims:ime,jms:jme) ) if(.not.allocated(glw_p) ) allocate(glw_p(ims:ime,jms:jme) ) if(.not.allocated(grdflx_p) ) allocate(grdflx_p(ims:ime,jms:jme) ) @@ -216,6 +217,7 @@ subroutine deallocate_lsm if(allocated(gsw_p) ) deallocate(gsw_p ) if(allocated(hfx_p) ) deallocate(hfx_p ) if(allocated(isltyp_p) ) deallocate(isltyp_p ) + if(allocated(isctyp_p) ) deallocate(isctyp_p ) if(allocated(ivgtyp_p) ) deallocate(ivgtyp_p ) if(allocated(lai_p) ) deallocate(lai_p ) if(allocated(lh_p) ) deallocate(lh_p ) @@ -279,7 +281,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) character(len=StrKIND),pointer:: config_microp_scheme, & config_convection_scheme - integer,dimension(:),pointer:: isltyp,ivgtyp + integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp real(kind=RKIND),dimension(:),pointer :: acsnom,acsnow,canwat,chs,chs2,chklowq,cpm,cqs2,glw, & grdflx,gsw,hfx,lai,lh,noahres,potevp,qfx,qgh,qsfc, & @@ -336,6 +338,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) call mpas_pool_get_array(diag_physics,'znt' ,znt ) call mpas_pool_get_array(sfc_input,'isltyp' ,isltyp ) + call mpas_pool_get_array(sfc_input,'isctyp' ,isctyp ) call mpas_pool_get_array(sfc_input,'ivgtyp' ,ivgtyp ) call mpas_pool_get_array(sfc_input,'shdmin' ,shdmin ) call mpas_pool_get_array(sfc_input,'shdmax' ,shdmax ) @@ -409,6 +412,7 @@ subroutine lsm_from_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) znt_p(i,j) = znt(i) isltyp_p(i,j) = isltyp(i) + isctyp_p(i,j) = isctyp(i) ivgtyp_p(i,j) = ivgtyp(i) shdmin_p(i,j) = shdmin(i) shdmax_p(i,j) = shdmax(i) @@ -490,7 +494,7 @@ subroutine lsm_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) !local pointers: character(len=StrKIND),pointer:: config_microp_scheme - integer,dimension(:),pointer:: isltyp,ivgtyp + integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp real(kind=RKIND),dimension(:),pointer :: acsnom,acsnow,canwat,chs,chs2,chklowq,cpm,cqs2,glw, & grdflx,gsw,hfx,lai,lh,noahres,potevp,qfx,qgh,qsfc, & @@ -544,6 +548,7 @@ subroutine lsm_to_MPAS(configs,mesh,diag_physics,sfc_input,its,ite) call mpas_pool_get_array(diag_physics,'znt' ,znt ) call mpas_pool_get_array(sfc_input,'isltyp' ,isltyp ) + call mpas_pool_get_array(sfc_input,'isctyp' ,isctyp ) call mpas_pool_get_array(sfc_input,'ivgtyp' ,ivgtyp ) call mpas_pool_get_array(sfc_input,'shdmin' ,shdmin ) call mpas_pool_get_array(sfc_input,'shdmax' ,shdmax ) 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 8bbec89..72147f2 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F +++ b/src/core_atmosphere/physics/mpas_atmphys_driver_lsm_noahmp.F @@ -58,7 +58,7 @@ subroutine lsm_noahmp_fromMPAS(configs,mesh,diag,diag_physics,diag_physics_noahm integer:: i,its,ite integer:: n,ns,nsoil,nsnow,nzsnow - integer,dimension(:),pointer:: isltyp,ivgtyp + integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp real(kind=RKIND),dimension(:),pointer:: latCell,lonCell real(kind=RKIND),dimension(:),pointer:: shdmax,shdmin,vegfra,tmn,xice,xland diff --git a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F index 07481fd..bce73fe 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F +++ b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahinit.F @@ -91,7 +91,7 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input) character(len=StrKIND),pointer:: mminlu,mminsl integer,pointer:: nCells,nSoilLevels - integer,dimension(:),pointer:: ivgtyp,isltyp + integer,dimension(:),pointer:: ivgtyp,isltyp,isctyp real(kind=RKIND),dimension(:),pointer:: snoalb,snow,snowh real(kind=RKIND),dimension(:,:),pointer:: tslb,smois,sh2o @@ -118,6 +118,7 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input) call mpas_pool_get_dimension(mesh,'nSoilLevels',nSoilLevels) call mpas_pool_get_array(sfc_input,'isltyp', isltyp) + call mpas_pool_get_array(sfc_input,'isctyp', isctyp) call mpas_pool_get_array(sfc_input,'ivgtyp', ivgtyp) call mpas_pool_get_array(sfc_input,'sh2o' , sh2o ) call mpas_pool_get_array(sfc_input,'smois' , smois ) @@ -149,6 +150,20 @@ subroutine lsminit(dminfo,mesh,configs,diag_physics,sfc_input) call physics_error_fatal("module_sf_noahlsm.F: lsminit: out of range value "// & "of ISLTYP. Is this field in the input?" ) + ! Make sure all cells have reasonable soil colour indices. + errflag = 0 + do iCell = 1, nCells + if(isctyp(iCell) < 1) then + errflag = 1 + write(err_message,*) "module_sf_noahlsm.F: lsminit: out of range ISCTYP ", & + iCell,isctyp(iCell),isltyp(iCell),ivgtyp(iCell) + call physics_message(err_message) + endif + end do + if (errflag == 1) & + call physics_error_fatal("module_sf_noahlsm.F: lsminit: out of range value "// & + "of ISCTYP. Is this field in the input?" ) + !initializes soil liquid water content SH2O: do iCell = 1, nCells bx = bb(isltyp(iCell)) diff --git a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F index e720166..13a28e2 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F +++ b/src/core_atmosphere/physics/mpas_atmphys_lsm_noahmpinit.F @@ -231,7 +231,7 @@ subroutine noahmp_init(configs,mesh,diag_physics,diag_physics_noahmp,output_noah logical,pointer:: urban_physics integer,pointer:: nsoilcomps - integer,dimension(:),pointer:: isltyp,ivgtyp + integer,dimension(:),pointer:: isltyp,isctyp,ivgtyp integer,dimension(:),pointer:: isnowxy integer,dimension(:),pointer:: irnumsi,irnummi,irnumfi @@ -319,11 +319,14 @@ subroutine noahmp_init(configs,mesh,diag_physics,diag_physics_noahmp,output_noah !--- initialization of time-invariant surface variables needed in subroutine NoahmpInitMain: call mpas_pool_get_array(sfc_input,'dzs' ,dzs ) call mpas_pool_get_array(sfc_input,'isltyp',isltyp) + call mpas_pool_get_array(sfc_input,'isctyp',isctyp) call mpas_pool_get_array(sfc_input,'ivgtyp',ivgtyp) do i = its, ite mpas_noahmp%isltyp(i) = isltyp(i) + mpas_noahmp%isctyp(i) = isctyp(i) mpas_noahmp%ivgtyp(i) = ivgtyp(i) + enddo do ns = 1, nsoil mpas_noahmp%dzs(ns) = dzs(ns,its) diff --git a/src/core_atmosphere/physics/mpas_atmphys_update_surface.F b/src/core_atmosphere/physics/mpas_atmphys_update_surface.F index 6e2057d..4d6a59e 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_update_surface.F +++ b/src/core_atmosphere/physics/mpas_atmphys_update_surface.F @@ -151,7 +151,7 @@ subroutine physics_update_sst(dminfo,config_frac_seaice,mesh,sfc_input,diag_phys call mpas_pool_get_array(sfc_input,'isice' ,isice ) call mpas_pool_get_array(sfc_input,'iswater' ,iswater ) - call mpas_pool_get_array(sfc_input,'isltyp' ,isltyp ) + call mpas_pool_get_array(sfc_input,'isltyp' ,isltyp ) call mpas_pool_get_array(sfc_input,'ivgtyp' ,ivgtyp ) call mpas_pool_get_array(sfc_input,'landmask' ,landmask ) call mpas_pool_get_array(sfc_input,'vegfra' ,vegfra ) diff --git a/src/core_atmosphere/physics/mpas_atmphys_vars.F b/src/core_atmosphere/physics/mpas_atmphys_vars.F index 28a2364..f69f5b1 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_vars.F +++ b/src/core_atmosphere/physics/mpas_atmphys_vars.F @@ -822,6 +822,7 @@ module mpas_atmphys_vars integer,dimension(:,:),allocatable:: & isltyp_p, &!dominant soil type category [-] + isctyp_p, &!dominant soil colour category [-] ivgtyp_p !dominant vegetation category [-] real(kind=RKIND),dimension(:),allocatable:: & diff --git a/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md b/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md index ae45c39..6427e6b 100644 --- a/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md +++ b/src/core_atmosphere/physics/physics_noahmp/RELEASE_NOTES.md @@ -400,7 +400,7 @@ - but in the file listed in the namelist as: HRLDAS_SETUP_FILE = " - The initialization fields are: SNOW,CANWAT,TSK,TSLB,SMOIS - - This file also contains the static grid/domain information: XLAT,XLONG,TMN,HGT,SEAICE,MAPFAC_MX,MAPFAC_MY,SHDMAX,SHDMIN,XLAND,IVGTYP,ISLTYP,DZS,ZS + - This file also contains the static grid/domain information: XLAT,XLONG,TMN,HGT,SEAICE,MAPFAC_MX,MAPFAC_MY,SHDMAX,SHDMIN,XLAND,IVGTYP,ISLTYP,ISCTYP,DZS,ZS - This file can also contains some optional fields: LAI - NOTE: a WRF input file can be used as a HRLDAS_SETUP_FILE 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 2de35ed..1cba671 100644 --- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/ConfigVarInTransferMod.F90 @@ -62,7 +62,7 @@ subroutine ConfigVarInTransfer(noahmp, NoahmpIO) ! config domain variable noahmp%config%domain%SurfaceType = 1 noahmp%config%domain%NumSwRadBand = 2 - noahmp%config%domain%SoilColor = 4 + noahmp%config%domain%SoilColor = NoahmpIO%ISCTYP(I) noahmp%config%domain%NumCropGrowStage = 8 noahmp%config%domain%FlagSoilProcess = NoahmpIO%calculate_soil noahmp%config%domain%NumSoilTimeStep = NoahmpIO%soil_update_steps diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90 index 12df9b1..9aba978 100644 --- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarFinalizeMod.F90 @@ -37,6 +37,7 @@ subroutine NoahmpIOVarFinalizeDefault(NoahmpIO) if ( allocated (NoahmpIO%zsoil) ) deallocate ( NoahmpIO%zsoil ) ! depth to soil interfaces [m] if ( allocated (NoahmpIO%ivgtyp) ) deallocate ( NoahmpIO%ivgtyp ) ! vegetation type if ( allocated (NoahmpIO%isltyp) ) deallocate ( NoahmpIO%isltyp ) ! soil type + if ( allocated (NoahmpIO%isctyp) ) deallocate ( NoahmpIO%isctyp ) ! soil colour class if ( allocated (NoahmpIO%vegfra) ) deallocate ( NoahmpIO%vegfra ) ! vegetation fraction [] if ( allocated (NoahmpIO%tmn) ) deallocate ( NoahmpIO%tmn ) ! deep soil temperature [K] if ( allocated (NoahmpIO%xland) ) deallocate ( NoahmpIO%xland ) ! =2 ocean; =1 land/seaice diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90 index ada853d..84891f3 100644 --- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarInitMod.F90 @@ -41,6 +41,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO) if ( .not. allocated (NoahmpIO%zsoil) ) allocate ( NoahmpIO%zsoil (1:nsoil ) ) ! depth to soil interfaces [m] if ( .not. allocated (NoahmpIO%ivgtyp) ) allocate ( NoahmpIO%ivgtyp (its:ite ) ) ! vegetation type if ( .not. allocated (NoahmpIO%isltyp) ) allocate ( NoahmpIO%isltyp (its:ite ) ) ! soil type + if ( .not. allocated (NoahmpIO%isctyp) ) allocate ( NoahmpIO%isctyp (its:ite ) ) ! soil colour class if ( .not. allocated (NoahmpIO%vegfra) ) allocate ( NoahmpIO%vegfra (its:ite ) ) ! vegetation fraction [] if ( .not. allocated (NoahmpIO%tmn) ) allocate ( NoahmpIO%tmn (its:ite ) ) ! deep soil temperature [K] if ( .not. allocated (NoahmpIO%xland) ) allocate ( NoahmpIO%xland (its:ite ) ) ! =2 ocean; =1 land/seaice @@ -467,6 +468,7 @@ subroutine NoahmpIOVarInitDefault(NoahmpIO) NoahmpIO%ice = undefined_int NoahmpIO%ivgtyp = undefined_int NoahmpIO%isltyp = undefined_int + NoahmpIO%isctyp = undefined_int NoahmpIO%isnowxy = undefined_int NoahmpIO%coszen = undefined_real NoahmpIO%xlat = undefined_real diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90 index 0a3cd93..33f6765 100644 --- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpIOVarType.F90 @@ -68,6 +68,7 @@ module NoahmpIOVarType integer :: soil_update_steps ! number of model time steps to update soil process integer, allocatable, dimension(:) :: ivgtyp ! vegetation type integer, allocatable, dimension(:) :: isltyp ! soil type + integer, allocatable, dimension(:) :: isctyp ! soil colour class real(kind=kind_noahmp), allocatable, dimension(:) :: coszen ! cosine zenith angle real(kind=kind_noahmp), allocatable, dimension(:) :: xlat ! latitude [rad] real(kind=kind_noahmp), allocatable, dimension(:,:) :: dz8w ! thickness of atmo layers [m] diff --git a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90 b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90 index eb01ceb..7339ba2 100644 --- a/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90 +++ b/src/core_atmosphere/physics/physics_noahmp/drivers/mpas/NoahmpReadTableMod.F90 @@ -28,7 +28,7 @@ subroutine NoahmpReadTable(NoahmpIO) !------------------------------------------------------- integer, parameter :: MVT = 27 ! number of vegetation types integer, parameter :: MBAND = 2 ! number of radiation bands - integer, parameter :: MSC = 8 ! number of soil texture + integer, parameter :: MSC = 21 ! number of soil colour types integer, parameter :: MAX_SOILTYP = 30 ! max number of soil types integer, parameter :: NCROP = 5 ! number of crop types integer, parameter :: NSTAGE = 8 ! number of crop growth stages diff --git a/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL b/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL index c9d37c5..4ab3b66 100644 --- a/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL +++ b/src/core_atmosphere/physics/physics_noahmp/parameters/NoahmpTable.TBL @@ -405,13 +405,19 @@ / &noahmp_rad_parameters - !------------------------------------------------------------------------------ - ! soil color: 1 2 3 4 5 6 7 8 soil color index for soil albedo - !------------------------------------------------------------------------------ - ALBSAT_VIS = 0.15, 0.11, 0.10, 0.09, 0.08, 0.07, 0.06, 0.05 ! saturated soil albedo at visible band - ALBSAT_NIR = 0.30, 0.22, 0.20, 0.18, 0.16, 0.14, 0.12, 0.10 ! saturated soil albedo at NIR band - ALBDRY_VIS = 0.27, 0.22, 0.20, 0.18, 0.16, 0.14, 0.12, 0.10 ! dry soil albedo at visible band - ALBDRY_NIR = 0.54, 0.44, 0.40, 0.36, 0.32, 0.28, 0.24, 0.20 ! dry soil albedo at NIR band + !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ + ! soil color index for soil albedo + ! soil color: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 + !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ + ! saturated soil albedo at visible band + ALBSAT_VIS = 0.25, 0.23, 0.21, 0.20, 0.19, 0.18, 0.17, 0.16, 0.15, 0.14, 0.13, 0.12, 0.11, 0.10, 0.09, 0.08, 0.07, 0.06, 0.05, 0.04, 0.00 + ! saturated soil albedo at NIR band + ALBSAT_NIR = 0.50, 0.46, 0.42, 0.40, 0.38, 0.36, 0.34, 0.32, 0.30, 0.28, 0.26, 0.24, 0.22, 0.20, 0.18, 0.16, 0.14, 0.12, 0.10, 0.08, 0.00 + ! dry soil albedo at visible band + ALBDRY_VIS = 0.36, 0.34, 0.32, 0.31, 0.30, 0.29, 0.28, 0.27, 0.26, 0.25, 0.24, 0.23, 0.22, 0.20, 0.18, 0.16, 0.14, 0.12, 0.10, 0.08, 0.00 + ! dry soil albedo at NIR band + ALBDRY_NIR = 0.61, 0.57, 0.53, 0.51, 0.49, 0.48, 0.45, 0.43, 0.41, 0.39, 0.37, 0.35, 0.33, 0.31, 0.29, 0.27, 0.25, 0.23, 0.21, 0.16, 0.00 + !------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ ALBICE = 0.80, 0.55 ! albedo land ice: 1=vis, 2=nir ALBLAK = 0.60, 0.40 ! albedo frozen lakes: 1=vis, 2=nir OMEGAS = 0.8 , 0.4 ! two-stream parameter omega for snow diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml index ce1dc87..6b10574 100644 --- a/src/core_init_atmosphere/Registry.xml +++ b/src/core_init_atmosphere/Registry.xml @@ -420,6 +420,7 @@ + @@ -525,6 +526,7 @@ + @@ -796,6 +798,9 @@ + + diff --git a/src/core_init_atmosphere/mpas_geotile_manager.F b/src/core_init_atmosphere/mpas_geotile_manager.F index 04f9c60..b6b1663 100644 --- a/src/core_init_atmosphere/mpas_geotile_manager.F +++ b/src/core_init_atmosphere/mpas_geotile_manager.F @@ -98,6 +98,7 @@ module mpas_geotile_manager !> * isice = 24 !> * isurban = 1 !> * isoilwater = 14 + !> * islcolwater = 21 ! !----------------------------------------------------------------------- function mpas_geotile_mgr_init(mgr, path) result(ierr) @@ -119,7 +120,7 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr) integer, pointer :: tile_z_start, tile_z_end integer, pointer :: signed integer, pointer :: tile_bdr - integer, pointer :: iswater, islake, isice, isurban, isoilwater + integer, pointer :: iswater, islake, isice, isurban, isoilwater, islcolwater integer, pointer :: category_min, category_max integer :: err_level real (kind=RKIND), pointer :: dx ! Grid spacing in the x-direction @@ -188,7 +189,7 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr) ! ! If this is a categorical field, then check to see if it has category_max and category_min, - ! and then set the defaults of iswater, islake, isice, isurban and isoilwater + ! and then set the defaults of iswater, islake, isice, isurban, isoilwater and islcolwater ! call mpas_pool_get_config(mgr % pool, 'type', fieldType) if (fieldType == 'categorical') then @@ -219,12 +220,14 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr) isice => null() isurban => null() isoilwater => null() + islcolwater => null() call mpas_pool_get_config(mgr % pool, 'iswater', iswater) call mpas_pool_get_config(mgr % pool, 'islake', islake) call mpas_pool_get_config(mgr % pool, 'isice', isice) call mpas_pool_get_config(mgr % pool, 'isurban', isurban) call mpas_pool_get_config(mgr % pool, 'isoilwater', isoilwater) + call mpas_pool_get_config(mgr % pool, 'islcolwater', islcolwater) if (.not. associated(iswater)) then call mpas_pool_add_config(mgr % pool, 'iswater', 16) @@ -245,6 +248,10 @@ function mpas_geotile_mgr_init(mgr, path) result(ierr) if (.not. associated(isoilwater)) then call mpas_pool_add_config(mgr % pool, 'isoilwater', 14) endif + + if (.not. associated(isoilwater)) then + call mpas_pool_add_config(mgr % pool, 'islcolwater', 21) + end if endif ! diff --git a/src/core_init_atmosphere/mpas_init_atm_static.F b/src/core_init_atmosphere/mpas_init_atm_static.F index 09fe4c5..cceab7b 100644 --- a/src/core_init_atmosphere/mpas_init_atm_static.F +++ b/src/core_init_atmosphere/mpas_init_atm_static.F @@ -93,6 +93,7 @@ end subroutine interp_accumulation_function real (kind=RKIND) :: soilcomp_msgval = 255.0_RKIND ! Modified later based on index file for soilcomp integer, dimension(:), pointer :: lu_index integer, dimension(:), pointer :: soilcat_top + integer, dimension(:), pointer :: soilcol_idx integer, dimension(:), pointer :: nhs integer, dimension(:,:), allocatable:: ncat ! Landmask is used by the accumulation function for maxsnoalb and soilcomp, @@ -152,6 +153,7 @@ subroutine init_atm_static(mesh, dims, configs) integer, pointer :: isice_lu, iswater_lu integer :: iswater_soil + integer :: iswater_slcol integer, pointer :: nCells, nCellsSolve, nEdges, nVertices, maxEdges logical, pointer :: on_a_sphere real (kind=RKIND), pointer :: sphere_radius @@ -182,6 +184,7 @@ subroutine init_atm_static(mesh, dims, configs) real (kind=RKIND), pointer :: missing_value integer, dimension(:), pointer :: lu_index integer, dimension(:), pointer :: soilcat_top + integer, dimension(:), pointer :: soilcol_idx integer, dimension(:), pointer :: landmask integer, dimension(:), pointer :: bdyMaskCell character(len=StrKIND), pointer :: mminlu @@ -264,6 +267,7 @@ subroutine init_atm_static(mesh, dims, configs) call mpas_pool_get_array(mesh, 'isice_lu', isice_lu) call mpas_pool_get_array(mesh, 'iswater_lu', iswater_lu) call mpas_pool_get_array(mesh, 'soilcat_top', soilcat_top) + call mpas_pool_get_array(mesh, 'soilcol_idx', soilcol_idx) call mpas_pool_get_array(mesh, 'landmask', landmask) call mpas_pool_get_array(mesh, 'snoalb', snoalb) call mpas_pool_get_array(mesh, 'greenfrac', greenfrac) @@ -418,17 +422,28 @@ subroutine init_atm_static(mesh, dims, configs) supersample_fac=supersample_fac_30s) call mpas_log_write('--- end interpolate SOILCAT_TOP') +! +! Interpolate SOILCOL_IDX +! + geog_sub_path = 'soilcolour_30s/' + + call mpas_log_write('--- start interpolate SOILCOL_IDX') + call interp_soilcol(mesh, tree, trim(geog_data_path)//trim(geog_sub_path), iswater_slcol, & + supersample_fac=supersample_fac_30s) + call mpas_log_write('--- end interpolate SOILCOL_IDX') + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! KLUDGE TO FIX SOIL TYPE OVER ANTARCTICA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! where (lu_index == isice_lu) soilcat_top = 16 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! -! CORRECT INCONSISTENT SOIL AND LAND USE DATA +! CORRECT INCONSISTENT SOIL TYPE, SOIL COLOUR AND LAND USE DATA !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do iCell = 1,nCells - if (lu_index(iCell) == iswater_lu .or. & - soilcat_top(iCell) == iswater_soil) then + if (lu_index(iCell) == iswater_lu .or. & + soilcat_top(iCell) == iswater_soil .or. & + soilcol_idx(iCell) == iswater_slcol ) then if (lu_index(iCell) /= iswater_lu) then call mpas_log_write('Turning lu_index into water at $i', intArgs=(/iCell/)) lu_index(iCell) = iswater_lu @@ -437,6 +452,10 @@ subroutine init_atm_static(mesh, dims, configs) call mpas_log_write('Turning soilcat_top into water at $i', intArgs=(/iCell/)) soilcat_top(iCell) = iswater_soil end if + if (soilcol_idx(iCell) /= iswater_slcol) then + call mpas_log_write('Turning soilcol_idx into water at $i', intArgs=(/iCell/)) + soilcol_idx(iCell) = iswater_slcol + end if end if end do @@ -1907,6 +1926,89 @@ subroutine interp_soilcat(mesh, kdtree, geog_data_path, iswater_soil, supersampl end subroutine interp_soilcat + !*********************************************************************** + ! + ! routine interp_soilcol + ! + !> \brief Interpolate soil colour indices + !> \author Marcos Longo + !> \date 15 October 2025 + !> \details + !> This sub-routine is based on interp_soilcat. The procedure interpolates soil + !> colour indices by using the init_atm_map_static_data routine and then by + !> accumulating the pixel values into each cell using + !> category_interp_accumulation. + !> + !> Variable mesh should be an mpas_pool that contains nCells and soilcol_idx; + !> kdtree should be an initialised mpas_kd_type tree with dimensions + !> (xCell, yCell, zCell); and geog_data_path should be the path to the + !> soil colour index data set. The values used by this data set to flag + !> pixels with no data (e.g., pixels fully covered with water or ice). Some + !> data sets may not have the no data flag, but we harmonise this data set + !> with the land use class mapping. + !----------------------------------------------------------------------- + subroutine interp_soilcol(mesh, kdtree, geog_data_path, iswater_slcol, supersample_fac) + + implicit none + + ! Input variables + type (mpas_pool_type), intent(inout) :: mesh + type (mpas_kd_type), pointer, intent(in) :: kdtree + character (len=*), intent(in) :: geog_data_path + integer, intent(out) :: iswater_slcol + integer, intent(in), optional :: supersample_fac + + ! Local variables + type (mpas_geotile_mgr_type) :: mgr + integer, pointer :: nCells + integer, pointer :: iswater_slcol_ptr + + real (kind=RKIND), pointer :: scalefactor + + integer :: iCell + integer :: ierr + + ierr = mgr % init(trim(geog_data_path)) + if (ierr /= 0) then + call mpas_log_write("Error occured initalising interpolation for "//trim(geog_data_path), messageType=MPAS_LOG_CRIT) + return + end if + + call mpas_pool_get_dimension(mesh, 'nCells', nCells) + call mpas_pool_get_array(mesh, 'soilcol_idx', soilcol_idx) + call mpas_pool_get_config(mgr % pool, 'scale_factor', scalefactor ) + call mpas_pool_get_config(mgr % pool, 'category_min', category_min ) + call mpas_pool_get_config(mgr % pool, 'category_max', category_max ) + call mpas_pool_get_config(mgr % pool, 'islcolwater' , iswater_slcol_ptr) + + iswater_slcol = iswater_slcol_ptr + + allocate(ncat(category_min:category_max, nCells)) + ncat(:,:) = 0 + + call init_atm_map_static_data(mesh, mgr, kdtree, categorical_interp_criteria, categorical_interp_accumulation, & + supersample_fac=supersample_fac) + + do iCell = 1, nCells + ! Variable ncat will tally the classes occurring in the grid cell. We select the class that has the + ! highest count as the value representing the grid cell. We use maxloc for this, however, maxloc returns + ! the index corresponding to the maximum value starting from 1. We shift the result to make sure the + ! category value correctly maps between category_min and category_max. + soilcol_idx(iCell) = maxloc(ncat(:,iCell), dim=1) - 1 + category_min + end do + deallocate(ncat) + + ierr = mgr % finalize() + if (ierr /= 0) then + call mpas_log_write("Error occured finalising interpolation for "//trim(geog_data_path), messageType=MPAS_LOG_CRIT) + return + end if + + nullify(category_min) + nullify(category_max) + + end subroutine interp_soilcol + !*********************************************************************** ! ! routine derive_landmask diff --git a/src/framework/mpas_io.F b/src/framework/mpas_io.F index 09514a3..bf5a898 100644 --- a/src/framework/mpas_io.F +++ b/src/framework/mpas_io.F @@ -868,7 +868,7 @@ subroutine MPAS_io_inq_var(handle, fieldname, fieldtype, ndims, dimnames, dimsiz if (present(ierr)) ierr = MPAS_IO_ERR_BACKEND deallocate(new_fieldlist_node % fieldhandle) deallocate(new_fieldlist_node) - call mpas_log_write('Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) + call mpas_log_write('PIO_SUPPORT: Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) return end if !call mpas_log_write('Inquired about variable ID $i.', intArgs=(/new_fieldlist_node % fieldhandle % fieldid/) ) @@ -907,7 +907,7 @@ subroutine MPAS_io_inq_var(handle, fieldname, fieldtype, ndims, dimnames, dimsiz if (present(ierr)) ierr = MPAS_IO_ERR_BACKEND deallocate(new_fieldlist_node % fieldhandle) deallocate(new_fieldlist_node) - call mpas_log_write('Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) + call mpas_log_write('SMIOL_SUPPORT_A: Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) return end if new_fieldlist_node % fieldhandle % field_type = smiol_type @@ -955,7 +955,7 @@ subroutine MPAS_io_inq_var(handle, fieldname, fieldtype, ndims, dimnames, dimsiz deallocate(new_fieldlist_node % fieldhandle) deallocate(new_fieldlist_node) - call mpas_log_write('Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) + call mpas_log_write('SMIOL_SUPPORT_B: Variable ' // trim(fieldname) // ' not in input file.', MPAS_LOG_WARN) return end if new_fieldlist_node% fieldhandle % ndims = smiol_ndims