text:
c     Subroutine to calc basic stats for assembled vectors
c     AJ_Kettle, Oct26/2017

      SUBROUTINE calc_basic_stats_sd1(l_mlent,j_sd,j_day,f_ndflag,
     +  ft_sd_pres_hpa,ft_sd_airt_c,ft_sd_wetb_c,
     +  ft_sd_vprs_hpa,ft_sd_relh_pc,ft_sd_relhreg_pc,
     +  ft_sd_wdir_code,ft_sd_wspd_bft,ft_sd_ccov_okta,
     +  ft_sd_visi_code,ft_sd_grdcnd_code,ft_sd_ppt_mm,

     +  ft_pres_dayavg_hpa,ft_airt_daymax_c,ft_airt_daymin_c,
     +  ft_airt_range_c,ft_boden_daymin_c,ft_airt_dayavg_c,
     +  ft_vprs_dayavg_hpa,ft_relh_dayavg_pc,ft_wspd_dayavg_bft,
     +  ft_ccov_dayavg_okta,ft_sundur_daytot_h,ft_ppt_daytot_mm,
     +  ft_snowcover_cm,ft_snow_daytot_cm,ft_wspd_daymax_ms,
     +  ft_snowmelted_cm,ft_we_snowmelted_mm,ft_we_snowcover_mm,

     +  ft_sd_pres_min_hpa,ft_sd_pres_max_hpa,ft_sd_pres_avg_hpa,
     +    ft_sd_pres_ngood,ft_sd_pres_nbad,
     +  ft_sd_airt_min_c,ft_sd_airt_max_c,ft_sd_airt_avg_c,
     +    ft_sd_airt_ngood,ft_sd_airt_nbad,
     +  ft_sd_wetb_min_c,ft_sd_wetb_max_c,ft_sd_wetb_avg_c,
     +    ft_sd_wetb_ngood,ft_sd_wetb_nbad,
     +  ft_sd_vprs_min_hpa,ft_sd_vprs_max_hpa,ft_sd_vprs_avg_hpa,
     +    ft_sd_vprs_ngood,ft_sd_vprs_nbad,
     +  ft_sd_relh_min_pc,ft_sd_relh_max_pc,ft_sd_relh_avg_pc,
     +    ft_sd_relh_ngood,ft_sd_relh_nbad,
     +  ft_sd_relhreg_min_pc,ft_sd_relhreg_max_pc,
     +    ft_sd_relhreg_avg_pc,
     +    ft_sd_relhreg_ngood,ft_sd_relhreg_nbad,
     +  ft_sd_wdir_min_code,ft_sd_wdir_max_code,ft_sd_wdir_avg_code,
     +    ft_sd_wdir_ngood,ft_sd_wdir_nbad,
     +  ft_sd_wspd_min_bft,ft_sd_wspd_max_bft,ft_sd_wspd_avg_bft,
     +    ft_sd_wspd_ngood,ft_sd_wspd_nbad,
     +  ft_sd_ccov_min_okta,ft_sd_ccov_max_okta,
     +    ft_sd_ccov_avg_okta,
     +    ft_sd_ccov_ngood,ft_sd_ccov_nbad,
     +  ft_sd_visi_min_code,ft_sd_visi_max_code,
     +    ft_sd_visi_avg_code,
     +    ft_sd_visi_ngood,ft_sd_visi_nbad,
     +  ft_sd_grdcnd_min_code,ft_sd_grdcnd_max_code,
     +    ft_sd_grdcnd_avg_code,
     +    ft_sd_grdcnd_ngood,ft_sd_grdcnd_nbad,
     +  ft_sd_ppt_min_mm,ft_sd_ppt_max_mm,ft_sd_ppt_avg_mm,
     +    ft_sd_ppt_ngood,ft_sd_ppt_nbad,

     +  ft_day_pres_min_hpa,ft_day_pres_max_hpa,
     +    ft_day_pres_avg_hpa,
     +    ft_day_pres_ngood,ft_day_pres_nbad,
     +  ft_day_airt_daymax_min_c,ft_day_airt_daymax_max_c,
     +    ft_day_airt_daymax_avg_c,
     +    ft_day_airt_daymax_ngood,ft_day_airt_daymax_nbad,
     +  ft_day_airt_daymin_min_c,ft_day_airt_daymin_max_c,
     +    ft_day_airt_daymin_avg_c,
     +    ft_day_airt_daymin_ngood,ft_day_airt_daymin_nbad,
     +  ft_day_airt_range_min_c,ft_day_airt_range_max_c,
     +    ft_day_airt_range_avg_c,
     +    ft_day_airt_range_ngood,ft_day_airt_range_nbad,
     +  ft_day_boden_daymin_min_c,ft_day_boden_daymin_max_c,
     +    ft_day_boden_daymin_avg_c,
     +    ft_day_boden_daymin_ngood,ft_day_boden_daymin_nbad,
     +  ft_day_airt_dayavg_min_c,ft_day_airt_dayavg_max_c,
     +    ft_day_airt_dayavg_avg_c,
     +    ft_day_airt_dayavg_ngood,ft_day_airt_dayavg_nbad,
     +  ft_day_vprs_dayavg_min_hpa,ft_day_vprs_dayavg_max_hpa,
     +    ft_day_vprs_dayavg_avg_hpa,
     +    ft_day_vprs_dayavg_ngood,ft_day_vprs_dayavg_nbad,
     +  ft_day_relh_dayavg_min_pc,ft_day_relh_dayavg_max_pc,
     +    ft_day_relh_dayavg_avg_pc,
     +    ft_day_relh_dayavg_ngood,ft_day_relh_dayavg_nbad,
     +  ft_day_wspd_dayavg_min_bft,ft_day_wspd_dayavg_max_bft,
     +    ft_day_wspd_dayavg_avg_bft,
     +    ft_day_wspd_dayavg_ngood,ft_day_wspd_dayavg_nbad,
     +  ft_day_ccov_dayavg_min_okta,ft_day_ccov_dayavg_max_okta,
     +    ft_day_ccov_dayavg_avg_okta,
     +    ft_day_ccov_dayavg_ngood,ft_day_ccov_dayavg_nbad,
     +  ft_day_sundur_daytot_min_h,ft_day_sundur_daytot_max_h,
     +    ft_day_sundur_daytot_avg_h,
     +    ft_day_sundur_daytot_ngood,ft_day_sundur_daytot_nbad,
     +  ft_day_ppt_daytot_min_mm,ft_day_ppt_daytot_max_mm,
     +    ft_day_ppt_daytot_avg_mm,
     +    ft_day_ppt_daytot_ngood,ft_day_ppt_daytot_nbad,
     +  ft_day_snowcover_min_cm,ft_day_snowcover_max_cm,
     +    ft_day_snowcover_avg_cm,
     +    ft_day_snowcover_ngood,ft_day_snowcover_nbad,
     +  ft_day_snow_daytot_min_cm,ft_day_snow_daytot_max_cm,
     +    ft_day_snow_daytot_avg_cm,
     +    ft_day_snow_daytot_ngood,ft_day_snow_daytot_nbad,
     +  ft_day_wspd_daymax_min_ms,ft_day_wspd_daymax_max_ms,
     +    ft_day_wspd_daymax_avg_ms,
     +    ft_day_wspd_daymax_ngood,ft_day_wspd_daymax_nbad,
     +  ft_day_snowmelted_min_cm,ft_day_snowmelted_max_cm,
     +    ft_day_snowmelted_avg_cm,
     +    ft_day_snowmelted_ngood,ft_day_snowmelted_nbad,
     +  ft_day_we_snowmelted_min_mm,ft_day_we_snowmelted_max_mm,
     +    ft_day_we_snowmelted_avg_mm,
     +    ft_day_we_snowmelted_ngood,ft_day_we_snowmelted_nbad,
     +  ft_day_we_snowcover_min_mm,ft_day_we_snowcover_max_mm,
     +    ft_day_we_snowcover_avg_mm,
     +    ft_day_we_snowcover_ngood,ft_day_we_snowcover_nbad)

      IMPLICIT NONE
c************************************************************************
      REAL                :: f_ndflag

      INTEGER             :: l_mlent
      INTEGER             :: j_sd
      REAL                :: ft_sd_pres_hpa(l_mlent)
      REAL                :: ft_sd_airt_c(l_mlent)
      REAL                :: ft_sd_wetb_c(l_mlent)
      REAL                :: ft_sd_vprs_hpa(l_mlent)
      REAL                :: ft_sd_relh_pc(l_mlent)
      REAL                :: ft_sd_relhreg_pc(l_mlent)
      REAL                :: ft_sd_wdir_code(l_mlent)
      REAL                :: ft_sd_wspd_bft(l_mlent)
      REAL                :: ft_sd_ccov_okta(l_mlent)
      REAL                :: ft_sd_visi_code(l_mlent)
      REAL                :: ft_sd_grdcnd_code(l_mlent)
      REAL                :: ft_sd_ppt_mm(l_mlent)

      INTEGER             :: j_day
      REAL                :: ft_pres_dayavg_hpa(l_mlent)
      REAL                :: ft_airt_daymax_c(l_mlent)
      REAL                :: ft_airt_daymin_c(l_mlent)
      REAL                :: ft_airt_range_c(l_mlent)
      REAL                :: ft_boden_daymin_c(l_mlent)
      REAL                :: ft_airt_dayavg_c(l_mlent)
      REAL                :: ft_vprs_dayavg_hpa(l_mlent)
      REAL                :: ft_relh_dayavg_pc(l_mlent)
      REAL                :: ft_wspd_dayavg_bft(l_mlent)
      REAL                :: ft_ccov_dayavg_okta(l_mlent)
      REAL                :: ft_sundur_daytot_h(l_mlent)
      REAL                :: ft_ppt_daytot_mm(l_mlent)
      REAL                :: ft_snowcover_cm(l_mlent)
      REAL                :: ft_snow_daytot_cm(l_mlent)
      REAL                :: ft_wspd_daymax_ms(l_mlent)
      REAL                :: ft_snowmelted_cm(l_mlent)
      REAL                :: ft_we_snowmelted_mm(l_mlent)
      REAL                :: ft_we_snowcover_mm(l_mlent)
c****
c     Output variables

      REAL                :: ft_sd_pres_min_hpa
      REAL                :: ft_sd_pres_max_hpa
      REAL                :: ft_sd_pres_avg_hpa
      REAL                :: ft_sd_pres_ngood
      REAL                :: ft_sd_pres_nbad

      REAL                :: ft_sd_airt_min_c
      REAL                :: ft_sd_airt_max_c
      REAL                :: ft_sd_airt_avg_c
      REAL                :: ft_sd_airt_ngood
      REAL                :: ft_sd_airt_nbad

      REAL                :: ft_sd_wetb_min_c
      REAL                :: ft_sd_wetb_max_c
      REAL                :: ft_sd_wetb_avg_c
      REAL                :: ft_sd_wetb_ngood
      REAL                :: ft_sd_wetb_nbad

      REAL                :: ft_sd_vprs_min_hpa
      REAL                :: ft_sd_vprs_max_hpa
      REAL                :: ft_sd_vprs_avg_hpa
      REAL                :: ft_sd_vprs_ngood
      REAL                :: ft_sd_vprs_nbad

      REAL                :: ft_sd_relh_min_pc
      REAL                :: ft_sd_relh_max_pc
      REAL                :: ft_sd_relh_avg_pc
      REAL                :: ft_sd_relh_ngood
      REAL                :: ft_sd_relh_nbad

      REAL                :: ft_sd_relhreg_min_pc
      REAL                :: ft_sd_relhreg_max_pc
      REAL                :: ft_sd_relhreg_avg_pc
      REAL                :: ft_sd_relhreg_ngood
      REAL                :: ft_sd_relhreg_nbad

      REAL                :: ft_sd_wdir_min_code
      REAL                :: ft_sd_wdir_max_code
      REAL                :: ft_sd_wdir_avg_code
      REAL                :: ft_sd_wdir_ngood
      REAL                :: ft_sd_wdir_nbad

      REAL                :: ft_sd_wspd_min_bft
      REAL                :: ft_sd_wspd_max_bft
      REAL                :: ft_sd_wspd_avg_bft
      REAL                :: ft_sd_wspd_ngood
      REAL                :: ft_sd_wspd_nbad

      REAL                :: ft_sd_ccov_min_okta
      REAL                :: ft_sd_ccov_max_okta
      REAL                :: ft_sd_ccov_avg_okta
      REAL                :: ft_sd_ccov_ngood
      REAL                :: ft_sd_ccov_nbad

      REAL                :: ft_sd_visi_min_code
      REAL                :: ft_sd_visi_max_code
      REAL                :: ft_sd_visi_avg_code
      REAL                :: ft_sd_visi_ngood
      REAL                :: ft_sd_visi_nbad

      REAL                :: ft_sd_grdcnd_min_code
      REAL                :: ft_sd_grdcnd_max_code
      REAL                :: ft_sd_grdcnd_avg_code
      REAL                :: ft_sd_grdcnd_ngood
      REAL                :: ft_sd_grdcnd_nbad

      REAL                :: ft_sd_ppt_min_mm
      REAL                :: ft_sd_ppt_max_mm
      REAL                :: ft_sd_ppt_avg_mm
      REAL                :: ft_sd_ppt_ngood
      REAL                :: ft_sd_ppt_nbad
