!
!
!
SUBROUTINE thermcell_flux(ngrid,nlay,ptimestep,masse,                         &
                          lalim,lmin,lmax,alim_star,entr_star,detr_star,      &
                          f,rhobarz,zlev,zw2,fm,entr,                         &
                          detr,zqla,lev_out,lunout1,igout)
      
      
!==============================================================================
! thermcell_flux: deduction des flux
!==============================================================================
      
      USE print_control_mod, ONLY: prt_level
      USE thermcell_mod
      
      IMPLICIT NONE
      
      
!==============================================================================
! Declaration
!==============================================================================
      
!      inputs:
!      ------- 
      
      INTEGER ngrid
      INTEGER nlay
      INTEGER igout
      INTEGER lev_out
      INTEGER lunout1
      INTEGER lmin(ngrid)
      INTEGER lmax(ngrid)
      INTEGER lalim(ngrid)
      
      REAL alim_star(ngrid,nlay)
      REAL entr_star(ngrid,nlay)
      REAL detr_star(ngrid,nlay)
      REAL zw2(ngrid,nlay+1)
      REAL zlev(ngrid,nlay+1)
      REAL masse(ngrid,nlay)
      REAL ptimestep
      REAL rhobarz(ngrid,nlay)
      REAL f(ngrid)
      REAL zqla(ngrid,nlay)
      REAL zmax(ngrid)
      
!      outputs:
!      --------      
      
      REAL entr(ngrid,nlay)
      REAL detr(ngrid,nlay)
      REAL fm(ngrid,nlay+1)
      
!      local:
!      ------
      
      INTEGER ig,l,k
      INTEGER lout
      
      REAL zfm                            ! mass flux such as updraft fraction is equal to its maximal authorized value alphamax
      REAL f_old                          ! save initial value of mass flux
      REAL eee0                           ! save initial value of entrainment
      REAL ddd0                           ! save initial value of detrainment
      REAL eee                            ! eee0 - layer mass * maximal authorized mass fraction / dt
      REAL ddd                            ! ddd0 - eee
      REAL zzz                            ! temporary variable set to mass flux
      REAL fact
      REAL test
      
      INTEGER ncorecfm1
      INTEGER ncorecfm2
      INTEGER ncorecfm3
      INTEGER ncorecfm4
      INTEGER ncorecfm5
      INTEGER ncorecfm6
      INTEGER ncorecalpha
      
      LOGICAL check_debug
      LOGICAL labort_physic
      
      CHARACTER (len=20) :: modname='thermcell_flux'
      CHARACTER (len=80) :: abort_message
      
!==============================================================================
! Initialization
!==============================================================================
      
      ncorecfm1 = 0
      ncorecfm2 = 0
      ncorecfm3 = 0
      ncorecfm4 = 0
      ncorecfm5 = 0
      ncorecfm6 = 0
      
      ncorecalpha = 0
      
      entr(:,:) = 0.
      detr(:,:) = 0.
      fm(:,:)   = 0.
      
      IF (prt_level.ge.10) THEN
         write(lunout1,*) 'Dans thermcell_flux 0'
         write(lunout1,*) 'flux base ', f(igout)
         write(lunout1,*) 'lmax ', lmax(igout)
         write(lunout1,*) 'lalim ', lalim(igout)
         write(lunout1,*) 'ig= ', igout
         write(lunout1,*) ' l E*    A*     D*  '
         write(lunout1,'(i4,3e15.5)') (l,entr_star(igout,l),                  &
         &                             alim_star(igout,l),detr_star(igout,l), &
         &                             l=1,lmax(igout))
      ENDIF
      
!==============================================================================
! Calcul de l'entrainement, du detrainement et du flux de masse
!==============================================================================
      
!------------------------------------------------------------------------------
! Multiplication par la norme issue de la relation de fermeture
!------------------------------------------------------------------------------
      
      DO l=1,nlay
         entr(:,l) = f(:) * (entr_star(:,l) + alim_star(:,l))
         detr(:,l) = f(:) * detr_star(:,l)
      ENDDO
      
!------------------------------------------------------------------------------
! Calcul du flux de masse
!------------------------------------------------------------------------------
      
      DO l=1,nlay
         DO ig=1,ngrid
            IF (l.lt.lmax(ig) .and. l.ge.lmin(ig)) THEN
               fm(ig,l+1) = fm(ig,l) + entr(ig,l) - detr(ig,l)
            ELSEIF (l.eq.lmax(ig)) THEN
               fm(ig,l+1) = 0.
               detr(ig,l) = fm(ig,l) + entr(ig,l)
            ELSE
               fm(ig,l+1) = 0.
               entr(ig,l) = 0.
               detr(ig,l) = 0.
            ENDIF
         ENDDO
      ENDDO
      
!==============================================================================
! Tests de validite
!==============================================================================
      
      DO l=1,nlay
         
!------------------------------------------------------------------------------
! Verification de la positivite du flux de masse entrant
!------------------------------------------------------------------------------
         
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : I remove the correction and replace it by an uncompromising test.
!      According to the previous derivations, we MUST have positive mass flux
!      everywhere! Indeed, as soon as fm becomes negative, the plume stops.
!      Then the only value which can be negative is the mass flux at the top of
!      the plume. However, this value was set to zero a few lines above.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         
         labort_physic=.false.
         
         DO ig=1,ngrid
            IF (fm(ig,l).lt.0.) THEN
               labort_physic = .true.
               igout = ig
               lout = l
            ENDIF
         ENDDO
         
         IF (labort_physic) THEN
            print *, 'ERROR: mass flux has negative value(s)!'
            print *, 'ig,l,fm',igout, lout, fm(igout,lout)
            print *, 'lmin,lmax', lmin(igout), lmax(igout)
            print *, '-------------------------------'
            print *, 'fm+', fm(igout,lout+1)
            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
            print *, 'fm ', fm(igout,lout)
            print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
            print *, 'fm-', fm(igout,lout-1)
            abort_message = 'Negative incoming fm in thermcell_flux'
            CALL abort_physic(modname,abort_message,1)
         ENDIF
         
!------------------------------------------------------------------------------
! Test entrainment positif
!------------------------------------------------------------------------------
         
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : According to the previous derivations, we MUST have positive entrainment
!      and detrainment everywhere!
!      Indeed, they are set to max between zero and a computed value.
!      Then tests are uncompromising.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         
         labort_physic=.false.
         
         DO ig=1,ngrid
            IF (entr(ig,l).lt.0.) THEN
               labort_physic = .true.
               igout = ig
               lout = l
            ENDIF
         ENDDO
         
         IF (labort_physic) THEN
            print *, 'ERROR: entrainment has negative value(s)!'
            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
            print *, 'lmin,lmax', lmin(igout), lmax(igout)
            print *, '-------------------------------'
            print *, 'fm+', fm(igout,lout+1)
            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
            print *, 'fm ', fm(igout,lout)
            print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
            print *, 'fm-', fm(igout,lout-1)
            abort_message = 'Negative entrainment in thermcell_flux'
            CALL abort_physic(modname,abort_message,1)
         ENDIF
         
!------------------------------------------------------------------------------
! Test detrainment positif
!------------------------------------------------------------------------------
         
         DO ig=1,ngrid
            IF (detr(ig,l).lt.0.) THEN
               labort_physic = .true.
               igout = ig
               lout = l
            ENDIF
         ENDDO
         
         IF (labort_physic) THEN
            print *, 'ERROR: detrainment has negative value(s)!'
            print *, 'ig,l,detr', igout, lout, detr(igout,lout)
            print *, 'lmin,lmax', lmin(igout), lmax(igout)
            print *, '-------------------------------'
            print *, 'fm+', fm(igout,lout+1)
            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
            print *, 'fm ', fm(igout,lout)
            print *, 'entr,detr', entr(igout,lout-1), detr(igout,lout-1)
            print *, 'fm-', fm(igout,lout-1)
            abort_message = 'Negative detrainment in thermcell_flux'
            CALL abort_physic(modname,abort_message,1)
         ENDIF
         
!------------------------------------------------------------------------------
! Test sur fraca constante ou croissante au-dessus de lalim
!------------------------------------------------------------------------------
         
! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0
         IF (iflag_thermals_optflux==0) THEN
            DO ig=1,ngrid
               IF (l.ge.lalim(ig).and.l.le.lmax(ig).and.(zw2(ig,l+1).gt.1.e-10).and.(zw2(ig,l).gt.1.e-10) ) THEN
                  zzz = fm(ig,l) * rhobarz(ig,l+1) * zw2(ig,l+1)              &
                  &   / (rhobarz(ig,l) * zw2(ig,l))
                  
                  IF (fm(ig,l+1).gt.zzz) THEN
                     detr(ig,l) = detr(ig,l) + fm(ig,l+1) - zzz
                     fm(ig,l+1) = zzz
                     ncorecfm4  = ncorecfm4 + 1
                  ENDIF
               ENDIF
            ENDDO
         ENDIF
         
!------------------------------------------------------------------------------
! Test sur flux de masse constant ou decroissant au-dessus de lalim
!------------------------------------------------------------------------------
         
! AB : Do we have to decree that? If so, set iflag_thermals_optflux to 0
         IF (iflag_thermals_optflux==0) THEN
            DO ig=1,ngrid
               IF ((fm(ig,l+1).gt.fm(ig,l)).and.(l.gt.lalim(ig))) then
                  f_old = fm(ig,l+1)
                  fm(ig,l+1) = fm(ig,l)
                  detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1)
                  ncorecfm5  = ncorecfm5 + 1
               ENDIF
            ENDDO
         ENDIF
         
!------------------------------------------------------------------------------
! Test detrainment inferieur au flux de masse
!------------------------------------------------------------------------------
         
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : Even if fm has no negative value, it can be lesser than detr.
!      That's not suitable because when we will mix the plume with the
!      environment, it will detrain more mass than it is physically able to do.
!      When it occures, that imply that entr + fm is greater than detr,
!      otherwise we get a negative outgoing mass flux (cf. above).
!      That is why we correct the detrainment as follows.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         
         DO ig=1,ngrid
            IF (detr(ig,l).gt.fm(ig,l)) THEN
               detr(ig,l) = fm(ig,l)
               entr(ig,l) = fm(ig,l+1)
               
               IF (prt_level.ge.10) THEN
                  print *, 'WARNING: Detrainment is modified in thermcell_flux!'
                  print *, 'ig,l,lmax', ig, l, lmax(ig)
               ENDIF
               
               IF (prt_level.ge.100) THEN
                  print *, 'fm-', fm(ig,l)
                  print *, 'entr,detr', entr(ig,l), detr(ig,l)
                  print *, 'fm+', fm(ig,l+1)
               ENDIF
               
               ncorecfm6  = ncorecfm6 + 1
               
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! Dans le cas ou on est au dessus de la couche d'alimentation et que le
! detrainement est plus fort que le flux de masse, on stope le thermique.
! test : on commente
!
! AB : Do we have to stop the plume here?
! AB : Attention, if lmax is modified, be sure that fm, zw2, entr and detr are
!      set to zero above this level!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!               IF (l.gt.lalim(ig)) THEN
!                  lmax(ig)=l
!                  fm(ig,l+1)=0.
!                  entr(ig,l)=0.
!               ELSE
!                  ncorecfm2=ncorecfm2+1
!               ENDIF
            ENDIF
!            
!            IF (l.gt.lmax(ig)) THEN
!               detr(ig,l) = 0.
!               fm(ig,l+1) = 0.
!               entr(ig,l) = 0.
!            ENDIF
         ENDDO
         
!         labort_physic=.false.
!         
!         DO ig=1,ngrid
!            IF (entr(ig,l).lt.0.) THEN
!               labort_physic=.true.
!               igout=ig
!            ENDIF
!         ENDDO
!         
!         IF (labort_physic) THEN
!            print *, 'ERROR: detrainment is greater than mass flux and'
!            print *, '       entrainment cannot compensate it!'
!            print *, 'ig,l,lmax(ig)', igout, l, lmax(igout)
!            print *, '--------------------'
!            print *, 'fm+', fm(igout,lout+1)
!            print *, 'entr,detr', entr(igout,lout), detr(igout,lout)
!            print *, 'fm-', fm(igout,lout)
!            abort_message = 'problem in thermcell_flux'
!            CALL abort_physic (modname,abort_message,1)
!         ENDIF
         
!------------------------------------------------------------------------------
! Test fraction masse entrainee inferieure a fomass_max
!------------------------------------------------------------------------------
         
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : Entrainment is bigger than the maximal authorized value.
!      If we consider that the excess entrainement is in fact plume air which
!      is not detrained then we compensate it by decreasing detr.
!      If it's not enough, we can increase entr in the layer above and decrease
!      the outgoing mass flux in the current layer.
!      If it's still unsufficient, code will abort (now commented).
!      Else we reset entr to its intial value and we renormalize entrainment,
!      detrainment and mass flux profiles such as we do not exceed the maximal
!      authorized entrained mass.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
         
         labort_physic=.false.
         
         DO ig=1,ngrid
            eee0 = entr(ig,l)
            ddd0 = detr(ig,l)
            eee  = entr(ig,l) - masse(ig,l) * fomass_max / ptimestep
            ddd  = detr(ig,l) - eee
            
            IF (eee.gt.0.) THEN
               entr(ig,l) = entr(ig,l) - eee
               ncorecfm3  = ncorecfm3 + 1
               
               IF (prt_level.ge.10) THEN
                  print *, 'WARNING: Entrainment is modified in thermcell_flux!'
                  print *, 'ig,l,lmax', ig, l, lmax(ig)
               ENDIF
               
               IF (ddd.gt.0.) THEN                          ! on reduit le detr dans la couche
                  detr(ig,l) = ddd
               ELSEIF (l.eq.lmax(ig)) THEN                  ! on ajuste le detr dans la couche
                  detr(ig,l) = fm(ig,l) + entr(ig,l)
               ELSEIF (entr(ig,l+1).gt.ABS(ddd)) THEN       ! on reduit le detr dans la couche et l'entr au-dessus
                  detr(ig,l) = 0.
                  fm(ig,l+1) = fm(ig,l) + entr(ig,l)
                  entr(ig,l+1) = entr(ig,l+1) + ddd
!               ELSE
!                  entr(ig,l) = entr(ig,l) + eee
!                  igout=ig
!                  lout=l
!                  labort_physic=.true.
               ELSE
                  fact = (eee0 - eee) / eee0
                  entr(ig,l) = eee0
                  detr(ig,:) = detr(ig,:) * fact
                  entr(ig,:) = entr(ig,:) * fact
                  fm(ig,:)   = fm(ig,:) * fact
                  
                  IF (prt_level.ge.1) THEN
                     print *, 'WARNING: Normalisation is modified in thermcell_flux!'
                     print *, 'old f, new f :', f(ig), fact * f(ig)
                  ENDIF
               ENDIF
            ENDIF
         ENDDO
         
         IF (labort_physic) THEN
            print *, 'ERROR: Entrainment is greater than its maximal authorized value!'
            print *, '       Nor detrainment neither entrainment can compensate it!'
            print *, 'ig,l,entr', igout, lout, entr(igout,lout)
            print *, 'lmin,lmax', lmin(igout), lmax(igout)
            print *, '--------------------'
            print *, 'e_max :', masse(igout,lout)*fomass_max/ptimestep
            print *, '   fomass_max :', fomass_max
            print *, '   masse      :', masse(igout,lout)
            print *, '   ptimestep  :', ptimestep
            print *, 'norm  :', f(igout)
            print *, 'alim* :', alim_star(igout,lout)
            print *, 'entr* :', entr_star(igout,lout)
            print *, '--------------------'
            print *, 'fm l+1   ', fm(igout,lout+2)
            print *, 'entr <-- ', entr(igout,lout+1)
            print *, 'detr  -->', detr(igout,lout+1)
            print *, 'fm l     ', fm(igout,lout+1)
            print *, 'entr <-- ', entr(igout,lout)
            print *, 'detr  -->', detr(igout,lout)
            print *, 'fm l-1   ', fm(igout,lout)
            call flush
            abort_message = 'problem in thermcell_flux'
            CALL abort_physic (modname,abort_message,1)
         ENDIF
         
!------------------------------------------------------------------------------
! Test fraction couverte inferieure a alphamax
!------------------------------------------------------------------------------
         
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! FH : Partie a revisiter.
!      Il semble qu'etaient codees ici deux optiques dans le cas F/(rho*w) > 1
!      - soit limiter la hauteur du thermique en considerant que c'est 
!        la derniere chouche
!      - soit limiter F a rho w.
!      Dans le second cas, il faut en fait limiter a un peu moins que ca parce
!      qu'on a des 1/(1-alpha) un peu plus loin dans thermcell_main et qu'il
!      semble de toutes facons deraisonable d'avoir des fractions de 1...
!      Ci-dessous, et dans l'etat actuel, le premier des deux if est sans doute inutile.  
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~                 
         
         DO ig=1,ngrid
            zfm = rhobarz(ig,l+1) * zw2(ig,l+1) * alphamax
            IF (fm(ig,l+1) .gt. zfm) THEN
               f_old = fm(ig,l+1)
               fm(ig,l+1) = zfm
               detr(ig,l) = detr(ig,l) + f_old - fm(ig,l+1)
               
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : The previous change doesn't observe the equation df/dz = e - d.
!      To avoid this issue, what is better to do? I choose to increase the
!      entrainment in the layer above when l<lmax. If l=lmax then fm(l+1)=0 and
!      there are never any problems.
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
               IF (l.lt.lmax(ig)) THEN
                  entr(ig,l+1) = entr(ig,l+1) + f_old - fm(ig,l+1)
               ENDIF
               
               IF (prt_level.ge.10) THEN
                  print *, 'Mass flux is modified in thermcell_flux'
                  print *, 'ig,l,lmax', ig, l, lmax(ig)
               ENDIF
               
               IF (prt_level.ge.100) THEN
                  print *, 'fm-', fm(ig,l)
                  print *, 'entr,detr', entr(ig,l), detr(ig,l)
                  print *, 'fm+', fm(ig,l+1)
               ENDIF
               
               ncorecalpha = ncorecalpha + 1
               
            ENDIF
         ENDDO
         
      ENDDO
      
!------------------------------------------------------------------------------
! We check if fm, entr and detr vanish correctly in layer lmax
!------------------------------------------------------------------------------
      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : test added to avoid problem when lmax = 0, which is not a fortran index.
!      Theoretically, we always get lmax >= lmin >= linf > 0
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      DO ig=1,ngrid
         IF (lmax(ig).gt.0) THEN
            detr(ig,lmax(ig)) = fm(ig,lmax(ig)) + entr(ig,lmax(ig))
            fm(ig,lmax(ig)+1) = 0.
            entr(ig,lmax(ig)) = 0.
         ENDIF
      ENDDO
      
!------------------------------------------------------------------------------
! Impression des bidouilles utilisees pour corriger les problemes
!------------------------------------------------------------------------------
      
      IF (prt_level.ge.1) THEN
         
         IF (ncorecfm4.ge.1) THEN
            print *, 'WARNING: Deacreasing updraft fraction above lalim!'
            print *, 'in', ncorecfm4, 'point(s)'
         ENDIF
         
         IF (ncorecfm5.ge.1) THEN
            print *, 'WARNING: Increasing mass flux above lalim!'
            print *, 'in', ncorecfm5, 'point(s)'
         ENDIF
         
         IF (ncorecfm6.ge.1) THEN
            print *, 'WARNING: Detrainment is greater than mass flux!'
            print *, 'in', ncorecfm6, 'point(s)'
         ENDIF
         
! AB : those outputs are commented because corresponding test is also commented.
!         IF (ncorecfm2.ge.1) THEN
!            print *, 'WARNING: Detrainment is greater than mass flux!'
!            print *, 'in', ncorecfm2, 'point(s)'
!         ENDIF
         
         IF (ncorecfm3.ge.1) THEN
            print *, 'WARNING: Entrained mass is greater than maximal authorized value!'
            print *, 'in', ncorecfm3, 'point(s)'
         ENDIF
         
         IF (ncorecalpha.ge.1) THEN
            print *, 'WARNING: Updraft fraction is greater than maximal authorized value!'
            print *, 'in', ncorecalpha, 'point(s)'
         ENDIF
         
      ENDIF
      
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! AB : temporary test added to check the validity of eq. df/dz = e - d
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!      DO l=1,nlay
!         DO ig=1,ngrid
!            test = abs(fm(ig,l) + entr(ig,l) - detr(ig,l) - fm(ig,l+1))
!            IF (test.gt.1.e-10) THEN
!               print *, 'WARNING: df/dz != entr - detr'
!               print *, 'ig,l', ig, l
!               print *, 'fm+  :', fm(ig,l+1)
!               print *, 'entr,detr', entr(ig,l), detr(ig,l)
!               print *, 'fm   :', fm(ig,l)
!               print *, 'entr,detr', entr(ig,l-1), detr(ig,l-1)
!               print *, 'fm-  :', fm(ig,l-1)
!               print *, 'err. :', test
!            ENDIF
!         ENDDO
!      ENDDO
      
      
RETURN
END