text:
! Direct Inversion for 3-D Azimuthal Anisotropy
! V1.1 (2017) based on Mineos to calculate kernel
! V2.0 (2019) based on CPS-tregn96 to calcualte kernel
! Copyright:
!    Author: Chuanming Liu (USTC)
!     Email: Chuanming.liu@colorado.edu
! Reference:
! Liu, C., Yao, H., Yang, H., Shen, W., Fang, H., Hu, S., Qiao, L., 2019. Direct
! inversion for three-dimensional shear wavespeed azimuthal anisotropy based on
! surface-wave ray tracing: methodology and application to Yunnan, southwest China.
! Journal of Geophysics Research: Solid Earth. doi: 10.1029/2018JB016920

program SurfAniso
        use lsmrModule, only:lsmr
        use lsmrblasInterface, only : dnrm2
        use omp_lib
        implicit none

        character inputfile*100
        character logfile*100
        character outmodel*100
        character outsyn*100
        logical ex
        character dummy*40
        character datafile*80
        integer nx,ny,nz
        real goxd,gozd
        real dvxd,dvzd
        integer nsrc,nrc
        real weightVs,weight0,weight1,weightGcs,weight
        real damp,dampAA,dampVs
        real minthk
        integer kmax,kmaxRc
        real*8,dimension(:),allocatable:: tRc
        real,dimension(:),allocatable:: depz
        integer itn
        integer nout
        integer localSize
        real mean,std_devs,balances,balanceb
        integer msurf
        real,dimension(:),allocatable:: obst,dsyn,cbst,dist
        real,dimension(:),allocatable:: pvall
      	real sta1_lat,sta1_lon,sta2_lat,sta2_lon
      	real dist1
        integer dall
       	integer istep
       	real,parameter :: pi=3.1415926535898
       	integer checkstat
       	integer ii,jj,kk
        real, dimension (:,:), allocatable :: scxf,sczf
        real, dimension (:,:,:), allocatable :: rcxf,rczf
        integer,dimension(:,:),allocatable::wavetype,igrt,nrc1
       	integer,dimension(:),allocatable::nsrc1
       	integer,dimension(:,:),allocatable::periods
       	real,dimension(:),allocatable::rw
       	integer, dimension(:),allocatable::iw,col
        real,dimension(:),allocatable::dv
        real,dimension(:,:,:),allocatable::vsf
      	character strf
      	integer veltp,wavetp
      	real velvalue
      	integer knum,knumo,err
      	integer istep1,istep2
      	integer period
      	integer knumi,srcnum,count1
      	integer HorizonType,VerticalType
      	character line*200
        integer iter,maxiter
        ! integer*4/8
        integer maxnar
        real acond,anorm,arnorm,rnorm,xnorm
        character str1
        real atol,btol
        real conlim
        integer istop
        integer itnlim
        integer lenrw,leniw
        integer nar,nar_tmp,nars
        integer count3,nvz,nvx
        integer m,maxvp,n
        integer i,j,k
        real spfra
        real noiselevel
        integer ifsyn
        integer writepath
        real averdws
        real maxnorm

        character SenGTruefile*200, SenGfile*200
        integer lay
        integer maxm
        real, dimension(:), allocatable:: sigmaT,resSigma,datweight
        real, dimension(:), allocatable:: Tdata,resbst,fwdT,fwdTvs,fwdTaa,Tref,RefTaa,resbst_iso

        real meandeltaT,stddeltaT
        real,dimension(:,:),allocatable:: GVs,GGs,GGc
        real,dimension(:,:,:), allocatable :: Lsen_Gsc
        real gaussian
        external gaussian
        real synT
        !real,dimension(:,:,:),allocatable::vsRela,vsAbs,vsreal
        real,dimension(:,:,:),allocatable:: gcf,gsf
        integer nar1,nar2
        real vsref
        real res2Norm,Mnorm2,Enorm2,Vsnorm2,Gcnorm2,Gsnorm2
        real resNSigma2Norm
        real Enorm2Lame,Vsnorm2Lame,Gcnorm2Lame,Gsnorm2Lame
        real MwNorm2

        real*8,dimension(:,:),allocatable:: tRcV
        real sumObs
        integer Refwritepath,realwritepath
        real VariGc,VariGs,VariVs
        integer count4,para
        real LambdaVs,LambdaGc,LambdaGs
        real PreRes,meanAbs
        integer iter_mod,rmax
        character(len=2) :: id='00'
        character(len=60) :: filename
        real*8 :: startT,endT
        logical:: isTikh
        real :: summ, threshold0,thresholdVs,threshold
        real :: pertV
        logical :: iso_inv, iso_mod
        real Minvel,MaxVel
        real,dimension(:),allocatable:: norm
        real,dimension(:),allocatable:: GcInv, GsInv, VsInv
        real mindVs, maxdVs, meadVs
        real minGs, maxGs, meaGs
        real minGc, maxGc, meaGc
        real lameGcs, lameVs
        real mean2
        integer Nre, NreVs, narVs
        ! open output for LSMR
        startT=OMP_get_wtime()
        nout=36
        open(nout,file='lsmr.txt')
        ! terminal output
        write(*,*)
        write(*,*) '                       DAzimSurfTomo'
        write(*,*)

        ! read contral file
        if (iargc() < 1) then
            write(*,*) 'input file [para.in (Default)]:'
            read(*,'(a)') inputfile
            if (len_trim(inputfile) <=1 ) then
                inputfile='para.in'
            else
                inputfile=inputfile(1:len_trim(inputfile))
            endif
        else
            call getarg(1, inputfile)
        endif
        inquire(file=inputfile, exist=ex)
        if (.not. ex)   stop 'unable to open the inputfile'

        open(10,file=inputfile,status='old',action='read')
        read(10,'(a30)') dummy
        read(10,'(a30)') dummy
        read(10,'(a30)') dummy
        ! read(10,*) SenGfile
        read(10,*) datafile
        read(10,*) nx,ny,nz
        read(10,*) goxd,gozd
        read(10,*) dvxd,dvzd
        read(10,*) minthk
        read(10,*) Minvel,Maxvel
        read(10,*) nsrc
        read(10,*) spfra
        read(10,*) maxiter
        read(10,*) iso_mod
        read(10,'(a30)') dummy
        read(10,*) weightVs
        read(10,*) weightGcs
        read(10,*) damp
        ! read(10,*) isTikh
        ! read(10,*) LcorrXY,LcorrZZ
        ! write(*,*) 'integrated  sensitivity kernel for Vs and Gcs file based on MOD  '
        ! write(*,'(a)')SenGfile
        write(*,*)'input Rayleigh wave phase velocity data file:'
        write(*,'(a)') datafile
        write(*,*)  'model origin:latitude,longitue'
        write(*,'(2f10.4)') goxd,gozd
        write(*,*) 'grid spacing:latitude,longitue'
        write(*,'(2f10.4)') dvxd,dvzd
        write(*,*) 'model dimension:nx,ny,nz'
        write(*,'(3i5)') nx,ny,nz
        write(*,*)'depth refined interval layer '
        write(*,'(f8.1)')minthk
        write(*,*)'smoothing for dVsv '
        write(*,'(f8.1)')weightVs
        write(*,*)'smoothing for Gc, Gs '
        write(*,'(f8.1)')weightGcs
        write(*,*)'damping'
        write(*,'(f8.1)')damp

        ! write(*,*)' Regularization Type: (T) 1st order Tikhonov ;(F) Gaussian'
        ! write(*,*) 'correlation length: XY, ZZ (km)'
        ! write(*,'(f8.1)')LcorrXY, LcorrZZ

        if (nz.LE.1)  stop 'error nz value.'
        read(10,'(a30)')dummy
        read(10,*) kmaxRc
        write(*,*) 'number of period'
        write(*,'(i6)') kmaxRc

        if(kmaxRc.gt.0)then
           allocate(tRc(kmaxRc),stat=checkstat)
           if (checkstat > 0) stop 'error allocating RP'
           read(10,*)(tRc(i),i=1,kmaxRc)
        else
            stop 'Can only deal with Rayleigh wave phase velocity data!'
        endif

	    write(logfile,'(a,a)') trim(inputfile),'_inv.log'
        !open(66,file=logfile,action='write')
        open(66, file=logfile)
        write(66,*)
        write(66,*) '                  DAzimSurfTomo'
        write(66,*)
        write(66,*) 'model origin:latitude,longitue'
        write(66,'(2f10.4)') goxd,gozd
        write(66,*) 'grid spacing:latitude,longitue'
        write(66,'(2f10.4)') dvxd,dvzd
        write(66,*) 'model dimension:nx,ny,nz'
        write(66,'(3i5)') nx,ny,nz

        write(*,*)'Rayleigh wave phase velocity used,periods:(s)'
        write(*,'(50f6.1)')(tRc(i),i=1,kmaxRc)
        write(66,*)'Rayleigh wave phase velocity used,periods:(s)'
        write(66,'(50f6.1)')(tRc(i),i=1,kmaxRc)

        ! close(10)
        nrc=nsrc
        kmax=kmaxRc
        lay=nz-1
        !-----------------------------------------------------------------------!
        ! read measurements
        inquire(file=datafile, exist=ex)
        if (.not. ex) then
            write(66,'(a)') 'unable to open the datafile'
            close(66)
            stop 'unable to open the datafile'
        endif
        write(*,*) 'begin load data file.....'

        open(unit=87,file=datafile,status='old')

        allocate(scxf(nsrc,kmax),sczf(nsrc,kmax),&
        rcxf(nrc,nsrc,kmax),rczf(nrc,nsrc,kmax),stat=checkstat)
        if(checkstat > 0)  write(6,*)'scxf error with allocate','nsrc=',nsrc,'kmax=',kmax

        allocate(periods(nsrc,kmax),wavetype(nsrc,kmax),nrc1(nsrc,kmax),nsrc1(kmax),&
        igrt(nsrc,kmax),stat=checkstat)
        if(checkstat > 0) write(6,*)'error with allocate'

        allocate(obst(nrc*nsrc*kmax),dist(nrc*nsrc*kmax),stat=checkstat)
        obst=0
        dist=0
        if(checkstat > 0)  write(6,*)'error with allocate'

        allocate(pvall(nrc*nsrc*kmax),stat=checkstat)
        if(checkstat > 0)    write(6,*)'error with allocate'

        istep=0
        istep2=0
        dall=0
        knumo=12345
	    knum=0
	    istep1=0
	    pvall=0
        do
            read(87,'(a)',iostat=err) line
            if(err.eq.0) then
                if(line(1:1).eq.'#') then
                   read(line,*) str1, sta1_lat, sta1_lon, period, wavetp, veltp
                   if(wavetp.eq.2.and.veltp.eq.0) knum=period
                   if(wavetp.eq.2.and.veltp.eq.1) stop 'can not deal with Rayleigh wave group data'
                   if(wavetp.eq.1.and.veltp.eq.0) stop 'can not deal with Love wave phase data'
                   if(wavetp.eq.1.and.veltp.eq.1) stop 'can not deal with Love wave group data'
                   if(knum.ne.knumo) then
                       istep=0
                       istep2=istep2+1
                   endif
                   istep=istep+1
                   istep1=0
                   sta1_lat=(90.0-sta1_lat)*pi/180.0
                   sta1_lon=sta1_lon*pi/180.0
                   scxf(istep,knum)=sta1_lat
                   sczf(istep,knum)=sta1_lon
                   periods(istep,knum)=period
                   wavetype(istep,knum)=wavetp
                   igrt(istep,knum)=veltp
                   nsrc1(knum)=istep
                   knumo=knum
                else
                   read(line,*) sta2_lat,sta2_lon,velvalue
                   istep1=istep1+1
                   dall=dall+1
                   sta2_lat=(90.0-sta2_lat)*pi/180.0
                   sta2_lon=sta2_lon*pi/180.0
                   rcxf(istep1,istep,knum)=sta2_lat
                   rczf(istep1,istep,knum)=sta2_lon
                   call delsph(sta1_lat,sta1_lon,sta2_lat,sta2_lon,dist1)
                   dist(dall)=dist1
                   obst(dall)=dist1/velvalue
                   pvall(dall)=velvalue
                   nrc1(istep,knum)=istep1
                endif
            else
                exit
            endif
        enddo
        close(87)
        write(*,'(a,i7)') ' Number of all measurements', dall
        !----------------------------------------------------------------------!
        !  Initialization
        allocate(depz(nz), stat=checkstat)
        allocate(vsf(nx,ny,nz), stat=checkstat)

        maxvp=(nx-2)*(ny-2)*(nz-1)
        maxm=(nx-2)*(ny-2)*(nz-1)*3
        maxnar = spfra*dall*nx*ny*nz*3 !sparsity fraction
        allocate(dv(maxm), stat=checkstat)
        allocate(GcInv(maxvp), GsInv(maxvp), VsInv(maxvp), stat=checkstat)
        allocate(rw(maxnar), stat=checkstat)
        allocate(iw(2*maxnar+1), stat=checkstat)
        allocate(col(maxnar), stat=checkstat)
        allocate(norm(maxvp), stat=checkstat)
        !allocate(cbst(dall+maxm*maxm),dsyn(dall),stat=checkstat)
        allocate(cbst(dall+maxm*3), dsyn(dall), stat=checkstat)
        allocate(dsyn(dall), stat=checkstat)
        allocate(sigmaT(dall), Tref(dall), datweight(dall), stat=checkstat)
        allocate(Tdata(dall), resbst(dall), fwdT(dall), fwdTvs(dall), fwdTaa(dall), RefTaa(dall), stat=checkstat)
        allocate(resbst_iso(dall), stat=checkstat)
        allocate(resSigma(dall), stat=checkstat)
        allocate(Lsen_Gsc(nx*ny,kmaxRc,nz-1), stat=checkstat)
        allocate(GVs(dall,maxvp), GGc(dall,maxvp), GGs(dall,maxvp), stat=checkstat)
        allocate(gcf(nx-2,ny-2,nz-1),gsf(nx-2,ny-2,nz-1), stat=checkstat)
        allocate( tRcV((nx-2)*(ny-2),kmaxRc), stat=checkstat)
        !----------------------------------------------------------------------!
        !  read reference isotropic model: vsf
        open(11, file='MOD', status='old')
        vsf=0
        read(11,*) (depz(i),i=1,nz)
        do k = 1,nz
            do j = 1,ny
                read(11,*)(vsf(i,j,k),i=1,nx)
            enddo
        enddo
        close(11)
        write(*,*) ' grid points in depth direction:(km)'
        write(*,'(50f7.2)') depz

        call  CalRmax(nz,depz,minthk,rmax)
        !----------------------------------------------------------------------!
        !                     iteration part
        !----------------------------------------------------------------------!
        RefTaa=0
        ! iso_mod=.false.
        open(34, file='IterVel.out')
        do iter=1, maxiter
        iter_mod=mod(iter, 2)
        write(6,*)  ' -----------------------------------------------------------'
        write(66,*) ' -----------------------------------------------------------'

        if (iso_mod) then
           iso_inv=.true.
           write(66,*)iter,'th iteration, invert for isotropic Vs para.'
           write(6,*) iter,'th iteration, invert for isotropic Vs para.'
           maxm =  maxvp
        else
            iso_inv=.false.
            write(66,*)iter,'th iteration, invert for dVs, Gc, Gs '
            write(6,*) iter,'th iteration, invert for dVs, Gc, Gs '
            ! maxm =  maxvp*2 ! joint
            maxm =  maxvp*3
        endif
        write(6,*)  ' -----------------------------------------------------------'
        write(66,*) ' -----------------------------------------------------------'
        !----------------------------------------------------------------------!
        ! compute G matrix based on the sensitivity kernel and ray-tracing
        ! forward calculation of traveltime misfit based on the reference iso- or aniso- model.
        Refwritepath = 0 ! switch of write ray path
        dsyn = 0
        GGc = 0
        GGs = 0
        GVs = 0
        tRcV = 0
        iw = 0
        rw = 0.0
        col = 0
        nar = 0
        Lsen_Gsc=0.0
        if (iso_inv) then
            call CalSurfG(nx,ny,nz,maxvp,vsf,iw,rw,col,dsyn,&
            GVs, dall,&
            goxd, gozd, dvxd, dvzd, kmaxRc, tRc, periods, depz, minthk,&
            scxf, sczf, rcxf, rczf, nrc1, nsrc1, kmax, nsrc, nrc, nar)
        else
            call CalSurfGAnisoJoint(nx, ny, nz, maxvp, vsf, iw, rw, col, dsyn, &
            GVs, GGc, GGs, Lsen_Gsc, dall, rmax, tRcV, &
            goxd, gozd, dvxd, dvzd, kmaxRc, tRc, periods,depz, minthk,&
            scxf, sczf, rcxf, rczf, nrc1, nsrc1, kmax, nsrc, nrc, nar, Refwritepath)
        endif
        ! write(*,*) ' Finish G matrix calculation.'
        !----------------------------------------------------------------------!
        ! output the corresponding isotropic phase velocity map of the MOD.
        if (iter.eq.1)then
            open(77,file='period_phaseVMOD.dat')
            call WTPeriodPhaseV(nx, ny, gozd, goxd, dvzd, dvxd, kmaxRc, tRc, tRcV, 77)
        end if
        !----------------------------------------------------------------------!
        ! calculate data residual based on the reference iso- or aniso- model.
        ! cbst: travel time residual.
        ! dsyn: predicated isotropic travel time
        ! Tdata: traveltime residual, used in the inversion
        cbst=0
        Tdata=0
        sigmaT=0
        Tref=0
        ! For both joint and iso inversion. Data fitting: obs(i) - Tref-iso(i)
        do i=1,dall
            Tref(i)=dsyn(i)
            cbst(i)=obst(i) - Tref(i)
            Tdata(i)=cbst(i)
        enddo

        ! Statistic of the data residual based on the reference model.
        mean=sum(cbst(1:dall))/dall
        std_devs=sqrt(sum((cbst(1:dall)-mean)**2)/dall)
        meanAbs=sum(abs(cbst(1:dall)))/dall

        write(6 ,'(a, f12.4,a,f10.2,a,f10.2,a)') '  Before Inversion: abs mean, std, RMS of Res:', meanAbs,' s ',&
        std_devs,' s ', dnrm2(dall,cbst,1)/sqrt(real(dall)),' s'
        write(66,'(a, f12.4,a,f10.2,a,f10.2,a)') '  Before Inversion: abs mean, std, RMS of Res:', meanAbs,' s ',&
        std_devs,' s ', dnrm2(dall,cbst,1)/sqrt(real(dall)),' s'

        !----------------------------------------------------------------------!
        ! data weight
        ! datweight=0.0
        ! if (iso_inv) then
        !     thresholdVs=0.2
        !     do i = 1,dall
        !     datweight(i) = 1.0
        !     if(abs(cbst(i)) > thresholdVs)  datweight(i) = exp(-(abs(cbst(i))-thresholdVs))
        !     cbst(i) = cbst(i)*datweight(i)
        !     enddo
        ! else
        !     ! thresholdVs=0.5
        !     threshold=std_devs*2.
        !     do i = 1,dall
        !     datweight(i) = 1.0
        !     if(abs(cbst(i)) > threshold)  datweight(i) = exp(-(abs(cbst(i))-threshold))
        !     cbst(i) = cbst(i)*datweight(i)
        !     enddo
        ! For Yunnan.
        call CalDdatSigma(dall, obst, cbst, sigmaT, meandeltaT)
	    do i=1,dall
            datweight(i)=1/sigmaT(i)
            cbst(i)=cbst(i)*datweight(i)
        enddo

        do i = 1,nar
            rw(i) = rw(i)*datweight(iw(1+i))
        enddo

        meanAbs=sum(abs(cbst(1:dall)))/dall
        write(6 ,'(a, f8.3, a, f8.3,a, f7.3, a)') '  mean data weight:',sum(datweight(1:dall))/dall, &
        ' |  abs data mean with weight:',meanAbs,'s  |  dt/t0:',  meandeltaT*100, ' %'
        write(66,'(a, f8.3, a, f8.3,a, f7.3, a)') '  mean data weight:',sum(datweight(1:dall))/dall, &
        ' |  abs data mean with weight:',meanAbs,'s  |  dt/t0:',  meandeltaT*100, ' %'
        if (iso_mod) then
            norm=0
            do i=1,nar
            norm(col(i))=norm(col(i))+abs(rw(i))
            enddo
        endif
        ! write out res for 1-th and final inversion
        ! if ((iter .eq. 1) .or. (iter .eq. maxiter)) then
        !     write(id,'(I2.2)') iter
        !     filename = 'Traveltime_use_'//TRIM(id)//'th.dat'
        !     open(88, file=filename)
        !     write(88,'(a)') 'Dist(km)        T_obs(s)       T_forward(s)         Res(s)        weight        W_Res'
        !     do i=1,dall
        !         write(88,*) dist(i), obst(i), dsyn(i), Tdata(i), datweight(i), cbst(i)
        !     enddo
        !     close(88)
        ! endif
        !----------------------------------------------------------------------!
        ! Setting of iteration.
        !----------------------------------------------------------------------!
        ! A*x = b
        ! m       input      m, the number of rows in A. (data)
        ! n       input      n, the number of columns in A. (model)
        ! count3: number of increase of row number.
        ! nar: number of G which value is not zero, including augmented matrix.
        !      if Gn*m does not contain zeros,  nar =n*m
        ! nar1: nar before regularization
        ! rw(i) gives value of the i-th non-zero value of G
        ! iw(i+1) means this value is in the iw(i+1)-th row.
        ! col(i) means this value is in the col(i)-th column. = iw(i+1+nar)
        !----------------------------------------------------------------------!
        nar1=nar
        lameGcs=weightGcs
        ! lameVs=dnrm2(dall,cbst,1)**2/dall*weightVs
        lameVs=weightVs
        count3=0
        if (iso_inv) then
            call TikhonovRegularization(nx, ny, nz, maxvp, dall, nar, rw, iw, col, count3, iso_inv, lameGcs, lameVs)
        else
            call        TikhRegul_joint(nx, ny, nz, maxvp, dall, nar, rw, iw, col, narVs, count3, lameGcs, lameVs)
        endif

        forall(i=1:count3)
           cbst(dall+i)=0
        end forall
        !----------------------------------------------------------------------!
        !----------------------------------------------------------------------!
        if (nar > maxnar) stop 'increase sparsity fraction(spfra)'
        write(*  ,'(a,3f8.2)') '  damp,  lamebda Gsc, lamebda Vs: ', damp, lameGcs, lameVs
        write(66 ,'(a,3f8.2)') '  damp,  lamebda Gsc, lamebda Vs: ', damp, lameGcs, lameVs

        m = dall + count3 ! data number, rows number of G
        n = maxm          ! model number, columns of G
        iw(1) = nar
        do i=1,nar
           iw(1+nar+i)=col(i)
        enddo
        leniw = 2*nar+1
        lenrw = nar
        dv = 0
        ! Control of LSQR
        ! atol:   estimate of error to G, accurate to about 4 digits, set btol = 1.0e-4. (default AA: 1e-4)
        ! btol:   estimate of error to data, accurate to about 4 digits, set btol = 1.0e-4. (default AA: 1e-4)
        ! conlim: the apparent condition number of the matrix Abar, controls the amp of Gc,s
        !         conlim and damp may be used separately or together to regularize ill-conditioned systems.
        ! itnlim: an upper limit on the number of iterations.
        if (iso_inv) then
            atol = 1e-3
            btol = 1e-3
            conlim = 1200
            itnlim = 1000
            localSize = n/4
        else
            atol = 1e-5
            btol = 1e-4
            conlim = 200
            itnlim = 500
            localSize = 10
        endif

        istop = 0
        anorm = 0.0
        acond = 0.0
        arnorm = 0.0
        xnorm = 0.0

        call LSMR(m, n, leniw, lenrw, iw, rw, cbst, damp,&
        atol, btol, conlim, itnlim, localSize, nout,&
        dv, istop, itn, anorm, acond, rnorm, arnorm, xnorm)
        if(istop==3) THEN
            write(* , '(a)') '  istop = 3, large condition number, LSMR failed'
            write(66, '(a)') '  istop = 3, large condition number, LSMR failed'
        endif
        write(*,'(a)') '  Finish LSMR.......'
        write(*, '(a, i7)')  '  itn=               ',itn
        write(66,'(a, i7)')  '  itn=               ',itn
        write(*,'(a, f7.1)') '  L2 norm of A=      ',anorm
        write(*, '(a, f7.1)')'  Condition NO. of A=',acond
        write(66,'(a, f7.1)')'  Condition NO. of A=',acond
        write(*,'(a, f7.1)') '  rnorm=             ',rnorm
        write(*,'(a, f7.1)') '  arnorm=            ',arnorm
        write(*,'(a, f7.3)') '  norm of dv =       ',xnorm
        !----------------------------------------------------------------------!
        ! Based on perturbation, construct isotropic or anisotropic velocity model
        ! For vsRela will be used in subroutine: CalSurfGAniso, it will be important to get a right value.
        ! vsf initially comes from MOD, then can be updated after each iteration.
        if (iso_inv) then
            do k=1,nz-1
            do j=1,ny-2
            do i=1,nx-2
            pertV = dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)
            if (pertV.ge.0.500) pertV=0.500
            if (pertV.le.-0.500) pertV=-0.500
            if (abs(pertV)<1e-5) pertV=0.0
            dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i) = pertV
            vsf(i+1,j+1,k)=vsf(i+1,j+1,k)+dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)

            if(vsf(i+1,j+1,k).lt.Minvel) vsf(i+1,j+1,k)=Minvel
            if(vsf(i+1,j+1,k).gt.Maxvel) vsf(i+1,j+1,k)=Maxvel
            enddo
            enddo
            enddo
        else
            gcf=0
            gsf=0
            do k=1,nz-1
            do j=1,ny-2
            do i=1,nx-2

            pertV = dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)
            if (pertV.ge. 0.5)   pertV=0.5
            if (pertV.le.-0.5)   pertV=-0.5
            if (abs(pertV)<1e-5) pertV=0.0
            dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i) = pertV
            vsf(i+1,j+1,k)=vsf(i+1,j+1,k)+dv((k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)

            if(vsf(i+1,j+1,k).lt.Minvel) vsf(i+1,j+1,k)=Minvel
            if(vsf(i+1,j+1,k).gt.Maxvel) vsf(i+1,j+1,k)=Maxvel

            gcf(i,j,k)=dv(maxvp  +(k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)
            gsf(i,j,k)=dv(maxvp*2+(k-1)*(nx-2)*(ny-2)+(j-1)*(nx-2)+i)
            enddo
            enddo
            enddo
        endif
        !----------------------------------------------------------------------!
        ! result analysis
        if (iso_inv)then
            mindVs=minval(dv(1:maxvp))
            maxdVs=maxval(dv(1:maxvp))
            meadVs=sum(abs(dv(1:maxvp)))/maxvp
            write(6 ,'(a,3f10.4)')  '  min  max and abs mean  dVs (km/s)', mindVs, maxdVs, meadVs
            write(66,'(a,3f10.4)')  '  min  max and abs mean  dVs (km/s)', mindVs, maxdVs, meadVs
            do k=1,nz-1
            VariVs=sum(abs(dv((k-1)*(nx-2)*(ny-2)+1:k*(nx-2)*(ny-2) )))/((nx-2)*(ny-2))
            write(66,'(a,f5.1,a,f5.1,a,f10.4)') '  Z ',depz(k),' - ',depz(k+1),&
            ' km  abs mean dVs (km/s)',VariVs
            write(6,'(a,f5.1,a,f5.1,a,f10.4)')  '  Z ',depz(k),' - ',depz(k+1),&
            ' km  abs mean dVs (km/s)',VariVs
            enddo
        else
            mindVs=minval(dv(1:maxvp))
            maxdVs=maxval(dv(1:maxvp))
            meadVs=sum(abs(dv(1:maxvp)))/maxvp
            write(6 ,'(a,3f10.4)')  '  min  max and abs mean  dVs (km/s)', mindVs, maxdVs, meadVs
            write(66,'(a,3f10.4)')  '  min  max and abs mean  dVs (km/s)', mindVs, maxdVs, meadVs
            minGc=minval(dv(maxvp+1:maxvp+maxvp))*100
            maxGc=maxval(dv(maxvp+1:maxvp+maxvp))*100
            meaGc=sum(abs(dv(maxvp+1:maxvp+maxvp)))/maxvp*100
            minGs=minval(dv(maxvp*2+1:maxvp+maxvp*2))*100
            maxGs=maxval(dv(maxvp*2+1:maxvp+maxvp*2))*100
            meaGs=sum(abs(dv(maxvp*2+1:maxvp+maxvp*2)))/maxvp*100
            write(6 ,'(a,3f10.4)')  '  min  max and abs mean   Gc/L (%) ', minGc, maxGc, meaGc
            write(66,'(a,3f10.4)')  '  min  max and abs mean   Gc/L (%) ', minGc, maxGc, meaGc
            write(6 ,'(a,3f10.4)')  '  min  max and abs mean   Gs/L (%) ', minGs, maxGs, meaGs
            write(66,'(a,3f10.4)')  '  min  max and abs mean   Gs/L (%) ', minGs, maxGs, meaGs
            do k=1, nz-1
            VariVs=sum(abs(dv((k-1)*(nx-2)*(ny-2)+1:k*(nx-2)*(ny-2) )))/((nx-2)*(ny-2))
            VariGc=sum(abs(gcf(1:nx-2,1:ny-2,k)))/((nx-2)*(ny-2))
            VariGs=sum(abs(gsf(1:nx-2,1:ny-2,k)))/((nx-2)*(ny-2))
            write(66,'(a, f5.1, a, f5.1, a, 2f10.3, f9.4)')'  Z ',depz(k),' - ',depz(k+1),' km  Abs Mean &
               Gc (%)  Gs (%)   dVs (km/s)', VariGc*100, VariGs*100, VariVs
            write(6 ,'(a, f5.1, a, f5.1, a, 2f10.3, f9.4)')'  Z ',depz(k),' - ',depz(k+1),' km  Abs Mean &
               Gc (%)  Gs (%)   dVs (km/s)', VariGc*100, VariGs*100, VariVs
            enddo
        endif
        !----------------------------------------------------------------------!
        ! Calculate ||Lm||2
        Mnorm2 = 0
        MwNorm2 = 0
        if (iso_inv) then
            Nre=nar-nar1
            call Calmodel2Norm(nar1, nar, maxvp, Nre, rw, col, dv, lameVs)
        else
            Nre=nar-nar1
            NreVs=narVs-nar1
            call Calmodel2NormJoint(nar1, nar, maxvp, NreVs, Nre, rw, col, dv, lameGcs, lameVs)
        endif
        !----------------------------------------------------------------------!
        ! Calculate ||W(Gm-d)||2
        ! fist forward calculate the traveltime.
        ! Tdata: inversion used data: tobs- tref(iso)
        ! Calculate deltaT=deltaTvs+deltaTaa
        ! Here restT=deltaT_in-deltaT_out
        resbst=0
        fwdTvs=0
        fwdTaa=0
        if (iso_inv) then
            call CalVsReslNorm(maxvp,dall, GVs, dv(1:maxvp), datweight, Tdata, fwdTvs, resbst)
        else
            VsInv=dv(1:maxvp)
            GcInv=dv(1+maxvp:2*maxvp)
            GsInv=dv(1+maxvp*2:3*maxvp)
            call  CalReslNormJoint(maxvp, dall, GVs, GGc, GGs, VsInv, GcInv, GsInv, datweight, Tdata, fwdTvs, fwdTaa, resbst)
            ! update T_AA. useless.
            RefTaa=0
            forall(i=1:dall)
                RefTaa(i)=fwdTaa(i)
            end forall
        endif
        !----------------------------------------------------------------------!
        !  Residual Traveltime Analsysi
        !  dsyn: ref-iso traveltime from G-matrix based on vsf model.
        !  Tdata: inversion used data: tobs- tref(iso)
        !  resbst: residual, without weight(res2Norm with weight)
        if ((iter .eq. 1) .or. (iter .eq. maxiter)) then
            filename = 'Traveltime_statis_'//TRIM(id)//'th.dat'
            open(88,file=TRIM(filename))
            if (iso_mod) then
                write(88,'(7a)')'   Dist(km)   T_obs(s)  T_ref_iso   Res(in)   dT(dvs)   Res(out)'
                do i=1,dall
                    write(88,'(3f10.3, 3e12.3)') dist(i), obst(i), dsyn(i), Tdata(i), fwdTvs(i), resbst(i)
                enddo
            else
                write(88,'(7a)')'          Dist(km)       T_obs(s)        T_ref-iso        Res(in)   &
                dT(aa)        dT(dvs)        Res(out)'
                do i=1,dall
                    write(88,'(3f10.4, 4e12.3)') dist(i), obst(i), dsyn(i), Tdata(i), fwdTvs(i), fwdTaa(i), resbst(i)
                enddo
            endif
            close(88)
        endif
        !----------------------------------------------------------------------!
        ! Statistic for the  data residual.
        mean = sum(resbst(1:dall))/dall
        meanAbs=sum(abs(resbst(1:dall)))/dall
        std_devs = sqrt(sum((resbst(1:dall)-mean)**2)/dall)

        write(6 ,'(a,f12.4,a,f10.2,a,f10.2,a)')  '  After Inversion: abs mean, std, RMS of Res :', meanAbs,' s ',&
        std_devs,' s ', dnrm2(dall,resbst,1)/sqrt(real(dall)),' s'
        write(66,'(a,f12.4,a,f10.2,a,f10.2,a)')  '  After Inversion: abs mean, std, RMS of Res :', meanAbs,' s ',&
        std_devs,' s ', dnrm2(dall,resbst,1)/sqrt(real(dall)),' s'

        ! mean=sum(abs(RefTaa(1:dall)))/dall
        ! mean2=sum(abs(fwdTvs(1:dall)))/dall
        ! write(66,'(a,f12.4,a,f12.4)')'  ABS Mean Taa:',mean,'s    ABS Mean T(dv)', mean2
        !----------------------------------------------------------------------!
        if (iso_mod) then
        write(34,*)',OUTPUT S VELOCITY AT ITERATION',iter
        do k=1,nz
        do j=1,ny
        write(34,'(100f7.3)') (vsf(i,j,k),i=1,nx)
        enddo
        enddo
        write(34,*)',OUTPUT DWS AT ITERATION',iter
        do k=1,nz-1
        do j=2,ny-1
        write(34,'(100f10.3)') (norm((k-1)*(ny-2)*(nx-2)+(j-2)*(nx-2)+i-1),i=2,nx-1)
        enddo
        enddo
        endif
        !----------------------------------------------------------------------!
        write(66,'(a)') ' '
        write(6 ,'(a)') '  '
        enddo ! iteration

        ! if (iso_inv) then
        filename = 'MOD_Ref'
        open(11, file=filename)
        do k=1,nz
            write(11,'(f7.1)',advance='no') depz(k)
        enddo
        do k = 1,nz
        do j = 1,ny
        do i=1,nx
        if (i.eq.1) then
        write(11,'(/f8.4)',advance='no')vsf(i,j,k)
        else
        write(11,'(f8.4)',advance='no')vsf(i,j,k)
        endif
        enddo
        enddo
        enddo
        close(11)
        ! endif
        !----------------------------------------------------------------------!
        ! output inversion result.
        open(63,file='DSurfTomo.inv')
        open(73,file='Gc_Gs_model.inv')
        call writeVsmodel(nx,ny,nz,gozd,goxd,dvzd,dvxd,depz,vsf,63)
        call writeAzimuthal(nx,ny,nz,gozd,goxd,dvzd,dvxd,depz,gcf,gsf,vsf,73)
        close(63)
        close(73)
        !----------------------------------------------------------------------!
        open(77,file='phaseV_FWD.dat')
        call WTPeriodPhaseV(nx,ny,gozd,goxd,dvzd,dvxd,kmaxRc,tRc,tRcV,77)
        !----------------------------------------------------------------------!
        write(*,'(a)')'  Begin forward calculate period azimuthal A1, A2.'
        open(42,file='period_Azm_tomo.inv',status='replace',action='write')
        call FwdAzimuthalAniMap(nx,ny,nz,maxvp,&
            goxd,gozd,dvxd,dvzd,kmaxRc,tRc,&
           gcf,gsf,Lsen_Gsc,tRcV)
        !----------------------------------------------------------------------!
        ! output result model
        ! note: for the input Vs model is grid model (nz point)
        ! BUT, the inversion kernel and result is layer (nz-1 layer) model
        ! for Gc,Gs, the input and result are both layer models, so I set the lower boundary depth as the output
        ! for Vs, input is grid model, while the inversion result is layered model.
        !  I set the mid depth between the points as the depth index, and the output model uses the average Vs.
        write(66,*)'  -----------------------------------------------------------'
        write(*,*)  '  Program finishes successfully'
        write(66,*) '  Program finishes successfully'

        write(*,*)  '  Output inverted shear velocity model: Vs_model_Syn.rela  Vs_model_Syn.abs'
        write(66,*) '  Output inverted shear velocity model: Vs_model_Syn.rela  Vs_model_Syn.abs'

        endT=OMP_get_wtime()
        write(* ,'(a, f13.1, a)') '   All time cost= ',endT-startT,"s"
        write(66,'(a, f13.1, a)') '   All time cost= ',endT-startT,"s"

        close(10)   ! surfaniso.in
        close(nout) ! close lsmr.txt
        close(66)   ! close surf_tomo.log
        close(34)   ! IterVel.out
        deallocate(obst)
        deallocate(dsyn)
        deallocate(dist)
        deallocate(pvall)
        deallocate(depz)
        deallocate(norm)
        deallocate(scxf,sczf)
        deallocate(rcxf,rczf)
        deallocate(wavetype,igrt,nrc1)
        deallocate(nsrc1,periods)
        deallocate(rw)
        deallocate(iw,col)
        deallocate(cbst)
        deallocate(dv)
        deallocate(vsf)
        deallocate(sigmaT,Tdata,resbst,fwdT,resSigma,resbst_iso)
        deallocate(fwdTvs,fwdTaa)
        deallocate(RefTaa,Tref)
        deallocate(GVs,GGs,GGc,Lsen_Gsc)

        deallocate(gcf,gsf)
        if(kmaxRc.gt.0) then
            deallocate(tRc)
        endif
        deallocate(tRcV)
end program

subroutine writeVsmodel(nx, ny, nz, gozd, goxd, dvzd, dvxd, depz, vsf, Idout)
        integer nx,ny,nz
        real goxd,gozd
        real dvxd,dvzd
        real depz(nz)
        ! real vsAbs(nx-2,ny-2,nz-1), vsRela(nx-2,ny-2,nz-1)
        integer Idout
        real:: vsf(nx,ny,nz)
        integer:: k,j,i
        real :: vsref
        ! Because the (gozd,goxd) point is NE point of the inversion region.---> NE move 1 point.
        do k=1,nz
              do j=1,ny
                    do i=1,nx
                        write(Idout,'(5f8.4)') gozd+(j-2)*dvzd,goxd-(i-2)*dvxd,depz(k),vsf(i,j,k)
                    enddo
                enddo
            enddo
end subroutine writeVsmodel

subroutine writeAzimuthal(nx, ny, nz, gozd, goxd, dvzd, dvxd, depz, Gc, Gs, vsf, Idout)
        integer nx,ny,nz
        real goxd,gozd
        real dvxd,dvzd
        real depz(nz)
        real Gc(nx-2,ny-2,nz-1),Gs(nx-2,ny-2,nz-1)
        real:: vsf(nx,ny,nz)
        integer Idout
        integer:: k,j,i
        real cosTmp,sinTmp
        real vsref
        real*8 :: pi=3.1415926535898
        do k=1,nz-1
            do j=1,ny-2
                do i=1,nx-2
                    cosTmp=Gc(i,j,k)
                    sinTmp=Gs(i,j,k)
                    AzimAmp=0.5*sqrt(cosTmp**2+sinTmp**2)
                    AzimAng=atan2(sinTmp,cosTmp)/pi*180
                    if (AzimAng.LT.0.0) AzimAng=AzimAng+360
                    AzimAng=0.5*AzimAng
                    ! format: lon lat  depth(lower) Vs(mid-depth)  Angle  Amp(%/100) Gc/L(%)  Gs/L(%)
                    vsref=(vsf(i+1,j+1,k)+vsf(i+1,j+1,k+1))/2
                    write(Idout,'(8f10.4)') gozd+(j-1)*dvzd,goxd-(i-1)*dvxd,depz(k+1),vsref,AzimAng,AzimAmp,&
                           Gc(i,j,k)*100,Gs(i,j,k)*100
                enddo
            enddo
        enddo
end subroutine writeAzimuthal

subroutine WTPeriodPhaseV(nx,ny,gozd,goxd,dvzd,dvxd,kmaxRc,tRc,pvRc,Idout)
        implicit none
        integer nx,ny
        real goxd,gozd
        real dvxd,dvzd
        integer kmaxRc
        real*8 tRc(kmaxRc)
        real*8 pvRc((nx-2)*(ny-2),kmaxRc)
        integer Idout
        integer:: tt,k,jj,ii
        ! note:  pvRc((nx-2)*(ny-2),kmaxRc) ----different pvRc in
        do tt=1,kmaxRc
            do jj=1,ny-2
                do ii=1,nx-2
                    write(Idout,'(5f10.4)') gozd+(jj-1)*dvzd,goxd-(ii-1)*dvxd,tRc(tt),pvRc((jj-1)*(nx-2)+ii,tt)
                enddo
            enddo
        enddo
        close(Idout)
end subroutine WTPeriodPhaseV

subroutine CalRmax(nz, depz, minthk0, rmax)
        implicit none
        integer nz
        real depz(nz)
        real minthk0, minthk
        integer rmax
        integer k,i
        integer NL
        parameter (NL=200)
        integer nsublay(NL)
        real thk
        rmax=0
        do i=1, nz-1
            thk = depz(i+1)-depz(i)
            minthk = thk/minthk0
            nsublay(i) = int((thk+1.0e-4)/minthk) + 1
        enddo
        rmax=sum(nsublay(1:nz-1))
end subroutine CalRmax

max_stars_repo_name: core-man/DAzimSurfTomo

max_issues_repo_name: core-man/DAzimSurfTomo

max_forks_repo_name: core-man/DAzimSurfTomo


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
PROGRAM testBinaryTrees
#include "UnitTest.h"
USE UnitTest
USE IntrType
USE BinaryTrees

IMPLICIT NONE

TYPE(BinaryTreeType_Integer),POINTER :: t => NULL()
TYPE(BinaryTreeType_Index),POINTER :: t2 => NULL()
INTEGER(SIK),ALLOCATABLE :: emptyi(:)
LOGICAL(SBK) :: bool

CREATE_TEST('Binary Trees')

WRITE(*,*) '==================================================='
WRITE(*,*) 'TESTING BINARY TREES...'
WRITE(*,*) '==================================================='

WRITE(*,*) 'TESTING INTEGER BINARY TREE'

ALLOCATE(emptyi(1:0))
CALL CreateBinaryTree(t,emptyi)
CALL BurnBinaryTree(t)
CALL CreateBinaryTreeRmDup(t,(/4,8,4,6,7,2,1/))
ASSERT(t%val == 4, 'CreateBinaryTreeRmDup(t,...)')
WRITE(*,*) '  Passed: CALL CreateBinaryTreeRmDup(t,...)'
CALL BurnBinaryTree(t)
CALL CreateBinaryTree(t,(/8,3,4,6,7,2,1/))
ASSERT(t%val == 8, 'CreateBinaryTree(t,...)')
WRITE(*,*) '  Passed: CALL CreateBinaryTree(t,...)'
CALL CreateBinaryTree(t,(/8,3,4,6,7,2,1/))
bool = .NOT.(.NOT.t%match(2) .AND. t%match(10) .AND. .NOT.t%match(4))
ASSERT(bool, 't%match(...)')
WRITE(*,*) '  Passed: t%match(...)'

CALL BurnBinaryTree(t)
ASSERT(.NOT.ASSOCIATED(t), 'BurnBinaryTree(t)')
WRITE(*,*) '  Passed: CALL BurnBinaryTree(t)'
WRITE(*,*) '---------------------------------------------------'

WRITE(*,*) 'TESTING INDEXED BINARY TREE'
CALL CreateBinaryTreeRmDup(t2,(/4,8,4,6,7,2,1/),1,7)
ASSERT(t2%val == 4, 'CreateBinaryTreeRmDup(t,...)')
WRITE(*,*) '  Passed: CALL CreateBinaryTreeRmDup(t,...)'
CALL BurnBinaryTree(t2)
CALL CreateBinaryTree(t2,(/8,3,4,6,7,2,1/))
ASSERT(t2%val == 8, 'CreateBinaryTree(t2,...)')
WRITE(*,*) '  Passed: CALL CreateBinaryTree(t2,...)'
CALL CreateBinaryTree(t2,(/8,3,4,6,7,2,1/),1,7)
bool = .NOT.(.NOT.t2%match(2) .AND. t2%match(10) .AND. .NOT.t2%match(4))
ASSERT(bool, 't2%match(...)')
WRITE(*,*) '  Passed: t2%match(...)'
bool = .NOT.(t2%findIndex(4) /= 3 .AND. t2%findIndex(2) /= 6 .AND. &
      t2%findIndex(10) /= -1)
ASSERT(bool, 't2%findIndex(...)')
WRITE(*,*) '  Passed: t2%findIndex(...)'

CALL BurnBinaryTree(t2)
ASSERT(.NOT.ASSOCIATED(t), 'BurnBinaryTree(t2)')
WRITE(*,*) '  Passed: CALL BurnBinaryTree(t2)'

CALL SortedBinaryTree(t2,1,(/1,2,3,4,5,6/))
bool = .NOT.((t2%val /= 4 .AND. t2%index /= 4) .AND. &
      (t2%left%val /= 2 .AND. t2%left%index /= 2) .AND. &
      (t2%right%val /= 6 .AND. t2%right%index /= 6) .AND. &
      (t2%left%left%val /= 1 .AND. t2%left%left%index /= 1) .AND. &
      (t2%left%right%val /= 3 .AND. t2%left%right%index /= 3) .AND. &
      (t2%right%left%val /= 5 .AND. t2%right%left%index /= 5))
ASSERT(bool, 'SortedBinaryTree(...)')
WRITE(*,*) '  Passed: CALL SortedBinaryTree(...)'
CALL BurnBinaryTree(t2)
WRITE(*,*) '---------------------------------------------------'

DEALLOCATE(emptyi)

FINALIZE_TEST()

ENDPROGRAM testBinaryTrees

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: wgurecky/Futility


text:
subroutine dwnvel(nmmax     ,kmax      ,icx       , &
                & kcs       ,kfu       ,kfv       , &
                & kcu       ,kcv       ,s1        ,dps       ,u0eul     , &
                & v0eul     ,uuu       ,vvv       ,umod      ,zumod     , &
                & sig       ,hu        ,hv        ,kfsed     ,deltau    , &
                & deltav    ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: dwnvel.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/compute_sediment/dwnvel.f90 $
!!--description-----------------------------------------------------------------
! Calculate velocity components and magnitude at the
! zeta points based on velocity in the bottom
! computational layer
!
! Note: uses downwind velocity at any internal point,
! uses internal velocity at any open boundary, uses
! half of internal velocity in direction of any
! closed boundary or dry point.
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    !
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    real(fp)                             , pointer :: eps
    type (mornumericstype)               , pointer :: mornum
    logical                              , pointer :: v2dwbl
!
! Global variables
!
    integer                                           , intent(in)  :: icx
    integer                                           , intent(in)  :: kmax  !  Description and declaration in esm_alloc_int.f90
    integer                                           , intent(in)  :: nmmax !  Description and declaration in dimens.igs
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: kcs   !  Description and declaration in esm_alloc_int.f90
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: kcu   !  Description and declaration in esm_alloc_int.f90
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: kcv   !  Description and declaration in esm_alloc_int.f90
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: kfsed !  Description and declaration in esm_alloc_int.f90
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: kfu   !  Description and declaration in esm_alloc_int.f90
    integer   , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: kfv   !  Description and declaration in esm_alloc_int.f90
    real(prec), dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: dps   !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: deltau !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: deltav !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: hu    !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: hv    !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(in)  :: s1    !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub, kmax), intent(in)  :: u0eul
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub, kmax), intent(in)  :: v0eul
    real(fp)  , dimension(kmax)                       , intent(in)  :: sig   !  Description and declaration in esm_alloc_real.f90
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(out) :: umod
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(out) :: uuu
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(out) :: vvv
    real(fp)  , dimension(gdp%d%nmlb:gdp%d%nmub)      , intent(out) :: zumod
!
! Local variables
!
    integer  :: k
    integer  :: kmx
    integer  :: ndm
    integer  :: ndmd
    integer  :: ndmu
    integer  :: nm
    integer  :: nm_u1
    integer  :: nm_u2
    integer  :: nm_v1
    integer  :: nm_v2
    integer  :: nmd
    integer  :: nmu
    integer  :: num
    integer  :: numd
    real(fp) :: cc
    real(fp) :: fact
    real(fp) :: maxdepfrac
    real(fp) :: h1
    real(fp) :: ufac
    real(fp) :: uu
    real(fp) :: vfac
    real(fp) :: vv
    real(fp) :: deltas
!
!! executable statements -------------------------------------------------------
!
    eps                 => gdp%gdconst%eps
    mornum              => gdp%gdmorpar%mornum
    v2dwbl              => gdp%gdnumeco%v2dwbl
    !
    do nm = 1, nmmax
       !
       if (kfsed(nm) == 0) then
          uuu  (nm) = 0.0_fp
          vvv  (nm) = 0.0_fp
          umod (nm) = 0.0_fp
          zumod(nm) = 0.0_fp
          cycle
       endif
       !
       nmd  = nm  - icx
       numd = nmd + 1
       ndmd = nmd - 1
       num  = nm  + 1
       ndm  = nm  - 1
       nmu  = nm  + icx
       ndmu = nmu - 1
       !
       uu = 0.0_fp
       vv = 0.0_fp
       h1 = s1(nm) + real(dps(nm),fp)
       !
       if (v2dwbl) then
          fact   = max(kfu(nm) + kfu(nmd) + kfv(nm) + kfv(ndm), 1)
          deltas = (deltau(nm) + deltau(nmd) + deltav(nm) + deltav(ndm)) / fact
          maxdepfrac = 0.5_fp
       else
          deltas = 0.05_fp
          maxdepfrac = 0.05_fp
       endif
       !
       do k = kmax, 1, -1
          cc  = (1.0 + sig(k))*h1
          kmx = k
          if (cc>=maxdepfrac*h1 .or. cc>=deltas) then
             exit
          endif
       enddo
       !
       ufac = 0.5_fp
       vfac = 0.5_fp
       if (abs(kcs(nm)) == 1) then
          !
          ! Internal point
          ! Set velocity in U direction.
          !
          nm_u1 = nm
          nm_u2 = nmd
          !
          ! Set velocity in V direction.
          !
          nm_v1 = nm
          nm_v2 = ndm
       elseif (kcu(nm) + kcu(nmd) == 1) then
          !
          ! Open boundary (kcs(nm)==2) in v-direction
          !
          ! Set velocity in U direction.
          !
          nm_u1 = nm
          nm_u2 = nmd
          ufac  = 1.0_fp
          !
          ! Set velocity in V direction.
          !
          if (kcu(nm) == 1) then
             !
             ! Open boundary at left-hand side
             !
             nm_v1 = nmu
             nm_v2 = ndmu
          else
             !
             ! Open boundary at right-hand side
             !
             nm_v1 = nmd
             nm_v2 = ndmd
          endif
       else
          !
          ! Open boundary (kcs(nm)==2) in u-direction
          !
          ! Set velocity in U direction.
          !
          if (kcv(nm) == 1) then
             !
             ! Open boundary at lower side
             !
             nm_u1 = num
             nm_u2 = numd
          else
             !
             ! Open boundary at upper side
             !
             nm_u1 = ndm
             nm_u2 = ndmd
          endif
          !
          ! Set velocity in V direction.
          !
          nm_v1 = nm
          nm_v2 = ndm
          vfac  = 1.0_fp
       endif
       !
       uu = ufac * (  kfu(nm_u1)*u0eul(nm_u1, kmx)*hu(nm_u1) &
          &         + kfu(nm_u2)*u0eul(nm_u2, kmx)*hu(nm_u2)  )
       vv = vfac * (  kfv(nm_v1)*v0eul(nm_v1, kmx)*hv(nm_v1) &
          &         + kfv(nm_v2)*v0eul(nm_v2, kmx)*hv(nm_v2)  )
       !
       if (mornum%maximumwaterdepth) then
          !
          ! prevent any increase in velocity due to a decrease in water depth
          !
          if (kfu(nm_u1)==1) h1 = max(h1,hu(nm_u1))
          if (kfu(nm_u2)==1) h1 = max(h1,hu(nm_u2))
          if (kfv(nm_v1)==1) h1 = max(h1,hv(nm_v1))
          if (kfv(nm_v2)==1) h1 = max(h1,hv(nm_v2))
       endif
       !
       ! h1 will not be zero, because of cycle statement when kfsed==0
       !
       uu = uu/h1
       vv = vv/h1
       !
       if (abs(uu) < eps) then
          uu = 0.0_fp
       endif
       if (abs(vv) < eps) then
          vv = 0.0_fp
       endif
       !
       ! Calculate resultant velocity magnitude and height
       !
       umod(nm) = (uu*uu + vv*vv)**0.5
       if (kmax > 1) then
          zumod(nm) = h1 * (1.0+sig(kmx))
       else
          !
          ! (1+h1/z0rou)**(z0rou/h1)*exp(-1)*(z0rou+h1)-z0rou
          ! equals approximately exp(-1)*h1 if h1>>z0rou
          ! exp(-1) is approx 0.368
          !
          zumod(nm) = h1*0.368
       endif
       uuu(nm) = uu
       vvv(nm) = vv
    enddo
end subroutine dwnvel

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


!STARTOFREGISTRYGENERATEDFILE 'IfW_TSFFWind_Types.f90'
!
! WARNING This file is generated automatically by the FAST registry.
! Do not edit.  Your changes to this file will be lost.
!
! FAST Registry (v3.02.00, 22-Jun-2016)
!*********************************************************************************************************************************
! IfW_TSFFWind_Types
!.................................................................................................................................
! This file is part of IfW_TSFFWind.
!
! Copyright (C) 2012-2016 National Renewable Energy Laboratory
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
!
! W A R N I N G : This file was automatically generated from the FAST registry.  Changes made to this file may be lost.
!
!*********************************************************************************************************************************
!> This module contains the user-defined types needed in IfW_TSFFWind. It also contains copy, destroy, pack, and
!! unpack routines associated with each defined data type. This code is automatically generated by the FAST Registry.
MODULE IfW_TSFFWind_Types
!---------------------------------------------------------------------------------------------------------------------------------
USE NWTC_Library
IMPLICIT NONE
! =========  IfW_TSFFWind_InitInputType  =======
  TYPE, PUBLIC :: IfW_TSFFWind_InitInputType
    CHARACTER(1024)  :: WindFileName      !< Name of the wind file to use [-]
    INTEGER(IntKi)  :: SumFileUnit      !< Unit number for the summary file (-1 for none).  Provided by IfW. [-]
  END TYPE IfW_TSFFWind_InitInputType
! =======================
! =========  IfW_TSFFWind_InitOutputType  =======
  TYPE, PUBLIC :: IfW_TSFFWind_InitOutputType
    TYPE(ProgDesc)  :: Ver      !< Version information off FFWind submodule [-]
  END TYPE IfW_TSFFWind_InitOutputType
! =======================
! =========  IfW_TSFFWind_MiscVarType  =======
  TYPE, PUBLIC :: IfW_TSFFWind_MiscVarType
    INTEGER(IntKi)  :: TimeIndex = 0      !< An Index into the TData array [-]
  END TYPE IfW_TSFFWind_MiscVarType
! =======================
! =========  IfW_TSFFWind_ParameterType  =======
  TYPE, PUBLIC :: IfW_TSFFWind_ParameterType
    REAL(DbKi)  :: DT      !< Time step for cont. state integration & disc. state update [seconds]
    LOGICAL  :: TowerDataExist = .FALSE.      !< Data exists for the tower [-]
    LOGICAL  :: Periodic = .FALSE.      !< Flag to indicate if the wind file is periodic [-]
    REAL(SiKi) , DIMENSION(:,:,:,:), ALLOCATABLE  :: FFData      !< Array of FF data [-]
    REAL(SiKi) , DIMENSION(:,:,:), ALLOCATABLE  :: FFTower      !< Array of data along tower, below FF array [-]
    REAL(ReKi)  :: FFDTime = 0      !< Delta time [seconds]
    REAL(ReKi)  :: FFRate = 0      !< Data rate (1/FFTime) [Hertz]
    REAL(ReKi)  :: FFYHWid = 0      !< Half the grid width [meters]
    REAL(ReKi)  :: FFZHWid = 0      !< Half the grid height [meters]
    REAL(ReKi)  :: RefHt = 0      !< Reference (hub) height of the grid [meters]
    REAL(ReKi)  :: GridBase = 0      !< the height of the bottom of the grid [meters]
    REAL(ReKi)  :: InitXPosition = 0      !< the initial x position of grid (distance in FF is offset) [meters]
    REAL(ReKi)  :: InvFFYD = 0      !< reciprocal of delta y [1/meters]
    REAL(ReKi)  :: InvFFZD = 0      !< reciprocal of delta z [1/meters]
    REAL(ReKi)  :: InvMFFWS = 0      !< reciprocal of mean wind speed (MeanFFWS) [seconds/meter]
    REAL(ReKi)  :: MeanFFWS = 0      !< Mean wind speed (as defined in FF file), not necessarily of the portion used [meters/second]
    REAL(ReKi)  :: TotalTime = 0      !< The total time of the simulation [seconds]
    INTEGER(IntKi)  :: NFFComp = 3      !< Number of wind components [-]
    INTEGER(IntKi)  :: NFFSteps = 0      !< Number of time steps in the FF array [-]
    INTEGER(IntKi)  :: NYGrids = 0      !< Number of points in the lateral (y) direction of the grids [-]
    INTEGER(IntKi)  :: NZGrids = 0      !< Number of points in the vertical (z) direction of the grids [-]
    INTEGER(IntKi)  :: NTGrids = 0      !< Number of points in the vertical (z) direction on the tower (below the grids) [-]
    INTEGER(IntKi)  :: WindFileFormat      !< Binary file format description number [-]
  END TYPE IfW_TSFFWind_ParameterType
! =======================
CONTAINS
 SUBROUTINE IfW_TSFFWind_CopyInitInput( SrcInitInputData, DstInitInputData, CtrlCode, ErrStat, ErrMsg )
   TYPE(IfW_TSFFWind_InitInputType), INTENT(IN) :: SrcInitInputData
   TYPE(IfW_TSFFWind_InitInputType), INTENT(INOUT) :: DstInitInputData
   INTEGER(IntKi),  INTENT(IN   ) :: CtrlCode
   INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
   CHARACTER(*),    INTENT(  OUT) :: ErrMsg
! Local
   INTEGER(IntKi)                 :: i,j,k
   INTEGER(IntKi)                 :: i1, i1_l, i1_u  !  bounds (upper/lower) for an array dimension 1
   INTEGER(IntKi)                 :: i2, i2_l, i2_u  !  bounds (upper/lower) for an array dimension 2
   INTEGER(IntKi)                 :: i3, i3_l, i3_u  !  bounds (upper/lower) for an array dimension 3
   INTEGER(IntKi)                 :: i4, i4_l, i4_u  !  bounds (upper/lower) for an array dimension 4
   INTEGER(IntKi)                 :: ErrStat2
   CHARACTER(ErrMsgLen)           :: ErrMsg2
   CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_CopyInitInput'
!
   ErrStat = ErrID_None
   ErrMsg  = ""
    DstInitInputData%WindFileName = SrcInitInputData%WindFileName
    DstInitInputData%SumFileUnit = SrcInitInputData%SumFileUnit
 END SUBROUTINE IfW_TSFFWind_CopyInitInput

 SUBROUTINE IfW_TSFFWind_DestroyInitInput( InitInputData, ErrStat, ErrMsg )
  TYPE(IfW_TSFFWind_InitInputType), INTENT(INOUT) :: InitInputData
  INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
  CHARACTER(*),    INTENT(  OUT) :: ErrMsg
  CHARACTER(*),    PARAMETER :: RoutineName = 'IfW_TSFFWind_DestroyInitInput'
  INTEGER(IntKi)                 :: i, i1, i2, i3, i4, i5
!
  ErrStat = ErrID_None
  ErrMsg  = ""
 END SUBROUTINE IfW_TSFFWind_DestroyInitInput

 SUBROUTINE IfW_TSFFWind_PackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, SizeOnly )
  REAL(ReKi),       ALLOCATABLE, INTENT(  OUT) :: ReKiBuf(:)
  REAL(DbKi),       ALLOCATABLE, INTENT(  OUT) :: DbKiBuf(:)
  INTEGER(IntKi),   ALLOCATABLE, INTENT(  OUT) :: IntKiBuf(:)
  TYPE(IfW_TSFFWind_InitInputType),  INTENT(IN) :: InData
  INTEGER(IntKi),   INTENT(  OUT) :: ErrStat
  CHARACTER(*),     INTENT(  OUT) :: ErrMsg
  LOGICAL,OPTIONAL, INTENT(IN   ) :: SizeOnly
    ! Local variables
  INTEGER(IntKi)                 :: Re_BufSz
  INTEGER(IntKi)                 :: Re_Xferred
  INTEGER(IntKi)                 :: Db_BufSz
  INTEGER(IntKi)                 :: Db_Xferred
  INTEGER(IntKi)                 :: Int_BufSz
  INTEGER(IntKi)                 :: Int_Xferred
  INTEGER(IntKi)                 :: i,i1,i2,i3,i4,i5
  LOGICAL                        :: OnlySize ! if present and true, do not pack, just allocate buffers
  INTEGER(IntKi)                 :: ErrStat2
  CHARACTER(ErrMsgLen)           :: ErrMsg2
  CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_PackInitInput'
 ! buffers to store subtypes, if any
  REAL(ReKi),      ALLOCATABLE   :: Re_Buf(:)
  REAL(DbKi),      ALLOCATABLE   :: Db_Buf(:)
  INTEGER(IntKi),  ALLOCATABLE   :: Int_Buf(:)

  OnlySize = .FALSE.
  IF ( PRESENT(SizeOnly) ) THEN
    OnlySize = SizeOnly
  ENDIF
    !
  ErrStat = ErrID_None
  ErrMsg  = ""
  Re_BufSz  = 0
  Db_BufSz  = 0
  Int_BufSz  = 0
      Int_BufSz  = Int_BufSz  + 1*LEN(InData%WindFileName)  ! WindFileName
      Int_BufSz  = Int_BufSz  + 1  ! SumFileUnit
  IF ( Re_BufSz  .GT. 0 ) THEN
     ALLOCATE( ReKiBuf(  Re_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating ReKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF ( Db_BufSz  .GT. 0 ) THEN
     ALLOCATE( DbKiBuf(  Db_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating DbKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF ( Int_BufSz  .GT. 0 ) THEN
     ALLOCATE( IntKiBuf(  Int_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating IntKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF(OnlySize) RETURN ! return early if only trying to allocate buffers (not pack them)

  Re_Xferred  = 1
  Db_Xferred  = 1
  Int_Xferred = 1

        DO I = 1, LEN(InData%WindFileName)
          IntKiBuf(Int_Xferred) = ICHAR(InData%WindFileName(I:I), IntKi)
          Int_Xferred = Int_Xferred   + 1
        END DO ! I
      IntKiBuf ( Int_Xferred:Int_Xferred+(1)-1 ) = InData%SumFileUnit
      Int_Xferred   = Int_Xferred   + 1
 END SUBROUTINE IfW_TSFFWind_PackInitInput

 SUBROUTINE IfW_TSFFWind_UnPackInitInput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg )
  REAL(ReKi),      ALLOCATABLE, INTENT(IN   ) :: ReKiBuf(:)
  REAL(DbKi),      ALLOCATABLE, INTENT(IN   ) :: DbKiBuf(:)
  INTEGER(IntKi),  ALLOCATABLE, INTENT(IN   ) :: IntKiBuf(:)
  TYPE(IfW_TSFFWind_InitInputType), INTENT(INOUT) :: OutData
  INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
  CHARACTER(*),    INTENT(  OUT) :: ErrMsg
    ! Local variables
  INTEGER(IntKi)                 :: Buf_size
  INTEGER(IntKi)                 :: Re_Xferred
  INTEGER(IntKi)                 :: Db_Xferred
  INTEGER(IntKi)                 :: Int_Xferred
  INTEGER(IntKi)                 :: i
  LOGICAL                        :: mask0
  LOGICAL, ALLOCATABLE           :: mask1(:)
  LOGICAL, ALLOCATABLE           :: mask2(:,:)
  LOGICAL, ALLOCATABLE           :: mask3(:,:,:)
  LOGICAL, ALLOCATABLE           :: mask4(:,:,:,:)
  LOGICAL, ALLOCATABLE           :: mask5(:,:,:,:,:)
  INTEGER(IntKi)                 :: i1, i1_l, i1_u  !  bounds (upper/lower) for an array dimension 1
  INTEGER(IntKi)                 :: i2, i2_l, i2_u  !  bounds (upper/lower) for an array dimension 2
  INTEGER(IntKi)                 :: i3, i3_l, i3_u  !  bounds (upper/lower) for an array dimension 3
  INTEGER(IntKi)                 :: i4, i4_l, i4_u  !  bounds (upper/lower) for an array dimension 4
  INTEGER(IntKi)                 :: ErrStat2
  CHARACTER(ErrMsgLen)           :: ErrMsg2
  CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_UnPackInitInput'
 ! buffers to store meshes, if any
  REAL(ReKi),      ALLOCATABLE   :: Re_Buf(:)
  REAL(DbKi),      ALLOCATABLE   :: Db_Buf(:)
  INTEGER(IntKi),  ALLOCATABLE   :: Int_Buf(:)
    !
  ErrStat = ErrID_None
  ErrMsg  = ""
  Re_Xferred  = 1
  Db_Xferred  = 1
  Int_Xferred  = 1
      DO I = 1, LEN(OutData%WindFileName)
        OutData%WindFileName(I:I) = CHAR(IntKiBuf(Int_Xferred))
        Int_Xferred = Int_Xferred   + 1
      END DO ! I
      OutData%SumFileUnit = IntKiBuf( Int_Xferred )
      Int_Xferred   = Int_Xferred + 1
 END SUBROUTINE IfW_TSFFWind_UnPackInitInput

 SUBROUTINE IfW_TSFFWind_CopyInitOutput( SrcInitOutputData, DstInitOutputData, CtrlCode, ErrStat, ErrMsg )
   TYPE(IfW_TSFFWind_InitOutputType), INTENT(IN) :: SrcInitOutputData
   TYPE(IfW_TSFFWind_InitOutputType), INTENT(INOUT) :: DstInitOutputData
   INTEGER(IntKi),  INTENT(IN   ) :: CtrlCode
   INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
   CHARACTER(*),    INTENT(  OUT) :: ErrMsg
! Local
   INTEGER(IntKi)                 :: i,j,k
   INTEGER(IntKi)                 :: ErrStat2
   CHARACTER(ErrMsgLen)           :: ErrMsg2
   CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_CopyInitOutput'
!
   ErrStat = ErrID_None
   ErrMsg  = ""
      CALL NWTC_Library_Copyprogdesc( SrcInitOutputData%Ver, DstInitOutputData%Ver, CtrlCode, ErrStat2, ErrMsg2 )
         CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg,RoutineName)
         IF (ErrStat>=AbortErrLev) RETURN
 END SUBROUTINE IfW_TSFFWind_CopyInitOutput

 SUBROUTINE IfW_TSFFWind_DestroyInitOutput( InitOutputData, ErrStat, ErrMsg )
  TYPE(IfW_TSFFWind_InitOutputType), INTENT(INOUT) :: InitOutputData
  INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
  CHARACTER(*),    INTENT(  OUT) :: ErrMsg
  CHARACTER(*),    PARAMETER :: RoutineName = 'IfW_TSFFWind_DestroyInitOutput'
  INTEGER(IntKi)                 :: i, i1, i2, i3, i4, i5
!
  ErrStat = ErrID_None
  ErrMsg  = ""
  CALL NWTC_Library_Destroyprogdesc( InitOutputData%Ver, ErrStat, ErrMsg )
 END SUBROUTINE IfW_TSFFWind_DestroyInitOutput

 SUBROUTINE IfW_TSFFWind_PackInitOutput( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, SizeOnly )
  REAL(ReKi),       ALLOCATABLE, INTENT(  OUT) :: ReKiBuf(:)
  REAL(DbKi),       ALLOCATABLE, INTENT(  OUT) :: DbKiBuf(:)
  INTEGER(IntKi),   ALLOCATABLE, INTENT(  OUT) :: IntKiBuf(:)
  TYPE(IfW_TSFFWind_InitOutputType),  INTENT(IN) :: InData
  INTEGER(IntKi),   INTENT(  OUT) :: ErrStat
  CHARACTER(*),     INTENT(  OUT) :: ErrMsg
  LOGICAL,OPTIONAL, INTENT(IN   ) :: SizeOnly
    ! Local variables
  INTEGER(IntKi)                 :: Re_BufSz
  INTEGER(IntKi)                 :: Re_Xferred
  INTEGER(IntKi)                 :: Db_BufSz
  INTEGER(IntKi)                 :: Db_Xferred
  INTEGER(IntKi)                 :: Int_BufSz
  INTEGER(IntKi)                 :: Int_Xferred
  INTEGER(IntKi)                 :: i,i1,i2,i3,i4,i5
  LOGICAL                        :: OnlySize ! if present and true, do not pack, just allocate buffers
  INTEGER(IntKi)                 :: ErrStat2
  CHARACTER(ErrMsgLen)           :: ErrMsg2
  CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_PackInitOutput'
 ! buffers to store subtypes, if any
  REAL(ReKi),      ALLOCATABLE   :: Re_Buf(:)
  REAL(DbKi),      ALLOCATABLE   :: Db_Buf(:)
  INTEGER(IntKi),  ALLOCATABLE   :: Int_Buf(:)

  OnlySize = .FALSE.
  IF ( PRESENT(SizeOnly) ) THEN
    OnlySize = SizeOnly
  ENDIF
    !
  ErrStat = ErrID_None
  ErrMsg  = ""
  Re_BufSz  = 0
  Db_BufSz  = 0
  Int_BufSz  = 0
   ! Allocate buffers for subtypes, if any (we'll get sizes from these)
      Int_BufSz   = Int_BufSz + 3  ! Ver: size of buffers for each call to pack subtype
      CALL NWTC_Library_Packprogdesc( Re_Buf, Db_Buf, Int_Buf, InData%Ver, ErrStat2, ErrMsg2, .TRUE. ) ! Ver
        CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
        IF (ErrStat >= AbortErrLev) RETURN

      IF(ALLOCATED(Re_Buf)) THEN ! Ver
         Re_BufSz  = Re_BufSz  + SIZE( Re_Buf  )
         DEALLOCATE(Re_Buf)
      END IF
      IF(ALLOCATED(Db_Buf)) THEN ! Ver
         Db_BufSz  = Db_BufSz  + SIZE( Db_Buf  )
         DEALLOCATE(Db_Buf)
      END IF
      IF(ALLOCATED(Int_Buf)) THEN ! Ver
         Int_BufSz = Int_BufSz + SIZE( Int_Buf )
         DEALLOCATE(Int_Buf)
      END IF
  IF ( Re_BufSz  .GT. 0 ) THEN
     ALLOCATE( ReKiBuf(  Re_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating ReKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF ( Db_BufSz  .GT. 0 ) THEN
     ALLOCATE( DbKiBuf(  Db_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating DbKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF ( Int_BufSz  .GT. 0 ) THEN
     ALLOCATE( IntKiBuf(  Int_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating IntKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF(OnlySize) RETURN ! return early if only trying to allocate buffers (not pack them)

  Re_Xferred  = 1
  Db_Xferred  = 1
  Int_Xferred = 1

      CALL NWTC_Library_Packprogdesc( Re_Buf, Db_Buf, Int_Buf, InData%Ver, ErrStat2, ErrMsg2, OnlySize ) ! Ver
        CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
        IF (ErrStat >= AbortErrLev) RETURN

      IF(ALLOCATED(Re_Buf)) THEN
        IntKiBuf( Int_Xferred ) = SIZE(Re_Buf); Int_Xferred = Int_Xferred + 1
        IF (SIZE(Re_Buf) > 0) ReKiBuf( Re_Xferred:Re_Xferred+SIZE(Re_Buf)-1 ) = Re_Buf
        Re_Xferred = Re_Xferred + SIZE(Re_Buf)
        DEALLOCATE(Re_Buf)
      ELSE
        IntKiBuf( Int_Xferred ) = 0; Int_Xferred = Int_Xferred + 1
      ENDIF
      IF(ALLOCATED(Db_Buf)) THEN
        IntKiBuf( Int_Xferred ) = SIZE(Db_Buf); Int_Xferred = Int_Xferred + 1
        IF (SIZE(Db_Buf) > 0) DbKiBuf( Db_Xferred:Db_Xferred+SIZE(Db_Buf)-1 ) = Db_Buf
        Db_Xferred = Db_Xferred + SIZE(Db_Buf)
        DEALLOCATE(Db_Buf)
      ELSE
        IntKiBuf( Int_Xferred ) = 0; Int_Xferred = Int_Xferred + 1
      ENDIF
      IF(ALLOCATED(Int_Buf)) THEN
        IntKiBuf( Int_Xferred ) = SIZE(Int_Buf); Int_Xferred = Int_Xferred + 1
        IF (SIZE(Int_Buf) > 0) IntKiBuf( Int_Xferred:Int_Xferred+SIZE(Int_Buf)-1 ) = Int_Buf
        Int_Xferred = Int_Xferred + SIZE(Int_Buf)
        DEALLOCATE(Int_Buf)
      ELSE
        IntKiBuf( Int_Xferred ) = 0; Int_Xferred = Int_Xferred + 1
      ENDIF
 END SUBROUTINE IfW_TSFFWind_PackInitOutput

 SUBROUTINE IfW_TSFFWind_UnPackInitOutput( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg )
  REAL(ReKi),      ALLOCATABLE, INTENT(IN   ) :: ReKiBuf(:)
  REAL(DbKi),      ALLOCATABLE, INTENT(IN   ) :: DbKiBuf(:)
  INTEGER(IntKi),  ALLOCATABLE, INTENT(IN   ) :: IntKiBuf(:)
  TYPE(IfW_TSFFWind_InitOutputType), INTENT(INOUT) :: OutData
  INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
  CHARACTER(*),    INTENT(  OUT) :: ErrMsg
    ! Local variables
  INTEGER(IntKi)                 :: Buf_size
  INTEGER(IntKi)                 :: Re_Xferred
  INTEGER(IntKi)                 :: Db_Xferred
  INTEGER(IntKi)                 :: Int_Xferred
  INTEGER(IntKi)                 :: i
  LOGICAL                        :: mask0
  LOGICAL, ALLOCATABLE           :: mask1(:)
  LOGICAL, ALLOCATABLE           :: mask2(:,:)
  LOGICAL, ALLOCATABLE           :: mask3(:,:,:)
  LOGICAL, ALLOCATABLE           :: mask4(:,:,:,:)
  LOGICAL, ALLOCATABLE           :: mask5(:,:,:,:,:)
  INTEGER(IntKi)                 :: ErrStat2
  CHARACTER(ErrMsgLen)           :: ErrMsg2
  CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_UnPackInitOutput'
 ! buffers to store meshes, if any
  REAL(ReKi),      ALLOCATABLE   :: Re_Buf(:)
  REAL(DbKi),      ALLOCATABLE   :: Db_Buf(:)
  INTEGER(IntKi),  ALLOCATABLE   :: Int_Buf(:)
    !
  ErrStat = ErrID_None
  ErrMsg  = ""
  Re_Xferred  = 1
  Db_Xferred  = 1
  Int_Xferred  = 1
      Buf_size=IntKiBuf( Int_Xferred )
      Int_Xferred = Int_Xferred + 1
      IF(Buf_size > 0) THEN
        ALLOCATE(Re_Buf(Buf_size),STAT=ErrStat2)
        IF (ErrStat2 /= 0) THEN
           CALL SetErrStat(ErrID_Fatal, 'Error allocating Re_Buf.', ErrStat, ErrMsg,RoutineName)
           RETURN
        END IF
        Re_Buf = ReKiBuf( Re_Xferred:Re_Xferred+Buf_size-1 )
        Re_Xferred = Re_Xferred + Buf_size
      END IF
      Buf_size=IntKiBuf( Int_Xferred )
      Int_Xferred = Int_Xferred + 1
      IF(Buf_size > 0) THEN
        ALLOCATE(Db_Buf(Buf_size),STAT=ErrStat2)
        IF (ErrStat2 /= 0) THEN
           CALL SetErrStat(ErrID_Fatal, 'Error allocating Db_Buf.', ErrStat, ErrMsg,RoutineName)
           RETURN
        END IF
        Db_Buf = DbKiBuf( Db_Xferred:Db_Xferred+Buf_size-1 )
        Db_Xferred = Db_Xferred + Buf_size
      END IF
      Buf_size=IntKiBuf( Int_Xferred )
      Int_Xferred = Int_Xferred + 1
      IF(Buf_size > 0) THEN
        ALLOCATE(Int_Buf(Buf_size),STAT=ErrStat2)
        IF (ErrStat2 /= 0) THEN
           CALL SetErrStat(ErrID_Fatal, 'Error allocating Int_Buf.', ErrStat, ErrMsg,RoutineName)
           RETURN
        END IF
        Int_Buf = IntKiBuf( Int_Xferred:Int_Xferred+Buf_size-1 )
        Int_Xferred = Int_Xferred + Buf_size
      END IF
      CALL NWTC_Library_Unpackprogdesc( Re_Buf, Db_Buf, Int_Buf, OutData%Ver, ErrStat2, ErrMsg2 ) ! Ver
        CALL SetErrStat(ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName)
        IF (ErrStat >= AbortErrLev) RETURN

      IF(ALLOCATED(Re_Buf )) DEALLOCATE(Re_Buf )
      IF(ALLOCATED(Db_Buf )) DEALLOCATE(Db_Buf )
      IF(ALLOCATED(Int_Buf)) DEALLOCATE(Int_Buf)
 END SUBROUTINE IfW_TSFFWind_UnPackInitOutput

 SUBROUTINE IfW_TSFFWind_CopyMisc( SrcMiscData, DstMiscData, CtrlCode, ErrStat, ErrMsg )
   TYPE(IfW_TSFFWind_MiscVarType), INTENT(IN) :: SrcMiscData
   TYPE(IfW_TSFFWind_MiscVarType), INTENT(INOUT) :: DstMiscData
   INTEGER(IntKi),  INTENT(IN   ) :: CtrlCode
   INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
   CHARACTER(*),    INTENT(  OUT) :: ErrMsg
! Local
   INTEGER(IntKi)                 :: i,j,k
   INTEGER(IntKi)                 :: ErrStat2
   CHARACTER(ErrMsgLen)           :: ErrMsg2
   CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_CopyMisc'
!
   ErrStat = ErrID_None
   ErrMsg  = ""
    DstMiscData%TimeIndex = SrcMiscData%TimeIndex
 END SUBROUTINE IfW_TSFFWind_CopyMisc

 SUBROUTINE IfW_TSFFWind_DestroyMisc( MiscData, ErrStat, ErrMsg )
  TYPE(IfW_TSFFWind_MiscVarType), INTENT(INOUT) :: MiscData
  INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
  CHARACTER(*),    INTENT(  OUT) :: ErrMsg
  CHARACTER(*),    PARAMETER :: RoutineName = 'IfW_TSFFWind_DestroyMisc'
  INTEGER(IntKi)                 :: i, i1, i2, i3, i4, i5
!
  ErrStat = ErrID_None
  ErrMsg  = ""
 END SUBROUTINE IfW_TSFFWind_DestroyMisc

 SUBROUTINE IfW_TSFFWind_PackMisc( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, SizeOnly )
  REAL(ReKi),       ALLOCATABLE, INTENT(  OUT) :: ReKiBuf(:)
  REAL(DbKi),       ALLOCATABLE, INTENT(  OUT) :: DbKiBuf(:)
  INTEGER(IntKi),   ALLOCATABLE, INTENT(  OUT) :: IntKiBuf(:)
  TYPE(IfW_TSFFWind_MiscVarType),  INTENT(IN) :: InData
  INTEGER(IntKi),   INTENT(  OUT) :: ErrStat
  CHARACTER(*),     INTENT(  OUT) :: ErrMsg
  LOGICAL,OPTIONAL, INTENT(IN   ) :: SizeOnly
    ! Local variables
  INTEGER(IntKi)                 :: Re_BufSz
  INTEGER(IntKi)                 :: Re_Xferred
  INTEGER(IntKi)                 :: Db_BufSz
  INTEGER(IntKi)                 :: Db_Xferred
  INTEGER(IntKi)                 :: Int_BufSz
  INTEGER(IntKi)                 :: Int_Xferred
  INTEGER(IntKi)                 :: i,i1,i2,i3,i4,i5
  LOGICAL                        :: OnlySize ! if present and true, do not pack, just allocate buffers
  INTEGER(IntKi)                 :: ErrStat2
  CHARACTER(ErrMsgLen)           :: ErrMsg2
  CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_PackMisc'
 ! buffers to store subtypes, if any
  REAL(ReKi),      ALLOCATABLE   :: Re_Buf(:)
  REAL(DbKi),      ALLOCATABLE   :: Db_Buf(:)
  INTEGER(IntKi),  ALLOCATABLE   :: Int_Buf(:)

  OnlySize = .FALSE.
  IF ( PRESENT(SizeOnly) ) THEN
    OnlySize = SizeOnly
  ENDIF
    !
  ErrStat = ErrID_None
  ErrMsg  = ""
  Re_BufSz  = 0
  Db_BufSz  = 0
  Int_BufSz  = 0
      Int_BufSz  = Int_BufSz  + 1  ! TimeIndex
  IF ( Re_BufSz  .GT. 0 ) THEN
     ALLOCATE( ReKiBuf(  Re_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating ReKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF ( Db_BufSz  .GT. 0 ) THEN
     ALLOCATE( DbKiBuf(  Db_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating DbKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF ( Int_BufSz  .GT. 0 ) THEN
     ALLOCATE( IntKiBuf(  Int_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating IntKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF(OnlySize) RETURN ! return early if only trying to allocate buffers (not pack them)

  Re_Xferred  = 1
  Db_Xferred  = 1
  Int_Xferred = 1

      IntKiBuf ( Int_Xferred:Int_Xferred+(1)-1 ) = InData%TimeIndex
      Int_Xferred   = Int_Xferred   + 1
 END SUBROUTINE IfW_TSFFWind_PackMisc

 SUBROUTINE IfW_TSFFWind_UnPackMisc( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg )
  REAL(ReKi),      ALLOCATABLE, INTENT(IN   ) :: ReKiBuf(:)
  REAL(DbKi),      ALLOCATABLE, INTENT(IN   ) :: DbKiBuf(:)
  INTEGER(IntKi),  ALLOCATABLE, INTENT(IN   ) :: IntKiBuf(:)
  TYPE(IfW_TSFFWind_MiscVarType), INTENT(INOUT) :: OutData
  INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
  CHARACTER(*),    INTENT(  OUT) :: ErrMsg
    ! Local variables
  INTEGER(IntKi)                 :: Buf_size
  INTEGER(IntKi)                 :: Re_Xferred
  INTEGER(IntKi)                 :: Db_Xferred
  INTEGER(IntKi)                 :: Int_Xferred
  INTEGER(IntKi)                 :: i
  LOGICAL                        :: mask0
  LOGICAL, ALLOCATABLE           :: mask1(:)
  LOGICAL, ALLOCATABLE           :: mask2(:,:)
  LOGICAL, ALLOCATABLE           :: mask3(:,:,:)
  LOGICAL, ALLOCATABLE           :: mask4(:,:,:,:)
  LOGICAL, ALLOCATABLE           :: mask5(:,:,:,:,:)
  INTEGER(IntKi)                 :: ErrStat2
  CHARACTER(ErrMsgLen)           :: ErrMsg2
  CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_UnPackMisc'
 ! buffers to store meshes, if any
  REAL(ReKi),      ALLOCATABLE   :: Re_Buf(:)
  REAL(DbKi),      ALLOCATABLE   :: Db_Buf(:)
  INTEGER(IntKi),  ALLOCATABLE   :: Int_Buf(:)
    !
  ErrStat = ErrID_None
  ErrMsg  = ""
  Re_Xferred  = 1
  Db_Xferred  = 1
  Int_Xferred  = 1
      OutData%TimeIndex = IntKiBuf( Int_Xferred )
      Int_Xferred   = Int_Xferred + 1
 END SUBROUTINE IfW_TSFFWind_UnPackMisc

 SUBROUTINE IfW_TSFFWind_CopyParam( SrcParamData, DstParamData, CtrlCode, ErrStat, ErrMsg )
   TYPE(IfW_TSFFWind_ParameterType), INTENT(IN) :: SrcParamData
   TYPE(IfW_TSFFWind_ParameterType), INTENT(INOUT) :: DstParamData
   INTEGER(IntKi),  INTENT(IN   ) :: CtrlCode
   INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
   CHARACTER(*),    INTENT(  OUT) :: ErrMsg
! Local
   INTEGER(IntKi)                 :: i,j,k
   INTEGER(IntKi)                 :: i1, i1_l, i1_u  !  bounds (upper/lower) for an array dimension 1
   INTEGER(IntKi)                 :: i2, i2_l, i2_u  !  bounds (upper/lower) for an array dimension 2
   INTEGER(IntKi)                 :: i3, i3_l, i3_u  !  bounds (upper/lower) for an array dimension 3
   INTEGER(IntKi)                 :: i4, i4_l, i4_u  !  bounds (upper/lower) for an array dimension 4
   INTEGER(IntKi)                 :: ErrStat2
   CHARACTER(ErrMsgLen)           :: ErrMsg2
   CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_CopyParam'
!
   ErrStat = ErrID_None
   ErrMsg  = ""
    DstParamData%DT = SrcParamData%DT
    DstParamData%TowerDataExist = SrcParamData%TowerDataExist
    DstParamData%Periodic = SrcParamData%Periodic
IF (ALLOCATED(SrcParamData%FFData)) THEN
  i1_l = LBOUND(SrcParamData%FFData,1)
  i1_u = UBOUND(SrcParamData%FFData,1)
  i2_l = LBOUND(SrcParamData%FFData,2)
  i2_u = UBOUND(SrcParamData%FFData,2)
  i3_l = LBOUND(SrcParamData%FFData,3)
  i3_u = UBOUND(SrcParamData%FFData,3)
  i4_l = LBOUND(SrcParamData%FFData,4)
  i4_u = UBOUND(SrcParamData%FFData,4)
  IF (.NOT. ALLOCATED(DstParamData%FFData)) THEN
    ALLOCATE(DstParamData%FFData(i1_l:i1_u,i2_l:i2_u,i3_l:i3_u,i4_l:i4_u),STAT=ErrStat2)
    IF (ErrStat2 /= 0) THEN
      CALL SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%FFData.', ErrStat, ErrMsg,RoutineName)
      RETURN
    END IF
  END IF
    DstParamData%FFData = SrcParamData%FFData
ENDIF
IF (ALLOCATED(SrcParamData%FFTower)) THEN
  i1_l = LBOUND(SrcParamData%FFTower,1)
  i1_u = UBOUND(SrcParamData%FFTower,1)
  i2_l = LBOUND(SrcParamData%FFTower,2)
  i2_u = UBOUND(SrcParamData%FFTower,2)
  i3_l = LBOUND(SrcParamData%FFTower,3)
  i3_u = UBOUND(SrcParamData%FFTower,3)
  IF (.NOT. ALLOCATED(DstParamData%FFTower)) THEN
    ALLOCATE(DstParamData%FFTower(i1_l:i1_u,i2_l:i2_u,i3_l:i3_u),STAT=ErrStat2)
    IF (ErrStat2 /= 0) THEN
      CALL SetErrStat(ErrID_Fatal, 'Error allocating DstParamData%FFTower.', ErrStat, ErrMsg,RoutineName)
      RETURN
    END IF
  END IF
    DstParamData%FFTower = SrcParamData%FFTower
ENDIF
    DstParamData%FFDTime = SrcParamData%FFDTime
    DstParamData%FFRate = SrcParamData%FFRate
    DstParamData%FFYHWid = SrcParamData%FFYHWid
    DstParamData%FFZHWid = SrcParamData%FFZHWid
    DstParamData%RefHt = SrcParamData%RefHt
    DstParamData%GridBase = SrcParamData%GridBase
    DstParamData%InitXPosition = SrcParamData%InitXPosition
    DstParamData%InvFFYD = SrcParamData%InvFFYD
    DstParamData%InvFFZD = SrcParamData%InvFFZD
    DstParamData%InvMFFWS = SrcParamData%InvMFFWS
    DstParamData%MeanFFWS = SrcParamData%MeanFFWS
    DstParamData%TotalTime = SrcParamData%TotalTime
    DstParamData%NFFComp = SrcParamData%NFFComp
    DstParamData%NFFSteps = SrcParamData%NFFSteps
    DstParamData%NYGrids = SrcParamData%NYGrids
    DstParamData%NZGrids = SrcParamData%NZGrids
    DstParamData%NTGrids = SrcParamData%NTGrids
    DstParamData%WindFileFormat = SrcParamData%WindFileFormat
 END SUBROUTINE IfW_TSFFWind_CopyParam

 SUBROUTINE IfW_TSFFWind_DestroyParam( ParamData, ErrStat, ErrMsg )
  TYPE(IfW_TSFFWind_ParameterType), INTENT(INOUT) :: ParamData
  INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
  CHARACTER(*),    INTENT(  OUT) :: ErrMsg
  CHARACTER(*),    PARAMETER :: RoutineName = 'IfW_TSFFWind_DestroyParam'
  INTEGER(IntKi)                 :: i, i1, i2, i3, i4, i5
!
  ErrStat = ErrID_None
  ErrMsg  = ""
IF (ALLOCATED(ParamData%FFData)) THEN
  DEALLOCATE(ParamData%FFData)
ENDIF
IF (ALLOCATED(ParamData%FFTower)) THEN
  DEALLOCATE(ParamData%FFTower)
ENDIF
 END SUBROUTINE IfW_TSFFWind_DestroyParam

 SUBROUTINE IfW_TSFFWind_PackParam( ReKiBuf, DbKiBuf, IntKiBuf, Indata, ErrStat, ErrMsg, SizeOnly )
  REAL(ReKi),       ALLOCATABLE, INTENT(  OUT) :: ReKiBuf(:)
  REAL(DbKi),       ALLOCATABLE, INTENT(  OUT) :: DbKiBuf(:)
  INTEGER(IntKi),   ALLOCATABLE, INTENT(  OUT) :: IntKiBuf(:)
  TYPE(IfW_TSFFWind_ParameterType),  INTENT(IN) :: InData
  INTEGER(IntKi),   INTENT(  OUT) :: ErrStat
  CHARACTER(*),     INTENT(  OUT) :: ErrMsg
  LOGICAL,OPTIONAL, INTENT(IN   ) :: SizeOnly
    ! Local variables
  INTEGER(IntKi)                 :: Re_BufSz
  INTEGER(IntKi)                 :: Re_Xferred
  INTEGER(IntKi)                 :: Db_BufSz
  INTEGER(IntKi)                 :: Db_Xferred
  INTEGER(IntKi)                 :: Int_BufSz
  INTEGER(IntKi)                 :: Int_Xferred
  INTEGER(IntKi)                 :: i,i1,i2,i3,i4,i5
  LOGICAL                        :: OnlySize ! if present and true, do not pack, just allocate buffers
  INTEGER(IntKi)                 :: ErrStat2
  CHARACTER(ErrMsgLen)           :: ErrMsg2
  CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_PackParam'
 ! buffers to store subtypes, if any
  REAL(ReKi),      ALLOCATABLE   :: Re_Buf(:)
  REAL(DbKi),      ALLOCATABLE   :: Db_Buf(:)
  INTEGER(IntKi),  ALLOCATABLE   :: Int_Buf(:)

  OnlySize = .FALSE.
  IF ( PRESENT(SizeOnly) ) THEN
    OnlySize = SizeOnly
  ENDIF
    !
  ErrStat = ErrID_None
  ErrMsg  = ""
  Re_BufSz  = 0
  Db_BufSz  = 0
  Int_BufSz  = 0
      Db_BufSz   = Db_BufSz   + 1  ! DT
      Int_BufSz  = Int_BufSz  + 1  ! TowerDataExist
      Int_BufSz  = Int_BufSz  + 1  ! Periodic
  Int_BufSz   = Int_BufSz   + 1     ! FFData allocated yes/no
  IF ( ALLOCATED(InData%FFData) ) THEN
    Int_BufSz   = Int_BufSz   + 2*4  ! FFData upper/lower bounds for each dimension
      Re_BufSz   = Re_BufSz   + SIZE(InData%FFData)  ! FFData
  END IF
  Int_BufSz   = Int_BufSz   + 1     ! FFTower allocated yes/no
  IF ( ALLOCATED(InData%FFTower) ) THEN
    Int_BufSz   = Int_BufSz   + 2*3  ! FFTower upper/lower bounds for each dimension
      Re_BufSz   = Re_BufSz   + SIZE(InData%FFTower)  ! FFTower
  END IF
      Re_BufSz   = Re_BufSz   + 1  ! FFDTime
      Re_BufSz   = Re_BufSz   + 1  ! FFRate
      Re_BufSz   = Re_BufSz   + 1  ! FFYHWid
      Re_BufSz   = Re_BufSz   + 1  ! FFZHWid
      Re_BufSz   = Re_BufSz   + 1  ! RefHt
      Re_BufSz   = Re_BufSz   + 1  ! GridBase
      Re_BufSz   = Re_BufSz   + 1  ! InitXPosition
      Re_BufSz   = Re_BufSz   + 1  ! InvFFYD
      Re_BufSz   = Re_BufSz   + 1  ! InvFFZD
      Re_BufSz   = Re_BufSz   + 1  ! InvMFFWS
      Re_BufSz   = Re_BufSz   + 1  ! MeanFFWS
      Re_BufSz   = Re_BufSz   + 1  ! TotalTime
      Int_BufSz  = Int_BufSz  + 1  ! NFFComp
      Int_BufSz  = Int_BufSz  + 1  ! NFFSteps
      Int_BufSz  = Int_BufSz  + 1  ! NYGrids
      Int_BufSz  = Int_BufSz  + 1  ! NZGrids
      Int_BufSz  = Int_BufSz  + 1  ! NTGrids
      Int_BufSz  = Int_BufSz  + 1  ! WindFileFormat
  IF ( Re_BufSz  .GT. 0 ) THEN
     ALLOCATE( ReKiBuf(  Re_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating ReKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF ( Db_BufSz  .GT. 0 ) THEN
     ALLOCATE( DbKiBuf(  Db_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating DbKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF ( Int_BufSz  .GT. 0 ) THEN
     ALLOCATE( IntKiBuf(  Int_BufSz  ), STAT=ErrStat2 )
     IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating IntKiBuf.', ErrStat, ErrMsg,RoutineName)
       RETURN
     END IF
  END IF
  IF(OnlySize) RETURN ! return early if only trying to allocate buffers (not pack them)

  Re_Xferred  = 1
  Db_Xferred  = 1
  Int_Xferred = 1

      DbKiBuf ( Db_Xferred:Db_Xferred+(1)-1 ) = InData%DT
      Db_Xferred   = Db_Xferred   + 1
      IntKiBuf ( Int_Xferred:Int_Xferred+1-1 ) = TRANSFER( InData%TowerDataExist , IntKiBuf(1), 1)
      Int_Xferred   = Int_Xferred   + 1
      IntKiBuf ( Int_Xferred:Int_Xferred+1-1 ) = TRANSFER( InData%Periodic , IntKiBuf(1), 1)
      Int_Xferred   = Int_Xferred   + 1
  IF ( .NOT. ALLOCATED(InData%FFData) ) THEN
    IntKiBuf( Int_Xferred ) = 0
    Int_Xferred = Int_Xferred + 1
  ELSE
    IntKiBuf( Int_Xferred ) = 1
    Int_Xferred = Int_Xferred + 1
    IntKiBuf( Int_Xferred    ) = LBOUND(InData%FFData,1)
    IntKiBuf( Int_Xferred + 1) = UBOUND(InData%FFData,1)
    Int_Xferred = Int_Xferred + 2
    IntKiBuf( Int_Xferred    ) = LBOUND(InData%FFData,2)
    IntKiBuf( Int_Xferred + 1) = UBOUND(InData%FFData,2)
    Int_Xferred = Int_Xferred + 2
    IntKiBuf( Int_Xferred    ) = LBOUND(InData%FFData,3)
    IntKiBuf( Int_Xferred + 1) = UBOUND(InData%FFData,3)
    Int_Xferred = Int_Xferred + 2
    IntKiBuf( Int_Xferred    ) = LBOUND(InData%FFData,4)
    IntKiBuf( Int_Xferred + 1) = UBOUND(InData%FFData,4)
    Int_Xferred = Int_Xferred + 2

      IF (SIZE(InData%FFData)>0) ReKiBuf ( Re_Xferred:Re_Xferred+(SIZE(InData%FFData))-1 ) = PACK(InData%FFData,.TRUE.)
      Re_Xferred   = Re_Xferred   + SIZE(InData%FFData)
  END IF
  IF ( .NOT. ALLOCATED(InData%FFTower) ) THEN
    IntKiBuf( Int_Xferred ) = 0
    Int_Xferred = Int_Xferred + 1
  ELSE
    IntKiBuf( Int_Xferred ) = 1
    Int_Xferred = Int_Xferred + 1
    IntKiBuf( Int_Xferred    ) = LBOUND(InData%FFTower,1)
    IntKiBuf( Int_Xferred + 1) = UBOUND(InData%FFTower,1)
    Int_Xferred = Int_Xferred + 2
    IntKiBuf( Int_Xferred    ) = LBOUND(InData%FFTower,2)
    IntKiBuf( Int_Xferred + 1) = UBOUND(InData%FFTower,2)
    Int_Xferred = Int_Xferred + 2
    IntKiBuf( Int_Xferred    ) = LBOUND(InData%FFTower,3)
    IntKiBuf( Int_Xferred + 1) = UBOUND(InData%FFTower,3)
    Int_Xferred = Int_Xferred + 2

      IF (SIZE(InData%FFTower)>0) ReKiBuf ( Re_Xferred:Re_Xferred+(SIZE(InData%FFTower))-1 ) = PACK(InData%FFTower,.TRUE.)
      Re_Xferred   = Re_Xferred   + SIZE(InData%FFTower)
  END IF
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%FFDTime
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%FFRate
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%FFYHWid
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%FFZHWid
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%RefHt
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%GridBase
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%InitXPosition
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%InvFFYD
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%InvFFZD
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%InvMFFWS
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%MeanFFWS
      Re_Xferred   = Re_Xferred   + 1
      ReKiBuf ( Re_Xferred:Re_Xferred+(1)-1 ) = InData%TotalTime
      Re_Xferred   = Re_Xferred   + 1
      IntKiBuf ( Int_Xferred:Int_Xferred+(1)-1 ) = InData%NFFComp
      Int_Xferred   = Int_Xferred   + 1
      IntKiBuf ( Int_Xferred:Int_Xferred+(1)-1 ) = InData%NFFSteps
      Int_Xferred   = Int_Xferred   + 1
      IntKiBuf ( Int_Xferred:Int_Xferred+(1)-1 ) = InData%NYGrids
      Int_Xferred   = Int_Xferred   + 1
      IntKiBuf ( Int_Xferred:Int_Xferred+(1)-1 ) = InData%NZGrids
      Int_Xferred   = Int_Xferred   + 1
      IntKiBuf ( Int_Xferred:Int_Xferred+(1)-1 ) = InData%NTGrids
      Int_Xferred   = Int_Xferred   + 1
      IntKiBuf ( Int_Xferred:Int_Xferred+(1)-1 ) = InData%WindFileFormat
      Int_Xferred   = Int_Xferred   + 1
 END SUBROUTINE IfW_TSFFWind_PackParam

 SUBROUTINE IfW_TSFFWind_UnPackParam( ReKiBuf, DbKiBuf, IntKiBuf, Outdata, ErrStat, ErrMsg )
  REAL(ReKi),      ALLOCATABLE, INTENT(IN   ) :: ReKiBuf(:)
  REAL(DbKi),      ALLOCATABLE, INTENT(IN   ) :: DbKiBuf(:)
  INTEGER(IntKi),  ALLOCATABLE, INTENT(IN   ) :: IntKiBuf(:)
  TYPE(IfW_TSFFWind_ParameterType), INTENT(INOUT) :: OutData
  INTEGER(IntKi),  INTENT(  OUT) :: ErrStat
  CHARACTER(*),    INTENT(  OUT) :: ErrMsg
    ! Local variables
  INTEGER(IntKi)                 :: Buf_size
  INTEGER(IntKi)                 :: Re_Xferred
  INTEGER(IntKi)                 :: Db_Xferred
  INTEGER(IntKi)                 :: Int_Xferred
  INTEGER(IntKi)                 :: i
  LOGICAL                        :: mask0
  LOGICAL, ALLOCATABLE           :: mask1(:)
  LOGICAL, ALLOCATABLE           :: mask2(:,:)
  LOGICAL, ALLOCATABLE           :: mask3(:,:,:)
  LOGICAL, ALLOCATABLE           :: mask4(:,:,:,:)
  LOGICAL, ALLOCATABLE           :: mask5(:,:,:,:,:)
  INTEGER(IntKi)                 :: i1, i1_l, i1_u  !  bounds (upper/lower) for an array dimension 1
  INTEGER(IntKi)                 :: i2, i2_l, i2_u  !  bounds (upper/lower) for an array dimension 2
  INTEGER(IntKi)                 :: i3, i3_l, i3_u  !  bounds (upper/lower) for an array dimension 3
  INTEGER(IntKi)                 :: i4, i4_l, i4_u  !  bounds (upper/lower) for an array dimension 4
  INTEGER(IntKi)                 :: ErrStat2
  CHARACTER(ErrMsgLen)           :: ErrMsg2
  CHARACTER(*), PARAMETER        :: RoutineName = 'IfW_TSFFWind_UnPackParam'
 ! buffers to store meshes, if any
  REAL(ReKi),      ALLOCATABLE   :: Re_Buf(:)
  REAL(DbKi),      ALLOCATABLE   :: Db_Buf(:)
  INTEGER(IntKi),  ALLOCATABLE   :: Int_Buf(:)
    !
  ErrStat = ErrID_None
  ErrMsg  = ""
  Re_Xferred  = 1
  Db_Xferred  = 1
  Int_Xferred  = 1
      OutData%DT = DbKiBuf( Db_Xferred )
      Db_Xferred   = Db_Xferred + 1
      OutData%TowerDataExist = TRANSFER( IntKiBuf( Int_Xferred ), mask0 )
      Int_Xferred   = Int_Xferred + 1
      OutData%Periodic = TRANSFER( IntKiBuf( Int_Xferred ), mask0 )
      Int_Xferred   = Int_Xferred + 1
  IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN  ! FFData not allocated
    Int_Xferred = Int_Xferred + 1
  ELSE
    Int_Xferred = Int_Xferred + 1
    i1_l = IntKiBuf( Int_Xferred    )
    i1_u = IntKiBuf( Int_Xferred + 1)
    Int_Xferred = Int_Xferred + 2
    i2_l = IntKiBuf( Int_Xferred    )
    i2_u = IntKiBuf( Int_Xferred + 1)
    Int_Xferred = Int_Xferred + 2
    i3_l = IntKiBuf( Int_Xferred    )
    i3_u = IntKiBuf( Int_Xferred + 1)
    Int_Xferred = Int_Xferred + 2
    i4_l = IntKiBuf( Int_Xferred    )
    i4_u = IntKiBuf( Int_Xferred + 1)
    Int_Xferred = Int_Xferred + 2
    IF (ALLOCATED(OutData%FFData)) DEALLOCATE(OutData%FFData)
    ALLOCATE(OutData%FFData(i1_l:i1_u,i2_l:i2_u,i3_l:i3_u,i4_l:i4_u),STAT=ErrStat2)
    IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%FFData.', ErrStat, ErrMsg,RoutineName)
       RETURN
    END IF
    ALLOCATE(mask4(i1_l:i1_u,i2_l:i2_u,i3_l:i3_u,i4_l:i4_u),STAT=ErrStat2)
    IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating mask4.', ErrStat, ErrMsg,RoutineName)
       RETURN
    END IF
    mask4 = .TRUE.
      IF (SIZE(OutData%FFData)>0) OutData%FFData = REAL( UNPACK(ReKiBuf( Re_Xferred:Re_Xferred+(SIZE(OutData%FFData))-1 ), mask4, 0.0_ReKi ), SiKi)
      Re_Xferred   = Re_Xferred   + SIZE(OutData%FFData)
    DEALLOCATE(mask4)
  END IF
  IF ( IntKiBuf( Int_Xferred ) == 0 ) THEN  ! FFTower not allocated
    Int_Xferred = Int_Xferred + 1
  ELSE
    Int_Xferred = Int_Xferred + 1
    i1_l = IntKiBuf( Int_Xferred    )
    i1_u = IntKiBuf( Int_Xferred + 1)
    Int_Xferred = Int_Xferred + 2
    i2_l = IntKiBuf( Int_Xferred    )
    i2_u = IntKiBuf( Int_Xferred + 1)
    Int_Xferred = Int_Xferred + 2
    i3_l = IntKiBuf( Int_Xferred    )
    i3_u = IntKiBuf( Int_Xferred + 1)
    Int_Xferred = Int_Xferred + 2
    IF (ALLOCATED(OutData%FFTower)) DEALLOCATE(OutData%FFTower)
    ALLOCATE(OutData%FFTower(i1_l:i1_u,i2_l:i2_u,i3_l:i3_u),STAT=ErrStat2)
    IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating OutData%FFTower.', ErrStat, ErrMsg,RoutineName)
       RETURN
    END IF
    ALLOCATE(mask3(i1_l:i1_u,i2_l:i2_u,i3_l:i3_u),STAT=ErrStat2)
    IF (ErrStat2 /= 0) THEN
       CALL SetErrStat(ErrID_Fatal, 'Error allocating mask3.', ErrStat, ErrMsg,RoutineName)
       RETURN
    END IF
    mask3 = .TRUE.
      IF (SIZE(OutData%FFTower)>0) OutData%FFTower = REAL( UNPACK(ReKiBuf( Re_Xferred:Re_Xferred+(SIZE(OutData%FFTower))-1 ), mask3, 0.0_ReKi ), SiKi)
      Re_Xferred   = Re_Xferred   + SIZE(OutData%FFTower)
    DEALLOCATE(mask3)
  END IF
      OutData%FFDTime = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%FFRate = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%FFYHWid = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%FFZHWid = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%RefHt = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%GridBase = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%InitXPosition = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%InvFFYD = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%InvFFZD = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%InvMFFWS = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%MeanFFWS = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%TotalTime = ReKiBuf( Re_Xferred )
      Re_Xferred   = Re_Xferred + 1
      OutData%NFFComp = IntKiBuf( Int_Xferred )
      Int_Xferred   = Int_Xferred + 1
      OutData%NFFSteps = IntKiBuf( Int_Xferred )
      Int_Xferred   = Int_Xferred + 1
      OutData%NYGrids = IntKiBuf( Int_Xferred )
      Int_Xferred   = Int_Xferred + 1
      OutData%NZGrids = IntKiBuf( Int_Xferred )
      Int_Xferred   = Int_Xferred + 1
      OutData%NTGrids = IntKiBuf( Int_Xferred )
      Int_Xferred   = Int_Xferred + 1
      OutData%WindFileFormat = IntKiBuf( Int_Xferred )
      Int_Xferred   = Int_Xferred + 1
 END SUBROUTINE IfW_TSFFWind_UnPackParam

END MODULE IfW_TSFFWind_Types
!ENDOFREGISTRYGENERATEDFILE


text:
*     MB03XP EXAMPLE PROGRAM TEXT
*     Copyright (c) 2002-2020 NICONET e.V.
*
*     .. Parameters ..
      DOUBLE PRECISION ZERO, ONE
      PARAMETER        ( ZERO = 0.0D0, ONE = 1.0D0 )
      INTEGER          NIN, NOUT
      PARAMETER        ( NIN = 5, NOUT = 6 )
      INTEGER          NMAX
      PARAMETER        ( NMAX = 200 )
      INTEGER          LDA, LDB, LDQ, LDRES, LDZ, LDWORK
      PARAMETER        ( LDA = NMAX, LDB = NMAX, LDQ = NMAX,
     $                   LDRES = NMAX, LDWORK = NMAX, LDZ = NMAX )
*     .. Local Scalars ..
      INTEGER          I, IHI, ILO, INFO, J, N
*     .. Local Arrays ..
      DOUBLE PRECISION A(LDA,NMAX), ALPHAI(NMAX), ALPHAR(NMAX),
     $                 B(LDA,NMAX), BETA(NMAX), DWORK(LDWORK),
     $                 Q(LDQ,NMAX), RES(LDRES,3*NMAX), Z(LDZ,NMAX)
*     .. External Functions ..
      DOUBLE PRECISION DLANGE
      EXTERNAL         DLANGE
*     .. External Subroutines ..
      EXTERNAL         DGEMM, MB03XP
*     .. Executable Statements ..
      WRITE ( NOUT, FMT = 99999 )
*     Skip the heading in the data file and read the data.
      READ ( NIN, FMT = '()' )
      READ ( NIN, FMT = * )  N, ILO, IHI
      IF( N.LE.0 .OR. N.GT.NMAX ) THEN
         WRITE ( NOUT, FMT = 99990 ) N
      ELSE
         READ ( NIN, FMT = * ) ( ( A(I,J), J = 1,N ), I = 1,N )
         CALL DLACPY( 'All', N, N, A, LDA, RES(1,N+1), LDRES )
         READ ( NIN, FMT = * ) ( ( B(I,J), J = 1,N ), I = 1,N )
         CALL DLACPY( 'All', N, N, B, LDB, RES(1,2*N+1), LDRES )
         CALL MB03XP( 'S', 'I', 'I', N, ILO, IHI, A, LDA, B, LDB, Q,
     $                LDQ, Z, LDZ, ALPHAR, ALPHAI, BETA, DWORK, LDWORK,
     $                INFO )
         IF ( INFO.NE.0 ) THEN
            WRITE ( NOUT, FMT = 99998 ) INFO
         ELSE
            WRITE ( NOUT, FMT = 99996 )
            DO 10  I = 1, N
               WRITE (NOUT, FMT = 99991) ( A(I,J), J = 1,N )
10          CONTINUE
            CALL DGEMM( 'No Transpose', 'No Transpose', N, N, N, ONE,
     $                  RES(1,N+1), LDRES, Z, LDZ, ZERO, RES, LDRES )
            CALL DGEMM( 'No Transpose', 'No Transpose', N, N, N, -ONE,
     $                  Q, LDQ, A, LDA, ONE, RES, LDRES )
            WRITE ( NOUT, FMT = 99989 ) DLANGE( 'Frobenius', N, N, RES,
     $                                          LDRES, DWORK )
            WRITE ( NOUT, FMT = 99995 )
            DO 20  I = 1, N
               WRITE (NOUT, FMT = 99991) ( B(I,J), J = 1,N )
20          CONTINUE
            CALL DGEMM( 'No Transpose', 'No Transpose', N, N, N, ONE,
     $                  RES(1,2*N+1), LDRES, Q, LDQ, ZERO, RES, LDRES )
            CALL DGEMM( 'No Transpose', 'No Transpose', N, N, N, -ONE,
     $                  Z, LDZ, B, LDB, ONE, RES, LDRES )
            WRITE ( NOUT, FMT = 99988 ) DLANGE( 'Frobenius', N, N, RES,
     $                                          LDRES, DWORK )
            WRITE ( NOUT, FMT = 99994 )
            DO 30  I = 1, N
               WRITE (NOUT, FMT = 99991) ( Q(I,J), J = 1,N )
30          CONTINUE
            CALL DGEMM( 'Transpose', 'No Transpose', N, N, N, ONE, Q,
     $                  LDQ, Q, LDQ, ONE, RES, LDRES )
            DO 40  I = 1, N
               RES(I,I) = RES(I,I) - ONE
40          CONTINUE
            WRITE ( NOUT, FMT = 99987 ) DLANGE( 'Frobenius', N, N, RES,
     $                                          LDRES, DWORK )
            WRITE ( NOUT, FMT = 99993 )
            DO 50  I = 1, N
               WRITE (NOUT, FMT = 99991) ( Z(I,J), J = 1,N )
50          CONTINUE
            CALL DGEMM( 'Transpose', 'No Transpose', N, N, N, ONE, Z,
     $                  LDZ, Z, LDZ, ONE, RES, LDRES )
            DO 60 I = 1, N
               RES(I,I) = RES(I,I) - ONE
60          CONTINUE
            WRITE ( NOUT, FMT = 99986 ) DLANGE( 'Frobenius', N, N, RES,
     $                                          LDRES, DWORK )
            WRITE ( NOUT, FMT = 99992 )
            DO 70  I = 1, N
               WRITE ( NOUT, FMT = 99991 )
     $                 ALPHAR(I), ALPHAI(I), BETA(I)
70          CONTINUE
         END IF
      END IF
*
      STOP
*
99999 FORMAT (' MB03XP EXAMPLE PROGRAM RESULTS',/1X)
99998 FORMAT (' INFO on exit from MB03XP = ',I2)
99996 FORMAT (' The reduced matrix A is ')
99995 FORMAT (/' The reduced matrix B is ')
99994 FORMAT (/' The orthogonal factor Q is ')
99993 FORMAT (/' The orthogonal factor Z is ')
99992 FORMAT (/4X,'ALPHAR',4X,'ALPHAI',4X,'BETA')
99991 FORMAT (1000(1X,F9.4))
99990 FORMAT (/' N is out of range.',/' N = ',I5)
99989 FORMAT (/' Residual: || A*Z - Q*S ||_F = ',G7.2)
99988 FORMAT (/' Residual: || B*Q - Z*T ||_F = ',G7.2)
99987 FORMAT (/' Orthogonality of Q: || Q''*Q - I ||_F = ',G7.2)
99986 FORMAT (/' Orthogonality of Z: || Z''*Z - I ||_F = ',G7.2)
      END

ext: f

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference