c Code to calculate sea ice lead width distribution in lead numbers or
c lead number distribution (LND) from lead width distribution in lead
c areas (LWD)
      parameter(nx1=360,ny1=120,nx=nx1,ny=ny1,imt=nx1,jmt=ny1,nt=1)
      parameter(mx1=360,my1=120,mx=mx1,my=my1)
      parameter(m1=12,idelt=4,nregion=4,cutoff=0.10,npts=8)
      parameter(m2=2*6,m2p=m2+1,m2n=m2-1,m2p1=m2+2)
      parameter(m3=2*6,m3p=m3+1,m3n=m3-1,m3p1=m3+2)
      parameter(nyear1=1980,nyear2=2023,nmon=12)

      dimension heff(imt,jmt),h(m1),xn(nx1,ny1),yn(nx1,ny1)
      dimension heff1(imt,jmt), therm(imt,jmt), tair(imt,jmt)
      dimension heffm(imt,jmt,12),gice(imt,jmt,m1),aiday(imt,jmt)
      dimension tair1(imt,jmt),tairm(imt,jmt,365)
      dimension clon(imt,jmt),clat(imt,jmt),kmt(imt,jmt)
      DIMENSION imink(npts),jmink(npts),clat_k(npts), clon_k(npts)
      DIMENSION iyeark(npts),idayk(npts)
      dimension mask(nx1,ny1,nregion),count(nregion)
     &,count_overice(nregion), hmon(nregion), hyr(nregion)
     &,vmon(nregion), vyr(nregion)
     &,month(12),month_end(13)

      real gfsd(imt,jmt,m2),cfsd(imt,jmt,m2),gfnd(imt,jmt,m2)
     &,deff(imt,jmt),Lmesh(m2)

      real glwd(imt,jmt,m3),clwd(imt,jmt,m3),glnd(imt,jmt,m3)
     &,weff(imt,jmt),Wmesh(m3)

      dimension ulat(imt,jmt),ulon(imt,jmt),HTN(imt,jmt),HTE(imt,jmt)
     &,HUS(imt,jmt),HUW(imt,jmt),angle(imt,jmt),dxt(imt,jmt)
     &,dyt(imt,jmt)

      character *80 fopen(5), f1,f2,f3
      character *4 cyear(1900:2100),cyear1(1900:2100)
      character *12 char
      integer SLEN

      data month/31,28,31,30,31,30,31,31,30,31,30,31/
      data month_end/0,31,59,90,120,151,181,212,243,273,304,334,365/

      f1='glwd.H'
      f2='area.H'
      f3='glnd.H'

c for ice thickness distribution
      H(1)=0.0E+00
      if(m1.eq.12) then
c M1=2*6
        C1=0.20E+00 
        C2=110.0E+00
        C2=c2*1.
        C3=7.0E+00
      else if(m1.eq.8) then
c M1=2*4
        C1=0.20E+00
        C2=110.0E+00
        C2=c2*1.
        C3=14.2E+00
      end if

      DO 90 M=2,M1
      OM=FLOAT(M-1)
      OMS=OM**2
      H(M)=H(M-1)+C1+C3*(1.0E+00-EXP(-OMS/C2))
 90   CONTINUE

c M2=2*6, for floe size distribution
      IF(m2.EQ.12) THEN
        GC1=0.2
        GC2=110.0*100.0
        GC2=GC2*1.0
        GC3=7.0*7792.0
      END IF

      Lmesh(1)=0.0
      DO 91 M=2,m2
      OM=M-1
      OMS=OM**2
      Lmesh(M)=Lmesh(M-1)+GC1+GC3*(1.0-EXP(-OMS/GC2))
 91   CONTINUE

