  MODULE sfcphys_module

  implicit none

  CONTAINS

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine getcecd(u1,v1,s1,za,u10,v10,s10,xland,znt,ust,cd,ch,cq,ubar,vbar,sbar,use_avg_sfc)
      use input
      use constants
      use mpi
      implicit none

      real, intent(inout), dimension(ib:ie,jb:je) :: u1,v1,s1
      real, intent(in), dimension(ib:ie,jb:je) :: za
      real, intent(inout), dimension(ibl:iel,jbl:jel) :: u10,v10,s10
      real, intent(in), dimension(ib:ie,jb:je) :: xland
      real, intent(inout), dimension(ib:ie,jb:je) :: znt,ust,cd,ch,cq
      double precision, intent(inout) :: ubar,vbar,sbar
      logical, intent(in) :: use_avg_sfc

      integer :: i,j,i1,i2,j1,j2,n,nmax
      real :: wsp,wlast,var,rznt,usave,vsave,ssave
      double precision :: temd

      real, parameter :: dcd1  =  1.0e-3
      real, parameter :: dcd2  =  2.4e-3
      real, parameter :: dwsp1 =  5.0
      real, parameter :: dwsp2 = 25.0

      real, parameter :: dfac = (dcd2-dcd1)/(dwsp2-dwsp1)

!-----------------------------------------------------------------------
!
!  This subroutine determines several important variables at the surface
!  (eg, drag coefficient, roughness length, friction velocity).
!
!  Note:  This is a simple scheme that assumes neutral stability
!
!         For stability effects, use sfcmodel = 2 or 3
!
!  180213:  Rearrange into different sections (no need to iterate in most cases)
!
!-----------------------------------------------------------------------

