c Code to calculate sea ice mean lead width using lead width
c distribution in lead areas (LWD)
      parameter (imt=  360, jmt= 120, km= 30)
      parameter (nx1=imt, ny1=jmt, nx=imt,ny=jmt)
      parameter(m1=2*6,m1p=m1+1,m1n=m1-1,m1p1=m1+2)
      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(nregion=4,cutoff=0.15)
      parameter(nyear1=1980,nyear2=2023,nmon1=1,nmon2=12)

      dimension heff(imt,jmt),h(m1),xn(nx1,ny1),yn(nx1,ny1),hi(imt,jmt)
      dimension clon(imt,jmt),clat(imt,jmt),kmt(imt,jmt)
      dimension mask(nx1,ny1,nregion),count(nregion)
     &,count_overice(nregion), hmon(nregion), hyr(nregion)
     &,vmon(nregion),vyr(nregion),tmon(nregion),tyr(nregion)
     &,KY1(0:12)

      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)

      dimension temp(imt,jmt,12), temp1(imt,jmt), work(imt,jmt)
      real Lmesh(m2),gfsd(imt,jmt,m2),deff(imt,jmt)

      real Wmesh(m3),glwd(imt,jmt,m3),weff(imt,jmt)

      DIMENSION wx(nx1,ny1),wy(nx1,ny1),area(nx1,ny1,4)
     &,avg_n(4,2,12),avg_s(4,2,12)
     &,dxtice(imt),dytice(jmt),g(nx1,ny1,m1)
      character *80 fopen(5), f1,f2,f12,f22
      character *4 cyear(1900:2100),cyear1(1900:2100)
      character *12 char
      integer SLEN
      DATA KY1/0,31,59,90,120,151,181,212,243,273,304,334,365/

      f2= 'glwd.H'
      f22='weff.H'

      open(3,file=
     &'z')
      open(4,file=
     &'z1')

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
        C3=20.0E+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 width distribution
      IF(m3.EQ.12) THEN
        GC1=0.2
        GC2=110.0*100.0
        GC2=GC2*1.0
        GC3=7.0*7600.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)

      c0=0.0
      p5=0.5
c DXT and DYT are in km
c      call s_shift(WORK,HTN)
      WORK=eoshift(HTN,dim=2,shift=-1)
      DXT = p5*(HTN + WORK)
c      call w_shift(WORK,HTE)
      WORK=eoshift(HTE,dim=1,shift=-1)

      DYT = p5*(HTE + WORK)
c      where (DXT == c0) DXT=c1
c      where (DYT == c0) DYT=c1

        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 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

!mask(3) is arctic basin + barents sea
      do i=99+1,265+1
      do j=33+1,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
        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 whole model domain or Arctic circle
      do i=2,nx1-1
      do j=2,ny1-1
c      if(kmt(i,j).gt.0) then ! whole model domain
      if(kmt(i,j).gt.0.and.clat(i,j).gt.66.5) then ! arctic circle
      mask(i,j,4)=1
      end if
      end do
      end do

      do 999 iyear=nyear1,nyear2

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

      do imon=1,12
      do i=1,imt
      do j=1,jmt
      temp(i,j,imon)=0.0
      end do
      end do
      end do

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

      do imon=nmon1,nmon2

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

c below is to calculate mean lead width in model sub-regions
      do i=1,imt
      do j=1,jmt
        heff(i,j)=weff(i,j)
      end do
      end do

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

      do i=1,imt-0
      do j=2,jmt-1
      do ia=1,nregion
      if(mask(i,j,ia).eq.1) then
        vmon(ia)=vmon(ia)+heff(i,j)*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)+heff(i,j)*dxt(i,j)*dyt(i,j)
      tmon(ia)=tmon(ia)+temp(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*2
      tmon(ia)=tmon(ia)/count_overice(ia)
      hyr(ia)=hyr(ia)+hmon(ia)
      vyr(ia)=vyr(ia)+vmon(ia)
      tyr(ia)=tyr(ia)+tmon(ia)
      end do

      write(3,'(2i5,8f8.1)')iyear,imon,hmon(1),hmon(2),hmon(3),hmon(4)
     &,                                vmon(1),vmon(2),vmon(3),vmon(4)

      end do

      do ia=1,nregion
      hyr(ia)=hyr(ia)/(nmon2-nmon1+1.0)
      vyr(ia)=vyr(ia)/(nmon2-nmon1+1.0)
      end do

      write(4,'(i5,8f8.1)')iyear,hyr(1),hyr(2),hyr(3),hyr(4)
     &,                          vyr(1),vyr(2),vyr(3),vyr(4)
      write(*,'(i5,8f8.1)')iyear,hyr(1),hyr(2),hyr(3),hyr(4)
     &,                          vyr(1),vyr(2),vyr(3),vyr(4)

      close(2)
      close(22)
999   continue
      
      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