c M3=2*6, for lead wisth distribution
      IF(m3.EQ.12) THEN
        GC1=0.2
        GC2=414.0
        GC2=gc2*1.0
        GC3=7.0
      END IF

      Wmesh(1)=0.0
      DO 92 M=2,m3
      OM=M-1
      OMS=OM**2
      Wmesh(M)=Wmesh(M-1)+GC1+GC3*(1.0-EXP(-OMS/GC2))
 92   CONTINUE

      open(20,file='grid.dat.rot')
      read(20,'(10f8.2)') ((clon(i,j),i=1,nx1),j=1,ny1)
      read(20,'(10f8.2)') ((clat(i,j),i=1,nx1),j=1,ny1)
      close(20)

      open(24,file='grid.dat.pop')
        read(24,'(10f8.2)') ((ulat(i,j),i=1,nx),j=1,ny)
        read(24,'(10f8.2)') ((ulon(i,j),i=1,nx),j=1,ny)
        read(24,'(10f8.2)') ((HTN  (i,j),i=1,nx),j=1,ny)
        read(24,'(10f8.2)') ((HTE  (i,j),i=1,nx),j=1,ny)
        read(24,'(10f8.2)') ((HUS  (i,j),i=1,nx),j=1,ny)
        read(24,'(10f8.2)') ((HUW  (i,j),i=1,nx),j=1,ny)
        read(24,'(10f8.2)') ((angle(i,j),i=1,nx),j=1,ny)
      close(24)

      do i=2,nx1-1
      do j=2,ny1-1
        dxt(i,j)=0.5*(HTN(i,j)+HTN(i,j-1))  ! in km
        dyt(i,j)=0.5*(HTE(i,j)+HTE(i-1,j))  ! in km
      end do
      end do

      do j=2,ny1-1
        dxt(1,j)=dxt(2,j)  ! in km
        dyt(1,j)=dyt(2,j)  ! in km
        dxt(nx1,j)=dxt(nx1-1,j)  ! in km
        dyt(nx1,j)=dyt(nx1-1,j)  ! in km
      end do

      do i=1,nx1
        dxt(i,1)=dxt(i,2)  ! in km
        dyt(i,1)=dyt(i,2)  ! in km
        dxt(i,ny1)=dxt(i,ny1-1)  ! in km
        dyt(i,ny1)=dyt(i,ny1-1)  ! in km
      end do

        pi=3.1416
        radius=6370.
        do 8 i=1,nx1
        do 8 j=1,ny1
        rr=2*radius*sind((45-0.5*clat(i,j)))
        xn(i,j)=rr*cosd((clon(i,j)-35))
        yn(i,j)=rr*sind((clon(i,j)-35))
8       continue

      open(20,file=
     &'io.dat_360_120.output',form='formatted')
      read(20,'(360i2)') kmt
c      read(20) kmt
      close(20)

c get model sub-regions
      do i=1,nx1
      do j=1,ny1
      do k=1,nregion
      mask(i,j,k)=0
      end do
      end do
      end do

c mask(3) is arctic basin + barents sea
      do i=99+1,265+1
      do j=33+1+(ny1-my1),ny1-3+1
      if(kmt(i,j).gt.0) mask(i,j,3)=1
      end do
      end do
      do i=96+1,98+1
      do j=ny1-23,2,-1
      if(kmt(i,j).gt.0) mask(i,j,3)=1
      end do
      end do

      k=0
      jk=0
      do j=ny1-29,ny1-29-25,-1
        k=k+1
        if (mod(k,2) .eq. 1) jk=jk+1
        do i=92-jk+1, 92-jk+50+1
          if(kmt(i,j).gt.0) mask(i,j,3)=1
        end do
      end do

      do i=129-50+1,130+1
      do j=1+1,64+1+(ny1-my1)
        if(kmt(i,j).gt.0) mask(i,j,3)=1
      end do
      end do

c mask(1) is GIN Sea + South of GIN Sea
      do i=2,nx1/2
      do j=2,ny1-1
      if(kmt(i,j).gt.0.and.mask(i,j,3).eq.0) then
      mask(i,j,1)=1
      end if
      end do
      end do

c mask(2) is Bering Sea
      do i=nx1/2-30,nx1-80
      do j=2,33+1
      if(kmt(i,j).gt.0.and.mask(i,j,3).eq.0) then
      mask(i,j,2)=1
      end if
      end do
      end do

c mask(4) for globe ocean
      do i=2,nx1-1
      do j=2,ny1-1
      if(kmt(i,j).gt.0) then
      mask(i,j,4)=1
      end if
      end do
      end do

      do 999 iyear=nyear1,nyear2

c      write(*,*) iyear

      write(unit=cyear(iyear),fmt='(i4)') iyear
      i=slen(f1)
      open(1,file=f1(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')

      i=slen(f2)
      open(2,file=f2(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='old')

      i=slen(f3)
      open(12,file=f3(1:i)//cyear(iyear)
     &,access='direct',form='unformatted',recl=nx1*ny1*4
     &,status='unknown')

      do i=1,nregion
      hyr(i)=0.0
      vyr(i)=0.0
      end do

      iday=0
      do i=1,imt
      do j=1,jmt
      tair1(i,j)=0.0
      end do
      end do

      do imon=1,nmon

c calculate mean lead width
      weff=0.0
      do m=1,m3
      read(1) ((glwd(i,j,m),i=1,NX1),j=1,NY1)
      weff=weff+glwd(:,:,m)*Wmesh(m)
      end do

      read(2) ((aiday(i,j),i=1,NX1),j=1,NY1)

c calculate lead width distribution in lead numbers per area of 10km x 10km

      cell_area=10000.0*10000.0 ! 10km x 10km
      cell_length=10000.0

      glnd=0.0
      DO J=1,jmt
      DO I=1,imt
        if(aiday(i,j).gt.0.15) then
          do m=2,m3
           glnd(i,j,m)=glwd(i,j,m)*cell_area/(cell_length*Wmesh(m))
          end do
        end if
      end do
      end do

      do m=1,m3
c output of lead width distribution in lead numbers per 10km x 10km area
        write(12) ((glnd(i,j,m),i=1,NX1),j=1,NY1)
      end do

c below is to calculate mean lead width over model sub-regions
      tair=weff
      tairm(:,:,imon)=tair

      do i=1,nregion
      hmon(i)=0.0
      vmon(i)=0.0
      count_overice(i)=0.0
      end do

      do i=2,imt-1
      do j=2,jmt-1
      do ia=1,nregion
      if(mask(i,j,ia).eq.1) then
      vmon(ia)=vmon(ia)+tairm(i,j,imon)*dxt(i,j)*dyt(i,j)
      end if
      if(mask(i,j,ia).eq.1) then
      count_overice(ia)=count_overice(ia)+dxt(i,j)*dyt(i,j)
      hmon(ia)=hmon(ia)+tairm(i,j,imon)*dxt(i,j)*dyt(i,j)
      end if
      end do
      end do
      end do

      do ia=1,nregion
      count_overice(ia)=max(count_overice(ia),1.0)
      end do

      do ia=1,nregion
      hmon(ia)=hmon(ia)/count_overice(ia)
      vmon(ia)=vmon(ia)*1000.0*1000.0*1.0E-12  ! 10**12 m**3
      hyr(ia)=hyr(ia)+hmon(ia)
      vyr(ia)=vyr(ia)+vmon(ia)
      end do

      end do

      do ia=1,nregion
      hyr(ia)=hyr(ia)/12.0
      vyr(ia)=vyr(ia)/12.0
      end do

      close(1)
      close(2)
      close(12)
999   continue

      close(31)
      close(32)
      
      stop  
      end

      INTEGER FUNCTION slen (string)
C ---
C --- this function computes the length of a character string less
C --- trailing blanks
C --- slen > 0, length of string less trailing blanks
C ---      = 0, character string is blank
C ---
      CHARACTER*(*) string
      CHARACTER*1 cblank
      INTEGER i
      DATA cblank/' '/
C ---
      DO 50 i = LEN(string), 1, -1
         IF (string(i:i) .NE. ' ')  GO TO 100
50    CONTINUE
      i = 0
100   CONTINUE
      slen = i
      RETURN
      END