!!!    print *,'  za = ',za(1,1),za(ni/2,1),za(ni,1)

    nmax = 0

    i1 = 1
    i2 = ni

    j1 = 1
    j2 = nj

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

    if( use_avg_sfc )then

      i1 = 1
      i2 = 1
      j1 = 1
      j2 = 1

      ubar = 0.0
      vbar = 0.0
      sbar = 0.0

      do j=1,nj
      do i=1,ni
        ubar = ubar + u1(i,j)
        vbar = vbar + v1(i,j)
        sbar = sbar + s1(i,j)
      enddo
      enddo

      call MPI_ALLREDUCE(MPI_IN_PLACE,ubar,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
      call MPI_ALLREDUCE(MPI_IN_PLACE,vbar,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
      call MPI_ALLREDUCE(MPI_IN_PLACE,sbar,1,MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)

      temd = 1.0/dble(nx*ny)
      ubar = ubar*temd
      vbar = vbar*temd
      sbar = sbar*temd

      usave = u1(1,1)
      vsave = v1(1,1)
      ssave = s1(1,1)

      u1(1,1) = ubar
      v1(1,1) = vbar
      s1(1,1) = sbar

    endif

!-----------------------------------------------------------------------
!  cd,z0,ust:

  sfc_options:  &
  IF( set_znt.eq.1 )THEN
    ! specify z0 (roughness length):

    !$omp parallel do default(shared)   &
    !$omp private(i,j,wlast,n,var,rznt)
    do j=j1,j2
    do i=i1,i2
      znt(i,j) = max(1.0e-12,cnst_znt)
      rznt = 1.0/znt(i,j)
      cd(i,j) = ( karman/alog((10.0+znt(i,j))*rznt) )**2
      var = alog((10.0+znt(i,j))*rznt)/alog((za(i,j)+znt(i,j))*rznt)
      u10(i,j) = u1(i,j)*var
      v10(i,j) = v1(i,j)*var
      s10(i,j) = s1(i,j)*var
      ust(i,j) = max( s1(i,j)*karman/alog((za(i,j)+znt(i,j))*rznt) , 1.0e-6 )
    enddo
    enddo

  ELSEIF( set_ust.eq.1 )THEN  sfc_options
      ! specify ustar (friction velocity):

    !$omp parallel do default(shared)   &
    !$omp private(i,j,wlast,n,var,rznt)
    do j=j1,j2
    do i=i1,i2
      ust(i,j) = cnst_ust
      znt(i,j) = za(i,j)/(exp(max(1.0e-10,s1(i,j))*karman/max(1.0e-10,ust(i,j)))-1.0)
      rznt = 1.0/znt(i,j)
      cd(i,j) = ( karman/alog((10.0+znt(i,j))*rznt) )**2
      var = alog((10.0+znt(i,j))*rznt)/alog((za(i,j)+znt(i,j))*rznt)
      u10(i,j) = u1(i,j)*var
      v10(i,j) = v1(i,j)*var
      s10(i,j) = s1(i,j)*var
    enddo
    enddo

  ELSEIF( cecd.eq.1 )THEN  sfc_options
    ! specify Cd (drag coefficient):

    !$omp parallel do default(shared)   &
    !$omp private(i,j,wlast,n,var,rznt)
    do j=j1,j2
    do i=i1,i2
      cd(i,j) = max(1.0e-10,cnstcd)
      znt(i,j) = 10.0/(exp(karman/sqrt(cd(i,j)))-1.0)
      rznt = 1.0/znt(i,j)
      var = alog((10.0+znt(i,j))*rznt)/alog((za(i,j)+znt(i,j))*rznt)
      u10(i,j) = u1(i,j)*var
      v10(i,j) = v1(i,j)*var
      s10(i,j) = s1(i,j)*var
      ust(i,j) = max( s1(i,j)*karman/alog((za(i,j)+znt(i,j))*rznt) , 1.0e-6 )
    enddo
    enddo

  ELSE  sfc_options
    ! Cd,z0,ust are functions of 10-m windspeed over water ... need to iterate:

    !$omp parallel do default(shared)   &
    !$omp private(i,j,wlast,n,var,rznt)
    do j=j1,j2
    do i=i1,i2
      IF(xland(i,j).gt.1.5)THEN
        !-----------------------------------------
        ! water:  roughness length (z0) is a function of windspeed
        ! use last known z0 for first guess:
        rznt = 1.0/znt(i,j)
        var = alog((10.0+znt(i,j))*rznt)/alog((za(i,j)+znt(i,j))*rznt)
        s10(i,j) = s1(i,j)*var
        wlast = -1000.0
        n = 0
        do while( abs(s10(i,j)-wlast).gt.0.001 )
          n = n + 1
          wlast = s10(i,j)
          IF(cecd.eq.2)THEN
            ! Deacon's formula:  see Rotunno and Emanuel (1987, JAS, p. 547)
            cd(i,j) = 1.1e-3+(4.0e-5*s10(i,j))
          ELSEIF(cecd.eq.3)THEN
            ! based on Fairall et al (2003, JClim) at low wind speeds
            ! based on Donelan et al (2004, GRL) at high wind speeds
            cd(i,j) = dcd1+(s10(i,j)-dwsp1)*dfac
            cd(i,j) = min(cd(i,j),dcd2)
            cd(i,j) = max(cd(i,j),dcd1)
          ENDIF
          znt(i,j) = max( 1.0e-20 , 10.0/(exp(karman/sqrt(cd(i,j)))-1.0) )
          rznt = 1.0/znt(i,j)
          var = alog((10.0+znt(i,j))*rznt)/alog((za(i,j)+znt(i,j))*rznt)
          s10(i,j) = s1(i,j)*var
          if(n.gt.10) print *,'  getcecd:  myid,n,s10 = ',myid,n,s10(i,j)
          if(n.gt.20)then
            call stopcm1
          endif
        enddo
!!!        nmax = max(nmax,n)
        ! end water
        !-----------------------------------------
      ELSE
        !-----------------------------------------
        ! land:  roughness length (z0) is specified ... no need to iterate
        rznt = 1.0/znt(i,j)
        cd(i,j) = ( karman/alog((10.0+znt(i,j))*rznt) )**2
        var = alog((10.0+znt(i,j))*rznt)/alog((za(i,j)+znt(i,j))*rznt)
        s10(i,j) = s1(i,j)*var
        ! end land
        !-----------------------------------------
      ENDIF
      u10(i,j) = u1(i,j)*var
      v10(i,j) = v1(i,j)*var
      ust(i,j) = max( s1(i,j)*karman/alog((za(i,j)+znt(i,j))*rznt) , 1.0e-6 )
    enddo
    enddo

  ENDIF  sfc_options

!-----------------------------------------------------------------------
!  ch,cq:

    IF(isfcflx.eq.1)THEN
      ! surface fluxes of heat/moisture are included ... get exchange coefficients:

      !$omp parallel do default(shared)   &
      !$omp private(i,j,wlast,n,var,rznt)
      do j=j1,j2
      do i=i1,i2
        IF(xland(i,j).gt.1.5)THEN
          !-------
          ! water:
          IF(cecd.eq.1)THEN
            ! constant value (from namelist.input):
            ch(i,j) = cnstce
            cq(i,j) = cnstce
          ELSEIF(cecd.eq.2)THEN
            ! Deacon's formula:  see Rotunno and Emanuel (1987, JAS, p. 547)
            ch(i,j) = 1.1e-3+(4.0e-5*s10(i,j))
            cq(i,j) = 1.1e-3+(4.0e-5*s10(i,j))
          ELSEIF(cecd.eq.3)THEN
            ! Constant, based on Drennan et al. (2007, JAS, p. 1103)
            ch(i,j) = 1.20e-3
            cq(i,j) = 1.20e-3
          ENDIF
          ! end water
          !-------
        ELSE
          !-------
          ! land ... just set Ce to Cd (for now):
          ch(i,j) = cd(i,j)
          cq(i,j) = cd(i,j)
          ! end land
          !-------
        ENDIF
      enddo
      enddo

    ENDIF

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

    if( use_avg_sfc )then

      do j=1,nj
      do i=1,ni
        u10(i,j) = u10(1,1)
        v10(i,j) = v10(1,1)
        s10(i,j) = s10(1,1)
        znt(i,j) = znt(1,1)
        ust(i,j) = ust(1,1)
        cd(i,j) = cd(1,1)
        ch(i,j) = ch(1,1)
        cq(i,j) = cq(1,1)
      enddo
      enddo

      u1(1,1) = usave
      v1(1,1) = vsave
      s1(1,1) = ssave

    endif

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

!!!      print *,'  nmax = ',nmax

      if(timestats.ge.1) time_sfcphys=time_sfcphys+mytime()

      end subroutine getcecd

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine sfcflux(dt,ruh,xf,rvh,pi0s,ch,cq,pi0,thv0,th0,rf0,tsk,thflux,qvflux,mavail,   &
                         rho,rf,u1,v1,s1,ppi,tha,qva,qbsfc,psfc,u10,v10,s10,qsfc,znt,rtime)
      use input
      use constants
      use cm1libs , only : rslf
      implicit none

      real, intent(in) :: dt
      real, intent(in), dimension(ib:ie) :: ruh
      real, intent(in), dimension(ib:ie+1) :: xf
      real, intent(in), dimension(jb:je) :: rvh
      real, intent(in), dimension(ib:ie,jb:je) :: pi0s,ch,cq
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: pi0,thv0,th0,rf0
      real, intent(in), dimension(ib:ie,jb:je) :: tsk
      real, intent(inout), dimension(ib:ie,jb:je) :: psfc,thflux,qvflux
      real, intent(in), dimension(ibl:iel,jbl:jel) :: mavail
      real, intent(in), dimension(ib:ie,jb:je) :: u1,v1,s1
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: rho,rf,ppi,tha
      real, intent(in), dimension(ibm:iem,jbm:jem,kbm:kem) :: qva
      double precision, intent(inout) :: qbsfc
      real, intent(inout), dimension(ibl:iel,jbl:jel) :: u10,v10,s10,qsfc
      real, intent(inout), dimension(ib:ie,jb:je) :: znt
      real, intent(in) :: rtime

      integer :: i,j
      real :: pisfc,qvsat,tem,shf
      double precision, dimension(nj) :: bud1

!-----------------------------------------------------------------------
!
!  This subroutine calculates surface fluxes of heat and moisture.
!
!-----------------------------------------------------------------------

  !c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c!

  setting_flux:  &
  IF( set_flx.eq.1 )THEN

    ! specified heat fluxes:

    !$omp parallel do default(shared)   &
    !$omp private(i,j)
    DO j=1,nj
    do i=1,ni
      thflux(i,j) = cnst_shflx
    enddo
    ENDDO

    IF( imoist.eq.1 )THEN

      !$omp parallel do default(shared)   &
      !$omp private(i,j)
      DO j=1,nj
      do i=1,ni
        qvflux(i,j) = cnst_lhflx
      enddo
      ENDDO

    ENDIF

  !c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c!

  ELSE

    ! normal code:

!$omp parallel do default(shared)   &
!$omp private(i,j,pisfc,qvsat)
    DO j=1,nj

      !  sensible heat flux:
      do i=1,ni
        pisfc = (psfc(i,j)*rp00)**rovcp
        thflux(i,j)=ch(i,j)*s10(i,j)*(tsk(i,j)/pisfc-th0(i,j,1)-tha(i,j,1))
      enddo

      !  latent heat flux:
      IF(imoist.eq.1)THEN
        do i=1,ni
          qvsat=rslf(psfc(i,j),tsk(i,j))
          qsfc(i,j)=qvsat
          qvflux(i,j)=cq(i,j)*s10(i,j)*(qvsat-qva(i,j,1))*mavail(i,j)
        enddo
      ENDIF

    ENDDO

  ENDIF  setting_flux

  !c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c-c!

  !-------------------------------------------------------------------------
    ! for certain cases, turn off surface heat flux after specified time:

    IF( testcase.eq.2 )THEN

      ! shear-driven PBL:
      ! heat flux for t < 3000 s to spin-up turbulence
      ! (Moeng and Sullivan, 1994, JAS)
      if( rtime.gt.3000.0 ) thflux = 0.0

    ELSEIF( testcase.eq.6 )THEN

      ! hurricane PBL  (heat flux for 3000 s to spin-up turbulence)
      ! (Bryan et al, 2017, BLM)
      if( rtime.gt.3000.0 ) thflux = 0.0

    ELSEIF( testcase.eq.7 )THEN

      ! Replace normal fluxes ... use lowest-model-level wind speed:
      ! (VanZanten et al 2011, JAMES)

      !$omp parallel do default(shared)   &
      !$omp private(i,j,pisfc,qvsat)
      DO j=1,nj
        do i=1,ni
          pisfc = (psfc(i,j)*rp00)**rovcp
          thflux(i,j) = 0.001094*s1(i,j)*(tsk(i,j)/pisfc-th0(i,j,1)-tha(i,j,1))
        enddo
        IF(imoist.eq.1)THEN
          do i=1,ni
            qvsat=rslf(psfc(i,j),tsk(i,j))
            qsfc(i,j)=qvsat
            qvflux(i,j) = 0.001133*s1(i,j)*(qvsat-qva(i,j,1))*mavail(i,j)
          enddo
        ENDIF
      ENDDO

    ENDIF

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

    IF(imoist.eq.1)THEN
      ! some budget calculations:
!$omp parallel do default(shared)   &
!$omp private(i,j)
      do j=1,nj
        bud1(j)=0.0d0
        do i=1,ni
          bud1(j)=bud1(j)+qvflux(i,j)*ruh(i)*rvh(j)*rf(i,j,1)
        enddo
      enddo
      tem = dt*dx*dy
      do j=1,nj
        qbsfc=qbsfc+bud1(j)*tem
      enddo
    ENDIF

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

      if(timestats.ge.1) time_sfcphys=time_sfcphys+mytime()

      end subroutine sfcflux

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine sfcdiags(tsk,thflux,qvflux,cd,ch,cq,u1,v1,s1,wspd,        &
                          xland,psfc,qsfc,u10,v10,hfx,qfx,cda,znt,gz1oz0,  &
                          psim,psih,br,zol,mol,hpbl,dsxy,th2,t2,q2,fm,fh,  &
                          zs,za,pi0s,pi0,th0,ppi,tha,rho,rf,qa)
      use input
      use constants
      use cm1libs , only : rslf
      implicit none

      real, intent(in), dimension(ib:ie,jb:je) :: tsk,thflux,qvflux,   &
                                                  cd,ch,cq,u1,v1,s1,xland,psfc
      real, intent(inout), dimension(ibl:iel,jbl:jel) :: qsfc,u10,v10,hfx,qfx,wspd, &
                                    cda,gz1oz0,psim,psih,br,zol,mol,hpbl,dsxy,th2,t2,q2,fm,fh
      real, intent(inout), dimension(ib:ie,jb:je) :: znt
      real, intent(in), dimension(ib:ie,jb:je) :: zs
      real, intent(in), dimension(ib:ie,jb:je) :: za
      real, intent(in), dimension(ib:ie,jb:je) :: pi0s
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: pi0,th0,ppi,tha,rho,rf,qa

      integer :: i,j
      real :: pisfc,thgb,thx,thvx,tskv,govrth,dthvdz,vconv,vsgd,dthvm,   &
              val,fluxc,rznt

      REAL    , PARAMETER ::  VCONVC=1.
      REAL    , PARAMETER ::  CZO=0.0185
      REAL    , PARAMETER ::  OZO=1.59E-5
      real :: ep1

      EP1 = rv/rd - 1.0

      ! surface layer diagnostics:

!$omp parallel do default(shared)   &
!$omp private(i,j,pisfc,thgb,thx,thvx,tskv,govrth,dthvdz,vconv,vsgd,   &
!$omp dthvm,val,fluxc,rznt)
      do j=1,nj
      do i=1,ni
        pisfc = (psfc(i,j)*rp00)**rovcp
        thgb = tsk(i,j)/pisfc
        thx = th0(i,j,1)+tha(i,j,1)
        thvx = thx*(1.+EP1*qa(i,j,1))
        qsfc(i,j) = rslf(psfc(i,j),tsk(i,j))
        tskv = thgb*(1.0+ep1*qsfc(i,j))
        govrth = g/thx
        rznt = 1.0/znt(i,j)
        gz1oz0(i,j) = alog((za(i,j)+znt(i,j))*rznt)
        DTHVDZ = THVX-TSKV
        ! cm1r18:  use same formulation over land and water:
!!!        if (xland(i,j).lt.1.5) then
          ! land:
          fluxc = max(thflux(i,j) + ep1*tskv*qvflux(i,j),0.)
          VCONV = vconvc*(g/tsk(i,j)*hpbl(i,j)*fluxc)**.33
!!!        else
!!!          ! ocean:
!!!          IF(-DTHVDZ.GE.0)THEN
!!!            DTHVM=-DTHVDZ
!!!          ELSE
!!!            DTHVM=0.
!!!          ENDIF
!!!          VCONV = 2.*SQRT(DTHVM)
!!!        endif
! Mahrt and Sun low-res correction
        VSGD = 0.32 * (max(dsxy(i,j)/5000.-1.,0.))**.33
        wspd(i,j) = sqrt( s1(i,j)*s1(i,j) + vconv*vconv + vsgd*vsgd )
        wspd(i,j) = max(0.1,wspd(i,j))
        br(i,j) = govrth*za(i,j)*DTHVDZ/(wspd(i,j)**2)
        hfx(i,j) = thflux(i,j)*cp*rf(i,j,1)
        qfx(i,j) = qvflux(i,j)*rf(i,j,1)
        cda(i,j) = cd(i,j)
        ! impose neutral sfc layer:
        psim(i,j) = 0.0
        psih(i,j) = 0.0
        zol(i,j) = 0.0
        mol(i,j) = 0.0
        fm(i,j) = GZ1OZ0(i,j)
        fh(i,j) = GZ1OZ0(i,j)
        ! get 2-m th/q/t:
        val = alog((2.0+znt(i,j))*rznt)/alog((za(i,j)+znt(i,j))*rznt)
        th2(i,j) = thgb+(thx-thgb)*val
        q2(i,j) = qsfc(i,j)+(qa(i,j,1)-qsfc(i,j))*val
        t2(i,j) = th2(i,j)*pisfc
      enddo
      enddo

      if(timestats.ge.1) time_sfcphys=time_sfcphys+mytime()
      end subroutine sfcdiags

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine gethpbl(zh,th0,tha,qa,hpbl)
      use input
      use constants
      implicit none

      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: zh,th0,tha,qa
      real, intent(inout), dimension(ibl:iel,jbl:jel) :: hpbl

      integer :: i,j,kk
      real :: thx,thvx,thv,thvlast,thcrit,ep1

      EP1 = rv/rd - 1.0

      ! (NEEDED BY SFCLAY ... THIS IS A ROUGH ESTIMATE ONLY)
      ! (ONLY NEEDED WHEN IPBL=0)
      ! (USE WITH CAUTION)
      ! extraordinarily simple calculation:  define pbl depth as
      ! level where thv is first greater than thv at lowest model level
      ! 110104:  add 0.5 K, for the sake of slightly stable PBLs

!$omp parallel do default(shared)   &
!$omp private(i,j,kk,thx,thvx,thv,thvlast,thcrit)
      do j=1,nj
      do i=1,ni
        hpbl(i,j) = 0.0
        kk = 1
        thx = th0(i,j,1)+tha(i,j,1)
        thvx = thx*(1.+EP1*qa(i,j,1))
        thvlast = thvx
        thcrit = thvx+0.5
        do while( hpbl(i,j).lt.1.0e-12 .and. kk.lt.nk )
          kk = kk + 1
          thv = (th0(i,j,kk)+tha(i,j,kk))*(1.0+EP1*qa(i,j,kk))
          if( thv.ge.thcrit )then
            hpbl(i,j) = zh(i,j,kk-1)+(zh(i,j,kk)-zh(i,j,kk-1))   &
                                    *(thcrit-thvlast)/(thv-thvlast)
          endif
          thvlast = thv
        enddo
        if( kk.gt.(nk-1) .or. hpbl(i,j).lt.1.0e-12 ) hpbl(i,j) = 0.0
      enddo
      enddo

      if(timestats.ge.1) time_sfcphys=time_sfcphys+mytime()
      end subroutine gethpbl

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine gethpbl2(psfc,qsfc,thflux,qvflux,ust,tsk,zh,th0,tha,qa,ua,va,thv1,thvsfc,govthv,brilast,hpbl,riout)
      use input
      use constants
      use cm1libs , only : rslf
      implicit none

      real, intent(in), dimension(ib:ie,jb:je) :: psfc
      real, intent(inout), dimension(ibl:iel,jbl:jel) :: qsfc
      real, intent(in), dimension(ib:ie,jb:je) :: thflux,qvflux,ust,tsk
      real, intent(in), dimension(ib:ie,jb:je,kb:ke) :: zh,th0,tha,qa
      real, intent(in), dimension(ib:ie+1,jb:je,kb:ke) :: ua
      real, intent(in), dimension(ib:ie,jb:je+1,kb:ke) :: va
      real, intent(inout), dimension(ib:ie,jb:je) :: thv1,thvsfc,govthv,brilast
      real, intent(inout), dimension(ibl:iel,jbl:jel) :: hpbl
      real, intent(inout), dimension(ib:ie,jb:je,kb:ke) :: riout

      integer :: i,j,k,ntot,nloop
      real :: thv,bri,uavg,vavg,thvflux,tmpprt,thf1,thf2,gamfac,ws
      logical :: doit

      ! critical (ie, threshold) bulk Richardson number:
      real, parameter :: bric = 0.125

      ! cm1r19.6: follows formulation in YSU scheme.

      if( imoist.eq.1 )then
        do j=1,nj
        do i=1,ni
          qsfc(i,j) = rslf(psfc(i,j),tsk(i,j))
        enddo
        enddo
      else
        do j=1,nj
        do i=1,ni
          qsfc(i,j) = 0.0
        enddo
        enddo
      endif

      do j=1,nj
      do i=1,ni
        thv = (th0(i,j,1)+tha(i,j,1))*(1.0+repsm1*qa(i,j,1))
        govthv(i,j) = g/thv
      enddo
      enddo

    hloop:  &
    DO nloop = 1 , 2

      k = 1

    if( nloop.eq.1 )then
      do j=1,nj
      do i=1,ni
        thv = (th0(i,j,1)+tha(i,j,1))*(1.0+repsm1*qa(i,j,1))
        uavg = 0.5*(ua(i,j,k)+ua(i+1,j,k))
        vavg = 0.5*(va(i,j,k)+va(i,j+1,k))
        thvsfc(i,j) = tsk(i,j)*((p00/psfc(i,j))**rovcp)*(1.0+repsm1*qsfc(i,j))
        brilast(i,j) = govthv(i,j)*zh(i,j,1)*(thv-thvsfc(i,j))/max( 1.0 , uavg**2 + vavg**2 )
        thv1(i,j) = thv
      enddo
      enddo
    else
      do j=1,nj
      do i=1,ni
        thv = (th0(i,j,1)+tha(i,j,1))*(1.0+repsm1*qa(i,j,1))
        uavg = 0.5*(ua(i,j,k)+ua(i+1,j,k))
        vavg = 0.5*(va(i,j,k)+va(i,j+1,k))
        thf1 = (1.0+repsm1*qa(i,j,1))
        thf2 = repsm1*(th0(i,j,1)+tha(i,j,1))
        thvflux = max( 0.0 , thf1*thflux(i,j) + thf2*qvflux(i,j) )
        ws = ( ust(i,j)**3 + 8.0*karman*0.5*govthv(i,j)*thvflux*hpbl(i,j) )**0.3333333
        ws = min( ws , ust(i,j)*16.0 )
        ws = max( ws , ust(i,j)/5.0 )
        ws = max( ws , 0.0001 )
        gamfac = 6.8/ws
        tmpprt = max( 0.0 , min( gamfac*thflux(i,j) , 3.0 )*thf1    &
                           +min( gamfac*qvflux(i,j) , 2.0e-3 )*thf2 )
        thv1(i,j) = thv + tmpprt
        brilast(i,j) = govthv(i,j)*zh(i,j,1)*(thv-thvsfc(i,j))/max( 1.0 , uavg**2 + vavg**2 )
        riout(i,j,1) = brilast(i,j)
        riout(i,j,nk) = thv1(i,j)
        riout(i,j,nk-1) = ws
      enddo
      enddo
    endif

      doit = .true.
      ntot = 0

      hpbl = -1.0

      do k=2,nk
        if( doit )then
          do j=1,nj
          do i=1,ni
            if( hpbl(i,j).lt.0.0 )then
              thv = (th0(i,j,k)+tha(i,j,k))*(1.0+repsm1*qa(i,j,k))
              uavg = 0.5*(ua(i,j,k)+ua(i+1,j,k))
              vavg = 0.5*(va(i,j,k)+va(i,j+1,k))
              bri = govthv(i,j)*zh(i,j,k)*(thv-thv1(i,j))/max( 1.0 , uavg**2 + vavg**2 )
              riout(i,j,k) = bri
              if( bri.gt.bric )then
                if( abs(bri-brilast(i,j)).lt.0.000001 .or. k.eq.1 )then
                  hpbl(i,j) = zh(i,j,k)
                else
                  hpbl(i,j) = zh(i,j,k-1)+(zh(i,j,k)-zh(i,j,k-1))  &
                                         *(bric-brilast(i,j))  &
                                         /(bri-brilast(i,j))
                endif
                ntot = ntot+1
              endif
              brilast(i,j) = bri
            endif
          enddo
          enddo
        endif
        if( ntot.eq.(ni*nj) ) doit = .false.
      enddo

    ENDDO  hloop

      if(timestats.ge.1) time_sfcphys=time_sfcphys+mytime()
      end subroutine gethpbl2

!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

  END MODULE sfcphys_module


text:
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!

C RCS file, release, date & time of last delta, author, state, [and locker]
C $Header: /project/work/rep/STENEX/src/se_snl/se_comm_info_ext.f,v 1.1 2004/03/26 16:16:47 yoj Exp $

C what(1) key, module and SID; SCCS file; date and time of last delta:
C %W% %P% %G% %U%

C --------------------------------------------------------------------------
C Note: all these variables with prefix se_ are for stencil exchange library
C       only
C
C to define communication info variables:
C
C   se_ngb_pe    -- an array to indicate a communication with a certain
C                   processor is required base upon near-neighbour
C                   communication pattern: -1 denotes no communication is
C                   needed, and a non -1 number denotes processor number with
C                   which communication is formed
C   se_numdim    -- dimensionality of a data structure which requires
C                   communication
C   se_decompstr -- indicator of which dimenion(s) of data is/are decomposed,
C                   0 (not decomposed), 1 (decomposed)
C --------------------------------------------------------------------------

        module se_comm_info_ext

          integer :: se_ngb_pe(8)
          integer :: se_numdim
          character (len=10) :: se_decompstr

          integer :: se_twoway_npcol, se_twoway_nprow

        end module se_comm_info_ext

ext: f

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed


text:
      subroutine ire_mk(ctfm,ctfmp,n_A,n_Athth,n_Axx,n_B,n_K,n_Sx,n_psi,
     &x,Nx,hx,myzero,phys_bdy,res)
      implicit none
      integer i
      integer Nx
      real*8 hx
      real*8 myzero
      real*8 ctfm(Nx)
      real*8 ctfmp(Nx)
      real*8 n_A(Nx)
      real*8 n_Athth(Nx)
      real*8 n_Axx(Nx)
      real*8 n_B(Nx)
      real*8 n_K(Nx)
      real*8 n_Sx(Nx)
      real*8 n_psi(Nx)
      real*8 x(Nx)
      integer phys_bdy(2)
      real*8 res(Nx)
      real*8 qb
      if (phys_bdy(1) .eq. 1) then
      do i=1, 1, 1
      qb = -0.1D1 * n_psi(i) ** 6 / n_A(i) * n_Sx(i)
      res(i)=qb
      end do
      endif
      do i=2, Nx-1, 1
      qb = -0.1D1 * n_psi(i) ** 6 / n_A(i) * n_Sx(i) + (0.2D1 * n_psi(i)
     # ** 6 / n_A(i) ** 2 * n_Axx(i) - 0.2D1 * n_psi(i) ** 6 / n_A(i) /
     #n_B(i) * n_Athth(i)) / ctfm(i) + (-0.5000000000000000D0 * n_psi(i)
     # ** 6 / n_A(i) ** 3 * n_Axx(i) * (-0.1D1 * n_A(i - 1) + n_A(i + 1)
     #) + 0.5000000000000000D0 * n_psi(i) ** 6 / n_A(i) ** 2 * (-0.1D1 *
     # n_Axx(i - 1) + n_Axx(i + 1)) + 0.5000000000000000D0 * n_psi(i) **
     # 6 / n_A(i) ** 2 / n_B(i) * (-0.1D1 * n_B(i - 1) + n_B(i + 1)) * n
     #_Axx(i) - 0.5000000000000000D0 * n_psi(i) ** 6 / n_A(i) / n_B(i) *
     #* 2 * (-0.1D1 * n_B(i - 1) + n_B(i + 1)) * n_Athth(i) + 0.33333333
     #33333333D0 * n_psi(i) ** 6 / n_A(i) * (n_K(i - 1) - 0.1D1 * n_K(i
     #+ 1)) + 0.3D1 * n_psi(i) ** 5 / n_A(i) ** 2 * (-0.1D1 * n_psi(i -
     #1) + n_psi(i + 1)) * n_Axx(i)) / ctfmp(i) / hx
      res(i)=qb
      end do
      if (phys_bdy(2) .eq. 1) then
      do i=Nx, Nx, 1
      qb = myzero * x(i)
      res(i)=qb
      end do
      endif
      END


text:
module Mod_memor
!-----------------------------------------------------------------------
!
!    This module contains memory allocation class
!
!-----------------------------------------------------------------------
  use typre
  implicit none
  private
  public :: MemoryMan, TotalMemo, TotalMax!, alloc

  integer(8) :: TotalMemo = 0, TotalMax

  type MemoryMan
      private
      !public :: alloc
      integer(8) :: MaxMemo = 0
      integer(8) :: CurrentMemo = 0
      integer(ip) :: lumem = 0
      integer(ip) :: luerr = 0

contains
     !private
     !public :: alloc
     procedure :: allocrp1
     procedure :: allocrp2
     procedure :: allocrp3
     procedure :: allocrp4
     procedure :: allocip1
     procedure :: allocip2
     procedure :: allocip3          ! integers
     procedure :: pallocip1
     procedure :: pallocip2
     procedure :: pallocip3          ! integers pointers
     procedure :: pallocip4
     procedure :: alloci1p1                        ! i1p types
     procedure :: alloci2p1                        ! i2p types
     procedure :: allocr1p1
     procedure :: allocr1p2                 ! r1p types
     procedure :: allocr2p1
     procedure :: allocr2p2                 ! r2p types
     procedure :: allocr3p1
     procedure :: allocr4p1
     procedure :: allocr3p2                   ! r3p types
     procedure :: pallocrp1
     procedure :: pallocrp2
     procedure :: pallocrp3
     procedure :: pallocrp4
     procedure :: alloclg1
     procedure :: alloclg2
     procedure :: allocch1
!      procedure :: allocPostprFile

     generic :: alloc => allocrp1, allocrp2,allocrp3,allocrp4, &
            allocip1,allocip2,allocip3,&          ! integers
            alloci1p1,&                        ! i1p types
            alloci2p1,&                        ! i2p types
            allocr1p1,allocr1p2,&                 ! r1p types
            allocr2p1,allocr2p2,&                 ! r2p types
            allocr3p1,allocr3p2,&                   ! r3p types
            alloclg1,alloclg2,& ! logicals
            allocch1,allocr4p1 ! characters

     procedure :: allocObj

     generic :: palloc => pallocip1,pallocip2,pallocip3,pallocip4,&        ! integers
                          pallocrp1,pallocrp2,pallocrp3,pallocrp4,&
                          palloclp1

     procedure :: deallocrp1
     procedure :: deallocrp2
     procedure :: deallocrp3
     procedure :: deallocrp4
     procedure :: deallocip1
     procedure :: deallocip2
     procedure :: deallocip3          ! integers
     procedure :: pdeallocip1
     procedure :: pdeallocip2
     procedure :: pdeallocip3          ! integers
     procedure :: pdeallocip4
     procedure :: dealloci1p1                        ! i1p types
     procedure :: dealloci2p1                        ! i2p types
     procedure :: deallocr1p1
     procedure :: deallocr1p2                 ! r1p types
     procedure :: deallocr2p1
     procedure :: deallocr2p2                 ! r2p types
     procedure :: deallocr3p1
     procedure :: deallocr4p1

     procedure :: deallocr3p2                   ! r3p types
     procedure :: pdeallocrp1
     procedure :: pdeallocrp2
     procedure :: pdeallocrp3
     procedure :: pdeallocrp4

     procedure :: palloclp1
     procedure :: pdealloclp1

     procedure :: dealloclg1
     procedure :: dealloclg2

     procedure :: deallocch1

     generic :: dealloc => deallocrp1, deallocrp2,deallocrp3,deallocrp4, &
         deallocip1,deallocip2,deallocip3,&          ! integers
         dealloci1p1,&                        ! i1p types
         dealloci2p1,&                        ! i2p types
         deallocr1p1,deallocr1p2,&            ! r1p types
         deallocr2p1,deallocr2p2,&            ! r2p types
         deallocr3p1,deallocr3p2,&            ! r3p types
         dealloclg1,dealloclg2,&              ! logicals
         deallocch1,deallocr4p1                           ! characters

     generic :: pdealloc => pdeallocip1,pdeallocip2,pdeallocip3,pdeallocip4, &
                            pdeallocrp1,pdeallocrp2,pdeallocrp3,pdeallocrp4,&
                            pdealloclp1

     procedure :: deallocObj

     !Realloc
     procedure :: reallocip1
     procedure :: reallocrp1
     procedure :: reallocrp2
     procedure :: reallocip2
     procedure :: reallocip3
     procedure :: reallocrp3
     procedure :: reallocrp4

     generic :: realloc => reallocip1,reallocrp1,reallocrp2,reallocip2,reallocip3,reallocrp3,reallocrp4

     procedure :: preallocip1
     procedure :: preallocrp1
     procedure :: preallocrp2
     generic :: prealloc => preallocip1, preallocrp1,preallocrp2

     procedure :: init
     procedure :: GetValue
     procedure :: GetLumem

  end type

contains

  subroutine init(a,lumem,luerr)
    implicit none
    class(MemoryMan) :: a
    integer(ip), intent(in) :: lumem, luerr

    a%MaxMemo = 0
    a%CurrentMemo = 0
    a%lumem = lumem
    a%luerr = luerr
  end subroutine

  subroutine GetValue(a,Current,MaxMemo,TotalMem,TotalMa)
    implicit none
    class(MemoryMan) :: a
    integer(8) :: Current, Maxmemo, TotalMem,TotalMa

    MaxMemo = a%MaxMemo
    Current = a%CurrentMemo
    TotalMem = TotalMemo
    TotalMa  = TotalMax
  end subroutine

  subroutine memctr(CurrentMemo,MaxMemo,lbyts,vanam,lumem,vacal)
      use typre
      implicit none
      integer(ip), intent(in)    :: lbyts
      character(*), intent(in)   :: vanam
      integer(ip), intent(in)    :: lumem
      character(*), intent(in)   :: vacal       ! Calling routine
      integer(8), intent(inout) :: MaxMemo,CurrentMemo
      character(6)               :: lbyte, lbyte2, lbyte3
      real(rp) :: rMaxMemo, rCurrentMemo,rTotalMemo,rMemory,rdummy

      CurrentMemo = CurrentMemo + lbyts
      TotalMemo   = TotalMemo   + lbyts
      MaxMemo = max(MaxMemo,CurrentMemo)
      TotalMax = max(TotalMemo,TotalMax)

      if(lumem>0) then
         call Mem2String(CurrentMemo,rCurrentMemo,rdummy,lbyte)
         call Mem2String(TotalMemo,rTotalMemo,rdummy,lbyte2)
         call Mem2String(int(lbyts,8),rMemory,rdummy,lbyte3)

         write(lumem,'(a12,2x,a8,f8.1,1x,a6,2x,a13,f8.1,1x,a6,2x,a12,2x,f8.1,1x,a6,1x,a10)') &
            trim(vanam),'Memory: ', rMemory,lbyte3,'Module memo: ', rCurrentMemo,lbyte, 'Total memo: ', rTotalMemo,lbyte2,vacal
      else
         !write(*,*) 'hey check this: betterh than runend, some buffers do not have output channel'

      endif

end subroutine memctr

!-------------------------------------------------------------------
! General Objects (unknown but fixed size)
!-------------------------------------------------------------------

  subroutine allocObj(a,istat,vanam,vacal,lbyts)
    implicit none
    class(MemoryMan)            :: a
    integer(ip)                 :: istat
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine
    integer(ip), intent(in)     :: lbyts

    if(istat==0) then
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

  end subroutine

  subroutine deallocObj(a,istat,vanam,vacal,elbyts)
    implicit none
    class(MemoryMan)            :: a
    integer(ip)                 :: istat
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine
    integer(ip), intent(in)     :: elbyts

    integer(ip) :: lbyts

    if(istat==0) then
      lbyts=-elbyts
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

!-----------------------------------------------------------------------
! Real matrices
!---------------------------------------------------------------------

  subroutine allocrp1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    real(rp), allocatable           :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1
      varia=0.0_rp
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocrp1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    real(rp), allocatable           :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine reallocrp1(a,newndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip), intent(in)     :: newndim1

    real(rp), allocatable           :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    real(rp), allocatable :: aux_varia(:)

    integer(ip) :: oldndim1

    allocate(aux_varia(newndim1),stat=istat)
    if(istat==0) then
      lbyts=rp*newndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    oldndim1 = size(varia)

    aux_varia(1:min(newndim1,oldndim1)) = varia(1:min(newndim1,oldndim1))
    call move_alloc(aux_varia,varia)

    lbyts = -rp*oldndim1
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine allocrp2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    real(rp), allocatable           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2
      varia=0.0_rp
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocrp2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    real(rp), allocatable           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine reallocrp2(a,newndim1,newndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip), intent(in)     :: newndim1,newndim2

    real(rp), allocatable           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    real(rp), allocatable :: aux_varia(:,:)

    integer(ip) :: oldndim1,oldndim2

    allocate(aux_varia(newndim1,newndim2),stat=istat)
    if(istat==0) then
      lbyts=rp*newndim1*newndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    oldndim1 = size(varia,1)
    oldndim2 = size(varia,2)

    lbyts = -rp*size(varia)
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    aux_varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2)) = varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2))
    call move_alloc(aux_varia,varia)

  end subroutine

  subroutine allocrp3(a,ndim1,ndim2,ndim3,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3       ! Size of the variable
    real(rp), allocatable           :: varia(:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2,ndim3),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2*ndim3
      varia=0.0_rp
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocrp3(a,ndim1,ndim2,ndim3,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3       ! Size of the variable
    real(rp), allocatable           :: varia(:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2*ndim3
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine reallocrp3(a,newndim1,newndim2,newndim3,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip), intent(in)     :: newndim1,newndim2,newndim3

    real(rp), allocatable           :: varia(:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    real(rp), allocatable :: aux_varia(:,:,:)

    integer(ip) :: oldndim1,oldndim2,oldndim3

    allocate(aux_varia(newndim1,newndim2,newndim3),stat=istat)
    if(istat==0) then
      lbyts=rp*newndim1*newndim2*newndim3
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    oldndim1 = size(varia,1)
    oldndim2 = size(varia,2)
    oldndim3 = size(varia,3)

    lbyts = -rp*oldndim1*oldndim2*oldndim3
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    aux_varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2),1:min(newndim3,oldndim3)) = varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2),1:min(newndim3,oldndim3))
    call move_alloc(aux_varia,varia)

  end subroutine

  subroutine allocrp4(a,ndim1,ndim2,ndim3,ndim4,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3,ndim4       ! Size of the variable
    real(rp), allocatable           :: varia(:,:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2,ndim3,ndim4),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2*ndim3*ndim4
      varia=0.0_rp
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocrp4(a,ndim1,ndim2,ndim3,ndim4,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3,ndim4       ! Size of the variable
    real(rp), allocatable           :: varia(:,:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2*ndim3*ndim4
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine reallocrp4(a,newndim1,newndim2,newndim3,newndim4,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip), intent(in)     :: newndim1,newndim2,newndim3,newndim4

    real(rp), allocatable           :: varia(:,:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    real(rp), allocatable :: aux_varia(:,:,:,:)

    integer(ip) :: oldndim1,oldndim2,oldndim3,oldndim4

    allocate(aux_varia(newndim1,newndim2,newndim3,newndim4),stat=istat)
    if(istat==0) then
      lbyts=rp*newndim1*newndim2*newndim3*newndim4
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    oldndim1 = size(varia,1)
    oldndim2 = size(varia,2)
    oldndim3 = size(varia,3)
    oldndim4 = size(varia,4)

    lbyts = -rp*oldndim1*oldndim2*oldndim3*oldndim4
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    aux_varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2),1:min(newndim3,oldndim3),1:min(newndim4,oldndim4)) = varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2),1:min(newndim3,oldndim3),1:min(newndim4,oldndim4))
    call move_alloc(aux_varia,varia)

  end subroutine

