      subroutine actions (ob_cur, ob_num, idtbl)
      use conditional_module
      use climate_module
      use time_module
      use hru_module, only : fertno3, fertnh3, fertorgn, fertorgp, fertsolp,   &
        ihru, ipl, isol, ndeat, phubase, sol_sumno3, sol_sumsolp, hru, yield
      use soil_module
      use plant_module
      use plant_data_module
      use mgt_operations_module
      use tillage_data_module
      use reservoir_module
      use sd_channel_module
      use hru_lte_module
      use basin_module
      use organic_mineral_mass_module
      use hydrograph_module
      use tiles_data_module
      use output_landscape_module
      use conditional_module
      use constituent_mass_module
      use calibration_data_module

      implicit none

      integer, intent (in)  :: ob_cur      !none     |sequential number of individual objects
      integer, intent (in)  :: ob_num      !none     |sequential number for all objects
      integer, intent (in)  :: idtbl       !none     |
      integer :: icom                      !none     |
      integer :: iac                       !none     |counter
      integer :: ial                       !none     |counter
      integer :: jj                        !none     |counter
      integer :: i                         !none     |counter
      integer :: iburn                     !none     |burn type from fire data base
      integer :: idtill                    !none     |tillage type
      integer :: ifertop                   !         |surface application fraction from chem app data base
      integer :: ifrt                      !         |fertilizer type from fert data base
      integer :: ipestop                   !         |surface application fraction from chem app data base
      integer :: ipst                      !         |pesticide type from pest data base
      integer :: iharvop                   !         |harvest operation type
      integer :: iihru                     !         |
      integer :: ilu                       !         |landuse type
      integer :: j                         !none     |counter
      integer :: iob
      integer :: idp                       !         |
      integer :: istr                      !         |
      integer :: iob_out
      integer :: inhyd                     !         |
      integer :: ihyd_in                   !         |
      integer :: icon                      !         |
      integer :: iplt_bsn
      integer :: irrop                     !         |
      integer :: igr
      integer :: ireg                      !         |
      integer :: ilum
      real :: hiad1                        !         |
      real :: irrig_m3                     !         |
      real :: amt_mm                       !         |
      real :: biomass                      !         |
      real :: frt_kg
      real :: wur                          !         |
      real :: frac                         !         |
      real :: rto                          !         |
      real :: rto1                         !         |
      real :: pest_kg                      !kg/ha    |amount of pesticide applied
      character(len=1) :: action           !         |

      do iac = 1, d_tbl%acts
        action = "n"
        do ial = 1, d_tbl%alts
          if (d_tbl%act_hit(ial) == "y" .and. d_tbl%act_outcomes(iac,ial) == "y") then
            action = "y"
            exit
          end if
        end do

        if (action == "y") then
          select case (d_tbl%act(iac)%typ)

          !irrigation demand - hru action
          case ("irr_demand")
            ipl = 1
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            irrop = d_tbl%act_typ(iac)
            irrig(j)%demand = irrop_db(irrop)%amt_mm * hru(j)%area_ha / 1000.       ! ha-m = mm * ha / 1000.

            !! if unlimited source, set irrigation applied directly to hru
            if (d_tbl%act(iac)%file_pointer == "unlim") then
              irrig(j)%applied = irrop_db(irrop)%amt_mm * irrop_db(irrop)%eff * (1. - irrop_db(irrop)%surq)
              irrig(j)%runoff = irrop_db(irrop)%amt_mm * irrop_db(irrop)%eff * irrop_db(irrop)%surq
              !set organics and constituents from irr.ops ! irrig(j)%water =  cs_irr(j) =
              if (pco%mgtout == "y") then
                write (2612, *) j, time%yrc, time%mo, time%day, "        ", "IRRIGATE", phubase(j),  &
                    pcom(j)%plcur(ipl)%phuacc, soil(j)%sw,pl_mass(j)%tot(ipl)%m, rsd1(j)%tot(ipl)%m, &
                    sol_sumno3(j), sol_sumsolp(j), irrig(j)%applied
              end if
            else
              !! set demand for irrigation from channel, reservoir or aquifer
              if (pco%mgtout == "y") then
                write (2612, *) j, time%yrc, time%mo, time%day, "        ", "IRRIG_DMD", phubase(j), &
                    pcom(j)%plcur(ipl)%phuacc, soil(j)%sw,pl_mass(j)%tot(ipl)%m, rsd1(j)%tot(ipl)%m, &
                    sol_sumno3(j), sol_sumsolp(j), irrig(j)%demand
              end if
            end if

          !irrigate - hru action
          case ("irrigate")
            ipl = 1
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            irrop = mgt%op4                        !irrigation amount (mm) from irr.ops data base
            irrig(j)%applied = irrop_db(irrop)%amt_mm * irrop_db(irrop)%eff * (1. - irrop_db(irrop)%surq)
            irrig(j)%runoff = irrop_db(irrop)%amt_mm * irrop_db(irrop)%surq

            irrig(j)%demand = mgt%op3 * hru(j)%area_ha / 1000.        ! ha-m = mm * ha / 1000.

            !select object type
            iob = d_tbl%act(iac)%const
            select case (d_tbl%act(iac)%option)
            case ("aqu")
              irrig_m3 = amin1 (ch_stor(j)%flo, irrig(j)%demand)
              rto = irrig_m3 / ch_stor(j)%flo                               ! ratio of water removed from aquifer volume
              rto1 = (1. - rto)
              irrig(j)%applied = irrig(j)%applied * hru(j)%area_ha / 10.    ! convert to mm for hru application
              irrig(j)%water = rto * ch_stor(j)                             ! organics in irrigation water
              aqu(j) = rto1 * ch_stor(j)                                    ! remainder stays in aquifer
              cs_irr(j) = rto * cs_aqu(j)                                   ! constituents in irrigation water
              cs_aqu(j) = rto1 * cs_aqu(j)                                  ! remainder stays in aquifer

            case ("cha")
              irrig_m3 = amin1 (ch_stor(j)%flo, irrig(j)%demand)
              rto = irrig_m3 / ch_stor(j)%flo                               ! ratio of water removed from channel volume
              rto1 = (1. - rto)
              irrig(j)%applied = irrig(j)%applied * hru(j)%area_ha / 10.    ! convert to mm for hru application
              irrig(j)%water = rto * ch_stor(j)                             ! organics in irrigation water
              !ch_stor(j) = rto1 * ch_water(j)                              ! remainder stays in channel
              cs_irr(j) = rto * ch_water(j)                                 ! constituents in irrigation water
              ch_water(j) = rto1 * ch_water(j)                              ! remainder stays in channel

            case ("res")
              irrig_m3 = amin1 (res(j)%flo, irrig(j)%demand)
              rto = irrig_m3 / res(j)%flo                                   ! ratio of water removed from res volume
              rto1 = (1. - rto)
              irrig(j)%applied = irrig(j)%applied * hru(j)%area_ha / 10.    ! convert to mm for hru application
              irrig(j)%water = rto * res(j)                                 ! organics in irrigation water
              res(j) = rto1 * res(j)                                        ! remainder stays in reservoir
              cs_irr(j) = rto * res_water(j)                                ! constituents in irrigation water
              res_water(j) = rto1 * res_water(j)                            ! remainder stays in reservoir

            end select

            if (pco%mgtout == "y") then
              write (2612, *) j, time%yrc, time%mo, time%day, "        ", "IRRIGATE", phubase(j),  &
                  pcom(j)%plcur(ipl)%phuacc, soil(j)%sw,pl_mass(j)%tot(ipl)%m, rsd1(j)%tot(ipl)%m, &
                  sol_sumno3(j), sol_sumsolp(j), irrig(j)%demand
            end if

          !fertilize
          case ("fertilize")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            if (pcom(j)%dtbl(idtbl)%num_actions(iac) <= Int(d_tbl%act(iac)%const2)) then
              ipl = 1
              ifrt = d_tbl%act_typ(iac)               !fertilizer type from fert data base
              frt_kg = d_tbl%act(iac)%const           !amount applied in kg/ha
              ifertop = d_tbl%act_app(iac)            !surface application fraction from chem app data base
              call pl_fert (j, ifrt, frt_kg, ifertop)

              if (pco%mgtout == "y") then
                !write (2612, *) j, time%yrc, time%mo, time%day, chemapp_db(mgt%op4)%name, "    FERT", &
                write (2612, *) j, time%yrc, time%mo, time%day, chemapp_db(j)%name, "    FERT",       &
                  phubase(j),pcom(j)%plcur(ipl)%phuacc, soil(j)%sw, pl_mass(j)%tot(ipl)%m,            &
                  rsd1(j)%tot(ipl)%m, sol_sumno3(j), sol_sumsolp(j), frt_kg, fertno3, fertnh3,        &
                  fertorgn, fertsolp, fertorgp
              endif
              pcom(j)%dtbl(idtbl)%num_actions(iac) = pcom(j)%dtbl(idtbl)%num_actions(iac) + 1
            end if

          !tillage
          case ("till")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            if (pcom(j)%dtbl(idtbl)%num_actions(iac) <= Int(d_tbl%act(iac)%const2)) then
              idtill = d_tbl%act_app(iac)
              ipl = 1
              call mgt_newtillmix(j, 0., idtill)

              if (pco%mgtout == "y") then
                write (2612, *) j, time%yrc, time%mo, time%day, tilldb(idtill)%tillnm, "TILLAGE",    &
                    phubase(j), pcom(j)%plcur(ipl)%phuacc, soil(j)%sw, pl_mass(j)%tot(ipl)%m,        &
                    rsd1(j)%tot(ipl)%m, sol_sumno3(j), sol_sumsolp(j), tilldb(idtill)%effmix
              end if
              pcom(j)%dtbl(idtbl)%num_actions(iac) = pcom(j)%dtbl(idtbl)%num_actions(iac) + 1
            end if

          !plant
          case ("plant")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur
            icom = pcom(j)%pcomdb
            pcom(j)%days_plant = 1       !reset days since last planting

            if (pcom(j)%dtbl(idtbl)%num_actions(iac) <= Int(d_tbl%act(iac)%const2)) then
              do ipl = 1, pcom(j)%npl
                idp = pcomdb(icom)%pl(ipl)%db_num
                if (d_tbl%act(iac)%option == pcomdb(icom)%pl(ipl)%cpnm) then
                  pcom(j)%plcur(ipl)%gro = "y"
                  pcom(j)%plcur(ipl)%idorm = "n"
                if (pco%mgtout == "y") then
                  write (2612, *) j, time%yrc, time%mo, time%day, pldb(idp)%plantnm, pcomdb(icom)%name,  &
                      phubase(j), pcom(j)%plcur(ipl)%phuacc, soil(ihru)%sw,                              &
                      pl_mass(j)%tot(ipl)%m, rsd1(j)%tot_com%m, sol_sumno3(j),                           &
                      sol_sumsolp(j), pcom(j)%plg(ipl)%lai, pcom(j)%plcur(ipl)%laimx_pop
                  end if
                end if
              end do
              pcom(j)%dtbl(idtbl)%num_actions(iac) = pcom(j)%dtbl(idtbl)%num_actions(iac) + 1
            end if

          !harvest only
          case ("harvest")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            if (pcom(j)%dtbl(idtbl)%num_actions(iac) <= Int(d_tbl%act(iac)%const2)) then
              iharvop = d_tbl%act_typ(iac)
              icom = pcom(j)%pcomdb
              pcom(j)%days_harv = 1       !reset days since last harvest

              do ipl = 1, pcom(j)%npl
                biomass = pl_mass(j)%tot(ipl)%m
                if (d_tbl%act(iac)%option == pcomdb(icom)%pl(ipl)%cpnm .or. d_tbl%act(iac)%option == "all") then

                  !harvest specific type
                  select case (harvop_db(iharvop)%typ)
                  case ("biomass")
                    call mgt_harvbiomass (j, ipl, iharvop)
                  case ("grain")
                    call mgt_harvgrain (j, ipl, iharvop)
                  case ("residue")
                  case ("tree")
                  case ("tuber")
                    call mgt_harvtuber (j, ipl, iharvop)
                  end select

                  !! sum yield and num. of harvest to calc ave yields
                  pl_mass(j)%yield_tot(ipl) = pl_mass(j)%yield_tot(ipl) + pl_yield
                  pcom(j)%plcur(ipl)%harv_num = pcom(j)%plcur(ipl)%harv_num + 1

                  !! sum basin crop yields and area harvested
                  iplt_bsn = pcom(j)%plcur(ipl)%bsn_num
                  bsn_crop_yld(iplt_bsn)%area_ha = bsn_crop_yld(iplt_bsn)%area_ha + hru(j)%area_ha
                  bsn_crop_yld(iplt_bsn)%yield = bsn_crop_yld(iplt_bsn)%yield + yield * hru(j)%area_ha / 1000.
                  !! sum regional crop yields for soft calibration
                  ireg = hru(j)%crop_reg
                  do ilum = 1, plcal(ireg)%lum_num
                    if (plcal(ireg)%lum(ilum)%meas%name == mgt%op_char) then
                      plcal(ireg)%lum(ilum)%ha = plcal(ireg)%lum(ilum)%ha + hru(j)%area_ha
                      plcal(ireg)%lum(ilum)%sim%yield = plcal(ireg)%lum(ilum)%sim%yield + pl_yield%m * hru(j)%area_ha / 1000.
                    end if
                  end do

                  idp = pcom(j)%plcur(ipl)%idplt
                  if (pco%mgtout == "y") then
                    write (2612, *) j, time%yrc, time%mo, time%day,  pldb(idp)%plantnm, "HARVEST",      &
                        phubase(j), pcom(j)%plcur(ipl)%phuacc, soil(j)%sw, biomass, rsd1(j)%tot(ipl)%m, &
                        sol_sumno3(j), sol_sumsolp(j), pl_yield%m, pcom(j)%plstr(ipl)%sum_n, &
                        pcom(j)%plstr(ipl)%sum_p, pcom(j)%plstr(ipl)%sum_tmp, pcom(j)%plstr(ipl)%sum_w, &
                        pcom(j)%plstr(ipl)%sum_a
                  end if
                end if
                pcom(j)%plcur(ipl)%phuacc = 0.
              end do
              pcom(j)%dtbl(idtbl)%num_actions(iac) = pcom(j)%dtbl(idtbl)%num_actions(iac) + 1
            end if

          !kill plant
          case ("kill")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            if (pcom(j)%dtbl(idtbl)%num_actions(iac) <= Int(d_tbl%act(iac)%const2)) then
              icom = pcom(j)%pcomdb
              do ipl = 1, pcom(j)%npl
                biomass = pl_mass(j)%tot(ipl)%m
                if (d_tbl%act(iac)%option == pcomdb(icom)%pl(ipl)%cpnm .or. d_tbl%act(iac)%option == "all") then

                  call mgt_killop (j, ipl)

                  idp = pcom(j)%plcur(ipl)%idplt
                  if (pco%mgtout == "y") then
                    write (2612, *) j, time%yrc, time%mo, time%day,  pldb(idp)%plantnm, "HARV/KILL",     &
                        phubase(j), pcom(j)%plcur(ipl)%phuacc, soil(j)%sw, biomass, rsd1(j)%tot(ipl)%m,  &
                        sol_sumno3(j), sol_sumsolp(j), yield, pcom(j)%plstr(ipl)%sum_n,                  &
                        pcom(j)%plstr(ipl)%sum_p, pcom(j)%plstr(ipl)%sum_tmp, pcom(j)%plstr(ipl)%sum_w,  &
                        pcom(j)%plstr(ipl)%sum_a
                  end if
                end if
                pcom(j)%plcur(ipl)%phuacc = 0.
                phubase(j) = 0.
              end do
              pcom(j)%dtbl(idtbl)%num_actions(iac) = pcom(j)%dtbl(idtbl)%num_actions(iac) + 1
            end if

          !harvest and kill
          case ("harvest_kill")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            if (pcom(j)%dtbl(idtbl)%num_actions(iac) <= Int(d_tbl%act(iac)%const2)) then
              iharvop = d_tbl%act_typ(iac)
              icom = pcom(j)%pcomdb
              pcom(j)%days_harv = 1       !reset days since last harvest

              do ipl = 1, pcom(j)%npl
                biomass = pl_mass(j)%tot(ipl)%m
                if (d_tbl%act(iac)%option == pcomdb(icom)%pl(ipl)%cpnm .or. d_tbl%act(iac)%option == "all") then

                  !harvest specific type
                  select case (harvop_db(iharvop)%typ)
                  case ("biomass")
                    call mgt_harvbiomass (j, ipl, iharvop)
                  case ("grain")
                    call mgt_harvgrain (j, ipl, iharvop)
                  case ("residue")
                  case ("tree")
                  case ("tuber")
                    call mgt_harvtuber (j, ipl, iharvop)
                  end select

                  call mgt_killop (j, ipl)

                  !! sum yield and num. of harvest to calc ave yields
                  pl_mass(j)%yield_tot(ipl) = pl_mass(j)%yield_tot(ipl) + pl_yield
                  pcom(j)%plcur(ipl)%harv_num = pcom(j)%plcur(ipl)%harv_num + 1

                  !! sum basin crop yields and area harvested
                  iplt_bsn = pcom(j)%plcur(ipl)%bsn_num
                  bsn_crop_yld(iplt_bsn)%area_ha = bsn_crop_yld(iplt_bsn)%area_ha + hru(j)%area_ha
                  bsn_crop_yld(iplt_bsn)%yield = bsn_crop_yld(iplt_bsn)%yield + yield * hru(j)%area_ha / 1000.
                  !! sum regional crop yields for soft calibration
                  ireg = hru(j)%crop_reg
                  do ilum = 1, plcal(ireg)%lum_num
                    if (plcal(ireg)%lum(ilum)%meas%name == d_tbl%act(iac)%option) then
                      plcal(ireg)%lum(ilum)%ha = plcal(ireg)%lum(ilum)%ha + hru(j)%area_ha
                      plcal(ireg)%lum(ilum)%sim%yield = plcal(ireg)%lum(ilum)%sim%yield + pl_yield%m * hru(j)%area_ha / 1000.
                    end if
                  end do

                  idp = pcom(j)%plcur(ipl)%idplt
                  if (pco%mgtout == "y") then
                    write (2612, *) j, time%yrc, time%mo, time%day,  pldb(idp)%plantnm, "HARV/KILL",        &
                        phubase(j), pcom(j)%plcur(ipl)%phuacc, soil(j)%sw, biomass, rsd1(j)%tot(ipl)%m,     &
                        sol_sumno3(j), sol_sumsolp(j), pl_yield%m, pcom(j)%plstr(ipl)%sum_n,   &
                        pcom(j)%plstr(ipl)%sum_p, pcom(j)%plstr(ipl)%sum_tmp, pcom(j)%plstr(ipl)%sum_w,     &
                        pcom(j)%plstr(ipl)%sum_a
                  end if
                end if
                pcom(j)%plcur(ipl)%phuacc = 0.
                phubase(j) = 0.
              end do
              pcom(j)%dtbl(idtbl)%num_actions(iac) = pcom(j)%dtbl(idtbl)%num_actions(iac) + 1
            end if

          !reset rotation year
          case ("rot_reset")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur
            if (d_tbl%act(iac)%const < 1) d_tbl%act(iac)%const = 1
            pcom(j)%rot_yr = d_tbl%act(iac)%const

          !apply pesticide
          case ("pest_apply")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            if (pcom(j)%dtbl(idtbl)%num_actions(iac) <= Int(d_tbl%act(iac)%const2)) then
              ipl = 1
              ipst = d_tbl%act_typ(iac)                                     !pesticide type from fert data base
              ipestop = d_tbl%act_app(iac)                                  !surface application fraction from chem app data base
              pest_kg = d_tbl%act(iac)%const * chemapp_db(ipestop)%app_eff  !amount applied in kg/ha
              call pest_apply (j, ipst, pest_kg, ipestop)

              if (pco%mgtout == "y") then
                write (2612, *) j, time%yrc, time%mo, time%day_mo, d_tbl%act(iac)%option, "    PEST ",        &
                 phubase(j), pcom(j)%plcur(ipl)%phuacc, soil(j)%sw,pl_mass(j)%tot(ipl)%m,           &
                 rsd1(j)%tot(ipl)%m, sol_sumno3(j), sol_sumsolp(j), pest_kg
              endif
              pcom(j)%dtbl(idtbl)%num_actions(iac) = pcom(j)%dtbl(idtbl)%num_actions(iac) + 1
            end if

          case ("graze")    !! grazing operation
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            igr = d_tbl%act_typ(iac)
            graze = grazeop_db(igr)
            call pl_graze

              !if (pco%mgtout == "y") then
              !  write (2612, *) j, time%yrc, time%mo, time%day, "         ", "    GRAZE",         &
              !    phubase(j), pcom(j)%plcur(ipl)%phuacc, soil(j)%sw,pl_mass(j)%tot(ipl)%m,        &
              !    rsd1(j)%tot(ipl)%m, sol_sumno3(j), sol_sumsolp(j), grazeop_db(igr)%eat, grazeop_db(igr)%manure
              !end if
              pcom(j)%dtbl(idtbl)%num_actions(iac) = pcom(j)%dtbl(idtbl)%num_actions(iac) + 1

          !initiate growing season for hru_lte
          case ("grow_init")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            hlt(j)%gro = "y"
            hlt(j)%g = 0.
            hlt(j)%alai = 0.
            hlt(j)%dm = 0.
            hlt(j)%hufh = 0.

          !end growing season for hru_lte
          case ("grow_end")
            !calculate yield - print lai, biomass and yield
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

              idp = hlt(j)%iplant
              if (hlt(j)%pet < 10.) then
                wur = 100.
              else
                wur = 100. * hlt(j)%aet / hlt(j)%pet
              endif
              hiad1 = (pldb(idp)%hvsti - pldb(idp)%wsyf) *                            &
                        (wur / (wur + Exp(6.13 - .0883 * wur))) + pldb(idp)%wsyf
              hiad1 = amin1 (hiad1, pldb(idp)%hvsti)
              yield = 0.8 * hlt(j)%dm * hiad1  ! * hlt(isd)%stress
              hlt(j)%yield = yield / 1000.
              hlt(j)%npp = hlt(j)%dm / 1000.
              hlt(j)%lai_mx = hlt(j)%alai
              !compute annual net primary productivity (npp) for perennial non-harvested?
              !use output.mgt print code
              !write() isd, time%day, time%yrc, pldb(iplt)%plantnm, hlt(isd)%alai, hlt(isd)%dm, yield
              hlt(j)%gro = "n"
              hlt(j)%g = 0.
              hlt(j)%alai = 0.
              hlt(j)%dm = 0.     !adjust for non-harvested perennials?
              hlt(j)%hufh = 0.
              hlt(j)%aet = 0.
              hlt(j)%pet = 0.

          !drainage water management
          case ("drain_control") !! set drain depth for drainage water management
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            if (pcom(j)%dtbl(idtbl)%num_actions(iac) <= Int(d_tbl%act(iac)%const2)) then
             istr = hru(j)%tiledrain
              hru(j)%lumv%sdr_dep = d_tbl%act(iac)%const * sdr(istr)%depth
              !if (hru(j)%lumv%sdr_dep > 0) then
              !  do jj = 1, soil(j)%nly
              !    if (hru(j)%lumv%sdr_dep < soil(j)%phys(jj)%d) hru(j)%lumv%ldrain = jj
              !    if (hru(j)%lumv%sdr_dep < soil(j)%phys(jj)%d) exit
              !  end do
              !else
              !  hru(j)%lumv%ldrain = 0
              !end if
              pcom(j)%dtbl(idtbl)%num_actions(iac) = pcom(j)%dtbl(idtbl)%num_actions(iac) + 1
            end if

          !tile flow control for saturated buffers
          case ("flow_control") !! set flow fractions to buffer tile and direct to channel
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur
            select case (d_tbl%act(iac)%option)
            case ("min_flo")
              if (hwb_d(j)%qtile < d_tbl%act(iac)%const) then
                frac = 1.
              else
                frac = d_tbl%act(iac)%const / hwb_d(j)%qtile
              end if
              ! set inflow hydrograph fraction of recieving objects - used for dtbl flow fractions
              ! set first object hyd fractin as defined in decision table
              inhyd = dtbl_flo(idtbl)%act(iac)%ob_num
              ihyd_in = ob(ob_num)%rcvob_inhyd(inhyd)
              iob_out = ob(ob_num)%obj_out(inhyd)
              ob(iob_out)%frac_in(ihyd_in) = frac

              ! set second hydrograph fraction
              if (inhyd < ob(ob_num)%src_tot .and. dtbl_flo(idtbl)%act(iac)%typ /= "irrigate_direct") then
                inhyd = inhyd + 1
                ihyd_in = ob(ob_num)%rcvob_inhyd(inhyd)
                iob_out = ob(ob_num)%obj_out(inhyd)
                ob(iob_out)%frac_in(ihyd_in) = 1. - frac
              end if

            case ("linear")

            case ("power")

            end select

          !tile flow control for saturated buffers
          case ("tile_control") !! set flow fractions to buffer tile and direct to channel
            icon = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur
            select case (d_tbl%act(iac)%option)
            case ("min_flo")
              if (hwb_d(j)%qtile < d_tbl%act(iac)%const) then
                frac = 1.
              else
                frac = d_tbl%act(iac)%const / hwb_d(j)%qtile
              end if
              ! set inflow hydrograph fraction of recieving objects - used for dtbl flow fractions
              ! set first object hyd fractin as defined in decision table
              inhyd = dtbl_flo(idtbl)%act(iac)%ob_num
              ihyd_in = ob(ob_num)%rcvob_inhyd(inhyd)
              iob_out = ob(ob_num)%obj_out(inhyd)
              ob(iob_out)%frac_in(ihyd_in) = frac

              ! set second hydrograph fraction
              if (inhyd < ob(ob_num)%src_tot .and. dtbl_flo(idtbl)%act(iac)%typ /= "irrigate_direct") then
                inhyd = inhyd + 1
                ihyd_in = ob(ob_num)%rcvob_inhyd(inhyd)
                iob_out = ob(ob_num)%obj_out(inhyd)
                ob(iob_out)%frac_in(ihyd_in) = 1. - frac
              end if

            case ("linear")

            case ("power")

            end select

          !land use change
          case ("lu_change")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur
            ilu = d_tbl%act_typ(iac)
            hru(j)%dbs%land_use_mgt = ilu
            hru(j)%land_use_mgt_c = d_tbl%act(iac)%file_pointer
            isol = hru(j)%dbs%soil
            call plant_init (1)

          !channel change
          case ("chan_change")
            ich = ob_cur
            !set new cover and name for calibration
            sd_ch(ich)%cov = d_tbl%act(iac)%const
            sd_ch(ich)%order = d_tbl%act(iac)%file_pointer

          ! burning
          case ("burn")
            j = d_tbl%act(iac)%ob_num
            if (j == 0) j = ob_cur

            if (pcom(j)%dtbl(idtbl)%num_actions(iac) <= Int(d_tbl%act(iac)%const2)) then
              iburn = d_tbl%act_typ(iac)           !burn type from fire data base
              do ipl = 1, pcom(j)%npl
                call pl_burnop (j, ipl, iburn)
              end do

              if (pco%mgtout == "y") then
                write (2612, *) j, time%yrc, time%mo, time%day, "        ", "    BURN", phubase(j),    &
                    pcom(j)%plcur(ipl)%phuacc, soil(j)%sw,pl_mass(j)%tot(ipl)%m, rsd1(j)%tot(ipl)%m,   &
                    sol_sumno3(j), sol_sumsolp(j)
              end if
              pcom(j)%dtbl(idtbl)%num_actions(iac) = pcom(j)%dtbl(idtbl)%num_actions(iac) + 1
            end if

          !herd management - move the herd
          case ("herd")

          !water rights decision to move water
          case ("water_rights")

          end select
        end if
      end do

      return
      end subroutine actions