c****
      REAL                :: ft_day_pres_min_hpa
      REAL                :: ft_day_pres_max_hpa
      REAL                :: ft_day_pres_avg_hpa
      REAL                :: ft_day_pres_ngood
      REAL                :: ft_day_pres_nbad

      REAL                :: ft_day_airt_daymax_min_c
      REAL                :: ft_day_airt_daymax_max_c
      REAL                :: ft_day_airt_daymax_avg_c
      REAL                :: ft_day_airt_daymax_ngood
      REAL                :: ft_day_airt_daymax_nbad

      REAL                :: ft_day_airt_daymin_min_c
      REAL                :: ft_day_airt_daymin_max_c
      REAL                :: ft_day_airt_daymin_avg_c
      REAL                :: ft_day_airt_daymin_ngood
      REAL                :: ft_day_airt_daymin_nbad

      REAL                :: ft_day_airt_range_min_c
      REAL                :: ft_day_airt_range_max_c
      REAL                :: ft_day_airt_range_avg_c
      REAL                :: ft_day_airt_range_ngood
      REAL                :: ft_day_airt_range_nbad

      REAL                :: ft_day_boden_daymin_min_c
      REAL                :: ft_day_boden_daymin_max_c
      REAL                :: ft_day_boden_daymin_avg_c
      REAL                :: ft_day_boden_daymin_ngood
      REAL                :: ft_day_boden_daymin_nbad

      REAL                :: ft_day_airt_dayavg_min_c
      REAL                :: ft_day_airt_dayavg_max_c
      REAL                :: ft_day_airt_dayavg_avg_c
      REAL                :: ft_day_airt_dayavg_ngood
      REAL                :: ft_day_airt_dayavg_nbad

      REAL                :: ft_day_vprs_dayavg_min_hpa
      REAL                :: ft_day_vprs_dayavg_max_hpa
      REAL                :: ft_day_vprs_dayavg_avg_hpa
      REAL                :: ft_day_vprs_dayavg_ngood
      REAL                :: ft_day_vprs_dayavg_nbad

      REAL                :: ft_day_relh_dayavg_min_pc
      REAL                :: ft_day_relh_dayavg_max_pc
      REAL                :: ft_day_relh_dayavg_avg_pc
      REAL                :: ft_day_relh_dayavg_ngood
      REAL                :: ft_day_relh_dayavg_nbad

      REAL                :: ft_day_wspd_dayavg_min_bft
      REAL                :: ft_day_wspd_dayavg_max_bft
      REAL                :: ft_day_wspd_dayavg_avg_bft
      REAL                :: ft_day_wspd_dayavg_ngood
      REAL                :: ft_day_wspd_dayavg_nbad

      REAL                :: ft_day_ccov_dayavg_min_okta
      REAL                :: ft_day_ccov_dayavg_max_okta
      REAL                :: ft_day_ccov_dayavg_avg_okta
      REAL                :: ft_day_ccov_dayavg_ngood
      REAL                :: ft_day_ccov_dayavg_nbad

      REAL                :: ft_day_sundur_daytot_min_h
      REAL                :: ft_day_sundur_daytot_max_h
      REAL                :: ft_day_sundur_daytot_avg_h
      REAL                :: ft_day_sundur_daytot_ngood
      REAL                :: ft_day_sundur_daytot_nbad

      REAL                :: ft_day_ppt_daytot_min_mm
      REAL                :: ft_day_ppt_daytot_max_mm
      REAL                :: ft_day_ppt_daytot_avg_mm
      REAL                :: ft_day_ppt_daytot_ngood
      REAL                :: ft_day_ppt_daytot_nbad

      REAL                :: ft_day_snowcover_min_cm
      REAL                :: ft_day_snowcover_max_cm
      REAL                :: ft_day_snowcover_avg_cm
      REAL                :: ft_day_snowcover_ngood
      REAL                :: ft_day_snowcover_nbad

      REAL                :: ft_day_snow_daytot_min_cm
      REAL                :: ft_day_snow_daytot_max_cm
      REAL                :: ft_day_snow_daytot_avg_cm
      REAL                :: ft_day_snow_daytot_ngood
      REAL                :: ft_day_snow_daytot_nbad

      REAL                :: ft_day_wspd_daymax_min_ms
      REAL                :: ft_day_wspd_daymax_max_ms
      REAL                :: ft_day_wspd_daymax_avg_ms
      REAL                :: ft_day_wspd_daymax_ngood
      REAL                :: ft_day_wspd_daymax_nbad

      REAL                :: ft_day_snowmelted_min_cm
      REAL                :: ft_day_snowmelted_max_cm
      REAL                :: ft_day_snowmelted_avg_cm
      REAL                :: ft_day_snowmelted_ngood
      REAL                :: ft_day_snowmelted_nbad

      REAL                :: ft_day_we_snowmelted_min_mm
      REAL                :: ft_day_we_snowmelted_max_mm
      REAL                :: ft_day_we_snowmelted_avg_mm
      REAL                :: ft_day_we_snowmelted_ngood
      REAL                :: ft_day_we_snowmelted_nbad

      REAL                :: ft_day_we_snowcover_min_mm
      REAL                :: ft_day_we_snowcover_max_mm
      REAL                :: ft_day_we_snowcover_avg_mm
      REAL                :: ft_day_we_snowcover_ngood
      REAL                :: ft_day_we_snowcover_nbad

c************************************************************************
      print*,'just inside calc_basic_stats_sd1'
c      print*,'l_mlent=',l_mlent

c     Find basic stats for single file
c     PRES
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_pres_hpa,

     +   ft_sd_pres_min_hpa,ft_sd_pres_max_hpa,
     +   ft_sd_pres_avg_hpa,
     +   ft_sd_pres_ngood,ft_sd_pres_nbad)

c     AIRT
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_airt_c,

     +   ft_sd_airt_min_c,ft_sd_airt_max_c,
     +   ft_sd_airt_avg_c,
     +   ft_sd_airt_ngood,ft_sd_airt_nbad)

c     WETB
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_wetb_c,

     +   ft_sd_wetb_min_c,ft_sd_wetb_max_c,
     +   ft_sd_wetb_avg_c,
     +   ft_sd_wetb_ngood,ft_sd_wetb_nbad)

c     VPRS
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_vprs_hpa,

     +   ft_sd_vprs_min_hpa,ft_sd_vprs_max_hpa,
     +   ft_sd_vprs_avg_hpa,
     +   ft_sd_vprs_ngood,ft_sd_vprs_nbad)

c     RELH
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_relh_pc,

     +   ft_sd_relh_min_pc,ft_sd_relh_max_pc,
     +   ft_sd_relh_avg_pc,
     +   ft_sd_relh_ngood,ft_sd_relh_nbad)

c     RELHREG
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_relhreg_pc,

     +   ft_sd_relhreg_min_pc,ft_sd_relhreg_max_pc,
     +   ft_sd_relhreg_avg_pc,
     +   ft_sd_relhreg_ngood,ft_sd_relhreg_nbad)

c     WDIR
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_wdir_code,

     +   ft_sd_wdir_min_code,ft_sd_wdir_max_code,
     +   ft_sd_wdir_avg_code,
     +   ft_sd_wdir_ngood,ft_sd_wdir_nbad)

c     WSPD
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_wspd_bft,

     +   ft_sd_wspd_min_bft,ft_sd_wspd_max_bft,
     +   ft_sd_wspd_avg_bft,
     +   ft_sd_wspd_ngood,ft_sd_wspd_nbad)

c     CCOV
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_ccov_okta,

     +   ft_sd_ccov_min_okta,ft_sd_ccov_max_okta,
     +   ft_sd_ccov_avg_okta,
     +   ft_sd_ccov_ngood,ft_sd_ccov_nbad)

c     VISI
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_visi_code,

     +   ft_sd_visi_min_code,ft_sd_visi_max_code,
     +   ft_sd_visi_avg_code,
     +   ft_sd_visi_ngood,ft_sd_visi_nbad)

c     GRDCND
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_grdcnd_code,

     +   ft_sd_grdcnd_min_code,ft_sd_grdcnd_max_code,
     +   ft_sd_grdcnd_avg_code,
     +   ft_sd_grdcnd_ngood,ft_sd_grdcnd_nbad)

c     PPT
      CALL dwd_subday_stat_pack(l_mlent,j_sd,f_ndflag,
     +   ft_sd_ppt_mm,

     +   ft_sd_ppt_min_mm,ft_sd_ppt_max_mm,
     +   ft_sd_ppt_avg_mm,
     +   ft_sd_ppt_ngood,ft_sd_ppt_nbad)
c************************************************************************
c************************************************************************
c     Find basic stats for single file
c     PRES_DAYAVG
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_pres_dayavg_hpa,

     +   ft_day_pres_min_hpa,ft_day_pres_max_hpa,
     +   ft_day_pres_avg_hpa,
     +   ft_day_pres_ngood,ft_day_pres_nbad)

c     AIRT_DAYMAX
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_airt_daymax_c,

     +   ft_day_airt_daymax_min_c,ft_day_airt_daymax_max_c,
     +   ft_day_airt_daymax_avg_c,
     +   ft_day_airt_daymax_ngood,ft_day_airt_daymax_nbad)

c     AIRT_DAYMIN
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_airt_daymin_c,

     +   ft_day_airt_daymin_min_c,ft_day_airt_daymin_max_c,
     +   ft_day_airt_daymin_avg_c,
     +   ft_day_airt_daymin_ngood,ft_day_airt_daymin_nbad)

c     AIRT_RANGE
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_airt_range_c,

     +   ft_day_airt_range_min_c,ft_day_airt_daymin_max_c,
     +   ft_day_airt_range_avg_c,
     +   ft_day_airt_range_ngood,ft_day_airt_range_nbad)

c     BODEN_DAYMIN
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_boden_daymin_c,

     +   ft_day_boden_daymin_min_c,ft_day_boden_daymin_max_c,
     +   ft_day_boden_daymin_avg_c,
     +   ft_day_boden_daymin_ngood,ft_day_boden_daymin_nbad)

c     AIRT_DAYAVG
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_airt_dayavg_c,

     +   ft_day_airt_dayavg_min_c,ft_day_airt_dayavg_max_c,
     +   ft_day_airt_dayavg_avg_c,
     +   ft_day_airt_dayavg_ngood,ft_day_airt_dayavg_nbad)

c     VPRS_DAYAVG
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_vprs_dayavg_hpa,

     +   ft_day_vprs_dayavg_min_hpa,ft_day_vprs_dayavg_max_hpa,
     +   ft_day_vprs_dayavg_avg_hpa,
     +   ft_day_vprs_dayavg_ngood,ft_day_vprs_dayavg_nbad)

c     RELH_DAYAVG
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_relh_dayavg_pc,

     +   ft_day_relh_dayavg_min_pc,ft_day_relh_dayavg_max_pc,
     +   ft_day_relh_dayavg_avg_pc,
     +   ft_day_relh_dayavg_ngood,ft_day_relh_dayavg_nbad)

c     WSPD_DAYAVG
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_wspd_dayavg_bft,

     +   ft_day_wspd_dayavg_min_bft,ft_day_wspd_dayavg_max_bft,
     +   ft_day_wspd_dayavg_avg_bft,
     +   ft_day_wspd_dayavg_ngood,ft_day_wspd_dayavg_nbad)

c     CCOV_DAYAVG
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_ccov_dayavg_okta,

     +   ft_day_ccov_dayavg_min_okta,ft_day_ccov_dayavg_max_okta,
     +   ft_day_ccov_dayavg_avg_okta,
     +   ft_day_ccov_dayavg_ngood,ft_day_ccov_dayavg_nbad)

c     SUNDUR_DAYTOT
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_sundur_daytot_h,

     +   ft_day_sundur_daytot_min_h,ft_day_sundur_daytot_max_h,
     +   ft_day_sundur_daytot_avg_h,
     +   ft_day_sundur_daytot_ngood,ft_day_sundur_daytot_nbad)

c     PPT_DAYTOT
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_ppt_daytot_mm,

     +   ft_day_ppt_daytot_min_mm,ft_day_ppt_daytot_max_mm,
     +   ft_day_ppt_daytot_avg_mm,
     +   ft_day_ppt_daytot_ngood,ft_day_ppt_daytot_nbad)

c     SNOWCOVER
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_snowcover_cm,

     +   ft_day_snowcover_min_cm,ft_day_snowcover_max_cm,
     +   ft_day_snowcover_avg_cm,
     +   ft_day_snowcover_ngood,ft_day_snowcover_nbad)

c     SNOW_DAYTOT
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_snow_daytot_cm,

     +   ft_day_snow_daytot_min_cm,ft_day_snow_daytot_max_cm,
     +   ft_day_snow_daytot_avg_cm,
     +   ft_day_snow_daytot_ngood,ft_day_snow_daytot_nbad)

c     WSPD_DAYMAX
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_wspd_daymax_ms,

     +   ft_day_wspd_daymax_min_ms,ft_day_wspd_daymax_max_ms,
     +   ft_day_wspd_daymax_avg_ms,
     +   ft_day_wspd_daymax_ngood,ft_day_wspd_daymax_nbad)

c     SNOWMELTED
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_snowmelted_cm,

     +   ft_day_snowmelted_min_cm,ft_day_snowmelted_max_cm,
     +   ft_day_snowmelted_avg_cm,
     +   ft_day_snowmelted_ngood,ft_day_snowmelted_nbad)

c     WE_SNOWMELTED
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_we_snowmelted_mm,

     +   ft_day_we_snowmelted_min_mm,ft_day_we_snowmelted_max_mm,
     +   ft_day_we_snowmelted_avg_mm,
     +   ft_day_we_snowmelted_ngood,ft_day_we_snowmelted_nbad)

c     WE_SNOWCOVER
      CALL dwd_subday_stat_pack(l_mlent,j_day,f_ndflag,
     +   ft_we_snowcover_mm,

     +   ft_day_we_snowcover_min_mm,ft_day_we_snowcover_max_mm,
     +   ft_day_we_snowcover_avg_mm,
     +   ft_day_we_snowcover_ngood,ft_day_we_snowcover_nbad)

c************************************************************************
      GOTO 18

