text:
!-------------------------------------------------------------------------------
!>
!! 3D Operator module
!!
!! @par Description
!!         This module contains the subroutines for differential oeprators using vertical metrics.
!!
!! @author  H.Tomita
!!
!! @par History
!! @li      2004-02-17 (H.Tomita)    Imported from igdc-4.33
!! @li      2011-09-27 (T.Seiki)     merge optimization by RIST and M.Terai
!!
!<
module mod_oprt3d
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use mod_precision
  use mod_debug
  use mod_adm, only: &
     ADM_LOG_FID,    &
     TI  => ADM_TI,  &
     TJ  => ADM_TJ,  &
     AI  => ADM_AI,  &
     AIJ => ADM_AIJ, &
     AJ  => ADM_AJ,  &
     K0  => ADM_KNONE
  use mod_gmtr, only: &
     P_RAREA => GMTR_P_RAREA, &
     T_RAREA => GMTR_T_RAREA, &
     W1      => GMTR_T_W1,    &
     W2      => GMTR_T_W2,    &
     W3      => GMTR_T_W3,    &
     HNX     => GMTR_A_HNX,   &
     HNY     => GMTR_A_HNY,   &
     HNZ     => GMTR_A_HNZ,   &
     HTX     => GMTR_A_HTX,   &
     HTY     => GMTR_A_HTY,   &
     HTZ     => GMTR_A_HTZ,   &
     TNX     => GMTR_A_TNX,   &
     TNY     => GMTR_A_TNY,   &
     TNZ     => GMTR_A_TNZ,   &
     TN2X    => GMTR_A_TN2X,  &
     TN2Y    => GMTR_A_TN2Y,  &
     TN2Z    => GMTR_A_TN2Z
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: OPRT3D_divdamp

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  !-----------------------------------------------------------------------------
contains
  !-----------------------------------------------------------------------------
  subroutine OPRT3D_divdamp( &
       ddivdx, ddivdx_pl, &
       ddivdy, ddivdy_pl, &
       ddivdz, ddivdz_pl, &
       rhogvx, rhogvx_pl, &
       rhogvy, rhogvy_pl, &
       rhogvz, rhogvz_pl, &
       rhogw,  rhogw_pl   )
    use mod_adm, only: &
       ADM_have_pl,  &
       ADM_have_sgp, &
       ADM_lall,     &
       ADM_lall_pl,  &
       ADM_gall,     &
       ADM_gall_pl,  &
       ADM_kall,     &
       ADM_gall_1d,  &
       ADM_gmin,     &
       ADM_gmax,     &
       ADM_gslf_pl,  &
       ADM_gmin_pl,  &
       ADM_gmax_pl,  &
       ADM_kmin,     &
       ADM_kmax
    use mod_grd, only: &
       GRD_rdgz
    use mod_oprt, only: &
       cinterp_TN,     &
       cinterp_TN_pl,  &
       cinterp_HN,     &
       cinterp_HN_pl,  &
       cinterp_TRA,    &
       cinterp_TRA_pl, &
       cinterp_PRA,    &
       cinterp_PRA_pl
    use mod_vmtr, only: &
       VMTR_RGAM,        &
       VMTR_RGAM_pl,     &
       VMTR_RGAMH,       &
       VMTR_RGAMH_pl,    &
       VMTR_RGSQRTH,     &
       VMTR_RGSQRTH_pl,  &
       VMTR_C2WfactGz,   &
       VMTR_C2WfactGz_pl
    implicit none

    real(RP), intent(out) :: ddivdx   (ADM_gall   ,ADM_kall,ADM_lall   ) ! tendency
    real(RP), intent(out) :: ddivdx_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: ddivdy   (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP), intent(out) :: ddivdy_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(out) :: ddivdz   (ADM_gall   ,ADM_kall,ADM_lall   )
    real(RP), intent(out) :: ddivdz_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: rhogvx   (ADM_gall   ,ADM_kall,ADM_lall   ) ! rho*vx { gam2 x G^1/2 }
    real(RP), intent(in)  :: rhogvx_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: rhogvy   (ADM_gall   ,ADM_kall,ADM_lall   ) ! rho*vy { gam2 x G^1/2 }
    real(RP), intent(in)  :: rhogvy_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: rhogvz   (ADM_gall   ,ADM_kall,ADM_lall   ) ! rho*vz { gam2 x G^1/2 }
    real(RP), intent(in)  :: rhogvz_pl(ADM_gall_pl,ADM_kall,ADM_lall_pl)
    real(RP), intent(in)  :: rhogw    (ADM_gall   ,ADM_kall,ADM_lall   ) ! rho*w  { gam2 x G^1/2 }
    real(RP), intent(in)  :: rhogw_pl (ADM_gall_pl,ADM_kall,ADM_lall_pl)

    real(RP) :: sclt         (ADM_gall   ,TI:TJ) ! scalar on the hexagon vertex
    real(RP) :: sclt_pl      (ADM_gall_pl)
    real(RP) :: sclt_rhogw
    real(RP) :: sclt_rhogw_pl

    real(RP) :: rhogvx_vm   (ADM_gall   )          ! rho*vx / vertical metrics
    real(RP) :: rhogvx_vm_pl(ADM_gall_pl)
    real(RP) :: rhogvy_vm   (ADM_gall   )          ! rho*vy / vertical metrics
    real(RP) :: rhogvy_vm_pl(ADM_gall_pl)
    real(RP) :: rhogvz_vm   (ADM_gall   )          ! rho*vz / vertical metrics
    real(RP) :: rhogvz_vm_pl(ADM_gall_pl)
    real(RP) :: rhogw_vm    (ADM_gall,   ADM_kall) ! rho*w  / vertical metrics
    real(RP) :: rhogw_vm_pl (ADM_gall_pl,ADM_kall)

    integer :: nstart1, nstart2, nend
    integer :: ij
    integer :: ip1j, ijp1, ip1jp1
    integer :: im1j, ijm1, im1jm1

    integer :: gall, gall_1d, gmin, kmin, kmax
    integer :: g, k, l, v, n
    !---------------------------------------------------------------------------

    call DEBUG_rapstart('OPRT3D_divdamp')

    gall    = ADM_gall
    gall_1d = ADM_gall_1d
    gmin    = ADM_gmin
    kmin    = ADM_kmin
    kmax    = ADM_kmax

    nstart1 = suf(ADM_gmin-1,ADM_gmin-1)
    nstart2 = suf(ADM_gmin  ,ADM_gmin  )
    nend    = suf(ADM_gmax  ,ADM_gmax  )

    do l = 1, ADM_lall
       do k = ADM_kmin+1, ADM_kmax
       do g = 1, ADM_gall
          rhogw_vm(g,k) = ( VMTR_C2WfactGz(g,k,1,l) * rhogvx(g,k  ,l) &
                          + VMTR_C2WfactGz(g,k,2,l) * rhogvx(g,k-1,l) &
                          + VMTR_C2WfactGz(g,k,3,l) * rhogvy(g,k  ,l) &
                          + VMTR_C2WfactGz(g,k,4,l) * rhogvy(g,k-1,l) &
                          + VMTR_C2WfactGz(g,k,5,l) * rhogvz(g,k  ,l) &
                          + VMTR_C2WfactGz(g,k,6,l) * rhogvz(g,k-1,l) &
                          ) * VMTR_RGAMH(g,k,l)                       & ! horizontal contribution
                        + rhogw(g,k,l) * VMTR_RGSQRTH(g,k,l)            ! vertical   contribution
       enddo
       enddo
       do g = 1, ADM_gall
          rhogw_vm(g,ADM_kmin  ) = 0.0_RP
          rhogw_vm(g,ADM_kmax+1) = 0.0_RP
       enddo

       !$omp parallel default(none), private(g,n,k,ij,ip1j,ip1jp1,ijp1,im1j,ijm1,im1jm1,sclt_rhogw), &
       !$omp shared(gall,gall_1d,gmin,kmin,kmax,nstart1,nstart2,nend,l,ADM_have_sgp, &
       !$omp GRD_rdgz,VMTR_RGAM,rhogvx,rhogvy,rhogvz,ddivdx,ddivdy,ddivdz,sclt, &
       !$omp rhogvx_vm,rhogvy_vm,rhogvz_vm,rhogw_vm,cinterp_TN,cinterp_HN,cinterp_PRA,cinterp_TRA)
       do k = kmin, kmax

          !$omp do
          do g = 1, gall
             rhogvx_vm(g) = rhogvx(g,k,l) * VMTR_RGAM(g,k,l)
             rhogvy_vm(g) = rhogvy(g,k,l) * VMTR_RGAM(g,k,l)
             rhogvz_vm(g) = rhogvz(g,k,l) * VMTR_RGAM(g,k,l)
          enddo
          !$omp end do

          !$omp do
          do n = nstart1, nend
             ij     = n
             ip1j   = n + 1
             ip1jp1 = n + 1 + gall_1d

             sclt_rhogw = ( ( rhogw_vm(ij,k+1) + rhogw_vm(ip1j,k+1) + rhogw_vm(ip1jp1,k+1) ) &
                          - ( rhogw_vm(ij,k  ) + rhogw_vm(ip1j,k  ) + rhogw_vm(ip1jp1,k  ) ) &
                          ) / 3.0_RP * GRD_rdgz(k)

             sclt(n,TI) = ( - (rhogvx_vm(ij    )+rhogvx_vm(ip1j  )) * cinterp_TN(ij  ,l,AI ,1) &
                            - (rhogvx_vm(ip1j  )+rhogvx_vm(ip1jp1)) * cinterp_TN(ip1j,l,AJ ,1) &
                            + (rhogvx_vm(ip1jp1)+rhogvx_vm(ij    )) * cinterp_TN(ij  ,l,AIJ,1) &
                            - (rhogvy_vm(ij    )+rhogvy_vm(ip1j  )) * cinterp_TN(ij  ,l,AI ,2) &
                            - (rhogvy_vm(ip1j  )+rhogvy_vm(ip1jp1)) * cinterp_TN(ip1j,l,AJ ,2) &
                            + (rhogvy_vm(ip1jp1)+rhogvy_vm(ij    )) * cinterp_TN(ij  ,l,AIJ,2) &
                            - (rhogvz_vm(ij    )+rhogvz_vm(ip1j  )) * cinterp_TN(ij  ,l,AI ,3) &
                            - (rhogvz_vm(ip1j  )+rhogvz_vm(ip1jp1)) * cinterp_TN(ip1j,l,AJ ,3) &
                            + (rhogvz_vm(ip1jp1)+rhogvz_vm(ij    )) * cinterp_TN(ij  ,l,AIJ,3) &
                          ) * 0.5_RP * cinterp_TRA(ij,l,TI) &
                        + sclt_rhogw
          enddo
          !$omp end do nowait

          !$omp do
          do n = nstart1, nend
             ij     = n
             ijp1   = n     + gall_1d
             ip1jp1 = n + 1 + gall_1d

             sclt_rhogw = ( ( rhogw_vm(ij,k+1) + rhogw_vm(ijp1,k+1) + rhogw_vm(ip1jp1,k+1) ) &
                          - ( rhogw_vm(ij,k  ) + rhogw_vm(ijp1,k  ) + rhogw_vm(ip1jp1,k  ) ) &
                          ) / 3.0_RP * GRD_rdgz(k)

             sclt(n,TJ) = ( - (rhogvx_vm(ij    )+rhogvx_vm(ip1jp1)) * cinterp_TN(ij  ,l,AIJ,1) &
                            + (rhogvx_vm(ip1jp1)+rhogvx_vm(ijp1  )) * cinterp_TN(ijp1,l,AI ,1) &
                            + (rhogvx_vm(ijp1  )+rhogvx_vm(ij    )) * cinterp_TN(ij  ,l,AJ ,1) &
                            - (rhogvy_vm(ij    )+rhogvy_vm(ip1jp1)) * cinterp_TN(ij  ,l,AIJ,2) &
                            + (rhogvy_vm(ip1jp1)+rhogvy_vm(ijp1  )) * cinterp_TN(ijp1,l,AI ,2) &
                            + (rhogvy_vm(ijp1  )+rhogvy_vm(ij    )) * cinterp_TN(ij  ,l,AJ ,2) &
                            - (rhogvz_vm(ij    )+rhogvz_vm(ip1jp1)) * cinterp_TN(ij  ,l,AIJ,3) &
                            + (rhogvz_vm(ip1jp1)+rhogvz_vm(ijp1  )) * cinterp_TN(ijp1,l,AI ,3) &
                            + (rhogvz_vm(ijp1  )+rhogvz_vm(ij    )) * cinterp_TN(ij  ,l,AJ ,3) &
                          ) * 0.5_RP * cinterp_TRA(ij,l,TJ) &
                        + sclt_rhogw
          enddo
          !$omp end do

          !$omp do
          do n = nstart2, nend
             ij     = n
             im1j   = n - 1
             ijm1   = n     - gall_1d
             im1jm1 = n - 1 - gall_1d

             ddivdx(n,k,l) = ( + ( sclt(ijm1,  TJ) + sclt(ij,    TI) ) * cinterp_HN(ij,    l,AI ,1) &
                               + ( sclt(ij,    TI) + sclt(ij,    TJ) ) * cinterp_HN(ij,    l,AIJ,1) &
                               + ( sclt(ij,    TJ) + sclt(im1j,  TI) ) * cinterp_HN(ij,    l,AJ ,1) &
                               - ( sclt(im1jm1,TJ) + sclt(im1j,  TI) ) * cinterp_HN(im1j,  l,AI ,1) &
                               - ( sclt(im1jm1,TI) + sclt(im1jm1,TJ) ) * cinterp_HN(im1jm1,l,AIJ,1) &
                               - ( sclt(ijm1  ,TJ) + sclt(im1jm1,TI) ) * cinterp_HN(ijm1,  l,AJ ,1) &
                             ) * 0.5_RP * cinterp_PRA(ij,l)

             ddivdy(n,k,l) = ( + ( sclt(ijm1,  TJ) + sclt(ij,    TI) ) * cinterp_HN(ij,    l,AI ,2) &
                               + ( sclt(ij,    TI) + sclt(ij,    TJ) ) * cinterp_HN(ij,    l,AIJ,2) &
                               + ( sclt(ij,    TJ) + sclt(im1j,  TI) ) * cinterp_HN(ij,    l,AJ ,2) &
                               - ( sclt(im1jm1,TJ) + sclt(im1j,  TI) ) * cinterp_HN(im1j,  l,AI ,2) &
                               - ( sclt(im1jm1,TI) + sclt(im1jm1,TJ) ) * cinterp_HN(im1jm1,l,AIJ,2) &
                               - ( sclt(ijm1  ,TJ) + sclt(im1jm1,TI) ) * cinterp_HN(ijm1,  l,AJ ,2) &
                             ) * 0.5_RP * cinterp_PRA(ij,l)

             ddivdz(n,k,l) = ( + ( sclt(ijm1,  TJ) + sclt(ij,    TI) ) * cinterp_HN(ij,    l,AI ,3) &
                               + ( sclt(ij,    TI) + sclt(ij,    TJ) ) * cinterp_HN(ij,    l,AIJ,3) &
                               + ( sclt(ij,    TJ) + sclt(im1j,  TI) ) * cinterp_HN(ij,    l,AJ ,3) &
                               - ( sclt(im1jm1,TJ) + sclt(im1j,  TI) ) * cinterp_HN(im1j,  l,AI ,3) &
                               - ( sclt(im1jm1,TI) + sclt(im1jm1,TJ) ) * cinterp_HN(im1jm1,l,AIJ,3) &
                               - ( sclt(ijm1  ,TJ) + sclt(im1jm1,TI) ) * cinterp_HN(ijm1,  l,AJ ,3) &
                             ) * 0.5_RP * cinterp_PRA(ij,l)
          enddo

          if ( ADM_have_sgp(l) ) then
             !$omp master
             n = suf(gmin,gmin)

             ij     = n
             im1j   = n - 1
             ijm1   = n     - gall_1d
             im1jm1 = n - 1 - gall_1d

             sclt(im1jm1,TI) = sclt(ijm1,TJ) ! copy

             ddivdx(n,k,l) = ( + ( sclt(ijm1,  TJ) + sclt(ij,    TI) ) * cinterp_HN(ij,    l,AI ,1) &
                               + ( sclt(ij,    TI) + sclt(ij,    TJ) ) * cinterp_HN(ij,    l,AIJ,1) &
                               + ( sclt(ij,    TJ) + sclt(im1j,  TI) ) * cinterp_HN(ij,    l,AJ ,1) &
                               - ( sclt(im1jm1,TJ) + sclt(im1j,  TI) ) * cinterp_HN(im1j,  l,AI ,1) &
                               - ( sclt(im1jm1,TI) + sclt(im1jm1,TJ) ) * cinterp_HN(im1jm1,l,AIJ,1) &
                             ) * 0.5_RP * cinterp_PRA(ij,l)

             ddivdy(n,k,l) = ( + ( sclt(ijm1,  TJ) + sclt(ij,    TI) ) * cinterp_HN(ij,    l,AI ,2) &
                               + ( sclt(ij,    TI) + sclt(ij,    TJ) ) * cinterp_HN(ij,    l,AIJ,2) &
                               + ( sclt(ij,    TJ) + sclt(im1j,  TI) ) * cinterp_HN(ij,    l,AJ ,2) &
                               - ( sclt(im1jm1,TJ) + sclt(im1j,  TI) ) * cinterp_HN(im1j,  l,AI ,2) &
                               - ( sclt(im1jm1,TI) + sclt(im1jm1,TJ) ) * cinterp_HN(im1jm1,l,AIJ,2) &
                             ) * 0.5_RP * cinterp_PRA(ij,l)

             ddivdz(n,k,l) = ( + ( sclt(ijm1,  TJ) + sclt(ij,    TI) ) * cinterp_HN(ij,    l,AI ,3) &
                               + ( sclt(ij,    TI) + sclt(ij,    TJ) ) * cinterp_HN(ij,    l,AIJ,3) &
                               + ( sclt(ij,    TJ) + sclt(im1j,  TI) ) * cinterp_HN(ij,    l,AJ ,3) &
                               - ( sclt(im1jm1,TJ) + sclt(im1j,  TI) ) * cinterp_HN(im1j,  l,AI ,3) &
                               - ( sclt(im1jm1,TI) + sclt(im1jm1,TJ) ) * cinterp_HN(im1jm1,l,AIJ,3) &
                             ) * 0.5_RP * cinterp_PRA(ij,l)
             !$omp end master
             !$omp barrier
          endif
       enddo
       !$omp end parallel

       do g = 1, ADM_gall
          ddivdx(g,ADM_kmin-1,l) = 0.0_RP
          ddivdx(g,ADM_kmax+1,l) = 0.0_RP
          ddivdy(g,ADM_kmin-1,l) = 0.0_RP
          ddivdy(g,ADM_kmax+1,l) = 0.0_RP
          ddivdz(g,ADM_kmin-1,l) = 0.0_RP
          ddivdz(g,ADM_kmax+1,l) = 0.0_RP
       enddo
    enddo

    if ( ADM_have_pl ) then
       do l = 1, ADM_lall_pl
          do k = ADM_kmin+1, ADM_kmax
          do g = 1, ADM_gall_pl
             rhogw_vm_pl(g,k) = ( VMTR_C2WfactGz_pl(g,k,1,l) * rhogvx_pl(g,k  ,l) &
                                + VMTR_C2WfactGz_pl(g,k,2,l) * rhogvx_pl(g,k-1,l) &
                                + VMTR_C2WfactGz_pl(g,k,3,l) * rhogvy_pl(g,k  ,l) &
                                + VMTR_C2WfactGz_pl(g,k,4,l) * rhogvy_pl(g,k-1,l) &
                                + VMTR_C2WfactGz_pl(g,k,5,l) * rhogvz_pl(g,k  ,l) &
                                + VMTR_C2WfactGz_pl(g,k,6,l) * rhogvz_pl(g,k-1,l) &
                                ) * VMTR_RGAMH_pl(g,k,l)                          & ! horizontal contribution
                              + rhogw_pl(g,k,l) * VMTR_RGSQRTH_pl(g,k,l)            ! vertical   contribution
          enddo
          enddo
          do g = 1, ADM_gall_pl
             rhogw_vm_pl(g,ADM_kmin  ) = 0.0_RP
             rhogw_vm_pl(g,ADM_kmax+1) = 0.0_RP
          enddo

          n = ADM_GSLF_PL

          do k = ADM_kmin, ADM_kmax
             do v = 1, ADM_gall_pl
                rhogvx_vm_pl(v) = rhogvx_pl(v,k,l) * VMTR_RGAM_pl(v,k,l)
                rhogvy_vm_pl(v) = rhogvy_pl(v,k,l) * VMTR_RGAM_pl(v,k,l)
                rhogvz_vm_pl(v) = rhogvz_pl(v,k,l) * VMTR_RGAM_pl(v,k,l)
             enddo

             do v = ADM_gmin_pl, ADM_gmax_pl
                ij   = v
                ijp1 = v + 1
                if( ijp1 > ADM_gmax_pl ) ijp1 = ADM_gmin_pl

                sclt_rhogw_pl = ( ( rhogw_vm_pl(n,k+1) + rhogw_vm_pl(ij,k+1) + rhogw_vm_pl(ijp1,k+1) ) &
                                - ( rhogw_vm_pl(n,k  ) + rhogw_vm_pl(ij,k  ) + rhogw_vm_pl(ijp1,k  ) ) &
                                ) / 3.0_RP * GRD_rdgz(k)

                sclt_pl(v) = ( + ( rhogvx_vm_pl(n   ) + rhogvx_vm_pl(ij  ) ) * cinterp_TN_pl(ij  ,l,1,1) &
                               + ( rhogvy_vm_pl(n   ) + rhogvy_vm_pl(ij  ) ) * cinterp_TN_pl(ij  ,l,1,2) &
                               + ( rhogvz_vm_pl(n   ) + rhogvz_vm_pl(ij  ) ) * cinterp_TN_pl(ij  ,l,1,3) &
                               + ( rhogvx_vm_pl(ij  ) + rhogvx_vm_pl(ijp1) ) * cinterp_TN_pl(ij  ,l,2,1) &
                               + ( rhogvy_vm_pl(ij  ) + rhogvy_vm_pl(ijp1) ) * cinterp_TN_pl(ij  ,l,2,2) &
                               + ( rhogvz_vm_pl(ij  ) + rhogvz_vm_pl(ijp1) ) * cinterp_TN_pl(ij  ,l,2,3) &
                               - ( rhogvx_vm_pl(ijp1) + rhogvx_vm_pl(n   ) ) * cinterp_TN_pl(ijp1,l,1,1) &
                               - ( rhogvy_vm_pl(ijp1) + rhogvy_vm_pl(n   ) ) * cinterp_TN_pl(ijp1,l,1,2) &
                               - ( rhogvz_vm_pl(ijp1) + rhogvz_vm_pl(n   ) ) * cinterp_TN_pl(ijp1,l,1,3) &
                             ) * 0.5_RP * cinterp_TRA_pl(ij,l) &
                           + sclt_rhogw_pl
             enddo

             ddivdx_pl(:,k,l) = 0.0_RP
             ddivdy_pl(:,k,l) = 0.0_RP
             ddivdz_pl(:,k,l) = 0.0_RP

             do v = ADM_gmin_pl, ADM_gmax_pl
                ij   = v
                ijm1 = v - 1
                if( ijm1 < ADM_gmin_pl ) ijm1 = ADM_gmax_pl ! cyclic condition

                ddivdx_pl(n,k,l) = ddivdx_pl(n,k,l) + ( sclt_pl(ijm1) + sclt_pl(ij) ) * cinterp_HN_pl(ij,l,1)
                ddivdy_pl(n,k,l) = ddivdy_pl(n,k,l) + ( sclt_pl(ijm1) + sclt_pl(ij) ) * cinterp_HN_pl(ij,l,2)
                ddivdz_pl(n,k,l) = ddivdz_pl(n,k,l) + ( sclt_pl(ijm1) + sclt_pl(ij) ) * cinterp_HN_pl(ij,l,3)
             enddo

             ddivdx_pl(n,k,l) = ddivdx_pl(n,k,l) * 0.5_RP * cinterp_PRA_pl(n,l)
             ddivdy_pl(n,k,l) = ddivdy_pl(n,k,l) * 0.5_RP * cinterp_PRA_pl(n,l)
             ddivdz_pl(n,k,l) = ddivdz_pl(n,k,l) * 0.5_RP * cinterp_PRA_pl(n,l)
          enddo

       enddo
    else
       ddivdx_pl(:,:,:) = 0.0_RP
       ddivdy_pl(:,:,:) = 0.0_RP
       ddivdz_pl(:,:,:) = 0.0_RP
    endif

    call DEBUG_rapend('OPRT3D_divdamp')

    return
  end subroutine OPRT3D_divdamp

  !-----------------------------------------------------------------------------
  integer function suf(i,j)
    use mod_adm, only: &
       ADM_gall_1d
    implicit none

    integer :: i, j
    !---------------------------------------------------------------------------

    suf = ADM_gall_1d * (j-1) + i

  end function suf