text:
      SUBROUTINE slDMAV (RMAT, AXVEC)
*+
*     - - - - - -
*      D M A V
*     - - - - - -
*
*  From a rotation matrix, determine the corresponding axial vector.
*  (double precision)
*
*  A rotation matrix describes a rotation about some arbitrary axis,
*  called the Euler axis.  The "axial vector" returned by this routine
*  has the same direction as the Euler axis, and its magnitude is the
*  amount of rotation in radians.  (The magnitude and direction can be
*  separated by means of the routine slDVN.)
*
*  Given:
*    RMAT   d(3,3)   rotation matrix
*
*  Returned:
*    AXVEC  d(3)     axial vector (radians)
*
*  The reference frame rotates clockwise as seen looking along
*  the axial vector from the origin.
*
*  If RMAT is null, so is the result.
*
*  Last revision:   26 November 2005
*
*  Copyright P.T.Wallace.  All rights reserved.
*
*  License:
*    This program is free software; you can redistribute it and/or modify
*    it under the terms of the GNU General Public License as published by
*    the Free Software Foundation; either version 2 of the License, or
*    (at your option) any later version.
*
*    This program is distributed in the hope that it will be useful,
*    but WITHOUT ANY WARRANTY; without even the implied warranty of
*    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*    GNU General Public License for more details.
*
*    You should have received a copy of the GNU General Public License
*    along with this program (see SLA_CONDITIONS); if not, write to the
*    Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
*    Boston, MA  02110-1301  USA
*
*  Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
*-

      IMPLICIT NONE

      DOUBLE PRECISION RMAT(3,3),AXVEC(3)

      DOUBLE PRECISION X,Y,Z,S2,C2,PHI,F

      X = RMAT(2,3)-RMAT(3,2)
      Y = RMAT(3,1)-RMAT(1,3)
      Z = RMAT(1,2)-RMAT(2,1)
      S2 = SQRT(X*X+Y*Y+Z*Z)
      IF (S2.NE.0D0) THEN
         C2 = RMAT(1,1)+RMAT(2,2)+RMAT(3,3)-1D0
         PHI = ATAN2(S2,C2)
         F = PHI/S2
         AXVEC(1) = X*F
         AXVEC(2) = Y*F
         AXVEC(3) = Z*F
      ELSE
         AXVEC(1) = 0D0
         AXVEC(2) = 0D0
         AXVEC(3) = 0D0
      END IF

      END

ext: f

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
C
C
*+ synch_fit.f

       subroutine synch_fit(ndata,freq_data,flux_data,err_data,
     *                      iopt,res_data,s                   )
C      --------------------------------------------------------
C
C Fit a theoretical synchrotron spectrum with one of a number of models
C
C Input:
C   Number of data points
       integer     ndata
C   frequency data points
       real*4      freq_data(ndata)
C   flux data points
       real*4      flux_data(ndata)
C   error data points
       real*4      err_data(ndata)
C   option code
       integer     iopt
C Returned:
C   results array
       real*4      res_data(*)
C   status
       integer     s
C
C
C The supplied data are fitted with a theoretical synchrotron spectrum
C which can take one of a number of forms depending on the value of
C iopt.
C
C iopt      spectral form and minimized function
C --------------------------------------------------------------------
C
C 11        power-law of the form:  S * x**(-alpha)    minimizing SoS
C 110                   "                "             minimizing Chi2
C 12        power-law with thermal component
C                      S x**(-alpha) + T x**(-0.1)     minimizing SoS
C 120                   "                              minimizing Chi2
C 13        power-law with thermal component; alpha0 given
C                      S x**(-alpha0) + T x**(-0.1)    minimizing SoS
C 130                   "                              minimizing Chi2
C 14        power-law with thermal component; alpha constrained
C                      S x**(-alpha0+a) + T x**(-0.1)  minimizing SoS
C 140                   "                              minimizing Chi2
C
C 20        current synchrotron spectrum               minimizing SoS
C 200                   "                              minimizing Chi2
C 21        current synchrotron spectrum + T x**(-0.1) minimizing SoS
C 210                   "                              minimizing Chi2
C 22        current synchrotron spectrum + T x**(-0.1)
C                      with fixed break frequency      minimizing SoS
C 220                   "                              minimizing Chi2
C
C 31        as 21 except break >= max_frequency for
C                      flat spectra                    minimizing SoS
C 310                   "                              minimizing Chi2
C 32        as 21 except use power-laws for flat spectra
C                                                      minimizing SoS
C 320                   "                              minimizing Chi2
C 33        as 21 except use power-laws for flat spectra
C                      and break >= max_frequency for
C                      flat, but steepening spectrum   minimizing SoS
C 330                   "                              minimizing Chi2
*-

C define local data variables passed in common to the internal routine
       real*8                    x_data(20), y_data(20), e_data(20)
       common /local_synch_fit/  x_data,     y_data,     e_data
       real*8                    con(20)
       common /local_synch_con/  con

C variables used by NAG -- X holds results
       integer       lrwork, liwork
       parameter    (lrwork=300, liwork=100)
       real*8        x(4), r_work(lrwork), fsumsq
       integer       ifail, i_work(liwork)
       real*8        synch_bfun8, aa, xs

C parameters describing spectral type
       integer       spectral_type
       real*4        injection_index_E, alpha

C loop counter, output unit number and number of model parameters
       integer       n, nn, iout, nx, nm

C flags to indiacte a local change of mode
       integer       ioptl

C flag used in ordering of the data
       integer       io(20)
       real*4        fmin

C logical variable used to test for input data which could be fitted by
C the spectrum as defined
       logical       test_OK, test_steep

C local variable used as a flag when fitwas not perfect or data do not
C obey tests before fitting -- this is set only when an error state
C does NOT occur
       integer                    local_warn
       common  /local_synch_warn/ local_warn

       include '/mrao/anmap/include/synch_errors.inc'

C check status on entry and initialise local warning flag
       if (s.ne.0) return
       if (ndata.le.1) then
         local_warn = -1
         return
       else
         local_warn = 0
       endif
       ioptl = iopt

C prepare data; scale so that values are within a sensible range
       test_OK = .true.
       do n=1,ndata
         if (freq_data(n).le.0.0 .or. flux_data(n).le.0.0) then
           test_OK = .false.
         endif
         if (iopt.gt.100) then
           if (err_data(n).le.0.0) then
             test_OK = .false.
           endif
         endif
         io(n) = n
       enddo
       if (.not.test_OK) then
         local_warn = -1
         return
       end if

C order and scale data
       do n=1,ndata
         fmin = 1.0E+30
         do nn=1,ndata
           if (io(nn).ne.0) then
             if (freq_data(nn).lt.fmin) then
               fmin = freq_data(nn)
               nx = nn
             endif
           endif
         enddo
         io(nx) = 0
         x_data(n) = freq_data(nx)
         y_data(n) = flux_data(nx)
         e_data(n) = err_data(nx)
         if (n.eq.1) then
           nm = nx
         endif
       enddo
       if (iopt.gt.100) then
         do n=1,ndata
           e_data(n) = e_data(n)/y_data(1)
         enddo
       endif
       do n=ndata,1,-1
         x_data(n) = x_data(n)/x_data(1)
         y_data(n) = y_data(n)/y_data(1)
       end do
       if (.not.test_ok) then
         local_warn = -1
         return
       end if

C find spectral range
       call synch_enqsp(injection_index_E,spectral_type)

C look for standard steep spectrum -- useful for all modes
       test_OK = .false.
       do n=1,ndata-1
         do nn=n+1,ndata
           alpha = - log(y_data(nn)/y_data(n))/
     *               log(x_data(nn)/x_data(n))
           test_OK = test_OK .or. (2.0*alpha+1.0).gt.injection_index_E
         end do
       end do
C look for steepening spectrum
       test_steep = .true.
       aa = 0.0
       do n=1,ndata-1
         nn = n + 1
         alpha = - log(y_data(nn)/y_data(n))/
     *             log(x_data(nn)/x_data(n))
         test_steep = test_steep .and. alpha.gt.aa
         aa = alpha
       enddo

C setup trial spectral index
       aa = -log(flux_data(2))/log(freq_data(2))