c     Output stats
      print*,'ft_sd_pres_ngood=', ft_sd_pres_ngood
      print*,'ft_sd_pres_nbad=',  ft_sd_pres_nbad
      print*,'ft_sd_pres_min_hpa',ft_sd_pres_min_hpa
      print*,'ft_sd_pres_max_hpa',ft_sd_pres_max_hpa
      print*,'ft_sd_pres_avg_hpa',ft_sd_pres_avg_hpa

      print*,'ft_sd_airt_ngood=', ft_sd_airt_ngood
      print*,'ft_sd_airt_nbad=',  ft_sd_airt_nbad
      print*,'ft_sd_airt_min_c',  ft_sd_airt_min_c
      print*,'ft_sd_airt_max_c',  ft_sd_airt_max_c
      print*,'ft_sd_airt_avg_hpa',ft_sd_airt_avg_c

      print*,'ft_sd_wetb_ngood=', ft_sd_wetb_ngood
      print*,'ft_sd_wetb_nbad=',  ft_sd_wetb_nbad
      print*,'ft_sd_wetb_min_c',  ft_sd_wetb_min_c
      print*,'ft_sd_wetb_max_c',  ft_sd_wetb_max_c
      print*,'ft_sd_wetb_avg_c',  ft_sd_wetb_avg_c

      print*,'ft_sd_vprs_ngood=',   ft_sd_vprs_ngood
      print*,'ft_sd_vprs_nbad=',    ft_sd_vprs_nbad
      print*,'ft_sd_vprs_min_hpa',  ft_sd_vprs_min_hpa
      print*,'ft_sd_vprs_max_hpa',  ft_sd_vprs_max_hpa
      print*,'ft_sd_vprs_avg_hpa',  ft_sd_vprs_avg_hpa

      print*,'ft_sd_relh_ngood=',   ft_sd_relh_ngood
      print*,'ft_sd_relh_nbad=',    ft_sd_relh_nbad
      print*,'ft_sd_relh_min_pc',   ft_sd_relh_min_pc
      print*,'ft_sd_relh_max_pc',   ft_sd_relh_max_pc
      print*,'ft_sd_relh_avg_pc',   ft_sd_relh_avg_pc

      print*,'ft_sd_relhreg_ngood=',ft_sd_relhreg_ngood
      print*,'ft_sd_relhreg_nbad=', ft_sd_relhreg_nbad
      print*,'ft_sd_relhreg_min_pc',ft_sd_relhreg_min_pc
      print*,'ft_sd_relhreg_max_pc',ft_sd_relhreg_max_pc
      print*,'ft_sd_relhreg_avg_pc',ft_sd_relhreg_avg_pc

      print*,'ft_sd_wdir_ngood=',   ft_sd_wdir_ngood
      print*,'ft_sd_wdir_nbad=',    ft_sd_wdir_nbad
      print*,'ft_sd_wdir_min_code', ft_sd_wdir_min_code
      print*,'ft_sd_wdir_max_code', ft_sd_wdir_max_code
      print*,'ft_sd_wdir_avg_code', ft_sd_wdir_avg_code

      print*,'ft_sd_wspd_ngood=',   ft_sd_wspd_ngood
      print*,'ft_sd_wspd_nbad=',    ft_sd_wspd_nbad
      print*,'ft_sd_wspd_min_bft',  ft_sd_wspd_min_bft
      print*,'ft_sd_wspd_max_bft',  ft_sd_wspd_max_bft
      print*,'ft_sd_wspd_avg_bft',  ft_sd_wspd_avg_bft

      print*,'ft_sd_ccov_ngood=',   ft_sd_ccov_ngood
      print*,'ft_sd_ccov_nbad=',    ft_sd_ccov_nbad
      print*,'ft_sd_ccov_min_okta', ft_sd_ccov_min_okta
      print*,'ft_sd_ccov_max_okta', ft_sd_ccov_max_okta
      print*,'ft_sd_ccov_avg_okta', ft_sd_ccov_avg_okta

      print*,'ft_sd_visi_ngood=',   ft_sd_visi_ngood
      print*,'ft_sd_visi_nbad=',    ft_sd_visi_nbad
      print*,'ft_sd_visi_min_code', ft_sd_visi_min_code
      print*,'ft_sd_visi_max_code', ft_sd_visi_max_code
      print*,'ft_sd_visi_avg_code', ft_sd_visi_avg_code

      print*,'ft_sd_ppt_ngood=',    ft_sd_ppt_ngood
      print*,'ft_sd_ppt_nbad=',     ft_sd_ppt_nbad
      print*,'ft_sd_ppt_min_mm',    ft_sd_ppt_min_mm
      print*,'ft_sd_ppt_max_mm',    ft_sd_ppt_max_mm
      print*,'ft_sd_ppt_avg_mm',    ft_sd_ppt_avg_mm
c*****
      print*,'ft_day_pres_min_hpa=',     ft_day_pres_min_hpa
      print*,'ft_day_pres_max_hpa=',     ft_day_pres_max_hpa
      print*,'ft_day_pres_avg_hpa=',     ft_day_pres_avg_hpa
      print*,'ft_day_pres_ngood=',       ft_day_pres_ngood
      print*,'ft_day_pres_nbad=',        ft_day_pres_nbad

      print*,'ft_day_airt_daymax_min_c=',ft_day_airt_daymax_min_c
      print*,'ft_day_airt_daymax_max_c=',ft_day_airt_daymax_max_c
      print*,'ft_day_airt_daymax_avg_c=',ft_day_airt_daymax_avg_c
      print*,'ft_day_airt_daymax_ngood=',ft_day_airt_daymax_ngood
      print*,'ft_day_airt_daymax_nbad=', ft_day_airt_daymax_nbad

      print*,'ft_day_airt_daymin_min_c=',ft_day_airt_daymin_min_c
      print*,'ft_day_airt_daymin_max_c=',ft_day_airt_daymin_max_c
      print*,'ft_day_airt_daymin_avg_c=',ft_day_airt_daymin_avg_c
      print*,'ft_day_airt_daymin_ngood=',ft_day_airt_daymin_ngood
      print*,'ft_day_airt_daymin_nbad=', ft_day_airt_daymin_nbad

      print*,'ft_day_airt_range_min_c=', ft_day_airt_range_min_c
      print*,'ft_day_airt_range_max_c=', ft_day_airt_range_max_c
      print*,'ft_day_airt_range_avg_c=', ft_day_airt_range_avg_c
      print*,'ft_day_airt_range_ngood=', ft_day_airt_range_ngood
      print*,'ft_day_airt_range_nbad=',  ft_day_airt_range_nbad

      print*,'ft_day_boden_daymin_min_c=', ft_day_boden_daymin_min_c
      print*,'ft_day_boden_daymin_max_c=', ft_day_boden_daymin_max_c
      print*,'ft_day_boden_daymin_avg_c=', ft_day_boden_daymin_avg_c
      print*,'ft_day_boden_daymin_ngood=', ft_day_boden_daymin_ngood
      print*,'ft_day_boden_daymin_nbad=',  ft_day_boden_daymin_nbad

      print*,'ft_day_airt_dayavg_min_c=',  ft_day_airt_dayavg_min_c
      print*,'ft_day_airt_dayavg_max_c=',  ft_day_airt_dayavg_max_c
      print*,'ft_day_airt_dayavg_avg_c=',  ft_day_airt_dayavg_avg_c
      print*,'ft_day_airt_dayavg_ngood=',  ft_day_airt_dayavg_ngood
      print*,'ft_day_airt_dayavg_nbad=',   ft_day_airt_dayavg_nbad

      print*,'ft_day_vprs_dayavg_min_hpa=',ft_day_vprs_dayavg_min_hpa
      print*,'ft_day_vprs_dayavg_max_hpa=',ft_day_vprs_dayavg_max_hpa
      print*,'ft_day_vprs_dayavg_avg_hpa=',ft_day_vprs_dayavg_avg_hpa
      print*,'ft_day_vprs_dayavg_ngood=',  ft_day_vprs_dayavg_ngood
      print*,'ft_day_vprs_dayavg_nbad=',   ft_day_vprs_dayavg_nbad

      print*,'ft_day_relh_dayavg_min_pc=', ft_day_relh_dayavg_min_pc
      print*,'ft_day_relh_dayavg_max_pc=', ft_day_relh_dayavg_max_pc
      print*,'ft_day_relh_dayavg_avg_pc=', ft_day_relh_dayavg_avg_pc
      print*,'ft_day_relh_dayavg_ngood=',  ft_day_relh_dayavg_ngood
      print*,'ft_day_relh_dayavg_nbad=',   ft_day_relh_dayavg_nbad

      print*,'ft_day_wspd_dayavg_min_bft=',ft_day_wspd_dayavg_min_bft
      print*,'ft_day_wspd_dayavg_max_bft=',ft_day_wspd_dayavg_max_bft
      print*,'ft_day_wspd_dayavg_avg_bft=',ft_day_wspd_dayavg_avg_bft
      print*,'ft_day_wspd_dayavg_ngood=',  ft_day_wspd_dayavg_ngood
      print*,'ft_day_wspd_dayavg_nbad=',   ft_day_wspd_dayavg_nbad

      print*,'ft_day_ccov_dayavg_min_okta=',ft_day_ccov_dayavg_min_okta
      print*,'ft_day_ccov_dayavg_max_okta=',ft_day_ccov_dayavg_max_okta
      print*,'ft_day_ccov_dayavg_avg_okta=',ft_day_ccov_dayavg_avg_okta
      print*,'ft_day_ccov_dayavg_ngood=',   ft_day_ccov_dayavg_ngood
      print*,'ft_day_ccov_dayavg_nbad=',    ft_day_ccov_dayavg_nbad

      print*,'ft_day_sundur_daytot_min_h=', ft_day_sundur_daytot_min_h
      print*,'ft_day_sundur_daytot_max_h=', ft_day_sundur_daytot_max_h
      print*,'ft_day_sundur_daytot_avg_h=', ft_day_sundur_daytot_avg_h
      print*,'ft_day_sundur_daytot_ngood=', ft_day_sundur_daytot_ngood
      print*,'ft_day_sundur_daytot_nbad=',  ft_day_sundur_daytot_nbad

      print*,'ft_day_ppt_daytot_min_mm=',   ft_day_ppt_daytot_min_mm
      print*,'ft_day_ppt_daytot_max_mm=',   ft_day_ppt_daytot_max_mm
      print*,'ft_day_ppt_daytot_avg_mm=',   ft_day_ppt_daytot_avg_mm
      print*,'ft_day_ppt_daytot_ngood=',    ft_day_ppt_daytot_ngood
      print*,'ft_day_sundur_daytot_nbad=',  ft_day_sundur_daytot_nbad

      print*,'ft_day_snowcover_min_cm=',    ft_day_snowcover_min_cm
      print*,'ft_day_snowcover_max_cm=',    ft_day_snowcover_max_cm
      print*,'ft_day_snowcover_avg_cm=',    ft_day_snowcover_max_cm
      print*,'ft_day_snowcover_ngood=',     ft_day_snowcover_ngood
      print*,'ft_day_snowcover_nbad=',      ft_day_snowcover_nbad

      print*,'ft_day_snow_daytot_min_cm=',  ft_day_snow_daytot_min_cm
      print*,'ft_day_snow_daytot_max_cm=',  ft_day_snow_daytot_max_cm
      print*,'ft_day_snow_daytot_avg_cm=',  ft_day_snow_daytot_avg_cm
      print*,'ft_day_snow_daytot_ngood=',   ft_day_snow_daytot_ngood
      print*,'ft_day_snow_daytot_nbad=',    ft_day_snow_daytot_nbad

      print*,'ft_day_wspd_daymax_min_ms=',  ft_day_wspd_daymax_min_ms
      print*,'ft_day_wspd_daymax_max_ms=',  ft_day_wspd_daymax_max_ms
      print*,'ft_day_wspd_daymax_avg_ms=',  ft_day_wspd_daymax_avg_ms
      print*,'ft_day_wspd_daymax_ngood=',   ft_day_wspd_daymax_ngood
      print*,'ft_day_wspd_daymax_nbad=',    ft_day_wspd_daymax_nbad

      print*,'ft_day_snowmelted_min_cm=',   ft_day_snowmelted_min_cm
      print*,'ft_day_snowmelted_max_cm=',   ft_day_snowmelted_min_cm
      print*,'ft_day_snowmelted_avg_cm=',   ft_day_snowmelted_avg_cm
      print*,'ft_day_snowmelted_ngood=',    ft_day_snowmelted_ngood
      print*,'ft_day_snowmelted_nbad=',     ft_day_snowmelted_nbad

      print*,'ft_day_we_snowmelted_min_mm=',ft_day_we_snowmelted_min_mm
      print*,'ft_day_we_snowmelted_max_mm=',ft_day_we_snowmelted_max_mm
      print*,'ft_day_we_snowmelted_avg_mm=',ft_day_we_snowmelted_avg_mm
      print*,'ft_day_we_snowmelted_ngood=', ft_day_we_snowmelted_ngood
      print*,'ft_day_we_snowmelted_nbad=',  ft_day_we_snowmelted_nbad

      print*,'ft_day_we_snowcover_min_mm=', ft_day_we_snowcover_min_mm
      print*,'ft_day_we_snowcover_max_mm=', ft_day_we_snowcover_max_mm
      print*,'ft_day_we_snowcover_avg_mm=', ft_day_we_snowcover_avg_mm
      print*,'ft_day_we_snowcover_ngood=',  ft_day_we_snowcover_ngood
      print*,'ft_day_we_snowcover_nbad=',   ft_day_we_snowcover_nbad

      print*,'j_day=',j_day
      print*,'j_sd=',j_sd,j_sd/3.0

18    CONTINUE
c************************************************************************

      RETURN
      END

max_stars_repo_name: ajkettle/glamod-nuim

max_issues_repo_name: ajkettle/glamod-nuim

max_forks_repo_name: ajkettle/glamod-nuim


text:
c     Copyright (C) 2001-7  Leo Breiman and Adele Cutler and Merck & Co, Inc.
c     This program is free software; you can redistribute it and/or
c     modify it under the terms of the GNU General Public License
c     as published by the Free Software Foundation; either version 2
c     of the License, or (at your option) any later version.

c     This program is distributed in the hope that it will be useful,
c     but WITHOUT ANY WARRANTY; without even the implied warranty of
c     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c     GNU General Public License for more details.
c
c     Modified by Andy Liaw and Matt Wiener:
c     The main program is re-written as a C function to be called from R.
c     All calls to the uniform RNG is replaced with R's RNG.  Some subroutines
c     not called are excluded.  Variables and arrays declared as double as
c     needed.  Unused variables are deleted.
c
c     SUBROUTINE BUILDTREE

      subroutine buildtree(a, b, cl, cat, maxcat, mdim, nsample,
     1     nclass, treemap, bestvar, bestsplit, bestsplitnext, tgini,
     1     nodestatus,nodepop, nodestart, classpop, tclasspop,
     1     tclasscat,ta,nrnodes, idmove, ndsize, ncase, mtry, iv,
     1     nodeclass, ndbigtree, win, wr, wl, mred, nuse, mind)

c     Buildtree consists of repeated calls to two subroutines, Findbestsplit
c     and Movedata.  Findbestsplit does just that--it finds the best split of
c     the current node.  Movedata moves the data in the split node right and
c     left so that the data corresponding to each child node is contiguous.
c     The buildtree bookkeeping is different from that in Friedman's original
c     CART program.  ncur is the total number of nodes to date.
c     nodestatus(k)=1 if the kth node has been split.  nodestatus(k)=2 if the
c     node exists but has not yet been split, and =-1 of the node is terminal.
c     A node is terminal if its size is below a threshold value, or if it is
c     all one class, or if all the x-values are equal.  If the current node k
c     is split, then its children are numbered ncur+1 (left), and
c     ncur+2(right), ncur increases to ncur+2 and the next node to be split is
c     numbered k+1.  When no more nodes can be split, buildtree returns to the
c     main program.

      implicit double precision(a-h,o-z)
      integer a(mdim, nsample), cl(nsample), cat(mdim),
     1     treemap(2,nrnodes), bestvar(nrnodes),
     1     bestsplit(nrnodes), nodestatus(nrnodes), ta(nsample),
     1     nodepop(nrnodes), nodestart(nrnodes),
     1     bestsplitnext(nrnodes), idmove(nsample),
     1     ncase(nsample), b(mdim,nsample),
     1     iv(mred), nodeclass(nrnodes), mind(mred)

      double precision tclasspop(nclass), classpop(nclass, nrnodes),
     1     tclasscat(nclass, 32), win(nsample), wr(nclass),
     1     wl(nclass), tgini(mdim), xrand
      integer msplit, ntie

      msplit = 0
      call zerv(nodestatus,nrnodes)
      call zerv(nodestart,nrnodes)
      call zerv(nodepop,nrnodes)
      call zermr(classpop,nclass,nrnodes)

      do j=1,nclass
         classpop(j, 1) = tclasspop(j)
      end do
      ncur = 1
      nodestart(1) = 1
      nodepop(1) = nuse
      nodestatus(1) = 2
c     start main loop
      do 30 kbuild = 1, nrnodes
         if (kbuild .gt. ncur) goto 50
         if (nodestatus(kbuild) .ne. 2) goto 30
c     initialize for next call to findbestsplit
         ndstart = nodestart(kbuild)
         ndend = ndstart + nodepop(kbuild) - 1
         do j = 1, nclass
            tclasspop(j) = classpop(j,kbuild)
         end do
         jstat = 0

         call findbestsplit(a,b,cl,mdim,nsample,nclass,cat,maxcat,
     1        ndstart, ndend,tclasspop,tclasscat,msplit, decsplit,
     1        nbest,ncase, jstat,mtry,win,wr,wl,mred,mind)
         if (jstat .eq. -1) then
            nodestatus(kbuild) = -1
            goto 30
         else
            bestvar(kbuild) = msplit
            iv(msplit) = 1
            if (decsplit .lt. 0.0) decsplit = 0.0
            tgini(msplit) = tgini(msplit) + decsplit
            if (cat(msplit) .eq. 1) then
               bestsplit(kbuild) = a(msplit,nbest)
               bestsplitnext(kbuild) = a(msplit,nbest+1)
            else
               bestsplit(kbuild) = nbest
               bestsplitnext(kbuild) = 0
            endif
         endif

         call movedata(a,ta,mdim,nsample,ndstart,ndend,idmove,ncase,
     1        msplit,cat,nbest,ndendl)

c     leftnode no.= ncur+1, rightnode no. = ncur+2.
         nodepop(ncur+1) = ndendl - ndstart + 1
         nodepop(ncur+2) = ndend - ndendl
         nodestart(ncur+1) = ndstart
         nodestart(ncur+2) = ndendl + 1

c     find class populations in both nodes
         do n = ndstart, ndendl
            nc = ncase(n)
            j=cl(nc)
            classpop(j,ncur+1) = classpop(j,ncur+1) + win(nc)
         end do
         do n = ndendl+1, ndend
            nc = ncase(n)
            j = cl(nc)
            classpop(j,ncur+2) = classpop(j,ncur+2) + win(nc)
         end do
c     check on nodestatus
         nodestatus(ncur+1) = 2
         nodestatus(ncur+2) = 2
         if (nodepop(ncur+1).le.ndsize) nodestatus(ncur+1) = -1
         if (nodepop(ncur+2).le.ndsize) nodestatus(ncur+2) = -1
         popt1 = 0
         popt2 = 0
         do j = 1, nclass
            popt1 = popt1 + classpop(j,ncur+1)
            popt2 = popt2 + classpop(j,ncur+2)
         end do

         do j=1,nclass
            if (classpop(j,ncur+1).eq.popt1) nodestatus(ncur+1) = -1
            if (classpop(j,ncur+2).eq.popt2) nodestatus(ncur+2) = -1
         end do

         treemap(1,kbuild) = ncur + 1
         treemap(2,kbuild) = ncur + 2
         nodestatus(kbuild) = 1
         ncur = ncur+2
         if (ncur.ge.nrnodes) goto 50

 30   continue
 50   continue

      ndbigtree = nrnodes
      do k=nrnodes, 1, -1
         if (nodestatus(k) .eq. 0) ndbigtree = ndbigtree - 1
         if (nodestatus(k) .eq. 2) nodestatus(k) = -1
      end do

c     form prediction in terminal nodes
      do kn = 1, ndbigtree
         if(nodestatus(kn) .eq. -1) then
            pp = 0
            ntie = 1
            do j = 1, nclass
               if (classpop(j,kn) .gt. pp) then
                  nodeclass(kn) = j
                  pp = classpop(j,kn)
               end if
c     Break ties at random:
               if (classpop(j,kn) .eq. pp) then
                  ntie = ntie + 1
                  call rrand(xrand)
                  if (xrand .lt. 1.0 / ntie) then
                     nodeclass(kn)=j
                     pp=classpop(j,kn)
                  end if
               end if
            end do
         end if
      end do

      end

c     SUBROUTINE FINDBESTSPLIT
c     For the best split, msplit is the variable split on. decsplit is the
c     dec. in impurity.  If msplit is numerical, nsplit is the case number
c     of value of msplit split on, and nsplitnext is the case number of the
c     next larger value of msplit.  If msplit is categorical, then nsplit is
c     the coding into an integer of the categories going left.
      subroutine findbestsplit(a, b, cl, mdim, nsample, nclass, cat,
     1     maxcat, ndstart, ndend, tclasspop, tclasscat, msplit,
     2     decsplit, nbest, ncase, jstat, mtry, win, wr, wl,
     3     mred, mind)
      implicit double precision(a-h,o-z)
      integer a(mdim,nsample), cl(nsample), cat(mdim),
     1     ncase(nsample), b(mdim,nsample), nn, j
      double precision tclasspop(nclass), tclasscat(nclass,32), dn(32),
     1     win(nsample), wr(nclass), wl(nclass), xrand
      integer mind(mred), ncmax, ncsplit,nhit, ntie
      ncmax = 10
      ncsplit = 512
c     compute initial values of numerator and denominator of Gini
      pno = 0.0
      pdo = 0.0
      do j = 1, nclass
         pno = pno + tclasspop(j) * tclasspop(j)
         pdo = pdo + tclasspop(j)
      end do
      crit0 = pno / pdo
      jstat = 0

c     start main loop through variables to find best split
      critmax = -1.0e25
      do k = 1, mred
         mind(k) = k
      end do
      nn = mred
c     sampling mtry variables w/o replacement.
      do mt = 1, mtry
         call rrand(xrand)
         j = int(nn * xrand) + 1
         mvar = mind(j)
         mind(j) = mind(nn)
         mind(nn) = mvar
         nn = nn - 1
         lcat = cat(mvar)
         if (lcat .eq. 1) then
c     Split on a numerical predictor.
            rrn = pno
            rrd = pdo
            rln = 0
            rld = 0
            call zervr(wl, nclass)
            do j = 1, nclass
               wr(j) = tclasspop(j)
            end do
            ntie = 1
            do nsp = ndstart, ndend-1
               nc = a(mvar, nsp)
               u = win(nc)
               k = cl(nc)
               rln = rln + u * (2 * wl(k) + u)
               rrn = rrn + u * (-2 * wr(k) + u)
               rld = rld + u
               rrd = rrd - u
               wl(k) = wl(k) + u
               wr(k) = wr(k) - u
               if (b(mvar, nc) .lt. b(mvar, a(mvar, nsp + 1))) then
c     If neither nodes is empty, check the split.
                  if (dmin1(rrd, rld) .gt. 1.0e-5) then
                     crit = (rln / rld) + (rrn / rrd)
                     if (crit .gt. critmax) then
                        nbest = nsp
                        critmax = crit
                        msplit = mvar
                     end if
c     Break ties at random:
                     if (crit .eq. critmax) then
                        ntie = ntie + 1
                        call rrand(xrand)
                        if (xrand .lt. 1.0 / ntie) then
                           nbest = nsp
                           critmax = crit
                           msplit = mvar
                        end if
                     end if
                  end if
               end if
            end do
         else
c     Split on a categorical predictor.  Compute the decrease in impurity.
            call zermr(tclasscat, nclass, 32)
            do nsp = ndstart, ndend
               nc = ncase(nsp)
               l = a(mvar, ncase(nsp))
               tclasscat(cl(nc), l) = tclasscat(cl(nc), l) + win(nc)
            end do
            nnz = 0
            do i = 1, lcat
               su = 0
               do j = 1, nclass
                  su = su + tclasscat(j, i)
               end do
               dn(i) = su
               if(su .gt. 0) nnz = nnz + 1
            end do
            nhit = 0
            if (nnz .gt. 1) then
               if (nclass .eq. 2 .and. lcat .gt. ncmax) then
                  call catmaxb(pdo, tclasscat, tclasspop, nclass,
     &                 lcat, nbest, critmax, nhit, dn)
               else
                  call catmax(pdo, tclasscat, tclasspop, nclass, lcat,
     &                 nbest, critmax, nhit, maxcat, ncmax, ncsplit)
               end if
               if (nhit .eq. 1) msplit = mvar
c            else
c               critmax = -1.0e25
            end if
         end if
      end do
      if (critmax .lt. -1.0e10 .or. msplit .eq. 0) jstat = -1
      decsplit = critmax - crit0
      return
      end

C     ==============================================================
c     SUBROUTINE MOVEDATA
c     This subroutine is the heart of the buildtree construction. Based on the
c     best split the data in the part of the a matrix corresponding to the
c     current node is moved to the left if it belongs to the left child and
c     right if it belongs to the right child.

      subroutine movedata(a,ta,mdim,nsample,ndstart,ndend,idmove,
     1     ncase,msplit,cat,nbest,ndendl)
      implicit double precision(a-h,o-z)
      integer a(mdim,nsample),ta(nsample),idmove(nsample),
     1     ncase(nsample),cat(mdim),icat(32)

c     compute idmove=indicator of case nos. going left

      if (cat(msplit).eq.1) then
         do nsp=ndstart,nbest
            nc=a(msplit,nsp)
            idmove(nc)=1
         end do
         do nsp=nbest+1, ndend
            nc=a(msplit,nsp)
            idmove(nc)=0
         end do
         ndendl=nbest
      else
         ndendl=ndstart-1
         l=cat(msplit)
         call myunpack(l,nbest,icat)
         do nsp=ndstart,ndend
            nc=ncase(nsp)
            if (icat(a(msplit,nc)).eq.1) then
               idmove(nc)=1
               ndendl=ndendl+1
            else
               idmove(nc)=0
            endif
         end do
      endif

c     shift case. nos. right and left for numerical variables.

      do 40 msh=1,mdim
         if (cat(msh).eq.1) then
            k=ndstart-1
            do 50 n=ndstart,ndend
               ih=a(msh,n)
               if (idmove(ih).eq.1) then
                  k=k+1
                  ta(k)=a(msh,n)
               endif
 50         continue
            do 60 n=ndstart,ndend
               ih=a(msh,n)
               if (idmove(ih).eq.0) then
                  k=k+1
                  ta(k)=a(msh,n)
               endif
 60         continue

            do 70 k=ndstart,ndend
               a(msh,k)=ta(k)
 70         continue
         endif

 40   continue
      ndo=0
      if (ndo.eq.1) then
         do 140 msh = 1, mdim
            if (cat(msh) .gt. 1) then
               k = ndstart - 1
               do 150 n = ndstart, ndend
                  ih = ncase(n)
                  if (idmove(ih) .eq. 1) then
                     k = k + 1
                     ta(k) = a(msh, ih)
                  endif
 150           continue
               do 160 n = ndstart, ndend
                  ih = ncase(n)
                  if (idmove(ih) .eq. 0) then
                     k = k + 1
                     ta(k) = a(msh,ih)
                  endif
 160           continue

               do 170 k = ndstart, ndend
                  a(msh,k) = ta(k)
 170           continue
            endif

 140     continue
      end if

c     compute case nos. for right and left nodes.

      if (cat(msplit).eq.1) then
         do 80 n=ndstart,ndend
            ncase(n)=a(msplit,n)
 80      continue
      else
         k=ndstart-1
         do 90 n=ndstart, ndend
            if (idmove(ncase(n)).eq.1) then
               k=k+1
               ta(k)=ncase(n)
            endif
 90      continue
         do 100 n=ndstart,ndend
            if (idmove(ncase(n)).eq.0) then
               k=k+1
               ta(k)=ncase(n)
            endif
 100     continue
         do 110 k=ndstart,ndend
            ncase(k)=ta(k)
 110     continue
      endif

      end

      subroutine myunpack(l,npack,icat)

c     npack is a long integer.  The sub. returns icat, an integer of zeroes and
c     ones corresponding to the coefficients in the binary expansion of npack.

      integer icat(32),npack
      do j=1,32
         icat(j)=0
      end do
      n=npack
      icat(1)=mod(n,2)
      do k=2,l
         n=(n-icat(k-1))/2
         icat(k)=mod(n,2)
      end do
      end

      subroutine zerv(ix,m1)
      integer ix(m1)
      do n=1,m1
         ix(n)=0
      end do
      end

      subroutine zervr(rx,m1)
      double precision rx(m1)
      do n=1,m1
         rx(n)=0.0d0
      end do
      end

      subroutine zerm(mx,m1,m2)
      integer mx(m1,m2)
      do i=1,m1
         do j=1,m2
            mx(i,j)=0
         end do
      end do
      end

      subroutine zermr(rx,m1,m2)
      double precision rx(m1,m2)
      do i=1,m1
         do j=1,m2
            rx(i,j)=0.0d0
         end do
      end do
      end

      subroutine zermd(rx,m1,m2)
      double precision rx(m1,m2)
      do i=1,m1
         do j=1,m2
            rx(i,j)=0.0d0
         end do
      end do
      end

ext: f

max_issues_repo_name: andudu/rhoana

max_forks_repo_name: andudu/rhoana


