diff --git a/src/core_atmosphere/Registry.xml b/src/core_atmosphere/Registry.xml index 6378596797..f36de2d693 100644 --- a/src/core_atmosphere/Registry.xml +++ b/src/core_atmosphere/Registry.xml @@ -547,6 +547,18 @@ #ifdef DO_PHYSICS + + + + + + + + + + + + @@ -821,6 +833,18 @@ + + + + + + + + + + + + @@ -2079,11 +2103,6 @@ description="logical for configuration of fractional sea-ice" possible_values=".true. for sea-ice between 0 or 1; .false. for sea-ice equal to 0 or 1 (flag)."/> - - + + + + + + diff --git a/src/core_atmosphere/physics/mpas_atmphys_manager.F b/src/core_atmosphere/physics/mpas_atmphys_manager.F index 1b88522786..60b5725b4a 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_manager.F +++ b/src/core_atmosphere/physics/mpas_atmphys_manager.F @@ -39,9 +39,16 @@ module mpas_atmphys_manager !defines alarm to update the surface boundary conditions: character(len=*), parameter:: sfcbdyAlarmID = 'sfcbdy' -!defines alarm to update the background surface albedo and the greeness fraction: +!defines alarm to update the background surface albedo: + character(len=*), parameter:: albedoAlarmID = 'albedo' + +!defines alarm to update the background greeness fraction: character(len=*), parameter:: greenAlarmID = 'green' +!defines alarm to update the background leaf area index: + character(len=*), parameter:: laiAlarmID = 'lai' + + !defines alarm to update the ozone path length,the trace gas path length,the total emissivity, !and the total absorptivity in the "CAM" long-wave radiation codes. The default time interval !between updates is 6 hours and is set with config_camrad_abs_update (00:30:00). @@ -149,11 +156,14 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) !local pointers: logical,pointer:: config_frac_seaice, & config_o3climatology, & - config_sfc_albedo, & config_sst_update, & config_sstdiurn_update, & config_deepsoiltemp_update + character(len=StrKIND),pointer:: config_albedo_update, & + config_lai_update, & + config_greeness_update + character(len=StrKIND),pointer:: config_convection_scheme, & config_radt_lw_scheme, & config_radt_sw_scheme @@ -194,11 +204,14 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) call mpas_pool_get_config(domain%blocklist%configs,'config_frac_seaice' ,config_frac_seaice ) call mpas_pool_get_config(domain%blocklist%configs,'config_o3climatology' ,config_o3climatology ) - call mpas_pool_get_config(domain%blocklist%configs,'config_sfc_albedo' ,config_sfc_albedo ) call mpas_pool_get_config(domain%blocklist%configs,'config_sst_update' ,config_sst_update ) call mpas_pool_get_config(domain%blocklist%configs,'config_sstdiurn_update' ,config_sstdiurn_update ) call mpas_pool_get_config(domain%blocklist%configs,'config_deepsoiltemp_update',config_deepsoiltemp_update) + call mpas_pool_get_config(domain%blocklist%configs,'config_albedo_update' ,config_albedo_update ) + call mpas_pool_get_config(domain%blocklist%configs,'config_lai_update' ,config_lai_update ) + call mpas_pool_get_config(domain%blocklist%configs,'config_greeness_update' ,config_greeness_update ) + !update the current julian day and current year: currTime = mpas_get_clock_time(clock,MPAS_NOW,ierr) @@ -228,12 +241,28 @@ subroutine physics_timetracker(domain,dt,clock,itimestep,xtime_s) call mpas_pool_get_subpool(block%structs,'diag_physics',diag_physics) - !update the background surface albedo and greeness of vegetation: interpolation of input + !update the background surface albedo: interpolation of input + !monthly values to current day: + if(mpas_is_alarm_ringing(clock,albedoAlarmID,ierr=ierr)) then + call mpas_reset_clock_alarm(clock,albedoAlarmID,ierr=ierr) + call mpas_log_write('--- time to update background surface albedo') + call physics_update_surface_alb(timeStamp,mesh,sfc_input) + endif + + !update the background surface leaf area index: interpolation of input + !monthly values to current day: + if(mpas_is_alarm_ringing(clock,laiAlarmID,ierr=ierr)) then + call mpas_reset_clock_alarm(clock,laiAlarmID,ierr=ierr) + call mpas_log_write('--- time to update background leaf area index.') + call physics_update_surface_lai(timeStamp,mesh,sfc_input,diag_physics) + endif + + !update the background surface greenness of vegetation: interpolation of input !monthly values to current day: if(mpas_is_alarm_ringing(clock,greenAlarmID,ierr=ierr)) then call mpas_reset_clock_alarm(clock,greenAlarmID,ierr=ierr) - call mpas_log_write('--- time to update background surface albedo, greeness fraction.') - call physics_update_surface(timeStamp,config_sfc_albedo,mesh,sfc_input) + call mpas_log_write('--- time to update background greenness fraction of vegetation') + call physics_update_surface_gvf(timeStamp,mesh,sfc_input) endif !update surface boundary conditions with input sea-surface temperatures and fractional @@ -398,6 +427,8 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) config_radtsw_interval, & config_bucket_update, & config_camrad_abs_update, & + config_albedo_update, & + config_lai_update, & config_greeness_update logical,pointer:: config_sst_update @@ -434,11 +465,14 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) call mpas_pool_get_config(configs,'config_radtsw_interval' ,config_radtsw_interval ) call mpas_pool_get_config(configs,'config_bucket_update' ,config_bucket_update ) call mpas_pool_get_config(configs,'config_camrad_abs_update',config_camrad_abs_update) - call mpas_pool_get_config(configs,'config_greeness_update' ,config_greeness_update ) call mpas_pool_get_config(configs,'config_sst_update' ,config_sst_update ) call mpas_pool_get_config(configs,'config_frac_seaice' ,config_frac_seaice ) call mpas_pool_get_config(configs,'config_microp_re' ,config_microp_re ) + call mpas_pool_get_config(configs,'config_albedo_update' ,config_albedo_update ) + call mpas_pool_get_config(configs,'config_lai_update' ,config_lai_update ) + call mpas_pool_get_config(configs,'config_greeness_update' ,config_greeness_update ) + call mpas_pool_get_config(configs,'config_dt',config_dt) @@ -546,12 +580,32 @@ subroutine physics_run_init(configs,mesh,state,clock,stream_manager) endif -!set alarm for updating the background surface albedo and the greeness fraction: - call mpas_set_timeInterval(alarmTimeStep,timeString=config_greeness_update,ierr=ierr) - alarmStartTime = startTime - call mpas_add_clock_alarm(clock,greenAlarmID,alarmStartTime,alarmTimeStep,ierr=ierr) +!set alarm for updating the background surface albedo: + if(trim(config_albedo_update) /= "none") then + call mpas_set_timeInterval(alarmTimeStep,timeString=config_albedo_update,ierr=ierr) + alarmStartTime = startTime + call mpas_add_clock_alarm(clock,albedoAlarmID,alarmStartTime,alarmTimeStep,ierr=ierr) + if(ierr /= 0) & + call physics_error_fatal('subroutine physics_init: error creating alarm albedo') + endif + +!set alarm for updating the background leaf area index: + if(trim(config_lai_update) /= "none") then + call mpas_set_timeInterval(alarmTimeStep,timeString=config_lai_update,ierr=ierr) + alarmStartTime = startTime + call mpas_add_clock_alarm(clock,laiAlarmID,alarmStartTime,alarmTimeStep,ierr=ierr) + if(ierr /= 0) & + call physics_error_fatal('subroutine physics_init: error creating alarm lai') + endif + +!set alarm for updating the background greenness pf vegetation: + if(trim(config_greeness_update) /= "none") then + call mpas_set_timeInterval(alarmTimeStep,timeString=config_greeness_update,ierr=ierr) + alarmStartTime = startTime + call mpas_add_clock_alarm(clock,greenAlarmID,alarmStartTime,alarmTimeStep,ierr=ierr) if(ierr /= 0) & - call physics_error_fatal('subroutine physics_init: error creating alarm greeness') + call physics_error_fatal('subroutine physics_init: error creating alarm greenness') + endif !set alarm for updating the surface boundary conditions: if (config_sst_update) then diff --git a/src/core_atmosphere/physics/mpas_atmphys_update_surface.F b/src/core_atmosphere/physics/mpas_atmphys_update_surface.F index 6e2057d5cb..eeb02de42d 100644 --- a/src/core_atmosphere/physics/mpas_atmphys_update_surface.F +++ b/src/core_atmosphere/physics/mpas_atmphys_update_surface.F @@ -19,7 +19,9 @@ module mpas_atmphys_update_surface private public:: physics_update_sst, & physics_update_sstskin, & - physics_update_surface, & + physics_update_surface_alb, & + physics_update_surface_lai, & + physics_update_surface_gvf, & physics_update_deepsoiltemp @@ -29,7 +31,9 @@ module mpas_atmphys_update_surface ! ! subroutines in mpas_atmphys_update_surface: ! ------------------------------------------- -! physics_update_surface : update the surface albedo and greeness fraction. +! physics_update_surface_alb : update the surface albedo. +! physics_update_surface_lai : update the surface leaf area index. +! physics_update_surface_gvf : update the surface greeness fraction. ! physics_update_sst : update the sea-surface temperatures. ! physics_update_sstskin : add a diurnal cycle to the sea-surface temperatures. ! physics_update_deepsoiltemp: update the deep soil temperatures. @@ -51,27 +55,23 @@ module mpas_atmphys_update_surface !================================================================================================================= - subroutine physics_update_surface(current_date,config_sfc_albedo,mesh,sfc_input) + subroutine physics_update_surface_alb(current_date,mesh,sfc_input) !================================================================================================================= !input variables: type(mpas_pool_type),intent(in):: mesh character(len=*),intent(in):: current_date - logical,intent(in):: config_sfc_albedo !inout variables: type(mpas_pool_type),intent(inout):: sfc_input !local pointers: -!logical,pointer:: config_sfc_albedo integer,pointer:: nCellsSolve integer,dimension(:),pointer:: landmask real(kind=RKIND),dimension(:) ,pointer:: sfc_albbck real(kind=RKIND),dimension(:,:),pointer:: albedo12m - real(kind=RKIND),dimension(:) ,pointer:: vegfra,shdmin,shdmax - real(kind=RKIND),dimension(:,:),pointer:: greenfrac !local variables: integer:: iCell @@ -84,23 +84,85 @@ subroutine physics_update_surface(current_date,config_sfc_albedo,mesh,sfc_input) call mpas_pool_get_array(sfc_input,'albedo12m' , albedo12m ) call mpas_pool_get_array(sfc_input,'sfc_albbck', sfc_albbck) - call mpas_pool_get_array(sfc_input,'greenfrac' , greenfrac ) - call mpas_pool_get_array(sfc_input,'vegfra' , vegfra ) - call mpas_pool_get_array(sfc_input,'shdmin' , shdmin ) - call mpas_pool_get_array(sfc_input,'shdmax' , shdmax ) - !updates the surface background albedo for the current date as a function of the monthly-mean -!surface background albedo valid on the 15th day of the month, if config_sfc_albedo is true: - if(config_sfc_albedo) then +!surface background albedo valid on the 15th day of the month: + call monthly_interp_to_date(nCellsSolve,current_date,albedo12m,sfc_albbck) - call monthly_interp_to_date(nCellsSolve,current_date,albedo12m,sfc_albbck) + do iCell = 1, nCellsSolve + sfc_albbck(iCell) = sfc_albbck(iCell) / 100. + if(landmask(iCell) .eq. 0) sfc_albbck(iCell) = 0.08 + enddo - do iCell = 1, nCellsSolve - sfc_albbck(iCell) = sfc_albbck(iCell) / 100. - if(landmask(iCell) .eq. 0) sfc_albbck(iCell) = 0.08 - enddo + end subroutine physics_update_surface_alb - endif +!================================================================================================================= + subroutine physics_update_surface_lai(current_date,mesh,sfc_input,diag_physics) +!================================================================================================================= + +!input variables: + type(mpas_pool_type),intent(in):: mesh + character(len=*),intent(in):: current_date + +!inout variables: + type(mpas_pool_type),intent(inout):: sfc_input + type(mpas_pool_type),intent(inout):: diag_physics + +!local pointers: + + integer,pointer:: nCellsSolve + + real(kind=RKIND),dimension(:,:),pointer:: lai12m + real(kind=RKIND),dimension(:) ,pointer:: lai + +!local variables: + integer:: iCell + +!----------------------------------------------------------------------------------------------------------------- + + call mpas_pool_get_dimension(mesh,'nCellsSolve',nCellsSolve) + + call mpas_pool_get_array(sfc_input,'lai12m' , lai12m ) + call mpas_pool_get_array(diag_physics,'lai ' , lai ) + +!updates the leaf area index for the current date as a function of the monthly-mean lai +!valid on the 15th day of the month. + call monthly_interp_to_date(nCellsSolve,current_date,lai12m,lai) + + do iCell = 1,nCellsSolve + lai(iCell) = lai(iCell) / 10. + enddo + + end subroutine physics_update_surface_lai + +!================================================================================================================= + subroutine physics_update_surface_gvf(current_date,mesh,sfc_input) +!================================================================================================================= + +!input variables: + type(mpas_pool_type),intent(in):: mesh + character(len=*),intent(in):: current_date + +!inout variables: + type(mpas_pool_type),intent(inout):: sfc_input + +!local pointers: + + integer,pointer:: nCellsSolve + + real(kind=RKIND),dimension(:) ,pointer:: vegfra,shdmin,shdmax + real(kind=RKIND),dimension(:,:),pointer:: greenfrac + +!local variables: + integer:: iCell + +!----------------------------------------------------------------------------------------------------------------- + + call mpas_pool_get_dimension(mesh,'nCellsSolve',nCellsSolve) + + call mpas_pool_get_array(sfc_input,'greenfrac' , greenfrac ) + call mpas_pool_get_array(sfc_input,'vegfra' , vegfra ) + call mpas_pool_get_array(sfc_input,'shdmin' , shdmin ) + call mpas_pool_get_array(sfc_input,'shdmax' , shdmax ) !updates the green-ness fraction for the current date as a function of the monthly-mean green- !ness valid on the 15th day of the month. get the min/max for each cell for the monthly green- @@ -108,7 +170,7 @@ subroutine physics_update_surface(current_date,config_sfc_albedo,mesh,sfc_input) call monthly_interp_to_date(nCellsSolve,current_date,greenfrac,vegfra) call monthly_min_max(nCellsSolve,greenfrac,shdmin,shdmax) - end subroutine physics_update_surface + end subroutine physics_update_surface_gvf !================================================================================================================= subroutine physics_update_sst(dminfo,config_frac_seaice,mesh,sfc_input,diag_physics) diff --git a/src/core_init_atmosphere/Registry.xml b/src/core_init_atmosphere/Registry.xml index b107d16953..2cd7e2d625 100644 --- a/src/core_init_atmosphere/Registry.xml +++ b/src/core_init_atmosphere/Registry.xml @@ -176,6 +176,11 @@ description="The climatological monthly vegetation fraction dataset to use (case 7 only)" possible_values="`MODIS' or `NCEP'"/> + + + @@ -537,6 +543,7 @@ + @@ -820,6 +827,9 @@ + + diff --git a/src/core_init_atmosphere/mpas_init_atm_static.F b/src/core_init_atmosphere/mpas_init_atm_static.F index b9d289d8b6..6eab67c9aa 100644 --- a/src/core_init_atmosphere/mpas_init_atm_static.F +++ b/src/core_init_atmosphere/mpas_init_atm_static.F @@ -126,6 +126,7 @@ subroutine init_atm_static(mesh, dims, configs) character(len=StrKIND), pointer :: config_vegfrac_data character(len=StrKIND), pointer :: config_albedo_data character(len=StrKIND), pointer :: config_maxsnowalbedo_data + character(len=StrKIND), pointer :: config_lai_data logical, pointer :: config_noahmp_static character(len=StrKIND+1) :: geog_data_path ! same as config_geog_data_path, but guaranteed to have a trailing slash character(len=StrKIND+1) :: geog_sub_path ! subdirectory names in config_geog_data_path, with trailing slash @@ -179,6 +180,8 @@ subroutine init_atm_static(mesh, dims, configs) real (kind=RKIND), dimension(:,:), pointer :: greenfrac real (kind=RKIND), dimension(:,:), pointer :: albedo12m integer (kind=I8KIND), dimension(:,:), pointer :: albedo12m_int + real (kind=RKIND), dimension(:,:), pointer :: lai12m + integer (kind=I8KIND), dimension(:,:), pointer :: lai12m_int real (kind=RKIND) :: fillval real (kind=RKIND), pointer :: missing_value integer, dimension(:), pointer :: lu_index @@ -217,6 +220,7 @@ subroutine init_atm_static(mesh, dims, configs) call mpas_pool_get_config(configs, 'config_topo_data', config_topo_data) call mpas_pool_get_config(configs, 'config_vegfrac_data', config_vegfrac_data) call mpas_pool_get_config(configs, 'config_albedo_data', config_albedo_data) + call mpas_pool_get_config(configs, 'config_lai_data', config_lai_data) call mpas_pool_get_config(configs, 'config_maxsnowalbedo_data', config_maxsnowalbedo_data) call mpas_pool_get_config(configs, 'config_supersample_factor', supersample_fac) call mpas_pool_get_config(configs, 'config_30s_supersample_factor', supersample_fac_30s) @@ -270,6 +274,7 @@ subroutine init_atm_static(mesh, dims, configs) call mpas_pool_get_array(mesh, 'snoalb', snoalb) call mpas_pool_get_array(mesh, 'greenfrac', greenfrac) call mpas_pool_get_array(mesh, 'albedo12m', albedo12m) + call mpas_pool_get_array(mesh, 'lai12m', lai12m) call mpas_pool_get_array(mesh, 'shdmin', shdmin) call mpas_pool_get_array(mesh, 'shdmax', shdmax) @@ -1225,6 +1230,150 @@ subroutine init_atm_static(mesh, dims, configs) supersample_fac=supersample_fac_30s) call mpas_log_write('--- end interpolate SOILCL4') + endif + +! +! Interpolate LAI12M +! + if (trim(config_lai_data) == 'MODIS') then + + call mpas_log_write('Using MODIS FPAR 30-arc-second data for climatological monthly leaf area index') + + geog_sub_path = 'lai_modis_30s/' + + ierr = mgr % init(trim(geog_data_path)//trim(geog_sub_path)) + if (ierr /= 0) then + call mpas_log_write('Error occurred when initalizing the interpolation of monthly leaf area index (lai12m)', & + messageType=MPAS_LOG_CRIT) + endif + + call mpas_pool_get_config(mgr % pool, 'tile_bdr', tile_bdr) + call mpas_pool_get_config(mgr % pool, 'tile_x', tile_nx) + call mpas_pool_get_config(mgr % pool, 'tile_y', tile_ny) + call mpas_pool_get_config(mgr % pool, 'tile_z', tile_nz) + call mpas_pool_get_config(mgr % pool, 'missing_value', missing_value) + + allocate(nhs(nCells)) + allocate(lai12m_int(tile_nz, nCells)) + nhs(:) = 0 + lai12m(:,:) = 0.0 + lai12m_int(:,:) = 0_I8KIND + fillval = 0.0 + + do iCell = 1, nCells + if (nhs(iCell) == 0) then + tile => null() + ierr = mgr % get_tile(latCell(iCell), lonCell(iCell), tile) + if (ierr /= 0 .or. .not. associated(tile)) then + call mpas_log_write('Could not get tile that contained cell $i', intArgs=(/iCell/), messageType=MPAS_LOG_CRIT) + end if + + ierr = mgr % push_tile(tile) + if (ierr /= 0) then + call mpas_log_write("Error pushing this tile onto the stack: "//trim(tile % fname), messageType=MPAS_LOG_CRIT) + end if + end if + + do while (.not. mgr % is_stack_empty()) + tile => mgr % pop_tile() + if (tile % is_processed) then + cycle + end if + + call mpas_log_write('Processing tile: '//trim(tile % fname)) + + all_pixels_mapped_to_halo_cells = .true. + + do j = tile_bdr + 1, tile_ny + tile_bdr, 1 + do i = tile_bdr + 1, tile_nx + tile_bdr, 1 + + call mgr % tile_to_latlon(tile, j, i, lat_pt, lon_pt) + call mpas_latlon_to_xyz(xPixel, yPixel, zPixel, sphere_radius, lat_pt, lon_pt) + call mpas_kd_search(tree, (/xPixel, yPixel, zPixel/), res, max_distance=max_kdtree_distance2) + + if (.not. associated(res)) cycle + + ! + ! This field only matters for land cells, and for all but the outermost boundary cells, + ! we can safely assume that the nearest model grid cell contains the pixel (else, a different + ! cell would be nearest) + ! + if (landMask(res % id) == 1 .and. bdyMaskCell(res % id) < nBdyLayers) then + do k = 1, tile_nz + if (tile % tile(i, j, k) == missing_value) then + i8val = int(fillval, kind=I8KIND) + else + i8val = int(tile % tile(i,j,k), kind=I8KIND) + end if + lai12m_int(k, res % id) = lai12m_int(k, res % id) + i8val + end do + nhs(res % id) = nhs(res % id) + 1 + + if (res % id <= nCellsSolve) then + all_pixels_mapped_to_halo_cells = .false. + end if + ! For outermost land cells, additional work is needed to verify that the pixel + ! actually lies within the nearest cell + else if (landMask(res % id) == 1) then + if (mpas_in_cell(xPixel, yPixel, zPixel, xCell(res % id), yCell(res % id), zCell(res % id), & + nEdgesOnCell(res % id), verticesOnCell(:,res % id), xVertex, yVertex, zVertex)) then + do k = 1, tile_nz + if (tile % tile(i, j, k) == missing_value) then + i8val = int(fillval, kind=I8KIND) + else + i8val = int(tile % tile(i,j,k), kind=I8KIND) + end if + lai12m_int(k, res % id) = lai12m_int(k, res % id) + i8val + end do + nhs(res % id) = nhs(res % id) + 1 + + if (res % id <= nCellsSolve) then + all_pixels_mapped_to_halo_cells = .false. + end if + end if + end if + end do + end do + + tile % is_processed = .true. + + if (.not. all_pixels_mapped_to_halo_cells) then + ierr = mgr % push_neighbors(tile) + if (ierr /= 0) then + call mpas_log_write("Error pushing the tile neighbors of: "//trim(tile%fname), messageType=MPAS_LOG_CRIT) + end if + end if + + end do + end do + + do iCell = 1, nCells + ! For land points that have no overlap with valid data, and for water points, + ! just use the fill value... + if (nhs(iCell) == 0) then + lai12m(:,iCell) = fillval + else + lai12m(:,iCell) = real(real(lai12m_int(:,iCell), kind=R8KIND) / real(nhs(iCell), kind=R8KIND), kind=RKIND) + end if + end do + + deallocate(nhs) + deallocate(lai12m_int) + + ierr = mgr % finalize() + if (ierr /= 0) then + call mpas_log_write('Error occurred when finalizing the interpolation of monthly leaf area index (lai)', & + messageType=MPAS_LOG_CRIT) + endif + else + + call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR) + call mpas_log_write('Invalid monthly leaf area index dataset '''//trim(config_albedo_data) & + //''' selected for config_lai_data', messageType=MPAS_LOG_ERR) + call mpas_log_write(' Possible options are: ''MODIS''', messageType=MPAS_LOG_ERR) + call mpas_log_write('*****************************************************************', messageType=MPAS_LOG_ERR) + call mpas_log_write('Please correct the namelist.', messageType=MPAS_LOG_CRIT) + end if !