C initialise for call to NAG fitting routine
       if (iopt.eq.11 .or. iopt.eq.110) then
         nx = 2
         x(1) = y_data(1)
         x(2) = aa
       elseif (iopt.eq.12 .or. iopt.eq.120) then
         nx = 3
         if (ndata.lt.3) then
           local_warn = -1
           return
         endif
         if (test_OK) then
           x(1) = y_data(1)
           x(3) = aa
           x(2) = 0.0
         else
           x(1) = y_data(1)
           x(3) = aa
           x(2) = y_data(ndata)
         endif
       elseif (iopt.eq.13 .or. iopt.eq.130) then
         nx = 2
         if (test_OK) then
           x(1) = y_data(1)
           x(2) = 0.0
         else
           x(1) = y_data(1)
           x(2) = y_data(ndata)
         endif
         con(1) = res_data(3)
       elseif (iopt.eq.14 .or. iopt.eq.140) then
         nx = 3
         if (ndata.lt.3) then
           local_warn = -1
           return
         endif
         if (test_OK) then
           x(1) = y_data(1)
           x(3) = aa
           x(2) = 0.0
         else
           x(1) = y_data(1)
           x(3) = aa
           x(2) = y_data(ndata)
         endif
         con(1) = res_data(3)
       elseif (iopt.eq.20 .or. iopt.eq.200) then
         call synch_enqlim( con(1) )
         con(1) = con(1)*0.8D+0
         nx = 2
         if (test_OK) then
           x(2) = x_data(ndata)
           xs = x_data(1)/x(2)
           x(1) = y_data(1)*synch_bfun8(1.0D+0,s)/
     *                      synch_bfun8(xs,s)
         else
           local_warn = 100
           return
         endif
       elseif (iopt.eq.21 .or. iopt.eq.210) then
         nx = 3
         call synch_enqlim( con(1) )
         con(1) = con(1) * 0.8D+0
         if (test_OK) then
           x(3) = x_data(ndata)
           xs = x_data(1)/x(3)
           x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                      synch_bfun8(xs,s)
           x(2) = 0.0D+0
         else
           x(3) = 2.0D+0*x_data(ndata)
           xs = x_data(1)/x(3)
           x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                          synch_bfun8(xs,s)
           x(2) = y_data(ndata)
         endif
       elseif (iopt.eq.22 .or. iopt.eq.220) then
         nx = 2
         call synch_enqlim( con(1) )
         con(1) = res_data(3)/freq_data(nm)
         if (test_OK) then
           xs = x_data(1)/con(1)
           x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                      synch_bfun8(xs,s)
           x(2) = 0.0D+0
         else
           xs = x_data(1)/con(1)
           x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                          synch_bfun8(xs,s)
           x(2) = y_data(ndata)
         endif

       elseif (iopt.eq.31 .or. iopt.eq.310) then
         nx = 3
         ioptl = 21
         call synch_enqlim( con(1) )
         con(1) = con(1) * 0.8D+0
         if (test_OK) then
           x(3) = x_data(ndata)
           xs = x_data(1)/x(3)
           x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                      synch_bfun8(xs,s)
           x(2) = 0.0D+0
         else
           con(1) = 1.0D+0
           x(3) = 2.0D+0*x_data(ndata)
           xs = x_data(1)/x(3)
           x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                          synch_bfun8(xs,s)
           x(2) = y_data(ndata)
         endif

       elseif (iopt.eq.32 .or. iopt.eq.320) then
         nx = 3
         call synch_enqlim( con(1) )
         con(1) = con(1) * 0.8D+0
         if (test_OK) then
           ioptl = 21
           x(3) = x_data(ndata)
           xs = x_data(1)/x(3)
           x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                      synch_bfun8(xs,s)
           x(2) = 0.0D+0
         else
           ioptl = 22
           con(1) = 100.0*x_data(ndata)
           xs = x_data(1)/con(1)
           x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                          synch_bfun8(xs,s)
           x(2) = y_data(ndata)
         endif

       elseif (iopt.eq.33 .or. iopt.eq.330) then
         nx = 3
         call synch_enqlim( con(1) )
         con(1) = con(1) * 0.8D+0
         if (test_OK) then
           ioptl = 21
           x(3) = x_data(ndata)
           xs = x_data(1)/x(3)
           x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                      synch_bfun8(xs,s)
           x(2) = 0.0D+0
         else
           if (test_steep) then
             ioptl = 21
             con(1) = 1.0D+0
             x(3) = 2.0D+0*x_data(ndata)
             xs = x_data(1)/x(3)
             x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                            synch_bfun8(xs,s)
             x(2) = y_data(ndata)
           else
             ioptl = 22
             con(1) = 100.0*x_data(ndata)
             xs = x_data(1)/con(1)
             x(1) = y_data(1)*synch_bfun8(1.0D+0,s) /
     *                            synch_bfun8(xs,s)
             x(2) = y_data(ndata)
           endif
         endif

       end if

C initialise fitting mode and do the fit
       call synch_fitmode( ioptl )
       ifail = 1
       call e04fdf(ndata,nx,x,fsumsq,i_work,
     *             liwork,r_work,lrwork,ifail)

C check IFAIL on exit
       local_warn = 0
       if (ifail.eq.1) then
         call io_enqout(iout)
         call io_wrout('***(SYNCH_FIT) Error calling e04fdf')
         write(iout,*) '***(SYNCH_FIT) non-zero IFAIL = ',ifail
         s = ill_synNAG
         return
       else if (ifail.eq.2) then
         local_warn = 3000
       else if (ifail.eq.3) then
         local_warn = 2000
       else if (ifail.ge.4) then
         local_warn = 1000
       end if
       do n=1,3
         res_data(n) = abs(x(n))
       enddo
       res_data(4) = sqrt(fsumsq)
       res_data(1) = res_data(1) * flux_data(nm)

       if (ioptl.eq.11 .or. ioptl.eq.110) then
         res_data(1) = res_data(1)*(freq_data(nm)**res_data(2))
         res_data(2) = x(2)
         res_data(3) = 0.0
       endif
       if (ioptl.eq.13 .or. ioptl.eq.130 ) then
         res_data(3) = con(1)
       endif
       if (ioptl.eq.14 .or. ioptl.eq.140 ) then
         res_data(3) = con(1) + res_data(3)
       endif
       if (ioptl.eq.12 .or. ioptl.eq.120 .or.
     *     ioptl.eq.13 .or. ioptl.eq.130  .or.
     *     ioptl.eq.14 .or. ioptl.eq.140 ) then
         res_data(1) = res_data(1)*(freq_data(nm)**res_data(3))
         res_data(2) = res_data(2)*flux_data(nm)*(freq_data(nm)**0.1)
       endif
       if (ioptl.eq.22 .or. ioptl.eq.220 ) then
         res_data(3) = con(1)
       endif
       if (ioptl.eq.21 .or. ioptl.eq.210 .or.
     *     ioptl.eq.22 .or. ioptl.eq.220 ) then
         res_data(2) = res_data(2)*flux_data(nm)*(freq_data(nm)**0.1)
         res_data(3) = res_data(3)*freq_data(nm)
       endif

       if (ioptl.eq.20 .or. ioptl.eq.200) then
         res_data(2) = res_data(2)*freq_data(nm)
       endif
C check error on exit
       call iocmd_err(s,'synch_fit','Failed')
       end
C
C
*+ lsfun1
C
       subroutine lsfun1(ndata,nres,x,fvecc)
C      -------------------------------------
C
C Local routine used by NAG for fitting purposes
C
       integer     ndata, nres, n, status
       real*8      fvecc(ndata), x(nres)
       real*8      x_data(20), y_data(20), e_data(20)
       real*8      x_point, y_point
       real*8      synch_bfun8
       integer     mode, icount
       common /local_synch_fit/   x_data, y_data, e_data
       common /local_lsfun1_mode/ mode, icount
       real*8                     con(20)
       common /local_synch_con/   con

       status = 0
       if (mode.eq.11 .or. mode.eq.110) then
         do n=1,ndata
           fvecc(n) =
     *        abs(x(1)) * (x_data(n)**(-x(2))) -
     *        y_data(n)
           if (mode.eq.110) then
             fvecc(n) = fvecc(n) / e_data(n)
           endif
         end do
         icount = icount + 1

       else if (mode.eq.12 .or. mode.eq.120) then
         do n=1,ndata
           fvecc(n) =
     *        abs(x(1)) * (x_data(n)**(-abs(x(3)))) +
     *        abs(x(2))*(x_data(n)**(-0.1D+0)) -
     *        y_data(n)
           if (mode.eq.120) then
             fvecc(n) = fvecc(n) / e_data(n)
           endif
         end do
         icount = icount + 1

       else if (mode.eq.13 .or. mode.eq.130) then
         do n=1,ndata
           fvecc(n) =
     *        abs(x(1)) * (x_data(n)**(-con(1))) +
     *        abs(x(2))*(x_data(n)**(-0.1D+0)) -
     *        y_data(n)
           if (mode.eq.130) then
             fvecc(n) = fvecc(n) / e_data(n)
           endif
         end do
         icount = icount + 1

       else if (mode.eq.14 .or. mode.eq.140) then
         do n=1,ndata
           fvecc(n) =
     *        abs(x(1)) * (x_data(n)**(-(con(1)+abs(x(3))))) +
     *        abs(x(2))*(x_data(n)**(-0.1D+0)) -
     *        y_data(n)
           if (mode.eq.140) then
             fvecc(n) = fvecc(n) / e_data(n)
           endif
         end do
         icount = icount + 1

       else if (mode.eq.20 .or. mode.eq.200) then
         do n=1,ndata
           x_point = x_data(n) / max(abs(x(2)),x_data(ndata)/con(1))
           y_point = abs(x(1)) * synch_bfun8(x_point,status) /
     *                           synch_bfun8(1.0D+0,status)
           fvecc(n)= log(y_point) - log(y_data(n))
           if (mode.eq.200) then
             fvecc(n) = fvecc(n) * y_data(n) / e_data(n)
           endif
         end do
         icount = icount + 1

       else if (mode.eq.21 .or. mode.eq.210) then
         do n=1,ndata
           x_point = x_data(n) / max(abs(x(3)),x_data(ndata)/con(1))
           y_point = abs(x(1)) * synch_bfun8(x_point,status) /
     *                           synch_bfun8(1.0D+0,status) +
     *               abs(x(2))*(x_data(n)**(-0.1D+0))
           fvecc(n)= log(y_point) - log(y_data(n))
           if (mode.eq.210) then
             fvecc(n) = fvecc(n) * y_data(n) / e_data(n)
           endif
         end do
         icount = icount + 1

       else if (mode.eq.22 .or. mode.eq.220) then
         do n=1,ndata
           x_point = x_data(n) / con(1)
           y_point = abs(x(1)) * synch_bfun8(x_point,status) /
     *                           synch_bfun8(1.0D+0,status) +
     *               abs(x(2))*(x_data(n)**(-0.1D+0))
           fvecc(n)= log(y_point) - log(y_data(n))
           if (mode.eq.220) then
             fvecc(n) = fvecc(n) * y_data(n) / e_data(n)
           endif
         end do
         icount = icount + 1
       end if
       call iocmd_err(status,'LSFUN1 -- synch_fit',' ')
       end
C
C
*+
       real function synch_fitfun( x, iopt, res, s )
C      ---------------------------------------------
C
C Return the value of the fitted synchrotron spectrum at x
C
C Given:
C   frequency
       real*4        x
C   spectral type
       integer       iopt
C   results data
       real*4        res(*)
C Updated:
C   error status
       integer       s
C
C-
       real*4      synch_bfun

       if (s.ne.0) return

       if (iopt.eq.11 .or. iopt.eq.110) then
         synch_fitfun = res(1)*(x**(-res(2)))

       elseif (iopt.eq.12 .or. iopt.eq.120 .or.
     *         iopt.eq.13 .or. iopt.eq.130 .or.
     *         iopt.eq.14 .or. iopt.eq.140 ) then
         synch_fitfun = res(1)*(x**(-res(3))) +
     *                  res(2)*(x**(-0.1))

       elseif (iopt.eq.20 .or. iopt.eq.200) then
         s = 1
         synch_fitfun = res(1)*synch_bfun(x/res(2),s)/
     *                         synch_bfun(1.0,s)

       elseif (iopt.eq.21 .or. iopt.eq.210 .or.
     *         iopt.eq.22 .or. iopt.eq.220 .or.
     *         iopt.eq.31 .or. iopt.eq.310 .or.
     *         iopt.eq.32 .or. iopt.eq.320 .or.
     *         iopt.eq.33 .or. iopt.eq.330 .or.
     *         iopt.eq.34 .or. iopt.eq.340 ) then
         s = 1
         synch_fitfun = res(1)*synch_bfun(x/res(3),s)/
     *                         synch_bfun(1.0,s) +
     *                  res(2)*(x**(-0.1))

       endif
       end

max_stars_repo_name: CavendishAstrophysics/anmap

max_issues_repo_name: CavendishAstrophysics/anmap

max_forks_repo_name: CavendishAstrophysics/anmap