text:
 PROGRAM SIMPLE

  !set up modules
  USE sps_vars; USE sps_utils
  IMPLICIT NONE

  !NB: the various structure types are defined in sps_vars.f90
  !    variables not explicitly defined here are defined in sps_vars.f90

  INTEGER :: i
  !define variable for SSP spectrum
  REAL(SP), DIMENSION(ntfull,nspec)  :: spec_ssp
  !define variables for Mass and Lbol info
  REAL(SP), DIMENSION(ntfull)    :: mass_ssp,lbol_ssp
  CHARACTER(100) :: file1='', file2=''
  !structure containing all necessary parameters
  TYPE(PARAMS) :: pset
  !define structure for CSP spectrum
  TYPE(COMPSPOUT), DIMENSION(ntfull) :: ocompsp
  REAL(SP) :: ssfr6,ssfr7,ssfr8,ave_age

  !---------------------------------------------------------------!
  !---------------------------------------------------------------!

  ! Lets compute an SSP, solar metallicity, with a Chabrier IMF
  ! with no dust, and the 'default' assumptions regarding the
  ! locations of the isochrones

  imf_type  = 0             !define the IMF (1=Chabrier 2003)
                            !see sps_vars.f90 for details of this var
  pset%zmet = 10            !define the metallicity (see the manual)
                            !20 = solar metallacity

  CALL SPS_SETUP(pset%zmet) !read in the isochrones and spectral libraries

  !define the parameter set.  These are the default values, specified
  !in sps_vars.f90, but are explicitly included here for transparency
  pset%sfh   = 0     !set SFH to "SSP"
  pset%const = 1.
  pset%zred  = 0.0   !redshift
  pset%dust1 = 0.0   !dust parameter 1
  pset%dust2 = 0.0   !dust parameter 2

  pset%dell  = 0.0   !shift in log(L) for TP-AGB stars
  pset%delt  = 0.0   !shift in log(Teff) for TP-AGB stars
  pset%fbhb  = 0.0   !fraction of blue HB stars
  pset%sbss  = 0.0   !specific frequency of BS stars

  !compute the SSP
  CALL SSP_GEN(pset,mass_ssp,lbol_ssp,spec_ssp)
  !compute mags and write out mags and spec for SSP
  file1 = 'SSP_BPASS.out'
  CALL COMPSP(3,1,file1,mass_ssp,lbol_ssp,spec_ssp,pset,ocompsp)

  ! Now lets compute a 1 Gyr tau model SFH with a Salpeter IMF
  ! with a simple dust model, at a particular time
  ! (rather than outputing all the time info)

  imf_type  = 0                !define the IMF (0=Salpeter)
                               !see sps_vars.f90 for details of this var

  !NB: you only need to re-run SPS_SETUP if you have changed the metallicity
  !    or, even better (but slower), you can call SPS_SETUP(-1) and this will
  !    set up all the metallicities at once.
  CALL SPS_SETUP(pset%zmet)    !read in the isochrones and spectral libraries

  !define the parameter set.
  pset%sfh   = 1     !set SFH to "CSP"; sfh=1 means a normal tau model
  pset%tau   = 2.0   !tau units are Gyr
  pset%dust1 = 1.0   !dust parameter 1
  pset%dust2 = 0.3   !dust parameter 2
  !pset%tage  = 12.5  !age at which we want the mags

  !compute the CSP
  CALL SSP_GEN(pset,mass_ssp,lbol_ssp,spec_ssp)
  !compute mags, and write out mags and spec for CSP
  file1 = 'CSP.out'
  CALL COMPSP(3,1,file1,mass_ssp,lbol_ssp,spec_ssp,pset,ocompsp)

  !compute basic SFH statistics for the last entry in the ocompsp array
  !results are returned in the variables ssfr6,...,ave_age
  CALL SFHSTAT(pset,ocompsp(1),ssfr6,ssfr7,ssfr8,ave_age)

END PROGRAM SIMPLE

max_stars_repo_name: christopherlovell/fsps

max_issues_repo_name: christopherlovell/fsps

max_forks_repo_name: christopherlovell/fsps


text:
!
! 3-noded Lagrange element for Poisson equation in 2D
!

MODULE LAGRANGEELEMENT_2DPOISSON_TRIA3NODE
IMPLICIT NONE
PUBLIC :: LagrangeElem2DPoissonTria3Node

TYPE LagrangeElem2DPoissonTria3Node
    PRIVATE

    !
    ! member variables
    !
    INTEGER :: elenum             ! element number
    INTEGER :: ndim=2             ! dimension
    INTEGER :: degree=1           ! degree of polynomial basis functions
    INTEGER :: npElem=3           ! nodes per element
    INTEGER :: ndof=1             ! number of DOF per node
    INTEGER :: nsize=3            ! number of DOF per element
    INTEGER :: elmType            ! element type
    INTEGER :: matType            ! material type
    INTEGER :: secType            ! section type
    INTEGER :: subdomId           ! subdomain/processor id for parallel code
    INTEGER :: elmTypeNameNum     ! element type name
    INTEGER :: matId              ! material id
    INTEGER :: finiteInt          ! flag for small/finite strain
    INTEGER :: sss                ! plane stress or plane strain or axisymmetric
    INTEGER :: nivGP=0            ! number of internal variables per Gauss point
    INTEGER :: nGP=1              ! number of Gauss points

    ! internal variables, for plasicity problems
    DOUBLE PRECISION, DIMENSION(:), allocatable :: intVar1, intVar2

    INTEGER :: nodeNums(3), forAssyVec(3), globalDOFnums(3)

    !
    ! member functions
    !
    CONTAINS

        !PROCEDURE ::  Constructor
        PROCEDURE ::  prepareElemData
        PROCEDURE ::  calcStiffnessAndResidual
        PROCEDURE ::  assembleElementMatrixAndVector

END TYPE LagrangeElem2DPoissonTria3Node

CONTAINS
    SUBROUTINE Constructor(this)
      class(LagrangeElem2DPoissonTria3Node) :: this
          INTEGER :: degree = 1
          INTEGER :: npElem = 3
          INTEGER :: nlbf   = 3
          INTEGER :: ndof   = 1
          INTEGER :: nsize  = 3
    END SUBROUTINE Constructor

    SUBROUTINE prepareElemData(this)
      class(LagrangeElem2DPoissonTria3Node) :: this
          INTEGER :: degree = 1
          INTEGER :: npElem = 3
          INTEGER :: nlbf   = 3
          INTEGER :: ndof   = 1
          INTEGER :: nsize  = 3
    END SUBROUTINE prepareElemData

    SUBROUTINE calcStiffnessAndResidual(this, vertices, Klocal, Flocal)
      IMPLICIT NONE
      class(LagrangeElem2DPoissonTria3Node) :: this

      DOUBLE PRECISION, DIMENSION(3,3) :: Klocal
      DOUBLE PRECISION, DIMENSION(3) :: Flocal
      DOUBLE PRECISION, DIMENSION(:,:) :: vertices
      DOUBLE PRECISION, DIMENSION(3,2) :: Bmat
      DOUBLE PRECISION, DIMENSION(2,3) :: BmatTrans
      DOUBLE PRECISION :: x1, x2, x3, y1, y2, y3
      DOUBLE PRECISION :: area
      DOUBLE PRECISION :: kx=1.0, ky=1.0

      x1 = vertices(this%nodeNums(1),1)
      y1 = vertices(this%nodeNums(1),2)
      x2 = vertices(this%nodeNums(2),1)
      y2 = vertices(this%nodeNums(2),2)
      x3 = vertices(this%nodeNums(3),1)
      y3 = vertices(this%nodeNums(3),2)

      area = 0.5*(x1*(y2-y3) + x2*(y3-y1) + x3*(y1-y2))
      write(*,*) " area = ", area

      Bmat(1,1) = y2-y3; Bmat(2,1) = y3-y1; Bmat(3,1) = y1-y2
      Bmat(1,2) = x3-x2; Bmat(2,2) = x1-x3; Bmat(3,2) = x2-x1
      Bmat = Bmat/(2.0*area)

      Klocal = 0.0
      Flocal = 0.0

      BmatTrans = TRANSPOSE(Bmat)

      Klocal = MATMUL(Bmat, BmatTrans)
      Klocal = area*Klocal

      ! Flocal -= Klocal*dispC;

    END SUBROUTINE calcStiffnessAndResidual

    SUBROUTINE assembleElementMatrixAndVector(this, firstIter, flag, Klocal, Flocal, Mat, rhs, reac)
      IMPLICIT NONE
      class(LagrangeElem2DPoissonTria3Node) :: this
      LOGICAL :: firstIter, flag
      DOUBLE PRECISION :: fact=1.0
      DOUBLE PRECISION, dimension(:,:) :: Mat, Klocal
      DOUBLE PRECISION, dimension(:) :: Flocal, rhs, reac
      INTEGER :: i, j, r, c

      do i=1,this%npElem
        !add reaction forces
        !reac(globalDOFnums(ii)) = reac(globalDOFnums(ii)) + Flocal(ii);
        r = this%forAssyVec(i)

        if(r .ne. 0) then
          rhs(r) = rhs(r) + Flocal(i);

          do j=1,this%npElem
            c = this%forAssyVec(j)
            if(c .ne. 0) then
              Mat(r, c) = Mat(r,c) + Klocal(i,j)
            end if
          end do
        end if
      end do

      if(firstIter .eqv. .true.) then
        do i=1, this%npElem
          r = this%forAssyVec(i)

          if(r .ne. 0) then
            !fact = var1applied(globalDOFnums(ii))

            do j=1,this%npElem
              c = this%forAssyVec(j)
              if(c .ne. 0) then
                rhs(r) = rhs(r) - Klocal(j,i) * fact
              end if
            end do
          end if
        end do
      end if
    END SUBROUTINE assembleElementMatrixAndVector

END MODULE LAGRANGEELEMENT_2DPOISSON_TRIA3NODE

max_stars_repo_name: chennachaos/ParallelFEM-Fortran

max_issues_repo_name: chennachaos/ParallelFEM-Fortran

max_forks_repo_name: chennachaos/ParallelFEM-Fortran


text:
c$Id: imgvalx.f,v 1.2 2007/04/25 05:15:34 avi Exp $
c$Log: imgvalx.f,v $
c Revision 1.2  2007/04/25  05:15:34  avi
c gcc v3 compliant
c
c Revision 1.1  1998/10/02  20:23:27  mcavoy
c Initial revision
c
c Revision 1.1  1997/11/14  17:16:24  tscull
c Initial revision
c
c Revision 1.3  1997/04/15  08:21:03  avi
c trap input image NaNs in subroutine imgvalx
c
c Revision 1.2  1996/12/28  01:40:42  avi
c minor changes - don't evaluate nonexistant memory if nz.eq.1
c
c Revision 1.1  1996/12/27  20:48:20  avi
c Initial revision
c
c revision 1.2  1993/10/25  17:37:33  avi

      subroutine imgvalx(imgt,nx,ny,nz,center,mmppix,x,v,lslice)
c     avi snyder 08-06-94
c     returns in v the interpolated value of imgt at locus x.
c     lslice returns the slice from which most of the data is taken.
c     lslice = 0 (and v = 0.) if imgt is undefined at locus x.

      real*4 imgt(nx,ny,nz)
      real*4 center(3),mmppix(3),x(3)
      logical*4 ldefined,lok
      real*4 wx(3),xf(3)
      integer*4 ix(3),mx(3)

      mx(1)=nx
      mx(2)=ny
      mx(3)=nz

      ldefined=.true.
      do 1 k=1,3
      center(k)=sign(center(k),mmppix(k))
      xf(k)=(center(k)+x(k))/mmppix(k)
      ix(k)=nint(xf(k)-0.5)
      wx(k)=xf(k)-float(ix(k))
      if((ix(k).eq.0.and.wx(k).gt.0.999).or.mx(k).eq.1)then
        ix(k)=1
        wx(k)=0.
      endif
      if(mx(k).gt.1.and.ix(k).eq.mx(k).and.wx(k).lt.0.001)then
        ix(k)=mx(k)-1
        wx(k)=1.
      endif
    1 ldefined=ldefined.and.((mx(k).eq.1).or.(ix(k).gt.0).and.(ix(k).lt.mx(k)))

      if(.not.ldefined)then
        v=0.
        lslice=0
        return
      endif

      lok=(imgt(ix(1)+0,ix(2)+0,ix(3)+0).ge.0.0.or.imgt(ix(1)+0,ix(2)+0,ix(3)+0).le.0.0).and.
     &    (imgt(ix(1)+1,ix(2)+0,ix(3)+0).ge.0.0.or.imgt(ix(1)+1,ix(2)+0,ix(3)+0).le.0.0).and.
     &    (imgt(ix(1)+0,ix(2)+1,ix(3)+0).ge.0.0.or.imgt(ix(1)+0,ix(2)+1,ix(3)+0).le.0.0).and.
     &    (imgt(ix(1)+1,ix(2)+1,ix(3)+0).ge.0.0.or.imgt(ix(1)+1,ix(2)+1,ix(3)+0).le.0.0).and.
     &    (imgt(ix(1)+0,ix(2)+0,ix(3)+1).ge.0.0.or.imgt(ix(1)+0,ix(2)+0,ix(3)+1).le.0.0).and.
     &    (imgt(ix(1)+1,ix(2)+0,ix(3)+1).ge.0.0.or.imgt(ix(1)+1,ix(2)+0,ix(3)+1).le.0.0).and.
     &    (imgt(ix(1)+0,ix(2)+1,ix(3)+1).ge.0.0.or.imgt(ix(1)+0,ix(2)+1,ix(3)+1).le.0.0).and.
     &    (imgt(ix(1)+1,ix(2)+1,ix(3)+1).ge.0.0.or.imgt(ix(1)+1,ix(2)+1,ix(3)+1).le.0.0)
      if(.not.lok)then
        v=0.
        lslice=0
        return
      endif

      v=
     & +(1.-wx(3))*((1.-wx(2))*(
     &               (1.-wx(1))*imgt(ix(1)+0,ix(2)+0,ix(3)+0)
     &                  +wx(1) *imgt(ix(1)+1,ix(2)+0,ix(3)+0))
     &               +wx(2)*(
     &               (1.-wx(1))*imgt(ix(1)+0,ix(2)+1,ix(3)+0)
     &                  +wx(1) *imgt(ix(1)+1,ix(2)+1,ix(3)+0)))
       if(mx(3).gt.1)v=v
     &      +wx(3)*((1.-wx(2))*(
     &               (1.-wx(1))*imgt(ix(1)+0,ix(2)+0,ix(3)+1)
     &                  +wx(1) *imgt(ix(1)+1,ix(2)+0,ix(3)+1))
     &               +wx(2)*(
     &               (1.-wx(1))*imgt(ix(1)+0,ix(2)+1,ix(3)+1)
     &                  +wx(1) *imgt(ix(1)+1,ix(2)+1,ix(3)+1)))
      lslice=ix(3)
      if(wx(3).gt.0.5)lslice=lslice+1
      return
      end

      subroutine peak_find(mode,imgt,nx,ny,nz,mmppix,centert,coord)
