From ee6283aa5034c3aabdc74809aba186bd04d2bb05 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Wed, 3 Aug 2022 09:11:43 -0400 Subject: [PATCH 01/43] add snow albedo to catch_params.nc4 --- .../Utils/Raster/make_bcs | 2 + .../Utils/Raster/mkCatchParam.F90 | 39 ++- .../Utils/Raster/mod_process_hres_data.F90 | 251 ++++++++++++++++++ .../Utils/Raster/rmTinyCatchParaMod.F90 | 9 + 4 files changed, 300 insertions(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs index a174aff89..2d7935fa9 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs @@ -105,6 +105,7 @@ echo " ${C2}ICA : Icarus (/discover/nobackup/ltakacs/bcs/Icarus/)" echo " ${C2}NL3 : Icarus-NLv3 (/discover/nobackup/ltakacs/bcs/Icarus-NLv3/)" echo " ${C2}NL4 : NLv4 [SMAP] (/discover/nobackup/projects/gmao/smap/bcs_NLv4/NLv4/)" echo " ${C2}NL5 : NLv5 [SMAP]" +echo " ${C2}NL6 : NLv6 [SMAP with snow albedo from MODIS]" echo " ${C2}DEV : Development version${CR}" echo " " if ( $HELPMODE != YES ) then @@ -132,6 +133,7 @@ if ( $HELPMODE != YES ) then $dummy == 'NL3' | \ $dummy == 'NL4' | \ $dummy == 'NL5' | \ + $dummy == 'NL6' | \ $dummy == 'DEV') then set lbcsv = $dummy else if ( $dummy == '' ) then diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index fe035a9e5..f6f26e7b0 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -564,7 +564,7 @@ PROGRAM mkCatchParam if( F25Tag) call soil_para_high (nc,nr,regrid,gridnamer,F25Tag=F25Tag) if(.not.F25Tag) call soil_para_high (nc,nr,regrid,gridnamer) endif - if(SOILBCS=='HWSD') call soil_para_hwsd (nc,nr,gridnamer) + if(SOILBCS=='HWSD') call soil_para_hwsd (nc,nr,gridnamer) write (log_file,'(a)')' Done.' else write (log_file,'(a,a)')' Using existing file.' @@ -653,6 +653,43 @@ PROGRAM mkCatchParam endif write (log_file,'(a)')' ' + ! Biljana + ! Creating snow_alb_param.dat + ! --------------------------------------------------------------------- + if(use_snow_albedo)then + tmpstring = 'Step 14: Snow albedo from MODIS' + ! here I should test if snow_alb is written into nc file + ! and call the soil_snow_alb subroutine only if it's not there + ! inquire(file='clsm/catch_params.nc4', exist=file_exists) + ! if (.not.file_exists) then + ! write (log_file,'(a)')' clsm/catch_params.nc4 is missing. Something is wrong. STOP!' + ! stop + ! else + ! inquire for snow albedo variable + + ! if not present, inititate + write (log_file,'(a)')' Loading snow albedo...' + call soil_snow_alb (nc,nr,gridnamer) + write (log_file,'(a)')' Done.' + + ! if present and loaded, do nothing + + ! endif + + !fname_tmp = 'clsm/snow_alb_param.dat' + !write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + !inquire(file=trim(fname_tmp), exist=file_exists) + !if (.not.file_exists) then + ! write (log_file,'(a)')' Creating file...' + ! call soil_snow_alb (nc,nr,gridnamer) + ! write (log_file,'(a)')' Done.' + !else + ! write (log_file,'(a,a)')' Using existing file.' + !endif + write (log_file,'(a)')' ' + endif ! if use_snow_albedo + ! Biljana -- done creating snow albedo file + ! inquire(file='clsm/irrig.dat', exist=file_exists) ! if (.not.file_exists) call create_irrig_params (nc,nr,gridnamer) ! write (log_file,'(a)')'Done computing irrigation model parameters ...............13' diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index a75036635..de11d794e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -26,6 +26,7 @@ MODULE process_hres_data use leap_year use MAPL_ConstantsMod use lsm_routines, ONLY: sibalb +USE IFPORT ! Biljana #if defined USE_EXTERNAL_FINDLOC use findloc_mod, only: findloc @@ -38,6 +39,7 @@ MODULE process_hres_data private public :: soil_para_hwsd,hres_lai,hres_gswp2, merge_lai_data, grid2tile_modis6 +public :: soil_snow_alb ! Biljana public :: modis_alb_on_tiles_high,modis_scale_para_high,hres_lai_no_gswp public :: histogram, create_mapping, esa2mosaic , esa2clm public :: grid2tile_ndep_t2m_alb, CREATE_ROUT_PARA_FILE, map_country_codes, get_country_codes @@ -2986,6 +2988,254 @@ END SUBROUTINE hres_gswp2 !---------------------------------------------------------------------- + SUBROUTINE soil_snow_alb (nx,ny,gfiler) ! Biljana + +! Implement snow albedo calculated from MOIDS 22-year climatology. Store snow albedo +! values in clsm/catch_params.nc4 +! Biljana Orescanin July 2022, SSAI@NASA + + implicit none + integer, intent (in) :: nx, ny + character(*) :: gfiler + integer,allocatable,target,dimension (:,:) :: tile_id + + character*200 :: fname + character*2 :: vv,hh + integer :: n,maxcat,i,j,k,ncid,status + real,allocatable,dimension(:) :: min_lon,max_lon,min_lat,max_lat,snw_alb + integer(kind=4),parameter :: xdim = 1200, ydim = 1200 + real,parameter :: alb_res=10.0/1200.0 + real,dimension(xdim) :: lon_alb + real,dimension(ydim) :: lat_alb + real,dimension(xdim,ydim) :: stch_snw_alb_tmp + real,dimension(36,18,xdim,ydim) :: stch_snw_alb + real :: minlon,maxlon,minlat,maxlat,pad_lon,pad_lat + real :: sno_alb_cnt, sno_alb_sum,sno_alb_cnt2,sno_alb_sum2 + integer :: vvtil_min,hhtil_min,vvtil_max,hhtil_max,hhtil,vvtil + integer :: tindex1,pfaf1 + integer(kind=4) :: dummy,VarID,varid1,varid2,varid3 + integer(kind=4) :: imin,imax,jmin,jmax + integer(kind=4) :: imin2,imax2,jmin2,jmax2,count_init_invalid + logical :: file_exists + + ! Read number of catchment-tiles (maxcat) from catchment.def file + fname='clsm/catchment.def' + open (10,file=fname,status='old',action='read',form='formatted') + read(10,*) maxcat + + ! read min/max lat/lons, so those can be used to locate + ! snow albedo grids in the stitched MODIS albedo file + allocate (min_lon(1:maxcat)) + allocate (min_lat(1:maxcat)) + allocate (max_lon(1:maxcat)) + allocate (max_lat(1:maxcat)) + allocate (snw_alb(1:maxcat)) + + ! before populating, set all snow albedo values to missing + snw_alb(:)=-9999.0 + + do n = 1, maxcat + read (10,*) tindex1,pfaf1,minlon,maxlon,minlat,maxlat + min_lon(n) = minlon + max_lon(n) = maxlon + min_lat(n) = minlat + max_lat(n) = maxlat + end do + + close (10,status='keep') + + ! Read tile-id raster file + allocate(tile_id(1:nx,1:ny)) + + fname=trim(gfiler)//'.rst' + open (10,file=fname,status='old',action='read', & + form='unformatted',convert='little_endian') + + do j=1,ny + read(10)tile_id(:,j) + end do + + close (10,status='keep') + + !------------ Get the information on snow albedo ----- + ! ----------- The information on snow albedo is stored in 10x10 30-arcsec files. Read in this + ! information. Then loop over tiles to find a corresponding snow albedo mean + + ! read in all 10x10deg snow albedo files into a single [36,18,1200,1200] array + do hhtil=1,36 ! loop over all horziontal input files + do vvtil=1,18 ! loop over all vertical input files + + write(vv,'(i2.2)') vvtil + write(hh,'(i2.2)') hhtil + + fname = '/discover/nobackup/borescan/tools/idl/01_snow_fraction/06_modis_nsidc/'// & + '/data/data_out/snow_alb_all_08_Top99th_percentile_MOD10A1.A_30arcsec_'// & + '2000_2022_H'//hh//'V'//vv//'.nc' + + ! Open the file. NF90_NOWRITE tells netCDF we want read-only access to the file. + status=NF_OPEN(trim(fname),NF_NOWRITE, ncid) ; VERIFY_(STATUS) + ! Get the varid of the data variable, based on its name. + status=NF_INQ_VARID(ncid,'Snow_Albedo',VarID1) ; VERIFY_(STATUS) + status=NF_INQ_VARID(ncid,'lon' ,VarID2) ; VERIFY_(STATUS) + status=NF_INQ_VARID(ncid,'lat' ,VarID3) ; VERIFY_(STATUS) + ! Read the data. + status=NF_GET_VARA_REAL(ncid,VarID1,(/1,1/),(/xdim,ydim/),stch_snw_alb_tmp) ; VERIFY_(STATUS) + status=NF_GET_VARA_REAL(ncid,VarID2,(/1/) ,(/xdim/) ,lon_alb) ; VERIFY_(STATUS) + status=NF_GET_VARA_REAL(ncid,VarID3,(/1/) ,(/ydim/) ,lat_alb) ; VERIFY_(STATUS) + ! Close the file, freeing all resources. + status=NF_CLOSE(ncid); VERIFY_(STATUS) + + ! store into large aray + stch_snw_alb(hhtil,vvtil,:,:)=stch_snw_alb_tmp + + enddo + enddo + + ! open the file to write snow albedo output in +! fname ='clsm/snow_alb_param.dat' +! open(11,file=trim(fname),form='formatted',status='unknown',action = 'write') + + ! loop over tiles + print*, 'Starting tile loop for snow albedo. Biljana' + count_init_invalid=0 ! counter for non-valid snow albedo after matching tile size + + do n = 1, maxcat ! loop over tile + + ! set the current tile snow albedo to missing. Then start calculations to see if not missing. + snw_alb(n)=-9999.0 + + ! set sums and counts to zero + sno_alb_sum=0. + sno_alb_cnt=0. + sno_alb_sum2=0. + sno_alb_cnt2=0. + + ! This tile has min/max lat/lon info. Use this info to identify which 10x10deg + ! snow albedo file(s) to read in (and then loop over these files). + ! Using ceiling and floor for max and min range, the "halo" approach is + ! implemented. + vvtil_min= floor((min_lat(n)+ 90.0)/10.) + hhtil_min= floor((min_lon(n)+180.0)/10.) + vvtil_max=ceiling((max_lat(n)+ 90.0)/10.) + hhtil_max=ceiling((max_lon(n)+180.0)/10.) + + ! make sure vv's and hh's are within the range + ! if min>max, swap them + if (vvtil_min .gt. vvtil_max) then + dummy =vvtil_min + vvtil_min=vvtil_max + vvtil_max=dummy + endif + if (hhtil_min .gt. hhtil_max) then + dummy =hhtil_min + hhtil_min=hhtil_max + hhtil_max=dummy + endif + + ! if beyond the range, bring them back + vvtil_min=max(vvtil_min,1) + vvtil_max=min(vvtil_max,18) + hhtil_min=max(hhtil_min,1) + hhtil_max=min(hhtil_max,36) + + do hhtil=hhtil_min,hhtil_max ! loop over all horziontal input files + do vvtil=vvtil_min,vvtil_max ! loop over all vertical input files + + ! find indices covered by the tile + imin=floor((min_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) + imax=floor((max_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) + jmin=floor((min_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) + jmax=floor((max_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) + ! make sure to stay within the range + imin=max(imin,1) + imax=min(imax,xdim) + jmin=max(jmin,1) + jmax=min(jmax,ydim) + + ! sum snow albedo values and counts for current tile corresponding indices + sno_alb_sum= sno_alb_sum + & + sum(stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax), & + stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & + stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax).le.1.0) + + sno_alb_cnt= sno_alb_cnt + & + count(stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & + stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax).le.1.0) + + end do ! vvtil + end do ! hhtil + ! get mean snow albedo over the tile + snw_alb(n) = sno_alb_sum / max(1.0,sno_alb_cnt) + if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=-9999.0 !1.E15 + + ! if no valid solition, and if tile size is smaller than the snow albedo grid box, + ! then expand the search areaby 1 tile in each direction + + ! size of a tile (in both directions) + pad_lon=(max_lon(n)-min_lon(n)) + pad_lat=(max_lat(n)-min_lat(n)) + + if (snw_alb(n) .le. 0.0 .and. (pad_lon .lt. alb_res .or. pad_lat .lt. alb_res)) then + + count_init_invalid=count_init_invalid+1 + + do hhtil=hhtil_min,hhtil_max ! loop over all horziontal input files + do vvtil=vvtil_min,vvtil_max ! loop over all vertical input files + + ! find indices of snow albedo array corresponding to the current tile + imin2=floor((min_lon(n)-pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) + imax2=floor((max_lon(n)+pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) + jmin2=floor((min_lat(n)-pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) + jmax2=floor((max_lat(n)+pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) + imin2=max(imin2,1) + imax2=min(imax2,xdim) + jmin2=max(jmin2,1) + jmax2=min(jmax2,ydim) + + ! sum snow albedo values and counts for current tile corresponding indices + sno_alb_sum2= sno_alb_sum2 + & + sum(stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2), & + stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & + stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2).le.1.0) + + sno_alb_cnt2= sno_alb_cnt2 + & + count(stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & + stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2).le.1.0) + + end do ! vvtil + end do ! hhtil + + snw_alb(n) = sno_alb_sum2 / max(1.0,sno_alb_cnt2) + if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=-9999.0 !1.E15 + + endif + + ! write the current tile value into the filr +! write (11,'(i10,i8,f13.4)') tindex1,pfaf1,snw_alb(n) + + end do ! n-loop over tiles + + ! write snow albedo into clsm/catch_params.nc4 + inquire(file='clsm/catch_params.nc4', exist=file_exists) + + if(file_exists) then + status = NF_OPEN ('clsm/catch_params.nc4', NF_WRITE, ncid ) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCID,NC_VarID(NCID,'SNOWALB'),(/1/),(/maxcat/),real(snw_alb)) ; VERIFY_(STATUS) + STATUS = NF_CLOSE (NCID) ; VERIFY_(STATUS) + endif + +! ! close the output file +! write (11,'(a)')' ' +! write (11,'(a)')'TileIndex PfafID snw_alb' +! close (11, status = 'keep') + + print*, 'Ended tile loop for snow albedo. Biljana' + print*, 'There has been ',count_init_invalid,' inital non-valid snow values (out of',maxcat,')' + + END SUBROUTINE soil_snow_alb ! Biljana + + !-------------------------------------------------------------------------------------- + SUBROUTINE soil_para_hwsd (nx,ny,gfiler) ! Processing NGDC-HWSD-STATSGO merged soil properties with Woesten Soil @@ -5895,6 +6145,7 @@ SUBROUTINE open_landparam_nc4_files(N_tile) call DEF_VAR ( NCCatOUTID, CellID1,'TSB2' ,'water_transfer_param_4' ,'1' ) call DEF_VAR ( NCCatOUTID, CellID1,'WPWET' ,'wetness_at_wilting_point' ,'1' ) call DEF_VAR ( NCCatOUTID, CellID1,'DP2BR' ,'depth_to_bedrock' ,'mm' ) + call DEF_VAR ( NCCatOUTID, CellID1,'SNOWALB' ,'snow_albedo' ,'1' ) ! Biljana call DEF_VAR ( NCVegOUTID, CellID3,'ITY' ,'vegetation_type' ,'1' ) call DEF_VAR ( NCVegOUTID, CellID3,'Z2CH' ,'vegetation_height' ,'m' ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index fd381b50c..43f6dee9c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -45,6 +45,7 @@ module rmTinyCatchParaMod INTEGER, PARAMETER, public:: SRTM_maxcat = 291284 logical, public, save :: use_PEATMAP = .true. logical, public, save :: jpl_height = .true. + logical, public, save :: use_snow_albedo = .false. ! Biljana character*8, public, save :: LAIBCS = 'MODGEO' character*4, public, save :: SOILBCS = 'HWSD' character*6, public, save :: MODALB = 'MODIS2' @@ -117,6 +118,14 @@ SUBROUTINE init_bcs_config (LBSV) use_PEATMAP = .true. jpl_height = .true. + case ("NL6") !! Biljana New boundary condition option + LAIBCS = 'MODGEO' + SOILBCS = 'HWSD' + MODALB = 'MODIS2' + use_PEATMAP = .true. + jpl_height = .true. + use_snow_albedo= .true. + case ("DEV") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' From cb9ca6560cca12e3a303921e1fb8450f5879bbb1 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Mon, 8 Aug 2022 09:14:32 -0400 Subject: [PATCH 02/43] fix log message --- .../GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index f6f26e7b0..466c355d6 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -658,6 +658,7 @@ PROGRAM mkCatchParam ! --------------------------------------------------------------------- if(use_snow_albedo)then tmpstring = 'Step 14: Snow albedo from MODIS' + write (log_file,'(a)') trim(tmpstring) ! here I should test if snow_alb is written into nc file ! and call the soil_snow_alb subroutine only if it's not there ! inquire(file='clsm/catch_params.nc4', exist=file_exists) From 14f2eac0ad4367684cf3ef4eb27041d33bd35bfa Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Wed, 10 Aug 2022 20:46:01 -0400 Subject: [PATCH 03/43] edit comments, fix gnu error --- .../Utils/Raster/mkCatchParam.F90 | 27 ------------------- .../Utils/Raster/mod_process_hres_data.F90 | 19 +++++++------ .../Utils/Raster/rmTinyCatchParaMod.F90 | 4 +-- 3 files changed, 11 insertions(+), 39 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index 466c355d6..f8d22c2c7 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -653,43 +653,16 @@ PROGRAM mkCatchParam endif write (log_file,'(a)')' ' - ! Biljana ! Creating snow_alb_param.dat ! --------------------------------------------------------------------- if(use_snow_albedo)then tmpstring = 'Step 14: Snow albedo from MODIS' write (log_file,'(a)') trim(tmpstring) - ! here I should test if snow_alb is written into nc file - ! and call the soil_snow_alb subroutine only if it's not there - ! inquire(file='clsm/catch_params.nc4', exist=file_exists) - ! if (.not.file_exists) then - ! write (log_file,'(a)')' clsm/catch_params.nc4 is missing. Something is wrong. STOP!' - ! stop - ! else - ! inquire for snow albedo variable - - ! if not present, inititate write (log_file,'(a)')' Loading snow albedo...' call soil_snow_alb (nc,nr,gridnamer) write (log_file,'(a)')' Done.' - - ! if present and loaded, do nothing - - ! endif - - !fname_tmp = 'clsm/snow_alb_param.dat' - !write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' - !inquire(file=trim(fname_tmp), exist=file_exists) - !if (.not.file_exists) then - ! write (log_file,'(a)')' Creating file...' - ! call soil_snow_alb (nc,nr,gridnamer) - ! write (log_file,'(a)')' Done.' - !else - ! write (log_file,'(a,a)')' Using existing file.' - !endif write (log_file,'(a)')' ' endif ! if use_snow_albedo - ! Biljana -- done creating snow albedo file ! inquire(file='clsm/irrig.dat', exist=file_exists) ! if (.not.file_exists) call create_irrig_params (nc,nr,gridnamer) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index de11d794e..857607e07 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -26,7 +26,6 @@ MODULE process_hres_data use leap_year use MAPL_ConstantsMod use lsm_routines, ONLY: sibalb -USE IFPORT ! Biljana #if defined USE_EXTERNAL_FINDLOC use findloc_mod, only: findloc @@ -39,7 +38,7 @@ MODULE process_hres_data private public :: soil_para_hwsd,hres_lai,hres_gswp2, merge_lai_data, grid2tile_modis6 -public :: soil_snow_alb ! Biljana +public :: soil_snow_alb public :: modis_alb_on_tiles_high,modis_scale_para_high,hres_lai_no_gswp public :: histogram, create_mapping, esa2mosaic , esa2clm public :: grid2tile_ndep_t2m_alb, CREATE_ROUT_PARA_FILE, map_country_codes, get_country_codes @@ -2988,7 +2987,7 @@ END SUBROUTINE hres_gswp2 !---------------------------------------------------------------------- - SUBROUTINE soil_snow_alb (nx,ny,gfiler) ! Biljana + SUBROUTINE soil_snow_alb (nx,ny,gfiler) ! Implement snow albedo calculated from MOIDS 22-year climatology. Store snow albedo ! values in clsm/catch_params.nc4 @@ -3068,9 +3067,9 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) ! Biljana write(vv,'(i2.2)') vvtil write(hh,'(i2.2)') hhtil - fname = '/discover/nobackup/borescan/tools/idl/01_snow_fraction/06_modis_nsidc/'// & - '/data/data_out/snow_alb_all_08_Top99th_percentile_MOD10A1.A_30arcsec_'// & - '2000_2022_H'//hh//'V'//vv//'.nc' + fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/'// & + 'albedo/snow/MODIS/v1/snow_alb_MOD10A1.061_30arcsec_'// & + '2000_2022_H'//hh//'V'//vv//'.nc' ! Open the file. NF90_NOWRITE tells netCDF we want read-only access to the file. status=NF_OPEN(trim(fname),NF_NOWRITE, ncid) ; VERIFY_(STATUS) @@ -3096,7 +3095,7 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) ! Biljana ! open(11,file=trim(fname),form='formatted',status='unknown',action = 'write') ! loop over tiles - print*, 'Starting tile loop for snow albedo. Biljana' + print*, 'Starting tile loop for snow albedo. ' count_init_invalid=0 ! counter for non-valid snow albedo after matching tile size do n = 1, maxcat ! loop over tile @@ -3229,10 +3228,10 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) ! Biljana ! write (11,'(a)')'TileIndex PfafID snw_alb' ! close (11, status = 'keep') - print*, 'Ended tile loop for snow albedo. Biljana' + print*, 'Ended tile loop for snow albedo. ' print*, 'There has been ',count_init_invalid,' inital non-valid snow values (out of',maxcat,')' - END SUBROUTINE soil_snow_alb ! Biljana + END SUBROUTINE soil_snow_alb !-------------------------------------------------------------------------------------- @@ -6145,7 +6144,7 @@ SUBROUTINE open_landparam_nc4_files(N_tile) call DEF_VAR ( NCCatOUTID, CellID1,'TSB2' ,'water_transfer_param_4' ,'1' ) call DEF_VAR ( NCCatOUTID, CellID1,'WPWET' ,'wetness_at_wilting_point' ,'1' ) call DEF_VAR ( NCCatOUTID, CellID1,'DP2BR' ,'depth_to_bedrock' ,'mm' ) - call DEF_VAR ( NCCatOUTID, CellID1,'SNOWALB' ,'snow_albedo' ,'1' ) ! Biljana + call DEF_VAR ( NCCatOUTID, CellID1,'SNOWALB' ,'snow_albedo' ,'1' ) call DEF_VAR ( NCVegOUTID, CellID3,'ITY' ,'vegetation_type' ,'1' ) call DEF_VAR ( NCVegOUTID, CellID3,'Z2CH' ,'vegetation_height' ,'m' ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index 43f6dee9c..91ccd2459 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -45,7 +45,7 @@ module rmTinyCatchParaMod INTEGER, PARAMETER, public:: SRTM_maxcat = 291284 logical, public, save :: use_PEATMAP = .true. logical, public, save :: jpl_height = .true. - logical, public, save :: use_snow_albedo = .false. ! Biljana + logical, public, save :: use_snow_albedo = .false. character*8, public, save :: LAIBCS = 'MODGEO' character*4, public, save :: SOILBCS = 'HWSD' character*6, public, save :: MODALB = 'MODIS2' @@ -118,7 +118,7 @@ SUBROUTINE init_bcs_config (LBSV) use_PEATMAP = .true. jpl_height = .true. - case ("NL6") !! Biljana New boundary condition option + case ("NL6") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' From 0c092d05c212ee4cba5edc2916612718ddf78455 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Wed, 10 Aug 2022 20:50:56 -0400 Subject: [PATCH 04/43] edit --- .../Utils/Raster/mod_process_hres_data.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index 857607e07..577846f74 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -3068,8 +3068,8 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) write(hh,'(i2.2)') hhtil fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/'// & - 'albedo/snow/MODIS/v1/snow_alb_MOD10A1.061_30arcsec_'// & - '2000_2022_H'//hh//'V'//vv//'.nc' + '/albedo/snow/MODIS/v1/snow_alb_MOD10A1.061_30arcsec_'// & + '2000_2022_H'//hh//'V'//vv//'.nc' ! Open the file. NF90_NOWRITE tells netCDF we want read-only access to the file. status=NF_OPEN(trim(fname),NF_NOWRITE, ncid) ; VERIFY_(STATUS) From 3039868d6ba45c9ab13b4611896a8bed8f9b601b Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Wed, 10 Aug 2022 20:56:21 -0400 Subject: [PATCH 05/43] edit#2 --- .../Utils/Raster/mod_process_hres_data.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index 577846f74..0f31beed6 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -3067,9 +3067,7 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) write(vv,'(i2.2)') vvtil write(hh,'(i2.2)') hhtil - fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/'// & - '/albedo/snow/MODIS/v1/snow_alb_MOD10A1.061_30arcsec_'// & - '2000_2022_H'//hh//'V'//vv//'.nc' + fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v1/snow_alb_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' ! Open the file. NF90_NOWRITE tells netCDF we want read-only access to the file. status=NF_OPEN(trim(fname),NF_NOWRITE, ncid) ; VERIFY_(STATUS) From 473a6323642ca07b2230ca5ffae5c03a131f52e4 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Thu, 11 Aug 2022 13:05:19 -0400 Subject: [PATCH 06/43] add SNOWALB to CatmentRst object --- .../Utils/mk_restarts/CatchmentRst.F90 | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index 0ea043b54..b2ece52ff 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -76,6 +76,7 @@ module CatchmentRstMod real, allocatable :: sndzn1(:) real, allocatable :: sndzn2(:) real, allocatable :: sndzn3(:) + real, allocatable :: snowalb(:) real, allocatable :: ch(:,:) real, allocatable :: cm(:,:) real, allocatable :: cq(:,:) @@ -427,6 +428,9 @@ subroutine write_shared_nc4(this, formatter, rc) if (this%meta%has_variable('WW')) then call MAPL_VarWrite(formatter,"WW",this%ww) endif + if (this%meta%has_variable('SNOWALB')) then + call MAPL_VarWrite(formatter,"SNOWALB",this%snowalb) + endif _RETURN(_SUCCESS) @@ -483,6 +487,9 @@ subroutine allocate_catch(this,rc) if (this%meta%has_variable('TSURF')) then allocate( this% tsurf(ntiles) ) endif + if (this%meta%has_variable('SNOWALB')) then + allocate( this% snowalb(ntiles) ) + endif allocate( this% wesnn1(ntiles) ) allocate( this% wesnn2(ntiles) ) @@ -509,6 +516,7 @@ subroutine allocate_catch(this,rc) if (this%meta%has_variable('WW')) then allocate( this% ww(ntiles,4) ) endif + _RETURN(_SUCCESS) end subroutine allocate_catch @@ -529,6 +537,8 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) logical :: file_exists type(NetCDF4_Fileformatter) :: CatchFmt + type(Variable) :: var + type(FileMetadata) :: meta_ character*256 :: Iam = "add_bcs" @@ -611,6 +621,18 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) call MAPL_VarRead ( CatchFmt ,'WPWET', this%WPWET, __RC__) call MAPL_VarRead ( CatchFmt ,'DP2BR', DP2BR, __RC__) call MAPL_VarRead ( CatchFmt ,'POROS', this%POROS, __RC__) + meta_ = CatchFmt%read(__RC__) + if ( meta_%has_variable('SNOWALB')) then + if ( .not. allocated(this%snowalb)) allocate(this%snowalb(ntiles)) + call MAPL_VarRead ( CatchFmt ,'SNOWALB', this%snowalb, __RC__) + if ( .not. this%meta%has_variable('SNOWALB')) then + var = Variable(type=pFIO_REAL32, dimensions='tile') + call var%add_attribute('long_name', 'snow_albedo') + call var%add_attribute('units', '1') + call this%meta%add_variable('SNOWALB', var) + endif + endif + call CatchFmt%close() else open(unit=21, file=trim(DataDir)//'/clsm/mosaic_veg_typs_fracs',form='formatted') @@ -1069,6 +1091,10 @@ subroutine re_tile(this, InTileFile, OutBcsDir, OutTileFile, surflay, rc) var_out = this%tsurf(id_glb(:)) this%tsurf = var_out endif + if (this%meta%has_variable('SNOWALB')) then + var_out = this%snowalb(id_glb(:)) + this%snowalb = var_out + endif ! CH CM CQ FR WW ! WW From 1e1faa551d0bfe732f61c9e4d0018e0bf12f578f Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Thu, 18 Aug 2022 19:39:22 -0400 Subject: [PATCH 07/43] use snowalb info + add to surface rc file --- .../GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 49 ++++++++++++++++++- .../Shared/GEOS_SurfaceGridComp.rc | 8 +++ 2 files changed, 55 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index de12911f0..981c0eb33 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -896,6 +896,16 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddInternalSpec(GC ,& + LONG_NAME = 'snow_albedo' ,& + UNITS = '1' ,& + SHORT_NAME = 'SNOWALB' ,& + FRIENDLYTO = trim(COMP_NAME) ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) call MAPL_AddInternalSpec(GC ,& LONG_NAME = 'wetness_at_wilting_point' ,& UNITS = '1' ,& @@ -3763,6 +3773,7 @@ subroutine Driver ( RC ) real, dimension(:), pointer :: psis real, dimension(:), pointer :: bee real, dimension(:), pointer :: poros + real, dimension(:), pointer :: snowalb real, dimension(:), pointer :: wpwet real, dimension(:), pointer :: cond real, dimension(:), pointer :: gnu @@ -4132,6 +4143,10 @@ subroutine Driver ( RC ) integer :: nv, nVars integer :: nDims,dimSizes(3) integer :: ldas_ens_id, ldas_first_ens_id + + integer :: SNOW_ALBEDO_INFO + character(len=ESMF_MAXSTR) :: SURFRC + type(ESMF_Config) :: SCF !#--- ! -------------------------------------------------------------------------- @@ -4354,8 +4369,8 @@ subroutine Driver ( RC ) call MAPL_GetPointer(INTERNAL,CM ,'CM' ,RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL,CQ ,'CQ' ,RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL,FR ,'FR' ,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,DCQ ,'DCQ' ,RC=STATUS); VERIFY_(STATUS) - call MAPL_GetPointer(INTERNAL,DCH ,'DCH' ,RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL,DCQ ,'DCQ' ,RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL,DCH ,'DCH' ,RC=STATUS); VERIFY_(STATUS) if (N_CONST_LAND4SNWALB /= 0) then call MAPL_GetPointer(INTERNAL,RDU001 ,'RDU001' , RC=STATUS); VERIFY_(STATUS) call MAPL_GetPointer(INTERNAL,RDU002 ,'RDU002' , RC=STATUS); VERIFY_(STATUS) @@ -4857,6 +4872,25 @@ subroutine Driver ( RC ) SNOVR, SNONR, SNOVF, SNONF, & ! instantaneous snow albedos on tiles RCONSTIT, UUU, TPSN1OUT1, DRPAR, DFPAR) + + call MAPL_GetResource(MAPL,SURFRC,label='SURFRC:',default='GEOS_SurfaceGridComp.rc',RC=STATUS) ; VERIFY_(STATUS) + SCF = ESMF_ConfigCreate(rc=status) ; VERIFY_(STATUS) + call ESMF_ConfigLoadFile(SCF,SURFRC,rc=status) ; VERIFY_(STATUS) + call MAPL_GetResource(SCF,SNOW_ALBEDO_INFO,Label="SNOW_ALBEDO_INFO:", DEFAULT=0, RC=STATUS) ; VERIFY_(STATUS) + + if (SNOW_ALBEDO_INFO == 1) then + + call MAPL_GetPointer(INTERNAL,SNOWALB ,'SNOWALB' ,RC=STATUS); VERIFY_(STATUS) + + where (SNOWALB > 0. .and. SNOWALB <= 1.) + SNOVR = SNOWALB + SNONR = SNOWALB + SNOVF = SNOWALB + SNONF = SNOWALB + endwhere + + endif ! if SNOW_ALBEDO_INFO + ! -------------------------------------------------------------------------- ! albedo/swnet partitioning ! -------------------------------------------------------------------------- @@ -5535,6 +5569,17 @@ subroutine Driver ( RC ) SNOVR, SNONR, SNOVF, SNONF, & ! instantaneous snow albedos on tiles RCONSTIT, UUU, TPSN1OUT1,DRPAR, DFPAR) + if (SNOW_ALBEDO_INFO == 1) then + + where (SNOWALB > 0. .and. SNOWALB <= 1.) + SNOVR = SNOWALB + SNONR = SNOWALB + SNOVF = SNOWALB + SNONF = SNOWALB + endwhere + + endif ! if SNOW_ALBEDO_INFO + ALBVR = ALBVR *(1.-ASNOW) + SNOVR *ASNOW ALBVF = ALBVF *(1.-ASNOW) + SNOVF *ASNOW ALBNR = ALBNR *(1.-ASNOW) + SNONR *ASNOW diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc index 98900a487..ec7a5f67c 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc @@ -234,3 +234,11 @@ # # GEOSagcm=>PRESCRIBE_DVG: 0 # GEOSldas=>PRESCRIBE_DVG: 0 + +# ---- SNOW ALBEDO +# +# 0 : Snow albedo look-up table (default) +# 1 : Snow albedo derived from MODIS Collection MOD10A1.006 (2000-2022) +# +# GEOSagcm=>SNOW_ALBEDO_INFO: 0 +# GEOSldas=>SNOW_ALBEDO_INFO: 0 From 322d5bce5f80b244d664aa8c49cbc47729898d2c Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Wed, 24 Aug 2022 11:44:19 -0400 Subject: [PATCH 08/43] minor cleanup --- .../GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 | 2 -- .../Utils/Raster/rmTinyCatchParaMod.F90 | 6 ++++++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index f8d22c2c7..4a2164fac 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -653,8 +653,6 @@ PROGRAM mkCatchParam endif write (log_file,'(a)')' ' - ! Creating snow_alb_param.dat - ! --------------------------------------------------------------------- if(use_snow_albedo)then tmpstring = 'Step 14: Snow albedo from MODIS' write (log_file,'(a)') trim(tmpstring) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index 91ccd2459..d7c84ac2a 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -89,6 +89,7 @@ SUBROUTINE init_bcs_config (LBSV) GNU = 2.17 use_PEATMAP = .false. jpl_height = .false. + use_snow_albedo= .false. case ("GM4", "ICA") LAIBCS = 'GSWP2' @@ -96,6 +97,7 @@ SUBROUTINE init_bcs_config (LBSV) MODALB = 'MODIS2' use_PEATMAP = .false. jpl_height = .false. + use_snow_albedo= .false. case ("NL3") LAIBCS = 'MODGEO' @@ -103,6 +105,7 @@ SUBROUTINE init_bcs_config (LBSV) MODALB = 'MODIS2' use_PEATMAP = .false. jpl_height = .false. + use_snow_albedo= .false. case ("NL4") LAIBCS = 'MODGEO' @@ -110,6 +113,7 @@ SUBROUTINE init_bcs_config (LBSV) MODALB = 'MODIS2' use_PEATMAP = .false. jpl_height = .true. + use_snow_albedo= .false. case ("NL5") LAIBCS = 'MODGEO' @@ -117,6 +121,7 @@ SUBROUTINE init_bcs_config (LBSV) MODALB = 'MODIS2' use_PEATMAP = .true. jpl_height = .true. + use_snow_albedo= .false. case ("NL6") LAIBCS = 'MODGEO' @@ -132,6 +137,7 @@ SUBROUTINE init_bcs_config (LBSV) MODALB = 'MODIS2' use_PEATMAP = .true. jpl_height = .true. + use_snow_albedo= .false. end select From edf37eab56b1b1412f404d06d3ae62e6512f24d3 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Mon, 29 Aug 2022 07:58:18 -0400 Subject: [PATCH 09/43] add flag so only V06 has snowalb info in catch_params --- .../GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 30 +++--- .../Utils/Raster/make_bcs | 10 +- .../Utils/Raster/mkCatchParam.F90 | 4 +- .../Utils/Raster/mod_process_hres_data.F90 | 102 ++++++++---------- .../Utils/Raster/rmTinyCatchParaMod.F90 | 14 ++- .../Utils/mk_restarts/CatchmentRst.F90 | 15 +-- 6 files changed, 88 insertions(+), 87 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index 981c0eb33..2a4a44659 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -906,6 +906,7 @@ subroutine SetServices ( GC, RC ) RESTART = MAPL_RestartRequired ,& RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddInternalSpec(GC ,& LONG_NAME = 'wetness_at_wilting_point' ,& UNITS = '1' ,& @@ -4872,22 +4873,21 @@ subroutine Driver ( RC ) SNOVR, SNONR, SNOVF, SNONF, & ! instantaneous snow albedos on tiles RCONSTIT, UUU, TPSN1OUT1, DRPAR, DFPAR) - call MAPL_GetResource(MAPL,SURFRC,label='SURFRC:',default='GEOS_SurfaceGridComp.rc',RC=STATUS) ; VERIFY_(STATUS) SCF = ESMF_ConfigCreate(rc=status) ; VERIFY_(STATUS) call ESMF_ConfigLoadFile(SCF,SURFRC,rc=status) ; VERIFY_(STATUS) - call MAPL_GetResource(SCF,SNOW_ALBEDO_INFO,Label="SNOW_ALBEDO_INFO:", DEFAULT=0, RC=STATUS) ; VERIFY_(STATUS) + call MAPL_GetResource(SCF,SNOW_ALBEDO_INFO,Label="SNOW_ALBEDO_INFO:",DEFAULT=0,RC=STATUS) ; VERIFY_(STATUS) if (SNOW_ALBEDO_INFO == 1) then - call MAPL_GetPointer(INTERNAL,SNOWALB ,'SNOWALB' ,RC=STATUS); VERIFY_(STATUS) + call MAPL_GetPointer(INTERNAL,SNOWALB,'SNOWALB',RC=STATUS); VERIFY_(STATUS) - where (SNOWALB > 0. .and. SNOWALB <= 1.) - SNOVR = SNOWALB - SNONR = SNOWALB - SNOVF = SNOWALB - SNONF = SNOWALB - endwhere + where (SNOWALB > 0. .and. SNOWALB <= 1.) + SNOVR = SNOWALB + SNONR = SNOWALB + SNOVF = SNOWALB + SNONF = SNOWALB + endwhere endif ! if SNOW_ALBEDO_INFO @@ -5571,12 +5571,12 @@ subroutine Driver ( RC ) if (SNOW_ALBEDO_INFO == 1) then - where (SNOWALB > 0. .and. SNOWALB <= 1.) - SNOVR = SNOWALB - SNONR = SNOWALB - SNOVF = SNOWALB - SNONF = SNOWALB - endwhere + where (SNOWALB > 0. .and. SNOWALB <= 1.) + SNOVR = SNOWALB + SNONR = SNOWALB + SNOVF = SNOWALB + SNONF = SNOWALB + endwhere endif ! if SNOW_ALBEDO_INFO diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs index 2d7935fa9..fa04ecff6 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs @@ -105,8 +105,9 @@ echo " ${C2}ICA : Icarus (/discover/nobackup/ltakacs/bcs/Icarus/)" echo " ${C2}NL3 : Icarus-NLv3 (/discover/nobackup/ltakacs/bcs/Icarus-NLv3/)" echo " ${C2}NL4 : NLv4 [SMAP] (/discover/nobackup/projects/gmao/smap/bcs_NLv4/NLv4/)" echo " ${C2}NL5 : NLv5 [SMAP]" -echo " ${C2}NL6 : NLv6 [SMAP with snow albedo from MODIS]" -echo " ${C2}DEV : Development version${CR}" +echo " ${C2}V06 : V06 [NLv3 + JPL Veg Height + Peatlands + snow albedo from MODIS]" +echo " ${C2}V07 : V07 [NLv3 + Peatlands]" +echo " ${C2}V08 : V08 [NLv3 + snow albedo from MODIS]${CR}" echo " " if ( $HELPMODE != YES ) then echo " NOTE: Due to compiler differences, code improvements and bug fixes that" @@ -133,8 +134,9 @@ if ( $HELPMODE != YES ) then $dummy == 'NL3' | \ $dummy == 'NL4' | \ $dummy == 'NL5' | \ - $dummy == 'NL6' | \ - $dummy == 'DEV') then + $dummy == 'V06' | \ + $dummy == 'V07' | \ + $dummy == 'V08') then set lbcsv = $dummy else if ( $dummy == '' ) then echo $lbcsv diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index 4a2164fac..2a506c897 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -109,7 +109,7 @@ PROGRAM mkCatchParam USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " USAGE(6) =" -e: EASE : This is optional if catchment.def file is available already or " USAGE(7) =" the til file format is pre-Fortuna-2. " - USAGE(8) =" -v LBCSV : Choose bcs version (F25, GM4, ICA, NL3, NL4, NL5, or DEV) " + USAGE(8) =" -v LBCSV : Choose bcs version (F25, GM4, ICA, NL3, NL4, NL5, V06, V07 and V08) " ! Process Arguments !------------------ @@ -249,7 +249,7 @@ PROGRAM mkCatchParam close (10, status = 'keep') inquire(file='clsm/catch_params.nc4', exist=file_exists) - if (.not.file_exists) CALL open_landparam_nc4_files(N_tile) + if (.not.file_exists) CALL open_landparam_nc4_files(N_tile,use_snow_albedo) ! Creating cti_stats.dat ! ---------------------- diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index 0f31beed6..184209d61 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -2989,8 +2989,8 @@ END SUBROUTINE hres_gswp2 SUBROUTINE soil_snow_alb (nx,ny,gfiler) -! Implement snow albedo calculated from MOIDS 22-year climatology. Store snow albedo -! values in clsm/catch_params.nc4 +! Implement snow albedo calculated from MOIDS 22-year climatology. +! Store snow albedo values in clsm/catch_params.nc4 ! Biljana Orescanin July 2022, SSAI@NASA implicit none @@ -3000,7 +3000,7 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) character*200 :: fname character*2 :: vv,hh - integer :: n,maxcat,i,j,k,ncid,status + integer :: n,maxcat,j,ncid,status real,allocatable,dimension(:) :: min_lon,max_lon,min_lat,max_lat,snw_alb integer(kind=4),parameter :: xdim = 1200, ydim = 1200 real,parameter :: alb_res=10.0/1200.0 @@ -3009,28 +3009,28 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) real,dimension(xdim,ydim) :: stch_snw_alb_tmp real,dimension(36,18,xdim,ydim) :: stch_snw_alb real :: minlon,maxlon,minlat,maxlat,pad_lon,pad_lat - real :: sno_alb_cnt, sno_alb_sum,sno_alb_cnt2,sno_alb_sum2 + real :: sno_alb_cnt,sno_alb_sum,sno_alb_cnt2,sno_alb_sum2 integer :: vvtil_min,hhtil_min,vvtil_max,hhtil_max,hhtil,vvtil integer :: tindex1,pfaf1 - integer(kind=4) :: dummy,VarID,varid1,varid2,varid3 + integer(kind=4) :: dummy,varid1,varid2,varid3 integer(kind=4) :: imin,imax,jmin,jmax integer(kind=4) :: imin2,imax2,jmin2,jmax2,count_init_invalid - logical :: file_exists + logical :: file_exists ! Read number of catchment-tiles (maxcat) from catchment.def file fname='clsm/catchment.def' open (10,file=fname,status='old',action='read',form='formatted') read(10,*) maxcat - ! read min/max lat/lons, so those can be used to locate - ! snow albedo grids in the stitched MODIS albedo file + ! Read min/max lat/lons to use when locating snow albedo grids in + ! the stitched MODIS albedo file allocate (min_lon(1:maxcat)) allocate (min_lat(1:maxcat)) allocate (max_lon(1:maxcat)) allocate (max_lat(1:maxcat)) allocate (snw_alb(1:maxcat)) - ! before populating, set all snow albedo values to missing + ! Start by setting all snow albedo values to missing snw_alb(:)=-9999.0 do n = 1, maxcat @@ -3057,21 +3057,21 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) close (10,status='keep') !------------ Get the information on snow albedo ----- - ! ----------- The information on snow albedo is stored in 10x10 30-arcsec files. Read in this - ! information. Then loop over tiles to find a corresponding snow albedo mean + ! ----------- The information on snow albedo is stored in 10x10deg 30-arcsec resolution files. + ! ----------- Read in this information, then loop over the tiles to find a corresponding snow albedo. - ! read in all 10x10deg snow albedo files into a single [36,18,1200,1200] array - do hhtil=1,36 ! loop over all horziontal input files - do vvtil=1,18 ! loop over all vertical input files + ! Read in all 10x10deg snow albedo files into a single [36,18,1200,1200] array + do hhtil=1,36 ! loop over input files - horizontal direction + do vvtil=1,18 ! loop over input files - vertical direction write(vv,'(i2.2)') vvtil write(hh,'(i2.2)') hhtil fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v1/snow_alb_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' - ! Open the file. NF90_NOWRITE tells netCDF we want read-only access to the file. - status=NF_OPEN(trim(fname),NF_NOWRITE, ncid) ; VERIFY_(STATUS) - ! Get the varid of the data variable, based on its name. + ! Open the file. (NF90_NOWRITE ensures read-only access to the file) + status=NF_OPEN(trim(fname),NF_NOWRITE, ncid) ; VERIFY_(STATUS) + ! Based on vars name, get the varids. status=NF_INQ_VARID(ncid,'Snow_Albedo',VarID1) ; VERIFY_(STATUS) status=NF_INQ_VARID(ncid,'lon' ,VarID2) ; VERIFY_(STATUS) status=NF_INQ_VARID(ncid,'lat' ,VarID3) ; VERIFY_(STATUS) @@ -3082,42 +3082,36 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) ! Close the file, freeing all resources. status=NF_CLOSE(ncid); VERIFY_(STATUS) - ! store into large aray + ! Store snow albedo values into a single 4D aray stch_snw_alb(hhtil,vvtil,:,:)=stch_snw_alb_tmp enddo enddo - ! open the file to write snow albedo output in -! fname ='clsm/snow_alb_param.dat' -! open(11,file=trim(fname),form='formatted',status='unknown',action = 'write') - ! loop over tiles print*, 'Starting tile loop for snow albedo. ' - count_init_invalid=0 ! counter for non-valid snow albedo after matching tile size + count_init_invalid=0 ! counter for non-valid snow albedo avalues (informational use only; not needed for) - do n = 1, maxcat ! loop over tile + do n = 1, maxcat ! loop over tiles - ! set the current tile snow albedo to missing. Then start calculations to see if not missing. + ! Start by setting snow albedo to missing snw_alb(n)=-9999.0 - ! set sums and counts to zero - sno_alb_sum=0. - sno_alb_cnt=0. + ! Set sums and counts to zero + sno_alb_sum =0. + sno_alb_cnt =0. sno_alb_sum2=0. sno_alb_cnt2=0. - ! This tile has min/max lat/lon info. Use this info to identify which 10x10deg - ! snow albedo file(s) to read in (and then loop over these files). - ! Using ceiling and floor for max and min range, the "halo" approach is - ! implemented. + ! Use tile's min/max lat/lon info to identify the 10x10deg input file(s) + ! and read in snow albedo value(s). The "ceiling" and "floor" implements the "halo". vvtil_min= floor((min_lat(n)+ 90.0)/10.) hhtil_min= floor((min_lon(n)+180.0)/10.) vvtil_max=ceiling((max_lat(n)+ 90.0)/10.) hhtil_max=ceiling((max_lon(n)+180.0)/10.) - ! make sure vv's and hh's are within the range - ! if min>max, swap them + ! Safety checks: + ! 1. Make sure vv's and hh's are within the range. If min>max, swap them. if (vvtil_min .gt. vvtil_max) then dummy =vvtil_min vvtil_min=vvtil_max @@ -3129,27 +3123,27 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) hhtil_max=dummy endif - ! if beyond the range, bring them back + ! 2. Keep within the range. vvtil_min=max(vvtil_min,1) vvtil_max=min(vvtil_max,18) hhtil_min=max(hhtil_min,1) hhtil_max=min(hhtil_max,36) - do hhtil=hhtil_min,hhtil_max ! loop over all horziontal input files - do vvtil=vvtil_min,vvtil_max ! loop over all vertical input files + do hhtil=hhtil_min,hhtil_max ! loop over input files - horzontal direction + do vvtil=vvtil_min,vvtil_max ! loop over input files - vertical direction - ! find indices covered by the tile + ! Find indices ranges corresponding to the current tile area. imin=floor((min_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) imax=floor((max_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) jmin=floor((min_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) jmax=floor((max_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) - ! make sure to stay within the range + ! Keep within the range. imin=max(imin,1) imax=min(imax,xdim) jmin=max(jmin,1) jmax=min(jmax,ydim) - ! sum snow albedo values and counts for current tile corresponding indices + ! Generate sums and counts using current tile corresponding indices sno_alb_sum= sno_alb_sum + & sum(stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax), & stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & @@ -3161,14 +3155,15 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) end do ! vvtil end do ! hhtil - ! get mean snow albedo over the tile + + ! Calculate snow albedo for the current tile snw_alb(n) = sno_alb_sum / max(1.0,sno_alb_cnt) if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=-9999.0 !1.E15 - ! if no valid solition, and if tile size is smaller than the snow albedo grid box, - ! then expand the search areaby 1 tile in each direction + ! If no valid solition found, and if the tile size smaller than the snow albedo resolution, + ! expand the search area by 1-tile padding. - ! size of a tile (in both directions) + ! Size of a tile (in both directions) pad_lon=(max_lon(n)-min_lon(n)) pad_lat=(max_lat(n)-min_lat(n)) @@ -3176,10 +3171,10 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) count_init_invalid=count_init_invalid+1 - do hhtil=hhtil_min,hhtil_max ! loop over all horziontal input files - do vvtil=vvtil_min,vvtil_max ! loop over all vertical input files + do hhtil=hhtil_min,hhtil_max ! loop over input files - horzontal direction + do vvtil=vvtil_min,vvtil_max ! loop over input files - vertical direction - ! find indices of snow albedo array corresponding to the current tile + ! Repeat the steps for extracting snow albedo value imin2=floor((min_lon(n)-pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) imax2=floor((max_lon(n)+pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) jmin2=floor((min_lat(n)-pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) @@ -3189,7 +3184,6 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) jmin2=max(jmin2,1) jmax2=min(jmax2,ydim) - ! sum snow albedo values and counts for current tile corresponding indices sno_alb_sum2= sno_alb_sum2 + & sum(stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2), & stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & @@ -3207,9 +3201,6 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) endif - ! write the current tile value into the filr -! write (11,'(i10,i8,f13.4)') tindex1,pfaf1,snw_alb(n) - end do ! n-loop over tiles ! write snow albedo into clsm/catch_params.nc4 @@ -3221,11 +3212,6 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) STATUS = NF_CLOSE (NCID) ; VERIFY_(STATUS) endif -! ! close the output file -! write (11,'(a)')' ' -! write (11,'(a)')'TileIndex PfafID snw_alb' -! close (11, status = 'keep') - print*, 'Ended tile loop for snow albedo. ' print*, 'There has been ',count_init_invalid,' inital non-valid snow values (out of',maxcat,')' @@ -6093,12 +6079,13 @@ END SUBROUTINE gimms_clim_ndvi ! -------------------------------------------------------------------------- - SUBROUTINE open_landparam_nc4_files(N_tile) + SUBROUTINE open_landparam_nc4_files(N_tile,use_snow_albedo) implicit none integer :: NCCatOUTID, NCCatCNOUTID, NCVegOUTID integer :: STATUS, CellID1, CellID2, CellID3, SubID integer, intent (in) :: N_tile + logical, intent (in) :: use_snow_albedo integer, dimension(8) :: date_time_values character (22) :: time_stamp character (100) :: MYNAME @@ -6142,6 +6129,7 @@ SUBROUTINE open_landparam_nc4_files(N_tile) call DEF_VAR ( NCCatOUTID, CellID1,'TSB2' ,'water_transfer_param_4' ,'1' ) call DEF_VAR ( NCCatOUTID, CellID1,'WPWET' ,'wetness_at_wilting_point' ,'1' ) call DEF_VAR ( NCCatOUTID, CellID1,'DP2BR' ,'depth_to_bedrock' ,'mm' ) + if (use_snow_albedo) & call DEF_VAR ( NCCatOUTID, CellID1,'SNOWALB' ,'snow_albedo' ,'1' ) call DEF_VAR ( NCVegOUTID, CellID3,'ITY' ,'vegetation_type' ,'1' ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index d7c84ac2a..4ad1a859d 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -123,7 +123,7 @@ SUBROUTINE init_bcs_config (LBSV) jpl_height = .true. use_snow_albedo= .false. - case ("NL6") + case ("V06") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' @@ -131,13 +131,21 @@ SUBROUTINE init_bcs_config (LBSV) jpl_height = .true. use_snow_albedo= .true. - case ("DEV") + case ("V07") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' use_PEATMAP = .true. - jpl_height = .true. + jpl_height = .false. use_snow_albedo= .false. + + case ("V08") + LAIBCS = 'MODGEO' + SOILBCS = 'HWSD' + MODALB = 'MODIS2' + use_PEATMAP = .false. + jpl_height = .false. + use_snow_albedo= .true. end select diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index b2ece52ff..eca9d5e66 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -485,10 +485,10 @@ subroutine allocate_catch(this,rc) allocate( this% ghtcnt6(ntiles) ) if (this%meta%has_variable('TSURF')) then - allocate( this% tsurf(ntiles) ) + allocate( this% tsurf(ntiles) ) endif if (this%meta%has_variable('SNOWALB')) then - allocate( this% snowalb(ntiles) ) + allocate( this% snowalb(ntiles) ) endif allocate( this% wesnn1(ntiles) ) @@ -537,10 +537,10 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) logical :: file_exists type(NetCDF4_Fileformatter) :: CatchFmt - type(Variable) :: var - type(FileMetadata) :: meta_ + type(Variable) :: var + type(FileMetadata) :: meta_ - character*256 :: Iam = "add_bcs" + character*256 :: Iam = "add_bcs" ntiles = this%ntiles @@ -621,8 +621,10 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) call MAPL_VarRead ( CatchFmt ,'WPWET', this%WPWET, __RC__) call MAPL_VarRead ( CatchFmt ,'DP2BR', DP2BR, __RC__) call MAPL_VarRead ( CatchFmt ,'POROS', this%POROS, __RC__) + meta_ = CatchFmt%read(__RC__) - if ( meta_%has_variable('SNOWALB')) then + + if (meta_%has_variable('SNOWALB')) then if ( .not. allocated(this%snowalb)) allocate(this%snowalb(ntiles)) call MAPL_VarRead ( CatchFmt ,'SNOWALB', this%snowalb, __RC__) if ( .not. this%meta%has_variable('SNOWALB')) then @@ -1091,6 +1093,7 @@ subroutine re_tile(this, InTileFile, OutBcsDir, OutTileFile, surflay, rc) var_out = this%tsurf(id_glb(:)) this%tsurf = var_out endif + if (this%meta%has_variable('SNOWALB')) then var_out = this%snowalb(id_glb(:)) this%snowalb = var_out From cc188a93f16ee1f44064cfa90fc421593a2d1f7d Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Fri, 2 Sep 2022 09:25:34 -0400 Subject: [PATCH 10/43] version to lower case v --- .../GEOSsurface_GridComp/Utils/Raster/make_bcs | 12 ++++++------ .../Utils/Raster/mkCatchParam.F90 | 2 +- .../Utils/Raster/rmTinyCatchParaMod.F90 | 6 +++--- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs index fa04ecff6..b6b46c6fd 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs @@ -105,9 +105,9 @@ echo " ${C2}ICA : Icarus (/discover/nobackup/ltakacs/bcs/Icarus/)" echo " ${C2}NL3 : Icarus-NLv3 (/discover/nobackup/ltakacs/bcs/Icarus-NLv3/)" echo " ${C2}NL4 : NLv4 [SMAP] (/discover/nobackup/projects/gmao/smap/bcs_NLv4/NLv4/)" echo " ${C2}NL5 : NLv5 [SMAP]" -echo " ${C2}V06 : V06 [NLv3 + JPL Veg Height + Peatlands + snow albedo from MODIS]" -echo " ${C2}V07 : V07 [NLv3 + Peatlands]" -echo " ${C2}V08 : V08 [NLv3 + snow albedo from MODIS]${CR}" +echo " ${C2}v06 : v06 [NLv3 + JPL Veg Height + Peatlands + snow albedo from MODIS]" +echo " ${C2}v07 : v07 [NLv3 + Peatlands]" +echo " ${C2}v08 : v08 [NLv3 + snow albedo from MODIS]${CR}" echo " " if ( $HELPMODE != YES ) then echo " NOTE: Due to compiler differences, code improvements and bug fixes that" @@ -134,9 +134,9 @@ if ( $HELPMODE != YES ) then $dummy == 'NL3' | \ $dummy == 'NL4' | \ $dummy == 'NL5' | \ - $dummy == 'V06' | \ - $dummy == 'V07' | \ - $dummy == 'V08') then + $dummy == 'v06' | \ + $dummy == 'v07' | \ + $dummy == 'v08') then set lbcsv = $dummy else if ( $dummy == '' ) then echo $lbcsv diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index 2a506c897..0d4ef824f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -109,7 +109,7 @@ PROGRAM mkCatchParam USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " USAGE(6) =" -e: EASE : This is optional if catchment.def file is available already or " USAGE(7) =" the til file format is pre-Fortuna-2. " - USAGE(8) =" -v LBCSV : Choose bcs version (F25, GM4, ICA, NL3, NL4, NL5, V06, V07 and V08) " + USAGE(8) =" -v LBCSV : Choose bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07 and v08) " ! Process Arguments !------------------ diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index 4ad1a859d..9e37252f4 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -123,7 +123,7 @@ SUBROUTINE init_bcs_config (LBSV) jpl_height = .true. use_snow_albedo= .false. - case ("V06") + case ("v06") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' @@ -131,7 +131,7 @@ SUBROUTINE init_bcs_config (LBSV) jpl_height = .true. use_snow_albedo= .true. - case ("V07") + case ("v07") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' @@ -139,7 +139,7 @@ SUBROUTINE init_bcs_config (LBSV) jpl_height = .false. use_snow_albedo= .false. - case ("V08") + case ("v08") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' From d58f3460d2335aa2277a2c887136e0e9a46943f5 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Fri, 2 Sep 2022 17:43:01 -0400 Subject: [PATCH 11/43] =?UTF-8?q?cleanup=20and=20improved=20documentation?= =?UTF-8?q?=20-=20clarified=20that=20MODIS-based=20snow=20albedo=20option?= =?UTF-8?q?=20falls=20back=20to=20look-up=20table=20values=20where=20MODIS?= =?UTF-8?q?=20data=20are=20not=20available=20-=20clarified=20that=20MODIS-?= =?UTF-8?q?based=20snow=20albedo=20is=20not=20yet=20available=20for=20Catc?= =?UTF-8?q?hCN=20-=20moved=20processing=20of=20=E2=80=9CSNOW=5FALBEDO=5FIN?= =?UTF-8?q?FO:=E2=80=9D=20into=20SetServices()=20-=20cleaned=20up=20lists?= =?UTF-8?q?=20of=20available=20bcs=20versions=20-=20changed=20name=20of=20?= =?UTF-8?q?subroutine=20that=20processes=20MODIS=20snow=20albedo=20to=20mo?= =?UTF-8?q?re=20intuitive=20=E2=80=9CMODIS=5Fsnow=5Falb()=E2=80=9D=20-=20f?= =?UTF-8?q?ixed=20typos=20in=20comments?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../GEOS_CatchCNGridComp.F90 | 10 +++- .../GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 57 ++++++++++--------- .../Shared/GEOS_SurfaceGridComp.rc | 16 +++--- .../Utils/Raster/make_bcs | 12 ++-- .../Utils/Raster/mkCatchParam.F90 | 20 +++---- .../Utils/Raster/mod_process_hres_data.F90 | 22 +++---- 6 files changed, 75 insertions(+), 62 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 index b569ac07c..fab9000a2 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 @@ -50,7 +50,7 @@ subroutine SetServices ( GC, RC ) character(len=ESMF_MAXSTR) :: CATCHCN_VERSION character(len=ESMF_MAXSTR) :: SURFRC type(ESMF_Config) :: SCF - integer :: DO_GOSWIM, LSM_CHOICE, ATM_CO2 + integer :: DO_GOSWIM, LSM_CHOICE, ATM_CO2, SNOW_ALBEDO_INFO ! Begin... ! -------- @@ -69,6 +69,14 @@ subroutine SetServices ( GC, RC ) SCF = ESMF_ConfigCreate(rc=status) ; VERIFY_(STATUS) call ESMF_ConfigLoadFile(SCF,SURFRC,rc=status) ; VERIFY_(STATUS) call ESMF_ConfigGetAttribute (SCF, label='ATM_CO2:', value=ATM_CO2, DEFAULT=2, RC=STATUS) ; VERIFY_(STATUS) + + ! SNOW ALBEDO -- so far, only parameterization based on look-up table is implemented for CatchCN + ! 0 : parameterization based on look-up table + ! 1 : MODIS-derived snow albedo (where available, elsewhere fall back to option 0) + call ESMF_ConfigGetAttribute (SCF, label='SNOW_ALBEDO_INFO:', value=SNOW_ALBEDO_INFO, DEFAULT=0, RC=STATUS) ; VERIFY_(STATUS) + + _ASSERT( SNOW_ALBEDO_INFO==0, "SNOW_ALBEDO_INFO must be 0 for CatchCN") + call ESMF_ConfigGetAttribute (SCF, label='N_CONST_LAND4SNWALB:' , value=DO_GOSWIM , DEFAULT=0, RC=STATUS); VERIFY_(STATUS) if ( LSM_CHOICE == 2 ) then diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index 2a4a44659..0795918ae 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -136,7 +136,8 @@ module GEOS_CatchGridCompMod end type CATCH_WRAP !#-- -integer :: USE_ASCATZ0, Z0_FORMULATION, AEROSOL_DEPOSITION, N_CONST_LAND4SNWALB,CHOOSEMOSFC +integer :: USE_ASCATZ0, Z0_FORMULATION, AEROSOL_DEPOSITION, N_CONST_LAND4SNWALB +integer :: CHOOSEMOSFC, SNOW_ALBEDO_INFO real :: SURFLAY ! Default (Ganymed-3 and earlier) SURFLAY=20.0 for Old Soil Params ! (Ganymed-4 and later ) SURFLAY=50.0 for New Soil Params real :: FWETC, FWETL @@ -227,15 +228,21 @@ subroutine SetServices ( GC, RC ) call MAPL_GetResource (SCF, FWETL, label='FWETL:', DEFAULT=0.025, __RC__ ) endif - ! GOSWIM ANOW_ALBEDO + ! SNOW ALBEDO + ! 0 : parameterization based on look-up table + ! 1 : MODIS-derived snow albedo (where available, elsewhere fall back to option 0) + call MAPL_GetResource (SCF, SNOW_ALBEDO_INFO, label='SNOW_ALBEDO_INFO:', DEFAULT=0, __RC__ ) + + ! GOSWIM SNOW_ALBEDO ! 0 : GOSWIM snow albedo scheme is turned off ! 9 : i.e. N_CONSTIT in Stieglitz to turn on GOSWIM snow albedo scheme - call MAPL_GetResource (SCF, N_CONST_LAND4SNWALB, label='N_CONST_LAND4SNWALB:', DEFAULT=0 , __RC__ ) + call MAPL_GetResource (SCF, N_CONST_LAND4SNWALB, label='N_CONST_LAND4SNWALB:', DEFAULT=0, __RC__ ) ! 1: Use all GOCART aerosol values, 0: turn OFF everythying, ! 2: turn off dust ONLY,3: turn off Black Carbon ONLY,4: turn off Organic Carbon ONLY ! __________________________________________ - call MAPL_GetResource (SCF, AEROSOL_DEPOSITION, label='AEROSOL_DEPOSITION:', DEFAULT=0 , __RC__ ) + call MAPL_GetResource (SCF, AEROSOL_DEPOSITION, label='AEROSOL_DEPOSITION:', DEFAULT=0, __RC__ ) + call ESMF_ConfigDestroy(SCF, __RC__) ! Set the Run entry points @@ -4145,9 +4152,6 @@ subroutine Driver ( RC ) integer :: nDims,dimSizes(3) integer :: ldas_ens_id, ldas_first_ens_id - integer :: SNOW_ALBEDO_INFO - character(len=ESMF_MAXSTR) :: SURFRC - type(ESMF_Config) :: SCF !#--- ! -------------------------------------------------------------------------- @@ -4869,27 +4873,24 @@ subroutine Driver ( RC ) RHOFS, & SNWALB_VISMAX, SNWALB_NIRMAX, SLOPE, & WESNN, HTSNNN, SNDZN, & - ALBVR, ALBNR, ALBVF, ALBNF, & ! instantaneous snow-free albedos on tiles - SNOVR, SNONR, SNOVF, SNONF, & ! instantaneous snow albedos on tiles + ALBVR, ALBNR, ALBVF, ALBNF, & ! instantaneous snow-free albedos on tiles + SNOVR, SNONR, SNOVF, SNONF, & ! instantaneous snow albedos on tiles RCONSTIT, UUU, TPSN1OUT1, DRPAR, DFPAR) - call MAPL_GetResource(MAPL,SURFRC,label='SURFRC:',default='GEOS_SurfaceGridComp.rc',RC=STATUS) ; VERIFY_(STATUS) - SCF = ESMF_ConfigCreate(rc=status) ; VERIFY_(STATUS) - call ESMF_ConfigLoadFile(SCF,SURFRC,rc=status) ; VERIFY_(STATUS) - call MAPL_GetResource(SCF,SNOW_ALBEDO_INFO,Label="SNOW_ALBEDO_INFO:",DEFAULT=0,RC=STATUS) ; VERIFY_(STATUS) - if (SNOW_ALBEDO_INFO == 1) then - - call MAPL_GetPointer(INTERNAL,SNOWALB,'SNOWALB',RC=STATUS); VERIFY_(STATUS) - - where (SNOWALB > 0. .and. SNOWALB <= 1.) - SNOVR = SNOWALB - SNONR = SNOWALB - SNOVF = SNOWALB - SNONF = SNOWALB - endwhere - - endif ! if SNOW_ALBEDO_INFO + + ! where available, use MODIS-derived snow albedo from bcs (via Catch restart) + + call MAPL_GetPointer(INTERNAL,SNOWALB,'SNOWALB',RC=STATUS); VERIFY_(STATUS) + + where (SNOWALB > 0. .and. SNOWALB <= 1.) + SNOVR = SNOWALB + SNONR = SNOWALB + SNOVF = SNOWALB + SNONF = SNOWALB + endwhere + + endif ! -------------------------------------------------------------------------- ! albedo/swnet partitioning @@ -5570,7 +5571,9 @@ subroutine Driver ( RC ) RCONSTIT, UUU, TPSN1OUT1,DRPAR, DFPAR) if (SNOW_ALBEDO_INFO == 1) then - + + ! where available, use MODIS-derived snow albedo from bcs (via Catch restart) + where (SNOWALB > 0. .and. SNOWALB <= 1.) SNOVR = SNOWALB SNONR = SNOWALB @@ -5578,7 +5581,7 @@ subroutine Driver ( RC ) SNONF = SNOWALB endwhere - endif ! if SNOW_ALBEDO_INFO + endif ALBVR = ALBVR *(1.-ASNOW) + SNOVR *ASNOW ALBVF = ALBVF *(1.-ASNOW) + SNOVF *ASNOW diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc index ec7a5f67c..cf36809ae 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc @@ -114,6 +114,14 @@ # # GEOSagcm=>MODIS_DVG: 0 # NOT yet used in GEOSldas +# ---- SNOW ALBEDO +# +# 0 : Snow albedo parameterization based on look-up table (default) +# 1 : Snow albedo derived from MODIS Collection MOD10A1.006 (2000-2022) [NOT YET AVAILABLE FOR CatchCN] +# where available, elsewhere fall back to option 0 +# +# GEOSagcm=>SNOW_ALBEDO_INFO: 0 +# GEOSldas=>SNOW_ALBEDO_INFO: 0 #--------------------------------------------------------# # GOSWIM aerosol deposition on surface snow # @@ -235,10 +243,4 @@ # GEOSagcm=>PRESCRIBE_DVG: 0 # GEOSldas=>PRESCRIBE_DVG: 0 -# ---- SNOW ALBEDO -# -# 0 : Snow albedo look-up table (default) -# 1 : Snow albedo derived from MODIS Collection MOD10A1.006 (2000-2022) -# -# GEOSagcm=>SNOW_ALBEDO_INFO: 0 -# GEOSldas=>SNOW_ALBEDO_INFO: 0 + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs index b6b46c6fd..28078e395 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs @@ -103,12 +103,12 @@ echo " ${C2}F25 : Fortuna-2_5" echo " ${C2}GM4 : Ganymed-4_0 (/discover/nobackup/ltakacs/bcs/Ganymed-4_0/)" echo " ${C2}ICA : Icarus (/discover/nobackup/ltakacs/bcs/Icarus/)" echo " ${C2}NL3 : Icarus-NLv3 (/discover/nobackup/ltakacs/bcs/Icarus-NLv3/)" -echo " ${C2}NL4 : NLv4 [SMAP] (/discover/nobackup/projects/gmao/smap/bcs_NLv4/NLv4/)" -echo " ${C2}NL5 : NLv5 [SMAP]" -echo " ${C2}v06 : v06 [NLv3 + JPL Veg Height + Peatlands + snow albedo from MODIS]" -echo " ${C2}v07 : v07 [NLv3 + Peatlands]" -echo " ${C2}v08 : v08 [NLv3 + snow albedo from MODIS]${CR}" -echo " " +echo " ${C2}NL4 : NLv4 [SMAP] {NLv3 + JPL veg height} (/discover/nobackup/projects/gmao/smap/bcs_NLv4/NLv4/)" +echo " ${C2}NL5 : NLv5 [SMAP] {NLv3 + JPL veg height + PEATMAP}" +echo " ${C2}v06 : v06 {NLv3 + JPL veg height + PEATMAP + MODIS snow alb}" +echo " ${C2}v07 : v07 {NLv3 + PEATMAP}" +echo " ${C2}v08 : v08 {NLv3 + MODIS snow alb}${CR}" +echo " " if ( $HELPMODE != YES ) then echo " NOTE: Due to compiler differences, code improvements and bug fixes that" echo " have taken place since the above archived BCs were created, some parameter" diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index 0d4ef824f..163b3d80b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -9,7 +9,7 @@ PROGRAM mkCatchParam ! -y: Size of latitude dimension of input raster. DEFAULT: 4320 ! -b: position of the dateline in the first box. DEFAULT: DC ! -g: Gridname (name of the .til or .rst file without file extension) -! -v: LBCSV : Choose bcs version (ICA, NL3, NL4, NL5, or development) +! -v: LBCSV : Choose bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08) ! -e: EASE : This is optional if catchment.def file is available already or ! the til file format is pre-Fortuna-2. ! @@ -21,7 +21,7 @@ PROGRAM mkCatchParam ! Sarith Mahanama - March 23, 2012 ! Email: sarith.p.mahanama@nasa.gov - use rmTinyCatchParaMod + use rmTinyCatchParaMod ! includes: use_snow_albedo use process_hres_data ! use module_irrig_params, ONLY : create_irrig_params @@ -109,7 +109,7 @@ PROGRAM mkCatchParam USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " USAGE(6) =" -e: EASE : This is optional if catchment.def file is available already or " USAGE(7) =" the til file format is pre-Fortuna-2. " - USAGE(8) =" -v LBCSV : Choose bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07 and v08) " + USAGE(8) =" -v LBCSV : Choose bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08) " ! Process Arguments !------------------ @@ -654,13 +654,13 @@ PROGRAM mkCatchParam write (log_file,'(a)')' ' if(use_snow_albedo)then - tmpstring = 'Step 14: Snow albedo from MODIS' - write (log_file,'(a)') trim(tmpstring) - write (log_file,'(a)')' Loading snow albedo...' - call soil_snow_alb (nc,nr,gridnamer) - write (log_file,'(a)')' Done.' - write (log_file,'(a)')' ' - endif ! if use_snow_albedo + tmpstring = 'Step 14: Snow albedo from MODIS' + write (log_file,'(a)') trim(tmpstring) + write (log_file,'(a)')' Creating file...' + call MODIS_snow_alb (nc,nr,gridnamer) + write (log_file,'(a)')' Done.' + write (log_file,'(a)')' ' + endif ! inquire(file='clsm/irrig.dat', exist=file_exists) ! if (.not.file_exists) call create_irrig_params (nc,nr,gridnamer) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index 184209d61..7a64b69cd 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -38,7 +38,7 @@ MODULE process_hres_data private public :: soil_para_hwsd,hres_lai,hres_gswp2, merge_lai_data, grid2tile_modis6 -public :: soil_snow_alb +public :: MODIS_snow_alb public :: modis_alb_on_tiles_high,modis_scale_para_high,hres_lai_no_gswp public :: histogram, create_mapping, esa2mosaic , esa2clm public :: grid2tile_ndep_t2m_alb, CREATE_ROUT_PARA_FILE, map_country_codes, get_country_codes @@ -2987,9 +2987,9 @@ END SUBROUTINE hres_gswp2 !---------------------------------------------------------------------- - SUBROUTINE soil_snow_alb (nx,ny,gfiler) + SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) -! Implement snow albedo calculated from MOIDS 22-year climatology. +! Implement snow albedo calculated from MODIS 22-year climatology. ! Store snow albedo values in clsm/catch_params.nc4 ! Biljana Orescanin July 2022, SSAI@NASA @@ -3129,8 +3129,8 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) hhtil_min=max(hhtil_min,1) hhtil_max=min(hhtil_max,36) - do hhtil=hhtil_min,hhtil_max ! loop over input files - horzontal direction - do vvtil=vvtil_min,vvtil_max ! loop over input files - vertical direction + do hhtil=hhtil_min,hhtil_max ! loop through input files - horizontal direction + do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction ! Find indices ranges corresponding to the current tile area. imin=floor((min_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) @@ -3160,8 +3160,8 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) snw_alb(n) = sno_alb_sum / max(1.0,sno_alb_cnt) if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=-9999.0 !1.E15 - ! If no valid solition found, and if the tile size smaller than the snow albedo resolution, - ! expand the search area by 1-tile padding. + ! If no valid solution found, and if tile size smaller than snow albedo resolution, + ! expand search area by 1-tile padding. ! Size of a tile (in both directions) pad_lon=(max_lon(n)-min_lon(n)) @@ -3171,8 +3171,8 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) count_init_invalid=count_init_invalid+1 - do hhtil=hhtil_min,hhtil_max ! loop over input files - horzontal direction - do vvtil=vvtil_min,vvtil_max ! loop over input files - vertical direction + do hhtil=hhtil_min,hhtil_max ! loop through input files - horizontal direction + do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction ! Repeat the steps for extracting snow albedo value imin2=floor((min_lon(n)-pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) @@ -3213,9 +3213,9 @@ SUBROUTINE soil_snow_alb (nx,ny,gfiler) endif print*, 'Ended tile loop for snow albedo. ' - print*, 'There has been ',count_init_invalid,' inital non-valid snow values (out of',maxcat,')' + print*, 'Initially found ', count_init_invalid, ' tiles with no-data values for snow albedo (out of ', maxcat,')' - END SUBROUTINE soil_snow_alb + END SUBROUTINE MODIS_snow_alb !-------------------------------------------------------------------------------------- From f9a05d8725b8360b9b9cba510cab0e38a1047d2c Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 7 Sep 2022 18:30:11 -0400 Subject: [PATCH 12/43] Additional documentation and cleanup of snow albedo changes: - Introduced string (SNOWALB) for snow albedo version. - Edited comments. --- .../Utils/Raster/mkCatchParam.F90 | 42 +++++---- .../Utils/Raster/mkSMAPTilesPara_v2.F90 | 6 +- .../Utils/Raster/mod_process_hres_data.F90 | 15 ++-- .../Utils/Raster/rmTinyCatchParaMod.F90 | 88 ++++++++++++------- 4 files changed, 86 insertions(+), 65 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index afbf98d4a..92a744dbb 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -4,14 +4,14 @@ PROGRAM mkCatchParam ! ! !ARGUMENTS: ! -! Usage = "mkCatchParam -x nx -y ny -g Gridname -b DL -v LBSV -e EASE" +! Usage = "mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV -e EASE" ! -x: Size of longitude dimension of input raster. DEFAULT: 8640 ! -y: Size of latitude dimension of input raster. DEFAULT: 4320 ! -b: position of the dateline in the first box. DEFAULT: DC ! -g: Gridname (name of the .til or .rst file without file extension) -! -v: LBCSV : Choose bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08) -! -e: EASE : This is optional if catchment.def file is available already or -! the til file format is pre-Fortuna-2. +! -v: LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08) +! -e: EASE : This is optional if catchment.def file is available already or +! the til file format is pre-Fortuna-2. ! ! ! This program is good to generate @@ -21,7 +21,7 @@ PROGRAM mkCatchParam ! Sarith Mahanama - March 23, 2012 ! Email: sarith.p.mahanama@nasa.gov - use rmTinyCatchParaMod ! includes: use_snow_albedo + use rmTinyCatchParaMod use process_hres_data ! use module_irrig_params, ONLY : create_irrig_params @@ -33,7 +33,7 @@ PROGRAM mkCatchParam ! NC and NR are typically overwritten through command-line arguments "-x nx -y ny". integer :: NC = i_raster, NR = j_raster - character*4 :: LBSV = 'DEF' + character*5 :: LBCSV = 'UNDEF' character*128 :: GridName = '' character*128 :: ARG, MaskFile character*256 :: CMD @@ -60,6 +60,7 @@ PROGRAM mkCatchParam character*200 :: tmpstring, tmpstring1, tmpstring2 character*200 :: fname_tmp, fname_tmp2, fname_tmp3, fname_tmp4 integer :: N_tile + logical :: process_snow_albedo = .false. ! --------- VARIABLES FOR *OPENMP* PARALLEL ENVIRONMENT ------------ ! @@ -102,14 +103,14 @@ PROGRAM mkCatchParam ! call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH') ! call execute_command_line('cd ..') - USAGE(1) ="Usage: mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV -e EASE " - USAGE(2) =" -x: Size of longitude dimension of input raster. DEFAULT: 8640 " - USAGE(3) =" -y: Size of latitude dimension of input raster. DEFAULT: 4320 " - USAGE(4) =" -g: Gridname (name of the .til or .rst file without file extension) " - USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " - USAGE(6) =" -e: EASE : This is optional if catchment.def file is available already or " - USAGE(7) =" the til file format is pre-Fortuna-2. " - USAGE(8) =" -v LBCSV : Choose bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08) " + USAGE(1) ="Usage: mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV -e EASE " + USAGE(2) =" -x: Size of longitude dimension of input raster. DEFAULT: 8640 " + USAGE(3) =" -y: Size of latitude dimension of input raster. DEFAULT: 4320 " + USAGE(4) =" -g: Gridname (name of the .til or .rst file without file extension) " + USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " + USAGE(6) =" -e: EASE : This is optional if catchment.def file is available already or " + USAGE(7) =" the til file format is pre-Fortuna-2. " + USAGE(8) =" -v LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08) " ! Process Arguments !------------------ @@ -153,9 +154,9 @@ PROGRAM mkCatchParam case ('g') GridName = trim(arg) case ('v') - LBSV = trim(arg) + LBCSV = trim(arg) if (trim(arg).eq."F25") F25Tag = .true. - call init_bcs_config (trim(LBSV)) + call init_bcs_config (trim(LBCSV)) ! get bcs details from version string case ('b') DL = trim(arg) case ('e') @@ -191,11 +192,14 @@ PROGRAM mkCatchParam if(use_PEATMAP) PEATSOURCE = 'PEATMAP' if(jpl_height) VEGZSOURCE = 'JPL' + if (trim(SNOWALB)=='MODC006')) process_snow_albedo=.true. + if(n_threads == 1) then write (log_file,'(a)')trim(LAIBCS) write (log_file,'(a)')trim(MODALB) write (log_file,'(a)')trim(SOILBCS) + write (log_file,'(a)')trim(SNOWALB) write (log_file,'(a)')trim(MaskFile) write (log_file,'(a)')trim(PEATSOURCE) write (log_file,'(a)')trim(VEGZSOURCE) @@ -249,7 +253,7 @@ PROGRAM mkCatchParam close (10, status = 'keep') inquire(file='clsm/catch_params.nc4', exist=file_exists) - if (.not.file_exists) CALL open_landparam_nc4_files(N_tile,use_snow_albedo) + if (.not.file_exists) CALL open_landparam_nc4_files(N_tile,process_snow_albedo) ! Creating cti_stats.dat ! ---------------------- @@ -653,8 +657,8 @@ PROGRAM mkCatchParam endif write (log_file,'(a)')' ' - if(use_snow_albedo)then - tmpstring = 'Step 14: Snow albedo from MODIS' + if(process_snow_albedo)then + tmpstring = 'Step 14: Static snow albedo from MODIS' write (log_file,'(a)') trim(tmpstring) write (log_file,'(a)')' Creating file...' call MODIS_snow_alb (nc,nr,gridnamer) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90 index 70025ec11..55c1abc37 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90 @@ -16,7 +16,7 @@ PROGRAM mkSMAPTilesPara_v2 use LogRectRasterizeMod implicit none - character*4 :: LBSV = 'DEF' + character*5 :: LBCSV = 'UNDEF' integer i,j,ig,jg,i0,iop,n,d1,d2,j1,j2,i1,i2,ix, jx,icount,pcount integer :: NC = i_raster, NR = j_raster, NT = 16330000, ND = 10000, ND_raster = 10000 @@ -95,7 +95,7 @@ PROGRAM mkSMAPTilesPara_v2 elseif ( trim(arg) == '-v' ) then i = i+1 - call get_command_argument(i,LBSV) + call get_command_argument(i,LBCSV) else ! stop for any other arguments @@ -753,7 +753,7 @@ PROGRAM mkSMAPTilesPara_v2 ! now run mkCatchParam ! -------------------- - tmpstring1 = '-e EASE -g '//trim(gfile)//' -v '//trim(LBSV) + tmpstring1 = '-e EASE -g '//trim(gfile)//' -v '//trim(LBCSV) write(tmpstring2,'(2(a2,x,i5,x))')'-x',nc,'-y',nr tmpstring = 'bin/mkCatchParam.x '//trim(tmpstring2)//' '//trim(tmpstring1) print *,trim(tmpstring) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index 8128b7586..92b3b161f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -2994,9 +2994,9 @@ END SUBROUTINE hres_gswp2 SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) -! Implement snow albedo calculated from MODIS 22-year climatology. -! Store snow albedo values in clsm/catch_params.nc4 -! Biljana Orescanin July 2022, SSAI@NASA + ! Process static snow albedo calculated from MODIS climatology and write into clsm/catch_params.nc4 + ! + ! Biljana Orescanin July 2022, SSAI@NASA implicit none integer, intent (in) :: nx, ny @@ -3099,9 +3099,6 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) do n = 1, maxcat ! loop over tiles - ! Start by setting snow albedo to missing - snw_alb(n)=-9999.0 - ! Set sums and counts to zero sno_alb_sum =0. sno_alb_cnt =0. @@ -6089,13 +6086,13 @@ END SUBROUTINE gimms_clim_ndvi ! -------------------------------------------------------------------------- - SUBROUTINE open_landparam_nc4_files(N_tile,use_snow_albedo) + SUBROUTINE open_landparam_nc4_files(N_tile,process_snow_albedo) implicit none integer :: NCCatOUTID, NCCatCNOUTID, NCVegOUTID integer :: STATUS, CellID1, CellID2, CellID3, SubID integer, intent (in) :: N_tile - logical, intent (in) :: use_snow_albedo + logical, intent (in) :: process_snow_albedo integer, dimension(8) :: date_time_values character (22) :: time_stamp character (100) :: MYNAME @@ -6139,7 +6136,7 @@ SUBROUTINE open_landparam_nc4_files(N_tile,use_snow_albedo) call DEF_VAR ( NCCatOUTID, CellID1,'TSB2' ,'water_transfer_param_4' ,'1' ) call DEF_VAR ( NCCatOUTID, CellID1,'WPWET' ,'wetness_at_wilting_point' ,'1' ) call DEF_VAR ( NCCatOUTID, CellID1,'DP2BR' ,'depth_to_bedrock' ,'mm' ) - if (use_snow_albedo) & + if (process_snow_albedo) & call DEF_VAR ( NCCatOUTID, CellID1,'SNOWALB' ,'snow_albedo' ,'1' ) call DEF_VAR ( NCVegOUTID, CellID3,'ITY' ,'vegetation_type' ,'1' ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index fb75e1966..2060ca57b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -42,14 +42,19 @@ module rmTinyCatchParaMod public tgen, sat_param,REFORMAT_VEGFILES,base_param,ts_param public :: Get_MidTime, Time_Interp_Fac, compute_stats, c_data public :: ascat_r0, jpl_canoph, NC_VarID, init_bcs_config + INTEGER, PARAMETER, public:: SRTM_maxcat = 291284 - logical, public, save :: use_PEATMAP = .true. - logical, public, save :: jpl_height = .true. - logical, public, save :: use_snow_albedo = .false. - character*8, public, save :: LAIBCS = 'MODGEO' - character*4, public, save :: SOILBCS = 'HWSD' - character*6, public, save :: MODALB = 'MODIS2' - REAL, public, save :: GNU = 1.0 + + ! The following variables define the details of the BCS version (data sources). + ! Initialize to dummy values here and set to desired values in init_bcs_config(). + + logical, public, save :: use_PEATMAP = .false. + logical, public, save :: jpl_height = .false. + character*8, public, save :: LAIBCS = 'UNDEF' + character*4, public, save :: SOILBCS = 'UNDEF' + character*6, public, save :: MODALB = 'UNDEF' + character*8, public, save :: SNOWALB = 'UNDEF' + REAL, public, save :: GNU = MAPL_UNDEF type :: mineral_perc real :: clay_perc @@ -59,93 +64,108 @@ module rmTinyCatchParaMod contains - SUBROUTINE init_bcs_config (LBSV) + SUBROUTINE init_bcs_config (LBCSV) + ! determine BCs details from land BCs version string (LBCSV) + ! + ! LAIBCS: Leaf-Area-Index data set. DEFAULT : MODGEO + ! GLASSA : 8-day AVHRR clim, 1981-2017, 7200x3600 grid + ! GLASSM : 8-day MODIS clim, 2000-2017, 7200x3600 grid + ! MODISV6 : 8-day clim, 2002.01-2016.10, 86400x43200 grid + ! MODGEO : MODIS with GEOLAND2 overlaid on South America, Africa, and Australia + ! GEOLAND2 : 10-day clim, 1999-2011, 40320x20160 grid + ! GSWP2 : Monthly clim, 1982-1998, 360x180 grid + ! MODIS : 8-day clim, 2000-2013, 43200x21600 grid + ! GSWPH : Monthly clim, 1982-1998, 43200x21600 grid + ! + ! MODALB: MODIS Albedo data (snow-free). DEFAULT : MODIS2 + ! MODIS1 : 16-day clim, 1'x1' (21600x10800) MODIS data, 2000-2004 + ! MODIS2 : 8-day clim, 30"x30"(43200x21600) MODIS data, 2001-2011 + ! + ! SNOWALB: Snow albedo data. DEFAULT : LUT + ! LUT : Parameterization based on look-up table values. + ! MODC006 : Static snow albedo derived from MODIS Collection 6 data where available, LUT elsewhere. + ! + ! SOILBCS: Soil parameter data. DEFAULT : HWSD + ! HWSD : Merged HWSD-STATSGO2 soil properties on 43200x21600 with Woesten et al. (1999) parameters + implicit none + + character(*), intent (in) :: LBCSV ! land BCs version - character(*), intent (in) :: LBSV ! LBSV = land BCs version (?) - -! LAIBCS: Choice of LAI data set. DEFAULT : MODGEO -! GLASSA : 8-day AVHRR climatology from the period 1981-2017 on 7200x3600 grid -! GLASSM : 8-day MODIS climatology from the period 2000-2017 on 7200x3600 grid -! MODISV6 : 8-day climatology from the period 2002.01-2016.10 on 86400x43200 grid -! MODGEO : MODIS with GEOLAND2 overlaid on South America, Afirca and Australia -! GEOLAND2: 10-day climatology from the period 1999-2011 on 40320x20160 grid -! GSWP2 : Monthly climatology from the period 1982-1998 on 360x180 grid -! MODIS : 8-day climatology from the period 2000-2013 on 43200x21600 grid -! GSWPH : Monthly climatology from the period 1982-1998 on 43200x21600 grid " -! MODALB: Choice of MODIS Albedo data. DEFAULT : MODIS2 -! MODIS1 : 16-day Climatology from 1'x1 (21600x10800) MODIS data from the period 2000-2004 -! MODIS2 : 8-day Climatology from 30"x30"(43200x21600) MODIS data from the period 2001-2011 -! SOILBCS:Choice of soil data. DEFAULT :HWSD -! HWSD : Merged HWSD-STATSGO2 soil properties on 43200x21600 with Woesten et al. (1999) Parameters - - select case (trim(LBSV)) + select case (trim(LBCSV)) case ("F25") LAIBCS = 'GSWP2' SOILBCS = 'NGDC' MODALB = 'MODIS1' + SNOWALB = 'LUT' GNU = 2.17 use_PEATMAP = .false. jpl_height = .false. - use_snow_albedo= .false. case ("GM4", "ICA") LAIBCS = 'GSWP2' SOILBCS = 'NGDC' MODALB = 'MODIS2' + SNOWALB = 'LUT' + GNU = 1.0 use_PEATMAP = .false. jpl_height = .false. - use_snow_albedo= .false. case ("NL3") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' + SNOWALB = 'LUT' + GNU = 1.0 use_PEATMAP = .false. jpl_height = .false. - use_snow_albedo= .false. case ("NL4") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' + SNOWALB = 'LUT' + GNU = 1.0 use_PEATMAP = .false. jpl_height = .true. - use_snow_albedo= .false. case ("NL5") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' + SNOWALB = 'LUT' + GNU = 1.0 use_PEATMAP = .true. jpl_height = .true. - use_snow_albedo= .false. case ("v06") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' + SNOWALB = 'MODC006' + GNU = 1.0 use_PEATMAP = .true. jpl_height = .true. - use_snow_albedo= .true. case ("v07") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' + SNOWALB = 'LUT' + GNU = 1.0 use_PEATMAP = .true. jpl_height = .false. - use_snow_albedo= .false. case ("v08") LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' + SNOWALB = 'MODC006' + GNU = 1.0 use_PEATMAP = .false. jpl_height = .false. - use_snow_albedo= .true. end select From 85c37079bafc02889b761e5976879fd28aa5db2a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Wed, 7 Sep 2022 22:27:17 -0400 Subject: [PATCH 13/43] fixed build error in last commit --- .../GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index 2060ca57b..e98bc7f61 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -10,6 +10,7 @@ module rmTinyCatchParaMod use date_time_util use leap_year use MAPL_ConstantsMod + use MAPL_Base, ONLY: MAPL_UNDEF use lsm_routines, ONLY: sibalb implicit none From 034aac09493c2948e681aefa3ed69f2b312fa8e9 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 8 Sep 2022 08:48:34 -0400 Subject: [PATCH 14/43] fixed another build error in previous commit --- .../GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index 92a744dbb..47878299e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -192,7 +192,7 @@ PROGRAM mkCatchParam if(use_PEATMAP) PEATSOURCE = 'PEATMAP' if(jpl_height) VEGZSOURCE = 'JPL' - if (trim(SNOWALB)=='MODC006')) process_snow_albedo=.true. + if (trim(SNOWALB)=='MODC006') process_snow_albedo=.true. if(n_threads == 1) then From 0718bc8db4ee23f0187fb7669b909f117d19036f Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Fri, 16 Sep 2022 11:18:53 -0400 Subject: [PATCH 15/43] fix case options --- .../GEOSsurface_GridComp/Utils/Raster/make_bcs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs index f4fe05525..ce59fa564 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs @@ -141,9 +141,9 @@ if ( $HELPMODE != YES ) then $dummy == 'NL3' | \ $dummy == 'NL4' | \ $dummy == 'NL5' | \ - $dummy == 'v06' | \ - $dummy == 'v07' | \ - $dummy == 'v08') then + $dummy == 'V06' | \ + $dummy == 'V07' | \ + $dummy == 'V08') then set lbcsv = $dummy else if ( $dummy == '' ) then echo $lbcsv From c7fb5a0f18e7d42f67f3e1f2bf7633c3b1810bb1 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Fri, 16 Sep 2022 11:21:21 -0400 Subject: [PATCH 16/43] fix MODIS collection info --- .../GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc | 2 +- .../GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 | 2 +- .../Utils/Raster/rmTinyCatchParaMod.F90 | 6 +++--- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc index cf36809ae..3491c4f5c 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc @@ -117,7 +117,7 @@ # ---- SNOW ALBEDO # # 0 : Snow albedo parameterization based on look-up table (default) -# 1 : Snow albedo derived from MODIS Collection MOD10A1.006 (2000-2022) [NOT YET AVAILABLE FOR CatchCN] +# 1 : Snow albedo derived from MODIS Collection MOD10A1.061 (Feb/2000 - Mar/2022) [NOT YET AVAILABLE FOR CatchCN] # where available, elsewhere fall back to option 0 # # GEOSagcm=>SNOW_ALBEDO_INFO: 0 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index 47878299e..b07b9ba3b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -192,7 +192,7 @@ PROGRAM mkCatchParam if(use_PEATMAP) PEATSOURCE = 'PEATMAP' if(jpl_height) VEGZSOURCE = 'JPL' - if (trim(SNOWALB)=='MODC006') process_snow_albedo=.true. + if (trim(SNOWALB)=='MODC061') process_snow_albedo=.true. if(n_threads == 1) then diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index e98bc7f61..992bb7310 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -85,7 +85,7 @@ SUBROUTINE init_bcs_config (LBCSV) ! ! SNOWALB: Snow albedo data. DEFAULT : LUT ! LUT : Parameterization based on look-up table values. - ! MODC006 : Static snow albedo derived from MODIS Collection 6 data where available, LUT elsewhere. + ! MODC061 : Static snow albedo derived from MODIS Collection 6.1 data where available, LUT elsewhere. ! ! SOILBCS: Soil parameter data. DEFAULT : HWSD ! HWSD : Merged HWSD-STATSGO2 soil properties on 43200x21600 with Woesten et al. (1999) parameters @@ -145,7 +145,7 @@ SUBROUTINE init_bcs_config (LBCSV) LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' - SNOWALB = 'MODC006' + SNOWALB = 'MODC061' GNU = 1.0 use_PEATMAP = .true. jpl_height = .true. @@ -163,7 +163,7 @@ SUBROUTINE init_bcs_config (LBCSV) LAIBCS = 'MODGEO' SOILBCS = 'HWSD' MODALB = 'MODIS2' - SNOWALB = 'MODC006' + SNOWALB = 'MODC061' GNU = 1.0 use_PEATMAP = .false. jpl_height = .false. From dd2f91a203a19c448e21922c61de8dba469a1d84 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Fri, 16 Sep 2022 12:11:49 -0400 Subject: [PATCH 17/43] better option for cases --- .../GEOSsurface_GridComp/Utils/Raster/make_bcs | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs index ce59fa564..2253db27f 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs @@ -133,7 +133,6 @@ if ( $HELPMODE != YES ) then set dummy = `echo $<` - set dummy = `echo $dummy | tr "[:lower:]" "[:upper:]"` if( $dummy == 'F25' | \ $dummy == 'GM4' | \ @@ -141,9 +140,9 @@ if ( $HELPMODE != YES ) then $dummy == 'NL3' | \ $dummy == 'NL4' | \ $dummy == 'NL5' | \ - $dummy == 'V06' | \ - $dummy == 'V07' | \ - $dummy == 'V08') then + $dummy == 'v06' | \ + $dummy == 'v07' | \ + $dummy == 'v08') then set lbcsv = $dummy else if ( $dummy == '' ) then echo $lbcsv From cff2c973f0665c5a411bd63c6ccfbaed3b1f86e4 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Fri, 16 Sep 2022 14:12:56 -0400 Subject: [PATCH 18/43] snowalb optional --- .../GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index 0795918ae..f02c76f6f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -903,17 +903,6 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'snow_albedo' ,& - UNITS = '1' ,& - SHORT_NAME = 'SNOWALB' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) - call MAPL_AddInternalSpec(GC ,& LONG_NAME = 'wetness_at_wilting_point' ,& UNITS = '1' ,& @@ -1464,6 +1453,17 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + call MAPL_AddInternalSpec(GC ,& + LONG_NAME = 'snow_albedo' ,& + UNITS = '1' ,& + SHORT_NAME = 'SNOWALB' ,& + FRIENDLYTO = trim(COMP_NAME) ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RESTART = MAPL_RestartOptional ,& + RC=STATUS ) + VERIFY_(STATUS) + !---------- GOSWIM snow impurity related variables ---------- if (N_CONST_LAND4SNWALB /= 0) then From 5f7aba2b341e14b25e3657d1cff7b9ab36fbad2d Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Fri, 16 Sep 2022 14:21:18 -0400 Subject: [PATCH 19/43] conditional on internal variable SNOWALB --- .../GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 23 +++++++++++-------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index 0795918ae..c60a2f7f4 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -903,16 +903,6 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'snow_albedo' ,& - UNITS = '1' ,& - SHORT_NAME = 'SNOWALB' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartRequired ,& - RC=STATUS ) - VERIFY_(STATUS) call MAPL_AddInternalSpec(GC ,& LONG_NAME = 'wetness_at_wilting_point' ,& @@ -1387,6 +1377,19 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) + if (SNOW_ALBEDO_INFO == 1) then + call MAPL_AddInternalSpec(GC ,& + LONG_NAME = 'snow_albedo' ,& + UNITS = '1' ,& + SHORT_NAME = 'SNOWALB' ,& + FRIENDLYTO = trim(COMP_NAME) ,& + DIMS = MAPL_DimsTileOnly ,& + VLOCATION = MAPL_VLocationNone ,& + RESTART = MAPL_RestartRequired ,& + RC=STATUS ) + VERIFY_(STATUS) + endif + call MAPL_AddInternalSpec(GC ,& LONG_NAME = 'surface_heat_exchange_coefficient',& UNITS = 'kg m-2 s-1' ,& From 7f796b38509dadcb640be72d49f336beea405198 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Fri, 16 Sep 2022 14:27:18 -0400 Subject: [PATCH 20/43] typo --- .../GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index 65c84e154..53f000b18 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -1466,17 +1466,6 @@ subroutine SetServices ( GC, RC ) RC=STATUS ) VERIFY_(STATUS) - call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'snow_albedo' ,& - UNITS = '1' ,& - SHORT_NAME = 'SNOWALB' ,& - FRIENDLYTO = trim(COMP_NAME) ,& - DIMS = MAPL_DimsTileOnly ,& - VLOCATION = MAPL_VLocationNone ,& - RESTART = MAPL_RestartOptional ,& - RC=STATUS ) - VERIFY_(STATUS) - !---------- GOSWIM snow impurity related variables ---------- if (N_CONST_LAND4SNWALB /= 0) then From 0da8643eed7b861b09015341ade81b4122b6f7dd Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Mon, 19 Sep 2022 07:55:58 -0400 Subject: [PATCH 21/43] change missing value from -9999 to MAPL_UNDEF --- .../GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 1 + .../Utils/Raster/mod_process_hres_data.F90 | 7 ++++--- .../Utils/mk_restarts/CatchmentRst.F90 | 3 +++ 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index 53f000b18..75d06bc1b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -1381,6 +1381,7 @@ subroutine SetServices ( GC, RC ) LONG_NAME = 'snow_albedo' ,& UNITS = '1' ,& SHORT_NAME = 'SNOWALB' ,& + DEFAULT = MAPL_UNDEF ,& FRIENDLYTO = trim(COMP_NAME) ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index 92b3b161f..fda5b3424 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -26,6 +26,7 @@ MODULE process_hres_data use leap_year use MAPL_ConstantsMod use lsm_routines, ONLY: sibalb +use MAPL_Base, ONLY: MAPL_UNDEF #if defined USE_EXTERNAL_FINDLOC use findloc_mod, only: findloc @@ -3036,7 +3037,7 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) allocate (snw_alb(1:maxcat)) ! Start by setting all snow albedo values to missing - snw_alb(:)=-9999.0 + snw_alb(:)=MAPL_UNDEF do n = 1, maxcat read (10,*) tindex1,pfaf1,minlon,maxlon,minlat,maxlat @@ -3160,7 +3161,7 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) ! Calculate snow albedo for the current tile snw_alb(n) = sno_alb_sum / max(1.0,sno_alb_cnt) - if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=-9999.0 !1.E15 + if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=MAPL_UNDEF !1.E15 ! If no valid solution found, and if tile size smaller than snow albedo resolution, ! expand search area by 1-tile padding. @@ -3199,7 +3200,7 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) end do ! hhtil snw_alb(n) = sno_alb_sum2 / max(1.0,sno_alb_cnt2) - if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=-9999.0 !1.E15 + if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=MAPL_UNDEF !1.E15 endif diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index eca9d5e66..d07c0629b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -3,6 +3,7 @@ module CatchmentRstMod use mk_restarts_getidsMod use MAPL + use MAPL_Base, ONLY: MAPL_UNDEF use mpi use LSM_ROUTINES, ONLY: & catch_calc_soil_moist, & @@ -14,6 +15,7 @@ module CatchmentRstMod DZGT => CATCH_DZGT, & PEATCLSM_POROS_THRESHOLD + implicit none #ifndef __GFORTRAN__ integer :: ftell @@ -631,6 +633,7 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) var = Variable(type=pFIO_REAL32, dimensions='tile') call var%add_attribute('long_name', 'snow_albedo') call var%add_attribute('units', '1') + call var%add_attribute('_FillValue', MAPL_UNDEF) call this%meta%add_variable('SNOWALB', var) endif endif From 58093dfd126652243584753806e442ac299e124c Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Mon, 19 Sep 2022 11:42:19 -0400 Subject: [PATCH 22/43] restarts have _NoFill = true we can't add fill value --- .../GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index d07c0629b..40b9b03ef 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -633,7 +633,6 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) var = Variable(type=pFIO_REAL32, dimensions='tile') call var%add_attribute('long_name', 'snow_albedo') call var%add_attribute('units', '1') - call var%add_attribute('_FillValue', MAPL_UNDEF) call this%meta%add_variable('SNOWALB', var) endif endif From 837949d6edcb42468fd6082ec5a10f7ae9eb29d3 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Mon, 19 Sep 2022 11:43:30 -0400 Subject: [PATCH 23/43] must be done some other way --- .../GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index 75d06bc1b..53f000b18 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -1381,7 +1381,6 @@ subroutine SetServices ( GC, RC ) LONG_NAME = 'snow_albedo' ,& UNITS = '1' ,& SHORT_NAME = 'SNOWALB' ,& - DEFAULT = MAPL_UNDEF ,& FRIENDLYTO = trim(COMP_NAME) ,& DIMS = MAPL_DimsTileOnly ,& VLOCATION = MAPL_VLocationNone ,& From 25534c71a4a73b5254f4b5dd2b85a5493f72c55c Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Wed, 21 Sep 2022 11:10:41 -0400 Subject: [PATCH 24/43] better choice for descriptive name --- .../GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index 53f000b18..1c421158b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -1378,7 +1378,7 @@ subroutine SetServices ( GC, RC ) if (SNOW_ALBEDO_INFO == 1) then call MAPL_AddInternalSpec(GC ,& - LONG_NAME = 'snow_albedo' ,& + LONG_NAME = 'effective_snow_albedo' ,& UNITS = '1' ,& SHORT_NAME = 'SNOWALB' ,& FRIENDLYTO = trim(COMP_NAME) ,& From ec240f87cd9e7d75e3d9d2070c1533776a20876d Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Wed, 26 Oct 2022 20:51:37 -0400 Subject: [PATCH 25/43] add no data values; update comments --- .../GEOS_CatchCNGridComp.F90 | 2 +- .../GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 2 +- .../Shared/GEOS_SurfaceGridComp.rc | 2 +- .../Utils/Raster/mod_process_hres_data.F90 | 46 +++++++++---------- .../Utils/Raster/rmTinyCatchParaMod.F90 | 2 +- 5 files changed, 27 insertions(+), 27 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 index f30b0328c..4bc2943fc 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 @@ -86,7 +86,7 @@ subroutine SetServices ( GC, RC ) ! SNOW ALBEDO -- so far, only parameterization based on look-up table is implemented for CatchCN ! 0 : parameterization based on look-up table - ! 1 : MODIS-derived snow albedo (where available, elsewhere fall back to option 0) + ! 1 : MODIS-derived snow albedo (where available, elsewhere no data values filled with global land average Snow_Albedo=0.56) call ESMF_ConfigGetAttribute (SCF, label='SNOW_ALBEDO_INFO:', value=SNOW_ALBEDO_INFO, DEFAULT=0, RC=STATUS) ; VERIFY_(STATUS) _ASSERT( SNOW_ALBEDO_INFO==0, "SNOW_ALBEDO_INFO must be 0 for CatchCN") diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index bf6c47ead..c4b6f33c9 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -230,7 +230,7 @@ subroutine SetServices ( GC, RC ) ! SNOW ALBEDO ! 0 : parameterization based on look-up table - ! 1 : MODIS-derived snow albedo (where available, elsewhere fall back to option 0) + ! 1 : MODIS-derived snow albedo (where available, elsewhere no data values filled with global land average Snow_Albedo=0.56) call MAPL_GetResource (SCF, SNOW_ALBEDO_INFO, label='SNOW_ALBEDO_INFO:', DEFAULT=0, __RC__ ) ! GOSWIM SNOW_ALBEDO diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc index 3491c4f5c..282985451 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc @@ -118,7 +118,7 @@ # # 0 : Snow albedo parameterization based on look-up table (default) # 1 : Snow albedo derived from MODIS Collection MOD10A1.061 (Feb/2000 - Mar/2022) [NOT YET AVAILABLE FOR CatchCN] -# where available, elsewhere fall back to option 0 +# where available, elsewhere no data values filled with global land average Snow_Albedo=0.56 # # GEOSagcm=>SNOW_ALBEDO_INFO: 0 # GEOSldas=>SNOW_ALBEDO_INFO: 0 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index fda5b3424..ebab7f8cc 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -3010,15 +3010,13 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) real,allocatable,dimension(:) :: min_lon,max_lon,min_lat,max_lat,snw_alb integer(kind=4),parameter :: xdim = 1200, ydim = 1200 real,parameter :: alb_res=10.0/1200.0 - real,dimension(xdim) :: lon_alb - real,dimension(ydim) :: lat_alb real,dimension(xdim,ydim) :: stch_snw_alb_tmp real,dimension(36,18,xdim,ydim) :: stch_snw_alb real :: minlon,maxlon,minlat,maxlat,pad_lon,pad_lat real :: sno_alb_cnt,sno_alb_sum,sno_alb_cnt2,sno_alb_sum2 integer :: vvtil_min,hhtil_min,vvtil_max,hhtil_max,hhtil,vvtil integer :: tindex1,pfaf1 - integer(kind=4) :: dummy,varid1,varid2,varid3 + integer(kind=4) :: dummy,varid1 integer(kind=4) :: imin,imax,jmin,jmax integer(kind=4) :: imin2,imax2,jmin2,jmax2,count_init_invalid logical :: file_exists @@ -3073,18 +3071,20 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) write(vv,'(i2.2)') vvtil write(hh,'(i2.2)') hhtil - fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v1/snow_alb_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' + ! MODIS-based climatology albedo raster files filled for NoData values using global land + ! average Snow Albedo (0.56). Note: the average excludes Antarctica and Greenland ice + ! sheets and is weighted by the grid-cell area. + fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v2/snow_alb_FillVal_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' + + ! MODIS-based climatology albedo raster files (NoData set to 1e+15) + !fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v1/snow_alb_noFill_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' ! Open the file. (NF90_NOWRITE ensures read-only access to the file) status=NF_OPEN(trim(fname),NF_NOWRITE, ncid) ; VERIFY_(STATUS) ! Based on vars name, get the varids. status=NF_INQ_VARID(ncid,'Snow_Albedo',VarID1) ; VERIFY_(STATUS) - status=NF_INQ_VARID(ncid,'lon' ,VarID2) ; VERIFY_(STATUS) - status=NF_INQ_VARID(ncid,'lat' ,VarID3) ; VERIFY_(STATUS) ! Read the data. status=NF_GET_VARA_REAL(ncid,VarID1,(/1,1/),(/xdim,ydim/),stch_snw_alb_tmp) ; VERIFY_(STATUS) - status=NF_GET_VARA_REAL(ncid,VarID2,(/1/) ,(/xdim/) ,lon_alb) ; VERIFY_(STATUS) - status=NF_GET_VARA_REAL(ncid,VarID3,(/1/) ,(/ydim/) ,lat_alb) ; VERIFY_(STATUS) ! Close the file, freeing all resources. status=NF_CLOSE(ncid); VERIFY_(STATUS) @@ -3096,7 +3096,7 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) ! loop over tiles print*, 'Starting tile loop for snow albedo. ' - count_init_invalid=0 ! counter for non-valid snow albedo avalues (informational use only; not needed for) + count_init_invalid=0 ! counter for non-valid snow albedo values (informational use only) do n = 1, maxcat ! loop over tiles @@ -3147,14 +3147,14 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) jmax=min(jmax,ydim) ! Generate sums and counts using current tile corresponding indices - sno_alb_sum= sno_alb_sum + & - sum(stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax), & - stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & - stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax).le.1.0) + sno_alb_sum= sno_alb_sum + & + sum(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax), & + stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & + stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).le.1.0) - sno_alb_cnt= sno_alb_cnt + & - count(stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & - stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin:imax,jmin:jmax).le.1.0) + sno_alb_cnt= sno_alb_cnt + & + count(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & + stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).le.1.0) end do ! vvtil end do ! hhtil @@ -3187,14 +3187,14 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) jmin2=max(jmin2,1) jmax2=min(jmax2,ydim) - sno_alb_sum2= sno_alb_sum2 + & - sum(stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2), & - stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & - stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2).le.1.0) + sno_alb_sum2= sno_alb_sum2 + & + sum(stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2), & + stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & + stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).le.1.0) - sno_alb_cnt2= sno_alb_cnt2 + & - count(stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & - stch_snw_alb(hhtil:hhtil,vvtil:vvtil,imin2:imax2,jmin2:jmax2).le.1.0) + sno_alb_cnt2= sno_alb_cnt2 + & + count(stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & + stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).le.1.0) end do ! vvtil end do ! hhtil diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index 992bb7310..ac191c7cb 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -52,7 +52,7 @@ module rmTinyCatchParaMod logical, public, save :: use_PEATMAP = .false. logical, public, save :: jpl_height = .false. character*8, public, save :: LAIBCS = 'UNDEF' - character*4, public, save :: SOILBCS = 'UNDEF' + character*6, public, save :: SOILBCS = 'UNDEF' character*6, public, save :: MODALB = 'UNDEF' character*8, public, save :: SNOWALB = 'UNDEF' REAL, public, save :: GNU = MAPL_UNDEF From 49e1f246e7860abcff58eb273963afcce1e3c860 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Sun, 30 Oct 2022 19:25:22 -0400 Subject: [PATCH 26/43] adding case for Scott's testing --- .../GEOSsurface_GridComp/Utils/Raster/make_bcs | 4 +++- .../Utils/Raster/mkCatchParam.F90 | 18 +++++++++--------- .../Utils/Raster/rmTinyCatchParaMod.F90 | 9 +++++++++ 3 files changed, 21 insertions(+), 10 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs index 2253db27f..af40fb04b 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs @@ -115,6 +115,7 @@ echo " ${C2} = NL3 + JPL veg height + PEATMAP${CR}" echo " ${C2}v06 = NL3 + JPL veg height + PEATMAP + MODIS snow alb${CR}" echo " ${C2}v07 = NL3 + PEATMAP${CR}" echo " ${C2}v08 = NL3 + MODIS snow alb${CR}" +echo " ${C2}v09 = NL3 + PEATMAP + MODIS snow alb${CR}" echo " " if ( $HELPMODE != YES ) then echo " ${C1} *BCs produced by this code will differ from BCs in archived directories\!\!\! ${CR}" @@ -142,7 +143,8 @@ if ( $HELPMODE != YES ) then $dummy == 'NL5' | \ $dummy == 'v06' | \ $dummy == 'v07' | \ - $dummy == 'v08') then + $dummy == 'v08' | \ + $dummy == 'v09') then set lbcsv = $dummy else if ( $dummy == '' ) then echo $lbcsv diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index b07b9ba3b..b88cd2fa2 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -9,7 +9,7 @@ PROGRAM mkCatchParam ! -y: Size of latitude dimension of input raster. DEFAULT: 4320 ! -b: position of the dateline in the first box. DEFAULT: DC ! -g: Gridname (name of the .til or .rst file without file extension) -! -v: LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08) +! -v: LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09) ! -e: EASE : This is optional if catchment.def file is available already or ! the til file format is pre-Fortuna-2. ! @@ -103,14 +103,14 @@ PROGRAM mkCatchParam ! call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH') ! call execute_command_line('cd ..') - USAGE(1) ="Usage: mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV -e EASE " - USAGE(2) =" -x: Size of longitude dimension of input raster. DEFAULT: 8640 " - USAGE(3) =" -y: Size of latitude dimension of input raster. DEFAULT: 4320 " - USAGE(4) =" -g: Gridname (name of the .til or .rst file without file extension) " - USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " - USAGE(6) =" -e: EASE : This is optional if catchment.def file is available already or " - USAGE(7) =" the til file format is pre-Fortuna-2. " - USAGE(8) =" -v LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08) " + USAGE(1) ="Usage: mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV -e EASE " + USAGE(2) =" -x: Size of longitude dimension of input raster. DEFAULT: 8640 " + USAGE(3) =" -y: Size of latitude dimension of input raster. DEFAULT: 4320 " + USAGE(4) =" -g: Gridname (name of the .til or .rst file without file extension) " + USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " + USAGE(6) =" -e: EASE : This is optional if catchment.def file is available already or " + USAGE(7) =" the til file format is pre-Fortuna-2. " + USAGE(8) =" -v LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09) " ! Process Arguments !------------------ diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index ac191c7cb..f7d47a432 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -168,6 +168,15 @@ SUBROUTINE init_bcs_config (LBCSV) use_PEATMAP = .false. jpl_height = .false. + case ("v09") + LAIBCS = 'MODGEO' + SOILBCS = 'HWSD' + MODALB = 'MODIS2' + SNOWALB = 'MODC061' + GNU = 1.0 + use_PEATMAP = .true. + jpl_height = .false. + end select END SUBROUTINE init_bcs_config From 7d22313d4cdd6b30bb8114a894d62e51c14fe280 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Tue, 1 Nov 2022 07:54:45 -0400 Subject: [PATCH 27/43] fix for CF high res bcs runs --- .../Utils/Raster/mod_process_hres_data.F90 | 23 +++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index ebab7f8cc..555120adb 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -3100,6 +3100,29 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) do n = 1, maxcat ! loop over tiles + ! note: if a tile is within 30 arcsec of the date line, snow albedo value + ! will remain missing (no values between the edge of the map and tile). + ! There are very few tiles that fall into this gap. To save the computational + ! time required for looping over h01 and h36 files, those tiles are 'moved' + ! ~30 arc sec away from the date line, which will result in the assignment of + ! "neighboring" snow albedo value + if (min_lon(n) .lt. -179.99 .or. max_lon(n) .lt. -179.99) then + min_lon(n)=min_lon(n)+0.01 + max_lon(n)=max_lon(n)+0.01 + endif + if (min_lon(n) .gt. 179.99 .or. max_lon(n) .gt. 179.99) then + min_lon(n)=min_lon(n)-0.01 + max_lon(n)=max_lon(n)-0.01 + endif + if (min_lat(n) .lt. -89.99 .or. max_lat(n) .lt. -89.99) then + min_lat(n)=min_lat(n)+0.01 + max_lat(n)=max_lat(n)+0.01 + endif + if (min_lat(n) .gt. 89.99 .or. max_lat(n) .gt. 89.99) then + min_lat(n)=min_lat(n)-0.01 + max_lat(n)=max_lat(n)-0.01 + endif + ! Set sums and counts to zero sno_alb_sum =0. sno_alb_cnt =0. From 4086fd3a92f51d90c8819e143c524fe467fb12d4 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Tue, 8 Nov 2022 08:40:23 -0500 Subject: [PATCH 28/43] fix indexing --- .../Utils/Raster/mod_process_hres_data.F90 | 52 ++++++------------- 1 file changed, 17 insertions(+), 35 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index 555120adb..ff39c6c7b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -3100,29 +3100,6 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) do n = 1, maxcat ! loop over tiles - ! note: if a tile is within 30 arcsec of the date line, snow albedo value - ! will remain missing (no values between the edge of the map and tile). - ! There are very few tiles that fall into this gap. To save the computational - ! time required for looping over h01 and h36 files, those tiles are 'moved' - ! ~30 arc sec away from the date line, which will result in the assignment of - ! "neighboring" snow albedo value - if (min_lon(n) .lt. -179.99 .or. max_lon(n) .lt. -179.99) then - min_lon(n)=min_lon(n)+0.01 - max_lon(n)=max_lon(n)+0.01 - endif - if (min_lon(n) .gt. 179.99 .or. max_lon(n) .gt. 179.99) then - min_lon(n)=min_lon(n)-0.01 - max_lon(n)=max_lon(n)-0.01 - endif - if (min_lat(n) .lt. -89.99 .or. max_lat(n) .lt. -89.99) then - min_lat(n)=min_lat(n)+0.01 - max_lat(n)=max_lat(n)+0.01 - endif - if (min_lat(n) .gt. 89.99 .or. max_lat(n) .gt. 89.99) then - min_lat(n)=min_lat(n)-0.01 - max_lat(n)=max_lat(n)-0.01 - endif - ! Set sums and counts to zero sno_alb_sum =0. sno_alb_cnt =0. @@ -3131,10 +3108,15 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) ! Use tile's min/max lat/lon info to identify the 10x10deg input file(s) ! and read in snow albedo value(s). The "ceiling" and "floor" implements the "halo". - vvtil_min= floor((min_lat(n)+ 90.0)/10.) - hhtil_min= floor((min_lon(n)+180.0)/10.) - vvtil_max=ceiling((max_lat(n)+ 90.0)/10.) - hhtil_max=ceiling((max_lon(n)+180.0)/10.) + vvtil_min= floor((min_lat(n)+ 90.0)/10.)+1 + hhtil_min= floor((min_lon(n)+180.0)/10.)+1 + + ! if tile crosses the edge of the snow albedo 10x10deg box, expand the + ! search area into the neighbouring 10x10deg box + hhtil_max=hhtil_min + vvtil_max=vvtil_min + if (floor(min_lon(n)/10) .ne. floor(max_lon(n)/10)) hhtil_max=hhtil_min+1 + if (floor(min_lat(n)/10) .ne. floor(max_lat(n)/10)) vvtil_max=vvtil_min+1 ! Safety checks: ! 1. Make sure vv's and hh's are within the range. If min>max, swap them. @@ -3159,10 +3141,10 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction ! Find indices ranges corresponding to the current tile area. - imin=floor((min_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) - imax=floor((max_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) - jmin=floor((min_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) - jmax=floor((max_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) + imin=floor((min_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) +1 + imax=floor((max_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) +1 + jmin=floor((min_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) +1 + jmax=floor((max_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) +1 ! Keep within the range. imin=max(imin,1) imax=min(imax,xdim) @@ -3201,10 +3183,10 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction ! Repeat the steps for extracting snow albedo value - imin2=floor((min_lon(n)-pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) - imax2=floor((max_lon(n)+pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) - jmin2=floor((min_lat(n)-pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) - jmax2=floor((max_lat(n)+pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) + imin2=floor((min_lon(n)-pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0))+1 + imax2=floor((max_lon(n)+pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0))+1 + jmin2=floor((min_lat(n)-pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0))+1 + jmax2=floor((max_lat(n)+pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0))+1 imin2=max(imin2,1) imax2=min(imax2,xdim) jmin2=max(jmin2,1) From 91d0d5d93e510e7f00617eec43035d29a9707e9a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 10 Nov 2022 06:12:26 -0500 Subject: [PATCH 29/43] additional clean-up of snow albedo changes - removed obsolete checks for unphysical snow albedo values from Catch GC - editedd comments - fixed indentation in subroutine MODIS_snow_alb() --- .../GEOS_CatchCNGridComp.F90 | 2 +- .../GEOScatch_GridComp/GEOS_CatchGridComp.F90 | 36 +- .../Shared/GEOS_SurfaceGridComp.rc | 6 +- .../Utils/Raster/mod_process_hres_data.F90 | 446 +++++++++--------- .../Utils/mk_restarts/CatchmentRst.F90 | 6 +- 5 files changed, 250 insertions(+), 246 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 index 4bc2943fc..a2bbcb990 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOS_CatchCNGridComp.F90 @@ -86,7 +86,7 @@ subroutine SetServices ( GC, RC ) ! SNOW ALBEDO -- so far, only parameterization based on look-up table is implemented for CatchCN ! 0 : parameterization based on look-up table - ! 1 : MODIS-derived snow albedo (where available, elsewhere no data values filled with global land average Snow_Albedo=0.56) + ! 1 : MODIS-derived snow albedo (backfilled with global land average snow albedo) call ESMF_ConfigGetAttribute (SCF, label='SNOW_ALBEDO_INFO:', value=SNOW_ALBEDO_INFO, DEFAULT=0, RC=STATUS) ; VERIFY_(STATUS) _ASSERT( SNOW_ALBEDO_INFO==0, "SNOW_ALBEDO_INFO must be 0 for CatchCN") diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 index c4b6f33c9..257ac7bc7 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatch_GridComp/GEOS_CatchGridComp.F90 @@ -230,7 +230,7 @@ subroutine SetServices ( GC, RC ) ! SNOW ALBEDO ! 0 : parameterization based on look-up table - ! 1 : MODIS-derived snow albedo (where available, elsewhere no data values filled with global land average Snow_Albedo=0.56) + ! 1 : MODIS-derived snow albedo (backfilled with global land average snow albedo) call MAPL_GetResource (SCF, SNOW_ALBEDO_INFO, label='SNOW_ALBEDO_INFO:', DEFAULT=0, __RC__ ) ! GOSWIM SNOW_ALBEDO @@ -4877,16 +4877,17 @@ subroutine Driver ( RC ) if (SNOW_ALBEDO_INFO == 1) then - ! where available, use MODIS-derived snow albedo from bcs (via Catch restart) - + ! use MODIS-derived snow albedo from bcs (via Catch restart) + ! + ! as a restart parameter from the bcs, snow albedo must not have no-data-values + ! (checks for unphysical values should be in the make_bcs package) + call MAPL_GetPointer(INTERNAL,SNOWALB,'SNOWALB',RC=STATUS); VERIFY_(STATUS) - where (SNOWALB > 0. .and. SNOWALB <= 1.) - SNOVR = SNOWALB - SNONR = SNOWALB - SNOVF = SNOWALB - SNONF = SNOWALB - endwhere + SNOVR = SNOWALB + SNONR = SNOWALB + SNOVF = SNOWALB + SNONF = SNOWALB endif @@ -5570,14 +5571,15 @@ subroutine Driver ( RC ) if (SNOW_ALBEDO_INFO == 1) then - ! where available, use MODIS-derived snow albedo from bcs (via Catch restart) - - where (SNOWALB > 0. .and. SNOWALB <= 1.) - SNOVR = SNOWALB - SNONR = SNOWALB - SNOVF = SNOWALB - SNONF = SNOWALB - endwhere + ! use MODIS-derived snow albedo from bcs (via Catch restart) + ! + ! as a restart parameter from the bcs, snow albedo must not have no-data-values + ! (checks for unphysical values should be in the make_bcs package) + + SNOVR = SNOWALB + SNONR = SNOWALB + SNOVF = SNOWALB + SNONF = SNOWALB endif diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc index 282985451..1c11a0be2 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc @@ -117,8 +117,10 @@ # ---- SNOW ALBEDO # # 0 : Snow albedo parameterization based on look-up table (default) -# 1 : Snow albedo derived from MODIS Collection MOD10A1.061 (Feb/2000 - Mar/2022) [NOT YET AVAILABLE FOR CatchCN] -# where available, elsewhere no data values filled with global land average Snow_Albedo=0.56 +# 1 : Snow albedo derived from MODIS Collection MOD10A1.061 (Feb/2000 - Mar/2022) +# - backfilled with global land average snow albedo where unavailable +# - must use compatible bcs version (e.g., v06, v08, v09) +# - NOT YET AVAILABLE FOR CatchCN # # GEOSagcm=>SNOW_ALBEDO_INFO: 0 # GEOSldas=>SNOW_ALBEDO_INFO: 0 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index ff39c6c7b..8630ca216 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -2995,238 +2995,238 @@ END SUBROUTINE hres_gswp2 SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) - ! Process static snow albedo calculated from MODIS climatology and write into clsm/catch_params.nc4 - ! - ! Biljana Orescanin July 2022, SSAI@NASA - - implicit none - integer, intent (in) :: nx, ny - character(*) :: gfiler - integer,allocatable,target,dimension (:,:) :: tile_id - - character*200 :: fname - character*2 :: vv,hh - integer :: n,maxcat,j,ncid,status - real,allocatable,dimension(:) :: min_lon,max_lon,min_lat,max_lat,snw_alb - integer(kind=4),parameter :: xdim = 1200, ydim = 1200 - real,parameter :: alb_res=10.0/1200.0 - real,dimension(xdim,ydim) :: stch_snw_alb_tmp - real,dimension(36,18,xdim,ydim) :: stch_snw_alb - real :: minlon,maxlon,minlat,maxlat,pad_lon,pad_lat - real :: sno_alb_cnt,sno_alb_sum,sno_alb_cnt2,sno_alb_sum2 - integer :: vvtil_min,hhtil_min,vvtil_max,hhtil_max,hhtil,vvtil - integer :: tindex1,pfaf1 - integer(kind=4) :: dummy,varid1 - integer(kind=4) :: imin,imax,jmin,jmax - integer(kind=4) :: imin2,imax2,jmin2,jmax2,count_init_invalid - logical :: file_exists - - ! Read number of catchment-tiles (maxcat) from catchment.def file - fname='clsm/catchment.def' - open (10,file=fname,status='old',action='read',form='formatted') - read(10,*) maxcat - - ! Read min/max lat/lons to use when locating snow albedo grids in - ! the stitched MODIS albedo file - allocate (min_lon(1:maxcat)) - allocate (min_lat(1:maxcat)) - allocate (max_lon(1:maxcat)) - allocate (max_lat(1:maxcat)) - allocate (snw_alb(1:maxcat)) - - ! Start by setting all snow albedo values to missing - snw_alb(:)=MAPL_UNDEF - - do n = 1, maxcat - read (10,*) tindex1,pfaf1,minlon,maxlon,minlat,maxlat - min_lon(n) = minlon - max_lon(n) = maxlon - min_lat(n) = minlat - max_lat(n) = maxlat - end do - - close (10,status='keep') - - ! Read tile-id raster file - allocate(tile_id(1:nx,1:ny)) - - fname=trim(gfiler)//'.rst' - open (10,file=fname,status='old',action='read', & - form='unformatted',convert='little_endian') - - do j=1,ny - read(10)tile_id(:,j) - end do - - close (10,status='keep') - - !------------ Get the information on snow albedo ----- - ! ----------- The information on snow albedo is stored in 10x10deg 30-arcsec resolution files. - ! ----------- Read in this information, then loop over the tiles to find a corresponding snow albedo. - - ! Read in all 10x10deg snow albedo files into a single [36,18,1200,1200] array - do hhtil=1,36 ! loop over input files - horizontal direction - do vvtil=1,18 ! loop over input files - vertical direction - - write(vv,'(i2.2)') vvtil - write(hh,'(i2.2)') hhtil - - ! MODIS-based climatology albedo raster files filled for NoData values using global land - ! average Snow Albedo (0.56). Note: the average excludes Antarctica and Greenland ice - ! sheets and is weighted by the grid-cell area. - fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v2/snow_alb_FillVal_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' - - ! MODIS-based climatology albedo raster files (NoData set to 1e+15) - !fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v1/snow_alb_noFill_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' - - ! Open the file. (NF90_NOWRITE ensures read-only access to the file) - status=NF_OPEN(trim(fname),NF_NOWRITE, ncid) ; VERIFY_(STATUS) - ! Based on vars name, get the varids. - status=NF_INQ_VARID(ncid,'Snow_Albedo',VarID1) ; VERIFY_(STATUS) - ! Read the data. - status=NF_GET_VARA_REAL(ncid,VarID1,(/1,1/),(/xdim,ydim/),stch_snw_alb_tmp) ; VERIFY_(STATUS) - ! Close the file, freeing all resources. - status=NF_CLOSE(ncid); VERIFY_(STATUS) - - ! Store snow albedo values into a single 4D aray - stch_snw_alb(hhtil,vvtil,:,:)=stch_snw_alb_tmp - + ! Process static snow albedo calculated from MODIS climatology and write into clsm/catch_params.nc4 + ! + ! Biljana Orescanin July 2022, SSAI@NASA + + implicit none + integer, intent (in) :: nx, ny + character(*) :: gfiler + integer,allocatable,target,dimension (:,:) :: tile_id + + character*200 :: fname + character*2 :: vv,hh + integer :: n,maxcat,j,ncid,status + real,allocatable,dimension(:) :: min_lon,max_lon,min_lat,max_lat,snw_alb + integer(kind=4),parameter :: xdim = 1200, ydim = 1200 + real,parameter :: alb_res=10.0/1200.0 + real,dimension(xdim,ydim) :: stch_snw_alb_tmp + real,dimension(36,18,xdim,ydim) :: stch_snw_alb + real :: minlon,maxlon,minlat,maxlat,pad_lon,pad_lat + real :: sno_alb_cnt,sno_alb_sum,sno_alb_cnt2,sno_alb_sum2 + integer :: vvtil_min,hhtil_min,vvtil_max,hhtil_max,hhtil,vvtil + integer :: tindex1,pfaf1 + integer(kind=4) :: dummy,varid1 + integer(kind=4) :: imin,imax,jmin,jmax + integer(kind=4) :: imin2,imax2,jmin2,jmax2,count_init_invalid + logical :: file_exists + + ! Read number of catchment-tiles (maxcat) from catchment.def file + fname='clsm/catchment.def' + open (10,file=fname,status='old',action='read',form='formatted') + read(10,*) maxcat + + ! Read min/max lat/lons to use when locating snow albedo grids in + ! the stitched MODIS albedo file + allocate (min_lon(1:maxcat)) + allocate (min_lat(1:maxcat)) + allocate (max_lon(1:maxcat)) + allocate (max_lat(1:maxcat)) + allocate (snw_alb(1:maxcat)) + + ! Start by setting all snow albedo values to missing + snw_alb(:)=MAPL_UNDEF + + do n = 1, maxcat + read (10,*) tindex1,pfaf1,minlon,maxlon,minlat,maxlat + min_lon(n) = minlon + max_lon(n) = maxlon + min_lat(n) = minlat + max_lat(n) = maxlat + end do + + close (10,status='keep') + + ! Read tile-id raster file + allocate(tile_id(1:nx,1:ny)) + + fname=trim(gfiler)//'.rst' + open (10,file=fname,status='old',action='read', & + form='unformatted',convert='little_endian') + + do j=1,ny + read(10)tile_id(:,j) + end do + + close (10,status='keep') + + !------------ Get the information on snow albedo ----- + ! ----------- The information on snow albedo is stored in 10x10deg 30-arcsec resolution files. + ! ----------- Read in this information, then loop over the tiles to find a corresponding snow albedo. + + ! Read in all 10x10deg snow albedo files into a single [36,18,1200,1200] array + do hhtil=1,36 ! loop over input files - horizontal direction + do vvtil=1,18 ! loop over input files - vertical direction + + write(vv,'(i2.2)') vvtil + write(hh,'(i2.2)') hhtil + + ! MODIS-based climatology albedo raster files, backfilled with global land + ! average snow albedo (=0.56; average excludes Antarctica and Greenland ice + ! sheets and is weighted by the grid-cell area). + fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v2/snow_alb_FillVal_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' + + ! MODIS-based climatology albedo raster files (NoData set to 1e+15) + !fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v1/snow_alb_noFill_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' + + ! Open the file. (NF90_NOWRITE ensures read-only access to the file) + status=NF_OPEN(trim(fname),NF_NOWRITE, ncid) ; VERIFY_(STATUS) + ! Based on vars name, get the varids. + status=NF_INQ_VARID(ncid,'Snow_Albedo',VarID1) ; VERIFY_(STATUS) + ! Read the data. + status=NF_GET_VARA_REAL(ncid,VarID1,(/1,1/),(/xdim,ydim/),stch_snw_alb_tmp) ; VERIFY_(STATUS) + ! Close the file, freeing all resources. + status=NF_CLOSE(ncid); VERIFY_(STATUS) + + ! Store snow albedo values into a single 4D aray + stch_snw_alb(hhtil,vvtil,:,:)=stch_snw_alb_tmp + + enddo enddo - enddo - - ! loop over tiles - print*, 'Starting tile loop for snow albedo. ' - count_init_invalid=0 ! counter for non-valid snow albedo values (informational use only) - - do n = 1, maxcat ! loop over tiles - - ! Set sums and counts to zero - sno_alb_sum =0. - sno_alb_cnt =0. - sno_alb_sum2=0. - sno_alb_cnt2=0. - - ! Use tile's min/max lat/lon info to identify the 10x10deg input file(s) - ! and read in snow albedo value(s). The "ceiling" and "floor" implements the "halo". - vvtil_min= floor((min_lat(n)+ 90.0)/10.)+1 - hhtil_min= floor((min_lon(n)+180.0)/10.)+1 - - ! if tile crosses the edge of the snow albedo 10x10deg box, expand the - ! search area into the neighbouring 10x10deg box - hhtil_max=hhtil_min - vvtil_max=vvtil_min - if (floor(min_lon(n)/10) .ne. floor(max_lon(n)/10)) hhtil_max=hhtil_min+1 - if (floor(min_lat(n)/10) .ne. floor(max_lat(n)/10)) vvtil_max=vvtil_min+1 - - ! Safety checks: - ! 1. Make sure vv's and hh's are within the range. If min>max, swap them. - if (vvtil_min .gt. vvtil_max) then - dummy =vvtil_min - vvtil_min=vvtil_max - vvtil_max=dummy - endif - if (hhtil_min .gt. hhtil_max) then - dummy =hhtil_min - hhtil_min=hhtil_max - hhtil_max=dummy - endif - - ! 2. Keep within the range. - vvtil_min=max(vvtil_min,1) - vvtil_max=min(vvtil_max,18) - hhtil_min=max(hhtil_min,1) - hhtil_max=min(hhtil_max,36) - - do hhtil=hhtil_min,hhtil_max ! loop through input files - horizontal direction - do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction - - ! Find indices ranges corresponding to the current tile area. - imin=floor((min_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) +1 - imax=floor((max_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) +1 - jmin=floor((min_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) +1 - jmax=floor((max_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) +1 - ! Keep within the range. - imin=max(imin,1) - imax=min(imax,xdim) - jmin=max(jmin,1) - jmax=min(jmax,ydim) - - ! Generate sums and counts using current tile corresponding indices - sno_alb_sum= sno_alb_sum + & - sum(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax), & - stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & - stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).le.1.0) - - sno_alb_cnt= sno_alb_cnt + & - count(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & - stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).le.1.0) - - end do ! vvtil - end do ! hhtil - - ! Calculate snow albedo for the current tile - snw_alb(n) = sno_alb_sum / max(1.0,sno_alb_cnt) - if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=MAPL_UNDEF !1.E15 - - ! If no valid solution found, and if tile size smaller than snow albedo resolution, - ! expand search area by 1-tile padding. - - ! Size of a tile (in both directions) - pad_lon=(max_lon(n)-min_lon(n)) - pad_lat=(max_lat(n)-min_lat(n)) - - if (snw_alb(n) .le. 0.0 .and. (pad_lon .lt. alb_res .or. pad_lat .lt. alb_res)) then - - count_init_invalid=count_init_invalid+1 - - do hhtil=hhtil_min,hhtil_max ! loop through input files - horizontal direction - do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction - - ! Repeat the steps for extracting snow albedo value - imin2=floor((min_lon(n)-pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0))+1 - imax2=floor((max_lon(n)+pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0))+1 - jmin2=floor((min_lat(n)-pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0))+1 - jmax2=floor((max_lat(n)+pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0))+1 - imin2=max(imin2,1) - imax2=min(imax2,xdim) - jmin2=max(jmin2,1) - jmax2=min(jmax2,ydim) - - sno_alb_sum2= sno_alb_sum2 + & - sum(stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2), & - stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & - stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).le.1.0) - - sno_alb_cnt2= sno_alb_cnt2 + & - count(stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & - stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).le.1.0) - - end do ! vvtil - end do ! hhtil - - snw_alb(n) = sno_alb_sum2 / max(1.0,sno_alb_cnt2) - if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=MAPL_UNDEF !1.E15 - - endif - - end do ! n-loop over tiles - - ! write snow albedo into clsm/catch_params.nc4 + + ! loop over tiles + print*, 'Starting tile loop for snow albedo. ' + count_init_invalid=0 ! counter for non-valid snow albedo values (informational use only) + + do n = 1, maxcat ! loop over tiles + + ! Set sums and counts to zero + sno_alb_sum =0. + sno_alb_cnt =0. + sno_alb_sum2=0. + sno_alb_cnt2=0. + + ! Use tile's min/max lat/lon info to identify the 10x10deg input file(s) + ! and read in snow albedo value(s). The "ceiling" and "floor" implements the "halo". + vvtil_min= floor((min_lat(n)+ 90.0)/10.)+1 + hhtil_min= floor((min_lon(n)+180.0)/10.)+1 + + ! if tile crosses the edge of the snow albedo 10x10deg box, expand the + ! search area into the neighbouring 10x10deg box + hhtil_max=hhtil_min + vvtil_max=vvtil_min + if (floor(min_lon(n)/10) .ne. floor(max_lon(n)/10)) hhtil_max=hhtil_min+1 + if (floor(min_lat(n)/10) .ne. floor(max_lat(n)/10)) vvtil_max=vvtil_min+1 + + ! Safety checks: + ! 1. Make sure vv's and hh's are within the range. If min>max, swap them. + if (vvtil_min .gt. vvtil_max) then + dummy =vvtil_min + vvtil_min=vvtil_max + vvtil_max=dummy + endif + if (hhtil_min .gt. hhtil_max) then + dummy =hhtil_min + hhtil_min=hhtil_max + hhtil_max=dummy + endif + + ! 2. Keep within the range. + vvtil_min=max(vvtil_min,1) + vvtil_max=min(vvtil_max,18) + hhtil_min=max(hhtil_min,1) + hhtil_max=min(hhtil_max,36) + + do hhtil=hhtil_min,hhtil_max ! loop through input files - horizontal direction + do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction + + ! Find indices ranges corresponding to the current tile area. + imin=floor((min_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) +1 + imax=floor((max_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) +1 + jmin=floor((min_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) +1 + jmax=floor((max_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) +1 + ! Keep within the range. + imin=max(imin,1) + imax=min(imax,xdim) + jmin=max(jmin,1) + jmax=min(jmax,ydim) + + ! Generate sums and counts using current tile corresponding indices + sno_alb_sum= sno_alb_sum + & + sum(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax), & + stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & + stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).le.1.0) + + sno_alb_cnt= sno_alb_cnt + & + count(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & + stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).le.1.0) + + end do ! vvtil + end do ! hhtil + + ! Calculate snow albedo for the current tile + snw_alb(n) = sno_alb_sum / max(1.0,sno_alb_cnt) + if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=MAPL_UNDEF !1.E15 + + ! If no valid solution found, and if tile size smaller than snow albedo resolution, + ! expand search area by 1-tile padding. + + ! Size of a tile (in both directions) + pad_lon=(max_lon(n)-min_lon(n)) + pad_lat=(max_lat(n)-min_lat(n)) + + if (snw_alb(n) .le. 0.0 .and. (pad_lon .lt. alb_res .or. pad_lat .lt. alb_res)) then + + count_init_invalid=count_init_invalid+1 + + do hhtil=hhtil_min,hhtil_max ! loop through input files - horizontal direction + do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction + + ! Repeat the steps for extracting snow albedo value + imin2=floor((min_lon(n)-pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0))+1 + imax2=floor((max_lon(n)+pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0))+1 + jmin2=floor((min_lat(n)-pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0))+1 + jmax2=floor((max_lat(n)+pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0))+1 + imin2=max(imin2,1) + imax2=min(imax2,xdim) + jmin2=max(jmin2,1) + jmax2=min(jmax2,ydim) + + sno_alb_sum2= sno_alb_sum2 + & + sum(stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2), & + stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & + stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).le.1.0) + + sno_alb_cnt2= sno_alb_cnt2 + & + count(stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & + stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).le.1.0) + + end do ! vvtil + end do ! hhtil + + snw_alb(n) = sno_alb_sum2 / max(1.0,sno_alb_cnt2) + if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=MAPL_UNDEF !1.E15 + + endif + + end do ! n-loop over tiles + + ! write snow albedo into clsm/catch_params.nc4 inquire(file='clsm/catch_params.nc4', exist=file_exists) - + if(file_exists) then status = NF_OPEN ('clsm/catch_params.nc4', NF_WRITE, ncid ) ; VERIFY_(STATUS) status = NF_PUT_VARA_REAL(NCID,NC_VarID(NCID,'SNOWALB'),(/1/),(/maxcat/),real(snw_alb)) ; VERIFY_(STATUS) STATUS = NF_CLOSE (NCID) ; VERIFY_(STATUS) endif - - print*, 'Ended tile loop for snow albedo. ' - print*, 'Initially found ', count_init_invalid, ' tiles with no-data values for snow albedo (out of ', maxcat,')' - + + print*, 'Ended tile loop for snow albedo. ' + print*, 'Initially found ', count_init_invalid, ' tiles with no-data values for snow albedo (out of ', maxcat,')' + END SUBROUTINE MODIS_snow_alb - + !-------------------------------------------------------------------------------------- - + SUBROUTINE soil_para_hwsd (nx,ny,gfiler) ! Processing NGDC-HWSD-STATSGO merged soil properties with Woesten Soil diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index 40b9b03ef..34406b852 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -336,7 +336,7 @@ subroutine read_shared_nc4(this, formatter, rc) call MAPL_VarRead(formatter,"WW",this%ww, __RC__) endif _RETURN(_SUCCESS) - end subroutine + end subroutine read_shared_nc4 subroutine write_nc4 (this, filename, rc) class(CatchmentRst), intent(inout):: this @@ -522,7 +522,7 @@ subroutine allocate_catch(this,rc) _RETURN(_SUCCESS) end subroutine allocate_catch - ! This subroutine reads BCs from BCSDIR and hydrological varable + ! This subroutine reads BCs from BCSDIR and hydrological variable (??) subroutine add_bcs_to_rst(this, surflay, DataDir, rc) class(CatchmentRst), intent(inout) :: this real, intent(in) :: surflay @@ -542,7 +542,7 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) type(Variable) :: var type(FileMetadata) :: meta_ - character*256 :: Iam = "add_bcs" + character*256 :: Iam = "add_bcs_to_rst" ntiles = this%ntiles From b9b84ffc94cb759ffa2336ee4534af98501037e7 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 10 Nov 2022 17:37:30 -0500 Subject: [PATCH 30/43] removed obsolete reading of tile_id from subroutine MODIS_snow_alb() --- .../Utils/Raster/mod_process_hres_data.F90 | 18 +----------------- 1 file changed, 1 insertion(+), 17 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index 8630ca216..bea879c9f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -2993,16 +2993,13 @@ END SUBROUTINE hres_gswp2 !---------------------------------------------------------------------- - SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) + SUBROUTINE MODIS_snow_alb ( ) ! Process static snow albedo calculated from MODIS climatology and write into clsm/catch_params.nc4 ! ! Biljana Orescanin July 2022, SSAI@NASA implicit none - integer, intent (in) :: nx, ny - character(*) :: gfiler - integer,allocatable,target,dimension (:,:) :: tile_id character*200 :: fname character*2 :: vv,hh @@ -3047,19 +3044,6 @@ SUBROUTINE MODIS_snow_alb (nx,ny,gfiler) close (10,status='keep') - ! Read tile-id raster file - allocate(tile_id(1:nx,1:ny)) - - fname=trim(gfiler)//'.rst' - open (10,file=fname,status='old',action='read', & - form='unformatted',convert='little_endian') - - do j=1,ny - read(10)tile_id(:,j) - end do - - close (10,status='keep') - !------------ Get the information on snow albedo ----- ! ----------- The information on snow albedo is stored in 10x10deg 30-arcsec resolution files. ! ----------- Read in this information, then loop over the tiles to find a corresponding snow albedo. From e96a4ce41d6d1e4a213fb1e2ba33cb27542f4425 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 10 Nov 2022 17:40:47 -0500 Subject: [PATCH 31/43] fixed call to subroutine MODIS_snow_alb() [see previous commit] --- .../GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index b88cd2fa2..e163fb2d5 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -661,7 +661,7 @@ PROGRAM mkCatchParam tmpstring = 'Step 14: Static snow albedo from MODIS' write (log_file,'(a)') trim(tmpstring) write (log_file,'(a)')' Creating file...' - call MODIS_snow_alb (nc,nr,gridnamer) + call MODIS_snow_alb ( ) write (log_file,'(a)')' Done.' write (log_file,'(a)')' ' endif From 7bff94341bdd0cf7138c1cf1cc52d829918add07 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Sun, 13 Nov 2022 21:22:09 -0500 Subject: [PATCH 32/43] updated to use just filled raster files --- .../Utils/Raster/mod_process_hres_data.F90 | 175 ++++++------------ 1 file changed, 58 insertions(+), 117 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index bea879c9f..eba91c7f5 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -3003,19 +3003,16 @@ SUBROUTINE MODIS_snow_alb ( ) character*200 :: fname character*2 :: vv,hh - integer :: n,maxcat,j,ncid,status + integer :: n,maxcat,ncid,status real,allocatable,dimension(:) :: min_lon,max_lon,min_lat,max_lat,snw_alb integer(kind=4),parameter :: xdim = 1200, ydim = 1200 - real,parameter :: alb_res=10.0/1200.0 real,dimension(xdim,ydim) :: stch_snw_alb_tmp real,dimension(36,18,xdim,ydim) :: stch_snw_alb - real :: minlon,maxlon,minlat,maxlat,pad_lon,pad_lat - real :: sno_alb_cnt,sno_alb_sum,sno_alb_cnt2,sno_alb_sum2 + real :: minlon,maxlon,minlat,maxlat + real :: sno_alb_cnt,sno_alb_sum integer :: vvtil_min,hhtil_min,vvtil_max,hhtil_max,hhtil,vvtil integer :: tindex1,pfaf1 - integer(kind=4) :: dummy,varid1 - integer(kind=4) :: imin,imax,jmin,jmax - integer(kind=4) :: imin2,imax2,jmin2,jmax2,count_init_invalid + integer(kind=4) :: imin,imax,jmin,jmax,varid1 logical :: file_exists ! Read number of catchment-tiles (maxcat) from catchment.def file @@ -3077,122 +3074,67 @@ SUBROUTINE MODIS_snow_alb ( ) enddo enddo + + if (minval(stch_snw_alb) .le. 0.0 .or. maxval(stch_snw_alb) .gt. 1.0) then + print*, 'There is a problem with snow albedo raster file. Non-physical values present. STOP!' + stop + endif ! loop over tiles - print*, 'Starting tile loop for snow albedo. ' - count_init_invalid=0 ! counter for non-valid snow albedo values (informational use only) + print*, 'Starting tile loop for snow albedo.' do n = 1, maxcat ! loop over tiles - - ! Set sums and counts to zero - sno_alb_sum =0. - sno_alb_cnt =0. - sno_alb_sum2=0. - sno_alb_cnt2=0. - - ! Use tile's min/max lat/lon info to identify the 10x10deg input file(s) - ! and read in snow albedo value(s). The "ceiling" and "floor" implements the "halo". - vvtil_min= floor((min_lat(n)+ 90.0)/10.)+1 - hhtil_min= floor((min_lon(n)+180.0)/10.)+1 - - ! if tile crosses the edge of the snow albedo 10x10deg box, expand the - ! search area into the neighbouring 10x10deg box - hhtil_max=hhtil_min - vvtil_max=vvtil_min - if (floor(min_lon(n)/10) .ne. floor(max_lon(n)/10)) hhtil_max=hhtil_min+1 - if (floor(min_lat(n)/10) .ne. floor(max_lat(n)/10)) vvtil_max=vvtil_min+1 - - ! Safety checks: - ! 1. Make sure vv's and hh's are within the range. If min>max, swap them. - if (vvtil_min .gt. vvtil_max) then - dummy =vvtil_min - vvtil_min=vvtil_max - vvtil_max=dummy - endif - if (hhtil_min .gt. hhtil_max) then - dummy =hhtil_min - hhtil_min=hhtil_max - hhtil_max=dummy - endif - - ! 2. Keep within the range. - vvtil_min=max(vvtil_min,1) - vvtil_max=min(vvtil_max,18) - hhtil_min=max(hhtil_min,1) - hhtil_max=min(hhtil_max,36) - - do hhtil=hhtil_min,hhtil_max ! loop through input files - horizontal direction - do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction - - ! Find indices ranges corresponding to the current tile area. - imin=floor((min_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) +1 - imax=floor((max_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) +1 - jmin=floor((min_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) +1 - jmax=floor((max_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) +1 - ! Keep within the range. - imin=max(imin,1) - imax=min(imax,xdim) - jmin=max(jmin,1) - jmax=min(jmax,ydim) - - ! Generate sums and counts using current tile corresponding indices - sno_alb_sum= sno_alb_sum + & - sum(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax), & - stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & - stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).le.1.0) - - sno_alb_cnt= sno_alb_cnt + & - count(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).gt.0.0 .and. & - stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax).le.1.0) - - end do ! vvtil - end do ! hhtil - - ! Calculate snow albedo for the current tile - snw_alb(n) = sno_alb_sum / max(1.0,sno_alb_cnt) - if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=MAPL_UNDEF !1.E15 - - ! If no valid solution found, and if tile size smaller than snow albedo resolution, - ! expand search area by 1-tile padding. - - ! Size of a tile (in both directions) - pad_lon=(max_lon(n)-min_lon(n)) - pad_lat=(max_lat(n)-min_lat(n)) - - if (snw_alb(n) .le. 0.0 .and. (pad_lon .lt. alb_res .or. pad_lat .lt. alb_res)) then - - count_init_invalid=count_init_invalid+1 + + ! Set sums and counts to zero + sno_alb_sum=0. + sno_alb_cnt=0. + + ! Use tile's min/max lat/lon info to identify the 10x10deg input file(s) + ! indexes + vvtil_min=floor((min_lat(n)+ 90.0)/10.)+1 + hhtil_min=floor((min_lon(n)+180.0)/10.)+1 + + ! if tile crosses the edge of the snow albedo 10x10deg box, expand the + ! search area into the neighbouring 10x10deg box + hhtil_max=hhtil_min + vvtil_max=vvtil_min + if (floor(min_lon(n)/10) .ne. floor(max_lon(n)/10)) hhtil_max=hhtil_min+1 + if (floor(min_lat(n)/10) .ne. floor(max_lat(n)/10)) vvtil_max=vvtil_min+1 + + ! Safety check; keep within the range + vvtil_min=max(vvtil_min,1) + vvtil_max=min(vvtil_max,18) + hhtil_min=max(hhtil_min,1) + hhtil_max=min(hhtil_max,36) + + do hhtil=hhtil_min,hhtil_max ! loop through input files - horizontal direction + do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction - do hhtil=hhtil_min,hhtil_max ! loop through input files - horizontal direction - do vvtil=vvtil_min,vvtil_max ! loop through input files - vertical direction - - ! Repeat the steps for extracting snow albedo value - imin2=floor((min_lon(n)-pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0))+1 - imax2=floor((max_lon(n)+pad_lon+180.0 - (hhtil-1)*10.0) * (xdim/10.0))+1 - jmin2=floor((min_lat(n)-pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0))+1 - jmax2=floor((max_lat(n)+pad_lat+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0))+1 - imin2=max(imin2,1) - imax2=min(imax2,xdim) - jmin2=max(jmin2,1) - jmax2=min(jmax2,ydim) - - sno_alb_sum2= sno_alb_sum2 + & - sum(stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2), & - stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & - stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).le.1.0) - - sno_alb_cnt2= sno_alb_cnt2 + & - count(stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).gt.0.0 .and. & - stch_snw_alb(hhtil,vvtil,imin2:imax2,jmin2:jmax2).le.1.0) - - end do ! vvtil - end do ! hhtil + ! Find indices ranges corresponding to the current tile area. + imin=floor((min_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) +1 + imax=floor((max_lon(n)+180.0 - (hhtil-1)*10.0) * (xdim/10.0)) +1 + jmin=floor((min_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) +1 + jmax=floor((max_lat(n)+ 90.0 - (vvtil-1)*10.0) * (ydim/10.0)) +1 + + ! if no matching grids, go to the next vv/hh box + if (imin .gt. xdim .or. jmin .gt. ydim .or. imax .lt. 1 .or. jmax .lt. 1) cycle + + ! Keep within the range, to include only the portion of the tile within this vv/hh box + imin=max(imin,1) + imax=min(imax,xdim) + jmin=max(jmin,1) + jmax=min(jmax,ydim) - snw_alb(n) = sno_alb_sum2 / max(1.0,sno_alb_cnt2) - if (snw_alb(n) .le. 0.0 .or. snw_alb(n) .gt. 1.0 ) snw_alb(n)=MAPL_UNDEF !1.E15 + ! Generate sums and counts using current tile corresponding indices + sno_alb_sum = sno_alb_sum + sum(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax)) + sno_alb_cnt = sno_alb_cnt + (imax-imin+1)*(jmax-jmin+1) - endif - + end do ! vvtil + end do ! hhtil + + ! If matching grids found, calculate snow albedo for the current tile + if (sno_alb_sum .ne. 0.0 .and. sno_alb_cnt .ne. 0) snw_alb(n)=min(1.0,max(0.0,sno_alb_sum/sno_alb_cnt)) + end do ! n-loop over tiles ! write snow albedo into clsm/catch_params.nc4 @@ -3205,7 +3147,6 @@ SUBROUTINE MODIS_snow_alb ( ) endif print*, 'Ended tile loop for snow albedo. ' - print*, 'Initially found ', count_init_invalid, ' tiles with no-data values for snow albedo (out of ', maxcat,')' END SUBROUTINE MODIS_snow_alb From 3befabea8ef35f0584f3f00b3815c77df96d74bd Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Tue, 15 Nov 2022 14:49:56 -0500 Subject: [PATCH 33/43] additional minor clean-up of subroutine MODIS_snow_alb() --- .../Utils/Raster/mod_process_hres_data.F90 | 40 ++++++++++--------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 index eba91c7f5..4909efa06 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mod_process_hres_data.F90 @@ -2995,7 +2995,11 @@ END SUBROUTINE hres_gswp2 SUBROUTINE MODIS_snow_alb ( ) - ! Process static snow albedo calculated from MODIS climatology and write into clsm/catch_params.nc4 + ! Map static, MODIS climatology-based snow albedo from 30-arcsec raster + ! grid to tile space and write into clsm/catch_params.nc4. + ! + ! Assumes that input snow albedo on raster grid is backfilled + ! (i.e., does not contain no-data values). ! ! Biljana Orescanin July 2022, SSAI@NASA @@ -3003,7 +3007,7 @@ SUBROUTINE MODIS_snow_alb ( ) character*200 :: fname character*2 :: vv,hh - integer :: n,maxcat,ncid,status + integer :: n,N_tile,ncid,status real,allocatable,dimension(:) :: min_lon,max_lon,min_lat,max_lat,snw_alb integer(kind=4),parameter :: xdim = 1200, ydim = 1200 real,dimension(xdim,ydim) :: stch_snw_alb_tmp @@ -3015,23 +3019,23 @@ SUBROUTINE MODIS_snow_alb ( ) integer(kind=4) :: imin,imax,jmin,jmax,varid1 logical :: file_exists - ! Read number of catchment-tiles (maxcat) from catchment.def file + ! Read number of catchment-tiles (N_tile) from catchment.def file fname='clsm/catchment.def' open (10,file=fname,status='old',action='read',form='formatted') - read(10,*) maxcat + read(10,*) N_tile ! Read min/max lat/lons to use when locating snow albedo grids in ! the stitched MODIS albedo file - allocate (min_lon(1:maxcat)) - allocate (min_lat(1:maxcat)) - allocate (max_lon(1:maxcat)) - allocate (max_lat(1:maxcat)) - allocate (snw_alb(1:maxcat)) + allocate (min_lon(1:N_tile)) + allocate (min_lat(1:N_tile)) + allocate (max_lon(1:N_tile)) + allocate (max_lat(1:N_tile)) + allocate (snw_alb(1:N_tile)) ! Start by setting all snow albedo values to missing snw_alb(:)=MAPL_UNDEF - do n = 1, maxcat + do n = 1, N_tile read (10,*) tindex1,pfaf1,minlon,maxlon,minlat,maxlat min_lon(n) = minlon max_lon(n) = maxlon @@ -3041,7 +3045,7 @@ SUBROUTINE MODIS_snow_alb ( ) close (10,status='keep') - !------------ Get the information on snow albedo ----- + ! ----------- Get the information on snow albedo ----- ! ----------- The information on snow albedo is stored in 10x10deg 30-arcsec resolution files. ! ----------- Read in this information, then loop over the tiles to find a corresponding snow albedo. @@ -3057,9 +3061,6 @@ SUBROUTINE MODIS_snow_alb ( ) ! sheets and is weighted by the grid-cell area). fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v2/snow_alb_FillVal_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' - ! MODIS-based climatology albedo raster files (NoData set to 1e+15) - !fname = '/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/albedo/snow/MODIS/v1/snow_alb_noFill_MOD10A1.061_30arcsec_H'//hh//'V'//vv//'.nc' - ! Open the file. (NF90_NOWRITE ensures read-only access to the file) status=NF_OPEN(trim(fname),NF_NOWRITE, ncid) ; VERIFY_(STATUS) ! Based on vars name, get the varids. @@ -3083,7 +3084,7 @@ SUBROUTINE MODIS_snow_alb ( ) ! loop over tiles print*, 'Starting tile loop for snow albedo.' - do n = 1, maxcat ! loop over tiles + do n = 1, N_tile ! loop over tiles ! Set sums and counts to zero sno_alb_sum=0. @@ -3126,14 +3127,15 @@ SUBROUTINE MODIS_snow_alb ( ) jmax=min(jmax,ydim) ! Generate sums and counts using current tile corresponding indices - sno_alb_sum = sno_alb_sum + sum(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax)) + sno_alb_sum = sno_alb_sum + sum(stch_snw_alb(hhtil,vvtil,imin:imax,jmin:jmax)) sno_alb_cnt = sno_alb_cnt + (imax-imin+1)*(jmax-jmin+1) end do ! vvtil end do ! hhtil - ! If matching grids found, calculate snow albedo for the current tile - if (sno_alb_sum .ne. 0.0 .and. sno_alb_cnt .ne. 0) snw_alb(n)=min(1.0,max(0.0,sno_alb_sum/sno_alb_cnt)) + ! If matching grids found, calculate snow albedo for the current tile; + ! ensure that resulting value is within physical range [0,1]. + if (sno_alb_cnt .ne. 0) snw_alb(n)=min(1.0,max(0.0,sno_alb_sum/sno_alb_cnt)) end do ! n-loop over tiles @@ -3142,7 +3144,7 @@ SUBROUTINE MODIS_snow_alb ( ) if(file_exists) then status = NF_OPEN ('clsm/catch_params.nc4', NF_WRITE, ncid ) ; VERIFY_(STATUS) - status = NF_PUT_VARA_REAL(NCID,NC_VarID(NCID,'SNOWALB'),(/1/),(/maxcat/),real(snw_alb)) ; VERIFY_(STATUS) + status = NF_PUT_VARA_REAL(NCID,NC_VarID(NCID,'SNOWALB'),(/1/),(/N_tile/),real(snw_alb)) ; VERIFY_(STATUS) STATUS = NF_CLOSE (NCID) ; VERIFY_(STATUS) endif From 3b15cac801fdaadebbacc25cf33bc0399025ed0a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Tue, 15 Nov 2022 15:41:59 -0500 Subject: [PATCH 34/43] added "only" specs to "use" statements in CatchmentRst.F90 and CatchmentCNRst.F90 --- .../GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 | 3 ++- .../GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 | 4 +++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 index ecab0906e..d28137c71 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentCNRst.F90 @@ -1,7 +1,8 @@ #include "MAPL_Generic.h" module CatchmentCNRstMod - use mk_restarts_getidsMod + use mk_restarts_getidsMod, ONLY: & + GetIds use mpi use MAPL use CatchmentRstMod, only : CatchmentRst diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index 34406b852..f0ac9b971 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -1,7 +1,9 @@ #include "MAPL_Generic.h" module CatchmentRstMod - use mk_restarts_getidsMod + use mk_restarts_getidsMod, ONLY: & + GetIds, & + ReadTileFile_RealLatLon use MAPL use MAPL_Base, ONLY: MAPL_UNDEF use mpi From 5aea1621481c430ed3b7a1f5a06020deaf06b0c0 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Wed, 16 Nov 2022 15:40:34 -0500 Subject: [PATCH 35/43] no re-tile for bcs parameters in catchment restart file --- .../Utils/mk_restarts/CatchmentRst.F90 | 208 +++--------------- 1 file changed, 35 insertions(+), 173 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index f0ac9b971..8e2d04b47 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -266,40 +266,14 @@ subroutine read_shared_nc4(this, formatter, rc) integer, optional, intent(out):: rc integer :: status - call MAPL_VarRead(formatter,"BF1",this%bf1, __RC__) - call MAPL_VarRead(formatter,"BF2",this%bf2, __RC__) - call MAPL_VarRead(formatter,"BF3",this%bf3, __RC__) + ! the four variables are used as scaler call MAPL_VarRead(formatter,"VGWMAX",this%vgwmax, __RC__) call MAPL_VarRead(formatter,"CDCR1",this%cdcr1, __RC__) call MAPL_VarRead(formatter,"CDCR2",this%cdcr2, __RC__) - call MAPL_VarRead(formatter,"PSIS",this%psis, __RC__) - call MAPL_VarRead(formatter,"BEE",this%bee, __RC__) call MAPL_VarRead(formatter,"POROS",this%poros, __RC__) - call MAPL_VarRead(formatter,"WPWET",this%wpwet, __RC__) - call MAPL_VarRead(formatter,"COND",this%cond, __RC__) - call MAPL_VarRead(formatter,"GNU",this%gnu, __RC__) - call MAPL_VarRead(formatter,"ARS1",this%ars1, __RC__) - call MAPL_VarRead(formatter,"ARS2",this%ars2, __RC__) - call MAPL_VarRead(formatter,"ARS3",this%ars3, __RC__) - call MAPL_VarRead(formatter,"ARA1",this%ara1, __RC__) - call MAPL_VarRead(formatter,"ARA2",this%ara2, __RC__) - call MAPL_VarRead(formatter,"ARA3",this%ara3, __RC__) - call MAPL_VarRead(formatter,"ARA4",this%ara4, __RC__) - call MAPL_VarRead(formatter,"ARW1",this%arw1, __RC__) - call MAPL_VarRead(formatter,"ARW2",this%arw2, __RC__) - call MAPL_VarRead(formatter,"ARW3",this%arw3, __RC__) - call MAPL_VarRead(formatter,"ARW4",this%arw4, __RC__) - call MAPL_VarRead(formatter,"TSA1",this%tsa1, __RC__) - call MAPL_VarRead(formatter,"TSA2",this%tsa2, __RC__) - call MAPL_VarRead(formatter,"TSB1",this%tsb1, __RC__) - call MAPL_VarRead(formatter,"TSB2",this%tsb2, __RC__) - call MAPL_VarRead(formatter,"ATAU",this%atau, __RC__) - call MAPL_VarRead(formatter,"BTAU",this%btau, __RC__) + call MAPL_VarRead(formatter,"TC",this%tc, __RC__) call MAPL_VarRead(formatter,"QC",this%qc, __RC__) -! -! call MAPL_VarRead(formatter,"OLD_ITY",this%ity, __RC__) -! call MAPL_VarRead(formatter,"CAPAC",this%capac, __RC__) call MAPL_VarRead(formatter,"CATDEF",this%catdef, __RC__) call MAPL_VarRead(formatter,"RZEXC",this%rzexc, __RC__) @@ -445,36 +419,11 @@ subroutine allocate_catch(this,rc) integer, optional, intent(out):: rc integer :: ntiles ntiles = this%ntiles - allocate( this% bf1(ntiles) ) - allocate( this% bf2(ntiles) ) - allocate( this% bf3(ntiles) ) allocate( this% vgwmax(ntiles) ) allocate( this% cdcr1(ntiles) ) allocate( this% cdcr2(ntiles) ) - allocate( this% psis(ntiles) ) - allocate( this% bee(ntiles) ) allocate( this% poros(ntiles) ) - allocate( this% wpwet(ntiles) ) - allocate( this% cond(ntiles) ) - allocate( this% gnu(ntiles) ) - allocate( this% ars1(ntiles) ) - allocate( this% ars2(ntiles) ) - allocate( this% ars3(ntiles) ) - allocate( this% ara1(ntiles) ) - allocate( this% ara2(ntiles) ) - allocate( this% ara3(ntiles) ) - allocate( this% ara4(ntiles) ) - allocate( this% arw1(ntiles) ) - allocate( this% arw2(ntiles) ) - allocate( this% arw3(ntiles) ) - allocate( this% arw4(ntiles) ) - allocate( this% tsa1(ntiles) ) - allocate( this% tsa2(ntiles) ) - allocate( this% tsb1(ntiles) ) - allocate( this% tsb2(ntiles) ) - allocate( this% atau(ntiles) ) - allocate( this% btau(ntiles) ) - allocate( this% ity(ntiles) ) + allocate( this% tc(ntiles,4) ) allocate( this% qc(ntiles,4) ) allocate( this% capac(ntiles) ) @@ -491,9 +440,6 @@ subroutine allocate_catch(this,rc) if (this%meta%has_variable('TSURF')) then allocate( this% tsurf(ntiles) ) endif - if (this%meta%has_variable('SNOWALB')) then - allocate( this% snowalb(ntiles) ) - endif allocate( this% wesnn1(ntiles) ) allocate( this% wesnn2(ntiles) ) @@ -553,38 +499,35 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) inquire(file = trim(DataDir)//'/clsm/catch_params.nc4', exist=file_exists) inquire(file = trim(DataDir)//"/clsm/CLM_veg_typs_fracs",exist=NewLand ) - if (size(this%ara1) /= this%ntiles ) then - ! it is just re-allocate - this%ity = DP2BR - this%ARA1 = DP2BR - this%ARA2 = DP2BR - this%ARA3 = DP2BR - this%ARA4 = DP2BR - this%ARS1 = DP2BR - this%ARS2 = DP2BR - this%ARS3 = DP2BR - this%ARW1 = DP2BR - this%ARW2 = DP2BR - this%ARW3 = DP2BR - this%ARW4 = DP2BR - this%ATAU = DP2BR - this%BTAU = DP2BR - this%PSIS = DP2BR - this%BEE = DP2BR - this%BF1 = DP2BR - this%BF2 = DP2BR - this%BF3 = DP2BR - this%TSA1 = DP2BR - this%TSA2 = DP2BR - this%TSB1 = DP2BR - this%TSB2 = DP2BR - this%COND = DP2BR - this%WPWET = DP2BR - this%POROS = DP2BR - this%VGWMAX = DP2BR - this%cdcr1 = DP2BR - this%cdcr2 = DP2BR - endif + this%ity = DP2BR + this%ARA1 = DP2BR + this%ARA2 = DP2BR + this%ARA3 = DP2BR + this%ARA4 = DP2BR + this%ARS1 = DP2BR + this%ARS2 = DP2BR + this%ARS3 = DP2BR + this%ARW1 = DP2BR + this%ARW2 = DP2BR + this%ARW3 = DP2BR + this%ARW4 = DP2BR + this%ATAU = DP2BR + this%BTAU = DP2BR + this%PSIS = DP2BR + this%BEE = DP2BR + this%BF1 = DP2BR + this%BF2 = DP2BR + this%BF3 = DP2BR + this%TSA1 = DP2BR + this%TSA2 = DP2BR + this%TSB1 = DP2BR + this%TSB2 = DP2BR + this%COND = DP2BR + this%WPWET = DP2BR + this%POROS = DP2BR + this%VGWMAX = DP2BR + this%cdcr1 = DP2BR + this%cdcr2 = DP2BR if(file_exists) then call CatchFmt%Open(trim(DataDir)//'/clsm/catch_params.nc4', pFIO_READ, __RC__) @@ -935,92 +878,15 @@ subroutine re_tile(this, InTileFile, OutBcsDir, OutTileFile, surflay, rc) var_out = this%poros(id_glb(:)) this%poros = var_out - var_out = this%cond(id_glb(:)) - this%cond = var_out - - var_out = this%psis(id_glb(:)) - this%psis = var_out - - var_out = this%bee(id_glb(:)) - this%bee = var_out - - var_out = this%wpwet(id_glb(:)) - this%wpwet = var_out - - var_out = this%gnu(id_glb(:)) - this%gnu = var_out - var_out = this%vgwmax(id_glb(:)) this%vgwmax = var_out - var_out = this%bf1(id_glb(:)) - this%bf1 = var_out - - var_out = this%bf2(id_glb(:)) - this%bf2 = var_out - - var_out = this%bf3(id_glb(:)) - this%bf3 = var_out - var_out = this%cdcr1(id_glb(:)) this%cdcr1 = var_out var_out = this%cdcr2(id_glb(:)) this%cdcr2 = var_out - var_out = this%ars1(id_glb(:)) - this%ars1 = var_out - - var_out = this%ars2(id_glb(:)) - this%ars2 = var_out - - var_out = this%ars3(id_glb(:)) - this%ars3 = var_out - - var_out = this%ara1(id_glb(:)) - this%ara1 = var_out - - var_out = this%ara2(id_glb(:)) - this%ara2 = var_out - - var_out = this%ara3(id_glb(:)) - this%ara3 = var_out - - var_out = this%ara4(id_glb(:)) - this%ara4 = var_out - - var_out = this%arw1(id_glb(:)) - this%arw1 = var_out - - var_out = this%arw2(id_glb(:)) - this%arw2 = var_out - - var_out = this%arw3(id_glb(:)) - this%arw3 = var_out - - var_out = this%arw4(id_glb(:)) - this%arw4 = var_out - - var_out = this%tsa1(id_glb(:)) - this%tsa1 = var_out - - var_out = this%tsa2(id_glb(:)) - this%tsa2 = var_out - - var_out = this%tsb1(id_glb(:)) - this%tsb1 = var_out - - var_out = this%tsb2(id_glb(:)) - this%tsb2 = var_out - - var_out = this%atau(id_glb(:)) - this%atau = var_out - - var_out = this%btau(id_glb(:)) - this%btau = var_out - - this%ity = [(k*1., k=1, out_ntiles)] - tmp2d = this%tc deallocate(this%tc) allocate(this%tc(out_ntiles, 4)) @@ -1098,11 +964,6 @@ subroutine re_tile(this, InTileFile, OutBcsDir, OutTileFile, surflay, rc) this%tsurf = var_out endif - if (this%meta%has_variable('SNOWALB')) then - var_out = this%snowalb(id_glb(:)) - this%snowalb = var_out - endif - ! CH CM CQ FR WW ! WW if(allocated(tmp2d)) deallocate(tmp2d) @@ -1298,12 +1159,13 @@ end subroutine re_scale subroutine set_scale_var(this ) class(CatchmentRst), intent(inout) :: this - this%old_catdef = this%catdef this%old_poros = this%poros this%old_cdcr1 = this%cdcr1 this%old_cdcr2 = this%cdcr2 - this%old_rzexc = this%rzexc this%old_vgwmax = this%vgwmax + + this%old_catdef = this%catdef + this%old_rzexc = this%rzexc this%old_sndzn3 = this%sndzn3 this%old_ghtcnt1 = this%ghtcnt1 this%old_ghtcnt2 = this%ghtcnt2 From 8020e32b0cabee822085555a6e50c2eb95f00c19 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Thu, 17 Nov 2022 08:51:35 -0500 Subject: [PATCH 36/43] add snowalb to restart when using regrid.pl --- .../Utils/mk_restarts/mk_CatchRestarts.F90 | 45 ++++++++++++++----- 1 file changed, 34 insertions(+), 11 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 index 4e6b2d94f..14e15320b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 @@ -267,7 +267,8 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) real, allocatable :: ARW1(:), ARW2(:), ARW3(:), ARW4(:) real, allocatable :: TSA1(:), TSA2(:), TSB1(:), TSB2(:) real, allocatable :: ATAU2(:), BTAU2(:), DP2BR(:), rity(:) - + real, allocatable :: snowalb(:) + real :: zdep1, zdep2, zdep3, zmet, term1, term2, rdum real, allocatable :: var1(:),var2(:,:) character*256 :: vname @@ -280,9 +281,10 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) integer, pointer :: Ido(:), idx(:), id(:) logical :: InIsOld type(NetCDF4_Fileformatter) :: InFmt,OutFmt,CatchFmt - type(FileMetadata) :: InCfg,OutCfg + type(FileMetadata) :: InCfg,OutCfg, meta_ type(StringVariableMap), pointer :: variables type(Variable), pointer :: myVariable + type(Variable) :: var type(StringVariableMapIterator) :: var_iter character(len=:), pointer :: var_name,dname type(StringVector), pointer :: var_dimensions @@ -303,15 +305,6 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) call InFmt%open(InRestart,pFIO_READ,__RC__) InCfg=InFmt%read(__RC__) - call MAPL_IOChangeRes(InCfg,OutCfg,(/'tile'/),(/ntiles/),__RC__) - i = index(InRestart,'/',back=.true.) - OutFileName = "OutData/"//trim(InRestart(i+1:)) - call OutFmt%create(OutFileName,__RC__) - call OutFmt%write(OutCfg,__RC__) - call MAPL_IOCountNonDimVars(OutCfg,nvars,__RC__) - - allocate(written(nvars)) - written=.false. else @@ -411,6 +404,13 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) call MAPL_VarRead ( catchFmt ,'WPWET', WPWET, __RC__) call MAPL_VarRead ( catchFmt ,'DP2BR', DP2BR, __RC__) call MAPL_VarRead ( catchFmt ,'POROS', POROS, __RC__) + + meta_ = CatchFmt%read(__RC__) + if (meta_%has_variable('SNOWALB')) then + allocate(snowalb(ntiles)) + call MAPL_VarRead ( CatchFmt ,'SNOWALB', snowalb, __RC__) + endif + call catchFmt%close(__RC__) else @@ -571,6 +571,23 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) endif HAVE if (filetype == 0) then + + call MAPL_IOChangeRes(InCfg,OutCfg,(/'tile'/),(/ntiles/),__RC__) + i = index(InRestart,'/',back=.true.) + OutFileName = "OutData/"//trim(InRestart(i+1:)) + call OutFmt%create(OutFileName,__RC__) + + if (allocated(snowalb)) then + var = Variable(type=pFIO_REAL32, dimensions='tile') + call var%add_attribute('long_name', 'snow_albedo') + call var%add_attribute('units', '1') + call OutCFG%add_variable('SNOWALB', var) + endif + + call OutFmt%write(OutCfg,__RC__) + + + call MAPL_VarWrite(OutFmt,names(1),BF1(Idx)) call MAPL_VarWrite(OutFmt,names(2),BF2(Idx)) call MAPL_VarWrite(OutFmt,names(3),BF3(Idx)) @@ -602,8 +619,13 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) call MAPL_VarWrite(OutFmt,names(29),BTAU2(Idx)) call MAPL_VarWrite(OutFmt,'OLD_ITY',rity(Idx)) + if (allocated(snowalb)) then + call MAPL_VarWrite(OutFmt,'SNOWALB',snowalb(Idx)) + endif call MAPL_IOCountNonDimVars(InCfg,nvars) + allocate(written(nvars)) + written=.false. variables => InCfg%get_variables() var_iter = variables%begin() @@ -616,6 +638,7 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) if ( trim(var_name) == trim(names(j)) ) written(i) = .true. enddo if (trim(var_name) == "OLD_ITY" ) written(i) = .true. + if (trim(var_name) == "SNOWALB" ) written(i) = .true. call var_iter%next() From bceedb3b4285624eceb3d7593523f34fde139677 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Thu, 17 Nov 2022 14:02:29 -0500 Subject: [PATCH 37/43] minor clarification in comments (CatchmentRst.F90) --- .../GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index 8e2d04b47..ae80ddd03 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -266,12 +266,13 @@ subroutine read_shared_nc4(this, formatter, rc) integer, optional, intent(out):: rc integer :: status - ! the four variables are used as scaler + ! these four (time-invariant) variables are used for rescaling of prognostic variables call MAPL_VarRead(formatter,"VGWMAX",this%vgwmax, __RC__) call MAPL_VarRead(formatter,"CDCR1",this%cdcr1, __RC__) call MAPL_VarRead(formatter,"CDCR2",this%cdcr2, __RC__) call MAPL_VarRead(formatter,"POROS",this%poros, __RC__) + ! Catchment model prognostic variables (and some diagnostics needed in Catch restart for GCM) call MAPL_VarRead(formatter,"TC",this%tc, __RC__) call MAPL_VarRead(formatter,"QC",this%qc, __RC__) call MAPL_VarRead(formatter,"CAPAC",this%capac, __RC__) @@ -1146,7 +1147,7 @@ subroutine re_scale(this, surflay, wemin_in, wemin_out, rc) endif - ! PEATCLSM - ensure low CATDEF on peat tiles where "old" restart is not also peat + ! PEATCLSM - ensure low CATDEF on peat tiles where "old" restart is not also peat ! ------------------------------------------------------------------------------- where ( (this%old_poros < PEATCLSM_POROS_THRESHOLD) .and. (this%poros >= PEATCLSM_POROS_THRESHOLD) ) From 91ac79a1f1f138bf3428f619ae787cb62f6f68f8 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Fri, 18 Nov 2022 10:10:23 -0500 Subject: [PATCH 38/43] remove snowalb when target bc does not have --- .../Utils/mk_restarts/CatchmentRst.F90 | 3 ++- .../Utils/mk_restarts/mk_CatchRestarts.F90 | 10 ++++++---- 2 files changed, 8 insertions(+), 5 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index ae80ddd03..fdf99f642 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -581,8 +581,9 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) call var%add_attribute('units', '1') call this%meta%add_variable('SNOWALB', var) endif + elseif (this%meta%has_variable('SNOWALB')) then + call this%meta%remove_variable('SNOWALB') endif - call CatchFmt%close() else open(unit=21, file=trim(DataDir)//'/clsm/mosaic_veg_typs_fracs',form='formatted') diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 index 14e15320b..60cac1687 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 @@ -578,10 +578,12 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) call OutFmt%create(OutFileName,__RC__) if (allocated(snowalb)) then - var = Variable(type=pFIO_REAL32, dimensions='tile') - call var%add_attribute('long_name', 'snow_albedo') - call var%add_attribute('units', '1') - call OutCFG%add_variable('SNOWALB', var) + if (.not. OutCFG%has_varibale('SNOWALB') then + var = Variable(type=pFIO_REAL32, dimensions='tile') + call var%add_attribute('long_name', 'snow_albedo') + call var%add_attribute('units', '1') + call OutCFG%add_variable('SNOWALB', var) + endif endif call OutFmt%write(OutCfg,__RC__) From d718c1f97b140016ab827811dcc182c0f948523d Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Fri, 18 Nov 2022 10:57:33 -0500 Subject: [PATCH 39/43] correct typo --- .../GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 index 60cac1687..9fe3f0815 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 @@ -578,7 +578,7 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) call OutFmt%create(OutFileName,__RC__) if (allocated(snowalb)) then - if (.not. OutCFG%has_varibale('SNOWALB') then + if (.not. OutCFG%has_varibale('SNOWALB')) then var = Variable(type=pFIO_REAL32, dimensions='tile') call var%add_attribute('long_name', 'snow_albedo') call var%add_attribute('units', '1') From f0efe1ee6bcda7e19495939d6478eb25c19a14ef Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Fri, 18 Nov 2022 11:35:48 -0500 Subject: [PATCH 40/43] fix typo again --- .../GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 index 9fe3f0815..1c0df8c92 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 @@ -578,7 +578,7 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) call OutFmt%create(OutFileName,__RC__) if (allocated(snowalb)) then - if (.not. OutCFG%has_varibale('SNOWALB')) then + if (.not. OutCFG%has_variable('SNOWALB')) then var = Variable(type=pFIO_REAL32, dimensions='tile') call var%add_attribute('long_name', 'snow_albedo') call var%add_attribute('units', '1') From e375bdfc9d08688180e731f4997a390d9b63d4a7 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Fri, 18 Nov 2022 15:12:43 -0500 Subject: [PATCH 41/43] allocate GNU . bug fix --- .../GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index fdf99f642..36b994176 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -523,6 +523,7 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) this%TSA2 = DP2BR this%TSB1 = DP2BR this%TSB2 = DP2BR + this%GNU = DP2BR this%COND = DP2BR this%WPWET = DP2BR this%POROS = DP2BR From 2fe381e04cd6d647d324939be6c4955f26b4f2c2 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Mon, 21 Nov 2022 10:11:49 -0500 Subject: [PATCH 42/43] remove snowalb when there is not catch_params.nc4 --- .../GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 index 36b994176..477e17d08 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/CatchmentRst.F90 @@ -624,6 +624,9 @@ subroutine add_bcs_to_rst(this, surflay, DataDir, rc) CLOSE (24, STATUS = 'KEEP') CLOSE (25, STATUS = 'KEEP') CLOSE (26, STATUS = 'KEEP') + + if (this%meta%has_variable('SNOWALB')) call this%meta%remove_variable('SNOWALB') + endif do n=1,ntiles From 84a3caff3edb0281e4d96907e20be1f711abb2da Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Tue, 22 Nov 2022 11:10:31 -0500 Subject: [PATCH 43/43] restore mk_CatchRestarts so regrid.pl would not work with new bcs with snowalb --- .../Utils/mk_restarts/mk_CatchRestarts.F90 | 45 +++++-------------- 1 file changed, 11 insertions(+), 34 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 index 14e15320b..4e6b2d94f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 @@ -267,8 +267,7 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) real, allocatable :: ARW1(:), ARW2(:), ARW3(:), ARW4(:) real, allocatable :: TSA1(:), TSA2(:), TSB1(:), TSB2(:) real, allocatable :: ATAU2(:), BTAU2(:), DP2BR(:), rity(:) - real, allocatable :: snowalb(:) - + real :: zdep1, zdep2, zdep3, zmet, term1, term2, rdum real, allocatable :: var1(:),var2(:,:) character*256 :: vname @@ -281,10 +280,9 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) integer, pointer :: Ido(:), idx(:), id(:) logical :: InIsOld type(NetCDF4_Fileformatter) :: InFmt,OutFmt,CatchFmt - type(FileMetadata) :: InCfg,OutCfg, meta_ + type(FileMetadata) :: InCfg,OutCfg type(StringVariableMap), pointer :: variables type(Variable), pointer :: myVariable - type(Variable) :: var type(StringVariableMapIterator) :: var_iter character(len=:), pointer :: var_name,dname type(StringVector), pointer :: var_dimensions @@ -305,6 +303,15 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) call InFmt%open(InRestart,pFIO_READ,__RC__) InCfg=InFmt%read(__RC__) + call MAPL_IOChangeRes(InCfg,OutCfg,(/'tile'/),(/ntiles/),__RC__) + i = index(InRestart,'/',back=.true.) + OutFileName = "OutData/"//trim(InRestart(i+1:)) + call OutFmt%create(OutFileName,__RC__) + call OutFmt%write(OutCfg,__RC__) + call MAPL_IOCountNonDimVars(OutCfg,nvars,__RC__) + + allocate(written(nvars)) + written=.false. else @@ -404,13 +411,6 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) call MAPL_VarRead ( catchFmt ,'WPWET', WPWET, __RC__) call MAPL_VarRead ( catchFmt ,'DP2BR', DP2BR, __RC__) call MAPL_VarRead ( catchFmt ,'POROS', POROS, __RC__) - - meta_ = CatchFmt%read(__RC__) - if (meta_%has_variable('SNOWALB')) then - allocate(snowalb(ntiles)) - call MAPL_VarRead ( CatchFmt ,'SNOWALB', snowalb, __RC__) - endif - call catchFmt%close(__RC__) else @@ -571,23 +571,6 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) endif HAVE if (filetype == 0) then - - call MAPL_IOChangeRes(InCfg,OutCfg,(/'tile'/),(/ntiles/),__RC__) - i = index(InRestart,'/',back=.true.) - OutFileName = "OutData/"//trim(InRestart(i+1:)) - call OutFmt%create(OutFileName,__RC__) - - if (allocated(snowalb)) then - var = Variable(type=pFIO_REAL32, dimensions='tile') - call var%add_attribute('long_name', 'snow_albedo') - call var%add_attribute('units', '1') - call OutCFG%add_variable('SNOWALB', var) - endif - - call OutFmt%write(OutCfg,__RC__) - - - call MAPL_VarWrite(OutFmt,names(1),BF1(Idx)) call MAPL_VarWrite(OutFmt,names(2),BF2(Idx)) call MAPL_VarWrite(OutFmt,names(3),BF3(Idx)) @@ -619,13 +602,8 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) call MAPL_VarWrite(OutFmt,names(29),BTAU2(Idx)) call MAPL_VarWrite(OutFmt,'OLD_ITY',rity(Idx)) - if (allocated(snowalb)) then - call MAPL_VarWrite(OutFmt,'SNOWALB',snowalb(Idx)) - endif call MAPL_IOCountNonDimVars(InCfg,nvars) - allocate(written(nvars)) - written=.false. variables => InCfg%get_variables() var_iter = variables%begin() @@ -638,7 +616,6 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) if ( trim(var_name) == trim(names(j)) ) written(i) = .true. enddo if (trim(var_name) == "OLD_ITY" ) written(i) = .true. - if (trim(var_name) == "SNOWALB" ) written(i) = .true. call var_iter%next()