text:
!/ ------------------------------------------------------------------- /
      MODULE W3IDATMD
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         29-May-2009 |
!/                  +-----------------------------------+
!/
!/    02-Apr-2004 : Origination.                        ( version 3.06 )
!/    19-Jul-2006 : Adding auxiliary grids.             ( version 3.10 )
!/    04-Oct-2006 : Add filter to array pointers.       ( version 3.10 )
!/    29-May-2009 : Preparing distribution version.     ( version 3.14 )
!/
!/    Copyright 2009 National Weather Service (NWS),
!/       National Oceanic and Atmospheric Administration.  All rights
!/       reserved.  WAVEWATCH III is a trademark of the NWS.
!/       No unauthorized use without permission.
!/
!  1. Purpose :
!
!     Define data structures to set up wave model input data for
!     several models simultaneously.
!
!  2. Variables and types :
!
!      Name      Type  Scope    Description
!     ----------------------------------------------------------------
!      NIDATA    Int.  Public   Number of models in array dim.
!      IIDATA    Int.  Public   Selected model for output, init. at -1.
!      INPUT     TYPE  Public   Basic data structure.
!      INPUTS    INPUT Public   Array of data structures.
!     ----------------------------------------------------------------
!
!     All elements of INPUT are aliased to pointers with the same
!     name. Some aditional pointer provide previous equivalenced
!     parameters. These pointers are defined as :
!
!      Name      Type  Scope    Description
!     ----------------------------------------------------------------
!      TLN       I.A.  Public   Time for water level field.
!      TC0/N     I.A.  Public   Times for current fields.
!      TW0/N     I.A.  Public   Times for wind fields.
!      TIN       I.A.  Public   Time for ice field.
!      TnN       I.A.  Public   Time for data types 1-3.
!      TDN       I.A.  Public   Time for next data.
!      TG0/N     I.A.  Public   Times for grid motion data.
!      TFN       I.A.  Public   Array consolidating most above times.
!      GA0/N     Real  Public   Norm of grid speed vector.
!      GD0/N     Real  Public   Direction of grid speed vector.
!      WX0/N     R.A.  Public   Cartesian X and Y wind components
!      WY0/N     R.A.  Public      for both times.
!      DT0/N     R.A.  Public   Corr. air-sea temperature differences.
!      CX0/N     R.A.  Public   Cartesian X and Y current components
!      CY0/N     R.A.  Public      for both times.
!      WLEV      R.A.  Public   Next water level field.
!      ICEI      R.A.  Public   Ice concentrations.
!      IINIT     Log.  Public   Flag for array initialization.
!      HML       R.A.  Public   Mixed layer depth, QL, 150525
!      FLLEV     Log.  Public   Flag for water level input.
!      FLCUR     Log.  Public   Flag for current input.
!      FLWIND    Log.  Public   Flag for wind input.
!      FLICE     Log.  Public   Flag for ice input.
!      FLHML     Log.  Public   Flag for mixed layer depth input.
!                               QL, 150525
!      FLAGS     L.A.  Public   Array consolidating the above four
!                               flags, as well as four additional
!                               data flags.
!     ----------------------------------------------------------------
!
!  3. Subroutines and functions :
!
!      Name      Type  Scope    Description
!     ----------------------------------------------------------------
!      W3NINP    Subr. Public   Set number of grids/models.
!      W3DIMI    Subr. Public   Set dimensions of arrays.
!      W3SETI    Subr. Public   Point to selected grid / model.
!     ----------------------------------------------------------------
!
!  4. Subroutines and functions used :
!
!      Name      Type  Module   Description
!     ----------------------------------------------------------------
!      W3SETG    Subr. W3GDATMD Point to proper model grid.
!      STRACE    Subr. W3SERVMD Subroutine tracing.
!      EXTCDE    Subr. W3SERVMD Abort program with exit code.
!     ----------------------------------------------------------------
!
!  5. Remarks :
!
!     - The number of grids is taken from W3GDATMD, and needs to be
!       set first with W3DIMG.
!
!  6. Switches :
!
!     !/S    Enable subroutine tracing.
!     !/T    Enable test output
!
!  7. Source code :
!
!/ ------------------------------------------------------------------- /
      PUBLIC
!/
!/ Conventional declarations
!/
      INTEGER                 :: NIDATA = -1, IIDATA = -1
!/
!/ Data structure INPUT
!/
      TYPE INPUT
        INTEGER               :: TFN(2,8), TC0(2), TW0(2), TDN(2),    &
                                 TG0(2)
        REAL                  :: GA0, GD0, GAN, GDN
        ! QL, 150525, HML
        REAL, POINTER         :: WX0(:,:), WY0(:,:), DT0(:,:),        &
                                 WXN(:,:), WYN(:,:), DTN(:,:),        &
                                 CX0(:,:), CY0(:,:), CXN(:,:),        &
                                 CYN(:,:), WLEV(:,:), ICEI(:,:),      &
                                 HML(:,:)
        LOGICAL               :: IINIT
        LOGICAL               :: FLAGS(8)
      END TYPE INPUT
!/
!/ Data storage
!/
      TYPE(INPUT), TARGET, ALLOCATABLE :: INPUTS(:)
!/
!/ Data aliasses for structure INPUT(S)
!/
      INTEGER, POINTER        :: TFN(:,:), TLN(:), TC0(:), TCN(:),    &
                                 TW0(:), TWN(:), TIN(:), T0N(:),      &
                                 T1N(:), T2N(:), TDN(:), TG0(:), TGN(:)
      REAL, POINTER           :: GA0, GD0, GAN, GDN
      ! QL, 150525, HML
      REAL, POINTER           :: WX0(:,:), WY0(:,:), DT0(:,:),        &
                                 WXN(:,:), WYN(:,:), DTN(:,:),        &
                                 CX0(:,:), CY0(:,:), CXN(:,:),        &
                                 CYN(:,:), WLEV(:,:), ICEI(:,:),      &
                                 HML(:,:)
      LOGICAL, POINTER        :: IINIT
      LOGICAL, POINTER        :: FLAGS(:)
      ! QL, 150525, FLHML
      LOGICAL, POINTER        :: FLLEV, FLCUR, FLWIND, FLICE, FLHML
!/
      CONTAINS
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3NINP ( NDSE, NDST )
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         19-Jul-2006 !
!/                  +-----------------------------------+
!/
!/    02-Apr-2004 : Origination.                        ( version 3.06 )
!/    19-Jul-2006 : Adding auxiliary grids.             ( version 3.10 )
!/
!  1. Purpose :
!
!     Set up the number of grids to be used.
!
!  2. Method :
!
!     Use data stored in NGRIDS in W3GDATMD.
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!       NDSE    Int.   I   Error output unit number.
!       NDST    Int.   I   Test output unit number.
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!     See module documentation.
!
!  5. Called by :
!
!     Any program that uses this grid structure.
!
!  6. Error messages :
!
!     - Error checks on previous setting of variable NGRIDS.
!
!  7. Remarks :
!
!  8. Structure :
!
!  9. Switches :
!
!     !/S    Enable subroutine tracing.
!     !/T    Enable test output
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE W3GDATMD, ONLY: NGRIDS, NAUXGR
      USE W3SERVMD, ONLY: EXTCDE
!
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      INTEGER, INTENT(IN)     :: NDSE, NDST
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: I
!/
!
! -------------------------------------------------------------------- /
! 1.  Test input and module status
!
      IF ( NGRIDS .EQ. -1 ) THEN
          WRITE (NDSE,1001) NGRIDS
          CALL EXTCDE (1)
        END IF
!
! -------------------------------------------------------------------- /
! 2.  Set variable and allocate arrays
!
      ALLOCATE ( INPUTS(-NAUXGR:NGRIDS) )
      NIDATA = NGRIDS
!
! -------------------------------------------------------------------- /
! 3.  Initialize parameters
!
      DO I=-NAUXGR, NGRIDS
        INPUTS(I)%TFN(1,:) = -1
        INPUTS(I)%TFN(2,:) =  0
        INPUTS(I)%TC0(1)   = -1
        INPUTS(I)%TC0(2)   =  0
        INPUTS(I)%TW0(1)   = -1
        INPUTS(I)%TW0(2)   =  0
        INPUTS(I)%TDN(1)   = -1
        INPUTS(I)%TDN(2)   =  0
        INPUTS(I)%TG0(1)   = -1
        INPUTS(I)%TG0(2)   =  0
        INPUTS(I)%IINIT    = .FALSE.
        INPUTS(I)%FLAGS    = .FALSE.
        END DO
!
      RETURN
!
! Formats
!
 1001 FORMAT (/' *** ERROR W3NINP : NGRIDS NOT YET SET *** '/         &
               '                    NGRIDS = ',I10/                   &
               '                    RUN W3NMOD FIRST'/)
!
!/
!/ End of W3NINP ----------------------------------------------------- /
!/
      END SUBROUTINE W3NINP
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3DIMI  ( IMOD, NDSE, NDST )
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         04-Oct-2006 !
!/                  +-----------------------------------+
!/
!/    02-Apr-2004 : Origination.                        ( version 3.06 )
!/    19-Jul-2006 : Adding auxiliary grids.             ( version 3.10 )
!/    04-Oct-2006 : Add filter to array pointers.       ( version 3.10 )
!/
!  1. Purpose :
!
!     Initialize an individual data grid at the proper dimensions.
!
!  2. Method :
!
!     Allocate directly into the structure array. Note that
!     this cannot be done through the pointer alias!
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!       IMOD    Int.   I   Model number to point to.
!       NDSE    Int.   I   Error output unit number.
!       NDST    Int.   I   Test output unit number.
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!     See module documentation.
!
!  5. Called by :
!
!     Main wave model drivers.
!
!  6. Error messages :
!
!     - Check on input parameters.
!     - Check on previous allocation.
!
!  7. Remarks :
!
!     - W3SETI needs to be called after allocation to point to
!       proper allocated arrays.
!
!  8. Structure :
!
!     See source code.
!
!  9. Switches :
!
!     !/S    Enable subroutine tracing.
!     !/T    Enable test output
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE W3GDATMD, ONLY: NGRIDS, NAUXGR, IGRID, W3SETG, NX, NY
      USE W3SERVMD, ONLY: EXTCDE
!
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      INTEGER, INTENT(IN)     :: IMOD, NDSE, NDST
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
      INTEGER                 :: JGRID
!/
!
! -------------------------------------------------------------------- /
! 1.  Test input and module status
!
      IF ( NGRIDS .EQ. -1 ) THEN
          WRITE (NDSE,1001)
          CALL EXTCDE (1)
        END IF
!
      IF ( IMOD.LT.-NAUXGR .OR. IMOD.GT.NIDATA ) THEN
          WRITE (NDSE,1002) IMOD, -NAUXGR, NIDATA
          CALL EXTCDE (2)
        END IF
!
      IF ( INPUTS(IMOD)%IINIT ) THEN
          WRITE (NDSE,1003)
          CALL EXTCDE (3)
        END IF
!
      JGRID  = IGRID
      IF ( JGRID .NE. IMOD ) CALL W3SETG ( IMOD, NDSE, NDST )
!
! -------------------------------------------------------------------- /
! 2.  Allocate arrays
!
      FLLEV  => INPUTS(IMOD)%FLAGS(1)
      FLCUR  => INPUTS(IMOD)%FLAGS(2)
      FLWIND => INPUTS(IMOD)%FLAGS(3)
      FLICE  => INPUTS(IMOD)%FLAGS(4)
      ! QL, 150525, from coupler, set in wav_comp_mct.F90
      FLHML  => INPUTS(IMOD)%FLAGS(5)
!
      IF ( FLLEV  ) ALLOCATE ( INPUTS(IMOD)%WLEV(NX,NY) )
!
      IF ( FLCUR  ) ALLOCATE ( INPUTS(IMOD)%CX0(NX,NY) ,              &
                               INPUTS(IMOD)%CY0(NX,NY) ,              &
                               INPUTS(IMOD)%CXN(NX,NY) ,              &
                               INPUTS(IMOD)%CYN(NX,NY) )
!
      IF ( FLWIND ) ALLOCATE ( INPUTS(IMOD)%WX0(NX,NY) ,              &
                               INPUTS(IMOD)%WY0(NX,NY) ,              &
                               INPUTS(IMOD)%DT0(NX,NY) ,              &
                               INPUTS(IMOD)%WXN(NX,NY) ,              &
                               INPUTS(IMOD)%WYN(NX,NY) ,              &
                               INPUTS(IMOD)%DTN(NX,NY) )
      IF ( FLWIND ) THEN
          INPUTS(IMOD)%DT0 = 0.
          INPUTS(IMOD)%DTN = 0.
      END IF
!
      IF ( FLICE  ) ALLOCATE ( INPUTS(IMOD)%ICEI(NX,NY) )
! QL, 150525, HML from coupler
      IF ( FLHML  ) ALLOCATE ( INPUTS(IMOD)%HML(NX,NY) )
!
      INPUTS(IMOD)%IINIT  = .TRUE.
!
! -------------------------------------------------------------------- /
! 3.  Point to allocated arrays
!
      CALL W3SETI ( IMOD, NDSE, NDST )
!
! -------------------------------------------------------------------- /
! 4.  Update counters in grid
!
! -------------------------------------------------------------------- /
! 5.  Restore previous grid setting if necessary
!
      IF ( JGRID .NE. IMOD ) CALL W3SETG ( JGRID, NDSE, NDST )
!
      RETURN
!
! Formats
!
 1001 FORMAT (/' *** ERROR W3DIMI : GRIDS NOT INITIALIZED *** '/      &
               '                    RUN W3NMOD FIRST '/)
 1002 FORMAT (/' *** ERROR W3DIMI : ILLEGAL MODEL NUMBER *** '/       &
               '                    IMOD   = ',I10/                   &
               '                    NAUXGR = ',I10/                   &
               '                    NIDATA = ',I10/)
 1003 FORMAT (/' *** ERROR W3DIMI : ARRAY(S) ALREADY ALLOCATED *** ')
!
!/
!/ End of W3DIMI ----------------------------------------------------- /
!/
      END SUBROUTINE W3DIMI
!/ ------------------------------------------------------------------- /
      SUBROUTINE W3SETI ( IMOD, NDSE, NDST )
!/
!/                  +-----------------------------------+
!/                  | WAVEWATCH III           NOAA/NCEP |
!/                  |           H. L. Tolman            |
!/                  |                        FORTRAN 90 |
!/                  | Last update :         04-Oct-2006 !
!/                  +-----------------------------------+
!/
!/    02-Apr-2004 : Origination.                        ( version 3.06 )
!/    19-Jul-2006 : Adding auxiliary grids.             ( version 3.10 )
!/    04-Oct-2006 : Add filter to array pointers.       ( version 3.10 )
!/
!  1. Purpose :
!
!     Select one of the WAVEWATCH III grids / models.
!
!  2. Method :
!
!     Point pointers to the proper variables in the proper element of
!     the GRIDS array.
!
!  3. Parameters :
!
!     Parameter list
!     ----------------------------------------------------------------
!       IMOD    Int.   I   Model number to point to.
!       NDSE    Int.   I   Error output unit number.
!       NDST    Int.   I   Test output unit number.
!     ----------------------------------------------------------------
!
!  4. Subroutines used :
!
!     See module documentation.
!
!  5. Called by :
!
!     Any subroutine.
!
!  6. Error messages :
!
!     Many subroutines in the WAVEWATCH system.
!
!  7. Remarks :
!
!  8. Structure :
!
!  9. Switches :
!
!     !/S    Enable subroutine tracing.
!     !/T    Enable test output
!
! 10. Source code :
!
!/ ------------------------------------------------------------------- /
      USE W3GDATMD, ONLY: NAUXGR
!
      USE W3SERVMD, ONLY: EXTCDE
!
      IMPLICIT NONE
!/
!/ ------------------------------------------------------------------- /
!/ Parameter list
!/
      INTEGER, INTENT(IN)     :: IMOD, NDSE, NDST
!/
!/ ------------------------------------------------------------------- /
!/ Local parameters
!/
!/
!
! -------------------------------------------------------------------- /
! 1.  Test input and module status
!
      IF ( NIDATA .EQ. -1 ) THEN
          WRITE (NDSE,1001)
          CALL EXTCDE (1)
        END IF
!
      IF ( IMOD.LT.-NAUXGR .OR. IMOD.GT.NIDATA ) THEN
          WRITE (NDSE,1002) IMOD, -NAUXGR, NIDATA
          CALL EXTCDE (2)
        END IF
!
! -------------------------------------------------------------------- /
! 2.  Set model numbers
!
      IIDATA = IMOD
!
! -------------------------------------------------------------------- /
! 3.  Set pointers
!
      TFN    => INPUTS(IMOD)%TFN
      TC0    => INPUTS(IMOD)%TC0
      TW0    => INPUTS(IMOD)%TW0
      TG0    => INPUTS(IMOD)%TG0
      TDN    => INPUTS(IMOD)%TDN
!
      TLN    => INPUTS(IMOD)%TFN(:,1)
      TCN    => INPUTS(IMOD)%TFN(:,2)
      TWN    => INPUTS(IMOD)%TFN(:,3)
      TIN    => INPUTS(IMOD)%TFN(:,4)
      T0N    => INPUTS(IMOD)%TFN(:,5)
      T1N    => INPUTS(IMOD)%TFN(:,6)
      T2N    => INPUTS(IMOD)%TFN(:,7)
      TGN    => INPUTS(IMOD)%TFN(:,8)
!
      GA0    => INPUTS(IMOD)%GA0
      GD0    => INPUTS(IMOD)%GD0
      GAN    => INPUTS(IMOD)%GAN
      GDN    => INPUTS(IMOD)%GDN
!
      IINIT  => INPUTS(IMOD)%IINIT
      FLAGS  => INPUTS(IMOD)%FLAGS
!
      FLLEV  => INPUTS(IMOD)%FLAGS(1)
      FLCUR  => INPUTS(IMOD)%FLAGS(2)
      FLWIND => INPUTS(IMOD)%FLAGS(3)
      FLICE  => INPUTS(IMOD)%FLAGS(4)
      ! QL, 150525, HML
      FLHML  => INPUTS(IMOD)%FLAGS(5)
!
      IF ( IINIT ) THEN
!
          IF ( FLLEV  ) THEN
              WLEV   => INPUTS(IMOD)%WLEV
            END IF
!
          IF ( FLCUR  ) THEN
              CX0    => INPUTS(IMOD)%CX0
              CY0    => INPUTS(IMOD)%CY0
              CXN    => INPUTS(IMOD)%CXN
              CYN    => INPUTS(IMOD)%CYN
            END IF
!
          IF ( FLWIND  ) THEN
              WX0    => INPUTS(IMOD)%WX0
              WY0    => INPUTS(IMOD)%WY0
              DT0    => INPUTS(IMOD)%DT0
              WXN    => INPUTS(IMOD)%WXN
              WYN    => INPUTS(IMOD)%WYN
              DTN    => INPUTS(IMOD)%DTN
            END IF
!
          IF ( FLICE  ) THEN
              ICEI   => INPUTS(IMOD)%ICEI
            END IF
! QL, 150525, HML
          IF ( FLHML  ) THEN
              HML    => INPUTS(IMOD)%HML
            END IF
!
        END IF
!
      RETURN
!
! Formats
!
 1001 FORMAT (/' *** ERROR W3SETI : GRIDS NOT INITIALIZED *** '/      &
               '                    RUN W3NMOD FIRST '/)
 1002 FORMAT (/' *** ERROR W3SETI : ILLEGAL MODEL NUMBER *** '/       &
               '                    IMOD   = ',I10/                   &
               '                    NAUXGR = ',I10/                   &
               '                    NIDATA = ',I10/)
!
!/
!/ End of W3SETI ----------------------------------------------------- /
!/
      END SUBROUTINE W3SETI
!/
!/ End of module W3IDATMD -------------------------------------------- /
!/
      END MODULE W3IDATMD

max_stars_repo_name: cbegeman/OceanNGD

max_stars_repo_licenses: zlib-acknowledgement RSA-MD FTL

max_issues_repo_name: cbegeman/OceanNGD

max_issues_repo_licenses: zlib-acknowledgement RSA-MD FTL

max_forks_repo_name: cbegeman/OceanNGD

max_forks_repo_licenses: zlib-acknowledgement RSA-MD FTL


text:
      MODULE initsn_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  10:47:04  03/09/06
      SUBROUTINE initsn (INDEPS, DIRSM, DIRSMH)
      USE vast_kind_param,ONLY: DOUBLE
      use cosmo_C, only : nppa
      integer :: INDEPS
      real(DOUBLE), DIMENSION(3,NPPA) :: DIRSM
      real(DOUBLE), DIMENSION(3,NPPA/3) :: DIRSMH
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


program hello
  print *, 'Hell, world!'
end program hello


text:
      SUBROUTINE BMG2_SymStd_relax_lines_x( &
     &                K, SO, QF, Q, SOR, B,&
     &                II, JJ, iGs, jGs,&
     &                NOG, NStncl, IRELAX_SYM, UPDOWN,&
     &                DATADIST, RWORK, NMSGr,&
     &                MPICOMM, XLINECOMM, YLINECOMM, halof&
     &                ) BIND(C,NAME='MPI_BMG2_SymStd_relax_lines_x')

! ======================================================================
!  --------------------
!   DESCRIPTION:
!  --------------------
!
!     Perform zebra-line relaxation in x.
!
! ======================================================================
! $license_flag$
! ======================================================================
!  --------------------
!   INPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   INPUT/OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   OUTPUT:
!  --------------------
!
!
! ======================================================================
!  --------------------
!   LOCAL:
!  --------------------
!
!
! ======================================================================
      USE ModInterface
      IMPLICIT NONE

! -----------------------------
!     Includes
!
      INCLUDE 'mpif.h'
      INCLUDE 'MSG_f90.h'

      INCLUDE 'BMG_workspace_f90.h'
      INCLUDE 'BMG_constants_f90.h'
      INCLUDE 'BMG_stencils_f90.h'
      INCLUDE 'BMG_parameters_f90.h'

! ----------------------------
!     Argument Declarations
!
      integer(len_t), value :: II, JJ, NMSGr
      integer(c_int), value :: nog, nstncl
      integer(len_t), value :: iGs, jGs
      integer(c_int), value :: IRELAX_SYM, K, UPDOWN
      integer, value :: MPICOMM
      integer(c_int), value :: XLINECOMM, YLINECOMM
      integer(len_t) :: DATADIST(2,*)
      real(real_t) :: B(II,JJ), Q(II,JJ), QF(II,JJ), SO(II+1,JJ+1,NStncl), &
           SOR(II,JJ,2), RWORK(NMSGr)
      type(c_ptr) :: halof

! ----------------------------
!     Local Declarations
!
      INTEGER I, I1, J, J1, JBEG, JEND
      INTEGER JBEG_START, JBEG_END, JBEG_STRIDE
      INTEGER ierror, myid, size

      INTEGER NP, IERR, OFFSET, NLINES
      INTEGER STATUS(MPI_STATUS_SIZE), MULT, TAG
      INTEGER ptrn
      INTEGER myidy, sizey
      INTEGER line_solve_comm_type

! ======================================================================

      line_solve_comm_type = BMG_LINE_SOLVE_COMM_TRADITIONAL
      TAG = 3

      J1=JJ-1
      I1=II-1
!
!     Local Declarations
!     on the way down we relax red lines and then black lines
!     on the way up we relax in the opposite order
!

      call MPI_Comm_Rank(XLINECOMM,myid,ierror)
      call MPI_Comm_Size(XLINECOMM,size,ierror)

      call MPI_Comm_Rank(YLINECOMM,myidy,ierror)
      call MPI_Comm_Size(YLINECOMM,sizey,ierror)
!
      IF ( UPDOWN.EQ.BMG_DOWN .OR. IRELAX_SYM.EQ.BMG_RELAX_NONSYM ) THEN
         !
         !  Red, black => (start,end,stride) = (3,2,-1)
         !  Black, red => (start,end,stride) = (2,3,1)
         !
         JBEG_START  = 2 + MOD(jGs,2)
         JBEG_END    = 2 + MOD(jGs+1,2)
         JBEG_STRIDE = MOD(jGs+1,2) - MOD(jGs,2)
         !
      ELSEIF ( IRELAX_SYM.EQ.BMG_RELAX_SYM ) THEN
         !
         JBEG_START  = 2 + MOD(jGs+1,2)
         JBEG_END    = 2 + MOD(jGs,2)
         JBEG_STRIDE = MOD(jGs,2) - MOD(jGs+1,2)
         !
      ENDIF

      DO JBEG=JBEG_START, JBEG_END, JBEG_STRIDE

         IF ( NStncl.EQ.5 ) THEN
            !
            DO J=JBEG,J1,2
               DO I=2,I1
                  !
                  Q(I,J) = QF(I,J)&
     &                   + SO(I,J,KS)*Q(I,J-1)&
     &                   + SO(I,J+1,KS)*Q(I,J+1)&
     &                   + SO(I,J,KSW)*Q(I-1,J-1)&
     &                   + SO(I+1,J,KNW)*Q(I+1,J-1)&
     &                   + SO(I,J+1,KNW)*Q(I-1,J+1)&
     &                   + SO(I+1,J+1,KSW)*Q(I+1,J+1)
                  !
               ENDDO
            ENDDO
            !
         ELSE
            !
            DO J=JBEG,J1,2
               DO I=2,I1
                  Q(I,J) = QF(I,J)&
     &                   + SO(I,J,KS)*Q(I,J-1)&
     &                   + SO(I,J+1,KS)*Q(I,J+1)
               ENDDO
            ENDDO
            !
         ENDIF

         ! store the index of the last line
         JEND = J-2

         ! ====================================================

         NP = II-2
         CALL BMG2_SymStd_LineSolve_A(SOR,Q,II,JJ,&
     &        JBEG, RWORK, NP, NLINES)

         ! =====================================================

         IF (SIZE.gt.1) THEN    ! if there is more than one processor in
            if (myid .eq. 0) then
               CALL MPI_GATHER(MPI_IN_PLACE,NLINES*8,&
                    &           MPI_DOUBLE_PRECISION,&
                    &           RWORK,NLINES*8,MPI_DOUBLE_PRECISION,&
                    &           0, XLINECOMM, IERR)
            else
               CALL MPI_GATHER(RWORK,NLINES*8,&
                    &           MPI_DOUBLE_PRECISION,&
                    &           0.0d0,NLINES*8,MPI_DOUBLE_PRECISION,&
                    &           0, XLINECOMM, IERR)
            endif
         ENDIF

         ! ======================================================

         MULT = 0
         DO J=JBEG,J1,2

            CALL BMG2_SymStd_LineSolve_B (SOR(2,J,1), SOR(3,J,2), &
     &           SOR(2,J,2), Q(2,J), &
     &           RWORK(SIZE*NLINES*8 + 1),&
     &           RWORK(MULT*8+1), &
     &           RWORK(MULT*8+1), &
     &           NP, size, DATADIST, myid, 8*NLINES, 8*NLINES)

            MULT = MULT+1

         ENDDO

         ! ====================================================

         IF (SIZE.gt.1) THEN
            if (myid .eq. 0) then
               CALL MPI_SCATTER(RWORK(MYID*NLINES*8+1),NLINES*8,&
                    &           MPI_DOUBLE_PRECISION,MPI_IN_PLACE, NLINES*8, &
                    &           MPI_DOUBLE_PRECISION,0,XLINECOMM,IERR)
            else
               CALL MPI_SCATTER(0.0d0,NLINES*8,&
                    &           MPI_DOUBLE_PRECISION,RWORK, NLINES*8, &
                    &           MPI_DOUBLE_PRECISION,0,XLINECOMM,IERR)
            endif
         ENDIF

         ! ====================================================

         MULT = 0
         DO J=JBEG,J1,2

            CALL BMG2_SymStd_LineSolve_C (SOR(2,J,1), SOR(3,J,2), &
     &           SOR(2,J,2), Q(1,J),  &
     &           RWORK(MULT*8 + 1),&
     &           RWORK(SIZE*NLINES*8 + 1),&
     &           NP, size, DATADIST,&
     &           myid)

            MULT = MULT+1

         ENDDO

         ! ====================================================

         IF (line_solve_comm_type .EQ.&
     &        BMG_LINE_SOLVE_COMM_TUNED)  THEN

            IF (SIZEY .GT. 1) THEN

               IF (MYIDY .EQ. 0) THEN
                  ! send to the top
                  IF (JEND .EQ. JJ-1) THEN
                     CALL MPI_Send(Q(1,JJ-1),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY+1,0,YLINECOMM,IERR)
                  ENDIF
               ENDIF

               IF (MYIDY .EQ. SIZEY-1) THEN
                  ! receive from the bottom
                  IF (JBEG .EQ. 3) THEN
                     CALL MPI_Recv(Q(1,1),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY-1,0,YLINECOMM,STATUS,IERR)
                  ENDIF
               ENDIF

               IF (MYIDY .GT. 0 .AND. MYIDY .LT. SIZEY-1) THEN
                  ! send to the top and receive from the bottom
                  IF (JBEG .EQ. 3 .AND. JEND .EQ. JJ-1) THEN
                     CALL MPI_Sendrecv(Q(1,JJ-1),II,&
     &                    MPI_DOUBLE_PRECISION,&
     &                    MYIDY+1,0,Q(1,1),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY-1,0,YLINECOMM,STATUS,IERR)

                  ELSE IF (JBEG .NE. 3 .AND. JEND .EQ. JJ-1) THEN
                     ! send to the top
                     CALL MPI_Send(Q(1,JJ-1),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY+1,0,YLINECOMM,IERR)

                  ELSE IF (JBEG .EQ. 3 .AND. JEND .NE. JJ-1) THEN
                     ! receive from the bottom
                     CALL MPI_Recv(Q(1,1),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY-1,0,YLINECOMM,STATUS,IERR)
                  ENDIF

               ENDIF

               IF (MYIDY .EQ. SIZEY-1) THEN
                  ! send to the bottom
                  IF (JBEG .EQ. 2) THEN
                     CALL MPI_Send(Q(1,2),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY-1,0,YLINECOMM,IERR)
                  ENDIF
               ENDIF

               IF (MYIDY .EQ. 0) THEN
                  ! receive from the top
                  IF (JEND .EQ. JJ-2) THEN
                     CALL MPI_Recv(Q(1,JJ),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY+1,0,YLINECOMM,STATUS,IERR)
                  ENDIF
               ENDIF

               IF (MYIDY .GT. 0 .AND. MYIDY .LT. SIZEY-1) THEN
                  ! send to the bottom and receive from the top
                  IF (JBEG .EQ. 2 .AND. JEND .EQ. JJ-2) THEN
                     CALL MPI_Sendrecv(Q(1,2),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY-1,0,Q(1,JJ),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY+1,0,YLINECOMM,STATUS,IERR)

                  ELSE IF (JBEG .EQ. 2 .AND. JEND .NE. JJ-2) THEN
                     ! send to the bottom
                     CALL MPI_Send(Q(1,2),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY-1,0,YLINECOMM,IERR)

                  ELSE IF (JBEG .NE. 2 .AND. JEND .EQ. JJ-2) THEN
                     ! receive from the top
                     CALL MPI_Recv(Q(1,JJ),II,MPI_DOUBLE_PRECISION,&
     &                    MYIDY+1,0,YLINECOMM,STATUS,IERR)
                  ENDIF

               ENDIF

            ENDIF

         ELSE  IF (line_solve_comm_type .EQ.&
              &           BMG_LINE_SOLVE_COMM_TRADITIONAL)  THEN

            call halo_exchange(K, Q, halof);
         ELSE

            WRITE(*,*) 'ERROR: invalid value for parameter'
            WRITE(*,*) '      BMG(id_BMG2_LINE_SOLVE_COMM_TYPE)=',&
     &           line_solve_comm_type
            STOP

         ENDIF

      ENDDO

! ======================================================================

      RETURN
      END

max_stars_repo_name: cedar-framework/cedar

max_issues_repo_name: cedar-framework/cedar

max_forks_repo_name: cedar-framework/cedar


parameter(k=3000000)
dimension x(k),y(k),q(k),iflag(k),ip(k),ngr(200)
double precision d,dr,x,y,p,dr1,dr2,dd,df,xr,yr,xx,yy

p=3.1415926535879/180.
p=p*1.d+0

open(80,file='nz.dat')
open(90,file='nz_bmap_500_new.dat')

neve=500

kk=0
do i=1,3000000
  read(80,*,end=99)ia,m,ig,io,mi,s,y(i),x(i),pp,q(i)   !,id   t,x(i),y(i),z,q(i)   !
  iflag(i)=0
end do
99 npt=i-1

print*,npt

ke=1
ntot=0
do ic=1,100000000
  ieve=int(rand()*float(npt))
  if(iflag(ieve).ne.0)cycle
  kk=kk+1
  if(kk==6626)exit
  xr=x(ieve)*p
  yr=y(ieve)*p
  d=10.d+0
  do iv=1,500
    do j=1,k
      ip(j)=0
    end do
    ne=0
    ds=0.
    do j=1,npt
      if(iflag(j).ne.0)cycle
       xx=x(j)*p
       yy=y(j)*p
       df=abs(yy-yr)
       dr=dacos(dsin(xx)*dsin(xr)+dcos(xx)*dcos(xr)*dcos(df))
       dr=dr*6370.d+0
      if(dr<d)then
        ds=ds+dr
        ne=ne+1
        if(ne>k)print*,ne,'*****'
        ip(ne)=j
      end if
    end do
    if(ne==0)exit
    dm=ds/float(ne)
    s=dm*.1
    if(ne==1)s=5.
    if(neve-50<=ne.and.ne<=neve+50)then
      dd=d
      d=10.d+0
      exit
    end if
    if(ne<neve)d=d+s !1.d-2
    if(ne>neve)d=d-s  !1.d-2
  end do
  if(ne<neve-50.or.ne>neve+50)cycle
  do j=1,200
    ngr(j)=0
  end do
  qmax=-20.
  do j=1,ne
    if(q(ip(j))>qmax)qmax=q(ip(j))
    ll=int(q(ip(j))*10.)
    if(ll>200)print*,ll,'ll'
    ngr(ll)=ngr(ll)+1
  end do
  nmax=0
  do j=1,200
    if(ngr(j)>nmax)then
      nmax=ngr(j)
      jmax=j
    end if
  end do
  qc=jmax*.1
  if(qmax-qc<2.)cycle
  print*,ic,iv,ne,ke,kk,ntot
  do j=1,ne
    iflag(ip(j))=ke
  end do
  ke=ke+1
  ntot=ntot+ne

  if(float(npt-ntot)/float(npt)<.01)exit
end do
print*,npt-ntot,float(npt-ntot)/float(npt)

do i=1,npt
  write(90,'(2(e10.4,1x),i4,1x,f3.1)')x(i),y(i),iflag(i),q(i)
end do

end


!----------------------------------------------------------------------!
subroutine SHExpandLSQ_Wrapper(cilm, d, lat, lon, nmax, lmax, chi2)
    ! This wrapper ensures that the input arrays have
    ! correct dimensions when calling SHExpandLSQ
    use SHTOOLS, only: SHExpandLSQ

    implicit none
    real*8, intent(in) :: d(nmax), lat(nmax), lon(nmax)
    real*8, intent(out):: cilm(2, lmax + 1, lmax + 1)
    integer,intent(in) :: nmax, lmax
    real*8, intent(out):: chi2
    integer :: norm, csphase, exitstatus
    norm = 4
    csphase = 1
    exitstatus = 0

    call SHExpandLSQ(cilm, d, lat, lon, nmax, lmax, norm, chi2, csphase, &
        exitstatus)

end subroutine SHExpandLSQ_Wrapper

!---------------------------------------------------------------!
subroutine MakeGridPoints_Wrapper(cilm, lmax, n, lat, lon, points, dealloc)
    ! This wrapper ensures that the input arrays have
    ! correct dimensions when calling MakeGridPoint
    use SHTOOLS, only: MakeGridPoint

    implicit none
    integer,intent(in) :: lmax, n, dealloc
    real*8, intent(in) :: cilm(2, lmax + 1, lmax + 1)
    real*8, intent(in) :: lat(n), lon(n)
    real*8, intent(out) :: points(n)
    integer :: norm, csphase, i
    norm = 4
    csphase = 1
    points = 0.0d0

    do i = 1, n
    points(i) = MakeGridPoint( cilm, lmax, lat(i), lon(i),&
        norm, csphase, dealloc)
    end do

end subroutine MakeGridPoints_Wrapper

!---------------------------------------------------------------------!
subroutine GLQGridCoord_Wrapper(latglq, longlq, lmax, nlat, nlong)
    use SHTOOLS, only: GLQGridCoord

    implicit none

    integer, intent(in) :: lmax
    integer, intent(out) :: nlat, nlong
    real*8, intent(out) :: latglq(lmax + 1), longlq(2*lmax + 1)

    call GLQGridCoord(latglq, longlq, lmax, nlat, nlong)

end subroutine GLQGridCoord_Wrapper

!----------------------------------------------------------------------!
subroutine SHGLQ_Wrapper(lmax, zero, w, plx)
    use SHTOOLS, only: SHGLQ

    implicit none

    integer, intent(in) :: lmax
    real*8, intent(out) :: zero(lmax + 1), w(lmax + 1)
    real*8, intent(out) :: plx(1:lmax+1, 1:(lmax+1)*(lmax+2)/2)
    integer :: norm, csphase, cnorm
    norm = 4
    csphase = 1
    cnorm = 0
    call SHGLQ(lmax, zero, w, plx, norm, csphase, cnorm)

end subroutine SHGLQ_Wrapper

!----------------------------------------------------------------------!
subroutine SHExpandGLQ_Wrapper(cilm, lmax, gridglq, w, plx)
    use SHTOOLS, only: SHExpandGLQ

    implicit none

    integer, intent(in) :: lmax
    real*8, intent(in) :: w(lmax + 1), gridglq(lmax + 1, 2*lmax + 1)
    real*8, intent(in) :: plx(1:lmax+1, 1:(lmax+1)*(lmax+2)/2)
    real*8, intent(out) :: cilm(2,lmax + 1,lmax + 1)
    integer :: norm, csphase
    norm = 4
    csphase = 1

    call SHExpandGLQ(cilm, lmax, gridglq, w, plx, norm = norm, &
        csphase = csphase)
end subroutine SHExpandGLQ_Wrapper

!----------------------------------------------------------------------!
subroutine MakeGridGLQ_Wrapper(gridglq, cilm, lmax, plx)
    use SHTOOLS, only: MakeGridGLQ

    implicit none

    integer, intent(in) :: lmax
    real*8, intent(in) :: cilm(2,lmax + 1,lmax + 1)
    real*8, intent(in) :: plx(1:lmax+1, 1:(lmax+1)*(lmax+2)/2)
    real*8, intent(out) :: gridglq(lmax + 1, 2*lmax + 1)
    integer :: norm, csphase
    norm = 4
    csphase = 1

    call MakeGridGLQ(gridglq, cilm, lmax, plx, norm = norm, csphase = csphase)

end subroutine MakeGridGLQ_Wrapper