c     iabs(mode) optimize
c     1          x
c     2          y
c     3          x y
c     4          z
c     5          x z
c     6          y z
c     7          x y z
c     sign of extremum found = sign of mode
      real*4 imgt(nx,ny,nz)
      real*4 mmppix(3),centert(3),coord(3)
      parameter (nterm=3)
      real*4 taram(nterm),curva(nterm)
      real*4 daram(nterm)/3*2./
      real*4 coef(0:2),array(3,3)
      parameter (ni=2)
      parameter (npts=2*ni+1)
      real*4 x(-ni:ni),y(-ni:ni)
      logical*4 liter,lfar

      call imgvalx(imgt,nx,ny,nz,centert,mmppix,coord,v,lslice)
      if(lslice.le.0)goto 99
      write(*,102)coord,v
  102 format(3f7.2,f10.1)
      sense=sign(1.,float(mode))
      jmask=iabs(mode)

      nnn=0
      niter=5
   83 liter=.false.
      lfar=.false.
      do 81 j=1,3
      if(iand(jmask,2**(j-1)).eq.0)goto 81
      taram(j)=coord(j)
      do 84 i=-ni,ni
      coord(j)=taram(j)+daram(j)*float(i)/float(ni)
      call imgvalx(imgt,nx,ny,nz,centert,mmppix,coord,v,lslice)
      if(lslice.le.0)goto 99
      x(i)=float(i)/float(ni)
   84 y(i)=-sense*v
      call polfit(x,y,npts,coef,array,3,chisqr)
      t=.5*coef(1)/coef(2)
      if(coef(2).lt.0..or.abs(t).gt..5)then
        t0=t
        t=sign(.5,coef(1))
c       write(*,"(' T ',f10.6,' ->',f10.6)")t0,t
        lfar=.true.
      endif
      liter=liter.or.abs(t).gt.0.02
      curva(j)=coef(2)
      coord(j)=taram(j)-daram(j)*t
   81 continue
      nnn=nnn+1
      if(nnn.gt.20)then
        write(*,"(' failed convergence')")
        goto 89
      endif
      if(.not.lfar)then
        niter=niter-1
      endif
      if(liter.and.niter.gt.0)goto 83

      call imgvalx(imgt,nx,ny,nz,centert,mmppix,coord,v,lslice)
      if(lslice.le.0)goto 99
      write(*,102)coord,v
   89 return
   99 write (*,"('peak_find: image boundary reached')")
      return
      end

      subroutine edge_trace(imgt,nx,ny,nz,mmppix,center,param,edge,npmax,np)
      real*4 imgt(nx,ny,nz)
      real*4 mmppix(3),center(3),edge(3,npmax)
      real*4 param(16)
c     param(01:03)	normal to osculating plane
c     param(04)		del arclength
c     param(05)		sigma eval
      real*4 n1(3)/0.,0.,1./		! for initial t, n guess
      real*4 n2(3)/0.,1.,0./		! for initial t, n guess
      real*4 t(3),n(3),b(3),x(3)
      real*4 coef(0:2),array(3,3)
      parameter (nk=3)
      parameter (npts=2*nk+1)
      real*4 u(-nk:nk),v(-nk-1:nk+1),dv(-nk:nk)
      parameter(nr=32)
      logical*4 ldebug/.false./

      pi=4.*atan(1.)
      a=sqrt(dot(param,param))		! compute binormal unit vector
      do i=1,3
        b(i)=param(i)/a
      enddo

c     initial guess for t, n
      if(abs(dot(b,n1)).lt.abs(dot(b,n2)))then
        do i=1,3
          n(i)=n1(i)
        enddo
      else
        do i=1,3
          n(i)=n2(i)
        enddo
      endif
      call cross(n,b,t)

      np=1
   11 s=0.
      c=0.
      if(ldebug)write(*, "('t ',3f10.4)")t
      if(ldebug)write(*, "('n ',3f10.4)")n
      if(ldebug)write(*, "('np, edge(np)',i6,3f7.2)")np,(edge(i,np),i=1,3)
      do 21 k=0,nr-1
      a=2.*pi*float(k)/float(nr)
      do i=1,3
        x(i)=edge(i,np)+param(05)*(cos(a)*n(i)+sin(a)*t(i))
      enddo
      call imgvalx(imgt,nx,ny,nz,center,mmppix,x,g,lslice)
      if(lslice.le.0)goto 99
      s=s+sin(a)*g
      c=c+cos(a)*g
c     if(ldebug)write(*, "(i6,3f7.2,f7.0)")k,x,g
   21 continue
      theta=atan2(s,c)
      if(ldebug)write(*, "('s,c,atan2(s,c) (deg)',3f10.4)")s,c,theta*45./atan(1.)
      do i=1,3
        n(i)=sin(theta)*t(i)+cos(theta)*n(i)
      enddo
      call cross(n,b,t)
      if(ldebug)write(*, "('t ',3f10.4)")t
      if(ldebug)write(*, "('n ',3f10.4)")n

      iter=0
   10 iter=iter+1
      if(ldebug)write(*, "('iter',i4)")iter
      do 9 k=-nk-1,nk+1			! sample image along normal
      do i=1,3
        x(i)=edge(i,np)+param(05)*(float(k)/float(nk))*n(i)
      enddo
      call imgvalx(imgt,nx,ny,nz,center,mmppix,x,v(k),lslice)
      if(lslice.le.0)goto 99
c     if(ldebug)write(*, "('imgvalx x,v',3f7.2,f9.1)")x,v(k)
    9 continue
      do 8 k=-nk,nk			! differentiate
      u(k)=float(k)/float(nk)
    8 dv(k)=v(k+1)-v(k-1)
      if(ldebug)write(*, "('dv',10f9.1)")dv

      call polfit(u,dv,npts,coef,array,3,chisqr)
      if(ldebug)write(*, "('coef ',3g12.4)")coef
      q=.5*coef(1)/coef(2)		! find gradient extremum
      if(abs(q).gt.1.5.and.np.gt.1)then
        write(*, "('edge_trace: edge end encountered')")
        goto 98
      endif
      if(abs(q).gt.0.5)q=sign(0.5,q)
      do i=1,3				! move current edge point onto gradient extremum
        edge(i,np)=edge(i,np)-param(5)*q*n(i)
      enddo
      if(abs(q).gt.0.05.and.iter.lt.20)goto 10
      write(*, "('np, edge(np)',i6,3f7.2)")np,(edge(i,np),i=1,3)

      if(np.ge.npmax)return

      do j=1,np-2			! edge closed?
        do i=1,3
          x(i)=edge(i,np)-edge(i,j)
        enddo
        a=sqrt(dot(x,x))
        if(a.lt.param(4))then
          write(*, "('edge_trace: edge closed')")
          return
        endif
      enddo

      do i=1,3				! advance trace
        edge(i,np+1)=edge(i,np)+param(04)*t(i)
      enddo
      np=np+1
      goto 11

   99 write(*, "('edge_trace: image boundary reached')")
   98 np=np-1
      return
      end

      subroutine cross(a,b,c)
      real*4 a(3),b(3),c(3)
      c(1)= a(2)*b(3)-a(3)*b(2)
      c(2)=-a(1)*b(3)+a(3)*b(1)
      c(3)= a(1)*b(2)-a(2)*b(1)
      return
      end

      function dot(a,b)
      real*4 a(3),b(3)
      dot=a(1)*b(1)+a(2)*b(2)+a(3)*b(3)
      return
      end

      function dist(a,b)
      real*4 a(3),b(3)
      dist=sqrt((a(1)-b(1))**2+(a(2)-b(2))**2+(a(3)-b(3))**2)
      return
      end

      subroutine reflect(a,b)
      real*4 a(3),b(3)
      b(1)=-a(1)
      b(2)=-a(2)
      b(3)=-a(3)
      return
      end

max_stars_repo_name: llevitis/PUP

max_issues_repo_name: llevitis/PUP

max_forks_repo_name: llevitis/PUP


text:
c f77 -e meandiff.f -L$PCMDI/lib -lezget -lcdms -ldrs -L/usr/local/netcdf-2.3.2/libsrc -lnetcdf -o meandiff
c f77 -e -g meandiff.f -L/zooks0/ktaylor/pcmdi/util/ezget -lezgetdebug -L$PCMDI/lib -lcdms -ldrs -L/usr/local/netcdf-2.3.2/libsrc -lnetcdf -o meandiff
c  IBM
c xlf meandiff.f  -L/u/williams/ezget -lezget -L/u/williams/devel/cdms/lib -lcdms -L/u/williams/drs/lib -ldrs -L/u/williams/netcdf/libsrc -lnetcdf -o meandiff
c   HP
c fort77  +U77 meandiff.f -L/afs/nersc.gov/u/ahs/u4160/ezget -lezget -L/afs/nersc.gov/software/pub/drs/lib -lcdms -ldrs -lm -L/afs/nersc.gov/projects/graphics/netcdf/2.3.2pl4/lib -lnetcdf -o meandiff
c   sun/SOLARIS
c f77 -e meandiff.f -L/pcmdi/cirrus/ktaylor/ezget/solaris -lezget -L/pcmdi/cirrus/drs_solaris -lcdms -ldrs -L/pcmdi/cirrus/netcdf-2.3.2/libsrc -lnetcdf -o meandiff
c   SGI
c f77 meandiff.f -L/pcmdi/ktaylor/pcmdi/util/ezget/sgi -lezget -L/usr/local/cdms/lib -lcdms -L/usr/local/drs/drslib -ldrs -L/usr/local/lib -lnetcdf -o meandiff

      program meandiff

c   This program
c
c    ** retrieves observed and model-simulated monthly surface air
c       temperature for the year 1979 (first year of AMIP run).
c    ** determines which grid cells are not missing any data and then
c       for these grid cells
c    ** computes the global average, annual mean difference between the
c       observed and model-simulated temperatures.
c
c       The data contributing to the averages are area-weighted and also
c       weighted by the number of days in each month.

      parameter (nlon=100, nlat=50, nmon=12, n4=0)

      real datamodl(nlon,nlat,nmon), wtsmodl(nlon,nlat,nmon),
     &     dataobs(nlon,nlat,nmon), wtsobs(nlon,nlat,nmon),
     &     wtobs(nlon,nlat)
      double precision asum, wtsum
      integer lons, lats, mons, i4, i, j, m, mm

c     -------------------------------------------
c     Initialize EzGet and define "missing" data value:

      call initget
      call defmisc('input missing value', 'real', 1.0e20)

c     -------------------------------------------
c     Define variable 1 as the observed surface air temperature and
c     indicate path/filename where data are stored.

      call defvar(1, 'tas', '/scratch/staff/lisa/amip_obs/nasa-amip_t')

c     -------------------------------------------
c     Define domain for variable 1.

      call defdim(1, 1, 'longitude', 'width', 'range',
     &                                             -180.0, 180.0, 360.0)
      call defdim(1, 2, 'latitude', 'cosine', 'range', -90.0, 90.0, 0.0)
      call defdim(1, 3, 'time',    'month', 'nearest', 109., 120., 0.0)

c     -------------------------------------------
c     Extract variable 1 (observed field).

      lons = 0
      lats = 0
      mons = 0
      i4 = 0
      call getdata(1,  nlon,nlat,nmon,n4,  lons,lats,mons,i4,
     &             wtsobs, dataobs)

c     -------------------------------------------
c     Define variable 2 as the modeled surface air temperature and
c     indicate path/filename where data are stored.

      call defvar(2, 'tas', '/amipsp/drs/tas/tas_bmr')

c     -------------------------------------------
c     Define domain for variable 2.

      call defdim(2, 1, 'longitude', 'width', 'range',
     &                                             0.0, 0.0, 360.0)
      call defdim(2, 2, 'latitude', 'gaussian', 'range',
     &                                              0.0, 0.0, 0.0)
      call defdim(2, 3, 'time',    'month', 'nearest', 109., 120., 0.0)

c     -------------------------------------------
c      Regrid model output to observed grid

      call defregrd(2, 'to', 1, 'area-weighted', 0,0.0,0.0, 0,0.0,0.0)

c     -------------------------------------------
c     Extract variable 2 (modeled field).

      call getdata(2,  nlon,nlat,nmon,n4,  lons,lats,mons,i4,
     &             wtsmodl, datamodl)

c       -----------------------------------------
c       Compute annual mean at each grid cell, but only if all 12 months
c       of data are available for both the model and the observations.

c       If data is available, weight by the area of the grid cell and
c       the number of days in a month.

        do 300 j=1,lats
          do 200 i=1,lons

            mm = 0
            do 100 m=1,mons
              if (wtsmodl(i,j,m)*wtsobs(i,j,m) .gt. 0.0) mm = mm + 1
  100       continue

            if (mm .eq. 12) then
              wtobs(i,j) = wtsobs(i,j,1)
            else
              wtobs(i,j) = 0.0
            endif

  200     continue
  300   continue

c     -------------------------------------------
c     Compute global average, annual mean difference

      asum = 0.0
      wtsum = 0.0

      do 600 m=1,nmon
        do 500 j=1,nlat
          do 400 i=1,nlon

            asum = asum + wtobs(i,j)*(datamodl(i,j,m)-dataobs(i,j,m))
            wtsum = wtsum + wtobs(i,j)

  400     continue
  500   continue
  600 continue

c     -----------------------------------------
c     Write area-weighted annual mean difference between model and
c     observed fields.

      if (wtsum .gt. 0.0) then

        asum = asum/wtsum

        write(*,'("annually-averaged fraction of globe with data: ",
     &            f12.7)') wtsum
        write(*,'("global, annual mean difference: ", f12.3)') asum

      else

        write(*,'(
     &  " data missing everywhere for at least 1 month of the year")')

      endif

c     -------------------------------------------
c     Close all files opened by EzGet.

      call closeget

      end

max_stars_repo_name: UV-CDAT/EzGet

max_issues_repo_name: UV-CDAT/EzGet

max_forks_repo_name: UV-CDAT/EzGet


text:
MODULE m_potl0
! ******************************************************************
! evaluate the xc-potential vxc for charge density and its
! gradients,dens,... only for nonmagnetic.
! ******************************************************************
CONTAINS
   SUBROUTINE potl0(xcpot,jspins,dx,rad,dens, &
                    vxc)

      USE m_grdchlh
      USE m_mkgl0
      USE m_types
      IMPLICIT NONE

      CLASS(t_xcpot),intent(in)::xcpot
      INTEGER, INTENT (IN) :: jspins
      REAL,    INTENT (IN) :: dx
      REAL,    INTENT (IN) :: rad(:),dens(:,:)
      REAL,    INTENT (OUT):: vxc(:,:)

!     .. previously untyped names ..
      INTEGER,PARAMETER :: ndvgrd=6

      TYPE(t_gradients)::grad

      INTEGER i,ispin,msh
      REAL, ALLOCATABLE :: drr(:,:),ddrr(:,:)

      REAL              :: vx(size(vxc,1),jspins)

      msh = size(rad)
      ALLOCATE ( drr(msh,jspins),ddrr(msh,jspins))
