SUBROUTINE MESH90
!use consts
use global
implicit none  
    
real(DP) :: zzmin,zzmax,jtzi,jbzi,dzloc,indz,jzz,jzsupp,dzcor,nzloc2,dz,dzzone,dzmesh
real(DP) :: nzsupp,dz1,temp
integer :: iz,i,j,nzleft,ind0,nzloci
save
temp = 0_DP
CALL MESGER90(MAXRLO,IDIM,R,2,IMUT,temp,RADIUS,DRMIN,FAC,NRLOC,NSERIE,IRLOC)

      IRR=IMUT-1

CALL MESGEN90(2,IDIM,R,IMUT,IMUT,RADIUS,temp,DRRMIN,DRMIN,DRMAX,FAC,FACL)

      IMAX=IMUT-1
      IMAX1=IMAX+1
!-----DP OCT 1996
      IF (IMAX1.GT.IDIM) THEN
        WRITE(LUW1,400) IMAX1,IDIM
        RETURN
      END IF
  400 FORMAT(' ***ERROR***IN TYPE 60',/, &
     ' THE NUMBER OF RADIAL GLOBAL MESHES (IMAX1=',I3,')',/, &
     ' EXCEEDS THE ARRAY DIMENSION FIXED TO IDIM=',I3,'.',/, &
     ' IF YOU STILL WANT TO RUN YOUR CASE, INCREASE THE VALUE OF',/, &
     ' "IDIM" IN THE "TRNVDST.INC" INCLUDE FILE AND RECOMPILE.')
!-----END DP OCT 1996
      DEEP=DEPTH+HEIGHT
!-----IF(DEPTH.EQ.0.) GO TO 10
      IF(DEPTH.EQ.0.) GO TO 8
      IF(DEPTH.GT.ZDMIN) GO TO 6
      temp = 0
      CALL MESGEN90(1,JDIM,Z,2,IMUT,temp,DEPTH,DRRMIN,DRMIN,DRMAX, &
      FAC,FACL)

      GO TO 8
       
    6 CALL MESGEN90(3,JDIM,Z,2,IMUT,DEPTH,temp,DRRMIN,DRMIN,DRMAX, &
      FAC,FACL)

    8 CONTINUE
!-----DP OCT 1996
      IF(DEPTH.EQ.0.) IMUT=2
      JTZ=IMUT
      ZZMIN=DEPTH
      DO 10  IZ=1,NLAYST
        ZZMAX=ZZMIN+THLST(IZ)
        JTZS(IZ)=IMUT
        CALL MESGEN90(1,JDIM,Z,IMUT,IMUT,ZZMIN,ZZMAX,DRRMIN,DRMIN, &
        DRMAX,FAC,FACL)

        JBZS(IZ)=IMUT-1
        ZZMIN=ZZMAX
   10 CONTINUE
!-----END DP OCT 1996
      JBZ=IMUT-1
      CALL MESGEN90(2,JDIM,Z,IMUT,IMUT,DEEP,temp,DRRMIN,DRMIN,DRMAX, &
      FAC,FACL)


      JMAX=IMUT-1
      JMAX1=JMAX+1
!-----DP OCT 1996
      IF (JMAX1.GT.JDIM) THEN
        WRITE(LUW1,401) JMAX1,JDIM
        RETURN
      END IF
  401 FORMAT(' ***ERROR***IN TYPE 60',/, &
     ' THE NUMBER OF VERTICAL GLOBAL MESHES (JMAX1=',I3,')',/, &
     ' EXCEEDS THE ARRAY DIMENSION FIXED TO JDIM=',I3,'.',/, &
     ' IF YOU STILL WANT TO RUN YOUR CASE, INCREASE THE VALUE OF',/, &
     ' "JDIM" IN THE "TRNVDST.INC" INCLUDE FILE AND RECOMPILE.')
!-----END DP OCT 1996
      CALL MESGEN90(4,IDIMR,RR,2,IMUT,R0,R1,DRRMIN,DRMIN,DRMAX, &
      FAC,FACL)

      IMAXR=IMUT-1
      IMAXR1=IMAXR+1
!-----DP OCT 1996
      IF (IMAXR1.GT.IDIMR) THEN
        WRITE(LUW1,402) IMAXR1,IDIMR
        RETURN
      END IF
  402 FORMAT(' ***ERROR***IN TYPE 60',/, &
     ' THE NUMBER OF RADIAL LOCAL MESHES (IMAXR1=',I3,')',/, &
     ' EXCEEDS THE ARRAY DIMENSION FIXED TO IDIMR=',I3,'.',/, &
     ' IF YOU STILL WANT TO RUN YOUR CASE, INCREASE THE VALUE OF',/, &
     ' "IDIMR" IN THE "TRNVDST.INC" INCLUDE FILE AND RECOMPILE.')