!----------------------------------------------------------------------!
subroutine SHExpandDH_Wrapper(grid, n, cilm, lmax)
    use SHTOOLS, only: SHExpandDH
    implicit none

    integer, intent(in) :: n
    real*8, intent(in) :: grid(n,n)
    real*8, intent(out) :: cilm(2, n/2, n/2)
    integer, intent(out) :: lmax
    integer :: norm, sampling, csphase, lmax_calc
    norm = 4
    sampling = 1
    csphase = 1

    call SHExpandDH(grid, n, cilm, lmax, norm, sampling, csphase)
end subroutine SHExpandDH_Wrapper

!----------------------------------------------------------------------!
subroutine MakeGridDH_Wrapper(grid, n, cilm, lmax)
    use SHTOOLS, only: MakeGridDH
    implicit none

    integer, intent(in) :: lmax
    real*8, intent(in) :: cilm(2,lmax + 1,lmax + 1)
    integer, intent(out) :: n
    real*8, intent(out) :: grid(2*lmax + 2,2*lmax + 2)
    integer :: norm, sampling, csphase
    norm = 4
    sampling = 1
    csphase = 1

    call MakeGridDH(grid, n, cilm, lmax, norm, sampling, csphase)

end subroutine MakeGridDH_Wrapper

!----------------------------------------------------------------------!
subroutine SHPowerSpectrum_Wrapper(cilm, lmax, pspectrum)
    use SHTOOLS, only: SHPowerSpectrum

    implicit none
    real*8, intent(in) :: cilm(2, lmax + 1, lmax + 1)
    integer, intent(in) :: lmax
    real*8, intent(out) :: pspectrum(lmax + 1)

    call SHPowerSpectrum(cilm, lmax, pspectrum)

end subroutine SHPowerSpectrum_Wrapper

!----------------------------------------------------------------------!
subroutine PlmON_Wrapper(p, lmax, z)
    use SHTools, only: PlmON

    implicit none
    integer, intent(in) :: lmax
    real*8, intent(out) :: p( (lmax + 1)*(lmax + 2)/2 )
    real*8, intent(in)  :: z
    integer :: csphase, cnorm
    csphase = 1
    cnorm = 0

    call PlmON(p, lmax, z, csphase, cnorm)

end subroutine PlmON_Wrapper


text:
module f90deviceio
  interface
    subroutine f90print(N)
      character(*) :: N
      !$omp declare target (f90print)
    end subroutine f90print
    subroutine f90printi(N,i)
      character(*) :: N
      integer :: i
      !$omp declare target (f90printi)
    end subroutine f90printi
    subroutine f90printf(N,f)
      character(*) :: N
      real(4) :: f
      !$omp declare target (f90printf)
    end subroutine f90printf
    subroutine f90printd(N,d)
      character(*) :: N
      real(8) :: d
      !$omp declare target (f90printd)
   end subroutine f90printd
  end interface
end module

max_stars_repo_name: raramakr/aomp

max_issues_repo_name: raramakr/aomp

max_forks_repo_name: raramakr/aomp


program main
    use mod_one
    implicit none
    call one()
end program main


text:
! make time available in a module

module time_module

  use bl_types

  implicit none

  real(dp_t), save :: time

  private

  public :: time

end module time_module

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO


text:
Bears dont exist in Davis, but there have been occasional sightings nearby:

 1910
A 700 lb Cinnamon Bear was killed about four miles west of Rumsey. In the article it was noted that 50 years prior grizzly bears were quite plentiful in this section and that a man or two was killed within three miles of Davis.

 1999
Graduate student Melanie Truan described the sighting of a 200pound black bear by four colleagues who were canoeing down Putah Creek on May 20 of this year. The bear swam away, but Truan went out to the scene hoping to make a plaster cast of a bear track.
To my great luck and delight I found tracks in the mud and made a cast, she said, passing the cast around the room in the Alumni and Visitors Center on the UCD campus. Putah Creek News, Winter 1999.

  On a hike along Putah creek in the spring 2010, I saw bear tracks and what sure looked to be an abandoned lair.  They were on the north side, just west of the University Airport.

Hominid bears locally include the big hairy Rainbow Community gay men who use that label to refer to themselves?

You might be looking for
  BEAR License
  Black Bear Diner  Actually, bears do exist in Davis. They sell teddy bears in a can here.
  Users/Bear Users/Bear

However, bears do exist in other places Adventures Outside of Davis you may go outside of Davis, so do use caution when you leave the protective Davis borders.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! MIT License
!
! Copyright (c) 2016 Anders Steen Christensen, Guido Falk von Rudorff
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

! prototype external calculator of the kernel matrix and solver for alpha coefficients
! makes use of a PBLAS implementation, here the MKL variant
! core idea is to distribute the kernel matrix across MPI ranks to allow out-of-core datasets to be processed
!
! CSCS:
!   load modules PrgEnv-cray/6.0.4 cray-mpich/7.6.0 daint-gpu intel/17.0.4.196
!   build with ifort driver_mpi.f90 -o driver_mpi  -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_scalapack_lp64 -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lmkl_blacs_openmpi_lp64 -lpthread -lm -ldl  -m64 -I${MKLROOT}/include -Og -g -Wall -Wextra -pedantic -fimplicit-none -fcheck=all -fbacktrace
program qml_driver
    use mpi
    implicit none

    double precision, allocatable, dimension(:,:) :: Q
    double precision, allocatable, dimension(:):: Y

    double precision :: sigma
    integer :: rep_size
    integer :: n_molecules

    double precision, allocatable, dimension(:) :: alphas

    integer :: info
    integer :: na, i, numroc
    integer :: local_i, local_j, global_i, global_j

    integer :: lwork
    double precision, dimension(:), allocatable :: work
    double precision, dimension(:,:), allocatable :: local_K, local_B

    ! pblacs data structures
    integer :: local_K_rows, local_K_cols
    integer :: local_B_rows, local_B_cols
    integer :: local_id, num_ranks, ranks_rows, ranks_cols, context
    integer :: local_rank_col, local_rank_row, block_size
    integer :: ierr
    integer, dimension(9) :: desca, descb
    integer, dimension(2) :: dims
    double precision, dimension(2) :: work_query

    ! init pblacs and set cartesian grid of MPI ranks
    call blacs_pinfo(local_id, num_ranks)

    dims = 0
    call MPI_Dims_create(num_ranks, 2, dims, ierr)
    ranks_rows = dims(1)
    ranks_cols = dims(2)

    block_size = 100

    ! create BLACS context
    call blacs_get(0, 0, context)
    call blacs_gridinit(context, 'R', ranks_rows, ranks_cols)
    call blacs_gridinfo(context, ranks_rows, ranks_cols, &
        local_rank_row, local_rank_col)

    ! Read hyperparameters and arrat suzes
    open(unit = 9, file = "parameters.fout", form="formatted")

        read(9,*) sigma, rep_size, n_molecules

        ! Allocate labels
        allocate(Y(n_molecules))

        read(9,*) Y(:n_molecules)

    close(9)

    ! allocate alphas
    allocate(alphas(n_molecules))
    ! Allocate representations
    allocate(Q(n_molecules, rep_size))

    ! Read representations
    open(unit = 9, file = "representations.fout", form="formatted")

        ! Read representaions for each molecule
        do i = 1, n_molecules
            read(9,*) Q(i,:rep_size)
        enddo

    close(9)

    ! Size of kernel
    na = size(Q, dim=1)

    ! Allocate kernel and output
    local_K_rows = numroc(na, block_size, local_rank_row, 0, ranks_rows)
    local_K_cols = numroc(na, block_size, local_rank_col, 0, ranks_cols)
    allocate(local_K(local_K_rows, local_K_cols))
    local_B_rows = numroc(na, block_size, local_rank_row, 0, ranks_rows)
    local_B_cols = numroc(1, block_size, local_rank_col, 0, ranks_cols)
    allocate(local_B(local_B_rows, local_B_cols))

    ! Calculate Laplacian kernel
    do local_j = 1, local_K_cols
        call l2g(local_j, local_rank_col, ranks_cols, &
                block_size, global_j)
        do local_i = 1, local_K_rows
            call l2g(local_i, local_rank_row, ranks_rows, &
                block_size, global_i)
            local_K(local_i, local_j) = exp(-sum(abs(Q(global_j,:) - Q(global_i,:)))/sigma)
        enddo
    enddo

    ! Setup variables for LAPACK
    call descinit(desca, na, na, block_size, block_size, 0, 0, context, MAX(1, local_K_rows), info)
    call descinit(descb, na, 1, block_size, block_size, 0, 0, context, MAX(1, local_B_rows), info)

    ! Allocate local work arrays
    call pdgels("N", na, na, 1, local_K, 1, 1, desca, local_B, 1, 1, DESCB, work_query, -1, info)
    lwork = INT(work_query(1))
    allocate(work(lwork))

    ! copy data
    local_B = 0.0d0
    if (local_B_cols .ne. 0) then
        do local_i = 1, local_B_rows
            call l2g(local_i, local_rank_row, ranks_rows, block_size, global_i)
            local_B(local_i, 1) = y(global_i)
        enddo
    end if

    ! Solver
    call pdgels("N", na, na, 1, local_K, 1, 1, desca, local_B, 1, 1, DESCB, work, lwork, info)

    ! Copy LAPACK output
    alphas = 0.0d0
    if (local_B_cols .ne. 0) then
        do local_i = 1, local_B_rows
            call l2g(local_i, local_rank_row, ranks_rows, block_size, global_i)
            alphas(global_i) = local_B(local_i, 1)
        enddo
    endif
    call DGSUM2D(context, "All", "1-tree", na, 1, alphas, 1, -1, -1)

    ! Save alphas to file
    if (local_id .eq. 0) then
        open(unit = 9, file = "alphas_mpi.fout", form="formatted")
            write(9,*) alphas(:)
        close(9)
    end if

    ! Tear down MPI
    call blacs_exit(0)

    ! Clean up
    deallocate(work)
    deallocate(local_B)
    deallocate(Q)
    deallocate(local_K)
    deallocate(Y)
    deallocate(alphas)

end program qml_driver
! convert local index to global index in block-cyclic distribution

   subroutine l2g(il,p,np,nb,i)

   implicit none
   integer :: il   ! local array index, input
   integer :: p    ! processor array index, input
   integer :: np   ! processor array dimension, input
   integer :: nb   ! block size, input
   integer :: i    ! global array index, output
   integer :: ilm1

   ilm1 = il-1
   i    = (((ilm1/nb) * np) + p)*nb + mod(ilm1,nb) + 1

   return
   end subroutine l2g

max_stars_repo_name: charnley/qml

max_issues_repo_name: charnley/qml

max_forks_repo_name: charnley/qml


      subroutine salt_cha_res_read

      use constituent_mass_module
      use input_file_module
      use maximum_data_module
      use channel_data_module
      use hydrograph_module
      use sd_channel_module
      use organic_mineral_mass_module

      implicit none

      character (len=80) :: titldum
      character (len=80) :: header
      integer :: eof, imax
      logical :: i_exist
      integer :: isalt
      integer :: isalti

      eof = 0

      !read all export coefficient data
      inquire (file=in_init%salt_water, exist=i_exist)
      if (i_exist .or. in_init%salt_water /= "null") then
        do
          open (107,file=in_init%salt_water)
          read (107,*,iostat=eof) titldum
          if (eof < 0) exit
          read (107,*,iostat=eof) header
          if (eof < 0) exit
          imax = 0
          do while (eof == 0)
            read (107,*,iostat=eof) titldum     !name
            if (eof < 0) exit
            do isalt = 1, cs_db%num_salts
              read (107,*,iostat=eof) titldum   !water
              if (eof < 0) exit
              read (107,*,iostat=eof) titldum   !benthic
              if (eof < 0) exit
            end do
            imax = imax + 1
          end do

          db_mx%saltw_ini = imax

          allocate (salt_water_ini(imax))
          allocate (salt_init_name(imax))

          do isalt = 1, imax
            allocate (salt_water_ini(isalt)%water(cs_db%num_salts+5))
            allocate (salt_water_ini(isalt)%benthic(cs_db%num_salts+5))
          end do

          rewind (107)
          read (107,*,iostat=eof) titldum
          if (eof < 0) exit

          do isalti = 1, imax
            read (107,*,iostat=eof) header
            if (eof < 0) exit
            read (107,*,iostat=eof) salt_init_name(isalti)
            if (eof < 0) exit
            read (107,*,iostat=eof) titldum, salt_water_ini(isalti)%water
            if (eof < 0) exit
            read (107,*,iostat=eof) titldum, salt_water_ini(isalti)%benthic
            if (eof < 0) exit
          end do
          close (107)
          exit
        end do
      end if

      return
      end subroutine salt_cha_res_read