c     slp_both_adjusted.f
c
c     Calculate adjusted SLP and SLP anomaly
c
c     Datasets are adjusted with the model's historical run 1
c     and ERA interim's mean and standard deviation
c
      use commondata, only : leapday,ensmem,threesixty,mx,my,seasons,
     & days,tmp_path,in_path,slp_mean_path,slp_std_path,mean_path,
     & era_slp_mean,era_slp_std,date,time

      implicit none

      real :: values(20000),results(20000)
      real :: p(200,100),q(20000)
      real :: ave(20000,4),std(20000,4)
      real :: eraave(20000,4),erastd(20000,4)
      integer :: vari(6),valid,leap,runs,mem
      integer :: iyear,imonth,iday,ihour,id,t,k,i,j,fid
      character(len=200) :: snd
      character(len=50) :: sne(3),snl(3),sno(2),snoanom(2)
      character(len=200) :: smem,input,shmodel,model

      data sne /'_hist_r1','_SLP_Mean1981-2000_','.txt'/
      data snl /'_hist_r1','_SLP_Std1981-2000_','.txt'/
      data sno /'_r','_SLP_Adjusted.dat'/
      data snoanom /'_r','_SLPanom_Adjusted.dat'/

c      call idate(date)
c      call itime(time)
c      write(*,1000) date(2),date(1),date(3)+2000,time
c1000  format('slp_both_adjusted start at: Date ', i2.2, '/', i2.2, 
c     &         '/', i4.4,'; time ',i2.2, ':', i2.2, ':', i2.2 )

      call get_command_argument(1,input)
      read(input,*) model
      call get_command_argument(2,input)
      read(input,*) runs
      call get_command_argument(3,input)
      read(input,*) vari(1)
      call get_command_argument(4,input)
      read(input,*) vari(2)
      call get_command_argument(5,input)
      read(input,*) vari(3)
      call get_command_argument(6,input)
      read(input,*) vari(4)
      call get_command_argument(7,input)
      read(input,*) vari(5)
      call get_command_argument(8,input)
      read(input,*) vari(6)
      call get_command_argument(9,input)
      read(input,*) snd

      if(vari(6).eq.2) days=threesixty()
      shmodel=model(1:index(model,'_')-1)

      valid=mx*(my+1)
      do mem=runs,runs
      smem=ensmem(mem)

      open(2,file=trim(tmp_path)//trim(model)//trim(sno(1))//
     & trim(smem)//trim(sno(2)),
     & form='unformatted',access='direct',recl=mx*(my+1)*4)
      open(3,file=trim(tmp_path)//trim(model)//trim(snoanom(1))//
     & trim(smem)//trim(snoanom(2)),
     & form='unformatted',access='direct',recl=mx*(my-1)*4)
      open(10,file=trim(in_path)//trim(snd),
     & form='unformatted',access='direct',recl=mx*(my+1)*4)

      id=0
      fid=0
C     if(model.eq.'NorESM1-M_hist') id=1
C     print*,trim(model)//trim(sno(1))//trim(smem)//trim(sno(2))

      do k=1,4
         open(16,file=trim(slp_std_path)//trim(era_slp_std(1))//
     & trim(seasons(k))//trim(era_slp_std(2)))
         read(16,'(10f9.3)') (erastd(i,k),i=1,valid)
         close(16)
         open(17,file=trim(slp_mean_path)//trim(era_slp_mean(1))//
     & trim(seasons(k))//trim(era_slp_mean(2)))
         read(17,'(10f9.3)') (eraave(i,k),i=1,valid)
         close(17)
         open(18,file=trim(mean_path)//trim(shmodel)//trim(snl(1))//
     & trim(snl(2))//trim(seasons(k))//trim(snl(3)))
         read(18,'(10f9.3)') (std(i,k),i=1,valid)
         close(18)
         open(19,file=trim(mean_path)//trim(shmodel)//trim(sne(1))//
     & trim(sne(2))//trim(seasons(k))//trim(sne(3)))
         read(19,'(10f9.3)') (ave(i,k),i=1,valid)
         close(19)
      enddo

      iyear=vari(2)
      imonth=vari(3)
      iday=vari(4)
      ihour=vari(5)-6
      leap=leapday(iyear,vari(6))

      do t=1,vari(1)
         id=id+1
         ihour=ihour+6
         if(ihour.gt.18) then
            ihour=0
            iday=iday+1
            if(leap.eq.1.and.imonth.eq.2) then
               if(iday.gt.days(imonth)+1) then
                  iday=1
                  imonth=imonth+1
               endif
            else
               if(iday.gt.days(imonth)) then
                  iday=1
                  imonth=imonth+1
                  if(imonth.gt.12) then
                     imonth=1
                     iyear=iyear+1
                     leap=leapday(iyear,vari(6))
                  endif
               endif
            endif
         endif

         read(10,rec=id)((p(i,j),i=1,mx),j=1,my+1)
         do i=1,mx
            do j=1,my+1
               if(p(i,j).lt.-899.99) then
                 p(i,j)=-900.99
               else
                 p(i,j)=p(i,j)
               endif
            enddo
         enddo

         do i=1,mx
            do j=1,my+1
               values((j-1)*mx+i)=p(i,j)
            enddo
         enddo

         if(imonth.lt.4) then
            k=1
         else if(imonth.lt.7) then
            k=2
         else if(imonth.lt.10) then
            k=3
         else
            k=4
         endif

         if(ihour.eq.0.or.ihour.eq.6.or.ihour.eq.12.or.ihour.eq.18) then
           fid=fid+1
           do i=1,valid
              if(values(i).lt.-899.99.or.ave(i,k).lt.-899.99.or.
     & std(i,k).lt.-899.99.or.erastd(i,k).lt.-899.99.or.
     & eraave(i,k).lt.-899.99) then
                 results(i)=-900.99
              else 
                 results(i)=(values(i)-ave(i,k))/std(i,k)*erastd(i,k)
     & +eraave(i,k)
              endif 
           enddo
           write(2,rec=fid) (results(i),i=1,valid)

           do i=1,valid
              if(values(i).lt.-899.99.or.ave(i,k).lt.-899.99.or.
     & std(i,k).lt.-899.99.or.erastd(i,k).lt.-899.99) then
                 results(i)=-900.99
              else
                 results(i)=(values(i)-ave(i,k))/std(i,k)*erastd(i,k)
              endif
           enddo
           do i=1,mx*(my-1)
              q(i)=results(i+mx)
           enddo
           write(3,rec=fid) (q(i),i=1,mx*(my-1))
         endif
      enddo

99    close(2)
      close(3)
      close(10)
      enddo

c      call idate(date)
c      call itime(time)
c      write(*,2000) date(2),date(1),date(3)+2000,time
c2000  format('slp_both_adjusted end at: Date ', i2.2, '/', i2.2,
c     &         '/', i4.4,'; time ',i2.2, ':', i2.2, ':', i2.2 )

      end
