diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 index 32173c29d..f554a0bbf 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM40_GridComp/GEOS_CatchCNCLM40GridComp.F90 @@ -173,7 +173,7 @@ module GEOS_CatchCNCLM40GridCompMod end type OFFLINE_WRAP integer :: RUN_IRRIG, USE_ASCATZ0, Z0_FORMULATION, IRRIG_METHOD, AEROSOL_DEPOSITION, N_CONST_LAND4SNWALB -integer :: ATM_CO2, PRESCRIBE_DVG, SCALE_ALBFPAR,CHOOSEMOSFC +integer :: ATM_CO2, PRESCRIBE_DVG,CHOOSEMOSFC 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 :: CO2 @@ -297,13 +297,6 @@ subroutine SetServices ( GC, RC ) ! 3--Estimated LAI/SAI using anomalies at the beginning of the foeecast and climatological LAI/SAI call MAPL_GetResource (SCF, PRESCRIBE_DVG, label='PRESCRIBE_DVG:', DEFAULT=0 , __RC__ ) - ! SCALE_ALBFPAR: Scale CATCHCN ALBEDO and FPAR - ! 0-- NO scaling is performed - ! 1-- Scale albedo to match interannually varying MODIS NIRDF and VISDF anomaly - ! 2-- Scale albedo to match CDFs of model fPAR to MODIS CDFs of fPAR - ! 3-- Pefform above both 1 and 2 scalings - call MAPL_GetResource (SCF, SCALE_ALBFPAR, label='SCALE_ALBFPAR:', DEFAULT=0 , __RC__ ) - ! Global mean CO2 call MAPL_GetResource (SCF, CO2, label='CO2:', DEFAULT=350.e-6, __RC__ ) call MAPL_GetResource (SCF, CO2_YEAR_IN, label='CO2_YEAR:', DEFAULT= -9999, __RC__ ) @@ -5170,24 +5163,9 @@ subroutine Driver ( RC ) real, allocatable, dimension(:) :: ALBVR_tmp, ALBNR_tmp, ALBVF_tmp, ALBNF_tmp real, allocatable, dimension(:) :: SNOVR_tmp, SNONR_tmp, SNOVF_tmp, SNONF_tmp - ! Variables for FPAR scaling - ! -------------------------- - - real, save,allocatable,dimension (:,:,:,:) :: Kappa, Lambda, Mu - real, save,allocatable,dimension (:,:,:) :: MnVal, MxVal - integer, save, allocatable, dimension (:) :: modis_tid, ThisMIndex - integer :: n_modis, NTCurrent, CDFfile, infos, comms - integer, allocatable, dimension (:,:) :: modis_index - integer, allocatable, dimension (:) :: modis2cat - real , allocatable, dimension (:) :: m_lons, m_lats - real , allocatable, dimension (:,:) :: scaled_fpar, parav, parzone, unscaled_fpar - REAL , PARAMETER :: TILEINT = 2. - integer, PARAMETER :: NOCTAD = 46, NSETS = 2 - real :: CLM4_fpar, CLM4_cdf, MODIS_fpar, tmparr(1,1,1,2), & - ThisK, ThisL, ThisM, ThisMin, ThisMax, tmparr2(1,1,1), ThisFPAR, ZFPAR - character (len=ESMF_MAXSTR) :: VISMEANFILE, VISSTDFILE, NIRMEANFILE, NIRSTDFILE, FPARMEANFILE, FPARSTDFILE - real, allocatable, dimension (:) :: MODISVISmean, MODISVISstd, MODISNIRmean, MODISNIRstd, MODELFPARmean, MODELFPARstd - logical, save :: first_fpar = .true. + ! Variables for FPAR + ! -------------------------- + real , allocatable, dimension (:,:) :: parzone IAm=trim(COMP_NAME)//"::RUN2::Driver" @@ -5702,123 +5680,6 @@ subroutine Driver ( RC ) ENDIF READ_CT_CO2 ENDIF - ! OPTIONAL FPAR SCALING -! --------------------- - - if (SCALE_ALBFPAR >= 2) then - IF (ntiles > 0) THEN - INTILALIZE_FPAR_PARAM : if(first_fpar) then - - ! Initialize FPAR MODIS scale parameters - ! -------------------------------------- - -! CALL ESMF_VMGet(vm, MPICOMMUNICATOR=comms, rc=status) -! VERIFY_(status) -! call MPI_Info_create(infos, STATUS) -! call MPI_Info_set(infos, "romio_cb_read", "automatic", STATUS) - - STATUS = NF_OPEN ('FPAR_CDF_Params-M09.nc4', NF_NOWRITE, CDFfile) - STATUS = NF_INQ_DIMID (CDFfile, 'tile10D', k); VERIFY_(STATUS) - STATUS = NF_INQ_DIMLEN (CDFfile, K, n_modis) ; VERIFY_(STATUS) - - allocate (m_lons (1 : n_modis)) - allocate (m_lats (1 : n_modis)) - - STATUS = NF_GET_VARA_REAL (CDFfile, VarID(CDFfile,'lon'), (/1/), (/n_modis/), m_lons);VERIFY_(STATUS) - STATUS = NF_GET_VARA_REAL (CDFfile, VarID(CDFfile,'lat'), (/1/), (/n_modis/), m_lats);VERIFY_(STATUS) - - allocate (modis_index (1: 360/nint(TILEINT), 1: 180/nint(TILEINT))) - modis_index = -9999 - - ! vector to grid 10x10 MODIS tiles - - do i = 1, n_modis - - k = NINT (((m_lons(i) + TILEINT/2.) + 180.) / TILEINT) - n = NINT (((m_lats(i) + TILEINT/2.) + 90.) / TILEINT) - modis_index (k, n) = i - - end do - - ! for each catchment-tile overlying MODIS 10x10 tile - - allocate (modis2cat (1: NTILES)) - allocate (modis_tid (1: NTILES)) - - modis_tid = -9999 - modis2cat = 0 - - do i = 1, NTILES - - k = NINT ((CEILING (lons(i)*90./MAPL_PI)*2 + 180.) / TILEINT) - n = NINT ((CEILING (lats(i)*90./MAPL_PI)*2 + 90.) / TILEINT) - if(k <= 3) k = 3 - if(k >= 178) k = 178 - modis2cat (i) = modis_index (k,n) - - end do - - K = count(modis2cat > 0) - - allocate (unq_mask(1:K )) - allocate (loc_int (1:K )) - - loc_int = pack(modis2cat ,mask = (modis2cat > 0)) - call MAPL_Sort (loc_int) - unq_mask = .true. - - do i = 2,K - unq_mask(i) = .not.(loc_int(i) == loc_int(i-1)) - end do - - NUNQ = count(unq_mask) - - allocate (ThisIndex (1:NUNQ)) - ThisIndex = pack(loc_int, mask = unq_mask ) - - allocate (Kappa (1: NUNQ, 1: NUMPFT, 1 : NOCTAD, 1 : 2)) - allocate (Lambda(1: NUNQ, 1: NUMPFT, 1 : NOCTAD, 1 : 2)) - allocate (Mu (1: NUNQ, 1: NUMPFT, 1 : NOCTAD, 1 : 2)) - allocate (MnVal (1: NUNQ, 1: NUMPFT, 1 : NOCTAD)) - allocate (MxVal (1: NUNQ, 1: NUMPFT, 1 : NOCTAD)) - - Kappa = -9999. - Lambda = -9999. - Mu = -9999. - - do i = 1, NUNQ - - where (modis2cat == ThisIndex(i)) modis_tid = i - - end do - - do i = 1, NUNQ - do K = 1,NOCTAD - do n = 1, NUMPFT - IF (ThisIndex(i) >= 1) THEN - STATUS = NF_GET_VARA_REAL(CDFFile, VARID(CDFFile,'Kappa' ),(/ThisIndex(i),N,K,1/), (/1,1,1,2/), tmparr);VERIFY_(STATUS) - Kappa (i,N,K,:) = tmparr (1,1,1,:) - STATUS = NF_GET_VARA_REAL(CDFFile, VARID(CDFFile,'Lambda'),(/ThisIndex(i),N,K,1/), (/1,1,1,2/), tmparr);VERIFY_(STATUS) - Lambda(i,N,K,:) = tmparr (1,1,1,:) - STATUS = NF_GET_VARA_REAL(CDFFile, VARID(CDFFile,'Mu' ),(/ThisIndex(i),N,K,1/), (/1,1,1,2/), tmparr);VERIFY_(STATUS) - Mu (i,N,K,:) = tmparr (1,1,1,:) - STATUS = NF_GET_VARA_REAL(CDFFile, VARID(CDFFile,'MinVal'),(/ThisIndex(i),N,K/), (/1,1,1/), tmparr2);VERIFY_(STATUS) - MnVal(i,N,K) = tmparr2 (1,1,1) - STATUS = NF_GET_VARA_REAL(CDFFile, VARID(CDFFile,'MaxVal'),(/ThisIndex(i),N,K/), (/1,1,1/), tmparr2);VERIFY_(STATUS) - MxVal(i,N,K) = tmparr2 (1,1,1) - ENDIF - end do - end do - end do - status = NF_CLOSE (CDFFile) - - deallocate ( modis2cat, unq_mask, loc_int, modis_index, m_lons, m_lats) - - first_fpar = .false. - - endif INTILALIZE_FPAR_PARAM - endif - end if ! -------------------------------------------------------------------------- ! ALLOCATE LOCAL POINTERS @@ -6352,11 +6213,7 @@ subroutine Driver ( RC ) allocate( car1(ntiles) ) allocate( car2(ntiles) ) allocate( car4(ntiles) ) - allocate( parzone(ntiles,nveg) ) allocate( para(ntiles) ) - allocate( parav(ntiles,nveg) ) - allocate (scaled_fpar (NTILES,NVEG)) - allocate (unscaled_fpar(NTILES,NVEG)) allocate ( totwat(ntiles) ) if(.not. allocated(npp )) allocate( npp(ntiles) ) if(.not. allocated(gpp )) allocate( gpp(ntiles) ) @@ -6382,6 +6239,7 @@ subroutine Driver ( RC ) allocate( psnsunx(ntiles,nveg) ) allocate( psnshax(ntiles,nveg) ) + allocate( parzone(ntiles,nveg) ) allocate( sifsunx(ntiles,nveg) ) allocate( sifshax(ntiles,nveg) ) allocate( laisunx(ntiles,nveg) ) @@ -6657,8 +6515,6 @@ subroutine Driver ( RC ) end do para(:) = 0. ! zero out absorbed PAR summing array - parav(:, :) = 0. ! - scaled_fpar = 1. do nz = 1,nzone @@ -6754,8 +6610,8 @@ subroutine Driver ( RC ) do nv = 1,nveg para(:) = para(:) + parzone(:,nv)*wtzone(:,nz)*fvez(:,nv) - parav(:,nv) = parav (:,nv) + parzone(:,nv)*wtzone(:,nz) end do + if(associated(BTRANT)) btrant(:) = btrant(:) + btran(:)*wtzone(:,nz) if(associated(SIF)) then do nv = 1,nveg @@ -6765,157 +6621,6 @@ subroutine Driver ( RC ) end do - do nv = 1,nveg - unscaled_fpar (:,nv) = parav (:,nv)/ (DRPAR(:) + DFPAR(:) + 1.e-20) - end do - - NTCurrent = CEILING (real (dofyr) / 8.) - - ! FPAR scaling to match MODIS CDF - ! ------------------------------- - - DO_FS1 : if (SCALE_ALBFPAR >= 2) then - - IF (ntiles > 0) THEN - - NT_LOOP1 : do n = 1,NTILES - - NV_LOOP1 : do nv = 1,nveg - - CLM4_fpar = parav (n,nv) / (DRPAR (n) + DFPAR (n) + 1.e-20) - K = -1 - - if(CLM4_fpar > 0.) then - - k = NINT(ITY(N,nv)) - if(minval(Kappa (modis_tid (n), k, NTCurrent, :)) < 0.) then - k = -1 - if(nv == 1) k = NINT(ITY(N,2)) - if(nv == 2) k = NINT(ITY(N,1)) - if(nv == 3) k = NINT(ITY(N,4)) - if(nv == 4) k = NINT(ITY(N,3)) - if(minval(Kappa (modis_tid (n), k, NTCurrent, :)) < 0.) k = -1 - if((K == -1).and.(nv > 2)) then - if(minval(Kappa (modis_tid (n), NINT(ITY(N,2)), NTCurrent, :)) > 0.) k = NINT(ITY(N,2)) - if(minval(Kappa (modis_tid (n), NINT(ITY(N,1)), NTCurrent, :)) > 0.) k = NINT(ITY(N,1)) - endif - endif - - endif - - if((K > 0).and.(CLM4_fpar > 0)) then - - ! Computing probability of CLM4 FPAR - - ThisK = Kappa (modis_tid (n), k, NTCurrent, 2) - ThisL = Lambda (modis_tid (n), k, NTCurrent, 2) - ThisM = Mu (modis_tid (n), k, NTCurrent, 2) - ThisMin = MnVal (modis_tid (n), k, NTCurrent) - ThisMax = MxVal (modis_tid (n), k, NTCurrent) - - if (CLM4_fpar < ThisMin) CLM4_fpar = ThisMin - if (CLM4_fpar > ThisMax) CLM4_fpar = ThisMax - if((ThisL == 0.).or.(ThisM == 0.)) print *,thisK,ThisL, ThisM, CLM4_fpar, ThisMin, ThisMax - if((ThisL == 0.).or.(ThisM == 0.)) print *,n,k,NTCurrent,modis_tid (n) - CLM4_cdf = ThisK * betai (ThisL, ThisM, (CLM4_fpar - ThisMin)/ThisMax) - - ! Computing corresponding MODIS FPAR for the same probability - - ThisK = Kappa (modis_tid (n), k, NTCurrent, 1) - ThisL = Lambda (modis_tid (n), k, NTCurrent, 1) - ThisM = Mu (modis_tid (n), k, NTCurrent, 1) - ThisMin = MnVal (modis_tid (n), k, NTCurrent) - ThisMax = MxVal (modis_tid (n), k, NTCurrent) - - scaled_fpar (n,nv) = cdf2fpar (CLM4_cdf, ThisK, ThisL, ThisM, ThisMin, ThisMax) - if((scaled_fpar (n,nv) > 1.).or.(scaled_fpar (n,nv) < 0.)) then - print *, 'PROB 1', CLM4_cdf, ThisK, ThisL, ThisM, ThisMin, ThisMax, scaled_fpar (n,nv) - endif - - scaled_fpar (n,nv) = scaled_fpar (n,nv) / (CLM4_fpar + 1.e-20) - - endif - end do NV_LOOP1 - - end do NT_LOOP1 - - para (:) = 0. ! zero out absorbed PAR summing array - parav = 0. - - if(associated(BTRANT)) btrant = 0. - if(associated(SIF)) sif = 0. - - do nz = 1,num_zon - - if(nz == 1) then - btran = btran1 - tcx = tx1 - qax = qx1 - endif - - if(nz == 2) then - btran = btran2 - tcx = tx2 - qax = qx2 - endif - - if(nz == 3) then - btran = btran3 - tcx = tx3 - qax = qx3 - endif - - do nv = 1,num_veg - elaz(:,nv) = elai(:,nv,nz) - esaz(:,nv) = esai(:,nv,nz) - ityz(:,nv) = ityp(:,nv,nz) - fvez(:,nv) = fveg(:,nv,nz) - end do - - do n = 1,NTILES - if(tp1(n) < (Tzero-0.01)) btran(n) = 0. ! no photosynthesis if ground fully frozen - end do - - call compute_rc(NTILES,nveg,TCx,QAx, & - TA, PS, ZTH,DRPAR,DFPAR, & - elaz,esaz,ityz,fvez,btran,fwet, & - RCx,RCxDT,RCxDQ,psnsunx,psnshax,laisunx,laishax, & - dayl_fac,co2v,dtc,dea,parzone,sifsunx,sifshax, & - fpar_sf = scaled_fpar ) - - rc00(:,nz) = rcx(:) - rcdt(:,nz) = rcxdt(:) - rcdq(:,nz) = rcxdq(:) - - psnsun(:,:,nz) = psnsunx(:,:) - psnsha(:,:,nz) = psnshax(:,:) - laisun(:,:,nz) = laisunx(:,:) - laisha(:,:,nz) = laishax(:,:) - - do nv = 1,nveg - para(:) = para(:) + parzone(:,nv)*wtzone(:,nz)*fvez(:,nv) - parav(:,nv) = parav (:,nv) + parzone(:,nv)*wtzone(:,nz) - end do - - if(associated(BTRANT)) btrant(:) = btrant(:) + btran(:)*wtzone(:,nz) - if(associated(SIF)) then - do nv = 1,nveg - sif(:) = sif(:) + wtzone(:,nz)*fvez(:,nv)*(sifsunx(:,nv)*laisunx(:,nv) + sifshax(:,nv)*laishax(:,nv)) - end do - endif - - end do - - endif - - endif DO_FS1 - - ! Below we are recycling the scaled_fpar array - from this point, it contains fpar scaled or otherwise - ! ---------------------------------------------------------------------------------------------------- - - do nv = 1,nveg - scaled_fpar (:,nv) = parav (:,nv)/ (DRPAR(:) + DFPAR(:) + 1.e-20) - end do if(associated(CNCO2)) CNCO2 = CO2V * 1e6 deallocate (co2v) @@ -6940,40 +6645,6 @@ subroutine Driver ( RC ) BGALBVR, BGALBVF, BGALBNR, BGALBNF, & ! gkw: MODIS soil background albedo ALBVR, ALBNR, ALBVF, ALBNF, MODIS_SCALE=.TRUE. ) ! instantaneous snow-free albedos on tiles - if ((SCALE_ALBFPAR == 1).OR.(SCALE_ALBFPAR == 3)) then - - if(.not.allocated (MODISVISmean )) allocate (MODISVISmean (1:NTILES)) - if(.not.allocated (MODISVISstd )) allocate (MODISVISstd (1:NTILES)) - if(.not.allocated (MODISNIRmean )) allocate (MODISNIRmean (1:NTILES)) - if(.not.allocated (MODISNIRstd )) allocate (MODISNIRstd (1:NTILES)) - if(.not.allocated (MODELFPARmean)) allocate (MODELFPARmean (1:NTILES)) - if(.not.allocated (MODELFPARstd )) allocate (MODELFPARstd (1:NTILES)) - - if(ntiles > 0) then - - call MAPL_GetResource(MAPL,VISMEANFILE , label = 'VISMEAN_FILE:' , default = 'MODISVISmean.dat' , RC=STATUS ) ; VERIFY_(STATUS) - call MAPL_GetResource(MAPL,VISSTDFILE , label = 'VISSTD_FILE:' , default = 'MODISVISstd.dat' , RC=STATUS ) ; VERIFY_(STATUS) - call MAPL_GetResource(MAPL,NIRMEANFILE , label = 'NIRMEAN_FILE:' , default = 'MODISNIRmean.dat' , RC=STATUS ) ; VERIFY_(STATUS) - call MAPL_GetResource(MAPL,NIRSTDFILE , label = 'NIRSTD_FILE:' , default = 'MODISNIRstd.dat' , RC=STATUS ) ; VERIFY_(STATUS) - - call MAPL_GetResource(MAPL,FPARMEANFILE , label = 'MODELFPARMEAN_FILE:', default = 'MODELFPARmean.dat' , RC=STATUS ) ; VERIFY_(STATUS) - call MAPL_GetResource(MAPL,FPARSTDFILE , label = 'MODELFPARSTD_FILE:' , default = 'MODELFPARstd.dat' , RC=STATUS ) ; VERIFY_(STATUS) - - call MAPL_ReadForcing(MAPL,'MODISVISmean' ,VISMEANFILE ,CURRENT_TIME,MODISVISmean ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - call MAPL_ReadForcing(MAPL,'MODISVISstd' ,VISSTDFILE ,CURRENT_TIME,MODISVISstd ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - call MAPL_ReadForcing(MAPL,'MODISNIRmean' ,NIRMEANFILE ,CURRENT_TIME,MODISNIRmean ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - call MAPL_ReadForcing(MAPL,'MODISNIRstd' ,NIRSTDFILE ,CURRENT_TIME,MODISNIRstd ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - call MAPL_ReadForcing(MAPL,'MODELFPARmean',FPARMEANFILE,CURRENT_TIME,MODELFPARmean ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - call MAPL_ReadForcing(MAPL,'MODELFPARstd' ,FPARSTDFILE ,CURRENT_TIME,MODELFPARstd ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - - do n = 1,NTILES - ThisFPAR = (unscaled_fpar(n,1)*FVG(N,1) + unscaled_fpar(n,2)*FVG(N,2))/(FVG(N,1) + FVG(N,2) + 1.e-20) - ZFPAR = (ThisFPAR - MODELFPARmean (n)) / MODELFPARstd (n) - ALBVF(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISVISstd(n) + MODISVISmean (n))) - ALBNF(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISNIRstd(n) + MODISNIRmean (n))) - end do - endif - endif call STIEGLITZSNOW_CALC_TPSNOW(NTILES, HTSNNN(1,:), WESNN(1,:), TPSN1OUT1, FICE1) TPSN1OUT1 = TPSN1OUT1 + Tzero @@ -6990,16 +6661,6 @@ subroutine Driver ( RC ) BGALBVR, BGALBVF, BGALBNR, BGALBNF, & ! gkw: MODIS soil background albedo ALBVR_tmp, ALBNR_tmp, ALBVF_tmp, ALBNF_tmp, MODIS_SCALE=.TRUE. ) ! instantaneous snow-free albedos on tiles - if ((SCALE_ALBFPAR == 1).OR.(SCALE_ALBFPAR == 3)) then - if(ntiles > 0) then - do n = 1,NTILES - ThisFPAR = (unscaled_fpar(n,3)*FVG(N,3) + unscaled_fpar(n,4)*FVG(N,4))/(FVG(N,3) + FVG(N,4) + 1.e-20) - ZFPAR = (ThisFPAR - MODELFPARmean (n)) / MODELFPARstd (n) - ALBVF_tmp(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISVISstd(n) + MODISVISmean (n))) - ALBNF_tmp(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISNIRstd(n) + MODISNIRmean (n))) - end do - endif - endif call SNOW_ALBEDO(NTILES,N_snow, N_CONST_LAND4SNWALB, VEG2, LAI2, ZTH, & RHOFS, & @@ -7643,13 +7304,6 @@ subroutine Driver ( RC ) BGALBVR, BGALBVF, BGALBNR, BGALBNF, & ! gkw: MODIS soil background albedo ALBVR, ALBNR, ALBVF, ALBNF, MODIS_SCALE=.TRUE. ) ! instantaneous snow-free albedos on tiles - if ((SCALE_ALBFPAR == 1).OR.(SCALE_ALBFPAR == 3)) then - do n = 1,NTILES - ThisFPAR = (unscaled_fpar(n,1)*FVG(N,1) + unscaled_fpar(n,2)*FVG(N,2))/(FVG(N,1) + FVG(N,2) + 1.e-20) - ZFPAR = (ThisFPAR - MODELFPARmean (n)) / MODELFPARstd (n) - ALBVF(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISVISstd(n) + MODISVISmean (n))) - end do - endif call STIEGLITZSNOW_CALC_TPSNOW(NTILES, HTSNNN(1,:), WESNN(1,:), TPSN1OUT1, FICE1) TPSN1OUT1 = TPSN1OUT1 + Tzero @@ -7666,17 +7320,6 @@ subroutine Driver ( RC ) BGALBVR, BGALBVF, BGALBNR, BGALBNF, & ! gkw: MODIS soil background albedo ALBVR_tmp, ALBNR_tmp, ALBVF_tmp, ALBNF_tmp, MODIS_SCALE=.TRUE. ) ! instantaneous snow-free albedos on tiles - if ((SCALE_ALBFPAR == 1).OR.(SCALE_ALBFPAR == 3)) then - if(ntiles > 0) then - do n = 1,NTILES - ThisFPAR = (unscaled_fpar(n,3)*FVG(N,3) + unscaled_fpar(n,4)*FVG(N,4))/(FVG(N,3) + FVG(N,4) + 1.e-20) - ZFPAR = (ThisFPAR - MODELFPARmean (n)) / MODELFPARstd (n) - ALBVF_tmp(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISVISstd(n) + MODISVISmean (n))) - ALBNF_tmp(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISNIRstd(n) + MODISNIRmean (n))) - end do - endif - if(allocated (MODISVISmean)) deallocate (MODISVISmean, MODISVISstd, MODISNIRmean, MODISNIRstd, MODELFPARmean, MODELFPARstd) - endif call SNOW_ALBEDO(NTILES,N_snow, N_CONST_LAND4SNWALB, VEG2, LAI2, ZTH, & RHOFS, & @@ -8023,11 +7666,7 @@ subroutine Driver ( RC ) deallocate( car1 ) deallocate( car2 ) deallocate( car4 ) - deallocate( parzone ) deallocate( para ) - deallocate( parav ) - deallocate(scaled_fpar) - deallocate(UNscaled_fpar) deallocate( totwat ) deallocate( dayl ) deallocate(dayl_fac ) @@ -8042,6 +7681,7 @@ subroutine Driver ( RC ) deallocate( psnsunx ) deallocate( psnshax ) deallocate( sifsunx ) + deallocate( parzone ) deallocate( sifshax ) deallocate( laisunx ) deallocate( laishax ) @@ -8070,128 +7710,109 @@ subroutine Driver ( RC ) end subroutine Driver - ! ----------------- routines for CDF scaling ------------------- - - REAL FUNCTION cdf2fpar (cdf, k,l, m, m1, m2) - - REAL, intent (in) :: cdf, k,l,m, m1, m2 - REAL :: x, ThisCDF, ThisFPAR - integer, parameter :: nBINS = 40 - - x = real (nBINS) - ThisCDF = 1. - - do while (ThisCDF >= cdf) - ThisFPAR = 1. - (real(nbins)-x)/real(nbins) - 1./2./real(nbins) - ThisCDF = K * betai (L, M, ThisFPAR) - x = x - 1. - if(x == 0) exit - end do - - cdf2fpar = ThisFPAR * m2 + m1 - if(cdf2fpar > m2) cdf2fpar = m2 - if(cdf2fpar < m1) cdf2fpar = m1 - return - - END FUNCTION cdf2fpar - ! --------------------------------------------------------- +! Commented out functions betai(), betacf(), and gammln(). +! These functions are not used and were reproduced identically in +! GEOS_CatchCNCLM40GridComp.F90 and in GEOS_CatchCNCLM45GridComp.F90. +! Another copy was in GEOScatchCN_GridComp/utils/math_routines.F90 but +! there function betai() was missing the restriction 0.0125 0.9875) x = 0.9875 - - if(x.lt.0..or.x.gt.1.)print *, 'bad argument x in betai',x - if(x.lt.0..or.x.gt.1.)stop - if(x.eq.0..or.x.eq.1.)then - bt=0. - else - bt=exp(gammln(a+b)-gammln(a)-gammln(b) & - +a*log(x)+b*log(1.-x)) - endif - - if(x.lt.(a+1.)/(a+b+2.))then - betai=bt*betacf(a,b,x)/a - return - else - betai=1.-bt*betacf(b,a,1.-x)/b - return - endif - - END FUNCTION betai - - ! ------------------------------------------------------- - - FUNCTION betacf(a,b,x) - - INTEGER MAXIT - REAL betacf,a,b,x,EPS,FPMIN - PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30) - INTEGER m,m2 - REAL aa,c,d,del,h,qab,qam,qap - - qab=a+b - qap=a+1. - qam=a-1. - c=1. - d=1.-qab*x/qap - - if(abs(d).lt.FPMIN)d=FPMIN - d=1./d - h=d - do m=1,MAXIT - m2=2*m - aa=m*(b-m)*x/((qam+m2)*(a+m2)) - d=1.+aa*d - if(abs(d).lt.FPMIN)d=FPMIN - c=1.+aa/c - if(abs(c).lt.FPMIN)c=FPMIN - d=1./d - h=h*d*c - aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2)) - d=1.+aa*d - if(abs(d).lt.FPMIN)d=FPMIN - c=1.+aa/c - if(abs(c).lt.FPMIN)c=FPMIN - d=1./d - del=d*c - h=h*del - if(abs(del-1.).lt.EPS)exit - enddo - betacf=h - return - - END FUNCTION betacf - - ! -------------------------------------------------------------- - - FUNCTION gammln(xx) - - REAL gammln,xx - INTEGER j - DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) - - SAVE cof,stp - DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, & - 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, & - -.5395239384953d-5,2.5066282746310005d0/ - x=xx - y=x - tmp=x+5.5d0 - tmp=(x+0.5d0)*log(tmp)-tmp - ser=1.000000000190015d0 - do j=1,6 - y=y+1.d0 - ser=ser+cof(j)/y - enddo - gammln=tmp+log(stp*ser/x) - return - - END FUNCTION gammln +!! FUNCTION betai(a,b,x) +!! REAL betai,a,b,x +!! REAL bt +!! !external gammln +!! +!! if (x < 0.0125) x = 0.0125 +!! if (x > 0.9875) x = 0.9875 +!! +!! if(x.lt.0..or.x.gt.1.)print *, 'bad argument x in betai',x +!! if(x.lt.0..or.x.gt.1.)stop +!! if(x.eq.0..or.x.eq.1.)then +!! bt=0. +!! else +!! bt=exp(gammln(a+b)-gammln(a)-gammln(b) & +!! +a*log(x)+b*log(1.-x)) +!! endif +!! +!! if(x.lt.(a+1.)/(a+b+2.))then +!! betai=bt*betacf(a,b,x)/a +!! return +!! else +!! betai=1.-bt*betacf(b,a,1.-x)/b +!! return +!! endif +!! +!! END FUNCTION betai +!! +!! ! ------------------------------------------------------- +!! +!! FUNCTION betacf(a,b,x) +!! +!! INTEGER MAXIT +!! REAL betacf,a,b,x,EPS,FPMIN +!! PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30) +!! INTEGER m,m2 +!! REAL aa,c,d,del,h,qab,qam,qap +!! +!! qab=a+b +!! qap=a+1. +!! qam=a-1. +!! c=1. +!! d=1.-qab*x/qap +!! +!! if(abs(d).lt.FPMIN)d=FPMIN +!! d=1./d +!! h=d +!! do m=1,MAXIT +!! m2=2*m +!! aa=m*(b-m)*x/((qam+m2)*(a+m2)) +!! d=1.+aa*d +!! if(abs(d).lt.FPMIN)d=FPMIN +!! c=1.+aa/c +!! if(abs(c).lt.FPMIN)c=FPMIN +!! d=1./d +!! h=h*d*c +!! aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2)) +!! d=1.+aa*d +!! if(abs(d).lt.FPMIN)d=FPMIN +!! c=1.+aa/c +!! if(abs(c).lt.FPMIN)c=FPMIN +!! d=1./d +!! del=d*c +!! h=h*del +!! if(abs(del-1.).lt.EPS)exit +!! enddo +!! betacf=h +!! return +!! +!! END FUNCTION betacf +!! +!! ! -------------------------------------------------------------- +!! +!! FUNCTION gammln(xx) +!! +!! REAL gammln,xx +!! INTEGER j +!! DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) +!! +!! SAVE cof,stp +!! DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, & +!! 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, & +!! -.5395239384953d-5,2.5066282746310005d0/ +!! x=xx +!! y=x +!! tmp=x+5.5d0 +!! tmp=(x+0.5d0)*log(tmp)-tmp +!! ser=1.000000000190015d0 +!! do j=1,6 +!! y=y+1.d0 +!! ser=ser+cof(j)/y +!! enddo +!! gammln=tmp+log(stp*ser/x) +!! return +!! +!! END FUNCTION gammln ! -------------------------------------------------------------- diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 index 5caaa0425..f1b24d30e 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/GEOScatchCNCLM45_GridComp/GEOS_CatchCNCLM45GridComp.F90 @@ -174,7 +174,7 @@ module GEOS_CatchCNCLM45GridCompMod end type OFFLINE_WRAP integer :: RUN_IRRIG, USE_ASCATZ0, Z0_FORMULATION, IRRIG_METHOD, AEROSOL_DEPOSITION, N_CONST_LAND4SNWALB -integer :: ATM_CO2, SCALE_ALBFPAR,CHOOSEMOSFC +integer :: ATM_CO2, CHOOSEMOSFC 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 :: CO2 @@ -291,12 +291,6 @@ subroutine SetServices ( GC, RC ) ! 4: import AGCM model CO2 (AGCM only) call MAPL_GetResource (SCF, ATM_CO2, label='ATM_CO2:', DEFAULT=2 , __RC__ ) - ! SCALE_ALBFPAR: Scale CATCHCN ALBEDO and FPAR - ! 0-- NO scaling is performed - ! 1-- Scale albedo to match interannually varying MODIS NIRDF and VISDF anomaly - ! 2-- Scale albedo to match CDFs of model fPAR to MODIS CDFs of fPAR - ! 3-- Pefform above both 1 and 2 scalings - call MAPL_GetResource (SCF, SCALE_ALBFPAR, label='SCALE_ALBFPAR:', DEFAULT=0 , __RC__ ) ! Global mean CO2 call MAPL_GetResource (SCF, CO2, label='CO2:', DEFAULT=350.e-6, __RC__ ) @@ -4453,7 +4447,6 @@ subroutine RUN2 ( GC, IMPORT, EXPORT, CLOCK, RC ) type(MAPL_MetaComp),pointer :: MAPL type(ESMF_Alarm) :: ALARM - integer :: IM,JM integer :: incl_Louis_extra_derivs @@ -5155,24 +5148,9 @@ subroutine Driver ( RC ) real, allocatable, dimension(:) :: ALBVR_tmp, ALBNR_tmp, ALBVF_tmp, ALBNF_tmp real, allocatable, dimension(:) :: SNOVR_tmp, SNONR_tmp, SNOVF_tmp, SNONF_tmp - ! Variables for FPAR scaling - ! -------------------------- - - real, save,allocatable,dimension (:,:,:,:) :: Kappa, Lambda, Mu - real, save,allocatable,dimension (:,:,:) :: MnVal, MxVal - integer, save, allocatable, dimension (:) :: modis_tid, ThisMIndex - integer :: n_modis, NTCurrent, CDFfile, infos, comms - integer, allocatable, dimension (:,:) :: modis_index - integer, allocatable, dimension (:) :: modis2cat - real , allocatable, dimension (:) :: m_lons, m_lats - real , allocatable, dimension (:,:) :: scaled_fpar, parav, parzone, unscaled_fpar - REAL , PARAMETER :: TILEINT = 2. - integer, PARAMETER :: NOCTAD = 46, NSETS = 2 - real :: CLM4_fpar, CLM4_cdf, MODIS_fpar, tmparr(1,1,1,2), & - ThisK, ThisL, ThisM, ThisMin, ThisMax, tmparr2(1,1,1), ThisFPAR, ZFPAR - character (len=ESMF_MAXSTR) :: VISMEANFILE, VISSTDFILE, NIRMEANFILE, NIRSTDFILE, FPARMEANFILE, FPARSTDFILE - real, allocatable, dimension (:) :: MODISVISmean, MODISVISstd, MODISNIRmean, MODISNIRstd, MODELFPARmean, MODELFPARstd - logical, save :: first_fpar = .true. + ! Variables for FPAR + ! -------------------------- + real , allocatable, dimension (:,:) :: parzone IAm=trim(COMP_NAME)//"::RUN2::Driver" @@ -5699,123 +5677,6 @@ subroutine Driver ( RC ) ENDIF READ_CT_CO2 ENDIF - ! OPTIONAL FPAR SCALING -! --------------------- - - if (SCALE_ALBFPAR >= 2) then - IF (ntiles > 0) THEN - INTILALIZE_FPAR_PARAM : if(first_fpar) then - - ! Initialize FPAR MODIS scale parameters - ! -------------------------------------- - -! CALL ESMF_VMGet(vm, MPICOMMUNICATOR=comms, rc=status) -! VERIFY_(status) -! call MPI_Info_create(infos, STATUS) -! call MPI_Info_set(infos, "romio_cb_read", "automatic", STATUS) - - STATUS = NF_OPEN ('FPAR_CDF_Params-M09.nc4', NF_NOWRITE, CDFfile) - STATUS = NF_INQ_DIMID (CDFfile, 'tile10D', k); VERIFY_(STATUS) - STATUS = NF_INQ_DIMLEN (CDFfile, K, n_modis) ; VERIFY_(STATUS) - - allocate (m_lons (1 : n_modis)) - allocate (m_lats (1 : n_modis)) - - STATUS = NF_GET_VARA_REAL (CDFfile, VarID(CDFfile,'lon'), (/1/), (/n_modis/), m_lons);VERIFY_(STATUS) - STATUS = NF_GET_VARA_REAL (CDFfile, VarID(CDFfile,'lat'), (/1/), (/n_modis/), m_lats);VERIFY_(STATUS) - - allocate (modis_index (1: 360/nint(TILEINT), 1: 180/nint(TILEINT))) - modis_index = -9999 - - ! vector to grid 10x10 MODIS tiles - - do i = 1, n_modis - - k = NINT (((m_lons(i) + TILEINT/2.) + 180.) / TILEINT) - n = NINT (((m_lats(i) + TILEINT/2.) + 90.) / TILEINT) - modis_index (k, n) = i - - end do - - ! for each catchment-tile overlying MODIS 10x10 tile - - allocate (modis2cat (1: NTILES)) - allocate (modis_tid (1: NTILES)) - - modis_tid = -9999 - modis2cat = 0 - - do i = 1, NTILES - - k = NINT ((CEILING (lons(i)*90./MAPL_PI)*2 + 180.) / TILEINT) - n = NINT ((CEILING (lats(i)*90./MAPL_PI)*2 + 90.) / TILEINT) - if(k <= 3) k = 3 - if(k >= 178) k = 178 - modis2cat (i) = modis_index (k,n) - - end do - - K = count(modis2cat > 0) - - allocate (unq_mask(1:K )) - allocate (loc_int (1:K )) - - loc_int = pack(modis2cat ,mask = (modis2cat > 0)) - call MAPL_Sort (loc_int) - unq_mask = .true. - - do i = 2,K - unq_mask(i) = .not.(loc_int(i) == loc_int(i-1)) - end do - - NUNQ = count(unq_mask) - - allocate (ThisIndex (1:NUNQ)) - ThisIndex = pack(loc_int, mask = unq_mask ) - - allocate (Kappa (1: NUNQ, 1: NUMPFT, 1 : NOCTAD, 1 : 2)) - allocate (Lambda(1: NUNQ, 1: NUMPFT, 1 : NOCTAD, 1 : 2)) - allocate (Mu (1: NUNQ, 1: NUMPFT, 1 : NOCTAD, 1 : 2)) - allocate (MnVal (1: NUNQ, 1: NUMPFT, 1 : NOCTAD)) - allocate (MxVal (1: NUNQ, 1: NUMPFT, 1 : NOCTAD)) - - Kappa = -9999. - Lambda = -9999. - Mu = -9999. - - do i = 1, NUNQ - - where (modis2cat == ThisIndex(i)) modis_tid = i - - end do - - do i = 1, NUNQ - do K = 1,NOCTAD - do n = 1, NUMPFT - IF (ThisIndex(i) >= 1) THEN - STATUS = NF_GET_VARA_REAL(CDFFile, VARID(CDFFile,'Kappa' ),(/ThisIndex(i),N,K,1/), (/1,1,1,2/), tmparr);VERIFY_(STATUS) - Kappa (i,N,K,:) = tmparr (1,1,1,:) - STATUS = NF_GET_VARA_REAL(CDFFile, VARID(CDFFile,'Lambda'),(/ThisIndex(i),N,K,1/), (/1,1,1,2/), tmparr);VERIFY_(STATUS) - Lambda(i,N,K,:) = tmparr (1,1,1,:) - STATUS = NF_GET_VARA_REAL(CDFFile, VARID(CDFFile,'Mu' ),(/ThisIndex(i),N,K,1/), (/1,1,1,2/), tmparr);VERIFY_(STATUS) - Mu (i,N,K,:) = tmparr (1,1,1,:) - STATUS = NF_GET_VARA_REAL(CDFFile, VARID(CDFFile,'MinVal'),(/ThisIndex(i),N,K/), (/1,1,1/), tmparr2);VERIFY_(STATUS) - MnVal(i,N,K) = tmparr2 (1,1,1) - STATUS = NF_GET_VARA_REAL(CDFFile, VARID(CDFFile,'MaxVal'),(/ThisIndex(i),N,K/), (/1,1,1/), tmparr2);VERIFY_(STATUS) - MxVal(i,N,K) = tmparr2 (1,1,1) - ENDIF - end do - end do - end do - status = NF_CLOSE (CDFFile) - - deallocate ( modis2cat, unq_mask, loc_int, modis_index, m_lons, m_lats) - - first_fpar = .false. - - endif INTILALIZE_FPAR_PARAM - endif - end if ! -------------------------------------------------------------------------- ! ALLOCATE LOCAL POINTERS @@ -6354,11 +6215,7 @@ subroutine Driver ( RC ) allocate( car1(ntiles) ) allocate( car2(ntiles) ) allocate( car4(ntiles) ) - allocate( parzone(ntiles,nveg) ) allocate( para(ntiles) ) - allocate( parav(ntiles,nveg) ) - allocate (scaled_fpar (NTILES,NVEG)) - allocate (unscaled_fpar(NTILES,NVEG)) allocate ( totwat(ntiles) ) if(.not. allocated(npp )) allocate( npp(ntiles) ) if(.not. allocated(gpp )) allocate( gpp(ntiles) ) @@ -6427,6 +6284,7 @@ subroutine Driver ( RC ) allocate( psnsunx(ntiles,nveg) ) allocate( psnshax(ntiles,nveg) ) allocate( sifsunx(ntiles,nveg) ) + allocate( parzone(ntiles,nveg) ) allocate( sifshax(ntiles,nveg) ) allocate( laisunx(ntiles,nveg) ) allocate( laishax(ntiles,nveg) ) @@ -6767,8 +6625,6 @@ subroutine Driver ( RC ) end do para(:) = 0. ! zero out absorbed PAR summing array - parav(:, :) = 0. ! - scaled_fpar = 1. do nz = 1,nzone @@ -6904,8 +6760,8 @@ subroutine Driver ( RC ) do nv = 1,nveg para(:) = para(:) + parzone(:,nv)*wtzone(:,nz)*fvez(:,nv) - parav(:,nv) = parav (:,nv) + parzone(:,nv)*wtzone(:,nz) end do + if(associated(BTRANT)) btrant(:) = btrant(:) + btran(:)*wtzone(:,nz) ! NOTE: btran here doesn't reflect the modification to btran for soybean (and nbrdlf_dcd_tmp_shrub if CNDV is on) in subroutine Photosynthesis. if(associated(SIF)) then do nv = 1,nveg @@ -6915,159 +6771,6 @@ subroutine Driver ( RC ) end do - do nv = 1,nveg - unscaled_fpar (:,nv) = parav (:,nv)/ (DRPAR(:) + DFPAR(:) + 1.e-20) - end do - - NTCurrent = CEILING (real (dofyr) / 8.) - - ! FPAR scaling to match MODIS CDF - ! ------------------------------- - - DO_FS1 : if (SCALE_ALBFPAR >= 2) then - - IF (ntiles > 0) THEN - - NT_LOOP1 : do n = 1,NTILES - - NV_LOOP1 : do nv = 1,nveg - - CLM4_fpar = parav (n,nv) / (DRPAR (n) + DFPAR (n) + 1.e-20) - K = -1 - - if(CLM4_fpar > 0.) then - - k = NINT(ITY(N,nv)) - if(minval(Kappa (modis_tid (n), k, NTCurrent, :)) < 0.) then - k = -1 - if(nv == 1) k = NINT(ITY(N,2)) - if(nv == 2) k = NINT(ITY(N,1)) - if(nv == 3) k = NINT(ITY(N,4)) - if(nv == 4) k = NINT(ITY(N,3)) - if(minval(Kappa (modis_tid (n), k, NTCurrent, :)) < 0.) k = -1 - if((K == -1).and.(nv > 2)) then - if(minval(Kappa (modis_tid (n), NINT(ITY(N,2)), NTCurrent, :)) > 0.) k = NINT(ITY(N,2)) - if(minval(Kappa (modis_tid (n), NINT(ITY(N,1)), NTCurrent, :)) > 0.) k = NINT(ITY(N,1)) - endif - endif - - endif - - if((K > 0).and.(CLM4_fpar > 0)) then - - ! Computing probability of CLM4 FPAR - - ThisK = Kappa (modis_tid (n), k, NTCurrent, 2) - ThisL = Lambda (modis_tid (n), k, NTCurrent, 2) - ThisM = Mu (modis_tid (n), k, NTCurrent, 2) - ThisMin = MnVal (modis_tid (n), k, NTCurrent) - ThisMax = MxVal (modis_tid (n), k, NTCurrent) - - if (CLM4_fpar < ThisMin) CLM4_fpar = ThisMin - if (CLM4_fpar > ThisMax) CLM4_fpar = ThisMax - if((ThisL == 0.).or.(ThisM == 0.)) print *,thisK,ThisL, ThisM, CLM4_fpar, ThisMin, ThisMax - if((ThisL == 0.).or.(ThisM == 0.)) print *,n,k,NTCurrent,modis_tid (n) - CLM4_cdf = ThisK * betai (ThisL, ThisM, (CLM4_fpar - ThisMin)/ThisMax) - - ! Computing corresponding MODIS FPAR for the same probability - - ThisK = Kappa (modis_tid (n), k, NTCurrent, 1) - ThisL = Lambda (modis_tid (n), k, NTCurrent, 1) - ThisM = Mu (modis_tid (n), k, NTCurrent, 1) - ThisMin = MnVal (modis_tid (n), k, NTCurrent) - ThisMax = MxVal (modis_tid (n), k, NTCurrent) - - scaled_fpar (n,nv) = cdf2fpar (CLM4_cdf, ThisK, ThisL, ThisM, ThisMin, ThisMax) - if((scaled_fpar (n,nv) > 1.).or.(scaled_fpar (n,nv) < 0.)) then - print *, 'PROB 1', CLM4_cdf, ThisK, ThisL, ThisM, ThisMin, ThisMax, scaled_fpar (n,nv) - endif - - scaled_fpar (n,nv) = scaled_fpar (n,nv) / (CLM4_fpar + 1.e-20) - - endif - end do NV_LOOP1 - - end do NT_LOOP1 - - para (:) = 0. ! zero out absorbed PAR summing array - parav = 0. - - if(associated(BTRANT)) btrant = 0. - if(associated(SIF)) sif = 0. - - do nz = 1,num_zon - - if(nz == 1) then - btran = btran1 - tcx = tx1 - qax = qx1 - endif - - if(nz == 2) then - btran = btran2 - tcx = tx2 - qax = qx2 - endif - - if(nz == 3) then - btran = btran3 - tcx = tx3 - qax = qx3 - endif - - do nv = 1,num_veg - elaz(:,nv) = elai(:,nv,nz) - esaz(:,nv) = esai(:,nv,nz) - ityz(:,nv) = ityp(:,nv,nz) - fvez(:,nv) = fveg(:,nv,nz) - end do - - do n = 1,NTILES - if(tp1(n) < (Tzero-0.01)) btran(n) = 0. ! no photosynthesis if ground fully frozen - end do - - call compute_rc(NTILES,nveg,TCx,QAx,T2M10D, & - TA, PS, ZTH,DRPAR,DFPAR,albdir,albdif, & - elaz,esaz,ityz,fvez,btran,fwet, & - RCx,RCxDT,RCxDQ,psnsunx,psnshax,laisunx,laishax, & - dayl_fac,co2v,dtc,dea,parzone,sifsunx,sifshax, & - lmrsunx,lmrshax,fpar_sf = scaled_fpar ) - - rc00(:,nz) = rcx(:) - rcdt(:,nz) = rcxdt(:) - rcdq(:,nz) = rcxdq(:) - - psnsun(:,:,nz) = psnsunx(:,:) - psnsha(:,:,nz) = psnshax(:,:) - laisun(:,:,nz) = laisunx(:,:) - laisha(:,:,nz) = laishax(:,:) - lmrsun(:,:,nz) = lmrsunx(:,:) - lmrsha(:,:,nz) = lmrshax(:,:) - - do nv = 1,nveg - para(:) = para(:) + parzone(:,nv)*wtzone(:,nz)*fvez(:,nv) - parav(:,nv) = parav (:,nv) + parzone(:,nv)*wtzone(:,nz) - end do - - if(associated(BTRANT)) btrant(:) = btrant(:) + btran(:)*wtzone(:,nz) - if(associated(SIF)) then - do nv = 1,nveg - sif(:) = sif(:) + wtzone(:,nz)*fvez(:,nv)*(sifsunx(:,nv)*laisunx(:,nv) + sifshax(:,nv)*laishax(:,nv)) - end do - endif - - end do - - endif - - endif DO_FS1 - - ! Below we are recycling the scaled_fpar array - from this point, it contains fpar scaled or otherwise - ! ---------------------------------------------------------------------------------------------------- - - do nv = 1,nveg - scaled_fpar (:,nv) = parav (:,nv)/ (DRPAR(:) + DFPAR(:) + 1.e-20) - end do if(associated(CNCO2)) CNCO2 = CO2V * 1e6 deallocate (co2v) @@ -7092,40 +6795,6 @@ subroutine Driver ( RC ) BGALBVR, BGALBVF, BGALBNR, BGALBNF, & ! gkw: MODIS soil background albedo ALBVR, ALBNR, ALBVF, ALBNF, MODIS_SCALE=.TRUE. ) ! instantaneous snow-free albedos on tiles - if ((SCALE_ALBFPAR == 1).OR.(SCALE_ALBFPAR == 3)) then - - if(.not.allocated (MODISVISmean )) allocate (MODISVISmean (1:NTILES)) - if(.not.allocated (MODISVISstd )) allocate (MODISVISstd (1:NTILES)) - if(.not.allocated (MODISNIRmean )) allocate (MODISNIRmean (1:NTILES)) - if(.not.allocated (MODISNIRstd )) allocate (MODISNIRstd (1:NTILES)) - if(.not.allocated (MODELFPARmean)) allocate (MODELFPARmean (1:NTILES)) - if(.not.allocated (MODELFPARstd )) allocate (MODELFPARstd (1:NTILES)) - - if(ntiles > 0) then - - call MAPL_GetResource(MAPL,VISMEANFILE , label = 'VISMEAN_FILE:' , default = 'MODISVISmean.dat' , RC=STATUS ) ; VERIFY_(STATUS) - call MAPL_GetResource(MAPL,VISSTDFILE , label = 'VISSTD_FILE:' , default = 'MODISVISstd.dat' , RC=STATUS ) ; VERIFY_(STATUS) - call MAPL_GetResource(MAPL,NIRMEANFILE , label = 'NIRMEAN_FILE:' , default = 'MODISNIRmean.dat' , RC=STATUS ) ; VERIFY_(STATUS) - call MAPL_GetResource(MAPL,NIRSTDFILE , label = 'NIRSTD_FILE:' , default = 'MODISNIRstd.dat' , RC=STATUS ) ; VERIFY_(STATUS) - - call MAPL_GetResource(MAPL,FPARMEANFILE , label = 'MODELFPARMEAN_FILE:', default = 'MODELFPARmean.dat' , RC=STATUS ) ; VERIFY_(STATUS) - call MAPL_GetResource(MAPL,FPARSTDFILE , label = 'MODELFPARSTD_FILE:' , default = 'MODELFPARstd.dat' , RC=STATUS ) ; VERIFY_(STATUS) - - call MAPL_ReadForcing(MAPL,'MODISVISmean' ,VISMEANFILE ,CURRENT_TIME,MODISVISmean ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - call MAPL_ReadForcing(MAPL,'MODISVISstd' ,VISSTDFILE ,CURRENT_TIME,MODISVISstd ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - call MAPL_ReadForcing(MAPL,'MODISNIRmean' ,NIRMEANFILE ,CURRENT_TIME,MODISNIRmean ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - call MAPL_ReadForcing(MAPL,'MODISNIRstd' ,NIRSTDFILE ,CURRENT_TIME,MODISNIRstd ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - call MAPL_ReadForcing(MAPL,'MODELFPARmean',FPARMEANFILE,CURRENT_TIME,MODELFPARmean ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - call MAPL_ReadForcing(MAPL,'MODELFPARstd' ,FPARSTDFILE ,CURRENT_TIME,MODELFPARstd ,ON_TILES=.true.,RC=STATUS) ; VERIFY_(STATUS) - - do n = 1,NTILES - ThisFPAR = (unscaled_fpar(n,1)*FVG(N,1) + unscaled_fpar(n,2)*FVG(N,2))/(FVG(N,1) + FVG(N,2) + 1.e-20) - ZFPAR = (ThisFPAR - MODELFPARmean (n)) / MODELFPARstd (n) - ALBVF(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISVISstd(n) + MODISVISmean (n))) - ALBNF(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISNIRstd(n) + MODISNIRmean (n))) - end do - endif - endif call STIEGLITZSNOW_CALC_TPSNOW(NTILES, HTSNNN(1,:), WESNN(1,:), TPSN1OUT1, FICE1) TPSN1OUT1 = TPSN1OUT1 + Tzero @@ -7142,16 +6811,6 @@ subroutine Driver ( RC ) BGALBVR, BGALBVF, BGALBNR, BGALBNF, & ! gkw: MODIS soil background albedo ALBVR_tmp, ALBNR_tmp, ALBVF_tmp, ALBNF_tmp, MODIS_SCALE=.TRUE. ) ! instantaneous snow-free albedos on tiles - if ((SCALE_ALBFPAR == 1).OR.(SCALE_ALBFPAR == 3)) then - if(ntiles > 0) then - do n = 1,NTILES - ThisFPAR = (unscaled_fpar(n,3)*FVG(N,3) + unscaled_fpar(n,4)*FVG(N,4))/(FVG(N,3) + FVG(N,4) + 1.e-20) - ZFPAR = (ThisFPAR - MODELFPARmean (n)) / MODELFPARstd (n) - ALBVF_tmp(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISVISstd(n) + MODISVISmean (n))) - ALBNF_tmp(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISNIRstd(n) + MODISNIRmean (n))) - end do - endif - endif call SNOW_ALBEDO(NTILES,N_snow, N_CONST_LAND4SNWALB, VEG2, LAI2, ZTH, & RHOFS, & @@ -7896,13 +7555,6 @@ subroutine Driver ( RC ) BGALBVR, BGALBVF, BGALBNR, BGALBNF, & ! gkw: MODIS soil background albedo ALBVR, ALBNR, ALBVF, ALBNF, MODIS_SCALE=.TRUE. ) ! instantaneous snow-free albedos on tiles - if ((SCALE_ALBFPAR == 1).OR.(SCALE_ALBFPAR == 3)) then - do n = 1,NTILES - ThisFPAR = (unscaled_fpar(n,1)*FVG(N,1) + unscaled_fpar(n,2)*FVG(N,2))/(FVG(N,1) + FVG(N,2) + 1.e-20) - ZFPAR = (ThisFPAR - MODELFPARmean (n)) / MODELFPARstd (n) - ALBVF(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISVISstd(n) + MODISVISmean (n))) - end do - endif call STIEGLITZSNOW_CALC_TPSNOW(NTILES, HTSNNN(1,:), WESNN(1,:), TPSN1OUT1, FICE1) TPSN1OUT1 = TPSN1OUT1 + Tzero @@ -7919,17 +7571,6 @@ subroutine Driver ( RC ) BGALBVR, BGALBVF, BGALBNR, BGALBNF, & ! gkw: MODIS soil background albedo ALBVR_tmp, ALBNR_tmp, ALBVF_tmp, ALBNF_tmp, MODIS_SCALE=.TRUE. ) ! instantaneous snow-free albedos on tiles - if ((SCALE_ALBFPAR == 1).OR.(SCALE_ALBFPAR == 3)) then - - do n = 1,NTILES - ThisFPAR = (unscaled_fpar(n,3)*FVG(N,3) + unscaled_fpar(n,4)*FVG(N,4))/(FVG(N,3) + FVG(N,4) + 1.e-20) - ZFPAR = (ThisFPAR - MODELFPARmean (n)) / MODELFPARstd (n) - ALBVF_tmp(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISVISstd(n) + MODISVISmean (n))) - ALBNF_tmp(n) = AMIN1 (1., AMAX1(0.001,ZFPAR * MODISNIRstd(n) + MODISNIRmean (n))) - end do - - if(allocated (MODISVISmean)) deallocate (MODISVISmean, MODISVISstd, MODISNIRmean, MODISNIRstd, MODELFPARmean, MODELFPARstd) - endif call SNOW_ALBEDO(NTILES,N_snow, N_CONST_LAND4SNWALB, VEG2, LAI2, ZTH, & RHOFS, & @@ -8281,11 +7922,7 @@ subroutine Driver ( RC ) deallocate( car1 ) deallocate( car2 ) deallocate( car4 ) - deallocate( parzone ) deallocate( para ) - deallocate( parav ) - deallocate (scaled_fpar) - deallocate (UNscaled_fpar) deallocate( totwat ) deallocate( nfire ) deallocate(som_closs) @@ -8342,6 +7979,7 @@ subroutine Driver ( RC ) deallocate( psnsunx ) deallocate( psnshax ) deallocate( sifsunx ) + deallocate( parzone ) deallocate( sifshax ) deallocate( laisunx ) deallocate( laishax ) @@ -8374,128 +8012,109 @@ subroutine Driver ( RC ) end subroutine Driver - ! ----------------- routines for CDF scaling ------------------- - - REAL FUNCTION cdf2fpar (cdf, k,l, m, m1, m2) - - REAL, intent (in) :: cdf, k,l,m, m1, m2 - REAL :: x, ThisCDF, ThisFPAR - integer, parameter :: nBINS = 40 - - x = real (nBINS) - ThisCDF = 1. - - do while (ThisCDF >= cdf) - ThisFPAR = 1. - (real(nbins)-x)/real(nbins) - 1./2./real(nbins) - ThisCDF = K * betai (L, M, ThisFPAR) - x = x - 1. - if(x == 0) exit - end do - - cdf2fpar = ThisFPAR * m2 + m1 - if(cdf2fpar > m2) cdf2fpar = m2 - if(cdf2fpar < m1) cdf2fpar = m1 - return - - END FUNCTION cdf2fpar - - ! --------------------------------------------------------- - FUNCTION betai(a,b,x) - REAL betai,a,b,x - REAL bt - !external gammln - - if (x < 0.0125) x = 0.0125 - if (x > 0.9875) x = 0.9875 - - if(x.lt.0..or.x.gt.1.)print *, 'bad argument x in betai',x - if(x.lt.0..or.x.gt.1.)stop - if(x.eq.0..or.x.eq.1.)then - bt=0. - else - bt=exp(gammln(a+b)-gammln(a)-gammln(b) & - +a*log(x)+b*log(1.-x)) - endif - - if(x.lt.(a+1.)/(a+b+2.))then - betai=bt*betacf(a,b,x)/a - return - else - betai=1.-bt*betacf(b,a,1.-x)/b - return - endif - - END FUNCTION betai - - ! ------------------------------------------------------- - - FUNCTION betacf(a,b,x) - - INTEGER MAXIT - REAL betacf,a,b,x,EPS,FPMIN - PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30) - INTEGER m,m2 - REAL aa,c,d,del,h,qab,qam,qap - - qab=a+b - qap=a+1. - qam=a-1. - c=1. - d=1.-qab*x/qap - - if(abs(d).lt.FPMIN)d=FPMIN - d=1./d - h=d - do m=1,MAXIT - m2=2*m - aa=m*(b-m)*x/((qam+m2)*(a+m2)) - d=1.+aa*d - if(abs(d).lt.FPMIN)d=FPMIN - c=1.+aa/c - if(abs(c).lt.FPMIN)c=FPMIN - d=1./d - h=h*d*c - aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2)) - d=1.+aa*d - if(abs(d).lt.FPMIN)d=FPMIN - c=1.+aa/c - if(abs(c).lt.FPMIN)c=FPMIN - d=1./d - del=d*c - h=h*del - if(abs(del-1.).lt.EPS)exit - enddo - betacf=h - return - - END FUNCTION betacf - - ! -------------------------------------------------------------- +! Commented out functions betai(), betacf(), and gammln(). +! These functions are not used and were reproduced identically in +! GEOS_CatchCNCLM40GridComp.F90 and in GEOS_CatchCNCLM45GridComp.F90. +! Another copy was in GEOScatchCN_GridComp/utils/math_routines.F90 but +! there function betai() was missing the restriction 0.0125 0.9875) x = 0.9875 +!! +!! if(x.lt.0..or.x.gt.1.)print *, 'bad argument x in betai',x +!! if(x.lt.0..or.x.gt.1.)stop +!! if(x.eq.0..or.x.eq.1.)then +!! bt=0. +!! else +!! bt=exp(gammln(a+b)-gammln(a)-gammln(b) & +!! +a*log(x)+b*log(1.-x)) +!! endif +!! +!! if(x.lt.(a+1.)/(a+b+2.))then +!! betai=bt*betacf(a,b,x)/a +!! return +!! else +!! betai=1.-bt*betacf(b,a,1.-x)/b +!! return +!! endif +!! +!! END FUNCTION betai +!! +!! ! ------------------------------------------------------- +!! +!! FUNCTION betacf(a,b,x) +!! +!! INTEGER MAXIT +!! REAL betacf,a,b,x,EPS,FPMIN +!! PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30) +!! INTEGER m,m2 +!! REAL aa,c,d,del,h,qab,qam,qap +!! +!! qab=a+b +!! qap=a+1. +!! qam=a-1. +!! c=1. +!! d=1.-qab*x/qap +!! +!! if(abs(d).lt.FPMIN)d=FPMIN +!! d=1./d +!! h=d +!! do m=1,MAXIT +!! m2=2*m +!! aa=m*(b-m)*x/((qam+m2)*(a+m2)) +!! d=1.+aa*d +!! if(abs(d).lt.FPMIN)d=FPMIN +!! c=1.+aa/c +!! if(abs(c).lt.FPMIN)c=FPMIN +!! d=1./d +!! h=h*d*c +!! aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2)) +!! d=1.+aa*d +!! if(abs(d).lt.FPMIN)d=FPMIN +!! c=1.+aa/c +!! if(abs(c).lt.FPMIN)c=FPMIN +!! d=1./d +!! del=d*c +!! h=h*del +!! if(abs(del-1.).lt.EPS)exit +!! enddo +!! betacf=h +!! return +!! +!! END FUNCTION betacf +!! +!! ! -------------------------------------------------------------- +!! +!! FUNCTION gammln(xx) +!! +!! REAL gammln,xx +!! INTEGER j +!! DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) +!! +!! SAVE cof,stp +!! DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, & +!! 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, & +!! -.5395239384953d-5,2.5066282746310005d0/ +!! x=xx +!! y=x +!! tmp=x+5.5d0 +!! tmp=(x+0.5d0)*log(tmp)-tmp +!! ser=1.000000000190015d0 +!! do j=1,6 +!! y=y+1.d0 +!! ser=ser+cof(j)/y +!! enddo +!! gammln=tmp+log(stp*ser/x) +!! return +!! +!! END FUNCTION gammln ! -------------------------------------------------------------- diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/utils/compute_FPAR_CDF_M09.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/utils/compute_FPAR_CDF_M09.F90 deleted file mode 100755 index 2e1314a3d..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/utils/compute_FPAR_CDF_M09.F90 +++ /dev/null @@ -1,801 +0,0 @@ -#include "Raster.h" - -PROGRAM comp_FPAR_CDF - - use math_routines - use MAPL_SortMod - use date_time_util, ONLY: & - date_time_type, augment_date_time - use ieee_arithmetic, only: isnan => ieee_is_nan - IMPLICIT NONE - - INCLUDE 'netcdf.inc' - INCLUDE 'mpif.h' - - integer :: comm_rank, comm_size, error, info, DOY,NCInID,NCOutID, NCOutID2, n, iargc, STATUS - integer :: stat(MPI_STATUS_SIZE), NTILES, NCatch, NDATA1, NDATA2,N_VIS_DATA, N_NIR_DATA - character*400 :: arg(2) - REAL, PARAMETER :: TILESIZE = 10., TILEINT = 2. - INTEGER, PARAMETER :: NPFT = 19, NOCTAD = 46, NSETS = 2 ,MXCNT = 500000 - logical, parameter :: onebin = .true. - INTEGER, PARAMETER :: YearB = 2003, YearE = 2016 - character*400, PARAMETER :: & - BCSDIR = 'SMAP_EASEv2_M09/', & - EXPDIR = '/discover/nobackup/fzeng/Catchment/SMAP_EASEv2_M09/e0004s_wet2/output/SMAP_EASEv2_M09_GLOBAL/', & - EXNAME = 'e0004s_wet2', & - OUTFIL = 'global_alb_mu_std/FPAR_CDF_Params-M09', & - LOGFIL = 'global_alb_mu_std/FPAR_CDF_Params-M09_log.', & - GFILE = 'SMAP_EASEv2_M09_3856x1624' - - logical :: file_exists, put_aux = .true. -! real, dimension (4) :: limits = (/20., -130., 60., -60./) - real, dimension (4) :: limits = (/-90., -180., 90., 180./) - real, dimension (NBINS) :: MODIS_BINS, CLM4_BINS - INTEGER :: I_INDEX(10),JM, IM, NC, NT, year, day, maxcat, ThisTile, req - type(date_time_type), dimension (YearE - YearB + 1) :: octad_time - CHARACTER*8 :: YYYYMMDD - CHARACTER*6 :: YYYYMM - CHARACTER*4 :: YYYY - CHARACTER*2 :: MM, DD,MMR, TSLICE - INTEGER :: i, j,k, pf, i1,i2,i3,i4, f1,f2,f3,f4, ND - INTEGER, DIMENSION (:,:), allocatable :: veg_index, catchs_all - integer, allocatable, dimension (:) :: ldas2bcs, Catchs, NCATCH_ALL - real :: yr,mn,dy,dum,yr1,mn1,dy1, lw, up, db, mean, std, skew, & - minv, maxv, minv1, maxv1, minv2, maxv2, VISmean, FPARmean, FPARstd, NIRmean, NIRstd, VISstd, r2, var1, var2 - real, dimension (:,:), allocatable :: modis_fpar, clm4_fpar, modis_visdf, modis_nirdf - real, dimension (:), allocatable :: modis_cdf, clm4_cdf, data_read, data_save, modis_hist, clm4_hist - real, dimension (:), allocatable :: vis_this, nir_this, fpar_this, est_this - real (kind=8), dimension(NBINS) :: dbins,dcdf - real (kind=8), dimension(3) :: modis_param, clm4_param - real , dimension(NBINS*2 + 13 + 6) :: tmp_real - character*300 :: tmpstring - - call MPI_Init(error) - call MPI_COMM_Size(MPI_COMM_WORLD,comm_size,error) - call MPI_COMM_Rank(MPI_COMM_WORLD,comm_rank,error) - - call MPI_Info_create(info, error) - call MPI_Info_set(info, "romio_cb_read", "automatic", error) - - write (TSLICE ,'(i2.2)') comm_rank + 1 - open (99,file=trim(logfil)//TSLICE, form ='formatted', action='write', status= 'unknown') - - ! STEP 1 check/create CDF params files - - inquire(file=trim(OUTFIL)//'.nc4', exist = file_exists) - - if(.not.file_exists) then - if(comm_rank == 0) then - call create_CDF_ParamFile - WRITE (99,*)'CREATED CDF PARAM FILE : ', trim(OUTFIL)//'.nc4' - endif - endif - - call MPI_BARRIER( MPI_COMM_WORLD, error) - - ! READ # OF 10x10 MODIS TILES AND SMAP_TILE_IDs THAT CONTRIBUTE TO EACH MODIS TILE - - STATUS = NF_OPEN (trim(OUTFIL)//'_aux.nc4', NF_NOWRITE,NCInID) ; VERIFY_(STATUS) - STATUS = NF_INQ_DIMID (NCInID, 'tile10D', K); VERIFY_(STATUS) - CALL HANDLE_ERR(STATUS, 'INQ_DIM') - STATUS = NF_INQ_DIMLEN (NCInID, K, NTILES); VERIFY_(STATUS) - CALL HANDLE_ERR(STATUS, 'DIMLEN_NTILES') - WRITE (99,*)'NOF 10D CELLS : ', NTILES - - allocate (NCatch_all (1:NTILES)) - allocate (Catchs_all (1:16000,1:NTILES)) - - STATUS = NF_GET_VARA_INT (NCInID, VarID(NCInID,'nSMAP'), (/1/), (/NTILES/), NCatch_ALL); VERIFY_(STATUS) - STATUS = NF_GET_VARA_INT (NCInID, VarID(NCInID,'SMAPID'), (/1, 1/), (/16000,NTILES/), Catchs_ALL); VERIFY_(STATUS) - status = NF_CLOSE (NCInID) - - call MPI_BARRIER( MPI_COMM_WORLD, error) - - if(comm_rank == 0) then - - ! ROOT PROCESSOR OPENS FILES TO UPDATE CDF PARAMS - - STATUS = NF_OPEN (trim(OUTFIL)//'.nc4' ,NF_WRITE,NCOutID ) - VERIFY_(STATUS) - STATUS = NF_OPEN (trim(OUTFIL)//'_aux.nc4',NF_WRITE,NCOutID2) - VERIFY_(STATUS) - - endif - - ! read maxcat, LDASsa tile order, veg types and define binvals - - open (10,file=trim(BCSDIR)//'clsm/catchment.def', status='old',action='read', & - form='formatted') - - read (10,*) maxcat - - close (10, status ='keep') - - allocate (veg_index (1: NPFT, 1: MAXCAT)) - allocate (ldas2bcs (1: maxcat)) - allocate (data_read (1: maxcat)) - allocate (data_save (1: maxcat)) - allocate (modis_fpar (1: maxcat, yearE - yearB + 1)) - allocate (clm4_fpar (1: maxcat, yearE - yearB + 1)) - allocate (modis_visdf (1: maxcat, yearE - yearB)) - allocate (modis_nirdf (1: maxcat, yearE - yearB)) - - clm4_fpar = 0. - modis_fpar = 0. - data_save = 0. - veg_index = -9999 - - open (10,file =trim(EXPDIR)//'rc_out/'//trim(EXNAME)//'.ldas_tilecoord.bin',status='old',form='unformatted',convert='big_endian') - read (10) i - read (10) LDAS2BCS - close (10, status = 'keep') - - open (10,file=trim(BCSDIR)//'clsm/CLM_veg_typs_fracs', status='old',action='read', & - form='formatted') - - do i = 1, maxcat - - read (10,*) pf, pf, i1,i2,i3,i4, f1,f2,f3,f4 - - if(f1 >= f2) then - veg_index (i1,i) = i - else - veg_index (i2,i) = i - endif - - end do - - close (10, status ='keep') - -! OCTAD_LOOP : DO NT = 1, NOCTAD - - NT = comm_rank + 1 - ND = 8 - - ! BEGIN READING MODIS FPAR - - MODIS_LOOP : DO year = YearB, YearE - - write (YYYY ,'(i4.4)') year - - WRITE (99,*)'FPAR, VISDF and NIRDF FILES : ' - WRITE (99,*) YYYY//'//fpar.dat' - - open (10, file = trim(BCSDIR)//'MODIS6/'//YYYY//'//fpar.dat', form = 'unformatted', action = 'read') - - if(year < yearE) then - open (11, file = trim(BCSDIR)//'MODIS6/'//YYYY//'//visdf.dat', form = 'unformatted', action = 'read') - open (12, file = trim(BCSDIR)//'MODIS6/'//YYYY//'//nirdf.dat', form = 'unformatted', action = 'read') - WRITE (99,*) YYYY//'//visdf.dat' - WRITE (99,*) YYYY//'//nirdf.dat' - WRITE (99,*) ' ' - endif - - do k = 0, nt ! your processor rank - - read (10) yr,mn,dy,dum,dum,dum,yr1,mn1,dy1 - read (10) modis_fpar (:, year - yearB + 1) - - if ((k > 0) .and.(year < yearE)) then - read (11) modis_visdf (:, year - yearB + 1) - read (12) modis_nirdf (:, year - yearB + 1) - endif - - IF (k == NT) WRITE (99,*) 'PROCESSING TIME SLICE : ', yr,mn,dy,yr1,mn1,dy1 - - end do - - close (10, status = 'keep') - - if(year < yearE) then - close (11, status = 'keep') - close (12, status = 'keep') - endif - - END DO MODIS_LOOP - - ! END READING MODIS FPAR and BEGIN READING CLM4 FPAR - - WRITE (99,*) ' ' - WRITE (99,*) 'READING CLM4 FPAR: ' - - CLM4_LOOP : DO year = YearB, YearE - - octad_time(year - yearB + 1)%year = year - 1 - octad_time(year - yearB + 1)%month = 12 - octad_time(year - yearB + 1)%day = 31 - octad_time(year - yearB + 1)%hour = 0 - octad_time(year - yearB + 1)%min = 0 - octad_time(year - yearB + 1)%sec = 0 - - do k = 1, nt - - if((K == NT).and.(NT == 46)) ND = 5 - - DO day = 1,ND - call augment_date_time( 86400, octad_time(year - yearB + 1)) - - if(k == nt) then - write (YYYY, '(i4.4)') octad_time(year - yearB + 1)%year - write (MM , '(i2.2)') octad_time(year - yearB + 1)%month - write (DD , '(i2.2)') octad_time(year - yearB + 1)%day - YYYYMMDD = YYYY//MM//DD - WRITE (99,*) trim(EXNAME)//'.ens_avg.ldas_tile_daily_out.'//YYYYMMDD//'.bin' - open (10, file = trim(EXPDIR)//'cat/ens_avg/Y'//YYYY//'/M'//MM//'/'// & - trim(EXNAME)//'.ens_avg.ldas_tile_daily_out.'//YYYYMMDD//'.bin', & - form = 'unformatted', convert='big_endian', action = 'read') - - do n = 1,3 - read (10) data_read - if(n == 2) data_save = data_read - end do - - clm4_fpar (:,year - yearB + 1) = clm4_fpar (:,year - yearB + 1) + data_save / (data_read + 1.e-20)/ real (ND) - - close (10, status = 'keep') - - endif - end do - - ND = 8 - - end do - - ! reoder to the order of BCs - - data_read = clm4_fpar (:,year - yearB + 1) - - do n = 1, maxcat - clm4_fpar (LDAS2BCS(n),year - yearB + 1) = data_read (n) - end do - - END DO CLM4_LOOP - - ! Now compute CDFs - ! ---------------- - - ! loop through tiles - - allocate (modis_hist (1:MXCNT)) - allocate (clm4_hist (1:MXCNT)) - allocate (vis_this (1:MXCNT)) - allocate (nir_this (1:MXCNT)) - allocate (fpar_this (1:MXCNT)) - allocate (est_this (1:MXCNT)) - allocate (modis_cdf (1:NBINS)) - allocate (clm4_cdf (1:NBINS)) - - TILE_LOOP : DO ThisTile = 1,NTILES - - NCATCH = NCATCH_ALL (ThisTile) - - if(NCatch < 1) then - WRITE (99,*) 'nSMAP problem ', NCatch, ThisTile - endif - - allocate (Catchs (1: NCatch)) - ! STATUS = NF_GET_VARA_INT (NCInID, VarID(NCInID,'SMAPID'), (/1, ThisTile/), (/NCatch,1/), Catchs) - Catchs (1: NCatch) = Catchs_ALL (1: NCatch, ThisTile) - - PFT_LOOP: DO k =1, NPFT - - modis_hist = 0. - clm4_hist = 0. - modis_cdf = 0. - clm4_cdf = 0. - NDATA1 = 1 - NDATA2 = 1 - N_VIS_DATA = 1 - N_NIR_DATA = 1 - MODIS_BINS = 0. - CLM4_BINS = 0. - - DO year = YearB, YearE - DO N = 1, NCatch - if((veg_index(k,Catchs(n)) > 0).and.(modis_fpar (veg_index(k,Catchs(n)), year - yearB + 1) > 0. )) then - modis_hist (NDATA1) = modis_fpar (veg_index(k,Catchs(n)), year - yearB + 1) - NDATA1 = NDATA1 + 1 - endif - - if((veg_index(k,Catchs(n)) > 0).and.(clm4_fpar (veg_index(k,Catchs(n)),year - yearB + 1) > 0. )) then - clm4_hist (NDATA2) = clm4_fpar (veg_index(k,Catchs(n)),year - yearB + 1) - NDATA2 = NDATA2 + 1 - endif - - if (year < yearE) then - if((veg_index(k,Catchs(n)) > 0).and.(modis_fpar (veg_index(k,Catchs(n)), year - yearB + 1) >= 0. ).and.( modis_visdf(veg_index(k,Catchs(n)), year - yearB + 1) > 0.) & - .and.( modis_nirdf(veg_index(k,Catchs(n)), year - yearB + 1) > 0.)) then - vis_this (N_VIS_DATA) = modis_visdf (veg_index(k,Catchs(n)), year - yearB + 1) - fpar_this(N_VIS_DATA) = modis_fpar (veg_index(k,Catchs(n)), year - yearB + 1) - nir_this (N_VIS_DATA) = modis_nirdf (veg_index(k,Catchs(n)), year - yearB + 1) - N_VIS_DATA = N_VIS_DATA + 1 - endif - -! if((veg_index(k,Catchs(n)) > 0).and.(modis_fpar (veg_index(k,Catchs(n)), year - yearB + 1) >= 0. ).and.( modis_nirdf(veg_index(k,Catchs(n)), year - yearB + 1) > 0.)) then -! nir_this (N_NIR_DATA) = modis_nirdf (veg_index(k,Catchs(n)), year - yearB + 1) -! N_NIR_DATA = N_NIR_DATA + 1 -! endif - endif - - if((ndata1 > MXCNT).or.(ndata2 > MXCNT).or.(N_VIS_DATA > MXCNT)) then - WRITE (99,*) 'NDATA1 or NDATA2 exceeded ',ndata1, ndata2, N_VIS_DATA, N_NIR_DATA - stop - endif - - END DO - END DO - - NDATA1 = NDATA1 - 1 - NDATA2 = NDATA2 - 1 - N_NIR_DATA = N_NIR_DATA - 1 - N_VIS_DATA = N_VIS_DATA - 1 - - MINV1 = -9999. - MAXV1 = -9999. - MINV2 = -9999. - MAXV2 = -9999. - - ! curve fitting - - modis_param = -9999. - clm4_param = -9999. - FPARmean = -9999. - FPARstd = -9999. - - if((NDATA1 > 10).and.(NDATA2 > 10)) then - - WRITE (99,*) '# of SMAP DATA CELLS ',ThisTile, K, NDATA1,NDATA2, N_NIR_DATA,N_VIS_DATA - - maxv = MAXVAL ((/MAXVAL(modis_hist (1: NDATA1)),MAXVAL(clm4_hist (1: NDATA2))/)) - minv = MINVAL ((/MINVAL(modis_hist (1: NDATA1)),MINVAL(clm4_hist (1: NDATA2))/)) - - if (maxv > minv) then - - modis_hist (1: NDATA1) = (modis_hist (1: NDATA1) - minv)/maxv - call prob_den_func (ndata1, modis_hist (1: NDATA1) , modis_cdf, MODIS_bins, lwval= 0.,upval = 1.) - - dbins = MODIS_bins - dcdf = modis_cdf - modis_param (1) = .5 - modis_param (2) = .5 - modis_param (3) = 0.9 - call optimiz (nbins,dbins,dcdf,modis_param) - - FPARmean = SUM (clm4_hist (1:NDATA2)) / real (NDATA2) - var1 = 0. - - do i = 1,NDATA2 - var1 = var1 + (clm4_hist(i) - FPARmean)*(clm4_hist(i) - FPARmean) - end do - - FPARstd = sqrt (var1/real(NDATA2 - 1)) - - clm4_hist (1: NDATA2) = (clm4_hist (1: NDATA2) - minv)/maxv - call prob_den_func (ndata2, clm4_hist (1: NDATA2) , clm4_cdf, CLM4_bins, lwval= 0.,upval = 1.) - - dbins = CLM4_bins - dcdf = clm4_cdf - clm4_param (1) = .5 - clm4_param (2) = .5 - clm4_param (3) = 0.9 - call optimiz (NBINS,dbins,dcdf,clm4_param) - endif - endif - - ! albedo parameters - - VISmean = -9999. - NIRmean = -9999. - VISstd = -9999. - NIRstd = -9999. - - if (N_VIS_DATA > 12) then - NIRmean = SUM (nir_this (1:N_VIS_DATA)) / real (N_VIS_DATA) - VISmean = SUM (vis_this (1:N_VIS_DATA)) / real (N_VIS_DATA) - - var1 = 0. - var2 = 0. - - do i = 1,N_VIS_DATA - var1 = var1 + (vis_this(i) - VISmean)*(vis_this(i) - VISmean) - var2 = var2 + (nir_this(i) - NIRmean)*(nir_this(i) - NIRmean) - end do - - VISstd = sqrt (var1/real(N_VIS_DATA - 1)) - NIRstd = sqrt (var2/real(N_VIS_DATA - 1)) - - endif - - do i = 1, comm_size - - tmp_real = -9999. - - if((I == 1).and.(comm_rank == 0)) then - - if(put_aux) then - STATUS = NF_PUT_VARA_REAL(NCOutID2,VARID(NCOutID2,'CDF' ),(/1,ThisTile,k,NT,1/), (/NBINS,1,1,1,1/), modis_cdf );VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID2,VARID(NCOutID2,'CDF' ),(/1,ThisTile,k,NT,2/), (/NBINS,1,1,1,1/), clm4_cdf );VERIFY_(STATUS) - endif - - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'Kappa' ),(/ThisTile,k,NT,1/), (/1,1,1,1/), REAL (modis_param(3))); VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'Lambda'),(/ThisTile,k,NT,1/), (/1,1,1,1/), REAL (modis_param(1))); VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'Mu' ),(/ThisTile,k,NT,1/), (/1,1,1,1/), REAL (modis_param(2))) ; VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'Kappa' ),(/ThisTile,k,NT,2/), (/1,1,1,1/), REAL (clm4_param(3))); VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'Lambda'),(/ThisTile,k,NT,2/), (/1,1,1,1/), REAL (clm4_param(1))); VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'Mu' ),(/ThisTile,k,NT,2/), (/1,1,1,1/), REAL (clm4_param(2))) ; VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MinVal'),(/ThisTile,k,NT/) , (/1,1,1/), MINV); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MaxVal'),(/ThisTile,k,NT/) , (/1,1,1/), MAXV); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODISVISmean' ),(/ThisTile,k,NT/), (/1,1,1/),VISmean ); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODISNIRmean' ),(/ThisTile,k,NT/), (/1,1,1/),NIRmean ); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODISVISstd' ),(/ThisTile,k,NT/), (/1,1,1/),VISstd ); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODISNIRstd' ),(/ThisTile,k,NT/), (/1,1,1/),NIRstd ); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODELFPARmean'),(/ThisTile,k,NT/), (/1,1,1/),FPARmean); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODELFPARstd' ),(/ThisTile,k,NT/), (/1,1,1/),FPARstd ); VERIFY_(STATUS);VERIFY_(STATUS) - WRITE (99,*) 'Writing Out OCTAD', ThisTile,k,NT - - else if (I > 1) then - - if(I-1 == comm_rank) then - ! print *, comm_rank, 'sending', ThisTile, k,nt - tmp_real (1) = real (ThisTile) - tmp_real (2) = real(k) - tmp_real (3) = real(NT) - tmp_real (4) = REAL (modis_param(3)) - tmp_real (5) = REAL (modis_param(1)) - tmp_real (6) = REAL (modis_param(2)) - tmp_real (7) = REAL (clm4_param(3)) - tmp_real (8) = REAL (clm4_param(1)) - tmp_real (9) = REAL (clm4_param(2)) - tmp_real(10) = MINV - tmp_real(11) = MINV - tmp_real(12) = MAXV - tmp_real(13) = MAXV - tmp_real(14) = VISmean - tmp_real(15) = NIRmean - tmp_real(16) = VISstd - tmp_real(17) = NIRstd - tmp_real(18) = FPARmean - tmp_real(19) = FPARstd - - NC = 19 - n = NC + NBINS - tmp_real (NC+1: N) = modis_cdf(:) - NC = n - n = NC + NBINS - tmp_real (NC+1: N) = clm4_cdf(:) - NC = n - - call MPI_ISend(tmp_real ,2*NBINS + 19,MPI_real,0,999,MPI_COMM_WORLD,req,status) - call MPI_WAIT (req,MPI_STATUS_IGNORE,status) - - else if (comm_rank == 0) then - - call MPI_RECV(tmp_real,2*NBINS + 19,MPI_real,I-1,999,MPI_COMM_WORLD,MPI_STATUS_IGNORE,status) - - IM = NINT (tmp_real (1)) - JM = NINT (tmp_real (2)) - I1 = NINT (tmp_real (3)) - - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID,'Kappa' ),(/IM,JM,I1,1/), (/1,1,1,1/), tmp_real (4)); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID,'Lambda'),(/IM,JM,I1,1/), (/1,1,1,1/), tmp_real (5)); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID,'Mu' ),(/IM,JM,I1,1/), (/1,1,1,1/), tmp_real (6)) ; VERIFY_(STATUS) ;VERIFY_(STATUS) - - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID,'Kappa' ),(/IM,JM,I1,2/), (/1,1,1,1/), tmp_real (7)); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID,'Lambda'),(/IM,JM,I1,2/), (/1,1,1,1/), tmp_real (8)); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID,'Mu' ),(/IM,JM,I1,2/), (/1,1,1,1/), tmp_real (9)); VERIFY_(STATUS);VERIFY_(STATUS) - - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID,'MinVal'),(/IM,JM,I1/) , (/1,1,1/) , tmp_real(10)) ; VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID,'MaxVal'),(/IM,JM,I1/) , (/1,1,1/) , tmp_real(12)) ; VERIFY_(STATUS) ;VERIFY_(STATUS) - - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODISVISmean' ),(/IM,JM,I1/), (/1,1,1/), tmp_real(14)); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODISNIRmean' ),(/IM,JM,I1/), (/1,1,1/), tmp_real(15)); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODISVISstd' ),(/IM,JM,I1/), (/1,1,1/), tmp_real(16)); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODISNIRstd' ),(/IM,JM,I1/), (/1,1,1/), tmp_real(17)); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODELFPARmean'),(/IM,JM,I1/), (/1,1,1/), tmp_real(18)); VERIFY_(STATUS);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID, VARID(NCOutID ,'MODELFPARstd' ),(/IM,JM,I1/), (/1,1,1/), tmp_real(19)); VERIFY_(STATUS);VERIFY_(STATUS) - - NC = 19 - n = NC + NBINS - if(put_aux) STATUS = NF_PUT_VARA_REAL(NCOutID2,VARID(NCOutID2,'CDF' ),(/1,IM,JM,I1,1/), (/NBINS,1,1,1,1/), tmp_real (nc+1: N)); VERIFY_(STATUS) - NC = n - n = NC + NBINS - if(put_aux) STATUS = NF_PUT_VARA_REAL(NCOutID2,VARID(NCOutID2,'CDF' ),(/1,IM,JM,I1,2/), (/NBINS,1,1,1,1/), tmp_real (nc+1: N)); VERIFY_(STATUS) - NC = n -! n = NC + NBINS -! if(put_aux) STATUS = NF_PUT_VARA_REAL(NCOutID2,VARID(NCOutID2,'BINS' ),(/1,IM,JM,I1/) , (/NBINS,1,1,1/) , tmp_real (nc+1: N)); VERIFY_(STATUS) -! NC = n -! n = NC + NBINS -! if(put_aux) STATUS = NF_PUT_VARA_REAL(NCOutID2,VARID(NCOutID2,'BINS' ),(/1,IM,JM,I1,2/), (/NBINS,1,1,1,1/), tmp_real (nc+1: N)); VERIFY_(STATUS) - - WRITE (99,*) ' RECEIVED Writing Out OCTAD', I-1,IM,JM,I1 - endif - endif - end do - - END DO PFT_LOOP - - deallocate (Catchs) - - END DO TILE_LOOP - - ! END DO OCTAD_LOOP - if(comm_rank == 0) then - status = NF_CLOSE (NCOutID ) - status = NF_CLOSE (NCOutID2) - endif - close (99, status = 'keep') - - call MPI_BARRIER( MPI_COMM_WORLD, error) - call MPI_Finalize(STATUS) - -STOP - -CONTAINS - -!________________________________________________________________________ - - -SUBROUTINE create_CDF_ParamFile - - implicit none - - integer :: NCFOutID, NCFOutID2, status, pid, tid, did, lid, bid, vid,CID - integer :: i,j, k, n, maxcat,ii,jj, tile_count, nplus, nb, NB_max - integer :: nc_rst = 43200, nr_rst = 21600, DIJ, MXT = 16000 - integer :: PID2, TID2, DID2, LID2,BID2, CID2 - real :: dxy, lw, up, db - real, allocatable, dimension (:) :: abins - integer, allocatable, target, dimension (:,:) :: tile_id - integer, pointer, dimension (:,:) :: tile_id_box - integer, allocatable, dimension (:) :: tile_id_vec, bcs2ldas - integer, allocatable, dimension (:) :: density, loc_int - integer, allocatable, dimension (:) :: loc_val - logical, allocatable, dimension (:) :: unq_mask - character (22) :: time_stamp - integer, dimension(8) :: date_time_values - - - status = NF_CREATE (trim(OUTFIL)//'.nc4' , NF_NETCDF4, NCFOutID );VERIFY_(STATUS) - status = NF_CREATE (trim(OUTFIL)//'_aux.nc4', NF_NETCDF4, NCFOutID2);VERIFY_(STATUS) - ! Define Dimensions - - status = NF_DEF_DIM(NCFOutID, 'pft' , NPFT , PID);VERIFY_(STATUS) - status = NF_DEF_DIM(NCFOutID, 'octad' , NOCTAD, TID);VERIFY_(STATUS) - status = NF_DEF_DIM(NCFOutID, 'data' , NSETS , DID);VERIFY_(STATUS) - status = NF_DEF_DIM(NCFOutID, 'tile10D',NF_UNLIMITED,LID);VERIFY_(STATUS) - - - status = NF_DEF_DIM(NCFOutID2, 'pft' , NPFT , PID2);VERIFY_(STATUS) - status = NF_DEF_DIM(NCFOutID2, 'octad' , NOCTAD, TID2);VERIFY_(STATUS) - status = NF_DEF_DIM(NCFOutID2, 'data' , NSETS , DID2);VERIFY_(STATUS) - status = NF_DEF_DIM(NCFOutID2, 'tile10D',NF_UNLIMITED,LID2);VERIFY_(STATUS) - status = NF_DEF_DIM(NCFOutID2, 'bin' , nbins, BID2);VERIFY_(STATUS) - status = NF_DEF_DIM(NCFOutID2, 'nCELLS' , MXT, CID2);VERIFY_(STATUS) - - ! Define variables - - status = NF_DEF_VAR(NCFOutID, 'lon' , NF_FLOAT ,1 ,(/LID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'lat' , NF_FLOAT ,1 ,(/LID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'Kappa' , NF_FLOAT ,4 ,(/LID,PID,TID,DID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'Lambda' , NF_FLOAT ,4 ,(/LID,PID,TID,DID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'Mu' , NF_FLOAT ,4 ,(/LID,PID,TID,DID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MinVal' , NF_FLOAT ,3 ,(/LID,PID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MaxVal' , NF_FLOAT ,3 ,(/LID,PID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODISVISmean', NF_FLOAT ,3 ,(/LID,PID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODISNIRmean', NF_FLOAT ,3 ,(/LID,PID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODISVISstd' , NF_FLOAT ,3 ,(/LID,PID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODISNIRstd' , NF_FLOAT ,3 ,(/LID,PID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODELFPARmean',NF_FLOAT ,3 ,(/LID,PID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODELFPARstd' ,NF_FLOAT ,3 ,(/LID,PID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID2, 'lon' , NF_FLOAT ,1 ,(/LID2/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID2, 'lat' , NF_FLOAT ,1 ,(/LID2/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID2, 'nSMAP' , NF_SHORT ,1 ,(/LID2/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID2, 'BINS' , NF_FLOAT ,1 ,(/BID2/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID2, 'SMAPID', NF_INT ,2 ,(/CID2, LID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID2, 'CDF' , NF_FLOAT ,5 ,(/BID2,LID2,PID2,TID2,DID2/), vid);VERIFY_(STATUS) - -! Global attributes -! - call date_and_time(VALUES=date_time_values) - - write (time_stamp,'(i4.4,a1,i2.2,a1,i2.2,1x,a2,1x,i2.2,a1,i2.2,a1,i2.2)') & - date_time_values(1),'-',date_time_values(2),'-',date_time_values(3),'at', & - date_time_values(5),':',date_time_values(6),':',date_time_values(7) - - status = NF_PUT_ATT_TEXT(NCFOutID, NF_GLOBAL, 'CreatedBy', LEN_TRIM("Sarith Mahanama"), & - trim("Sarith Mahanama")) - status = NF_PUT_ATT_TEXT(NCFOutID, NF_GLOBAL, 'Date' , LEN_TRIM(time_stamp),trim(time_stamp)) - - status = NF_PUT_ATT_TEXT(NCFOutID2, NF_GLOBAL, 'CreatedBy', LEN_TRIM("Sarith Mahanama"), & - trim("Sarith Mahanama")) - status = NF_PUT_ATT_TEXT(NCFOutID2, NF_GLOBAL, 'Date' , LEN_TRIM(time_stamp),trim(time_stamp)) - - status = NF_ENDDEF(NCFOutID ) - status = NF_ENDDEF(NCFOutID2) - - allocate (abins (1:nbins)) - - lw = 0. - up = 1. - db = (up - lw)/real(nbins) - do n = 1, nbins - abins (n) = lw + real(n)*db - db/2. - end do - - STATUS = NF_PUT_VARA_REAL(NCFOutID2,VARID(NCFOutID2, 'BINS'),(/1/) , (/NBINS/) , abins); VERIFY_(STATUS) - - ! create TILESIZE x TILESIZE tiles at TILEINT - - dij = nint (TILESIZE * nc_rst/360) - dxy = 360./nc_rst - - ! read maxcat from catchment.def - - open (10,file=trim(BCSDIR)//'clsm/catchment.def', status='old',action='read', & - form='formatted') - - read (10,*) maxcat - - close (10, status ='keep') - - ! read tilecoord for tile order in LDASsa - - open (10,file =trim(EXPDIR)//'rc_out/'//trim(EXNAME)//'.ldas_tilecoord.bin',status='old',form='unformatted',convert='big_endian') - read (10) i - if (i /= maxcat) then - print *,'NTILES BCs/LDASsa mismatch:', i,maxcat - stop - endif - - allocate (tile_id_vec (1: maxcat)) - allocate (bcs2ldas (1: maxcat)) - read (10) tile_id_vec - close (10, status = 'keep') - - ! indexing to the LDASsa order - - do i = 1, maxcat - BCS2LDAS(tile_id_vec(i)) = i - end do - - ! read tile_id raster and index according to the order of LDASsa - - open (10,file=trim(BCSDIR)//'rst/'//trim(GFILE)//'.rst',status='old',action='read', & - form='unformatted',convert='little_endian') - - ALLOCATE (tile_id (1:nc_rst,1:nr_rst)) - - DO j = 1, nr_rst - read (10) tile_id (:, j) - END DO - - close (10, status ='keep') - - deallocate (tile_id_vec) - - ! find catchment-tiles that contribute to each 10x10 tile - - tile_count = 0 - Nb_max = 0 - - DO J = FLOOR (limits(1) + TILESIZE/2), CEILING (limits (3) - TILESIZE/2), NINT(TILEINT) - DO I = FLOOR (limits(2) + TILESIZE/2), CEILING (limits (4) - TILESIZE/2), NINT(TILEINT) - if (associated (tile_id_box)) NULLIFY (tile_id_box) - jj = (j + 90)*nc_rst/360 - dij/2 - ii = (i + 180)*nc_rst/360 - dij/2 - tile_id_box => tile_id (ii + 1 : ii + dij, jj +1 : jj + dij) - - NPLUS = count(tile_id_box >= 1 .and. tile_id_box <= maxcat) - - if(NPLUS > 0) then - - allocate (loc_int (1:NPLUS)) - allocate (unq_mask(1:NPLUS)) - loc_int = pack(tile_id_box,mask = (tile_id_box >= 1 .and. tile_id_box <= maxcat)) - call MAPL_Sort (loc_int) - unq_mask = .true. - do n = 2,NPLUS - unq_mask(n) = .not.(loc_int(n) == loc_int(n-1)) - end do - NB = count(unq_mask) - tile_count = tile_count + 1 - allocate(loc_val (1:NB)) - loc_val = 1.*pack(loc_int,mask =unq_mask) - - IF(NB_MAX < NB) NB_MAX = NB - - if(NB > MXT) then - print *, 'NB EXCEEDED MXT', NB, tile_count - stop - endif - - STATUS = NF_PUT_VARA_INT (NCFOutID2,VARID(NCFOutID2,'nSMAP' ),(/tile_count/), (/1/), NB);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCFOutID2,VARID(NCFOutID2,'lat' ),(/tile_count/), (/1/), real(j));VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCFOutID2,VARID(NCFOutID2,'lon' ),(/tile_count/), (/1/), real(i));VERIFY_(STATUS) - STATUS = NF_PUT_VARA_INT (NCFOutID2,VARID(NCFOutID2,'SMAPID'),(/1,tile_count/),(/NB, 1/), loc_val);VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCFOutID ,VARID(NCFOutID, 'lat' ),(/tile_count/), (/1/), real(j));VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCFOutID ,VARID(NCFOutID, 'lon' ),(/tile_count/), (/1/), real(i));VERIFY_(STATUS) - -! do k = 1,NBINS -! print *,k, loc_val(k), BCS2LDAS(loc_val(k)) -! STATUS = NF_PUT_VARA_INT (NCFOutID,VARID(NCFOutID,'SMAPID' ),(/k,tile_count/), (/1, 1/), BCS2LDAS(loc_val(k)));VERIFY_(STATUS) -! print *, k,nbins,BCS2LDAS(loc_val(k)) -! end do - deallocate (loc_val,loc_int,unq_mask) - - endif - END DO - END DO - - PRINT *, 'NB_MAX :', NB_MAX - - status = NF_CLOSE (NCFOutID ) - status = NF_CLOSE (NCFOutID2) - - deallocate (abins) - -END SUBROUTINE create_CDF_ParamFile - -! ---------------------------------------------------------------------- - -integer function VarID (NCFID, VNAME) - - integer, intent (in) :: NCFID - character(*), intent (in) :: VNAME - integer :: status - - STATUS = NF_INQ_VARID (NCFID, trim(VNAME) ,VarID) - IF (STATUS .NE. NF_NOERR) & - CALL HANDLE_ERR(STATUS, trim(VNAME)) - -end function VarID - -! ----------------------------------------------------------------------- - -SUBROUTINE HANDLE_ERR(STATUS, Line) - - INTEGER, INTENT (IN) :: STATUS - CHARACTER(*), INTENT (IN) :: Line - - IF (STATUS .NE. NF_NOERR) THEN - PRINT *, trim(Line),': ',STATUS, NF_STRERROR(STATUS) - STOP 'Stopped' - ENDIF - -END SUBROUTINE HANDLE_ERR - -! ---------------------------------------------------------------------- - -SUBROUTINE HISTOGRAM (NLENS, NBINS, density, loc_val, x, BIN) - - implicit none - - integer, intent (in) :: NBINS, NLENS - real, intent (in) :: x (NLENS) - integer, intent (out):: density (NBINS) - real, intent (inout) :: loc_val (NBINS) - real, intent (in), optional :: bin - real :: xdum(NLENS), xl, xu, min_value - integer :: n - - if(present (bin)) min_value = real(floor(minval(x))) - - DO N = 1, NBINS - if(present (bin)) then - xl = (N - 1)*BIN + min_value - loc_val (n) = xl - xu = xl + bin - XDUM = 0. - where((x >= xl).and.(x < xu))XDUM = 1 - else - XDUM = 0. - where(x == loc_val (n)) XDUM = 1 - endif - density(n) = int(sum(XDUM)) - END DO - -END SUBROUTINE HISTOGRAM - -!---------------------------------------------------------------------- - -END PROGRAM comp_FPAR_CDF diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/utils/math_routines.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/utils/math_routines.F90 deleted file mode 100755 index 88a6ee604..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/GEOSland_GridComp/GEOScatchCN_GridComp/utils/math_routines.F90 +++ /dev/null @@ -1,1186 +0,0 @@ -MODULE math_routines - -implicit none -private -public :: gammln, betai, betacf, rsq, prob_den_func, nbins, & - N_RANDOM_YEARS, shuffle, OPTIMIZ, N_PARAMS !,init_MPI -real, save :: U(97), CS, CD, CM -integer, save :: I97, J97 -integer, parameter :: N_PARAMS=3 -integer, parameter :: nbins = 40, N_RANDOM_YEARS = 41 - -! initialize to non-MPI values - include 'mpif.h' -!integer,public :: myid=0, numprocs=1, mpierr, mpistatus(MPI_STATUS_SIZE) -!logical, public :: root_proc=.true. - -contains - -! -! ---------------------------------------------------------------- -! - -subroutine Shuffle(a) - integer, intent(inout) :: a(:) - integer :: i, randpos, temp - real :: r - - call init_random_seed () - - do i = size(a), 2, -1 - call random_number(r) - randpos = int(r * i) + 1 - temp = a(randpos) - a(randpos) = a(i) - a(i) = temp - end do - end subroutine Shuffle - -! -! ------------------------------------------- -! - - SUBROUTINE init_random_seed() - INTEGER :: i, n, clock - INTEGER, DIMENSION(:), ALLOCATABLE :: seed - - CALL RANDOM_SEED(size = n) - ALLOCATE(seed(n)) - - CALL SYSTEM_CLOCK(COUNT=clock) - - seed = clock + 37 * (/ (i - 1, i = 1, n) /) - CALL RANDOM_SEED(PUT = seed) - - DEALLOCATE(seed) - -END SUBROUTINE init_random_seed - -! -! --------------------------------------------------------------------- -! - - SUBROUTINE random_rain(eta) - -!====================================================================== -! NPAC $Id: math_routines.F90,v 1.1 2018/05/10 15:05:29 smahanam Exp $ -!====================================================================== -! _____________________________________________________________________ -! | -! Test program to generate a vector of Gaussian deviates using | -! the Box-Muller method and the system-supplied uniform RNG. | -! K.A.Hawick, 15 July 1994. | -! | -! H W Yau. 23rd of August, 1996. | -! Northeast Parallel Architectures Center. | -! Syracuse University. | -!_____________________________________________________________________| -! | -! Edit record: | -! --/Jul/1996: M.McMahon -- ALIGN statement and EXTRINSIC statements | -! PROCESSORS statement as well. | -! 23/Aug/1996: HWY. Cleaning up. | -! Removed superfluous definition of module. | -! 26/Aug/1996: Added `IMPLICIT NONE'. | -! Fixed HPF directives. Cleaned up interface block to | -! gaussian_vector(). | -! F90 version. | -!_____________________________________________________________________| -! - - IMPLICIT NONE - - INTEGER, PARAMETER :: n = N_RANDOM_YEARS - REAL ti, tf,tcalc,tcom,tend,start - REAL, DIMENSION(1:n) :: a - REAL, DIMENSION(1:N)::ETA -! -! The vector of Guassian deviates - INTEGER, DIMENSION(1:n) :: j -! -! Histograms to check the distribution - INTEGER, DIMENSION(1:nbins) :: bins - INTEGER i,less, more,nb - INTEGER :: nover = 0 -! -! INTERFACE -! SUBROUTINE timer(return_time, initial_time) -! REAL, INTENT(IN) :: initial_time -! REAL, INTENT(OUT) :: return_time -! END SUBROUTINE -! END INTERFACE -! -! INTERFACE -! SUBROUTINE gaussian_vector(x,n,tcalc2,tcomm2) -! REAL, DIMENSION(:) :: x -! REAL,INTENT(INOUT) :: tcalc2,tcomm2 -! INTEGER,INTENT(IN) :: n -! REAL ETA -! END SUBROUTINE gaussian_vector -! END INTERFACE -!______________________________________________________________________ -! -! Executable Code. -!______________________________________________________________________ -! -! start_timer - CALL timer(start,0.0) - tcom = 0.0 - tcalc = 0.0 - CALL timer(ti,0.0) - bins(1:nbins) = 0 - CALL timer(tend,ti) - tcalc = tcalc + tend -! -! Ask for n deviate - CALL gaussian_vector( a, n,tend,tcom ) - tcalc = tcalc + tend - CALL timer(ti,0.0) -!F95 FORALL(i=1:n)j(i) = a(i) * real(nbins/7) + nbins/2 + 1 - DO i=1,n - j(i) = a(i) * real(nbins/7) + nbins/2 + 1 - END DO - CALL timer(tend,ti) - tcalc = tcalc + tend -! - CALL timer(ti,0.0) - DO i=1,n - ETA(I)=(REAL(J(I))-10.-1.)/5. - ENDDO -! write(*,*)eta -! write(*,*)sum(j)/300. - DO nb = 1,nbins - bins(nb) = count(j.EQ.nb,1) -!CCCCCCCC write(6,*) nb, bins(nb) - ENDDO - CALL timer(tend,ti) - tcom = tcom + tend -! -! Stop timer. - CALL timer(tf,start) -! -! Correction for serial execution. - tcalc = tcalc + tcom - tcom = 0.0 -! WRITE(6,6001) 0,n, -! 1 tcom,tcalc,(tf-tcom-tcalc),tf -! -! Write histogram - DO nb=1,nbins -! WRITE(6,*) nb,(REAL(NB)-10.-1.)/5., bins(nb) - ENDDO -! -! STOP - RETURN -! - 6001 FORMAT('Number of Processors = ',I4/ & - 'Problem size = ',I6/ & - 'Communications = ',F9.3/ & - 'Compute = ',F9.3/ & - 'Others = ',F9.3/ & - 'Total time = ',F9.3) -! - END SUBROUTINE random_rain -! -! ------------------------------------------------------ -! - - SUBROUTINE gaussian_vector( x, n, tcalc2,tcomm2 ) - - IMPLICIT NONE - ! - ! Box-Muller Method - ! See Knuth, Vol 2, 2nd Edn, PP 117 - INTEGER,INTENT(IN) :: n - REAL, DIMENSION(:) :: x - REAL,INTENT(INOUT) :: tcalc2,tcomm2 - REAL, DIMENSION(:), ALLOCATABLE :: v1, v2, r, f - REAL ti,cal1,cal2,com1,com2 - LOGICAL, DIMENSION(:), ALLOCATABLE :: mask - ! - ! Accept/reject efficiency is about 1.27 - integer m, i, np - !______________________________________________________________________ - ! - ! Executable code. - !______________________________________________________________________ - ! - CALL timer(ti,0.0) - np = n - ! - ! avoid fluctuation problems - IF( np .lt. 10 ) np = 10 - ALLOCATE( v1(np), v2(np), r(np), f(np), mask(np) ) - ! - ! Generate two deviates: - CALL random_number( v1 ) - CALL random_number( v2 ) - v1 = 2.0 * v1 -1.0 - v2 = 2.0 * v2 -1.0 - r = v1**2 + v2**2 - ! - ! are they in the unit circle? - mask = r < 1.0 - r = merge( r, 0.5, mask ) - f = sqrt( -2.0 * log(r) / r) - v1 = v1 * f - v2 = v2 * f - ! - ! since pack is a new intrinsic, here is serial code to show - ! what is being done. - ! m = 0 - ! do i=1,np - ! if( mask(i) )then - ! m = m + 1 - ! r(m) = v1(i) - ! f(m) = v2(i) - ! endif - ! enddo - CALL timer(cal1,ti) - - CALL timer(ti,0.0) - ! - ! we now have 2 * m deviates - m = count( mask ) - CALL timer(com1,ti) - - CALL timer(ti,0.0) - ! - ! and to save space, we'll store them in r and f - ! r = pack( v1, mask ) - ! f = pack( v2, mask ) - ! - ! We now get performance at the expense of memory. - r = v1 - f = v2 - ! - ! Statistically, this should not happen for large n - IF( 2*m .lt. n )THEN - WRITE(6,*) 'Not enough deviates! Got: ', 2*m, ', Needed: ', n - ! WRITE(6,*) 'Increase accept reject allowance in', - ! ' xgaussian_vector' - STOP - ENDIF - ! - ! use the two result vectors to patch up enough as - IF( m .LT. n )THEN - x(1:m)=r(1:m) - CALL timer(cal2,ti) - CALL timer(ti,0.0) - x(m+1:n) = f(1:n-m) - CALL timer(com2,ti) - ELSE - x(1:n) = r(1:n) - ENDIF - - tcalc2 = cal1 + cal2 - tcomm2 = com1 + com2 - DEALLOCATE( v1, v2, r, f, mask ) - - RETURN - END SUBROUTINE gaussian_vector -! -! ------------------------------------------------------------------- -! - SUBROUTINE timer(return_time, initial_time) - implicit none - REAL, INTENT(IN) :: initial_time - REAL, INTENT(OUT) :: return_time - INTEGER finish,rate - CALL system_clock( COUNT=finish,COUNT_RATE=rate) - return_time = FLOAT(finish) / FLOAT(rate) - initial_time - RETURN - END SUBROUTINE timer -! -! ------------------------------------------------------------------- -! - -subroutine prob_den_func (ndata, datain, cdf, bins, & - mean, std, skew, pdf, lwval, upval) - -implicit none -integer, intent (in) :: ndata -real, dimension(ndata), intent (in) :: datain -real, intent(out), dimension(nbins) :: cdf, bins -real, optional, intent(out), dimension(nbins) :: pdf -real, optional, intent (out) :: mean, std, skew -real, optional :: lwval, upval -real :: lw,up,db,var1,var2, variance -integer :: i,n - -lw = minval (datain) -up = maxval (datain) - -if(present(upval)) then - up = upval - lw = lwval -endif - -db = (up - lw)/real(nbins) -cdf = 0. -if(present(pdf)) pdf =0. - -do n = 1, nbins - bins (n) = lw + real(n)*db - db/2. - do i = 1,ndata - if(datain(i) <= bins (n) + db/2.) cdf(n) = cdf(n) + 1. - if(present(pdf)) then - if(n==1) then - if((datain(i) >= bins (n) - db/2.).and. & - (datain(i) <= bins (n) + db/2.)) & - pdf(n) = pdf(n) + 1. - else - if((datain(i) > bins (n) - db/2.).and. & - (datain(i) <= bins (n) + db/2.)) & - pdf(n) = pdf(n) + 1. - endif - endif - end do -end do - -cdf = cdf/real(ndata) - -if(present(pdf)) pdf = pdf/real(ndata) - -if(present (mean)) mean = sum(datain)/real(ndata) - -if(present (std)) then - - var1 = 0. - mean = sum(datain)/real(ndata) - - do i = 1,ndata - var1 = var1 + (datain(i) - mean)*(datain(i) - mean) - end do - - std = sqrt (var1/real(ndata - 1)) - - if(present (skew)) then - var2 = 0. - do i = 1,ndata - var2 = var2 + ((datain(i) - mean)/std)* & - ((datain(i) - mean)/std)* & - ((datain(i) - mean)/std) - end do - - skew = var2/real(ndata - 1) - - endif - -endif - -END subroutine prob_den_func - -! -! ------------------------------------------------------- -! - -FUNCTION betai(a,b,x) -REAL betai,a,b,x -REAL bt -!external gammln - -if(x.lt.0..or.x.gt.1.)print *, 'bad argument x in betai',x -if(x.lt.0..or.x.gt.1.)stop -if(x.eq.0..or.x.eq.1.)then - bt=0. -else - bt=exp(gammln(a+b)-gammln(a)-gammln(b) & - +a*log(x)+b*log(1.-x)) -endif - -if(x.lt.(a+1.)/(a+b+2.))then - betai=bt*betacf(a,b,x)/a - return -else - betai=1.-bt*betacf(b,a,1.-x)/b - return -endif -END FUNCTION betai -! -! ------------------------------------------------------- -! -FUNCTION betacf(a,b,x) -INTEGER MAXIT -REAL betacf,a,b,x,EPS,FPMIN -PARAMETER (MAXIT=100,EPS=3.e-7,FPMIN=1.e-30) -INTEGER m,m2 -REAL aa,c,d,del,h,qab,qam,qap - -qab=a+b -qap=a+1. -qam=a-1. -c=1. -d=1.-qab*x/qap - -if(abs(d).lt.FPMIN)d=FPMIN -d=1./d -h=d -do m=1,MAXIT - m2=2*m - aa=m*(b-m)*x/((qam+m2)*(a+m2)) - d=1.+aa*d - if(abs(d).lt.FPMIN)d=FPMIN - c=1.+aa/c - if(abs(c).lt.FPMIN)c=FPMIN - d=1./d - h=h*d*c - aa=-(a+m)*(qab+m)*x/((a+m2)*(qap+m2)) - d=1.+aa*d - if(abs(d).lt.FPMIN)d=FPMIN - c=1.+aa/c - if(abs(c).lt.FPMIN)c=FPMIN - d=1./d - del=d*c - h=h*del - if(abs(del-1.).lt.EPS)exit -enddo - betacf=h -return -END FUNCTION betacf -! -! -------------------------------------------------------------- -! -FUNCTION gammln(xx) -REAL gammln,xx -INTEGER j -DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) - -SAVE cof,stp -DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, & - 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, & - -.5395239384953d-5,2.5066282746310005d0/ -x=xx -y=x -tmp=x+5.5d0 -tmp=(x+0.5d0)*log(tmp)-tmp -ser=1.000000000190015d0 -do j=1,6 - y=y+1.d0 - ser=ser+cof(j)/y -enddo -gammln=tmp+log(stp*ser/x) -return -END FUNCTION gammln - -! -! ------------------------------------------------ -! - -SUBROUTINE RSQ (NDATA, X, Y, R2, RMSE, limits, slope, intercept) - -implicit none -integer, intent (in) :: NDATA -real, dimension (ndata), intent(in) :: x,y -real, optional, dimension(4), intent (in) :: limits -real, optional, intent(out) :: slope, intercept, RMSE -real, intent(out) :: r2 - -integer ::n -real :: ic -real :: sumx, sumy,sumxy, sumx2,sumy2,sig2x,sig2y,sig2xy,error - -ic =0. -sumx=0. -sumy=0. -sumxy=0. -sumy2=0. -sumx2=0. -error=0. -do n = 1,ndata - - if(present (limits)) then - if((x(n) > limits(1)).and.(x(n) < limits(2)).and. & - (y(n) > limits(3)).and.(y(n) < limits(4))) then - ic = ic + 1. - sumx = sumx + x(n) - sumy = sumy + y(n) - sumxy= sumxy + x(n)*y(n) - sumy2= sumy2 + y(n)*y(n) - sumx2= sumx2 + x(n)*x(n) - error= error + (y(n) - x(n))*(y(n) - x(n)) - endif - else - ic = ic + 1. - sumx = sumx + x(n) - sumy = sumy + y(n) - sumxy= sumxy + x(n)*y(n) - sumy2= sumy2 + y(n)*y(n) - sumx2= sumx2 + x(n)*x(n) - error= error + (y(n) - x(n))*(y(n) - x(n)) - endif - -end do - -if (present(intercept)) intercept = -9999. -if (present(slope)) slope = -9999. -if (present(rmse)) rmse = -9999. -r2 = -9999. - -if(ic /= 0) then - if(ic*sumx2 /= sumx*sumx) then - - if (present(intercept)) intercept = & - (sumy*sumx2 - sumx*sumxy) / (ic*sumx2 - sumx*sumx) - if (present(slope)) slope = & - (ic*sumxy - sumx*sumy) / (ic*sumx2 - sumx*sumx) - if (present(rmse)) rmse = sqrt(error/real(ic)) - endif - - sumx =sumx/ic - sumy =sumy/ic - sumxy=sumxy/ic - sumy2=sumy2/ic - sumx2=sumx2/ic - sig2x=sumx2-sumx*sumx - sig2y=sumy2-sumy*sumy - sig2xy=(sumxy-sumx*sumy)*(sumxy-sumx*sumy) - - r2 = sig2xy/(sig2x*sig2y + 1.e-20) -endif - -END SUBROUTINE RSQ - -! -! ------------------------------------------ -! - SUBROUTINE OPTIMIZ(NDATA,wet,eff,X) - - implicit none - - integer, PARAMETER :: N = N_PARAMS, NEPS = 4 - integer, intent (in) :: ndata - REAL (kind = 8) :: LB(N), UB(N), X(N), XOPT(N), CON(N), VM(N), & - FSTAR(NEPS), XP(N), T, EPS, RT, FOPT, & - EFF(NDATA),WET(NDATA) - - INTEGER NACP(N), NS, NT, NFCNEV, IER, ISEED1, ISEED2, & - MAXEVL, IPRINT, NACC, NOBDS, I - - LOGICAL MAX - -! Set underflows to zero on IBM mainframes. -! CALL XUFLOW(0) - - MAX = .false. - EPS = 1.0D-6 - RT = .5 - ISEED1 = 1 - ISEED2 = 2 - NS = 20 - NT = 5 - MAXEVL = 100000 - IPRINT = 0 - - DO I=1,N - LB(I)=0.0001 - UB(I)=100. - CON(I) = 2.0 - END DO - UB(3)= 1. - LB(3)= 0.1 -! -! Set input values of the input/output parameters. -! - T = 5.0 - DO I = 1, N - VM(I) = 1.0 - END DO - -! WRITE(*,1000) N, MAX, T, RT, EPS, NS, NT, NEPS, MAXEVL, IPRINT, & -! ISEED1, ISEED2 -! -! CALL PRTVEC(X,N,'STARTING VALUES') -! CALL PRTVEC(VM,N,'INITIAL STEP LENGTH') -! CALL PRTVEC(LB,N,'LOWER BOUND') -! CALL PRTVEC(UB,N,'UPPER BOUND') -! CALL PRTVEC(C,N,'C VECTOR') -! WRITE(*,'(/,'' **** END OF DRIVER ROUTINE OUTPUT ****'' & -! /,'' **** BEFORE CALL TO SA. ****'')') - - CALL SA(X,MAX,RT,EPS,NS,NT,NEPS,MAXEVL,LB,UB,CON,IPRINT,ISEED1, & - ISEED2,T,VM,XOPT,FOPT,NACC,NFCNEV,NOBDS,IER, & - FSTAR,XP,NACP,NDATA,EFF,WET) - -! WRITE(*,'(/,'' **** RESULTS AFTER SA **** '')') -! CALL PRTVEC(XOPT,N,'SOLUTION') -! CALL PRTVEC(VM,N,'FINAL STEP LENGTH') -! WRITE(*,1001) FOPT, NFCNEV, NACC, NOBDS, T, IER - -1000 FORMAT(/,' SIMULATED ANNEALING EXAMPLE',/, & - /,' NUMBER OF PARAMETERS: ',I3,' MAXIMAZATION: ',L5, & - /,' INITIAL TEMP: ', G8.2, ' RT: ',G8.2, ' EPS: ',G8.2, & - /,' NS: ',I3, ' NT: ',I2, ' NEPS: ',I2, & - /,' MAXEVL: ',I10, ' IPRINT: ',I1, ' ISEED1: ',I4, & - ' ISEED2: ',I4) -1001 FORMAT(/,' OPTIMAL FUNCTION VALUE: ',G20.13 & - /,' NUMBER OF FUNCTION EVALUATIONS: ',I10, & - /,' NUMBER OF ACCEPTED EVALUATIONS: ',I10, & - /,' NUMBER OF OUT OF BOUND EVALUATIONS: ',I10, & - /,' FINAL TEMP: ', G20.13,' IER: ', I3) - - RETURN - END SUBROUTINE OPTIMIZ - -! -! --------------------------------------------------------------- -! - - SUBROUTINE SA(X,MAX,RT,EPS,NS,NT,NEPS,MAXEVL,LB,UB,CON,IPRINT, & - ISEED1,ISEED2,T,VM,XOPT,FOPT,NACC,NFCNEV,NOBDS,IER, & - FSTAR,XP,NACP,NDATA,EFF,WET) - -! Type all external variables. - - INTEGER,PARAMETER :: N=N_PARAMS - integer :: ndata - REAL (KIND = 8) X(N), LB(N), UB(N), CON(N), VM(N), FSTAR(N), & - XOPT(N), XP(N), T, EPS, RT, FOPT - REAL (KIND = 8) EFF(NDATA),WET(NDATA) - INTEGER NACP(N), NS, NT, NEPS, NACC, MAXEVL, IPRINT, & - NOBDS, IER, NFCNEV, ISEED1, ISEED2 - LOGICAL MAX - -! Type all internal variables. - REAL (KIND = 8) F, FP, P, PP, RATIO - INTEGER NUP, NDOWN, NREJ, NNEW, LNOBDS, H, I, J, M - LOGICAL QUIT - -! Type all functions. -! REAL (KIND = 8) EXPREP - -! Initialize the random number generator RANMAR. - CALL RMARIN(ISEED1,ISEED2) - -! Set initial values. - NACC = 0 - NOBDS = 0 - NFCNEV = 0 - IER = 99 - - DO I = 1, N - XOPT(I) = X(I) - NACP(I) = 0 - END DO - - DO I = 1, NEPS - FSTAR(I) = 1.0D+20 - END DO - -! If the initial temperature is not positive, notify the user and -! return to the calling routine. - IF (T .LE. 0.0) THEN - WRITE(*,'(/,'' THE INITIAL TEMPERATURE IS NOT POSITIVE. '' & - /,'' RESET THE VARIABLE T. ''/)') - IER = 3 - RETURN - END IF - -! If the initial value is out of bounds, notify the user and return -! to the calling routine. - DO I = 1, N - IF ((X(I) .GT. UB(I)) .OR. (X(I) .LT. LB(I))) THEN - CALL PRT1 - IER = 2 - RETURN - END IF - END DO - -! Evaluate the function with input X and return value as F. - CALL FCN(X,F,NDATA,EFF,WET) - -! If the function is to be minimized, switch the sign of the function. -! Note that all intermediate and final output switches the sign back -! to eliminate any possible confusion for the user. - IF(.NOT. MAX) F = -F - NFCNEV = NFCNEV + 1 - FOPT = F - FSTAR(1) = F - IF(IPRINT .GE. 1) CALL PRT2(MAX,N,X,F) - -! Start the main loop. Note that it terminates if (i) the algorithm -! succesfully optimizes the function or (ii) there are too many -! function evaluations (more than MAXEVL). -100 NUP = 0 - NREJ = 0 - NNEW = 0 - NDOWN = 0 - LNOBDS = 0 - - DO M = 1, NT - DO J = 1, NS - DO H = 1, N - -! Generate XP, the trial value of X. Note use of VM to choose XP. - DO I = 1, N - IF (I .EQ. H) THEN - XP(I) = X(I) + (RANMAR()*2.- 1.) * VM(I) - ELSE - XP(I) = X(I) - END IF - -! If XP is out of bounds, select a point in bounds for the trial. - IF((XP(I) .LT. LB(I)) .OR. (XP(I) .GT. UB(I))) THEN - XP(I) = LB(I) + (UB(I) - LB(I))*RANMAR() - LNOBDS = LNOBDS + 1 - NOBDS = NOBDS + 1 - IF(IPRINT .GE. 3) CALL PRT3(MAX,N,XP,X,FP,F) - END IF - END DO -! Evaluate the function with the trial point XP and return as FP. - CALL FCN(XP,FP,NDATA,EFF,WET) - IF(.NOT. MAX) FP = -FP - NFCNEV = NFCNEV + 1 - IF(IPRINT .GE. 3) CALL PRT4(MAX,N,XP,X,FP,F) - -! If too many function evaluations occur, terminate the algorithm. - IF(NFCNEV .GE. MAXEVL) THEN - CALL PRT5 - IF (.NOT. MAX) FOPT = -FOPT - IER = 1 - RETURN - END IF - -! Accept the new point if the function value increases. - IF(FP .GE. F) THEN - IF(IPRINT .GE. 3) THEN - WRITE(*,'('' POINT ACCEPTED'')') - END IF - DO I = 1, N - X(I) = XP(I) - END DO - F = FP - NACC = NACC + 1 - NACP(H) = NACP(H) + 1 - NUP = NUP + 1 - -! If greater than any other point, record as new optimum. - IF (FP .GT. FOPT) THEN - IF(IPRINT .GE. 3) THEN - WRITE(*,'('' NEW OPTIMUM'')') - END IF - DO I = 1, N - XOPT(I) = XP(I) - END DO - FOPT = FP - NNEW = NNEW + 1 - END IF - - ! If the point is lower, use the Metropolis criteria to decide on - ! acceptance or rejection. - ELSE - P = EXPREP((FP - F)/T) - PP = RANMAR() - IF (PP .LT. P) THEN - IF(IPRINT .GE. 3) CALL PRT6(MAX) - DO I = 1, N - X(I) = XP(I) - END DO - F = FP - NACC = NACC + 1 - NACP(H) = NACP(H) + 1 - NDOWN = NDOWN + 1 - ELSE - NREJ = NREJ + 1 - IF(IPRINT .GE. 3) CALL PRT7(MAX) - END IF - END IF - - END DO - END DO - -! Adjust VM so that approximately half of all evaluations are accepted. - DO I = 1, N - RATIO = DFLOAT(NACP(I)) /DFLOAT(NS) - IF (RATIO .GT. .6) THEN - VM(I) = VM(I)*(1. + CON(I)*(RATIO - .6)/.4) - ELSE IF (RATIO .LT. .4) THEN - VM(I) = VM(I)/(1. + CON(I)*((.4 - RATIO)/.4)) - END IF - IF (VM(I) .GT. (UB(I)-LB(I))) THEN - VM(I) = UB(I) - LB(I) - END IF - END DO - - IF(IPRINT .GE. 2) THEN - CALL PRT8(N,VM,XOPT,X) - END IF - - DO I = 1, N - NACP(I) = 0 - END DO - - END DO - - IF(IPRINT .GE. 1) THEN - CALL PRT9(MAX,N,T,XOPT,VM,FOPT,NUP,NDOWN,NREJ,LNOBDS,NNEW) - END IF - -! Check termination criteria. - QUIT = .FALSE. - FSTAR(1) = F - IF ((FOPT - FSTAR(1)) .LE. EPS) QUIT = .TRUE. - DO I = 1, NEPS - IF (ABS(F - FSTAR(I)) .GT. EPS) QUIT = .FALSE. - END DO - -! Terminate SA if appropriate. - IF (QUIT) THEN - DO I = 1, N - X(I) = XOPT(I) - END DO - IER = 0 - IF (.NOT. MAX) FOPT = -FOPT - IF(IPRINT .GE. 1) CALL PRT10 - RETURN - END IF - -! If termination criteria is not met, prepare for another loop. - T = RT*T - DO I = NEPS, 2, -1 - FSTAR(I) = FSTAR(I-1) - END DO - F = FOPT - DO I = 1, N - X(I) = XOPT(I) - END DO - -! Loop again. - GO TO 100 - - END SUBROUTINE SA -! -! --------------------------------------------------------- -! - - SUBROUTINE FCN(X,F,NDATA,EFF,WET) - implicit none - integer, parameter :: N = N_PARAMS - REAL (KIND = 8) X(n),F - INTEGER, intent (in) :: ndata - INTEGER Nd - REAL (KIND = 8) EFF(NDATA),WET(NDATA) - real :: a,b,xv, r2, rmse - real , dimension(ndata) :: yval, xval - - a = X(1) - b = X(2) - - do nd = 1, ndata - xv = wet(nd) - xval(nd) = eff(nd) - yval(nd) = X(3)*betai(a,b,xv) - end do - -! if(maxval(yval) > 1.) print *,maxval(yval), minval(yval) - - call rsq (ndata, xval, yval, r2, rmse) - F = rmse -! if(maxval(yval) > 1.) F = 1.d10 -! F = r2 -! if(F > 1.) F = -1.d-10 - - RETURN - END SUBROUTINE FCN - -! -! --------------------------------------------------------------- -! - FUNCTION EXPREP(RDUM) - implicit none - REAL (KIND = 8) RDUM, EXPREP - - IF (RDUM .GT. 174.) THEN - EXPREP = 3.69D+75 - ELSE IF (RDUM .LT. -180.) THEN - EXPREP = 0.0 - ELSE - EXPREP = EXP(RDUM) - END IF - - RETURN - END FUNCTION EXPREP -! -! --------------------------------------------- -! - subroutine RMARIN(IJ,KL) - - implicit none - integer, intent (in) :: ij,kl - integer :: i,j,k,l,m, ii,jj - real :: t,s - -! real U(97), C, CD, CM -! integer I97, J97 -! common /raset1/ U, C, CD, CM, I97, J97 - - if( IJ .lt. 0 .or. IJ .gt. 31328 .or. & - KL .lt. 0 .or. KL .gt. 30081 ) then - print '(A)', ' The first random number seed must have a value between 0 and 31328' - print '(A)',' The second seed must have a value between 0 and 30081' - stop - endif - - i = mod(IJ/177, 177) + 2 - j = mod(IJ , 177) + 2 - k = mod(KL/169, 178) + 1 - l = mod(KL, 169) - - do ii = 1, 97 - s = 0.0 - t = 0.5 - do jj = 1, 24 - m = mod(mod(i*j, 179)*k, 179) - i = j - j = k - k = m - l = mod(53*l+1, 169) - if (mod(l*m, 64) .ge. 32) then - s = s + t - endif - t = 0.5 * t - end do - U(ii) = s - end do - - - CS = 362436.0 / 16777216.0 - CD = 7654321.0 / 16777216.0 - CM = 16777213.0 /16777216.0 - I97 = 97 - J97 = 33 - return - end subroutine RMARIN - -! -! --------------------------------------- -! - - real function ranmar() - real :: uni -! real U(97), C, CD, CM -! integer I97, J97 -! common /raset1/ U, C, CD, CM, I97, J97 - uni = U(I97) - U(J97) - if( uni .lt. 0.0 ) uni = uni + 1.0 - U(I97) = uni - I97 = I97 - 1 - if(I97 .eq. 0) I97 = 97 - J97 = J97 - 1 - if(J97 .eq. 0) J97 = 97 - CS = CS - CD - if( CS .lt. 0.0 ) CS = CS + CM - uni = uni - CS - if( uni .lt. 0.0 ) uni = uni + 1.0 - RANMAR = uni - return - END function ranmar - - SUBROUTINE PRT1 - implicit none - WRITE(*,'(/,'' THE STARTING VALUE (X) IS OUTSIDE THE BOUNDS '' & - /,'' (LB AND UB). EXECUTION TERMINATED WITHOUT ANY'' & - /,'' OPTIMIZATION. RESPECIFY X, UB OR LB SO THAT '' & - /,'' LB(I) .LT. X(I) .LT. UB(I), I = 1, N. ''/)') - - RETURN - END SUBROUTINE PRT1 - - SUBROUTINE PRT2(MAX,N,X,F) - implicit none - REAL (KIND = 8) X(N), F - INTEGER N - LOGICAL MAX - - WRITE(*,'('' '')') - CALL PRTVEC(X,N,'INITIAL X') - IF (MAX) THEN - WRITE(*,'('' INITIAL F: '',/, G25.18)') F - ELSE - WRITE(*,'('' INITIAL F: '',/, G25.18)') -F - END IF - - RETURN - END SUBROUTINE PRT2 - - SUBROUTINE PRT3(MAX,N,XP,X,FP,F) - implicit none - REAL (KIND = 8) XP(N), X(N), FP, F - INTEGER N - LOGICAL MAX - - WRITE(*,'('' '')') - CALL PRTVEC(X,N,'CURRENT X') - IF (MAX) THEN - WRITE(*,'('' CURRENT F: '',G25.18)') F - ELSE - WRITE(*,'('' CURRENT F: '',G25.18)') -F - END IF - CALL PRTVEC(XP,N,'TRIAL X') - WRITE(*,'('' POINT REJECTED SINCE OUT OF BOUNDS'')') - - RETURN - END SUBROUTINE PRT3 - - SUBROUTINE PRT4(MAX,N,XP,X,FP,F) - implicit none - REAL (KIND = 8) XP(N), X(N), FP, F - INTEGER N - LOGICAL MAX - - WRITE(*,'('' '')') - CALL PRTVEC(X,N,'CURRENT X') - IF (MAX) THEN - WRITE(*,'('' CURRENT F: '',G25.18)') F - CALL PRTVEC(XP,N,'TRIAL X') - WRITE(*,'('' RESULTING F: '',G25.18)') FP - ELSE - WRITE(*,'('' CURRENT F: '',G25.18)') -F - CALL PRTVEC(XP,N,'TRIAL X') - WRITE(*,'('' RESULTING F: '',G25.18)') -FP - END IF - - RETURN - END SUBROUTINE PRT4 - - SUBROUTINE PRT5 - implicit none - WRITE(*,'(/,'' TOO MANY FUNCTION EVALUATIONS; CONSIDER '' & - /,'' INCREASING MAXEVL OR EPS, OR DECREASING '' & - /,'' NT OR RT. THESE RESULTS ARE LIKELY TO BE '' & - /,'' POOR.'',/)') - - RETURN - END SUBROUTINE PRT5 - - SUBROUTINE PRT6(MAX) - implicit none - LOGICAL MAX - - IF (MAX) THEN - WRITE(*,'('' THOUGH LOWER, POINT ACCEPTED'')') - ELSE - WRITE(*,'('' THOUGH HIGHER, POINT ACCEPTED'')') - END IF - - RETURN - END SUBROUTINE PRT6 - - SUBROUTINE PRT7(MAX) - implicit none - LOGICAL MAX - - IF (MAX) THEN - WRITE(*,'('' LOWER POINT REJECTED'')') - ELSE - WRITE(*,'('' HIGHER POINT REJECTED'')') - END IF - - RETURN - END SUBROUTINE PRT7 - - SUBROUTINE PRT8(N,VM,XOPT,X) - implicit none - REAL (KIND = 8) VM(N), XOPT(N), X(N) - INTEGER N - - WRITE(*,'(/,'' INTERMEDIATE RESULTS AFTER STEP LENGTH ADJUSTMENT'',/)') - CALL PRTVEC(VM,N,'NEW STEP LENGTH (VM)') - CALL PRTVEC(XOPT,N,'CURRENT OPTIMAL X') - CALL PRTVEC(X,N,'CURRENT X') - WRITE(*,'('' '')') - - RETURN - END SUBROUTINE PRT8 - - SUBROUTINE PRT9(MAX,N,T,XOPT,VM,FOPT,NUP,NDOWN,NREJ,LNOBDS,NNEW) - implicit none - REAL (KIND = 8) XOPT(N), VM(N), T, FOPT - INTEGER N, NUP, NDOWN, NREJ, LNOBDS, NNEW, TOTMOV - LOGICAL MAX - - TOTMOV = NUP + NDOWN + NREJ - - WRITE(*,'(/, '' INTERMEDIATE RESULTS BEFORE NEXT TEMPERATURE REDUCTION'',/)') - WRITE(*,'('' CURRENT TEMPERATURE: '',G12.5)') T - IF (MAX) THEN - WRITE(*,'('' MAX FUNCTION VALUE SO FAR: '',G25.18)') FOPT - WRITE(*,'('' TOTAL MOVES: '',I8)') TOTMOV - WRITE(*,'('' UPHILL: '',I8)') NUP - WRITE(*,'('' ACCEPTED DOWNHILL: '',I8)') NDOWN - WRITE(*,'('' REJECTED DOWNHILL: '',I8)') NREJ - WRITE(*,'('' OUT OF BOUNDS TRIALS: '',I8)') LNOBDS - WRITE(*,'('' NEW MAXIMA THIS TEMPERATURE:'',I8)') NNEW - ELSE - WRITE(*,'('' MIN FUNCTION VALUE SO FAR: '',G25.18)') -FOPT - WRITE(*,'('' TOTAL MOVES: '',I8)') TOTMOV - WRITE(*,'('' DOWNHILL: '',I8)') NUP - WRITE(*,'('' ACCEPTED UPHILL: '',I8)') NDOWN - WRITE(*,'('' REJECTED UPHILL: '',I8)') NREJ - WRITE(*,'('' TRIALS OUT OF BOUNDS: '',I8)') LNOBDS - WRITE(*,'('' NEW MINIMA THIS TEMPERATURE:'',I8)') NNEW - END IF - CALL PRTVEC(XOPT,N,'CURRENT OPTIMAL X') - CALL PRTVEC(VM,N,'STEP LENGTH (VM)') - WRITE(*,'('' '')') - - RETURN - END SUBROUTINE PRT9 - - SUBROUTINE PRT10 - implicit none - WRITE(*,'(/,'' SA ACHIEVED TERMINATION CRITERIA. IER = 0. '',/)') - - RETURN - END SUBROUTINE PRT10 - - SUBROUTINE PRTVEC(VECTOR,NCOLS,NAME) - implicit none - INTEGER NCOLS, LL,I,J, LINES - REAL (KIND = 8) VECTOR(NCOLS) - CHARACTER *(*) NAME - - WRITE(*,1001) NAME - - IF (NCOLS .GT. 10) THEN - LINES = INT(NCOLS/10.) - - DO I = 1, LINES - LL = 10*(I - 1) - WRITE(*,1000) (VECTOR(J),J = 1+LL, 10+LL) - END DO - - WRITE(*,1000) (VECTOR(J),J = 11+LL, NCOLS) - ELSE - WRITE(*,1000) (VECTOR(J),J = 1, NCOLS) - END IF - -1000 FORMAT( 10(G12.5,1X)) -1001 FORMAT(/,25X,A) - - RETURN - - END SUBROUTINE PRTVEC - - ! ***************************************************************************** - ! - ! subroutine init_MPI() - ! - ! ! initialize MPI - ! - ! call MPI_INIT(mpierr) - ! - ! call MPI_COMM_RANK( MPI_COMM_WORLD, myid, mpierr ) - ! call MPI_COMM_SIZE( MPI_COMM_WORLD, numprocs, mpierr ) - ! - ! if (myid .ne. 0) root_proc = .false. - ! -!! call init_MPI_types() - ! - ! write (*,*) "MPI process ", myid, " of ", numprocs, " is alive" - ! write (*,*) "MPI process ", myid, ": root_proc=", root_proc -! -! end subroutine init_MPI - - -END MODULE math_routines 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 42d3fe363..98900a487 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Shared/GEOS_SurfaceGridComp.rc @@ -234,14 +234,3 @@ # # GEOSagcm=>PRESCRIBE_DVG: 0 # GEOSldas=>PRESCRIBE_DVG: 0 - -# ---- Scale CATCHCN ALBEDO and FPAR -# -# 0 : NO scaling is performed (default) -# 1 : Scale albedo to match interannually varying MODIS NIRDF and VISDF anomaly -# 2 : Scale albedo to match CDFs of model fPAR to MODIS CDFs of fPAR -# 3 : Perform both 1 and 2 above -# -# GEOSagcm=>SCALE_ALBFPAR: 0 -# GEOSldas=>SCALE_ALBFPAR: 0 - diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt index 98956b0cd..21d1777ea 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt @@ -9,7 +9,6 @@ rasterize.F90 read_riveroutlet.F90 CubedSphere_GridMod.F90 rmTinyCatchParaMod.F90 -comp_CATCHCN_AlbScale_parameters.F90 zip.c util.c ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/comp_CATCHCN_AlbScale_parameters.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/comp_CATCHCN_AlbScale_parameters.F90 deleted file mode 100644 index 16bd9f0f2..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/comp_CATCHCN_AlbScale_parameters.F90 +++ /dev/null @@ -1,835 +0,0 @@ -#define VERIFY_(A) IF(A/=0)THEN;PRINT *,'ERROR AT LINE ', __LINE__;STOP;ENDIF -#define ASSERT_(A) if(.not.A)then;print *,'Error:',__FILE__,__LINE__;stop;endif - -MODULE comp_CATCHCN_AlbScale_parameters - - use date_time_util, ONLY: & - date_time_type, augment_date_time - - implicit none - INCLUDE 'netcdf.inc' - - private - - public :: albedo4catchcn - - character*400, PARAMETER :: & - InBCSDIR = '/discover/nobackup/smahanam/MERRA3/FPAR/SMAP_EASEv2_M09/', & - !EXPDIR = '/archive/u/smahanam/FPAR-ALB/e0004s_wet2/output/SMAP_EASEv2_M09_GLOBAL/', & - EXPDIR = '/discover/nobackup/borescan/BCS/add_to_l_land/e0004s_wet2/output/SMAP_EASEv2_M09_GLOBAL/', & - EXNAME = 'e0004s_wet2', & - InGFILE = 'SMAP_EASEv2_M09_3856x1624' - - ! character*400 :: GFILE = 'til/CF0180x6C_TM0720xTM0410-Pfafstetter.til' - real, parameter :: MAPL_PI = 3.14159265358979323846d0 - integer, parameter :: yearB = 2001, yearE = 2015, InNTILES = 1684725, NOCTAD = 46 - integer, parameter :: yearB1= 2002 - contains - - SUBROUTINE albedo4catchcn (gfile) - - implicit none - character (*), intent (in) :: gfile - integer :: NTILES - integer, dimension (:), allocatable :: id_loc - - call preprocess_m09 - open (10, file = 'clsm/catchment.def', form = 'formatted', status= 'old', action = 'read') - read (10, *) NTILES - close (10, status = 'keep') - - allocate (id_loc (1: NTILES)) - - call get_id_loc (NTILES, trim (GFILE)//'.til', id_loc) - call regrid_alb (NTILES, id_loc) - - end SUBROUTINE albedo4catchcn - - ! ------------------------------------------------------------------------------- - - SUBROUTINE get_id_loc (NT, gfile, id_loc) - - implicit none - - integer, intent (in) :: NT - integer, dimension (NT), intent (inout) :: id_loc - character(*), intent (in) :: gfile - integer :: n, i, nplus, t_count - real, dimension (:), allocatable :: lon, lat, m09_lon, m09_lat, tid_m09 - integer, allocatable, dimension (:) :: sub_tid - real , allocatable, dimension (:) :: sub_lon, sub_lat, rev_dist - real :: dw, dx, dy, min_lon, max_lon, min_lat, max_lat - logical :: tile_found - logical, allocatable, dimension(:) :: mask - integer, allocatable :: low_ind(:), upp_ind(:) - -! --------- VARIABLES FOR *OPENMP* PARALLEL ENVIRONMENT ------------ -! -! NOTE: "!$" is for conditional compilation -! -logical :: running_omp = .false. -! -!$ integer :: omp_get_thread_num, omp_get_num_threads -! -integer :: n_threads=1 -! -! -! ----------- OpenMP PARALLEL ENVIRONMENT ---------------------------- -! -! FIND OUT WHETHER -omp FLAG HAS BEEN SET DURING COMPILATION -! -!$ running_omp = .true. ! conditional compilation -! -! ECHO BASIC OMP VARIABLES -! -!$OMP PARALLEL DEFAULT(NONE) SHARED(running_omp,n_threads) -! -!$OMP SINGLE -! -!$ n_threads = omp_get_num_threads() -! -!$ write (*,*) 'running_omp = ', running_omp -!$ write (*,*) -!$ write (*,*) 'parallel OpenMP with ', n_threads, 'threads' -!$ write (*,*) -!$OMP ENDSINGLE -! -!$OMP CRITICAL -!$ write (*,*) 'thread ', omp_get_thread_num(), ' alive' -!$OMP ENDCRITICAL -! -!$OMP BARRIER -! -!$OMP ENDPARALLEL - - allocate (lon (1: NT)) - allocate (lat (1: NT)) - allocate (m09_lon (1: InNTILES)) - allocate (m09_lat (1: InNTILES)) - allocate (tid_m09 (1: InNTILES)) - - call ReadCNTilFile (trim(InBCSDIR)//trim(InGFILE)//'.til', InNTILES, m09_lon, m09_lat) - call ReadCNTilFile (trim(GFILE), NT, lon, lat) - - Id_loc = -9999 - do n = 1, InNTILES - tid_m09(n) = n - end do - - ! Domain decomposition - ! -------------------- - - allocate(low_ind(n_threads)) - allocate(upp_ind(n_threads)) - low_ind(1) = 1 - upp_ind(n_threads) = nt - - if (running_omp) then - do i=1,n_threads-1 - upp_ind(i) = low_ind(i) + (NT/n_threads) - 1 - low_ind(i+1) = upp_ind(i) + 1 - ! print *,i,low_ind(i),upp_ind(i) - end do - ! print *,i,low_ind(i),upp_ind(i) - end if - -!$OMP PARALLELDO DEFAULT(NONE) & -!$OMP SHARED( n_threads, low_ind, upp_ind, Id_loc, & -!$OMP lon, lat, m09_lon, m09_lat, tid_m09) & -!$OMP PRIVATE(n,i,t_count,min_lon, max_lon, min_lat, max_lat, & -!$OMP sub_tid, nplus, sub_lon, sub_lat, rev_dist, dw, & -!$OMP tile_found, mask) - - DO t_count = 1,n_threads - - allocate (mask (1: InNTILES)) - - OUT_TILES : do n = low_ind(t_count),upp_ind(t_count) - if(MOD(n,10000) == 0) print *,'In ID_LOC', t_count,n - dw = 0.25 - - ZOOMOUT : do - - tile_found = .false. - - ! Min/Max lon/lat of the working window - ! ------------------------------------- - - min_lon = MAX(lon (n) - dw, -180.) - max_lon = MIN(lon (n) + dw, 180.) - min_lat = MAX(lat (n) - dw, -90.) - max_lat = MIN(lat (n) + dw, 90.) - - mask = .false. - mask = ((m09_lat >= min_lat .and. m09_lat <= max_lat).and.(m09_lon >= min_lon .and. m09_lon <= max_lon)) - nplus = count(mask = mask) - - if(nplus < 0) then - dw = dw + 0.5 - CYCLE - endif - - allocate (sub_tid (1:nplus)) - allocate (sub_lon (1:nplus)) - allocate (sub_lat (1:nplus)) - allocate (rev_dist(1:nplus)) - - sub_tid = PACK (tid_m09 , mask= mask) - sub_lon = PACK (m09_lon , mask= mask) - sub_lat = PACK (m09_lat , mask= mask) - - ! compute distance from the tile - - sub_lat = sub_lat * MAPL_PI/180. - sub_lon = sub_lon * MAPL_PI/180. - - SEEK : if(Id_loc(n) < 0) then - - rev_dist = 1.e20 - - do i = 1,nplus - - rev_dist(i) = haversine(to_radian(lat(n)), to_radian(lon(n)), & - sub_lat(i), sub_lon(i)) - - end do - - FOUND : if(minval (rev_dist) < 1.e19) then - Id_loc(n) = sub_tid(minloc(rev_dist,1)) - tile_found = .true. - - if(Id_loc(n) ==0) then - print *, rev_dist - print *, sub_tid - print *, minval (rev_dist) - print *, minloc(rev_dist,1) - stop - endif - endif FOUND - - endif SEEK - - deallocate (sub_tid, sub_lon, sub_lat, rev_dist) - - if(tile_found) GO TO 100 - - ! if not increase the window size - dw = dw + 0.25 - - end do ZOOMOUT - -100 continue - - END do OUT_TILES - - deallocate (mask) - - end DO ! PARALLEL - -!$OMP ENDPARALLELDO - - END SUBROUTINE get_id_loc - - ! ***************************************************************************** - - function to_radian(degree) result(rad) - - ! degrees to radians - real,intent(in) :: degree - real :: rad - - rad = degree*MAPL_PI/180. - - end function to_radian - - ! ------------------------------------------------------------------------------------------------- - - SUBROUTINE regrid_alb (NTILES, id_loc) - - implicit none - integer, intent (in) :: NTILES - integer, dimension (NTILES), intent (in) :: id_loc - character*10 :: string - integer :: STATUS, ncid, NCOutID, t, time_slice, time_slice_next, yr, mn, dd, yr1, mn1, dd1, n_tslices - character (len=4), dimension (:), allocatable :: MMDD, MMDD_next - real, allocatable, dimension (:) :: varin, varout - - n_tslices = NOCTAD - - status = NF_OPEN('data/CATCH/MODIS-Albedo2/MCD43GF_wsa_H11V13.nc',NF_NOWRITE, ncid); VERIFY_(STATUS) - status = NF_INQ_DIM (ncid,3,string, n_tslices); VERIFY_(STATUS) - allocate (MMDD (0: n_tslices + 1)) - allocate (MMDD_next (0: n_tslices + 1)) - status = NF_GET_VARA_text(ncid, 3,(/1,1/),(/4,n_tslices/),MMDD(1:n_tslices)); VERIFY_(STATUS) - status = NF_CLOSE(ncid); VERIFY_(STATUS) - - mmdd(0) = mmdd(n_tslices) - mmdd(n_tslices + 1)= mmdd(1) - mmdd_next(0:n_tslices - 1) = mmdd(1:n_tslices) - mmdd_next(n_tslices: n_tslices + 1) = mmdd (1:2) - - allocate (varin (1:InNTILES)) - allocate (varout(1:NTILES)) - - STATUS = NF_OPEN('data/CATCH/CATCHCN_fPAR_Alb_stats2.nc4', NF_NOWRITE,NCOutID) ; VERIFY_(STATUS) - open (10, file = 'clsm/MODISVISmean.dat' ,form='unformatted',status='unknown',convert='little_endian') - open (11, file = 'clsm/MODISNIRmean.dat' ,form='unformatted',status='unknown',convert='little_endian') - open (12, file = 'clsm/MODISVISstd.dat' ,form='unformatted',status='unknown',convert='little_endian') - open (13, file = 'clsm/MODISNIRstd.dat' ,form='unformatted',status='unknown',convert='little_endian') - open (14, file = 'clsm/MODELFPARmean.dat' ,form='unformatted',status='unknown',convert='little_endian') - open (15, file = 'clsm/MODELFPARstd.dat' ,form='unformatted',status='unknown',convert='little_endian') - open (16, file = 'clsm/MODISFPARmean.dat' ,form='unformatted',status='unknown',convert='little_endian') - open (17, file = 'clsm/MODISFPARstd.dat' ,form='unformatted',status='unknown',convert='little_endian') - - do t =0,n_tslices+1 - - time_slice = t - yr = 1 - yr1= 1 - if(t == 0) then - time_slice = n_tslices - yr = 1 - 1 - endif - - if(t >= n_tslices) then - yr1 = 1 + 1 - if(t ==n_tslices + 1) then - time_slice = 1 - yr = 1 + 1 - endif - endif - - read(mmdd(t),'(i2.2,i2.2)') mn,dd - read(mmdd_next(t),'(i2.2,i2.2)') mn1,dd1 - - STATUS = NF_GET_VARA_REAL(NCOutID,VARID(NCOutID,'MODISVISmean' ),(/1,time_slice/), (/InNTILES,1/), varin ) ; VERIFY_(STATUS) - varout = varin (id_loc) - write(10) float((/yr,mn,dd,0,0,0,yr1,mn1,dd1,0,0,0,ntiles,1/)) - write(11) float((/yr,mn,dd,0,0,0,yr1,mn1,dd1,0,0,0,ntiles,1/)) - write(12) float((/yr,mn,dd,0,0,0,yr1,mn1,dd1,0,0,0,ntiles,1/)) - write(13) float((/yr,mn,dd,0,0,0,yr1,mn1,dd1,0,0,0,ntiles,1/)) - write(14) float((/yr,mn,dd,0,0,0,yr1,mn1,dd1,0,0,0,ntiles,1/)) - write(15) float((/yr,mn,dd,0,0,0,yr1,mn1,dd1,0,0,0,ntiles,1/)) - write(16) float((/yr,mn,dd,0,0,0,yr1,mn1,dd1,0,0,0,ntiles,1/)) - write(17) float((/yr,mn,dd,0,0,0,yr1,mn1,dd1,0,0,0,ntiles,1/)) - - STATUS = NF_GET_VARA_REAL(NCOutID,VARID(NCOutID,'MODISVISmean' ),(/1,time_slice/), (/InNTILES,1/), varin ) ; VERIFY_(STATUS) - varout = varin (id_loc) - write (10) varout - - STATUS = NF_GET_VARA_REAL(NCOutID,VARID(NCOutID,'MODISNIRmean' ),(/1,time_slice/), (/InNTILES,1/), varin ) ; VERIFY_(STATUS) - varout = varin (id_loc) - write (11) varout - - STATUS = NF_GET_VARA_REAL(NCOutID,VARID(NCOutID,'MODISVISstd' ),(/1,time_slice/) , (/InNTILES,1/), varin ) ; VERIFY_(STATUS) - varout = varin (id_loc) - write (12) varout - - STATUS = NF_GET_VARA_REAL(NCOutID,VARID(NCOutID,'MODISNIRstd' ),(/1,time_slice/) , (/InNTILES,1/), varin ) ; VERIFY_(STATUS) - varout = varin (id_loc) - write (13) varout - - STATUS = NF_GET_VARA_REAL(NCOutID,VARID(NCOutID,'MODELFPARmean' ),(/1,time_slice/), (/InNTILES,1/), varin ) ; VERIFY_(STATUS) - varout = varin (id_loc) - write (14) varout - - STATUS = NF_GET_VARA_REAL(NCOutID,VARID(NCOutID,'MODELFPARstd' ),(/1,time_slice/) , (/InNTILES,1/), varin ) ; VERIFY_(STATUS) - varout = varin (id_loc) - write (15) varout - - STATUS = NF_GET_VARA_REAL(NCOutID,VARID(NCOutID,'MODISFPARmean' ),(/1,time_slice/), (/InNTILES,1/), varin ) ; VERIFY_(STATUS) - varout = varin (id_loc) - write (16) varout - - STATUS = NF_GET_VARA_REAL(NCOutID,VARID(NCOutID,'MODISFPARstd' ),(/1,time_slice/) , (/InNTILES,1/), varin ) ; VERIFY_(STATUS) - varout = varin (id_loc) - write (17) varout - end do - - deallocate (varin, varout) - - close (10 , status = 'keep') - close (11 , status = 'keep') - close (12 , status = 'keep') - close (13 , status = 'keep') - close (14 , status = 'keep') - close (15 , status = 'keep') - close (16 , status = 'keep') - close (17 , status = 'keep') - - END SUBROUTINE regrid_alb - - ! ------------------------------------------------------------------------------------------------- - - SUBROUTINE preprocess_m09 - - implicit none - logical :: file_exists - INTEGER :: NT, ND, DAY, year, STATUS, NCOutID, k - CHARACTER*8 :: YYYYMMDD - CHARACTER*6 :: YYYYMM - CHARACTER*4 :: YYYY - CHARACTER*2 :: MM, DD - real :: yr,mn,dy,dum,yr1,mn1,dy1 - - real, allocatable, dimension (:,:) :: MODIS_VISDF, MODIS_NIRDF, MODEL_fPAR - real, allocatable, dimension (:) :: data_read, data_save, & - MODISVISmean, MODISNIRmean, MODISVISstd, MODISNIRstd, MODELFPARmean, MODELFPARstd, & - MODISFPARmean, MODISFPARstd - integer, allocatable, dimension (:) :: ldas2bcs - type(date_time_type), dimension (YearE - YearB + 1) :: octad_time - - - inquire(file='data/CATCH/CATCHCN_fPAR_Alb_stats2.nc4', exist=file_exists) - if(.not. file_exists) call create_stat_file - -! open (99,file='clsm/comp_CATCHCN_AlbScale_parameters.log', form ='formatted', action='write', status= 'unknown') - - STATUS = NF_OPEN ('data/CATCH/CATCHCN_fPAR_Alb_stats2.nc4', NF_WRITE,NCOutID) ; VERIFY_(STATUS) - - allocate (MODIS_VISDF (1:InNTILES, yearE - yearB + 1)) - allocate (MODIS_NIRDF (1:InNTILES, yearE - yearB + 1)) - allocate (MODEL_fPAR (1:InNTILES, yearE - yearB + 1)) - allocate (ldas2bcs (1:InNTILES)) - allocate (data_read (1:InNTILES)) - allocate (data_save (1:InNTILES)) - allocate (MODISVISmean (1:InNTILES)) - allocate (MODISNIRmean (1:InNTILES)) - allocate (MODISVISstd (1:InNTILES)) - allocate (MODISNIRstd (1:InNTILES)) - allocate (MODELFPARmean(1:InNTILES)) - allocate (MODELFPARstd (1:InNTILES)) - allocate (MODISFPARmean(1:InNTILES)) - allocate (MODISFPARstd (1:InNTILES)) - - open (10,file =trim(EXPDIR)//'rc_out/'//trim(EXNAME)//'.ldas_tilecoord.bin',status='old',form='unformatted',convert='big_endian') - read (10) k - read (10) LDAS2BCS - close(10, status = 'keep') - - OPEN_FILES1 : DO year = YearB, YearE - - write (YYYY ,'(i4.4)') year - open (10 + year - yearB, file = trim(InBCSDIR)//'MODIS6/'//YYYY//'//visdf.dat', form = 'unformatted', action = 'read') - open (30 + year - yearB, file = trim(InBCSDIR)//'MODIS6/'//YYYY//'//nirdf.dat', form = 'unformatted', action = 'read') -! WRITE (99,*)10 + year - yearB, YYYY//'//visdf.dat' -! WRITE (99,*)30 + year - yearB, YYYY//'//nirdf.dat' -! WRITE (99,*) ' ' - WRITE (*,*)10 + year - yearB, YYYY//'//visdf.dat' - WRITE (*,*)30 + year - yearB, YYYY//'//nirdf.dat' - WRITE (*,*) ' ' - octad_time(year - yearB + 1)%year = year - 1 - octad_time(year - yearB + 1)%month = 12 - octad_time(year - yearB + 1)%day = 31 - octad_time(year - yearB + 1)%hour = 0 - octad_time(year - yearB + 1)%min = 0 - octad_time(year - yearB + 1)%sec = 0 - - END DO OPEN_FILES1 - - ND = 8 - - OCTAD_LOOP : DO NT = 1, NOCTAD - - ! BEGIN READING MODIS VISDF/NIRDF - - MODIS_VISDF = 0. - MODIS_NIRDF = 0. - MODEL_fPAR = 0. - - if(NT == NOCTAD) ND = 5 -! WRITE (99,*) NT, ND, yearB, yearE - print *, NT, ND, yearB, yearE - - READ_YEARS : DO year = YearB, YearE - - read (10 + year - yearB) modis_visdf (:, year - yearB + 1) - read (30 + year - yearB) modis_nirdf (:, year - yearB + 1) - - DAILY_LOOP : DO day = 1,ND - - call augment_date_time(86400, octad_time(year - yearB + 1)) - - write (YYYY, '(i4.4)') octad_time(year - yearB + 1)%year - write (MM , '(i2.2)') octad_time(year - yearB + 1)%month - write (DD , '(i2.2)') octad_time(year - yearB + 1)%day - - YYYYMMDD = YYYY//MM//DD -! WRITE (99,*) trim(EXNAME)//'.ens_avg.ldas_tile_daily_out.'//YYYYMMDD//'.bin' - print *, trim(EXNAME)//'.ens_avg.ldas_tile_daily_out.'//YYYYMMDD//'.bin' - open (60, file = trim(EXPDIR)//'cat/ens_avg/Y'//YYYY//'/M'//MM//'/'// & - trim(EXNAME)//'.ens_avg.ldas_tile_daily_out.'//YYYYMMDD//'.bin', & - form = 'unformatted', convert='big_endian', action = 'read') - - do k = 1,3 - read (60) data_read - if(k == 2) data_save = data_read - end do - - MODEL_fPAR (:,year - yearB + 1) = MODEL_fPAR (:,year - yearB + 1) + data_save / (data_read + 1.e-20)/ real (ND) - close (60, status = 'keep') - - END DO DAILY_LOOP - - ! reoder to the order of BCs - - data_read = MODEL_fPAR (:,year - yearB + 1) - - do k = 1, InNTILES - MODEL_fPAR (LDAS2BCS(k),year - yearB + 1) = data_read (k) - end do - END DO READ_YEARS - - ! COMPUTE STATS - - CALL compute_stats (MODIS_VISDF, MODIS_NIRDF, MODEL_fPAR, & - MODISVISmean, MODISNIRmean, MODISVISstd, MODISNIRstd, MODELFPARmean, MODELFPARstd) - - STATUS = NF_PUT_VARA_REAL(NCOutID,VARID(NCOutID,'MODISVISmean' ),(/1,NT/), (/InNTILES,1/), MODISVISmean ) ; VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID,VARID(NCOutID,'MODISNIRmean' ),(/1,NT/), (/InNTILES,1/), MODISNIRmean ) ; VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID,VARID(NCOutID,'MODISVISstd' ),(/1,NT/), (/InNTILES,1/), MODISVISstd ) ; VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID,VARID(NCOutID,'MODISNIRstd' ),(/1,NT/), (/InNTILES,1/), MODISNIRstd ) ; VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID,VARID(NCOutID,'MODELFPARmean' ),(/1,NT/), (/InNTILES,1/), MODELFPARmean) ; VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID,VARID(NCOutID,'MODELFPARstd' ),(/1,NT/), (/InNTILES,1/), MODELFPARstd ) ; VERIFY_(STATUS) - - END DO OCTAD_LOOP - - CLOSE_FILES1 : DO year = YearB, YearE - - close (10 + year - yearB, status = 'keep') - close (30 + year - yearB, status = 'keep') - - END DO CLOSE_FILES1 - - ! MODIS FPAR - ! ---------- - - deallocate (MODIS_VISDF, MODIS_NIRDF, MODEL_fPAR) - allocate (MODIS_VISDF (1:InNTILES, yearE - yearB1 + 1)) - allocate (MODIS_NIRDF (1:InNTILES, yearE - yearB1 + 1)) - allocate (MODEL_fPAR (1:InNTILES, yearE - yearB1 + 1)) - - OPEN_FILES2 : DO year = YearB1, YearE - - write (YYYY ,'(i4.4)') year - open (10 + year - yearB1, file = trim(InBCSDIR)//'MODIS6/'//YYYY//'//visdf.dat', form = 'unformatted', action = 'read') - open (30 + year - yearB1, file = trim(InBCSDIR)//'MODIS6/'//YYYY//'//nirdf.dat', form = 'unformatted', action = 'read') - open (50 + year - yearB1, file = trim(InBCSDIR)//'MODIS6/'//YYYY//'//fpar.dat', form = 'unformatted', action = 'read') -! WRITE (99,*) " " -! WRITE (99,*) "MODIS FPAR" -! WRITE (99,*) "==========" -! WRITE (99,*) " " - -! WRITE (99,*)10 + year - yearB1, YYYY//'//visdf.dat' -! WRITE (99,*)30 + year - yearB1, YYYY//'//nirdf.dat' -! WRITE (99,*)50 + year - yearB1, YYYY//'//fpar.dat' -! WRITE (99,*) ' ' - WRITE (*,*)10 + year - yearB1, YYYY//'//visdf.dat' - WRITE (*,*)30 + year - yearB1, YYYY//'//nirdf.dat' - WRITE (*,*)50 + year - yearB1, YYYY//'//fpar.dat' - WRITE (*,*) ' ' - read (50 + year - yearB1) yr,mn,dy,dum,dum,dum,yr1,mn1,dy1 - WRITE (*,*) yr,mn,dy,yr1,mn1,dy1 - read (50 + year - yearB1) MODEL_FPAR (:, year - yearB1 + 1) - - octad_time(year - yearB1 + 1)%year = year - 1 - octad_time(year - yearB1 + 1)%month = 12 - octad_time(year - yearB1 + 1)%day = 31 - octad_time(year - yearB1 + 1)%hour = 0 - octad_time(year - yearB1 + 1)%min = 0 - octad_time(year - yearB1 + 1)%sec = 0 - - END DO OPEN_FILES2 - - ND = 8 - - OCTAD_LOOP2 : DO NT = 1, NOCTAD - - ! BEGIN READING MODIS VISDF/NIRDF/FPAR - - MODIS_VISDF = 0. - MODIS_NIRDF = 0. - MODEL_fPAR = 0. - - if(NT == NOCTAD) ND = 5 -! WRITE (99,*) NT, ND, yearB1, yearE - print *, NT, ND, yearB1, yearE - READ_YEARS2 : DO year = YearB1, YearE - - read (10 + year - yearB1) modis_visdf (:, year - yearB1 + 1) - read (30 + year - yearB1) modis_nirdf (:, year - yearB1 + 1) - read (50 + year - yearB1) yr,mn,dy,dum,dum,dum,yr1,mn1,dy1 - WRITE (*,*) yr,mn,dy,dum,dum,dum,yr1,mn1,dy1 - read (50 + year - yearB1) MODEL_FPAR (:, year - yearB1 + 1) - END DO READ_YEARS2 - - ! COMPUTE STATS - - CALL compute_stats (MODIS_VISDF, MODIS_NIRDF, MODEL_fPAR, & - MODISVISmean, MODISNIRmean, MODISVISstd, MODISNIRstd, MODELFPARmean, MODELFPARstd) - - STATUS = NF_PUT_VARA_REAL(NCOutID,VARID(NCOutID,'MODISFPARmean' ),(/1,NT/), (/InNTILES,1/), MODELFPARmean) ; VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCOutID,VARID(NCOutID,'MODISFPARstd' ),(/1,NT/), (/InNTILES,1/), MODELFPARstd ) ; VERIFY_(STATUS) - - END DO OCTAD_LOOP2 - - CLOSE_FILES2 : DO year = YearB1, YearE - - close (10 + year - yearB1, status = 'keep') - close (30 + year - yearB1, status = 'keep') - close (50 + year - yearB1, status = 'keep') - - END DO CLOSE_FILES2 - - STATUS = NF_CLOSE (NCOutID ) - - END SUBROUTINE preprocess_m09 - - ! ----------------------------------------------------------------- - - SUBROUTINE create_stat_file - - implicit none - - integer :: NCFOutID, STATUS, vid, tid, lid, n, k - character (22) :: time_stamp, tmpstr - integer, dimension(8) :: date_time_values - real, dimension (:), allocatable :: lons, lats - - STATUS = NF_CREATE ('data/CATCH/CATCHCN_fPAR_Alb_stats2.nc4', NF_NETCDF4, NCFOutID );VERIFY_(STATUS) - STATUS = NF_DEF_DIM(NCFOutID, 'octad', NOCTAD, TID) ; VERIFY_(STATUS) - STATUS = NF_DEF_DIM(NCFOutID, 'tiles', InNTILES, LID) ; VERIFY_(STATUS) - STATUS = NF_DEF_VAR(NCFOutID, 'lon' , NF_FLOAT ,1 ,(/LID/), vid);VERIFY_(STATUS) - STATUS = NF_DEF_VAR(NCFOutID, 'lat' , NF_FLOAT ,1 ,(/LID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODISVISmean', NF_FLOAT ,2 ,(/LID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODISNIRmean', NF_FLOAT ,2 ,(/LID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODISVISstd' , NF_FLOAT ,2 ,(/LID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODISNIRstd' , NF_FLOAT ,2 ,(/LID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODELFPARmean',NF_FLOAT ,2 ,(/LID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODELFPARstd' ,NF_FLOAT ,2 ,(/LID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODISFPARmean',NF_FLOAT ,2 ,(/LID,TID/), vid);VERIFY_(STATUS) - status = NF_DEF_VAR(NCFOutID, 'MODISFPARstd' ,NF_FLOAT ,2 ,(/LID,TID/), vid);VERIFY_(STATUS) - ! Global attributes - - call date_and_time(VALUES=date_time_values) - - write (time_stamp,'(i4.4,a1,i2.2,a1,i2.2,1x,a2,1x,i2.2,a1,i2.2,a1,i2.2)') & - date_time_values(1),'-',date_time_values(2),'-',date_time_values(3),'at', & - date_time_values(5),':',date_time_values(6),':',date_time_values(7) - - status = NF_PUT_ATT_TEXT(NCFOutID, NF_GLOBAL, 'CreatedBy', LEN_TRIM("Sarith Mahanama"), & - trim("Sarith Mahanama")) - status = NF_PUT_ATT_TEXT(NCFOutID, NF_GLOBAL, 'Date' , LEN_TRIM(time_stamp),trim(time_stamp)) - - status = NF_ENDDEF(NCFOutID ) - - ! Read and put lat/lon data - - allocate (lons (1:InNTILES)) - allocate (lats (1:InNTILES)) - - open (10, file = trim(InBCSDIR)//trim(InGFILE)//'.til', form = 'formatted', status = 'old') - - do n = 1,8 - read (10,*) tmpstr - end do - - do n = 1, InNTILES - read (10,*) k, k, lons(n), lats(n) - end do - - close (10, status = 'keep') - - STATUS = NF_PUT_VARA_REAL(NCFOutID,VARID(NCFOutID,'lon' ),(/1/), (/InNTILES/), lons) ; VERIFY_(STATUS) - STATUS = NF_PUT_VARA_REAL(NCFOutID,VARID(NCFOutID,'lat' ),(/1/), (/InNTILES/), lats) ; VERIFY_(STATUS) - - STATUS = NF_CLOSE (NCFOutID ) - - END SUBROUTINE create_stat_file - - ! ---------------------------------------------------------------------- - - SUBROUTINE compute_stats (MODIS_VISDF, MODIS_NIRDF, MODEL_fPAR, & - MODISVISmean, MODISNIRmean, MODISVISstd, MODISNIRstd, MODELFPARmean, MODELFPARstd) - - implicit none - real, dimension (:,:), intent (in) :: MODIS_VISDF, MODIS_NIRDF, MODEL_fPAR - real, dimension (:) , intent (inout) :: MODISVISmean, MODISNIRmean, MODISVISstd, & - MODISNIRstd, MODELFPARmean, MODELFPARstd - integer :: NX, NY, N, t - REAL :: MF, MV, MN, SF, SV, SN, ZV, ZN, CV, CN, CF - - NX = size (MODIS_VISDF,1) - NY = size (MODIS_VISDF,2) - print *,'Entered compute_stats', NX, NY - if (NX /= InNTILES) then - print *, 'NX NTILLES MISMAATCH : ', InNTILES, NX, NY - STOP - ENDIF - - DO N = 1, NX - -! MF = SUM (MODEL_fPAR (N,:)) / REAL (NY) -! MV = SUM (MODIS_VISDF (N,:)) / REAL (NY) -! MN = SUM (MODIS_nirDF (N,:)) / REAL (NY) - MF = 0. - MV = 0. - MN = 0. - CV = 0. - CN = 0. - CF = 0. - - do T = 1, NY - if ((MODEL_fPAR (N,T) >= 0.).AND.(MODEL_fPAR (N,T) <= 1.)) then - MF = MF + MODEL_fPAR (N,T) - CF = CF + 1. - endif - if ((MODIS_VISDF (N,T) >= 0.).AND.(MODIS_VISDF (N,T) <= 1.)) then - MV = MV + MODIS_VISDF (N,T) - CV = CV + 1. - endif - if ((MODIS_NIRDF (N,T) >= 0.).AND.(MODIS_NIRDF (N,T) <= 1.)) then - MN = MN + MODIS_NIRDF (N,T) - CN = CN + 1. - endif - end do - - IF(CF > 0) MF = MF / CF - IF(CV > 0) MV = MV / CV - IF(CN > 0) MN = MN / CN - - ! STANDARD DEVIATION - - SF = 1.e-15 - SV = 1.e-15 - SN = 1.e-15 - - do T = 1, NY - if ((MODEL_fPAR (N,T) >= 0.).AND.(MODEL_fPAR (N,T) <= 1.)) SF = SF + (MODEL_fPAR (N,t) - MF)*(MODEL_fPAR (N,t) - MF) - if ((MODIS_VISDF (N,T) >= 0.).AND.(MODIS_VISDF (N,T) <= 1.)) SV = SV + (MODIS_VISDF (N,t) - MV)*(MODIS_VISDF (N,t) - MV) - if ((MODIS_NIRDF (N,T) >= 0.).AND.(MODIS_NIRDF (N,T) <= 1.)) SN = SN + (MODIS_NIRDF (N,t) - MN)*(MODIS_NIRDF (N,t) - MN) - end do - - IF(CF > 0) SF = SQRT (SF / CF) - IF(CV > 0) SV = SQRT (SV / CV) - IF(CN > 0) SN = SQRT (SN / CN) - - ! CORRELATION - - ZV = 0. - ZN = 0. - - DO T = 1, NY - ZV = ZV + (MODEL_fPAR (N,t) - MF)*(MODIS_VISDF (N,t) - MV)/SF/SV - ZN = ZN + (MODEL_fPAR (N,t) - MF)*(MODIS_NIRDF (N,t) - MN)/SF/SV - END DO - - ZV = ZV / REAL (NY) - ZN = ZN / REAL (NY) - - MODISVISmean (N) = MV - MODISNIRmean (N) = MN - MODISVISstd (N) = SV - MODISNIRstd (N) = SN - MODELFPARmean(N) = MF - MODELFPARstd (N) = SF - - if(ZV < 0.) MODISVISstd (N) = -1. * MODISVISstd (N) - if(Zn < 0.) MODISnirstd (N) = -1. * MODISnirstd (N) - - END DO - - print *,'Leaving compute_stats' - - END SUBROUTINE compute_stats - - ! ---------------------------------------------------------------------- - - integer function VarID (NCFID, VNAME) - - integer, intent (in) :: NCFID - character(*), intent (in) :: VNAME - integer :: status - - STATUS = NF_INQ_VARID (NCFID, trim(VNAME) ,VarID) - IF (STATUS .NE. NF_NOERR) & - CALL HANDLE_ERR(STATUS, trim(VNAME)) - - end function VarID - - ! ----------------------------------------------------------------------- - - SUBROUTINE HANDLE_ERR(STATUS, Line) - - INTEGER, INTENT (IN) :: STATUS - CHARACTER(*), INTENT (IN) :: Line - - IF (STATUS .NE. NF_NOERR) THEN - PRINT *, trim(Line),': ',STATUS, NF_STRERROR(STATUS) - STOP 'Stopped' - ENDIF - - END SUBROUTINE HANDLE_ERR - - ! ***************************************************************************** - - subroutine ReadCNTilFile (InCNTileFile, nt, xlon, xlat) - - implicit none - character(*), intent (in) :: InCNTileFile - integer , intent (in) :: nt - real, dimension (nt), intent(inout) :: xlon, xlat - integer :: n,icnt,ityp - real :: xval,yval, pf - - open(11,file=InCNTileFile, & - form='formatted',action='read',status='old') - - do n = 1,8 ! skip header - read(11,*) - end do - - icnt = 0 - ityp = 100 - - do while (ityp == 100) ! loop over land tiles - read(11,*) ityp,pf,xval,yval - if(ityp == 100) then - icnt = icnt + 1 - xlon(icnt) = xval - xlat(icnt) = yval - endif - end do - - close(11) - - end subroutine ReadCNTilFile - - ! ***************************************************************************** - - real function haversine(deglat1,deglon1,deglat2,deglon2) - ! great circle distance -- adapted from Matlab - real,intent(in) :: deglat1,deglon1,deglat2,deglon2 - real :: a,c, dlat,dlon,lat1,lat2 - real,parameter :: radius = 6371.0E3 - -! dlat = to_radian(deglat2-deglat1) -! dlon = to_radian(deglon2-deglon1) - ! lat1 = to_radian(deglat1) -! lat2 = to_radian(deglat2) - dlat = deglat2-deglat1 - dlon = deglon2-deglon1 - lat1 = deglat1 - lat2 = deglat2 - a = (sin(dlat/2))**2 + cos(lat1)*cos(lat2)*(sin(dlon/2))**2 - if(a>=0. .and. a<=1.) then - c = 2*atan2(sqrt(a),sqrt(1-a)) - haversine = radius*c / 1000. - else - haversine = 1.e20 - endif - end function - - ! ***************************************************************************** - - END MODULE comp_CATCHCN_AlbScale_parameters 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 e21883729..e614afbf1 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs @@ -780,14 +780,6 @@ cd clsm.${IM}x${JM} /bin/mv green.dat green_clim_${RS}_DC.data /bin/mv lnfm.dat lnfm_clim_${RS}_DC.data /bin/mv ndvi.dat ndvi_clim_${RS}_DC.data - /bin/mv MODELFPARmean.dat MODELFPARmean_${RS}_DC.dat - /bin/mv MODELFPARstd.dat MODELFPARstd_${RS}_DC.dat - /bin/mv MODISFPARmean.dat MODISFPARmean_${RS}_DC.dat - /bin/mv MODISFPARstd.dat MODISFPARstd_${RS}_DC.dat - /bin/mv MODISNIRmean.dat MODISNIRmean_${RS}_DC.dat - /bin/mv MODISNIRstd.dat MODISNIRstd_${RS}_DC.dat - /bin/mv MODISVISmean.dat MODISVISmean_${RS}_DC.dat - /bin/mv MODISVISstd.dat MODISVISstd_${RS}_DC.dat /bin/rm -f sedfile cat > sedfile << EOF @@ -1010,14 +1002,6 @@ cd clsm.C${NC} /bin/mv green.dat green_clim_${RC}.data /bin/mv lnfm.dat lnfm_clim_${RC}.data /bin/mv ndvi.dat ndvi_clim_${RC}.data - /bin/mv MODELFPARmean.dat MODELFPARmean_${RC}.dat - /bin/mv MODELFPARstd.dat MODELFPARstd_${RC}.dat - /bin/mv MODISFPARmean.dat MODISFPARmean_${RC}.dat - /bin/mv MODISFPARstd.dat MODISFPARstd_${RC}.dat - /bin/mv MODISNIRmean.dat MODISNIRmean_${RC}.dat - /bin/mv MODISNIRstd.dat MODISNIRstd_${RC}.dat - /bin/mv MODISVISmean.dat MODISVISmean_${RC}.dat - /bin/mv MODISVISstd.dat MODISVISstd_${RC}.dat /bin/rm -f sedfile if( $CUBED_SPHERE_OCEAN == TRUE ) then @@ -1175,14 +1159,6 @@ cd clsm.C${NC} /bin/mv green.dat green_clim_${RC}.data /bin/mv lnfm.dat lnfm_clim_${RC}.data /bin/mv ndvi.dat ndvi_clim_${RC}.data - /bin/mv MODELFPARmean.dat MODELFPARmean_${RC}.dat - /bin/mv MODELFPARstd.dat MODELFPARstd_${RC}.dat - /bin/mv MODISFPARmean.dat MODISFPARmean_${RC}.dat - /bin/mv MODISFPARstd.dat MODISFPARstd_${RC}.dat - /bin/mv MODISNIRmean.dat MODISNIRmean_${RC}.dat - /bin/mv MODISNIRstd.dat MODISNIRstd_${RC}.dat - /bin/mv MODISVISmean.dat MODISVISmean_${RC}.dat - /bin/mv MODISVISstd.dat MODISVISstd_${RC}.dat /bin/rm -f sedfile if( $CUBED_SPHERE_OCEAN == TRUE ) then @@ -1357,14 +1333,6 @@ cd clsm.${IM}x${JM} /bin/mv green.dat green_clim_${RS}_DE.data /bin/mv lnfm.dat lnfm_clim_${RS}_DE.data /bin/mv ndvi.dat ndvi_clim_${RS}_DE.data - /bin/mv MODELFPARmean.dat MODELFPARmean_${RS}_DE.dat - /bin/mv MODELFPARstd.dat MODELFPARstd_${RS}_DE.dat - /bin/mv MODISFPARmean.dat MODISFPARmean_${RS}_DE.dat - /bin/mv MODISFPARstd.dat MODISFPARstd_${RS}_DE.dat - /bin/mv MODISNIRmean.dat MODISNIRmean_${RS}_DE.dat - /bin/mv MODISNIRstd.dat MODISNIRstd_${RS}_DE.dat - /bin/mv MODISVISmean.dat MODISVISmean_${RS}_DE.dat - /bin/mv MODISVISstd.dat MODISVISstd_${RS}_DE.dat cd ../ 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 877a3c7ec..7a32853c7 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -23,7 +23,6 @@ PROGRAM mkCatchParam use rmTinyCatchParaMod use process_hres_data - use comp_CATCHCN_AlbScale_parameters, ONLY : albedo4catchcn ! use module_irrig_params, ONLY : create_irrig_params implicit none @@ -665,7 +664,6 @@ PROGRAM mkCatchParam ! if (.not.file_exists) call create_irrig_params (nc,nr,gridnamer) ! write (log_file,'(a)')'Done computing irrigation model parameters ...............13' - ! call albedo4catchcn (gridnamet) write (log_file,'(a)')'============================================================' write (log_file,'(a)')'DONE creating CLSM data files...............................'