!
!-->  evaluate gradients of dens.
!
      CALL xcpot%alloc_gradients(msh,jspins,grad)
      DO ispin = 1, jspins
         CALL grdchlh(1,1,msh,dx,rad,dens(1:1,ispin),ndvgrd,&
                     drr(1:1,ispin),ddrr(1:1,ispin))
      ENDDO

      CALL mkgl0(jspins,rad,dens,drr,ddrr,&
                 grad)
!
! --> calculate the potential.
!
      CALL xcpot%get_vxc(jspins, dens(:msh,:), vxc, vx, grad)

      DEALLOCATE ( drr,ddrr )
   END SUBROUTINE potl0
END MODULE m_potl0

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text: #redirect Town History

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	SUBROUTINE DM_RCLH  ( iflno, ipos, iheadr, iret )
C************************************************************************
C* DM_RCLH								*
C*									*
C* This subroutine reads a column header from a DM file.		*
C*									*
C* DM_RCLH  ( IFLNO, IPOS, IHEADR, IRET )				*
C*									*
C* Input parameters:							*
C*	IFLNO		INTEGER		File number			*
C*	IPOS		INTEGER		Location			*
C*									*
C* Output parameters:							*
C*	IHEADR (*)	INTEGER		Header array			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -4 = file is not open		*
C*					 -9 = invalid column		*
C*					-11 = undefined header		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	4/87						*
C* m. gamazaychikov/CWS 04/11   Add code for A2DB connectivity          *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'dmcmn.cmn'
        INCLUDE         'dbcmn.cmn'
C
	INTEGER		iheadr (*)
C
        INTEGER         istid(2), intdtf (3)
C------------------------------------------------------------------------
C*      For A2DB requests - set the needed information.
C
        IF (  dbread ) THEN
           IF ( dbdatasrc .eq. 'grid' ) THEN
C
C*            For grid data - set the time.
C
              CALL ST_LSTR  ( dbtime, ldbtime, ier )
              IF ( ldbtime .gt. 0 ) THEN
                 CALL TG_CTOI  ( dbtime, intdtf, ier )
                 CALL TG_ITOF  ( intdtf, iheadr, ier )
              ELSE
                 iret = -9
                 RETURN
              END IF
           ELSE
C
C*            For point data - set the station id.
C
              CALL ST_LSTR (dbstid, lstr, ier)
              IF ( lstr .eq. 4 .and. dbstid(1:1) .eq. 'K' ) THEN
                 ist=2
              ELSE
                 ist = 1
              END IF
              IF ( dbdatasrc .eq. 'metar'  ) THEN
                CALL ST_STOI (dbstid(ist:lstr), 8, nv, istid, ier )
                iheadr(1) = istid(1)
                iheadr(2) = IMISSD
              ELSE IF ( dbdatasrc .eq. 'bufrua' ) THEN
                iheadr(1) = IMISSD
                CALL ST_NUMB (dbstid(ist:lstr), istid2, ier )
                iheadr(2) = istid2
              END IF
              iheadr(3) = dbstlt
              iheadr(4) = dbstln
              iheadr(5) = dbstel
              iheadr(8) = IMISSD
           END IF
           iret = 0
           RETURN
        END IF
C
C*	Check that the file number is valid.
C
	CALL DM_CHKF ( iflno, iret )
	IF  ( iret .ne. 0 ) RETURN
C
C*	Check for valid position.
C
	IF ( (ipos .le. 0) .or. (ipos .gt. kcol (iflno))) THEN
	    iret = -9
	    DO  i = 1, kckeys ( iflno )
		iheadr (i) = IMISSD
	    END DO
	  ELSE
C
C*	    Check that this header is defined.
C
	    jloc = ipos + krow ( iflno )
	    IF  ( kheadr ( 0, jloc, iflno ) .ne. IMISSD )  THEN
C
C*		Retrieve row header.
C
		DO  i = 1, kckeys (iflno)
		    iheadr (i) = kheadr ( i, jloc, iflno )
		END DO
C
C*		Set error return.
C
	      ELSE
		iret = -11
		DO  i = 1, kckeys ( iflno )
		    iheadr (i) = IMISSD
		END DO
	    END IF
	END IF
C*
	RETURN
	END

max_stars_repo_name: sgdecker/pygempak

max_issues_repo_name: sgdecker/pygempak

max_forks_repo_name: sgdecker/pygempak


text:
c******************************************************
c  Fitting output with cilinder. The cilinder radius
c  equals sig/2.
c******************************************************
      SUBROUTINE RES_FORM(lidd,iidd,nlidd,iin,nmm,ff,sig)
      IMPLICIT NONE
      include "tomo.h"
      include "line.h"
      integer*4 nmm
      real*4    ff(nmm),sig
      integer*4 lidd(nmm),iidd(nmm)
      integer*4 i,i1,i2,i3,iin,iou,nlidd,ilo,j
      real*8    s,s1,x1(3),del
c*
      nlidd=0
      do i=1,nmm
      ff(i)=0.0
      iidd(i)=0
      enddo
c* Get the central point coordinates, x1(i)
      iou=ioutr(iin)
      i1=icrpnt(iou)/1000000
      i3=icrpnt(iou)-i1*1000000
      i2=i3/1000
      i3=i3-i2*1000
      s1=dr(i1)*dr(i1)+dr(i2)*dr(i2)+dr(i3)*dr(i3)
      s=DSQRT(s1)
      x1(1)=dr(i1)/s
      x1(2)=dr(i2)/s
      x1(3)=dr(i3)/s
c more wide/close area
      do i=1,nall
        s1=0.0d0
        do j=1,3
          s1=s1+dpx(j,i)*dpx(j,iou)
        enddo
      if(DABS(s1).gt.1.0d0) s1=DSIGN(1.0d0,s1)
      del=DACOS(s1)*DBLE(R)
      if(del.le.sig) then
        ilo=ilocr(i)
        if(ilo.ne.0) then
              nlidd=nlidd+1
              lidd(nlidd)=ilo
              if(del.le.sig/2.0) iidd(nlidd)=1
        ff(ilo)=1.0d0
        endif
      endif
      enddo
      return
      end

max_stars_repo_name: NoiseCIEI/RayTomo

max_issues_repo_name: NoiseCIEI/RayTomo

max_forks_repo_name: NoiseCIEI/RayTomo


!
! NESDIS_MW_SnowEM_LUT
!
! Module containing the parameters related to microwave snow emissivity model
!
!
! CREATION HISTORY:
!       Written by:     Banghua Yan, 03-Jun-2005
!                       banghua.yan@noaa.gov
!                       Fuzhong Weng
!                       fuzhong.weng@noaa.gov
!

MODULE NESDIS_MW_SnowEM_LUT

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module use
  USE CSEM_Type_Kinds, ONLY: fp => CSEM_fp
  USE CSEM_Interpolation, ONLY: NPTS        , &
                                CSEM_LPoly_type  , &
                                CSEM_find_index  , &
                                CSEM_interp_1D   , &
                                CSEM_Clear_LPoly , &
                                CSEM_LPoly
  ! Disable implicit typing
  IMPLICIT NONE

  ! ------------
  ! Visibilities
  ! ------------
  PRIVATE

  ! -----------------
  ! Module parameters
  ! -----------------
  ! Version Id for the module
  CHARACTER(*), PARAMETER :: MODULE_VERSION_ID = &
  '$Id: NESDIS_SnowEM_Parameters.f90 21141 2012-09-14 17:40:43Z paul.vandelst@noaa.gov $'

  INTEGER, PUBLIC, PARAMETER :: N_MWSNOW_TYPES = 16
 ! Snow types
  INTEGER, PUBLIC, PARAMETER :: INVALID_SNOW_TYPE   = -999
  INTEGER, PUBLIC, PARAMETER :: WET_SNOW            =  1
  INTEGER, PUBLIC, PARAMETER :: GRASS_AFTER_SNOW    =  2
  INTEGER, PUBLIC, PARAMETER :: RS_SNOW_A           =  3
  INTEGER, PUBLIC, PARAMETER :: POWDER_SNOW         =  4
  INTEGER, PUBLIC, PARAMETER :: RS_SNOW_B           =  5
  INTEGER, PUBLIC, PARAMETER :: RS_SNOW_C           =  6
  INTEGER, PUBLIC, PARAMETER :: RS_SNOW_D           =  7
  INTEGER, PUBLIC, PARAMETER :: THIN_CRUST_SNOW     =  8
  INTEGER, PUBLIC, PARAMETER :: RS_SNOW_E           =  9
  INTEGER, PUBLIC, PARAMETER :: BOTTOM_CRUST_SNOW_A = 10
  INTEGER, PUBLIC, PARAMETER :: SHALLOW_SNOW        = 11
  INTEGER, PUBLIC, PARAMETER :: DEEP_SNOW           = 12
  INTEGER, PUBLIC, PARAMETER :: CRUST_SNOW          = 13
  INTEGER, PUBLIC, PARAMETER :: MEDIUM_SNOW         = 14
  INTEGER, PUBLIC, PARAMETER :: BOTTOM_CRUST_SNOW_B = 15
  INTEGER, PUBLIC, PARAMETER :: THICK_CRUST_SNOW    = 16

  CHARACTER(LEN=20), PUBLIC, PARAMETER, DIMENSION(N_MWSNOW_TYPES):: &
    NESDIS_SNOW_TYPE_LIST = (/ &
   'WET_SNOW            ', &
   'GRASS_AFTER_SNOW    ', &
   'RS_SNOW_A           ', &
   'POWDER_SNOW         ', &
   'RS_SNOW_B           ', &
   'RS_SNOW_C           ', &
   'RS_SNOW_D           ', &
   'THIN_CRUST_SNOW     ', &
   'RS_SNOW_E           ', &
   'BOTTOM_CRUST_SNOW_A ', &
   'SHALLOW_SNOW        ', &
   'DEEP_SNOW           ', &
   'CRUST_SNOW          ', &
   'MEDIUM_SNOW         ', &
   'BOTTOM_CRUST_SNOW_B ', &
   'THICK_CRUST_SNOW    ' /)

  INTEGER, PUBLIC, PARAMETER :: N_MWSNOW_FREQUENCY  = 10
  INTEGER, PUBLIC, PARAMETER :: N_AMSRE_SNOW_FREQ   = 7

  REAL(fp),PARAMETER, DIMENSION(N_AMSRE_SNOW_FREQ):: FREQUENCY_AMSRE = &
    (/ 6.925_fp, 10.65_fp, 18.7_fp,23.8_fp, 36.5_fp, 89.0_fp, 150._fp/)

 ! Define sixteen MW H-POL emissivity spectra for AMSRE ALGORITHMS
  REAL(fp),PARAMETER, DIMENSION(N_AMSRE_SNOW_FREQ,N_MWSNOW_TYPES) ::  &
  NESDIS_AMSRE_SnowEM_LUT_H = RESHAPE( (/ &
    0.93_fp, 0.92_fp, 0.93_fp, 0.94_fp, 0.93_fp, 0.93_fp, 0.90_fp, &
    0.91_fp, 0.90_fp, 0.90_fp, 0.90_fp, 0.91_fp, 0.90_fp, 0.85_fp, &
    0.85_fp, 0.85_fp, 0.84_fp, 0.84_fp, 0.82_fp, 0.80_fp, 0.80_fp, &
    0.90_fp, 0.90_fp, 0.92_fp, 0.92_fp, 0.90_fp, 0.80_fp, 0.79_fp, &
    0.82_fp, 0.81_fp, 0.77_fp, 0.76_fp, 0.74_fp, 0.74_fp, 0.74_fp, &
    0.84_fp, 0.83_fp, 0.80_fp, 0.78_fp, 0.77_fp, 0.75_fp, 0.69_fp, &
    0.77_fp, 0.77_fp, 0.76_fp, 0.75_fp, 0.73_fp, 0.71_fp, 0.71_fp, &
    0.95_fp, 0.94_fp, 0.95_fp, 0.94_fp, 0.89_fp, 0.75_fp, 0.65_fp, &
    0.73_fp, 0.73_fp, 0.74_fp, 0.72_fp, 0.71_fp, 0.68_fp, 0.67_fp, &
    0.88_fp, 0.87_fp, 0.86_fp, 0.85_fp, 0.80_fp, 0.68_fp, 0.63_fp, &
    0.86_fp, 0.84_fp, 0.80_fp, 0.78_fp, 0.72_fp, 0.62_fp, 0.57_fp, &
    0.87_fp, 0.85_fp, 0.83_fp, 0.80_fp, 0.77_fp, 0.68_fp, 0.62_fp, &
    0.82_fp, 0.78_fp, 0.74_fp, 0.71_fp, 0.67_fp, 0.64_fp, 0.64_fp, &
    0.90_fp, 0.90_fp, 0.89_fp, 0.88_fp, 0.83_fp, 0.53_fp, 0.48_fp, &
    0.87_fp, 0.85_fp, 0.84_fp, 0.82_fp, 0.74_fp, 0.53_fp, 0.49_fp, &
    0.85_fp, 0.84_fp, 0.83_fp, 0.81_fp, 0.79_fp, 0.51_fp, 0.46_fp  &
  /),(/N_AMSRE_SNOW_FREQ,N_MWSNOW_TYPES/))

 ! Define sixteen MW V-POL emissivity spectra for AMSRE ALGORITHMS

  REAL(fp),PARAMETER, DIMENSION(N_AMSRE_SNOW_FREQ,N_MWSNOW_TYPES) ::  &
  NESDIS_AMSRE_SnowEM_LUT_V = RESHAPE( (/ &
    0.96_fp, 0.94_fp, 0.96_fp, 0.95_fp, 0.94_fp, 0.94_fp, 0.91_fp, &
    0.96_fp, 0.94_fp, 0.95_fp, 0.96_fp, 0.96_fp, 0.92_fp, 0.87_fp, &
    0.99_fp, 0.97_fp, 0.96_fp, 0.96_fp, 0.93_fp, 0.87_fp, 0.87_fp, &
    0.98_fp, 0.97_fp, 0.99_fp, 0.98_fp, 0.96_fp, 0.84_fp, 0.83_fp, &
    0.97_fp, 0.95_fp, 0.93_fp, 0.92_fp, 0.89_fp, 0.84_fp, 0.84_fp, &
    1.00_fp, 0.97_fp, 0.96_fp, 0.94_fp, 0.91_fp, 0.84_fp, 0.78_fp, &
    0.99_fp, 0.96_fp, 0.93_fp, 0.90_fp, 0.86_fp, 0.80_fp, 0.80_fp, &
    0.98_fp, 0.97_fp, 0.98_fp, 0.97_fp, 0.92_fp, 0.77_fp, 0.67_fp, &
    0.98_fp, 0.95_fp, 0.90_fp, 0.86_fp, 0.82_fp, 0.74_fp, 0.73_fp, &
    0.96_fp, 0.95_fp, 0.95_fp, 0.93_fp, 0.87_fp, 0.71_fp, 0.66_fp, &
    0.97_fp, 0.95_fp, 0.94_fp, 0.90_fp, 0.84_fp, 0.68_fp, 0.63_fp, &
    0.96_fp, 0.94_fp, 0.92_fp, 0.90_fp, 0.85_fp, 0.77_fp, 0.71_fp, &
    0.98_fp, 0.96_fp, 0.93_fp, 0.90_fp, 0.81_fp, 0.71_fp, 0.71_fp, &
    0.99_fp, 0.97_fp, 0.98_fp, 0.96_fp, 0.92_fp, 0.57_fp, 0.52_fp, &
    1.00_fp, 0.97_fp, 0.97_fp, 0.95_fp, 0.86_fp, 0.58_fp, 0.54_fp, &
    0.98_fp, 0.96_fp, 0.96_fp, 0.94_fp, 0.89_fp, 0.56_fp, 0.51_fp  &
  /),(/N_AMSRE_SNOW_FREQ,N_MWSNOW_TYPES/))

 ! Define sixteen MW weighted emissivity spectra for AMSRE ALGORITHMS

  REAL(fp),PARAMETER, DIMENSION(N_AMSRE_SNOW_FREQ,N_MWSNOW_TYPES) :: &
  NESDIS_AMSRE_SnowEM_LUT_M  =   RESHAPE( (/ &
    0.91_fp, 0.93_fp, 0.94_fp, 0.95_fp, 0.95_fp, 0.93_fp, 0.93_fp, &
    0.91_fp, 0.92_fp, 0.91_fp, 0.90_fp, 0.91_fp, 0.91_fp, 0.91_fp, &
    0.90_fp, 0.89_fp, 0.88_fp, 0.87_fp, 0.86_fp, 0.82_fp, 0.82_fp, &
    0.92_fp, 0.93_fp, 0.94_fp, 0.94_fp, 0.92_fp, 0.80_fp, 0.80_fp, &
    0.87_fp, 0.86_fp, 0.83_fp, 0.80_fp, 0.79_fp, 0.77_fp, 0.77_fp, &
    0.89_fp, 0.88_fp, 0.85_fp, 0.84_fp, 0.83_fp, 0.79_fp, 0.79_fp, &
    0.84_fp, 0.83_fp, 0.82_fp, 0.80_fp, 0.78_fp, 0.72_fp, 0.72_fp, &
    0.95_fp, 0.96_fp, 0.96_fp, 0.95_fp, 0.91_fp, 0.75_fp, 0.75_fp, &
    0.80_fp, 0.80_fp, 0.80_fp, 0.79_fp, 0.75_fp, 0.70_fp, 0.70_fp, &
    0.91_fp, 0.90_fp, 0.89_fp, 0.87_fp, 0.82_fp, 0.69_fp, 0.69_fp, &
    0.90_fp, 0.89_fp, 0.85_fp, 0.82_fp, 0.76_fp, 0.65_fp, 0.65_fp, &
    0.89_fp, 0.88_fp, 0.86_fp, 0.83_fp, 0.78_fp, 0.70_fp, 0.70_fp, &
    0.88_fp, 0.86_fp, 0.80_fp, 0.75_fp, 0.69_fp, 0.67_fp, 0.67_fp, &
    0.96_fp, 0.97_fp, 0.92_fp, 0.87_fp, 0.72_fp, 0.50_fp, 0.50_fp, &
    0.93_fp, 0.94_fp, 0.89_fp, 0.85_fp, 0.74_fp, 0.48_fp, 0.48_fp, &
    0.88_fp, 0.88_fp, 0.87_fp, 0.85_fp, 0.77_fp, 0.52_fp, 0.52_fp  &
  /),(/N_AMSRE_SNOW_FREQ,N_MWSNOW_TYPES/))

  REAL(fp), PARAMETER, DIMENSION(N_MWSNOW_FREQUENCY) :: NESDIS_MWSNOW_FREQUENCY= (/  &
    4.90_fp,6.93_fp,10.65_fp,18.7_fp,23.8_fp,31.4_fp,50.3_fp,52.5_fp,89._fp,150._fp /)

  REAL(fp),PARAMETER :: NESDIS_MWSNOW_Emiss_LUT(N_MWSNOW_FREQUENCY,N_MWSNOW_TYPES) = &
    RESHAPE( (/ &
    0.87_fp,0.89_fp,0.91_fp,0.93_fp,0.94_fp,0.94_fp,0.94_fp,0.93_fp,0.92_fp,0.90_fp, &
    0.91_fp,0.91_fp,0.92_fp,0.91_fp,0.90_fp,0.90_fp,0.91_fp,0.91_fp,0.91_fp,0.86_fp, &
    0.90_fp,0.89_fp,0.88_fp,0.87_fp,0.86_fp,0.86_fp,0.85_fp,0.85_fp,0.82_fp,0.82_fp, &
    0.91_fp,0.91_fp,0.93_fp,0.93_fp,0.93_fp,0.93_fp,0.89_fp,0.88_fp,0.79_fp,0.79_fp, &
    0.90_fp,0.89_fp,0.88_fp,0.85_fp,0.84_fp,0.83_fp,0.83_fp,0.82_fp,0.79_fp,0.73_fp, &
    0.90_fp,0.89_fp,0.86_fp,0.82_fp,0.80_fp,0.79_fp,0.78_fp,0.78_fp,0.77_fp,0.77_fp, &
    0.88_fp,0.86_fp,0.85_fp,0.80_fp,0.78_fp,0.77_fp,0.77_fp,0.76_fp,0.72_fp,0.72_fp, &
    0.93_fp,0.94_fp,0.96_fp,0.96_fp,0.95_fp,0.93_fp,0.87_fp,0.86_fp,0.74_fp,0.65_fp, &
    0.87_fp,0.86_fp,0.84_fp,0.80_fp,0.76_fp,0.76_fp,0.75_fp,0.75_fp,0.70_fp,0.69_fp, &
    0.87_fp,0.86_fp,0.83_fp,0.77_fp,0.73_fp,0.68_fp,0.66_fp,0.66_fp,0.68_fp,0.67_fp, &
    0.89_fp,0.89_fp,0.88_fp,0.87_fp,0.86_fp,0.82_fp,0.77_fp,0.76_fp,0.69_fp,0.64_fp, &
    0.88_fp,0.87_fp,0.86_fp,0.83_fp,0.81_fp,0.77_fp,0.74_fp,0.73_fp,0.69_fp,0.64_fp, &
    0.86_fp,0.86_fp,0.86_fp,0.85_fp,0.82_fp,0.78_fp,0.69_fp,0.68_fp,0.51_fp,0.47_fp, &
    0.89_fp,0.88_fp,0.87_fp,0.83_fp,0.80_fp,0.75_fp,0.70_fp,0.70_fp,0.64_fp,0.60_fp, &
    0.91_fp,0.92_fp,0.93_fp,0.88_fp,0.84_fp,0.76_fp,0.66_fp,0.64_fp,0.48_fp,0.44_fp, &
    0.94_fp,0.95_fp,0.97_fp,0.91_fp,0.86_fp,0.74_fp,0.63_fp,0.63_fp,0.50_fp,0.45_fp  &
  /), (/N_MWSNOW_FREQUENCY, N_MWSNOW_TYPES/))

  ! --------------------------------------
  ! Structure definitions to hold forward
  ! variables across FWD, TL, and AD calls
  ! --------------------------------------
  TYPE :: iVar_type
    PRIVATE
    ! The interpolating polynomials
    TYPE(CSEM_LPoly_type) :: xlp
    ! The LUT interpolation indices
    INTEGER :: i1, i2
    ! The LUT interpolation boundary check
    LOGICAL :: x_outbound
    ! The interpolation input
    REAL(fp) :: x_int
    ! The data to be interpolated
    REAL(fp) :: x(NPTS)
  END TYPE iVar_type

  PUBLIC :: Load_MW_SnowEM_LUT

  INTERFACE Load_MW_SnowEM_LUT
    MODULE PROCEDURE Set_MW_SnowEM_LUT
    MODULE PROCEDURE Set_AMSRE_SnowEM_LUT
  END INTERFACE Load_MW_SnowEM_LUT

CONTAINS

  FUNCTION Set_MW_SnowEM_LUT( &
      frequency,         & ! INPUT, wavelength in micrometer
      emissivity)        & ! OUTPUT, surface emissivity (0 - 1)
      RESULT (Error_Status)

    REAL(fp), INTENT( IN )  :: frequency(:)
    REAL(fp), INTENT( OUT ) :: emissivity(:,:)

    INTEGER :: snow_type
    INTEGER :: i,Error_Status

    ! --------------------------------------------------  !
    !       internal variables                            !
    ! --------------------------------------------------  !
    Error_Status = 0
    DO i = 1, size(frequency)
      DO snow_type = 1, N_MWSNOW_TYPES
        Error_Status = Interp_MW_SnowEM_LUT( &
          frequency(i),               & ! INPUT, wavelength in micrometer
          emissivity(snow_type,i),    & ! OUTPUT, surface emissivity (0 - 1)
          snow_type )
      ENDDO
    ENDDO
  END FUNCTION  Set_MW_SnowEM_LUT

  FUNCTION Set_AMSRE_SnowEM_LUT( &
      frequency,           & ! INPUT, wavelength in micrometer
      emissivity_m,        & ! OUTPUT, surface emissivity (0 - 1)
      emissivity_h,        & ! OUTPUT, surface emissivity (0 - 1)
      emissivity_v)        & ! OUTPUT, surface emissivity (0 - 1)
    RESULT (Error_Status)

    REAL(fp), INTENT( IN )  :: frequency(:)
    REAL(fp), INTENT( OUT ) :: emissivity_m(:,:)
    REAL(fp), INTENT( OUT ) :: emissivity_h(:,:)
    REAL(fp), INTENT( OUT ) :: emissivity_v(:,:)

    INTEGER :: snow_type
    INTEGER :: i,Error_Status

    ! --------------------------------------------------  !
    !       internal variables                            !
    ! --------------------------------------------------  !
    Error_Status = 0

    DO i = 1, size(frequency)
      DO snow_type = 1, N_MWSNOW_TYPES
        Error_Status = Interp_AMSRE_SnowEM_LUT( &
          frequency(i),                 & ! INPUT, wavelength in micrometer
          emissivity_m(snow_type,i),    & ! OUTPUT, surface emissivity (0 - 1)
          emissivity_h(snow_type,i),    & ! OUTPUT, surface emissivity (0 - 1)
          emissivity_v(snow_type,i),    & ! OUTPUT, surface emissivity (0 - 1)
          snow_type )
      ENDDO
    ENDDO
  END FUNCTION  Set_AMSRE_SnowEM_LUT

  FUNCTION Interp_MW_SnowEM_LUT( &
      frequency,         & ! INPUT, wavelength in micrometer
      emissivity,        & ! OUTPUT, surface emissivity (0 - 1)
      snow_type )        & ! INPUT, surface type (1 - 16)
    RESULT (Error_Status)

    INTEGER , INTENT( IN )  :: snow_type
    REAL(fp), INTENT( IN )  :: frequency
    REAL(fp), INTENT( OUT ) :: emissivity
    INTEGER :: Error_Status
    TYPE(iVar_type) :: iVar

    ! --------------------------------------------------  !
    !       internal variables                            !
    ! --------------------------------------------------  !

    iVar%x_int = MAX(MIN(NESDIS_MWSNOW_FREQUENCY(N_MWSNOW_FREQUENCY),&
                 frequency), NESDIS_MWSNOW_FREQUENCY(1))
    CALL CSEM_find_index(NESDIS_MWSNOW_FREQUENCY, &
                 iVar%x_int, iVar%i1, iVar%i2, iVar%x_outbound)
    iVar%x = NESDIS_MWSNOW_FREQUENCY(iVar%i1:iVar%i2)

    ! Calculate the interpolating polynomial
    CALL CSEM_LPoly( iVar%x, iVar%x_int, & ! Input
                iVar%xlp            )      ! Output

    ! Perform Interpolation
    CALL CSEM_interp_1D( &
          NESDIS_MWSNOW_Emiss_LUT(iVar%i1:iVar%i2, Snow_Type), &
          iVar%xlp, Emissivity)
    Error_Status = 0
  END FUNCTION Interp_MW_SnowEM_LUT

  FUNCTION Interp_AMSRE_SnowEM_LUT( &
      frequency,         & ! INPUT, wavelength in micrometer
      emissivity_m,      & ! OUTPUT, surface emissivity (0 - 1)
      emissivity_h,      & ! OUTPUT, surface emissivity (0 - 1)
      emissivity_v,      & ! OUTPUT, surface emissivity (0 - 1)
      snow_type )        & ! INPUT, surface type (1 - 16)
    RESULT (Error_Status)

    INTEGER , INTENT( IN )  :: snow_type
    REAL(fp), INTENT( IN )  :: frequency
    REAL(fp), INTENT( OUT ) :: emissivity_m, emissivity_h, emissivity_v
    INTEGER :: Error_Status
    TYPE(iVar_type) :: iVar

    ! --------------------------------------------------  !
    !       internal variables                            !
    ! --------------------------------------------------  !

    iVar%x_int = MAX(MIN(FREQUENCY_AMSRE(N_AMSRE_SNOW_FREQ),&
                 frequency), FREQUENCY_AMSRE(1))
    CALL CSEM_find_index(FREQUENCY_AMSRE, &
                 iVar%x_int, iVar%i1, iVar%i2, iVar%x_outbound)
    iVar%x = FREQUENCY_AMSRE(iVar%i1:iVar%i2)

    ! Calculate the interpolating polynomial
    CALL CSEM_LPoly( iVar%x, iVar%x_int, & ! Input
                iVar%xlp            )      ! Output

    ! Perform Interpolation
    CALL CSEM_interp_1D( &
          NESDIS_AMSRE_SnowEM_LUT_M(iVar%i1:iVar%i2, Snow_Type), &
          iVar%xlp, Emissivity_m)
    CALL CSEM_interp_1D( &
          NESDIS_AMSRE_SnowEM_LUT_H(iVar%i1:iVar%i2, Snow_Type), &
          iVar%xlp, Emissivity_h)
    CALL CSEM_interp_1D( &
          NESDIS_AMSRE_SnowEM_LUT_V(iVar%i1:iVar%i2, Snow_Type), &
          iVar%xlp, Emissivity_v)
    Error_Status = 0
  END FUNCTION Interp_AMSRE_SnowEM_LUT

END MODULE NESDIS_MW_SnowEM_LUT