end module mod_oprt3d
!-------------------------------------------------------------------------------

max_stars_repo_name: hisashiyashiro/nicamdc

max_issues_repo_name: nicamdc-dev/nicamdc

max_forks_repo_name: nicamdc-dev/nicamdc


text:
Rocklin is a suburban town located about an hour east of Davis along Interstate 80. Rocklin has grown tremendously over the past 10 years. It is home to Sierra College, which is known for its attractive campus and nature trails. A number of Sierra College students end up transferring to UC Davis.

Rocklin doesnt have much in terms of largescale shopping.  This is mostly due to its residents complaints of adding traffic or changing the dynamic of the community.  To most of the teenage population and rational adults, the traffic complaint is idiotic because all the businesses that are rejected in Rocklin, set up in Roseville.  They are typically found very close to if not on the dividing line between Rocklin and Roseville.  The most well known being the Roseville Galleria.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
subroutine setgrid

! Initialize vertical grid

use vars
use params

implicit none

real latit, long
integer i,j,it,jt,k, kmax

if(nrestart.eq.0) then

 if(doconstdz) then

    z(1) = 0.5*dz
    do k=2,nz
     z(k)=z(k-1)+dz
    end do

 else

    open(8,file='./'//trim(case)//'/grd',status='old',form='formatted')
    do k=1,nz
      read(8,fmt=*,end=111) z(k)
      kmax=k
    end do
    goto 222
111 do k=kmax+1,nz
     z(k)=z(k-1)+(z(k-1)-z(k-2))
    end do
222 continue
    close (8)

 end if

end if

if(.not.doconstdz) dz = 0.5*(z(1)+z(2))

do k=2,nzm
   adzw(k) = (z(k)-z(k-1))/dz
end do
adzw(1) = 1.
adzw(nz) = adzw(nzm)
adz(1) = 1.
do k=2,nzm-1
   adz(k) = 0.5*(z(k+1)-z(k-1))/dz
end do
adz(nzm) = adzw(nzm)
zi(1) = 0.
do k=2,nz
   zi(k) = zi(k-1) + adz(k-1)*dz
end do

do k=1,nzm
  gamaz(k)=ggr/cp*z(k)
end do

if(dofplane) then

  if(fcor.eq.-999.) fcor = 4*pi/86400.*sin(latitude0*pi/180.)
  fcorz = 0.
  if(docoriolisz) fcorz =  sqrt(4.*(2*pi/(3600.*24.))**2-fcor**2)
  fcory(:) = fcor
  fcorzy(:) = fcorz
  longitude(:,:) = longitude0

else

 call task_rank_to_index(rank,it,jt)

 do j=0,ny
     latit=latitude0+dy*(j+jt-(ny_gl+YES3D-1)/2-1)*2.5e-8*360.
     fcory(j)= 4.*pi/86400.*sin(latit*pi/180.)
     fcorzy(j) = 0.
     if(j.ne.0.and.docoriolisz) fcorzy(j) = sqrt(4.*(2*pi/(3600.*24.))**2-fcory(j)**2)
 end do

end if ! dofplane

if (doradlat) then
 call task_rank_to_index(rank,it,jt)
 do j=1,ny
     latitude(:,j) = latitude0+dy*(j+jt-(ny_gl+YES3D-1)/2-1)*2.5e-8*360.
  end do
else
  latitude(:,:) = latitude0
end if

if (doradlon) then
 call task_rank_to_index(rank,it,jt)
 do i=1,nx
     longitude(i,:) = longitude0+dx/cos(latitude0*pi/180.)* &
                            (i+it-nx_gl/2-1)*2.5e-8*360.
  end do
else
  longitude(:,:) = longitude0
end if

end

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


text:
      subroutine sumrnf
c*******************************************************************
c                                                                  *
c   This subroutine is called from SR STMGET to sum up the number  *
c   of rainfall producing events and rainfall amounts during the   *
c   simulation period.                                             *
c                                                                  *
c*******************************************************************
c
      include 'pmxelm.inc'
      include 'pmxpln.inc'
      include 'pmxhil.inc'
c
c*******************************************************************
c                                                                  *
c   Common Blocks                                                  *
c                                                                  *
c*******************************************************************
c
      include 'cavloss.inc'
c
      include 'chydrol.inc'
c
      include 'cstmflg.inc'
c
      include 'cstruc.inc'
      include 'csumout.inc'
c
c*******************************************************************
c                                                                  *
c   sumout variables updated                                       *
c       nraint,nrainy,nrainm(13),traint,trainy,trainm(13)          *
c                                                                  *
c*******************************************************************
c
      nraint = nraint + 1
      nrainy = nrainy + 1
      nrainm = nrainm + 1
c
      traint = traint + (prcp*1000.0)
      trainy = trainy + (prcp*1000.0)
      trainm = trainm + (prcp*1000.0)
c
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
! :::::::::::::::::::::::::::: FILPATCH :::::::::::::::::::::::::;
!
!  fill the portion of valbig from rows  nrowst
!                             and  cols  ncolst
!  the patch can also be described by the corners (xlp,ybp) by (xrp,ytp).
!  vals are needed at time t, and level level,
!
!  first fill with  values obtainable from the level level
!  grids. if any left unfilled, then enlarge remaining rectangle of
!  unfilled values by 1 (for later linear interp), and recusively
!  obtain the remaining values from  coarser levels.
!
! :::::::::::::::::::::::::::::::::::::::;:::::::::::::::::::::::;
!
recursive subroutine filrecur(level,nvar,valbig,aux,naux,t,mitot,mjtot, &
                              nrowst,ncolst,ilo,ihi,jlo,jhi,patchOnly)

    use amr_module, only: hxposs, hyposs, xlower, ylower, xupper, yupper
    use amr_module, only: outunit, nghost, xperdom, yperdom, spheredom
    use amr_module, only: iregsz, jregsz, intratx, intraty, NEEDS_TO_BE_SET

    !for setaux timing
    use amr_module, only: timeSetaux, timeSetauxCPU

    implicit none

    ! Input
    integer, intent(in) :: level, nvar, naux, mitot, mjtot, nrowst, ncolst
    integer, intent(in) :: ilo,ihi,jlo,jhi
    real(kind=8), intent(in) :: t
    logical  :: patchOnly

    ! Output
    real(kind=8), intent(in out) :: valbig(nvar,mitot,mjtot)
    real(kind=8), intent(in out) :: aux(naux,mitot,mjtot)

    ! Local storage
    integer :: i_fine, j_fine, i_coarse, j_coarse, n, k
    integer  :: iplo, iphi, jplo, jphi
    integer :: mitot_patch, mjtot_patch, mitot_coarse, mjtot_coarse, nghost_patch, lencrse
    integer :: refinement_ratio_x, refinement_ratio_y
    integer :: unset_indices(4)
    real(kind=8) :: dx_fine, dy_fine, dx_coarse, dy_coarse
    real(kind=8) :: xlow_coarse,ylow_coarse, xlow_fine, ylow_fine, xhi_fine,yhi_fine

    !for timing
    integer :: clock_start, clock_finish, clock_rate
    real(kind=8) :: cpu_start, cpu_finish

    ! Interpolation variables
    real(kind=8) :: eta1, eta2, valp10, valm10, valc, valp01, valm01, dupc, dumc
    real(kind=8) :: ducc, du, fu, dvpc, dvmc, dvcc, dv, fv, valint

    ! Cell set tracking
    logical :: set
    integer(kind=1) :: flaguse(ihi-ilo+1, jhi-jlo+1)

    ! Scratch storage
    !  use stack-based scratch arrays instead of alloc, since dont really
    !  need to save beyond these routines, and to allow dynamic memory resizing
    !
    !     use 1d scratch arrays that are potentially the same size as
    !     current grid, since may not coarsen.
    !     need to make it 1d instead of 2 and do own indexing, since
    !     when pass it in to subroutines they treat it as having dierent
    !     dimensions than the max size need to allocate here

    !--      dimension valcrse((ihi-ilo+2)*(jhi-jlo+2)*nvar)  ! NB this is a 1D array
    !--      dimension auxcrse((ihi-ilo+2)*(jhi-jlo+2)*naux)  ! the +2 is to expand on coarse grid to enclose fine
    ! ### turns out you need 3 rows, forget offset of 1 plus one on each side
    ! the +3 is to expand on coarse grid to enclose fine
    real(kind=8) :: valcrse((ihi-ilo+3) * (jhi-jlo+3) * nvar)
    real(kind=8) :: auxcrse((ihi-ilo+3) * (jhi-jlo+3) * naux)

    ! We begin by filling values for grids at level level. If all values can be
    ! filled in this way, we return;
    mitot_patch = ihi-ilo + 1 ! nrowp
    mjtot_patch = jhi-jlo + 1

    dx_fine     = hxposs(level)
    dy_fine     = hyposs(level)

    ! Coordinates of edges of patch (xlp,xrp,ybp,ytp)
    xlow_fine = xlower + ilo * dx_fine
    ylow_fine = ylower + jlo * dy_fine

    ! Fill in the patch as much as possible using values at this level
    call intfil(valbig,mitot,mjtot,t,flaguse,nrowst,ncolst, ilo,  &
                ihi,jlo,jhi,level,nvar,naux)

    ! Trimbd returns set = true if all of the entries are filled (=1.).
    ! set = false, otherwise. If set = true, then no other levels are
    ! are required to interpolate, and we return.
    !
    ! Note that the used array is filled entirely in intfil, i.e. the
    ! marking done there also takes  into account the points filled by
    ! the boundary conditions. bc2amr will be called later, after all 4
    ! boundary pieces filled.
    call trimbd(flaguse,mitot_patch,mjtot_patch,set,unset_indices)
    ! il,ir,jb,jt = unset_indices(4)

    ! If set is .true. then all cells have been set and we can skip to setting
    ! the remaining boundary cells.  If it is .false. we need to interpolate
    ! some values from coarser levels, possibly calling this routine
    ! recursively.
    if (.not.set) then

        ! Error check
        if (level == 1) then
            write(outunit,*)" error in filrecur - level 1 not set"
            write(outunit,'("start at row: ",i4," col ",i4)') nrowst,ncolst
            print *," error in filrecur - level 1 not set"
            print *," should not need more recursion "
            print *," to set patch boundaries"
            print '("start at row: ",i4," col ",i4)', nrowst,ncolst
            stop
        endif

        ! We begin by initializing the level level arrays, so that we can use
        ! purely recursive formulation for interpolating.
        dx_coarse  = hxposs(level - 1)
        dy_coarse  = hyposs(level - 1)

        ! Adjust unset_indices to account for the patch
        ! isl, isr, jsb, jst
        unset_indices(1) = unset_indices(1) + ilo - 1
        unset_indices(2) = unset_indices(2) + ilo - 1
        unset_indices(3) = unset_indices(3) + jlo - 1
        unset_indices(4) = unset_indices(4) + jlo - 1

        ! Coarsened geometry
        refinement_ratio_x = intratx(level - 1)
        refinement_ratio_y = intraty(level - 1)

        ! New patch rectangle (after we have partially filled it in) but in the
        ! coarse patches [iplo,iphi,jplo,jphi]

        iplo = (unset_indices(1) - refinement_ratio_x + nghost * refinement_ratio_x) &
                                                 / refinement_ratio_x - nghost
        iphi = (unset_indices(2) + refinement_ratio_x) / refinement_ratio_x
        jplo = (unset_indices(3) - refinement_ratio_y + nghost * refinement_ratio_y) &
                                                / refinement_ratio_y - nghost
        jphi = (unset_indices(4) + refinement_ratio_y) / refinement_ratio_y

        xlow_coarse = xlower + iplo * dx_coarse
        ylow_coarse = ylower + jplo * dy_coarse

        ! Coarse grid number of spatial points (nrowc,ncolc)
        mitot_coarse   =  iphi - iplo + 1
        mjtot_coarse   =  jphi - jplo + 1

        ! Check to make sure we created big enough scratch arrays
        if (mitot_coarse > ihi - ilo + 3 .or. &
            mjtot_coarse > jhi - jlo + 3) then

            print *," did not make big enough work space in filrecur "
            print *," need coarse space with mitot_coarse,mjtot_coarse ",mitot_coarse,mjtot_coarse
            print *," made space for ilo,ihi,jlo,jhi ",ilo,ihi,jlo,jhi
            stop
        endif

        ! Set the aux array values for the coarse grid, this could be done
        ! instead in intfil using possibly already available bathy data from the
        ! grids
        if (naux > 0) then
            nghost_patch = 0
            lencrse = (ihi-ilo+3)*(jhi-jlo+3)*naux ! set 1 component, not all naux
            do k = 1, lencrse, naux
              auxcrse(k) = NEEDS_TO_BE_SET ! new system checks initialization before setting aux vals
            end do
            call system_clock(clock_start,clock_rate)
            call cpu_time(cpu_start)
            call setaux(nghost_patch, mitot_coarse,mjtot_coarse,  &
                       xlow_coarse, ylow_coarse,            &
                        dx_coarse,dy_coarse,naux,auxcrse)
            call system_clock(clock_finish,clock_rate)
            call cpu_time(cpu_finish)
            timeSetaux = timeSetaux + clock_finish - clock_start
            timeSetauxCPU = timeSetauxCPU + cpu_finish - cpu_start
        endif

        ! Fill in the edges of the coarse grid. for recursive calls, patch indices and
        ! 'coarse grid' indices are the same (no actual coarse grid here, so cant use mptr
        ! must pass indices. patchOnly argument  is thus true
        if ((xperdom .or. (yperdom .or. spheredom)) .and. sticksout(iplo,iphi,jplo,jphi)) then
            call prefilrecur(level-1,nvar,valcrse,auxcrse,naux,t,mitot_coarse,mjtot_coarse,1,1,   &
                             iplo,iphi,jplo,jphi,iplo,iphi,jplo,jphi,.true.)
        else
            call filrecur(level-1,nvar,valcrse,auxcrse,naux,t,mitot_coarse,mjtot_coarse,1,1,   &
                          iplo,iphi,jplo,jphi,.true.)
        endif

        do i_fine = 1,mitot_patch
            i_coarse = 2 + (i_fine - (unset_indices(1) - ilo) - 1) / refinement_ratio_x
            eta1 = (-0.5d0 + real(mod(i_fine - 1, refinement_ratio_x),kind=8)) &
                                / real(refinement_ratio_x,kind=8)

            do j_fine  = 1,mjtot_patch
                j_coarse = 2 + (j_fine - (unset_indices(3) - jlo) - 1) / refinement_ratio_y
                eta2 = (-0.5d0 + real(mod(j_fine - 1, refinement_ratio_y),kind=8)) &
                                    / real(refinement_ratio_y,kind=8)

                if (flaguse(i_fine,j_fine) == 0) then

                    do n=1,nvar
                       ! write(*,*)n,i_coarse+1,j_coarse,coarse_index(n,i_coarse + 1,j_coarse)
                        valp10 = valcrse(coarse_index(n,i_coarse + 1,j_coarse))
                        valm10 = valcrse(coarse_index(n,i_coarse - 1,j_coarse))
                        valc   = valcrse(coarse_index(n,i_coarse    ,j_coarse))
                        valp01 = valcrse(coarse_index(n,i_coarse    ,j_coarse + 1))
                        valm01 = valcrse(coarse_index(n,i_coarse    ,j_coarse - 1))

                        dupc = valp10 - valc
                        dumc = valc   - valm10
                        ducc = valp10 - valm10
                        du   = min(abs(dupc), abs(dumc))
                        du   = min(2.d0 * du, 0.5d0 * abs(ducc))
                        fu = max(0.d0, sign(1.d0, dupc * dumc))

                        dvpc = valp01 - valc
                        dvmc = valc   - valm01
                        dvcc = valp01 - valm01
                        dv   = min(abs(dvpc), abs(dvmc))
                        dv   = min(2.d0 * dv, 0.5d0 * abs(dvcc))
                        fv = max(0.d0,sign(1.d0, dvpc * dvmc))

                        valint = valc + eta1 * du * sign(1.d0, ducc) * fu &
                                      + eta2 * dv * sign(1.d0, dvcc) * fv

                        valbig(n,i_fine+nrowst-1,j_fine+ncolst-1) = valint

                    end do

                endif
            end do
        end do
    end if

    !  set bcs, whether or not recursive calls needed. set any part of patch that stuck out
    xhi_fine = xlower + (ihi + 1) * dx_fine
    yhi_fine = ylower + (jhi + 1) * dy_fine
    ! only call if a small coarser recursive patch
    ! otherwise whole grid bcs done from bound
    if (patchOnly) then
       call bc2amr(valbig,aux,mitot,mjtot,nvar,naux,dx_fine,dy_fine,level,t,xlow_fine,xhi_fine, &
                   ylow_fine,yhi_fine,xlower,ylower,xupper,yupper,xperdom,yperdom,spheredom)
    endif

contains

    integer pure function coarse_index(n,i,j)
        implicit none
        integer, intent(in) :: n, i, j
        coarse_index = n + nvar*(i-1)+nvar*mitot_coarse*(j-1)
    end function coarse_index

    logical pure function sticksout(iplo,iphi,jplo,jphi)
        implicit none
        integer, intent(in) :: iplo,iphi,jplo,jphi
        sticksout = (iplo < 0 .or. jplo < 0 .or. &
                     iphi >= iregsz(level - 1) .or. jphi >= jregsz(level - 1))
    end function sticksout

end subroutine filrecur

max_stars_repo_name: ian-r-rose/visualization

max_issues_repo_name: ian-r-rose/visualization

max_forks_repo_name: ian-r-rose/visualization


text:
Subroutine genei
  Implicit Double Precision (A-H, O-Z)
  Parameter (maxptn=400001)
  Common /para1/mul
  Common /para2/xmp, xmu, alpha, rscut2, cutof2
  Common /para3/nsevt, nevnt, nsbrun, ievt, isbrun
  Common /para5/iconfg, iordsc
  Common /prec1/gx0(maxptn), gy0(maxptn), gz0(maxptn), ft0(maxptn), px0(maxptn), py0(maxptn), pz0(maxptn), e0(maxptn), xmass0(maxptn), ityp0(maxptn)
  Common /prec5/eta(maxptn), rap(maxptn), tau(maxptn)
  Common /lor/enenew, pxnew, pynew, pznew
  Common /rndm3/iseedp
  Save
  External ran1
  iseed = iseedp
  incmul = 4000
  temp = 0.5D0
  etamin = -5D0
  etamax = 5D0
  r0 = 5D0
  tau0 = 0.1D0
  deta = etamax - etamin
  Do i = mul + 1, mul + incmul
     ityp0(i) = 21
     xmass0(i) = xmp
     Call energy(e, temp)
     Call momntm(px, py, pz, e)
     e = dsqrt(e**2+xmp**2)
     If (iconfg<=3) Then
        eta(i) = etamin + deta*ran1(iseed)
        bex = 0D0
        bey = 0D0
        bez = -tanh(eta(i))
        Call lorenz(e, px, py, pz, bex, bey, bez)
        px0(i) = pxnew
        py0(i) = pynew
        pz0(i) = pznew
        e0(i) = enenew
     Else
        px0(i) = px
        py0(i) = py
        pz0(i) = pz
        e0(i) = e
     End If
  End Do
  Do i = mul + 1, mul + incmul
     If (iconfg<=3) Then
        gz0(i) = tau0*sinh(eta(i))
        ft0(i) = tau0*cosh(eta(i))
        If (iconfg==1) Then
           Call posit1(x, y, r0)
           gx0(i) = x + px0(i)*ft0(i)/e0(i)
           gy0(i) = y + py0(i)*ft0(i)/e0(i)
        Else If (iconfg==2 .Or. iconfg==3) Then
           Call posit2(x, y)
           gx0(i) = x
           gy0(i) = y
        End If
     Else
        ft0(i) = 0D0
        Call posit3(x, y, z)
        gx0(i) = x
        gy0(i) = y
        gz0(i) = z
     End If
  End Do
  mul = mul + incmul
  If (mul>=maxptn .Or. mul==0) Then
     Print *, 'event', ievt, 'has', mul, 'number of gluon', 'adjusting counting is necessary'
     Stop 'adarr'
  End If
  Return
End Subroutine genei

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
      SUBROUTINE DMOMNT
C
C           'DMOMNT' CALCULATES FLEXIBLE ELEMENT ROOT BENDING MOMENTS
C            AND ELEMENT ROOT TENSIONS
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      COMMON/CANTNA/ A(10,3),ADOT(10,3),B(10,3),BDOT(10,3),DIN(10,3),
     .               DINDOT(10,3),DOUT(10,3),DOUTDT(10,3)
C
      COMMON/COMALP/ SZ02(10),SZ03(10),SZ04(10),SZ12(3,10),SZ13(3,10),
     .               SZ14(3,10),SZ15(3,10),SZ16(3,10),SZ21(9,10),
     .               SZ22(9,10),SZ23(9,10),SZ25(9,10),
     .               SZ26(9,10),SZ27(9,10),SZ28(9,10),SZ31(27,10),
     .               SZ32(27,10),SZ33(27,10),SZ34(27,10),SZ35(27,10),
     .               SZ41(81,10),SZ42(81,10),SZ43(81,10)
C
      COMMON/CSOLAR/ SAO(10),SKA(9),SKB(9),SKOA(10,3),SKOB(10,3),
     .               STMK(10),SKAA(10,9),SKBB(10,9)
C
      COMMON/DEBUG2/ IOUT,JOUT,KLUGE
C
      COMMON/DMMNT1/ ZKBM(6),EMAK(10),EMBK(10),ARTETA(3),CMTORK(3),
     .               ITORK,IBENDM,ITENSE,ITNS1
C
      COMMON/IPOOL1/ IGRAV,IDAMP,IK,K1,ITIM,IAB,IAPS,IBB,IBPS,NK(10),
     .               LK(10),LLK(10)
C
      COMMON/MRANG/ ADDOT(10,3),BDDOT(10,3),DOUTDD(10,3),DINDD(10,3)
C
      COMMON/PRCOM/ STORE(10,30),ILINE,ICOL,ICNT,IHD
C
      COMMON/RPOOL1/ RHOK(10),TIME,SA(3,3),FM1(3,3),ZLK(10),OMEG(3),
     .               ZLKP(10),ZLKDP(10),CMAT(3,3),GBAR(3,3),YBCM(3),
     .               ZBZK(3,10),FCM(3,3),DTO,PHID,PHI
C
      COMMON/RPOOL3/ ZMS,YIZM(3,2)
C
      COMMON/RPOOL5/ CKMAT(3,3,10),FM2(3,3)
C
      COMMON/RPOOL8/ SZ01(10),SZ11(3,10),SZ24(9,10)
C
      COMMON/SPRESX/ SPRES(10,3)
C
      COMMON/VARBLS/ DEPEND(150),DERIV(150)
C
      COMMON/TENSON/TSSO(10)
C
      DIMENSION VECTS(3),VECTT(3),VDOT(3),OMGD(3),OMEGDT(3),W(3),
     .          XMOM(3,3),DUM(3,3),RES(3,3),WDOT(3),VECTD(3),FM1T(3,3),
     .          FM2T(3,3),AT(3,3),PRD(3,3),SAP(3),RES2(3,3),
     .          XI(3),XID(3),SAT(3,3)
C
      NAMELIST/OOWWD/OMGD,OMEGDT,W,WDOT/D2/VECTS,VECTD,TDZ,VECTT,
     .               YBCM,ARTETA
C
C
      IF(IOUT.NE.1) WRITE(6,570)
      DO 10 L=1,10
      EMAK(L)=0.0D0
   10 EMBK(L)=0.0D0
      DO 11 J2=1,3
      IF(IOUT.NE.1) WRITE(6,600) (GBAR(J2,J1),J1=1,3)
      VECTT(J2)=0.0D0
      VDOT(J2)=0.0D0
      OMGD(J2)=0.0D0
      OMEGDT(J2)=0.0D0
      W(J2)=0.0D0
   11 WDOT(J2)=0.0D0
      IF(IBENDM.EQ.0) GO TO 35
      DO 50 K=1,IK
      LX=(LLK(K)-1)*3
      N=NK(K)
      IF(N.EQ.0) GO TO 50
   15 EIA=0.5D0*SKAA(K,1)
      EIB=0.5D0*SKBB(K,1)
      FLK=ZLK(K)
      DO 30 I=1,N
      ANT1=FUNA(K,K1,I)
      BNT1=FUNB(K,K1,I)
C
      AN0=SKOA(K,I)
      BN0=SKOB(K,I)
      IF(IOUT.NE.1) PRINT 25,AN0,BN0,ANT1,BNT1
   25 FORMAT(' AN0=',G20.12, ' BN0=',G20.12, ' ANT1=',G20.12,
     .  ' BNT1=',G20.12//)
C
      EMAK(K)=EMAK(K)+(ANT1-AN0)*ZKBM(LX+I)
      EMBK(K)=EMBK(K)+(BNT1-BN0)*ZKBM(LX+I)
   30 CONTINUE
      FLK2=FLK*FLK
      EMAK(K)=EIA/FLK2*EMAK(K)
      EMBK(K)=EIB/FLK2*EMBK(K)
   50 CONTINUE
   35 DO 150 K=1,IK
      IF(ITENSE.EQ.0) GO TO 533
      TDZ=RHOK(K)*ZLK(K)*2.0D0*SZ02(K)
C
      CALL XIMMT(K,XI,XID)
C
      DO 140 J3=1,3
      SAP(J3)=0.0D0
  140 CONTINUE
      DO 145 J2=1,3
      SAP(J2)=RHOK(K)*ZLK(K)*XI(J2)
      IF(IOUT.NE.1) PRINT 143,J2,SAP(J2)
  143 FORMAT(' VECTOR SAP ',I1, 2X, G20.12)
  145 CONTINUE
      SAPD=RHOK(K)*ZLK(K)*XID(1)
      IF(IDAMP.EQ.0.OR.K.GT.K1) GO TO 40
      IF(IDAMP.EQ.1.AND.K.LE.K1) GO TO 220
      PRINT 200,IDAMP,K,K1
  200 FORMAT ('0   DAMPER VALUE AND/OR MODES QUESTIONABLE, IDAMP = ',
     . I3, ' K = ', I3, ' K1 = ', I3)
      RETURN
  220 CONTINUE
      ICZ=1
      VECTS(1)=0.0D0
      VECTS(2)=DEPEND(11)
      VECTS(3)=0.0D0
      VECTD(1)=0.0D0
      VECTD(2)=DERIV(11)
      VECTD(3)=0.0D0
      DO 225 J4=1,3
      DO 225 J5=1,3
      FM1T(J4,J5)=FM1(J5,J4)
  225 CONTINUE
      GO TO 90
   40 CONTINUE
      ICZ=2
      DO 70 M1=1,3
      VECTS(M1)=0.0D0
   45 VECTD(M1)=0.0D0
      DO 70 L1=1,3
   70 FM1T(M1,L1)=FM2(L1,M1)
C
   90 CONTINUE
C
      DO 80 M1=1,3
      DO 80 L1=1,3
      XMOM(M1,L1)=CKMAT(L1,M1,K)
 80   CONTINUE
      CALL MPYMAT(XMOM,FM1T,DUM,1,1,RES,DUM)
C
      DO 250 M1=1,3
      DO 211 JX=1,3
      DO 211 JY=1,3
      RES2(JY,JX)=0.0D0
  211 CONTINUE
C
      OMGD(M1)=RES(M1,1)*DEPEND(7)+RES(M1,2)*DEPEND(8)
     .         +RES(M1,3)*DEPEND(9)
      OMEGDT(M1)=RES(M1,1)*DERIV(7)+RES(M1,2)*DERIV(8)
     .           +RES(M1,3)*DERIV(9)
C
      W(M1)=XMOM(M1,1)*VECTS(1)+XMOM(M1,2)*VECTS(2)
     .      +XMOM(M1,3)*VECTS(3)
      WDOT(M1)=XMOM(M1,1)*VECTD(1)+XMOM(M1,2)*VECTD(2)
     .         +XMOM(M1,3)*VECTD(3)
C
C          COMPUTE   1ST  TERM  -  ART
C
      VECTT(M1)=(RES(M1,1)*ARTETA(1)+RES(M1,2)*ARTETA(2)
     .          +RES(M1,3)*ARTETA(3))*TDZ
  250   CONTINUE
      SMULDT=2.0D0*RHOK(K)*ZLK(K)*((OMGD(2)+W(2))*XID(3)-(OMGD(3)
     .      +W(3))*XID(2))
       IF (IOUT .NE. 1) PRINT 147,SMULDT
  147  FORMAT(' SMULDT=', G20.12)
      DO  340  J2 =1, 3
  340 VECTS (J2) = 0.D0
       IF (IOUT .NE. 1) WRITE(6,OOWWD)
       T9=2.D0* SZ02 (K)
      DO  350   J2 = 1, 3
      VECTS(J2)=(XMOM(J2,1)*ZBZK(1,K)+XMOM(J2,2)*ZBZK(2,K)+XMOM(J2,3)*
     .          ZBZK(3,K))*T9+XI(J2)
  350  CONTINUE
      IF (IOUT .NE. 1) WRITE (6,D2)
      DO   300   J1 =  1, 3
      DO   300   J2 =  1, 3
C
       SAT (J1,J2) =SA (J2,J1)
  300 RES2  (J2, J1) =  0.D0
      DO 320 J1=1,3
      RES2(J1,1)=RES(J1,1)*YIZM(1,ICZ)+RES(J1,2)*YIZM(2,ICZ)+
     .           RES(J1,3)*YIZM(3,ICZ)
  320 CONTINUE
C
      SAPD1=TDZ*(-(OMGD(2)**2+OMGD(3)**2)*RES2(1,1)+(-OMEGDT(3)+
     .      OMGD(1)*OMGD(2))*RES2(2,1)+(OMEGDT(2)+OMGD(1)*OMGD(3))
     .      *RES2(3,1))
C
C
C
      SAP1=RHOK(K)*ZLK(K)*((-(OMGD(2)+W(2))**2-(OMGD(3)+W(3))**2)*
     .     VECTS(1)+(-OMEGDT(3)-WDOT(3)+OMGD(2)*(OMGD(1)+W(1))+(OMGD(2)
     .     +W(2))*W(1))*VECTS(2)+(OMEGDT(2)+WDOT(2)+OMGD(3)*(OMGD(1)
     .     +W(1))+(OMGD(3)+W(3))*W(1))*VECTS(3))
C
C
C
      TERMS=SAPD1+SAP1+SAPD+SMULDT+VECTT(1)+XID(1)
      IF(IOUT.NE.1)PRINT 999,TERMS
  999  FORMAT('0 INERTIA FORCES',G20.12)
       IF (IOUT .NE. 1) PRINT 610,SAPD1     ,SAP1     ,SAPD
  610  FORMAT(' SAPD1 ',G20.12, ' SAP1 ',G20.12,' SAPD ',G20.12)
C
C          GRAVITY   GRADIENT
C
      CALL  MPYMAT (RES, SAT,  DUM,  1, 1,  RES2, DUM)
      CALL MPYMAT (RES2, GBAR, DUM,  1, 1,  PRD, DUM)
      CALL MPYMAT (PRD,  SA,   DUM,  1, 1,  AT, DUM)
C
       IF (IOUT .EQ. 1)GO TO 407
       WRITE(6,403)
  403  FORMAT('0 AT'//)
       DO 406 I1=1, 3
       WRITE (6,405)(AT(I1,K2),K2=1,3)
  405  FORMAT('   ', 3G20.12)
  406  CONTINUE
  407  CONTINUE
      DO  410  J1 = 1, 3
      VECTS (J1) =  0.D0
      DO 410 J2 =1,3
      PRD (J1,J2) =  0.D0
      FM2T (J1,J2) =  FM1T(J2,J1)
  410 CONTINUE
C
C
C
      DO   420   J2 = 1, 3
      VECTS(J2)=((FM2T(J2,1)*ZBZK(1,K)+FM2T(J2,2)*ZBZK(2,K)+FM2T(J2,3)
     .          *ZBZK(3,K))+YIZM(J2,ICZ))*TDZ
  420 CONTINUE
C
C     COMBINE AT MATRIX X VECTS VECTOR
      DO  440  J1 =  1, 3
      PRD(J1,1)=AT(J1,1)*VECTS(1)+AT(J1,2)*VECTS(2)+AT(J1,3)*VECTS(3)
  440 CONTINUE
      IF (IOUT .EQ. 1) GO TO 450
       WRITE (6,445)
  445 FORMAT('0 PRD - 1ST EQUATION'//)
      DO 447 L7=1,3
       WRITE (6,405)(PRD(L7,L8),L8=1,3)
  447  CONTINUE
  450  CONTINUE
C
C       COMPUTE   2D  TERM  GRAVITY GRADIENT
C
      CALL  MPYMAT (AT, FM1T, DUM,  1, 2,  RES2, DUM)
      CALL  MPYMAT (RES2, XMOM, DUM, 1, 2, FM1T, DUM)
C
C   FINISH 2D LINE OF GRAVITY GRADIENT FORCES EQUATION
C
      DO   460   J1 =  1, 3
  460 FM2T(J1,1)=FM1T(J1,1)*SAP(1)+FM1T(J1,2)*SAP(2)+FM1T(J1,3)*SAP(3)
       IF (IOUT .EQ. 1) GO TO 470
       WRITE (6,471)
  471  FORMAT('0  FM2T - 2D EQUATION'//)
       DO 473 L7=1,3
  473  WRITE (6,405)(FM2T(L7,L8),L8=1,3)
  470  CONTINUE
      DO    490  J1 =  1, 3
  490 FM1T(J1,1)=(AT(J1,1)*YBCM(1)+AT(J1,2)*YBCM(2)+
     .           AT(J1,3)*YBCM(3))*(TDZ/ZMS)
       IF (IOUT .EQ. 1) GO TO 498
       WRITE (6, 492)
  492  FORMAT('0 FM1T - 3RD  EQUATION'//)
       DO 495 L7=1,3
       WRITE (6,405)(FM1T(L7,L8),L8=1,3)
  495  CONTINUE
  498  CONTINUE
C
C       GRAVITY   GRADIENT   FORMULA
C
      DO 510 J1=1,3
      RES2(J1,1)=PRD(J1,1)+FM2T(J1,1)-FM1T(J1,1)
  510 CONTINUE
       IF (IOUT .EQ. 1) GO TO 530
       WRITE (6,540)(RES2(J1,1), J1=1,3), SPRES(K,1)
  540  FORMAT(' GRAV ',3G20.12, ' SOLAR ',G20.12/)
  520  CONTINUE
  530 CONTINUE
C
       TSS1=-TERMS +RES2(1,1)+SPRES(K,1)
      TSSO(K)=TSS1
 533  CONTINUE
  150 CONTINUE
  570  FORMAT ('0  GBAR ARRAY'//)
  600  FORMAT('  ',3G20.12)
      RETURN
      END


text:
Faith and Social Justice is a multidenominational, Biblebased community that seeks to integrate its faith, politics, and lifestyle by studying, discussing, and taking action on social justice issues both locally and abroad. The group meets regularly for lively discussions of Scripture and contemporary books.   Particpants believe they have been called to care for the earth and to reach out to all people with love, nonviolence, and justice.

During the school year, the group meets monthly (often the 3rd Friday of each month from 79pm) in peoples home.  Email Terry to be added to the groups distribution list and receive 23 emails per month with news, upcoming events, and summaries of past meetings.

20070613 13:18:00 nbsp Hi!  Where do you all meet?  and what times generally?  or rather, how would one go about coming to a meeting to find out more? Users/RocksandDirt


text:
C->>> ---------------------------------------> ems_ca_iz_lg_bs_vr_st <<<
      subroutine ems_ca_iz_lg_bs_vr_st(ds, is)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'EMSMMGR.INC'
      include 'EMSMEM.INC'
      include 'EMSP.INC'
      include 'ICTVR.INC'
      integer is(0:is_n_en_m1)
      double precision ds(0:ds_n_en_m1)
      call ems_iz_lg_bs_vr_st(
     &     ds(p_lbc), ds(p_cbp), ds(p_ubc),
     &     is(p_st), ds(p_pr_act), ds(p_du_act), ds(p_scl),
     &     is(p_vr_in_r), is(p_vr_in_c))
      return
      end

C->>> ------------------------------------------> ems_iz_lg_bs_vr_st <<<
c     Initialise the status vector corresponding to primal logical
c     variables being basic and primal structural variables being
c     nonbasic and set all the primal activities to zero.
c
      subroutine ems_iz_lg_bs_vr_st(
     &     lbc, cbp, ubc,
     &     st, pr_act, du_act, scl,
     &     vr_in_r, vr_in_c)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'EMSMMGR.INC'
      include 'EMSMEM.INC'
      include 'EMSP.INC'
      include 'ICTVR.INC'
      include 'RLCTVR.INC'
      include 'EMSMSG.INC'
      double precision lbc(0:mx_n_c+mx_n_r)
      double precision cbp(0:mx_n_c+mx_n_r)
      double precision ubc(0:mx_n_c+mx_n_r)
      integer st(0:mx_n_c+mx_n_r)
      integer vr_in_r(0:n_r)
      integer vr_in_c(-vr_in_c_n_sn:n_c)
      double precision pr_act(0:mx_n_c+mx_n_r)
      double precision du_act(0:mx_n_c+mx_n_r)
      double precision scl(0:mx_n_c+mx_n_r)
      integer r_n, c_n, vr_n, sv_vr_in_c_0

      do 10, c_n = 1, n_c
         vr_in_c(c_n) = c_n
         st(c_n) = non_bc_vr_bs_st
         du_act(c_n) = zero
 10   continue
c
c     vr_in_c(0) indicates the dimension for which vr_in_c was defined
c     so preserve this value whilst using vr_in_c(0) to indicate the
c     number of entries currently in vr_in_c.
c
      sv_vr_in_c_0 = vr_in_c(0)
      vr_in_c(0) = n_c
      call ems_iz_st_bd_bt(vr_in_c(0), lbc, ubc, st)
      call ems_iz_vr_st_act(
     &     vr_in_c(0), lbc, cbp, ubc, st, pr_act, scl, tl_mx_iz_pr_act)
      vr_in_c(0) = sv_vr_in_c_0
      do 20, r_n = 1, n_r
         vr_n = mx_n_c+r_n
         vr_in_r(r_n) = vr_n
         st(vr_n) = bc_vr_bs_st + r_n
         du_act(vr_n) = zero
 20   continue
      vr_in_r(0) = n_r
      call ems_iz_st_bd_bt(vr_in_r, lbc, ubc, st)
      call ems_iz_vr_st_act(
     &     vr_in_r, lbc, cbp, ubc, st, pr_act, scl, tl_mx_iz_pr_act)
c
c     Indicate that the following are not correct for the model:
c
c     vr_in_c, INVERT, basic primal activities, edge weights and
c     row-wise representation of matrix columns being priced.
c
      ml_da_st_msk = ml_da_st_msk
     &     - iand(ml_da_st_msk, ml_da_st_vr_in_c)
     &     - iand(ml_da_st_msk, ml_da_st_inv)
     &     - iand(ml_da_st_msk, ml_da_st_bc_pr_act)
     &     - iand(ml_da_st_msk, ml_da_st_ed_wt)
     &     - iand(ml_da_st_msk, ml_da_st_r_mtx)
      return
      end

C->>> ---------------------------------------------> ems_iz_st_bd_bt <<<
      subroutine ems_iz_st_bd_bt(vr_ls, lbc, ubc, st)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      integer vr_ls(0:*)
      double precision lbc(0:mx_n_c+mx_n_r)
      double precision ubc(0:mx_n_c+mx_n_r)
      integer st(0:mx_n_c+mx_n_r)
      integer ls_en_n, vr_n, vr_st

      do 10, ls_en_n = 1, vr_ls(0)
         vr_n = vr_ls(ls_en_n)
         vr_st = st(vr_n)
         if (iand(vr_st, alt_bt) .eq. 0) then
            vr_st = vr_st - iand(st(vr_n), lb_ub)
            if (lbc(vr_n) .gt. -inf) vr_st = vr_st + lb_bt
            if (ubc(vr_n) .lt.  inf) vr_st = vr_st + ub_bt
            st(vr_n) = vr_st
         endif
 10   continue
      return
      end

C->>> --------------------------------------------> ems_iz_vr_st_act <<<
c     Initialises the up, down and ifs bits of the status, and the
c     primal activity for each of a list of (nonbasic) variables.
c     Assumes that the lower and upper bound bits are set correctly.
c     Passes through the primal activities, making sure that they have
c     not been set to excessively large/small values.
c
      subroutine ems_iz_vr_st_act(
     &     vr_ls, lbc, cbp, ubc, st, pr_act, scl, usr_tl_mx_iz_pr_act)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      integer vr_ls(0:*)
      double precision lbc(0:mx_n_c+mx_n_r)
      double precision cbp(0:mx_n_c+mx_n_r)
      double precision ubc(0:mx_n_c+mx_n_r)
      double precision pr_act(0:mx_n_c+mx_n_r)
      double precision scl(0:mx_n_c+mx_n_r)
      double precision usr_tl_mx_iz_pr_act
      integer st(0:mx_n_c+mx_n_r)
      integer ls_en_n, vr_n, vr_st
      double precision scl_v
      double precision pr_act_v
      logical scl_pr_act_v

      scl_v = one
      scl_pr_act_v = iand(ml_da_st_msk, ml_da_st_scl_ml_sol) .ne. 0
      do 10, ls_en_n = 1, vr_ls(0)
         vr_n = vr_ls(ls_en_n)
         vr_st = st(vr_n)
         if (scl_pr_act_v) scl_v = scl(vr_n)
c
c     Make sure that the up, dn and ifs bits are not set
c
         vr_st = vr_st - iand(vr_st, ifs_bt)
         vr_st = vr_st - iand(vr_st, up_dn)
         if (iand(vr_st, alt_bt) .eq. 0) then
c
c     Variable is standard
c
            if (iand(vr_st, ub_bt) .eq. 0) then
               if (iand(vr_st, lb_bt) .eq. 0) then
c
c     A FR variable: Set the primal activity to zero
c
c     Set the primal activity to zero and set the up and dn bits.
c
                  pr_act(vr_n) = zero
                  vr_st = vr_st + up_dn
               else
c
c     A LB variable: Set the primal activity to the lower bound
c
                  pr_act_v = lbc(vr_n)*scl_v
                  if (pr_act_v .ge. -usr_tl_mx_iz_pr_act) then
c
c     Set the primal activity to the lower bound and set the up bit...
c
                     pr_act(vr_n) = pr_act_v
                     vr_st = vr_st + up_bt
                  else
c
c     ... unless this is an excessively low value, in which case, set
c     the primal activity to the minimum allowed value and set the up
c     and dn bits.
c
                     pr_act(vr_n) = -usr_tl_mx_iz_pr_act
                     vr_st = vr_st + up_dn
                  endif
               endif
            else
               if (iand(vr_st, lb_bt) .eq. 0) then
c
c     A UB variable: Set the primal activity to the upper bound
c
                  pr_act_v = ubc(vr_n)*scl_v
                  if (pr_act_v .le. usr_tl_mx_iz_pr_act) then
c
c     Set the primal activity to the upper bound and set the dn bit...
c
                     pr_act(vr_n) = pr_act_v
                     vr_st = vr_st + dn_bt
                  else
c
c     ... unless this is an excessively high value, in which case, set
c     the primal activity to the maximum allowed value and set the up
c     and dn bits.
c
                     pr_act(vr_n) = usr_tl_mx_iz_pr_act
                     vr_st = vr_st + up_dn
                  endif
               else if (lbc(vr_n) .ne. ubc(vr_n)) then
c
c     A LB/UB variable: Set the variable to the lower bound
c
c                  vr_st = vr_st + up_bt
c                  pr_act(vr_n) = lbc(vr_n)*scl_v
c
c     Set the variable to the bound of smallest absolute value
c
                  if (abs(lbc(vr_n)) .le. abs(ubc(vr_n))) then
                     pr_act_v = lbc(vr_n)*scl_v
                     if (pr_act_v .ge. -usr_tl_mx_iz_pr_act) then
c
c     Set the primal activity to the lower bound and set the up bit...
c
                        pr_act(vr_n) = pr_act_v
                        vr_st = vr_st + up_bt
                     else
c
c     ... unless this is an excessively low value, in which case, set
c     the primal activity to the minimum allowed value and set the up
c     and dn bits.
c
                        pr_act(vr_n) = -usr_tl_mx_iz_pr_act
                        vr_st = vr_st + up_dn
                     endif
                  else
                     pr_act_v = ubc(vr_n)*scl_v
                     if (pr_act_v .le. usr_tl_mx_iz_pr_act) then
c
c     Set the primal activity to the upper bound and set the dn bit...
c
                        pr_act(vr_n) = pr_act_v
                        vr_st = vr_st + dn_bt
                     else
c
c     ... unless this is an excessively high value, in which case, set
c     the primal activity to the maximum allowed value and set the up
c     and dn bits.
c
                        pr_act(vr_n) = usr_tl_mx_iz_pr_act
                        vr_st = vr_st + up_dn
                     endif
                  endif
               else
c
c     A FX variable: Set the variable to the bound
c
                  pr_act_v = lbc(vr_n)*scl_v
                  if (pr_act_v .ge. -usr_tl_mx_iz_pr_act) then
c
c     Set the primal activity to the lower bound...
c
                     pr_act(vr_n) = pr_act_v
                  else
c
c     ... unless this is an excessively low value, in which case, set
c     the primal activity to the minimum allowed value and set the up
c     and dn bits.
c
                     pr_act(vr_n) = -usr_tl_mx_iz_pr_act
                     vr_st = vr_st + up_dn
                  endif
               endif
            endif
         else if (iand(vr_st, bp_bt) .ne. 0) then
c
c     Variable is BP: Set it at its break point [above/below] according
c     to the smaller absolute value of the [upper/lower] cost.
c
            vr_st = vr_st - iand(vr_st, lb_ub)
            pr_act_v = cbp(vr_n)*scl_v
            if (pr_act_v .ge. -usr_tl_mx_iz_pr_act .and.
     &           pr_act_v .le. usr_tl_mx_iz_pr_act) then
               pr_act(vr_n) = pr_act_v
               if (abs(lbc(vr_n)) .le. abs(ubc(vr_n))) then
c
c     Record the variable as being logically below its breakpoint
c
                  vr_st = vr_st + dn_bt + ub_bt
               else
c
c     Record the variable as being logically above its breakpoint
c
                  vr_st = vr_st + up_bt + lb_bt
               endif
            else if (pr_act_v .lt. -usr_tl_mx_iz_pr_act) then
c
c     ... unless this is an excessively low value... in which case, set
c     the primal activity to the minimum allowed value and indicate that
c     it is distinctly above its breakpoint
c
               pr_act(vr_n) = -usr_tl_mx_iz_pr_act
               vr_st = vr_st + up_dn + lb_bt
            else if (pr_act_v .lt. -usr_tl_mx_iz_pr_act) then
c
c     ... or high value, in which case, set the primal activity to the
c     maximum allowed value and indicate that it is distinctly below its
c     breakpoint.
c
               pr_act(vr_n) = usr_tl_mx_iz_pr_act
               vr_st = vr_st + up_dn + ub_bt
            endif
         else
c
c     Variable is PWL
c
         endif
         st(vr_n) = vr_st
 10   continue
      return
      end

C->>> -------------------------------------------------> ems_se_rsmi <<<
c     Sets the RSMI bounds and cost according to the type of variable
c     and the activity. Determines the basic variables, from the status,
c     maintaining any existing basis if possible.
c
      subroutine ems_se_rsmi(
     &     cbp, lbc, ubc, scl,
     &     st, pr_act, vr_in_r,
     &     rsmi_co, rsmi_lb, rsmi_ub)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      include 'RLCTVR.INC'
      include 'EMSMSG.INC'
      include 'EMSMSGN.INC'
      double precision cbp(0:mx_n_c+n_r)
      double precision lbc(0:mx_n_c+n_r)
      double precision ubc(0:mx_n_c+n_r)
      double precision scl(0:mx_n_c+n_r)
      integer st(0:mx_n_c+mx_n_r)
      double precision pr_act(0:mx_n_c+n_r)
      integer vr_in_r(0:n_r)
      double precision rsmi_co(0:mx_n_c+n_r)
      double precision rsmi_lb(0:mx_n_c+n_r)
      double precision rsmi_ub(0:mx_n_c+n_r)
      integer ix_n, vr_n, rsmi_vr_n
c      integer pwl_vr_n, pwl_vr_da_sa, pwl_vr_da_f_en
      integer r_n
      integer n_vr_in_r, n_vr_in_c
      double precision scl_v

      if (rq_inv .eq. rq_inv_no_rq_inv) then
c
c     If trying to avoid reinversion, check that the existing vr_in_r
c     consists of just basic variable.
c
         do 10, r_n = 1, n_r
            if (iand(st(vr_in_r(r_n)), bc_bt) .eq. 0) then
               rq_inv = rq_inv_nw_bs
               ml_da_st_msk =
     &              ml_da_st_msk - iand(ml_da_st_msk, ml_da_st_inv)
               go to 100
            end if
 10      continue
      end if
 100  continue
      n_vr_in_r = 0
      n_vr_in_c = 0
      do 110, ix_n = 1, n_r + n_c
         if (ix_n .le. n_c) then
            rsmi_vr_n = ix_n
         else
            rsmi_vr_n = ix_n + (mx_n_c-n_c)
         endif
         if (iand(st(rsmi_vr_n), su_vr_bs_bt) .ne. undn_vr_bs_st) then
            if (iand(st(rsmi_vr_n), bc_bt) .ne. 0) then
               if (n_vr_in_r .eq. n_r) then
c
c     Make the variable nonbasic if there is already a complete basis.
c
                  if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9900)
                  call ems_msg_wr_li(rsmi_msg_n)
                  rq_inv = rq_inv_nw_bs
                  ml_da_st_msk =
     &                 ml_da_st_msk - iand(ml_da_st_msk, ml_da_st_inv)
                  st(rsmi_vr_n) = st(rsmi_vr_n) - bc_bt
               endif
            else
               if (n_vr_in_c .eq. n_c) then
c
c     Make the variable basic if there is already a complete non-basis.
c
                  if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9910)
                  call ems_msg_wr_li(rsmi_msg_n)
                  rq_inv = rq_inv_nw_bs
                  ml_da_st_msk =
     &                 ml_da_st_msk - iand(ml_da_st_msk, ml_da_st_inv)
                  st(rsmi_vr_n) = st(rsmi_vr_n) + bc_bt
               endif
            endif
            if (iand(st(rsmi_vr_n), bc_bt) .ne. 0) then
c
c     The variable is basic.
c
               n_vr_in_r = n_vr_in_r + 1
               if (rq_inv .ne. rq_inv_no_rq_inv)
     &              vr_in_r(n_vr_in_r) = rsmi_vr_n
            else
c
c     The variable is nonbasic.
c
               n_vr_in_c = n_vr_in_c + 1
            end if
         endif
 110  continue
c
c     Second pass: Complete the sets of nonbasic and basic variables
c     (if necessary) using variables given the uncertain basic/nonbasic
c     status. Pass through all such variables now, starting with the
c     structurals, completing the set of nonbasic variables and then the
c     set of basic variables.
c
c     Initialise rsmi_lb/co/ub and corresponding status
c
      scl_v = one
      do 120, ix_n = 1, n_r + n_c
         if (ix_n .le. n_c) then
            vr_n =      ix_n
            rsmi_vr_n = ix_n
         else
            vr_n =      ix_n + (mx_n_c-n_c)
            rsmi_vr_n = ix_n + (mx_n_c-n_c)
         endif
c
c     Scale in forming the RSMI cost, lower bound and upper bound.
c     NB avoid scaling infinite bounds to be less than infinity.
c     If finite bounds are scaled to be infinite (unlikely) then this
c     shouldn't matter.
c
         if (iand(ml_da_st_msk, ml_da_st_scl_ml) .ne. 0)
     &        scl_v = scl(vr_n)
         if (iand(st(rsmi_vr_n), su_vr_bs_bt) .eq. undn_vr_bs_st) then
            st(rsmi_vr_n) = st(rsmi_vr_n) - undn_vr_bs_st
            if (n_vr_in_c .lt. n_c) then
c
c     Make the variable nonbasic
c
               n_vr_in_c = n_vr_in_c + 1
               st(rsmi_vr_n) = st(rsmi_vr_n) + non_bc_vr_bs_st
            else if (n_vr_in_r .lt. n_r) then
c
c     Make the variable basic
c
               n_vr_in_r = n_vr_in_r + 1
               vr_in_r(n_vr_in_r) = rsmi_vr_n
               st(rsmi_vr_n) = st(rsmi_vr_n) + bc_vr_bs_st
            else
               goto 8000
            endif
         endif
         if (iand(st(rsmi_vr_n), alt_bt) .eq. 0) then
c
c     Variable is Standard
c
            call ems_se_rsmi_std_vr(
     &           lbc(vr_n), cbp(vr_n), ubc(vr_n),
     &           rsmi_lb(rsmi_vr_n),
     &           rsmi_co(rsmi_vr_n),
     &           rsmi_ub(rsmi_vr_n),
     &           scl_v,
     &           st(rsmi_vr_n))
         else if (iand(st(rsmi_vr_n), bp_bt) .ne. 0) then
c
c     Variable is BP
c
            call ems_se_rsmi_bp_vr(
     &           lbc(vr_n), cbp(vr_n), ubc(vr_n),
     &           rsmi_lb(rsmi_vr_n),
     &           rsmi_co(rsmi_vr_n),
     &           rsmi_ub(rsmi_vr_n),
     &           pr_act(rsmi_vr_n),
     &           scl_v,
     &           st(rsmi_vr_n))
         else
c
c     Variable is PWL
c
         end if
 120  continue
c
c     Check that there is the right number of basic and nonbasic
c     variables.
c
      if (n_vr_in_c .ne. n_c .or. n_vr_in_r .ne. n_r) goto 8010
c
c     Record where the variable appears in vr_in_r.
c
      do 210, r_n = 1, n_r
         rsmi_vr_n = vr_in_r(r_n)
         st(rsmi_vr_n) = st(rsmi_vr_n) -
     &        iand(st(rsmi_vr_n), mx_mx_ml_a_dim) + r_n
 210  continue
      if (rq_inv .ne. rq_inv_no_rq_inv) then
c
c     Set av_eta_dse to be negative, in which case it is set to an
c     intelligent value after the first INVERT. NB Make sure it is set
c     to an intelligent value if a logical basis is identified and not
c     INVERTed.
c
         av_eta_dse = -one
      endif
c
c     Indicate that the model has vr_in_r correct.
c
      ml_da_st_msk = ior(ml_da_st_msk, ml_da_st_vr_in_r)
 7000 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
     &     n_vr_in_c, n_c,
     &     n_vr_in_r, n_r
      call ems_msg_wr_li(bug_msg_n)
      goto 7000
 8010 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9801)
     &     n_vr_in_c, n_c,
     &     n_vr_in_r, n_r
      call ems_msg_wr_li(bug_msg_n)
      goto 7000
 9900 format('Too many basic vars in status list')
 9910 format('Too many nonbasic vars in status list')
 9800 format('n_vr_in_c .ge. n_c .and. n_vr_in_r .ge. n_r',
     &     ' but there are still variables of uncertain status ',
     &     i9, i9, i9, i9)
 9801 format('n_vr_in_c .ne. n_c .or. n_vr_in_r .ne. n_r',
     &     i9, i9, i9, i9)
      end

C->>> ----------------------------------------------> ems_iz_vr_in_r <<<
c     Initialises vr_in_r from the status.
c     A variable is
c     Basic    if bc_bt and at least one of up_bt and dn_bt are set
c     Nonbasic if bc_bt is not set
c     Either   if bc_bt and  neither     of up_bt or  dn_bt is  set
c
c     After all the basic variables have been identified, an incomplete
c     basis is extended using variables which are neither Basic nor
c     Nonbasic
c
      subroutine ems_iz_vr_in_r(st, vr_in_r)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      include 'RLCTVR.INC'
      include 'EMSMSG.INC'
      include 'EMSMSGN.INC'
      integer st(0:mx_n_c+mx_n_r)
      integer vr_in_r(0:n_r)
      integer ix_n, vr_n, vr_st
      integer n_vr_in_r, n_vr_in_c

      n_vr_in_r = 0
      n_vr_in_c = 0
      do 110, ix_n = 1, n_r + n_c
         if (ix_n .le. n_c) then
            vr_n = ix_n
         else
            vr_n = ix_n + (mx_n_c-n_c)
         endif
         vr_st = st(vr_n)
         if (iand(vr_st, su_vr_bs_bt) .eq. undn_vr_bs_st) goto 110
         if (iand(vr_st, bc_bt) .ne. 0) then
            if (n_vr_in_r .eq. n_r) then
c
c     Make the variable nonbasic if there is already a complete basis.
c
               if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9900)
               call ems_msg_wr_li(rsmi_msg_n)
               rq_inv = rq_inv_nw_bs
               ml_da_st_msk =
     &              ml_da_st_msk - iand(ml_da_st_msk, ml_da_st_inv)
               vr_st = vr_st - bc_bt
            endif
         else
            if (n_vr_in_c .eq. n_c) then
c
c     Make the variable basic if there is already a complete non-basis.
c
               if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9910)
               call ems_msg_wr_li(rsmi_msg_n)
               rq_inv = rq_inv_nw_bs
               ml_da_st_msk =
     &              ml_da_st_msk - iand(ml_da_st_msk, ml_da_st_inv)
               vr_st = ior(vr_st, bc_vr_bs_st)
            endif
         endif
         if (iand(vr_st, bc_bt) .ne. 0) then
c
c     The variable is basic.
c
            n_vr_in_r = n_vr_in_r + 1
            vr_in_r(n_vr_in_r) = vr_n
            vr_st = vr_st - iand(vr_st, mx_mx_ml_a_dim) + n_vr_in_r
         else
c
c     The variable is nonbasic.
c
            n_vr_in_c = n_vr_in_c + 1
         end if
         st(vr_n) = vr_st
 110  continue
      if (n_vr_in_r+n_vr_in_c .ge. n_r+n_c) goto 1000
c
c     Second pass: Complete the sets of nonbasic and basic variables
c     (if necessary) using variables given the uncertain basic/nonbasic
c     status. Pass through all such variables now, starting with the
c     structurals, completing the set of nonbasic variables and then the
c     set of basic variables.
c
      do 210, ix_n = 1, n_r + n_c
         if (ix_n .le. n_c) then
            vr_n = ix_n
         else
            vr_n = ix_n + (mx_n_c-n_c)
         endif
         vr_st = st(vr_n)
         if (iand(vr_st, su_vr_bs_bt) .ne. undn_vr_bs_st) goto 210
         vr_st = vr_st - iand(vr_st, su_vr_bs_bt)
         if (n_vr_in_c .lt. n_c) then
c
c     Make the variable nonbasic
c
            n_vr_in_c = n_vr_in_c + 1
            vr_st = vr_st + non_bc_vr_bs_st
         else if (n_vr_in_r .lt. n_r) then
c
c     Make the variable basic
c
            n_vr_in_r = n_vr_in_r + 1
            vr_in_r(n_vr_in_r) = vr_n
            vr_st = vr_st + bc_vr_bs_st
         else
            goto 8000
         endif
 210  continue
 1000 continue
c
c     Check that there is the right number of basic and nonbasic
c     variables.
c
      if (n_vr_in_c .ne. n_c .or. n_vr_in_r .ne. n_r) goto 8010
c
c     Indicate that vr_in_r is correct for the model.
c
      ml_da_st_msk = ior(ml_da_st_msk, ml_da_st_vr_in_r)
c
c     Indicate that the following are not correct for the model:
c
c     vr_in_c, INVERT, status/primal activities, basic primal
c     activities, edge weights and row-wise representation of matrix
c     columns being priced.
c
      ml_da_st_msk = ml_da_st_msk
     &     - iand(ml_da_st_msk, ml_da_st_vr_in_c)
     &     - iand(ml_da_st_msk, ml_da_st_inv)
     &     - iand(ml_da_st_msk, ml_da_st_vr_st_fm_act)
     &     - iand(ml_da_st_msk, ml_da_st_bc_pr_act)
     &     - iand(ml_da_st_msk, ml_da_st_ed_wt)
     &     - iand(ml_da_st_msk, ml_da_st_r_mtx)
 7000 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
     &     n_vr_in_c, n_c,
     &     n_vr_in_r, n_r
      call ems_msg_wr_li(bug_msg_n)
      goto 7000
 8010 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9801)
     &     n_vr_in_c, n_c,
     &     n_vr_in_r, n_r
      call ems_msg_wr_li(bug_msg_n)
      goto 7000
 9900 format('Too many basic vars in status list')
 9910 format('Too many nonbasic vars in status list')
 9800 format('n_vr_in_c .ge. n_c .and. n_vr_in_r .ge. n_r',
     &     ' but there are still variables of uncertain status ',
     &     i9, i9, i9, i9)
 9801 format('n_vr_in_c .ne. n_c .or. n_vr_in_r .ne. n_r',
     &     i9, i9, i9, i9)
      end

C->>> ----------------------------------------------> ems_iz_rsmi_bd <<<
c     Initialise the RSMI bounds and cost according to the type of
c     variable and the activity.
c
      subroutine ems_iz_rsmi_bd(
     &     cbp, lbc, ubc, scl,
     &     st, pr_act, vr_in_r,
     &     rsmi_co, rsmi_lb, rsmi_ub)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      include 'RLCTVR.INC'
      include 'EMSMSG.INC'
      include 'EMSMSGN.INC'
      double precision cbp(0:mx_n_c+n_r)
      double precision lbc(0:mx_n_c+n_r)
      double precision ubc(0:mx_n_c+n_r)
      double precision scl(0:mx_n_c+n_r)
      integer st(0:mx_n_c+mx_n_r)
      double precision pr_act(0:mx_n_c+n_r)
      integer vr_in_r(0:n_r)
      double precision rsmi_co(0:mx_n_c+n_r)
      double precision rsmi_lb(0:mx_n_c+n_r)
      double precision rsmi_ub(0:mx_n_c+n_r)
      integer ix_n, vr_n, vr_st
      double precision scl_v, rcp_scl_v
      double precision pwl_lc, pwl_bp, pwl_uc
c
c     Initialise rsmi_lb/co/ub and corresponding status
c
      scl_v = one
      do 110, ix_n = 1, n_r + n_c
         if (ix_n .le. n_c) then
            vr_n = ix_n
         else
            vr_n = ix_n + (mx_n_c-n_c)
         endif
c
c     Scale in forming the RSMI cost, lower bound and upper bound.
c     NB avoid scaling infinite bounds to be less than infinity.
c     If finite bounds are scaled to be infinite (unlikely) then this
c     shouldn't matter.
c
         if (iand(ml_da_st_msk, ml_da_st_scl_ml) .ne. 0)
     &        scl_v = scl(vr_n)
         vr_st = st(vr_n)
         if (iand(vr_st, alt_bt) .eq. 0) then
c
c     Variable is Standard
c
            if (lbc(vr_n) .gt. -inf) then
               rsmi_lb(vr_n) = min(max(lbc(vr_n)*scl_v, -inf), inf)
            else
               rsmi_lb(vr_n) = -inf
            endif
            rsmi_co(vr_n) =    min(max(cbp(vr_n)/scl_v, -inf), inf)
            if (ubc(vr_n) .lt.  inf) then
               rsmi_ub(vr_n) = min(max(ubc(vr_n)*scl_v, -inf), inf)
            else
               rsmi_ub(vr_n) = inf
            endif
         else if (iand(vr_st, bp_bt) .ne. 0) then
c
c     Variable is BP
c
            rcp_scl_v = one/scl_v
            pwl_lc = min(max(lbc(vr_n)*rcp_scl_v, -inf), inf)
            pwl_bp = min(max(cbp(vr_n)*scl_v,     -inf), inf)
            pwl_uc = min(max(ubc(vr_n)*rcp_scl_v, -inf), inf)
            if (iand(vr_st, ub_bt) .ne. 0) then
               rsmi_lb(vr_n) = pwl_uc - pwl_lc
               rsmi_co(vr_n) = pwl_lc
               rsmi_ub(vr_n) = pwl_bp
            else
               rsmi_lb(vr_n) = pwl_bp
               rsmi_co(vr_n) = pwl_uc
               rsmi_ub(vr_n) = pwl_lc - pwl_uc
            endif
         else
c
c     Variable is PWL
c
         end if
         st(vr_n) = vr_st
 110  continue
      return
      end

C->>> ------------------------------------------> ems_se_rsmi_std_vr <<<
c     Set the values of rsmi_vr_lb, rsmi_vr_co, rsmi_vr_ub and the
c     LB/UB bits in the status.
c
      subroutine ems_se_rsmi_std_vr(
     &     ml_vr_lb, ml_vr_co, ml_vr_ub,
     &     rsmi_vr_lb, rsmi_vr_co, rsmi_vr_ub,
     &     scl_v,
     &     vr_st)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      integer vr_st
      double precision ml_vr_lb, ml_vr_co, ml_vr_ub
      double precision rsmi_vr_lb, rsmi_vr_co, rsmi_vr_ub
      double precision scl_v

      if (ml_vr_lb .gt. -inf) then
         rsmi_vr_lb = ml_vr_lb*scl_v
         vr_st = ior(vr_st, lb_bt)
      else
         rsmi_vr_lb = -inf
         vr_st = vr_st - iand(vr_st, lb_bt)
      endif
      rsmi_vr_co = ml_vr_co/scl_v
      if (ml_vr_ub .lt.  inf) then
         rsmi_vr_ub = ml_vr_ub*scl_v
         vr_st = ior(vr_st, ub_bt)
      else
         rsmi_vr_ub = inf
         vr_st = vr_st - iand(vr_st, ub_bt)
      endif
      return
      end

C->>> -------------------------------------------> ems_se_rsmi_bp_vr <<<
c     Set the values of rsmi_vr_lb, rsmi_vr_co, rsmi_vr_ub and the
c     LB/UB bits in the status.
c
      subroutine ems_se_rsmi_bp_vr(
     &     ml_vr_lb, ml_vr_co, ml_vr_ub,
     &     rsmi_vr_lb, rsmi_vr_co, rsmi_vr_ub,
     &     rsmi_vr_pr_act,
     &     scl_v,
     &     vr_st)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'RLCTVR.INC'
      integer vr_st
      double precision ml_vr_lb, ml_vr_co, ml_vr_ub
      double precision rsmi_vr_lb, rsmi_vr_co, rsmi_vr_ub
      double precision rsmi_vr_pr_act
      double precision scl_v
      double precision pwl_lc, pwl_bp, pwl_uc
      double precision rsdu, rcp_scl_v

      rcp_scl_v = one/scl_v
      pwl_lc = ml_vr_lb*rcp_scl_v
      pwl_bp = ml_vr_co*scl_v
      pwl_uc = ml_vr_ub*rcp_scl_v
      rsdu = rsmi_vr_pr_act - pwl_bp
      if (rsdu .lt. -tl_pr_ifs .or.
     &     rsdu .le. tl_pr_ifs .and. iand(vr_st, ub_bt) .ne. 0) then
c
c     Variable is below its break point or degenerate and logically
c     below its breakpoint:
c     Set lb:co:ub to uc-lc:lc:bp
c
         rsmi_vr_lb = pwl_uc - pwl_lc
         rsmi_vr_co = pwl_lc
         rsmi_vr_ub = pwl_bp
         vr_st = ior(vr_st, ub_bt) - iand(vr_st, lb_bt)
      else
c
c     Variable is above its break point or degenerate and logically
c     above its breakpoint:
c     Set lb:co:ub to bp:uc:lc-uc
c
         rsmi_vr_lb = pwl_bp
         rsmi_vr_co = pwl_uc
         rsmi_vr_ub = pwl_lc - pwl_uc
         vr_st = ior(vr_st, lb_bt) - iand(vr_st, ub_bt)
      endif
      return
      end

max_stars_repo_name: ERGO-Code/EMSOL

max_issues_repo_name: ERGO-Code/EMSOL

max_forks_repo_name: ERGO-Code/EMSOL


text:
C
C $Id: mdpset.f,v 1.10 2008-09-18 00:42:17 kennison Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE MDPSET (ARG1,ARG2,ARG3,ARG4,ARG5)
C
        CHARACTER*(*)    ARG1
        DOUBLE PRECISION ARG2(2),ARG3(2),ARG4(2),ARG5(2)
C
C Declare required common blocks.  See MAPBDX for descriptions of these
C common blocks and the variables in them.
C
        COMMON /MAPCM4/  GRDR,GRID,GRLA,GRLO,GRPO,OTOL,PDRE,PLA1,PLA2,
     +                   PLA3,PLA4,PLB1,PLB2,PLB3,PLB4,PLNO,PLTO,ROTA,
     +                   SRCH,XLOW,XROW,YBOW,YTOW,IDOT,IDSH,IDTL,ILCW,
     +                   ILTS,JPRJ,ELPF,INTF,LBLF,PRMF
        DOUBLE PRECISION GRDR,GRID,GRLA,GRLO,GRPO,OTOL,PDRE,PLA1,PLA2,
     +                   PLA3,PLA4,PLB1,PLB2,PLB3,PLB4,PLNO,PLTO,ROTA,
     +                   SRCH,XLOW,XROW,YBOW,YTOW
        INTEGER          IDOT,IDSH,IDTL,ILCW,ILTS,JPRJ
        LOGICAL          ELPF,INTF,LBLF,PRMF
        SAVE   /MAPCM4/
C
        COMMON /MAPCM5/  DDCT(5),DDCL(5),LDCT(6),LDCL(6),PDCT(19),
     +                   PDCL(19)
        CHARACTER*2      DDCT,DDCL,LDCT,LDCL,PDCT,PDCL
        SAVE   /MAPCM5/
C
C Declare local variables.
C
        INTEGER          I
C
C Check for an uncleared prior error.
C
        IF (ICFELL('MDPSET - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
C
C Transfer the parameters defining the map limits.
C
        I=IDICTL(ARG1,LDCT,6)
        IF (I.EQ.0) I=IDICTL(ARG1,LDCL,6)
        IF (I.EQ.0) GO TO 901
        ILTS=I
C
        PLA1=ARG2(1)
        PLA2=ARG3(1)
        PLA3=ARG4(1)
        PLA4=ARG5(1)
C
        IF (I.EQ.3) THEN
          PLB1=ARG2(2)
          PLB2=ARG3(2)
          PLB3=ARG4(2)
          PLB4=ARG5(2)
        END IF
C
C Set the flag to indicate that initialization is now required.
C
        INTF=.TRUE.
C
C Done.
C
        RETURN
C
C Error exit.
C
  901   CALL MDPCEM ('MDPSET - UNKNOWN MAP AREA SPECIFIER ',ARG1,2,1)
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
Graduate Students are students who are pursuing a Masters or Doctorate degree (or both).  There are a few thousand of them in Davis.  Since professional schools graduate school isnt really about taking a set number of units and requirements vary by department, they have been known to stick around for a long time.  That is, until they get sick of school or funding runs out.  Only 45% of UCD graduate students actually complete their program. To figure out when a grad student will graduate, retreat to a safe distance and ask him or her how long until theyre done.  Then multiply this value by 1.5 to get a good estimate.  Doctoral students may spout an acronym (A.B.D or E.B.D., for all but dissertation or everything but dissertation, respectively) when this question is asked; it means that they have completed all their doctoral coursework and passed their quals (qualifying exams), and only need to complete their dissertation in order to get the degree.  It is generally considered unwise to inquire into a grad students research: either its not going well, in which case the grad student will probably try to hurt you for reminding him or her of it; or its going well, in which case youll be treated to a rambling and impenetrable lecture on http://www.phdcomics.com/comics/archive.php?comicid976 Furrier Transformers or something.

Grad Students are members of Graduate Student Association GSA and are served by the Division of Graduate Studies. Many of them are also Teaching Assistants, and they can participate in community service through the Graduate Student Community Service Committee. They are often Q & A for Prospective Academics planning on becoming academics.

You can recognize them by the fact that they are Users/IDoNotExist too busy with work to describe themselves in the Davis Wiki.  Theyre generally more disheveled than undergraduate students, as theyre more concerned with academics than appearance, with the males often sporting the thesis beard. Another way to recognize graduate students is to go to campus on weekends and holidays.  Notice that youll see less people, but youll still see some.  To your surprise, some of those people are still acting like its still a school daybackpacks on, hurrying about like theyre about to miss a class.  If youve ever overheard someone making jokes about Foucault, bounded rationality, PCR machines, hegemony, or anything else that seem like something painfully awkward to be talking about outside of discussionyou can be 95% certain you are listening to a graduate student.

The natural habitat of a graduate student includes Mishkas, Sophias, and Delta of Venus. They are known to frequent trivia nights since free food is one of their three food groups (other twobeer and coffee). The courtship grounds of graduate students are largely unknown.  Grad students are nocturnal.  If encountered during the daytime, they are frequently irritable, and may attack without provocation.  Approach with care, and do not taunt.

http://www.phdcomics.com/ Piled Higher and Deeper, a popular Webcomics webcomic, can often be seen attached to office doors or bulletin boards throughout campus buildings, departments, and research labs.  A satire of life in graduate school, grad students everywhere can identify with it.  It pretty much accurately describes their lives and experiences here at Davis as well.

 And rewrite it, and include all of the picayune and frequently contradictory suggestions made by their committee, and defend it in front of a bunch of professors, and print it out so that it meets every last one of the University of Californias formatting requirements, and then rewrite it to account for new research thats come out since the first draft.

The other 5% are math majors. Sometimes math majors are even dorkier.

see also Graduate Student Social Event Calendar

Practical Advice

If you plan to go forward in academia, make sure that you network and volunteer.  As you move through junior faculty positions and write your first grant proposals, they will ask you to list synergistic efforts and activities not directly related to your graduate work.  Getting into grad school is not the end of your resume building. Once you obtain your advanced degree, to continue in academia you will have to defend your commitment to expanding knowledge (research), diversity, and the educational process.  Despite the work load, you should choose committees and services that you can volunteer on and try to forge relationships with others in your field.

Then there are those grad students too busy describing themselves on Davis Wiki to do their work. Thats why you dont get your papers back on time. Users/TonyMagagna

Im one of those grad students, too.  Users/BrentLaabs

I wish I was a dissheveled grad student ( too.  Users/MichaelGiardina Its were, Michael. Remember: the subjunctive follows the conditional. Dr. NO (a grad student who made it outBA to PhDin just under 7 and with only two breakdowns that required hospitalization)
  No you dont. Unless you are a masochist that is.  Users/EricKlein

Oh look, theres one of those dissheveled grad students now.  Were working weekends, holidays, AND nights. Users/DanComing

This must be a description of Grad Students in real fields like Science or Math.  English Grad students dont really do anything, unless theres something Im supposed to be doing that I dont know about.  Basically, were the people in the movie theater laughing at the movie because the plot is cliche or we see some reference to Salinger or Derrida that no one else gets.  Our job is to be pretentious and look good.  Thats it.  If theres work, no one has told me about it.  Users/NathanMilos
  Seriously.  Why would we be updating wiki pages about ourselves if we had anything better/more urgent/more exciting to do? Oh, the humanities.  Users/MelanieMadden

Does anyone know the location of the cute wooden apt. building that is shown on the Apartments page, next to the listings for Central Davis?  I want to live there!  There is a caption about how grad students may appreciate the smaller complexes......believe me, Im appreciatin and I want to appreciate my UHaul right on over there. Users/SignoraStella
  Its on K Street, between Fifth Street and Sixth Street Users/AlphaDog

Where is Users/SirPallas when you need him to talk about his theory of hazing? Users/KenBloom

20060911 11:04:44 nbsp Its generally considered rude to ask a PhD student how many years we have left. Its a little more acceptable to ask the question of a masters student.  Users/AmyGoogenspa
    I have not met anyone in my program who would consider that rude.  When are you getting out is pretty casual conversation, as I would imagine it is in prison.  If youre offended at someone asking that, youre taking yourself a little too seriously. Users/PhilSpear
        I would say the bigger problem is that most of them dont have an actual answer but feel obliged to make up some story, possibly unrelated to the question. Besides, whats the rush? Users/AlexMandel
           See http://www.phdcomics.com/comics/archive.php?comicid47 PhDComics.
    Ive never had any sort of negative reaction, and Ive been asking people from various departments at various universities around the country for many years.  Its basically part of the whats your research? and whats your group like? (or, if you know the department well, whos your advisor?) series of questions that usually pop up over lunch or at a diner.  Most people ramble on for a bit about it, talking about their planned series of papers or any hold ups they are dealing with. Users/JabberWokky jw

20061218 15:15:51 nbsp HEHEHE Users/HilarieLloyd

20080211 14:41:52 nbsp The hallmark of grad students  going home at 11pm so they can get up early on Saturday to work on their experiment! Users/pizzapieGuy

20080211 15:17:35 nbsp Grad students are often found munching at food sources that are open late.  We dont have time to actually eat something until 11 or 12 at night, by which time most food sources are no longer offering food. Users/IDoNotExist

Im fairly new to In n Out Burger since I didnt have it back at home in Utah. The fact that theyre open until 1 am is how I satisfy my hunger on those late nights! Users/RichardTrendall

20110406 12:13:51 nbsp Whoever wrote this page deserves a cookie! If you are still in Davis contact me and I will bake you some, if not thank you for the great laugh.  Users/DanielleC


program example_convection
    implicit none

    ! Number of cells
    integer, parameter :: nx=50, ny=50

    ! Variables
    real(kind=8),dimension(nx,ny) :: AE, AW, AN, AS, Su, Sp, AP ! Upwinding coefficients
    real(kind=8),dimension(nx,ny) :: A_E, A_W, A_N, A_S, S_u, S_p, A_P ! Second-order coefficients
    real(kind=8),dimension(0:nx+1,0:ny+1) :: phi
    real(kind=8) :: x_min, y_min, x_max, y_max, u, v, dx, dy, B, rho, m_e, m_w, m_n, m_s, phi_wbc, phi_sbc, EQ3, EQ2, phi_new
    real(kind=8) :: max_diff, new_diff, x, y
    integer :: i, j, iterations, max_iter
    character(100) :: filename

    ! Welcome message
    write(*,*)
    write(*,*) "Example Convection Problem"
    write(*,*) "--------------------------"

    ! Solver options
    B = 1.0
    max_iter = 10000000

    ! Domain
    x_min = 0.0
    x_max = 1.0
    y_min = 0.0
    y_max = 1.0

    ! Flow options
    u = 1.0
    v = 1.0
    rho = 1.0

    ! Boundary conditions
    phi_wbc = 1.0
    phi_sbc = 0.0

    ! Grid
    dx = (x_max-x_min)/nx
    dy = (y_max-y_min)/ny

    ! Initialize phi
    phi = 0.0
    phi(0,0:nx+1) = phi_sbc ! Southern boundary
    phi(0:ny+1,0) = phi_wbc ! Western boundary
    phi(0,0) = 0.5*(phi_sbc+phi_wbc) ! Southwest corner

    ! Mass flows
    m_e = rho*u*dy
    m_w = rho*u*dy
    m_n = rho*v*dx
    m_s = rho*v*dx

    ! Initialize all coefficients to interior values
    ! Upwind
    AE = 0.0
    AW = m_w
    AN = 0.0
    AS = m_s
    Su = 0.0
    Sp = 0.0

    ! Second-order
    A_E = -0.5*m_e
    A_W = 0.5*m_w
    A_N = -0.5*m_n
    A_S = 0.5*m_s
    S_u = 0.0
    S_p = 0.0

    ! West boundary; Dirichlet condition
    ! Upwind
    AE(2:ny-1,1) = 0.0
    AW(2:ny-1,1) = 0.0
    Sp(2:ny-1,1) = -m_e
    Su(2:ny-1,1) = m_w*phi_wbc

    ! Second-order
    A_E(2:ny-1,1) = -0.5*m_e
    A_W(2:ny-1,1) = 0.0
    S_p(2:ny-1,1) = -m_e
    S_u(2:ny-1,1) = m_w*phi_wbc

    ! East boundary; Neumann condition
    ! Upwind
    AE(2:ny-1,nx) = 0.0
    Aw(2:ny-1,nx) = m_w
    Sp(2:ny-1,nx) = -(m_e-m_w)

    ! Second-order
    A_E(2:ny-1,nx) = 0.0
    A_W(2:ny-1,nx) = 0.5*m_w
    S_p(2:ny-1,nx) = -(m_e-m_w)

    ! South boundary; Dirichlet condition
    ! Upwind
    AN(1,2:nx-1) = 0.0
    AS(1,2:nx-1) = 0.0
    Sp(1,2:nx-1) = -m_n
    Su(1,2:nx-1) = m_s*phi_sbc

    ! Second-order
    A_N(1,2:nx-1) = -0.5*m_n
    A_S(1,2:nx-1) = 0.0
    S_p(1,2:nx-1) = -m_n
    S_u(1,2:nx-1) = m_s*phi_sbc

    ! North boundary; Neumann condition
    ! Upwind
    AN(ny,2:nx-1) = 0.0
    AS(ny,2:nx-1) = m_s
    Sp(ny,2:nx-1) = -(m_n-m_s)

    ! Second-order
    A_N(ny,2:nx-1) = 0.0
    A_S(ny,2:nx-1) = 0.5*m_s
    S_p(ny,2:nx-1) = -(m_n-m_s)

    ! Northwest corner
    Su(ny,1) = m_w*phi_wbc
    Sp(ny,1) = -m_e-(m_n-m_s)
    AS(ny,1) = m_s
    AE(ny,1) = 0.0
    AW(ny,1) = 0.0
    AN(ny,1) = 0.0

    S_u(ny,1) = m_w*phi_wbc
    S_p(ny,1) = -m_e-(m_n-m_s)
    A_W(ny,1) = 0.0
    A_N(ny,1) = 0.0
    A_S(ny,1) = 0.5*m_s
    A_E(ny,1) = -0.5*m_e

    ! Southwest corner
    Su(1,1) = m_w*phi_wbc+m_s*phi_sbc
    Sp(1,1) = -m_e-m_n
    AE(1,1) = 0.0
    AW(1,1) = 0.0
    AN(1,1) = 0.0
    AS(1,1) = 0.0

    S_u(1,1) = m_w*phi_wbc+m_s*phi_sbc
    S_p(1,1) = -m_e-m_n
    A_W(1,1) = 0.0
    A_E(1,1) = -0.5*m_e
    A_S(1,1) = 0.0
    A_N(1,1) = -0.5*m_n

    ! Southeast corner
    Su(1,nx) = m_s*phi_sbc
    Sp(1,nx) = -m_n-(m_e-m_w)
    AW(1,nx) = m_w
    AN(1,nx) = 0.0
    AE(1,nx) = 0.0
    AS(1,nx) = 0.0

    S_u(1,nx) = m_s*phi_sbc
    S_p(1,nx) = -m_n-(m_e-m_w)
    A_E(1,nx) = 0.0
    A_W(1,nx) = 0.5*m_w
    A_S(1,nx) = 0.0
    A_N(1,nx) = -0.5*m_n

    ! Northeast corner
    Sp(ny,nx) = -(m_n-m_s)-(m_e-m_w)
    AN(ny,nx) = 0.0
    AE(ny,nx) = 0.0
    Aw(ny,nx) = m_w
    AS(ny,nx) = m_s

    S_p(ny,nx) = -(m_n-m_s)-(m_e-m_w)
    A_N(ny,nx) = 0.0
    A_S(ny,nx) = 0.5*m_s
    A_E(ny,nx) = 0.0
    A_W(ny,nx) = 0.5*m_w

    ! Calculate AP and A_P
    AP = AN+AE+AS+AW-SP+m_e-m_w+m_n-m_s
    A_P = A_N+A_E+A_S+A_W-S_P+m_e-m_w+m_n-m_s

    ! Iterate
    iterations = 0
    max_diff = 1.0
    write(*,*)
    write(*,*) "   Iteration      Max Correction"
    write(*,*) "--------------------------------"
    do while (.not. max_diff .eq. 0.0 .and. iterations<max_iter)

        ! Reset for iteration
        iterations = iterations+1
        max_diff = 0.0

        ! Loop through control points
        do i = 1,ny
            do j = 1,nx

                ! Upwind
                EQ3 = AE(i,j)*phi(i,j+1)+AW(i,j)*phi(i,j-1)+AN(i,j)*phi(i+1,j)+AS(i,j)*phi(i-1,j)+Su(i,j)

                ! Second order
                EQ2 = A_P(i,j)*phi(i,j)-A_E(i,j)*phi(i,j+1)-A_W(i,j)*phi(i,j-1)-A_N(i,j)*phi(i+1,j)-A_S(i,j)*phi(i-1,j)-S_u(i,j)

                ! Blend
                phi_new = EQ3/AP(i,j)-B/AP(i,j)*(EQ2-(AP(i,j)*phi(i,j)-EQ3))

                ! Check convergence
                new_diff = abs(phi_new-phi(i,j))
                if (new_diff > max_diff) then
                    max_diff = new_diff
                end if

                ! Assign
                phi(i,j) = phi_new

            end do
        end do
        write(*,*) iterations, max_diff
    end do

    ! Finalize north and east boundaries
    phi(ny+1,:) = phi(ny,:)
    phi(:,nx+1) = phi(:,nx)
    phi(ny+1,nx+1) = 0.5*(phi(ny+1,nx)+phi(ny,nx+1))

    ! Write results to file
    ! Open file
    filename = 'results.csv'
    open(1,file=filename)

        ! Header
        write(1,*) 'x,y,z,phi'

        ! Control points
        do i = 1,ny
            do j =1,nx
                x = x_min+0.5*dx*(2*j-1)
                y = y_min+0.5*dy*(2*i-1)
                write(1,*) x, ',', y, ',', 0.0, ',', phi(i,j)
            end do
        end do

        ! West and east boundaries
        do i = 1,ny
            y = y_min+0.5*dy*(2*i-1)
            write(1,*) x_min, ',', y, ',', 0.0, ',', phi(i,0) ! West
            write(1,*) x_max, ',', y, ',', 0.0, ',', phi(i,nx+1) ! East
        end do

        ! North and south boundaries
        do i = 1,nx
            x = x_min+0.5*dx*(2*i-1)
            write(1,*) x, ',', y_min, ',', 0.0, ',', phi(0,i) ! South
            write(1,*) x, ',', y_max, ',', 0.0, ',', phi(ny+1,i) ! North
        end do

        ! Corners
        write(1,*) x_min, ',', y_min, ',', 0.0, ',', phi(0,0) ! SW
        write(1,*) x_min, ',', y_max, ',', 0.0, ',', phi(ny+1,0) ! SE
        write(1,*) x_max, ',', y_min, ',', 0.0, ',', phi(0,nx+1) ! NW
        write(1,*) x_max, ',', y_max, ',', 0.0, ',', phi(ny+1,nx+1) ! NE

    close(1)

end program example_convection