c     gt_get_edge_missing.f
c
c     Get Gt from adjusted SLP with missing lat -90, lat 90
c
      use commondata, only :  leapday,ensmem,threesixty,mx,my,pi,days,
     & seasons,latlon,tmp_path,sur_path,date,time

      implicit none

      real :: p(200,100),q(200,100),dp2(200,100)
      real :: lat(200,100),lon(200,100),mis(200,100)
      real :: mp,phi,dpx1,dpx2,dpy1,dpy2,dp
      integer :: vari(6),fid,mem,runs,leap
      integer :: iyear,imonth,iday,ihour,id,t,i,j
      character(len=50) :: snd(2),sno(3),smem,input,model

      data snd /'_r','_SLP_Adjusted.dat'/
      data sno /'_r','_Gt_','.txt'/

c      call idate(date)
c      call itime(time)
c      write(*,1000) date(2),date(1),date(3)+2000,time
c1000  format('gt_get start at: Date ', i2.2, '/', i2.2, '/', i4.4,
c     &         '; 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)

      if(vari(6).eq.2) days=threesixty()

      open(2,file=trim(sur_path)//trim(latlon))
      do j=1,my+1
         do i=1,mx
            read(2,*) mp,lat(i,j),lon(i,j)
         enddo
         lon(mx+1,j)=lon(1,j)
         lat(mx+1,j)=lat(1,j)
      enddo
      close(2)

      do mem=runs,runs
      smem=ensmem(mem)

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

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

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

      open(21,file=trim(tmp_path)//trim(model)//trim(sno(1))//
     & trim(smem)//trim(sno(2))//trim(seasons(1))//trim(sno(3)))
      open(22,file=trim(tmp_path)//trim(model)//trim(sno(1))//
     & trim(smem)//trim(sno(2))//trim(seasons(2))//trim(sno(3)))
      open(23,file=trim(tmp_path)//trim(model)//trim(sno(1))//
     & trim(smem)//trim(sno(2))//trim(seasons(3))//trim(sno(3)))
      open(24,file=trim(tmp_path)//trim(model)//trim(sno(1))//
     & trim(smem)//trim(sno(2))//trim(seasons(4))//trim(sno(3)))

      do t=1,vari(1)/2
         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(3,rec=id) ((p(i,j),i=1,mx),j=1,my+1)
         do j=1,my+1
            p(mx+1,j)=p(1,j)
         enddo
         mis=0.0
         phi=pi*lat(1,2)/180.
         dpx1=(p(1,2)-p(mx,2))/cos(phi)
         dpx2=(p(2,2)-p(1,2))/cos(phi)
         dpy1=p(1,3)-p(1,2)
         dpy2=p(1,3)-p(1,2)
         dp=(dpx1**2+dpx2**2)/2.+(dpy1**2+dpy2**2)/2.
         dp2(1,2)=dp
         do j=3,my-1
            phi=pi*lat(1,j)/180.
            dpx1=(p(1,j)-p(mx,j))/cos(phi)
            dpx2=(p(2,j)-p(1,j))/cos(phi)
            dpy1=p(1,j)-p(1,j-1)
            dpy2=p(1,j+1)-p(1,j)
            dp=(dpx1**2+dpx2**2)/2.+(dpy1**2+dpy2**2)/2.
            dp2(1,j)=dp
            if(p(1,j).lt.-899.99.or.p(mx,j).lt.-899.99.or.
     &      p(2,j).lt.-899.99.or.p(1,j-1).lt.-899.99.or.
     &      p(1,j+1).lt.-899.99) then
               mis(1,j)=-900.99
            endif
         enddo
         phi=pi*lat(1,my)/180.
         dpx1=(p(1,my)-p(mx,my))/cos(phi)
         dpx2=(p(2,my)-p(1,my))/cos(phi)
         dpy1=p(1,my)-p(1,my-1)
         dpy2=p(1,my)-p(1,my-1)
         dp=(dpx1**2+dpx2**2)/2.+(dpy1**2+dpy2**2)/2.
         dp2(1,my)=dp
         do i=2,mx
            phi=pi*lat(i,2)/180.
            dpx1=(p(i,2)-p(i-1,2))/cos(phi)
            dpx2=(p(i+1,2)-p(i,2))/cos(phi)
            dpy1=p(i,3)-p(i,2)
            dpy2=p(i,3)-p(i,2)
            dp=(dpx1**2+dpx2**2)/2.+(dpy1**2+dpy2**2)/2.
            dp2(i,2)=dp
            do j=3,my-1
               phi=pi*lat(i,j)/180.
               dpx1=(p(i,j)-p(i-1,j))/cos(phi)
               dpx2=(p(i+1,j)-p(i,j))/cos(phi)
               dpy1=p(i,j)-p(i,j-1)
               dpy2=p(i,j+1)-p(i,j)
               dp=(dpx1**2+dpx2**2)/2.+(dpy1**2+dpy2**2)/2.
               dp2(i,j)=dp
               if(p(i,j).lt.-899.99.or.p(i-1,j).lt.-899.99.or.
     &         p(i+1,j).lt.-899.99.or.p(i,j-1).lt.-899.99.or.
     &         p(i,j+1).lt.-899.99) then
                  mis(i,j)=-900.99
               endif
            enddo
            phi=pi*lat(i,my)/180.
            dpx1=(p(i,my)-p(i-1,my))/cos(phi)
            dpx2=(p(i+1,my)-p(i,my))/cos(phi)
            dpy1=p(i,my)-p(i,my-1)
            dpy2=p(i,my)-p(i,my-1)
            dp=(dpx1**2+dpx2**2)/2.+(dpy1**2+dpy2**2)/2.
            dp2(i,my)=dp
         enddo
         do i=1,mx
            do j=1,my-1
               if(mis(i,j+1).lt.-899.99) then
c                  print*,iyear,imonth,iday,ihour,q(i,j)
                  q(i,j)=-900.99
               else
                  q(i,j)=dp2(i,j+1)
               endif
            enddo
         enddo

         if(imonth.lt.4) then
            fid=21
         else if(imonth.lt.7) then
            fid=22
         else if(imonth.lt.10) then
            fid=23
         else
            fid=24
         endif
 
         write(fid,'(i6,3i4.2)') iyear,imonth,iday,ihour
         write(fid,'(10f9.3)') ((q(i,j),i=1,mx),j=1,my-1)
      enddo

      close(3)
      close(21)
      close(22)
      close(23)
      close(24)
      enddo

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

      end
