!-----------------------BEGIN NOTICE -- DO NOT EDIT----------------------- ! NASA Goddard Space Flight Center Land Data Toolkit (LDT) v1.0 !-------------------------END NOTICE -- DO NOT EDIT----------------------- #include "LDT_misc.h" module LDT_climateBCs !BOP ! ! !MODULE: LDT_climateBCs ! ! !DESCRIPTION: ! This module is good to compute climatological BCs on the LIS_domain at output_interval ! using any pre-processed climatological BCs data from any temporal and latlon spatial resolutions. ! SUBROUTINE read_climate_bcs is the only PUBLIC subroutine that has : ! 4 input arguments - nest (nest #). and 3 parameters read via LDT config: SOURCE. out_interval and projection ! 1 oputput argument - CLIMDATA (LDT_rc%lnc(nest),LDT_rc%lnr(nest), time_steps) where time_step is 12 if out_interval iss "monthly" ! ! The module has 3 main subroutines and they have their own subroutines: ! 1) SUBROUTINE annual_cycle that contains FUNCTION MidTime and FUNCTION TimeFrac ! 2) SUBROUTINE init_cp ! 3) SUBROUTINE read_data that contains SUBROUTINE read_10deg_geos5tiles; SUBROUTINE regrid_to_lisgrid; and SUBROUTINE read_global_nc4 ! !REVISION HISTORY: ! ! 27 August 2020: Sarith Mahanama; Initial implementation use ESMF use CLSM_util, ONLY: NC_VarID, & c_data => G5_BCSDIR use LDT_timeMgrMod, only : LDT_calendar use LDT_gridmappingMod use LDT_coreMod use netcdf use LDT_gfracMod,only : LDT_gfrac_struc implicit none PRIVATE public :: read_clim_bcs type :: climate_param_type integer :: NTIMES, NX, NY type(ESMF_TimeInterval) :: data_dtstep, out_dtstep integer, allocatable, dimension (:) :: data_doys character*20 :: var1, var2 real :: sf integer, dimension(2) :: range logical :: stitch character*300 :: BCS_PATH real :: param_gridDesc(20) end type climate_param_type contains SUBROUTINE read_clim_bcs (nest, SOURCE, out_interval, project, CLIMDATA) ! SUBROUTINE read_clim_bcs (nest, SOURCE, CLIMDATA, maskarray) implicit none integer, intent (in) :: nest real, pointer, dimension(:,:,:), intent(inout) :: CLIMDATA real, optional, intent(inout) :: maskarray(LDT_rc%lnc(nest),LDT_rc%lnr(nest)) type(climate_param_type) :: CP call init_cp (SOURCE, CP) if (out_interval == "daily" ) call ESMF_TimeIntervalSet(CP%out_dtstep, h=24 , rc=status ) ; VERIFY_(STATUS) if (out_interval == "8day" ) call ESMF_TimeIntervalSet(CP%out_dtstep, h=24*8, rc=status ) ; VERIFY_(STATUS) if (out_interval == "monthly") call ESMF_TimeIntervalSet(CP%out_dtstep,mm=1 , rc=status ) ; VERIFY_(STATUS) if(trim (project) == 'latlon') CP%param_gridDesc(1) = 0. call annual_cycle (nest, project, CP, CLIMDATA) END SUBROUTINE read_clim_bcs !################################################################################### subroutine annual_cycle (nest, project, CP, CLIMDATA) implicit none integer, intent (in) :: nest type(climate_param_type), intent(inout) :: CP real, pointer, dimension(:,:,:), intent(inout) :: CLIMDATA character(*), intent (in) :: project type(ESMF_Time) :: CURRENT_TIME, OUT_RING, END_TIME, TIME1, TIME2 type(ESMF_TimeInterval) :: DAY_DT integer :: t1, t2, switch,out_time, time_count real :: subparam_gridDesc(20) ! Input parameter grid desc fgrd integer :: glpnc, glpnr ! Parameter (global) total columns and rows integer :: subpnc, subpnr ! Parameter subsetted columns and rows integer, allocatable :: lat_line(:,:), lon_line(:,:) real, allocatable, dimension (:,:) :: var_subset1, var_subset2, out_ave integer, parameter :: ref_year = 2002 logical :: last real :: time_frac ! ------------------------------------------------------------------- ! PREPARE SUBSETTED PARAMETER GRID FOR READING IN BCS DATA ! ------------------------------------------------------------------- !- Map Parameter Grid Info to LIS Target Grid/Projection Info -- subparam_gridDesc = 0. call LDT_RunDomainPts( nest, project, CP%param_gridDesc(:), & glpnc, glpnr, subpnc, subpnr, subparam_gridDesc, lat_line, lon_line) if((subpnc /= size(CLIMDATA,1)) .or.(subpnr /= size(CLIMDATA,2))) then write(LDT_logunit,*)'[ERR] LIS domain dimensions mismatch in READ_CLIM_BCS' write(LDT_logunit,*)' input dimensions : ', size(CLIMDATA,1), size(CLIMDATA,2) write(LDT_logunit,*)' domain dimensions: ', subpnc, subpnr VERIFY_(1) call LDT_endrun endif ! initialize time loop ! -------------------- call ESMF_TimeSet (CURRENT_TIME, yy=ref_year, mm=1, dd=1, rc=status) ; VERIFY_(STATUS) call ESMF_TimeSet (END_TIME , yy=ref_year + 1, mm=1, dd=1, rc=status) ; VERIFY_(STATUS) call ESMF_TimeIntervalSet (DAY_DT, h=24, rc=status) ; VERIFY_(STATUS) OUT_RING = CURRENT_TIME t1 = CP%NTIMES t2 = 1 out_time = 1 time_count= 1 switch = 0 last = .false. TIME1 = MidTime (CP%data_doys(t1),CP%data_doys(1),start=.true.) TIME2 = MidTime (CP%data_doys(1) ,CP%data_doys(2)) allocate (var_subset1 (subpnc,subpnr)) allocate (var_subset2 (subpnc,subpnr)) allocate (out_ave (subpnc,subpnr)) call read_data (CP, var_subset1, t1) call read_data (CP, var_subset2, t2) out_ave = 0. ! Annual loop to compute climatology at out_interval on the lis grid - this loops through non-leap ref_year do while (CURRENT_TIME <= END_TIME) UPDATE_BCS: if (CURRENT_TIME >= TIME2) then ! update time infor if(t2 == 1) t2 = 2 t1 = t2 t2 = t2 + 1 if(t2 > CP%NTIMES) then t2 = 1 t1 = CP%NTIMES switch = 1 endif TIME1 = TIME2 if(last) then TIME2 = MidTime (data_doys(1),data_doys(2), last = last) t2 = 1 else TIME2 = MidTime (data_doys(t1),data_doys(t2)) endif if (switch == 1) last = .true. ! read BCs data var_subset1 = var_subset2 call read_data (CP, var_subset2,t2) endif UPDATE_BCS ! sum of daily values to compute CLIMDATA data at out_interval time_frac = TimeFrac (TIME1, CURRENT_TIME, TIME2) out_ave = out_ave + (var_subset1*time_frac + var_subset2* (1. - time_frac)) PROC_OUTPPUT: if (CURRENT_TIME == OUT_RING) then if (out_time > size( CLIMDATA,3)) then write(LDT_logunit,*)'[ERR] out_interval mismatch in READ_CLIM_BCS' write(LDT_logunit,*)' out_interval : ', size(CLIMDATA,3) VERIFY_(1) call LDT_endrun endif OUT_RING = CURRENT_TIME + CP%out_dtstep CLIMDATA(:,:,out_time) = out_ave / real(time_count) out_time = out_time + 1 time_count = 0 out_ave = 0. end if PROC_OUTPPUT time_count = time_count + 1 CURRENT_TIME = CURRENT_TIME + DAY_DT end do contains function MidTime (DOY1, DOY2, start, last) result (MT) implicit none integer, intent (in) :: DOY1, DOY2 type(ESMF_Time) :: TIME1, MT, CYB type(ESMF_TimeInterval) :: DAY_SHIFT integer :: status, yyyy logical, intent(in), optional :: start, last if (present(start)) then yyyy = ref_year -1 else yyyy = ref_year endif call ESMF_TimeSet (CYB, yy=yyyy, mm=1, dd=1, rc=status) ; VERIFY_(STATUS) if(present(last)) then yyyy = ref_year + 1 call ESMF_TimeSet (CYB, yy=yyyy, mm=1, dd=1, rc=status) ; VERIFY_(STATUS) endif call ESMF_TimeIntervalSet (DAY_SHIFT, h=24*(DOY1-1), rc=status ) ; VERIFY_(STATUS) TIME1 = CYB + DAY_SHIFT if (DOY2 > DOY1) then call ESMF_TimeIntervalSet (DAY_SHIFT, h=24*(DOY2-DOY1)/2, rc=status ); VERIFY_(STATUS) MT = TIME1 + DAY_SHIFT else ! climatologies don't account for leap years) call ESMF_TimeIntervalSet (DAY_SHIFT, h=24*(366-DOY1)/2, rc=status ) ; VERIFY_(STATUS) MT = TIME1 + DAY_SHIFT if (DOY2 > 1) then call ESMF_TimeIntervalSet (DAY_SHIFT, h=24*(DOY2-1)/2, rc=status ); VERIFY_(STATUS) MT = MT + DAY_SHIFT endif endif end function MidTime ! --------------------------------------------------------------------------- real function TimeFrac (PTIME, CTIME, NTIME) result (tf) implicit none type(ESMF_Time),INTENT(IN) :: PTIME, CTIME, NTIME real(ESMF_KIND_R8) :: tdiff_21, tdiff_20, tfrac type(ESMF_TimeInterval) :: T21, T20 integer :: status T21 = NTIME - PTIME T20 = NTIME - CTIME call ESMF_TimeIntervalGet(T21, s_r8=tdiff_21,rc=status) ; VERIFY_(STATUS) call ESMF_TimeIntervalGet(T20,s_r8=tdiff_20,rc=status) ; VERIFY_(STATUS) tfrac = tdiff_20/tdiff_21 tf = tfrac end function TimeFrac end subroutine annual_cycle !################################################################################### SUBROUTINE init_cp (SOURCE, CP) implicit none type(climate_param_type), intent(inout) :: CP character(*), intent (in) :: SOURCE integer, parameter :: N_MONTHLY_DATES = 12 integer, parameter :: N_MODIS_DATES = 46 integer, dimension (N_MODIS_DATES), parameter :: & MODIS_DOYS = (/ & 1 , 9, 17, 25, 33, 41, 49, 57, 65, & 73 , 81, 89, 97,105,113,121,129,137, & 145,153,161,169,177,185,193,201,209, & 217,225,233,241,249,257,265,273,281, & 289,297,305,313,321,329,337,345,353,361/) integer, dimension (N_MONTHLY_DATES), parameter :: & MONTHLY_DOYS = (/ & 1 , 32, 60, 91,121,152, & 182,213,244,274,305,335/) integer :: status integer :: input_cols integer :: input_rows real :: IN_xres real :: IN_yres real :: subparam_gridDesc(20) ! Input parameter grid desc fgrd select case (SOURCE) case ('GSWPH') CP%NTIMES = N_MONTHLY_DATES CP%NX = 43200 CP%NY = 21600 CP%var1 = 'grnFrac' CP%sf = 0.001 CP%range(1) = 0 CP%range(2) = 1000 CP%BCS_PATH = trim(c_data)//'/GSWP2_30sec_VegParam/GSWP2_VegParam_' CP%stitch = .TRUE. allocate (CP%data_doys (1: CP%NTIMES)) CP%data_doys= MONTHLY_DOYS call ESMF_TimeIntervalSet(CP%data_dtstep, mm=1, rc=status ) ; VERIFY_(STATUS) case ('MCD15A2H') CP%NTIMES = N_MODIS_DATES CP%NX = 86400 CP%NY = 43200 CP%var1 = 'LAI_500m' CP%sf = 1. CP%range(1) = 0 CP%range(2) = 10 CP%BCS_PATH = 'LS_PARAMETERS/MODIS/MCD15A2H.006/MCD15A2H.006_LAI_YYYY' CP%stitch = .FALSE. allocate (CP%data_doys (1: CP%NTIMES)) CP%data_doys= MODIS_DOYS call ESMF_TimeIntervalSet(CP%data_dtstep, h=24*8, rc=status ) ; VERIFY_(STATUS) case ('MCD43GFv5') CP%NTIMES = N_MODIS_DATES CP%NX = 43200 CP%NY = 21600 CP%var1 = 'Alb_0.3_0.7' CP%var2 = 'Alb_0.7_5.0' CP%sf = 0.001 CP%range(1) = 0 CP%range(2) = 1000 CP%stitch = .TRUE. CP%BCS_PATH = trim(c_data)//'/MODIS-Albedo2/MCD43GF_wsa_' allocate (CP%data_doys (1: CP%NTIMES)) CP%data_doys= MODIS_DOYS call ESMF_TimeIntervalSet(CP%data_dtstep, h=24*8, rc=status ) ; VERIFY_(STATUS) case ('GLASSA') CP%NTIMES = N_MODIS_DATES CP%NX = 7200 CP%NY = 3600 CP%var1 = 'LAI' CP%sf = 0.01 CP%range(1) = 0 CP%range(2) = 1000 CP%stitch = .FALSE. CP%BCS_PATH = trim(c_data)//'/GLASS-LAI/AVHRR.v4/GLASS01B02.V04.AYYYY' allocate (CP%data_doys (1: CP%NTIMES)) CP%data_doys= MODIS_DOYS call ESMF_TimeIntervalSet(CP%data_dtstep, h=24*8, rc=status ) ; VERIFY_(STATUS) case ('GLASSM') CP%NTIMES = N_MODIS_DATES CP%NX = 7200 CP%NY = 3600 CP%var1 = 'LAI' CP%sf = 0.01 CP%range(1) = 0 CP%range(2) = 1000 CP%stitch = .FALSE. CP%BCS_PATH = trim(c_data)//'/GLASS-LAI/MODIS.v4/GLASS01B01.V04.AYYYY' allocate (CP%data_doys (1: CP%NTIMES)) CP%data_doys = MODIS_DOYS call ESMF_TimeIntervalSet(CP%data_dtstep, h=24*8, rc=status ) ; VERIFY_(STATUS) case default print *, 'Unknown climatological data label : ', trim (SOURCE) VERIFY_(1) end select input_cols = CP%NX input_rows = CP%NY IN_xres = 360./REAL(CP%NX) IN_yres = 180./REAL(CP%NY) CP%param_gridDesc(1) = 0. ! Latlon CP%param_gridDesc(2) = input_cols CP%param_gridDesc(3) = input_rows CP%param_gridDesc(4) = -60.0 + (IN_yres/2) ! LL lat CP%param_gridDesc(5) = -180.0 + (IN_xres/2) ! LL lon CP%param_gridDesc(6) = 128 CP%param_gridDesc(7) = 90.0 - (IN_yres/2) ! UR lat CP%param_gridDesc(8) = 180.0 - (IN_xres/2) ! UR lon CP%param_gridDesc(9) = IN_yres CP%param_gridDesc(10) = IN_xres CP%param_gridDesc(20) = 64 END SUBROUTINE init_cp !################################################################################### subroutine read_data implicit none contains subroutine read_10deg_geos5tiles implicit none integer :: jx, ix, status, time_slice, ncid,iLL,jLL, nc_10, nr_10, d_undef character*2 :: hh,vv real :: sf character*300 :: fname character*20 :: lai_name integer, allocatable, dimension (:,:) :: net_data1 do jx = 1,18 do ix = 1,36 write (vv,'(i2.2)')jx write (hh,'(i2.2)')ix fname = trim(c_data)//trim(lai_name)//'lai_clim.H'//hh//'V'//vv//'.nc' status = NF_OPEN(trim(fname),NF_NOWRITE, ncid) if(status == 0) then status = NF_GET_att_INT (ncid,NF_GLOBAL,'i_ind_offset_LL',iLL); VERIFY_(STATUS) status = NF_GET_att_INT (ncid,NF_GLOBAL,'j_ind_offset_LL',jLL); VERIFY_(STATUS) status = NF_GET_att_INT (ncid,4,'UNDEF',d_undef); VERIFY_(STATUS) status = NF_GET_att_REAL (ncid,4,'ScaleFactor',sf); VERIFY_(STATUS) status = NF_GET_VARA_INT (ncid, 4,(/1,1,time_slice/),(/nc_10,nr_10,1/),net_data1); VERIFY_(STATUS) status = NF_CLOSE(ncid) endif end do end do end subroutine read_10deg_geos5tiles ! --------------------------------------------------------------------- SUBROUTINE regrid_to_lisgrid (subpnc, subpnr, lon_line, lat_line) implicit none ! arrays read from LL integer, parameter :: nodata = -999 integer :: mi ! Total number of input param grid array points integer :: mo ! Total number of output LIS grid array points integer :: nc, nr, i integer, allocatable :: n11(:) ! Array that maps the location of each input grid ! point in the output grid. real, allocatable :: gi(:) ! Input parameter 1d grid logical*1,allocatable :: li(:) ! Input logical mask (to match gi) real :: go1(LDT_rc%lnc(n)*LDT_rc%lnr(n)) ! Output lis 1d grid logical*1 :: lo1(LDT_rc%lnc(n)*LDT_rc%lnr(n)) ! Output logical mask (to match go) do nr = 1, subpnr do nc = 1, subnc var_in (nc,nr) = data_in (lon_line(nc,1),lat_line (1,nr)) end do end do ! ------------------------------------------------------------------- ! AGGREGATING FINE-SCALE GRIDS TO COARSER LIS OUTPUT GRID ! ------------------------------------------------------------------- mi = subpnc*subpnr mo = LDT_rc%lnc(n)*LDT_rc%lnr(n) #if 0 if( mi .ne. mo .and. LDT_rc%irrigfrac_gridtransform(n) == "none" ) then write(LDT_logunit,*) "[ERR] Spatial transform, 'none', is selected, but number of" write(LDT_logunit,*) " input and output points do not match. Select other spatial" write(LDT_logunit,*) " option (e.g., average, neighbor, etc.)." write(LDT_logunit,*) "Program stopping ..." call LDT_endrun endif #endif allocate( li(mi), gi(mi), n11(mi) ) gi = nodata li = .false. lo1 = .false. !- Assign 2-D array to 1-D for aggregation routines: i = 0 do nr = 1, subpnr do nc = 1, subpnc i = i + 1 if( var_in(nc,nr) >= 0. ) then gi(i) = var_in(nc,nr) li(i) = .true. endif enddo enddo !- Create mapping between parameter domain and LIS grid domain: call upscaleByAveraging_input( subparam_gridDesc, & LDT_rc%gridDesc(n,:), mi, mo, n11 ) !- Calculate total counts for inland water type in each coarse gridcell: ! call upscaleByCnt( mi, mo, num_bins, noncrop, n11, li, gi, & call upscaleByCnt( mi, mo, num_bins, LDT_rc%udef, n11, li, gi, & lo2, go2 ) irrigfrac_cnt = 0. i = 0 do nr = 1, LDT_rc%lnr(n) do nc = 1, LDT_rc%lnc(n) i = i + 1 do t = 1, num_bins irrigfrac_cnt(nc,nr,t) = go2(i,t) end do enddo enddo !- Convert 3-D irrigration fraction to 2-D: irrigfrac = 0. do nr = 1, LDT_rc%lnr(n) do nc = 1, LDT_rc%lnc(n) ! Calculate gridpoint column total: isum = sum(irrigfrac_cnt(nc,nr,1:num_bins)) ! write(500,*) nc,nr, isum, irrigfrac_cnt(nc,nr,2),irrigfrac_cnt(nc,nr,3) ! Estimate 2-D fraction for just irrigation: if( isum > 0. ) then ! Irrigated layer only: irrigfrac(nc,nr) = (irrigfrac_cnt(nc,nr,2) / isum) * 100. ! Irrigated + rice patty layers: ! irrigfrac(nc,nr) = ((irrigfrac_cnt(nc,nr,2)+irrigfrac_cnt(nc,nr,3)) / isum) * 100. else irrigfrac(nc,nr) = 0. endif enddo end do fgrd = irrigfrac write(LDT_logunit,fmt=*) "[INFO] Done reading MODIS/GRIPC irrigation gridcell fractions" END SUBROUTINE regrid_to_lisgrid ! ------------------------------------------------------------------------- subroutine read_global_nc4 end subroutine read_global_nc4 end subroutine read_data end module LDT_climateBCs