!-----END DP OCT 1996
      DO 30 I=2,IMAX
   30 RM(I)=.5*(R(I)+R(I+1))
      RM(1)=R(1)
      RM(IMAX1)=R(IMAX1)
      DO 40 J=2,JMAX
   40 ZM(J)=.5*(Z(J)+Z(J+1))
      ZM(1)=Z(1)
      ZM(JMAX1)=Z(JMAX1)
      DO 50 I=2,IMAXR
   50 RMR(I)=.5*(RR(I)+RR(I+1))
      RMR(1)=RR(1)
      RMR(IMAXR1)=RR(IMAXR1)
      IRR1=IRR+1
      JBZ1=JBZ+1
!
!-Determine the vertical subregions
!-Initialise the array indice JZLoc(i)
      DO 60 I=1,MAXZLO
        JZLOC(I)=0
   60 CONTINUE
!-The number of subregions NZLOC can not exceed the number of meshes in
!-the store volume (=JBZ+1-JTZ), as the vertical division of the store
!-is kept unchanged for any NZLOC.
      IF (NZLOC.GT.(JBZ+1-JTZ)) NZLOC=(JBZ+1-JTZ)  
!-DP OCT 1996
!-CHECK THE VALUE OF NZLOC
      IF (NZLOC.GT.MZLOM1) THEN
        WRITE(LUW1,403) NZLOC,MZLOM1
        RETURN
      END IF
  403 FORMAT(' ***ERROR***IN TYPE 60',/, &
      ' THE NUMBER OF VERTICAL SUBREGIONS (NZLOC=',I3,')',/, &
      ' EXCEEDS THE MAXIMUM ALLOWED (MZLOM1=',I3,').',/, &
      ' IF YOU STILL WANT TO RUN YOUR CASE, INCREASE THE VALUE OF',/, &
      ' "MAXZLO" IN THE "TRNVDST.INC" INCLUDE FILE AND RECOMPILE.')
!-Distribution of the NZLOC subregions in the NLAYST ground layers
      DO 61 IZ=1,NLAYST
        NZLOCS(IZ)=1
   61 CONTINUE
      NZLEFT=NZLOC-NLAYST
      IZ=1
   62 CONTINUE
        IF ((NZLOCS(IZ).LT.(JBZS(IZ)+1-JTZS(IZ))).AND. &
          (NZLEFT.GE.1)) THEN
             NZLOCS(IZ)=NZLOCS(IZ)+1
             NZLEFT=NZLEFT-1
        END IF
!-------END IF ( .LT. ).AND.( .GE. )
        IZ=IZ+1
        IZ=MOD(IZ-1,NLAYST)+1
      IF (NZLEFT.GE.1) GOTO 62
      IND0=1
      DO 100 IZ=1,NLAYST
        NZLOCI= NZLOCS(IZ)
        JTZI  = JTZS(IZ)
        JBZI  = JBZS(IZ)
!-Recommended vertical extension of each vertical subregion.
        DZLOC=THLST(IZ)/NZLOCI
!-Number of meshes that can be put together.
        NZLEFT=(JBZI+1-JTZI)-NZLOCI
!-Distribution of the meshes in the different vertical subregions.
        JZLOC(IND0)=JTZI
        JZLOC(NZLOCI+IND0)=JBZI+1
        INDZ=JTZI
        JZZ=1
        JZSUPP=0
        DZCOR=0.
        NZLOC2=NZLOCI/2
        DO 90 I=1,NZLOC2
   70     CONTINUE
            DZ=Z(INDZ+JZZ)-Z(INDZ)
            IF (((DZLOC+DZCOR).LE.DZ ).OR.(NZLEFT.LT.2)) GO TO 80
            JZZ=JZZ+1
            NZLEFT=NZLEFT-2
            GO TO 70
   80     CONTINUE
          DZZONE= (DZLOC+DZCOR) - ( Z(INDZ+JZZ-1) - Z(INDZ) )
          DZMESH= Z(INDZ+JZZ)  -  Z(INDZ+JZZ-1)
          IF ( (JZZ.GT.1) .AND. (DZZONE.LT.(0.75*DZMESH)) ) THEN
            JZZ=JZZ-1
            NZLEFT=NZLEFT+2
          END IF
!-When NZLEFT=1, the number of meshes in the lower subregion is increased
!-by 1 if its vertical extension remain below 1.25*(DZLOC+DZCOR).
          NZSUPP=0
          IF (NZLEFT.EQ.1) THEN
            DZ1= Z(INDZ+JZZ+1) - Z(INDZ)
            IF ( DZ1 .LT. (1.25*(DZLOC+DZCOR)) ) NZSUPP=1
          ENDIF
          IF (NZSUPP.EQ.1)  THEN
            JZSUPP=1
            NZLEFT=NZLEFT-NZSUPP
          ENDIF
          JZLOC(NZLOCI+IND0-I)= JBZI+1 - (INDZ-JTZI) - JZZ - JZSUPP
          JZLOC(IND0+I)= INDZ + JZZ
          DZCOR= I*DZLOC - (Z(INDZ+JZZ)-Z(JTZI))
          INDZ= INDZ+JZZ
          JZZ= 1
   90   CONTINUE
        IND0=IND0+NZLOCI
  100 CONTINUE
  
RETURN
END subroutine mesh90