!-----------------------------------------------------------------------
! Integer Matrices
!-----------------------------------------------------------------------

  subroutine allocip1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    integer(ip), allocatable           :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=ip*ndim1
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocip1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    integer(ip), allocatable    :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ip*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine reallocip1(a,newndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip), intent(in)     :: newndim1

    integer(ip), allocatable           :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    integer(ip), allocatable :: aux_varia(:)

    integer(ip) :: oldndim1

    allocate(aux_varia(newndim1),stat=istat)
    if(istat==0) then
      lbyts=ip*newndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    oldndim1 = size(varia)

    aux_varia(1:min(newndim1,oldndim1)) = varia(1:min(newndim1,oldndim1))
    call move_alloc(aux_varia,varia)

    lbyts = -ip*oldndim1
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine allocip2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    integer(ip), allocatable           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2),stat=istat)
    if(istat==0) then
      lbyts=ip*ndim1*ndim2
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocip2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    integer(ip), allocatable           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ip*ndim1*ndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine reallocip2(a,newndim1,newndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip), intent(in)     :: newndim1,newndim2

    integer(ip), allocatable           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    integer(ip), allocatable :: aux_varia(:,:)

    integer(ip) :: oldndim1,oldndim2

    allocate(aux_varia(newndim1,newndim2),stat=istat)
    if(istat==0) then
      lbyts=ip*newndim1*newndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    oldndim1 = size(varia,1)
    oldndim2 = size(varia,2)

    aux_varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2)) = varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2))
    call move_alloc(aux_varia,varia)

    lbyts = -ip*oldndim1*oldndim2
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine allocip3(a,ndim1,ndim2,ndim3,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3       ! Size of the variable
    integer(ip), allocatable           :: varia(:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2,ndim3),stat=istat)
    if(istat==0) then
      lbyts=ip*ndim1*ndim2*ndim3
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocip3(a,ndim1,ndim2,ndim3,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3       ! Size of the variable
    integer(ip), allocatable           :: varia(:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ip*ndim1*ndim2*ndim3
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine reallocip3(a,newndim1,newndim2,newndim3,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip), intent(in)     :: newndim1,newndim2,newndim3

    integer(ip), allocatable           :: varia(:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    integer(ip), allocatable :: aux_varia(:,:,:)

    integer(ip) :: oldndim1,oldndim2,oldndim3

    allocate(aux_varia(newndim1,newndim2,newndim3),stat=istat)
    if(istat==0) then
      lbyts=ip*newndim1*newndim2*newndim3
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    oldndim1 = size(varia,1)
    oldndim2 = size(varia,2)
    oldndim3 = size(varia,3)

    aux_varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2),1:min(newndim3,oldndim3)) = varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2),1:min(newndim3,oldndim3))
    call move_alloc(aux_varia,varia)

    lbyts = -ip*oldndim1*oldndim2*oldndim3
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine allocip4(a,ndim1,ndim2,ndim3,ndim4,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3,ndim4       ! Size of the variable
    integer(ip), allocatable           :: varia(:,:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2,ndim3,ndim4),stat=istat)
    if(istat==0) then
      lbyts=ip*ndim1*ndim2*ndim3*ndim4
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocip4(a,ndim1,ndim2,ndim3,ndim4,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3,ndim4       ! Size of the variable
    integer(ip), allocatable           :: varia(:,:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ip*ndim1*ndim2*ndim3*ndim4
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

!----------------------------------------------------------------------
! Integer Pointers
!----------------------------------------------------------------------
  subroutine pallocip1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    integer(ip), pointer        :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=ip*ndim1
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pdeallocip1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    integer(ip), pointer        :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ip*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine preallocip1(a,newndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip), intent(in)     :: newndim1

    integer(ip), pointer        :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    integer(ip), allocatable, target :: aux_varia(:)

    integer(ip) :: oldndim1

    allocate(aux_varia(newndim1),stat=istat)
    if(istat==0) then
      lbyts=ip*newndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    oldndim1 = size(varia)

    aux_varia(1:min(newndim1,oldndim1)) = varia(1:min(newndim1,oldndim1))
    deallocate(varia)
    allocate(varia(newndim1))
    varia = aux_varia
    deallocate(aux_varia)

    lbyts = -ip*oldndim1
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pallocip2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    integer(ip), pointer        :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2),stat=istat)
    if(istat==0) then
      lbyts=ip*ndim1*ndim2
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pdeallocip2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    integer(ip), pointer        :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ip*ndim1*ndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pallocip3(a,ndim1,ndim2,ndim3,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3       ! Size of the variable
    integer(ip), pointer        :: varia(:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2,ndim3),stat=istat)
    if(istat==0) then
      lbyts=ip*ndim1*ndim2*ndim3
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pdeallocip3(a,ndim1,ndim2,ndim3,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3       ! Size of the variable
    integer(ip), pointer        :: varia(:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ip*ndim1*ndim2*ndim3
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pallocip4(a,ndim1,ndim2,ndim3,ndim4,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3,ndim4       ! Size of the variable
    integer(ip), pointer        :: varia(:,:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2,ndim3,ndim4),stat=istat)
    if(istat==0) then
      lbyts=ip*ndim1*ndim2*ndim3*ndim4
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pdeallocip4(a,ndim1,ndim2,ndim3,ndim4,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3,ndim4       ! Size of the variable
    integer(ip), pointer        :: varia(:,:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ip*ndim1*ndim2*ndim3*ndim4
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

!----------------------------------------------------------------------
! Real Pointers
!----------------------------------------------------------------------
  subroutine pallocrp1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    real(rp), pointer           :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pdeallocrp1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    real(rp), pointer           :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine preallocrp1(a,newndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip), intent(in)     :: newndim1

    real(rp), pointer           :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    real(rp), allocatable, target :: aux_varia(:)

    integer(ip) :: oldndim1

    allocate(aux_varia(newndim1),stat=istat)
    if(istat==0) then
      lbyts=rp*newndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    oldndim1 = size(varia)

    aux_varia(1:min(newndim1,oldndim1)) = varia(1:min(newndim1,oldndim1))
    deallocate(varia)
    allocate(varia(newndim1))
    varia = aux_varia
    deallocate(aux_varia)

    lbyts = -rp*oldndim1
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pallocrp2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    real(rp), pointer           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pdeallocrp2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    real(rp), pointer           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine preallocrp2(a,newndim1,newndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip), intent(in)     :: newndim1,newndim2

    real(rp), pointer           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    real(rp), pointer :: aux_varia(:,:) => NULL()

    integer(ip) :: oldndim1,oldndim2

    allocate(aux_varia(newndim1,newndim2),stat=istat)
    if(istat==0) then
      lbyts=rp*newndim1*newndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    oldndim1 = size(varia,1)
    oldndim2 = size(varia,2)

    lbyts = -rp*size(varia)
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)

    aux_varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2)) = varia(1:min(newndim1,oldndim1),1:min(newndim2,oldndim2))
    deallocate(varia)
    varia => aux_varia

  end subroutine

  subroutine pallocrp3(a,ndim1,ndim2,ndim3,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3       ! Size of the variable
    real(rp), pointer        :: varia(:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2,ndim3),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2*ndim3
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pdeallocrp3(a,ndim1,ndim2,ndim3,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3       ! Size of the variable
    real(rp), pointer           :: varia(:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2*ndim3
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pallocrp4(a,ndim1,ndim2,ndim3,ndim4,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3,ndim4       ! Size of the variable
    real(rp), pointer           :: varia(:,:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2,ndim3,ndim4),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2*ndim3*ndim4
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pdeallocrp4(a,ndim1,ndim2,ndim3,ndim4,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2,ndim3,ndim4       ! Size of the variable
    real(rp), pointer           :: varia(:,:,:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2*ndim3*ndim4
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

!-----------------------------------------------------------------------
! ip matrices
!---------------------------------------------------------------------

  subroutine alloci1p1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(i1p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,idime

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1
      do idime = 1,ndim1
         varia(idime)%l => NULL()
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine dealloci1p1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(i1p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine alloci1p2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    type(i1p), allocatable      :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,idime,jdime

    allocate(varia(ndim1,ndim2),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2
      do idime = 1,ndim2
         do jdime = 1,ndim1
            varia(jdime,idime)%l => NULL()
         enddo
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine dealloci1p2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    type(i1p), allocatable      :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

!-----------------------------------------------------------------------
! i2p matrices
!---------------------------------------------------------------------

  subroutine alloci2p1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(i2p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,idime

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1
      do idime = 1,ndim1
         varia(idime)%l => NULL()
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine dealloci2p1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(i2p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine alloci2p2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    type(i2p), allocatable      :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,idime,jdime

    allocate(varia(ndim1,ndim2),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2
      do idime = 1,ndim2
         do jdime = 1,ndim1
            varia(jdime,idime)%l => NULL()
         enddo
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine dealloci2p2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    type(i2p), allocatable      :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

!-----------------------------------------------------------------------
! r1p matrices
!---------------------------------------------------------------------

  subroutine allocr1p1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(r1p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,idime

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1
      do idime = 1,ndim1
         varia(idime)%a => NULL()
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocr1p1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(r1p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine allocr1p2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    type(r1p), allocatable      :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,idime,jdime

    allocate(varia(ndim1,ndim2),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2
      do idime = 1,ndim2
         do jdime = 1,ndim1
            varia(jdime,idime)%a => NULL()
         enddo
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocr1p2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    type(r1p), allocatable      :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

!-----------------------------------------------------------------------
! r2p matrices
!---------------------------------------------------------------------

  subroutine allocr2p1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(r2p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,jdime

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1
      do jdime = 1,ndim1
         varia(jdime)%a => NULL()
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocr2p1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(r2p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine allocr2p2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    type(r2p), allocatable      :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,idime,jdime

    allocate(varia(ndim1,ndim2),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2
      do idime = 1,ndim2
         do jdime = 1,ndim1
            varia(jdime,idime)%a => NULL()
         enddo
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocr2p2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    type(r2p), allocatable      :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

!-----------------------------------------------------------------------
! r3p matrices
!---------------------------------------------------------------------

  subroutine allocr3p1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(r3p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,jdime

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1
      do jdime = 1,ndim1
         varia(jdime)%a => NULL()
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocr3p1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(r3p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine allocr4p1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(r4p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,jdime

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1
      do jdime = 1,ndim1
         varia(jdime)%a => NULL()
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocr4p1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    type(r4p), allocatable      :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine allocr3p2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    type(r3p), allocatable      :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts,jdime,idime

    allocate(varia(ndim1,ndim2),stat=istat)
    if(istat==0) then
      lbyts=rp*ndim1*ndim2
      do idime = 1,ndim2
         do jdime = 1,ndim1
            varia(jdime,idime)%a => NULL()
         enddo
      enddo
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocr3p2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    type(r3p), allocatable      :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-rp*ndim1*ndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

!----------------------------------------------------------------------
! logicals
!----------------------------------------------------------------------

  subroutine alloclg1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    logical, allocatable           :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=ndim1
      varia=.false.
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine dealloclg1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    logical, allocatable           :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine alloclg2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    logical, allocatable           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1,ndim2),stat=istat)
    if(istat==0) then
      lbyts=ndim1*ndim2
      varia=.false.
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine dealloclg2(a,ndim1,ndim2,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1,ndim2       ! Size of the variable
    logical, allocatable           :: varia(:,:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ndim1*ndim2
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  !Logical pointers
    subroutine palloclp1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    logical, pointer            :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=1*ndim1
      varia=0_ip
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine pdealloclp1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    logical, pointer            :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-1*ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

!----------------------------------------------------------------------
! characters
!----------------------------------------------------------------------

  subroutine allocch1(a,ndim1,varia,vanam,vacal)
    implicit none

    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    character(*), allocatable :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    allocate(varia(ndim1),stat=istat)
    if(istat==0) then
      lbyts=ndim1
      varia=''
    else
      call memerr(0,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

  subroutine deallocch1(a,ndim1,varia,vanam,vacal)
    implicit none
    class(MemoryMan)            :: a
    integer(ip),  intent(in)    :: ndim1       ! Size of the variable
    character(*), allocatable :: varia(:)
    character(*), intent(in)    :: vanam       ! Variable name
    character(*), intent(in)    :: vacal       ! Calling routine

    integer(ip)                 :: istat,lbyts

    deallocate(varia,stat=istat)
    if(istat==0) then
      lbyts=-ndim1
    else
      call memerr(1,vanam,vacal,a%luerr)
    end if
    call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
  end subroutine

! !---------------------------------------------------------------------------
! ! ListenFile
! !---------------------------------------------------------------------------
!
!   subroutine allocListenFile(a,varia,vanam,vacal)
!     use Mod_Listen
!     implicit none
!
!     class(MemoryMan)            :: a
!     type(ListenFile), allocatable           :: varia
!     character(*), intent(in)    :: vanam       ! Variable name
!     character(*), intent(in)    :: vacal       ! Calling routine
!
!     integer(ip)                 :: istat,lbyts
!
!     allocate(varia,stat=istat)
!     if(istat==0) then
!       lbyts=maxwp*(5+5*rp)+8*ip+151
!     else
!       call memerr(0,vanam,vacal,a%luerr)
!     end if
!     call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
!   end subroutine
!
!   subroutine deallocListenFile(a,varia,vanam,vacal)
!     use Mod_Listen
!     implicit none
!     class(MemoryMan)            :: a
!     type(ListenFile), allocatable           :: varia(:)
!     character(*), intent(in)    :: vanam       ! Variable name
!     character(*), intent(in)    :: vacal       ! Calling routine
!
!     integer(ip)                 :: istat,lbyts
!
!     deallocate(varia,stat=istat)
!     if(istat==0) then
!       lbyts=-maxwp*(5+5*rp)+8*ip+151
!     else
!       call memerr(1,vanam,vacal,a%luerr)
!     end if
!     call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
!   end subroutine

!---------------------------------------------------------------------------
! PostprFile
!---------------------------------------------------------------------------

!   subroutine allocPostprFile(a,varia,vanam,vacal)
!     use Mod_postpr
!     implicit none
!
!     class(MemoryMan)            :: a
!     type(PostprFile), pointer           :: varia
!     character(*), intent(in)    :: vanam       ! Variable name
!     character(*), intent(in)    :: vacal       ! Calling routine
!
!     integer(ip)                 :: istat,lbyts
!
!     allocate(varia,stat=istat)
!     if(istat==0) then
!       lbyts=6*ip
!     else
!       call memerr(0,vanam,vacal,a%luerr)
!     end if
!     call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
!   end subroutine
!
!   subroutine deallocPostprFile(a,varia,vanam,vacal)
!     use Mod_postpr
!     implicit none
!     class(MemoryMan)            :: a
!     type(PostprFile), allocatable           :: varia(:)
!     character(*), intent(in)    :: vanam       ! Variable name
!     character(*), intent(in)    :: vacal       ! Calling routine
!
!     integer(ip)                 :: istat,lbyts
!
!     deallocate(varia,stat=istat)
!     if(istat==0) then
!       lbyts=-6*ip
!     else
!       call memerr(1,vanam,vacal,a%luerr)
!     end if
!     call memctr(a%CurrentMemo,a%MaxMemo,lbyts,vanam,a%lumem,vacal)
!   end subroutine

!-----------------------------------------------------------------------
! Count and error functions
!-----------------------------------------------------------------------

   subroutine memerr(itask,vanam,vacal,luerr)
      implicit none
      integer(ip), intent(in) :: itask,luerr
      character*(*)           :: vanam,vacal
      if(itask==0) then
         if(luerr>0) then
            write(luerr,'(a)') trim(vacal)//': MEMORY FOR '//trim(vanam)//' CANNOT BE ALLOCATED'
         else
            write(6,'(a)') trim(vacal)//': MEMORY FOR '//trim(vanam)//' CANNOT BE ALLOCATED'
         end if
         call runend(trim(vacal)//': MEMORY FOR '//trim(vanam)//' CANNOT BE ALLOCATED')
         stop
      else
         if(luerr>0) then
            write(luerr,'(a)') trim(vacal)//': MEMORY FOR '//trim(vanam)//' CANNOT BE DEALLOCATED'
         else
            write(6,'(a)') trim(vacal)//': MEMORY FOR '//trim(vanam)//' CANNOT BE DEALLOCATED'
         end if
         call runend(trim(vacal)//': MEMORY FOR '//trim(vanam)//' CANNOT BE DEALLOCATED')
         stop
      end if
   end subroutine memerr

   subroutine GetLumem(a,lumem)
      implicit none
      class(MemoryMan) :: a
      integer(ip) :: lumem

      lumem = a%lumem
   end subroutine

end module Mod_memor

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
!##########################################################################
        program fdstag
!##########################################################################
        use mpi
        use vars
        implicit none

        call init_parallelisation

        call read_mdmap

        call read_control

        call read_infodom

        call alloc_dom

        call localparameters

        call initial

	  IF (L_LSM)  CALL initial_LSM_3D_channel

        call initflowfield

        IF (LROUGH)  THEN								!Richard2015
          IF (.not.LRESTART) THEN
          call init_rough
          ELSE
          call rough_restart
         END IF
        END IF

        if (LIMB) call imb_initial

        if (LIMB) call PartLocMPI							!Pablo2015

        call iniflux

        if(.not.LRESTART) then
           if (time_averaging) then
              call update_mean
	      if (noise.gt.0.0) call add_noise(noise)
           end if

          ! if(SGS) then
           !   if(sgs_model.eq.1) then
           !      call eddyv_smag
           !   else if(sgs_model.eq.2) then
           !      call eddyv_wale
           !   else if(sgs_model.eq.3) then
           !      call eddyv_1eqn
           !   else if(sgs_model.eq.4) then
           !      call eddyv_keps
           !   end if
          ! end if
        end if

	  if (LPT) then						!Brunho2013
			if (myrank.eq.0) open(unit=202, file='particulillas')
			call init_particle
	  endif

        if ((solver.eq.2).and.(.not.L_LSM)) call coeff

        call MPI_BARRIER (MPI_COMM_WORLD,ierr)
        if(myrank.eq.0) then
           write (numfile,*) '============START ITERATIONS========='
           write (6,*) '============START ITERATIONS========='
        end if

        call flosol

        call end_parallelisation

        end program
!##########################################################################

ext: for


text:
*     TB01ZD EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2020 NICONET e.V.
*
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX, PMAX
      PARAMETER        ( NMAX = 20, PMAX = 20 )
      INTEGER          LDA, LDC, LDZ
      PARAMETER        ( LDA = NMAX, LDC = PMAX, LDZ = NMAX )
      INTEGER          LDWORK
      PARAMETER        ( LDWORK = MAX( NMAX, PMAX ) )
*     .. Local Scalars ..
      DOUBLE PRECISION TOL
      INTEGER          I, INFO, J, N, NCONT, P
      CHARACTER*1      JOBZ
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), B(NMAX), C(LDC,NMAX), DWORK(LDWORK),
     $                 TAU(NMAX), Z(LDZ,NMAX)
*     .. External Functions ..
      LOGICAL          LSAME
      EXTERNAL         LSAME
*     .. External Subroutines ..
      EXTERNAL         TB01ZD, DORGQR
*     .. Intrinsic Functions ..
      INTRINSIC        MAX
*     .. Executable Statements ..
*
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read in the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * ) N, P, TOL, JOBZ
      IF ( N.LE.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99993 ) N
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
         READ ( NIN, FMT = * ) ( B(I), I = 1,N )
         IF ( P.LE.0 .OR. P.GT.PMAX ) THEN
            WRITE ( NOUT, FMT = 99992 ) P
         ELSE
            READ ( NIN, FMT = * ) ( ( C(I,J), J = 1,N ), I = 1,P )
*           Find a controllable realization for the given system.
            CALL TB01ZD( JOBZ, N, P, A, LDA, B, C, LDC, NCONT, Z, LDZ,
     $                   TAU, TOL, DWORK, LDWORK, INFO )
*
            IF ( INFO.NE.0 ) THEN
               WRITE ( NOUT, FMT = 99998 ) INFO
            ELSE
               WRITE ( NOUT, FMT = 99997 ) NCONT
               DO 20 I = 1, NCONT
                  WRITE ( NOUT, FMT = 99994 ) ( A(I,J), J = 1,NCONT )
   20          CONTINUE
               WRITE ( NOUT, FMT = 99996 ) ( B(I), I = 1,NCONT )
               WRITE ( NOUT, FMT = 99991 )
               DO 30 I = 1, P
                  WRITE ( NOUT, FMT = 99994 ) ( C(I,J), J = 1,NCONT )
   30          CONTINUE
               IF ( LSAME( JOBZ, 'F' ) )
     $            CALL DORGQR( N, N, N, Z, LDZ, TAU, DWORK, LDWORK,
     $                         INFO )
               IF ( LSAME( JOBZ, 'F' ).OR.LSAME( JOBZ, 'I' ) ) THEN
                  WRITE ( NOUT, FMT = 99995 )
                  DO 40 I = 1, N
                     WRITE ( NOUT, FMT = 99994 ) ( Z(I,J), J = 1,N )
   40             CONTINUE
               END IF
            END IF
         END IF
      END IF
      STOP
*
99999 FORMAT (' TB01ZD EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from TB01ZD = ',I2)
99997 FORMAT (' The order of the controllable state-space representati',
     $       'on = ',I2,//' The state dynamics matrix A of a controlla',
     $       'ble realization is ')
99996 FORMAT (/' The input/state vector B of a controllable realizatio',
     $       'n is ',/(1X,F8.4))
99995 FORMAT (/' The similarity transformation matrix Z is ')
99994 FORMAT (20(1X,F8.4))
99993 FORMAT (/' N is out of range.',/' N = ',I5)
99992 FORMAT (/' P is out of range.',/' P = ',I5)
99991 FORMAT (/' The output/state matrix C of a controllable realizati',
     $       'on is ')
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference