source: trunk/SOURCES/ablation-0.2.f @ 10

!> \file ablation-0.2.f
!! Fichier contenant le subroutine Ablation 
!<
      
!> Subroutine ABLATION
!! Calcule l'ablation 
!! \author ...
!! \date ...
!! @note Used modules:
!! @note  - module3D_phy
!<


c********************************************************************

      subroutine ABLATION()

c     *********************************************************
c     (******** ABLATION ********)
c     *********************************************************


      use module3d_phy


      implicit none

c     (*** positive degree days ***)
      do i=1,nx
      do j=1,ny
!          tjuly(i,j)=min(tjuly(i,j),-4.0)
!          tann(i,j)=min(tann(i,j),tjuly(i,j))
        summ=0.0
        do nday=1,nyear   
          tt(nday)=tann(i,j)+(tjuly(i,j)-tann(i,j))*cos(py*nday) 
          temp=0.0 
          do k=1,50
            if (temp.gt.tt(nday)+2.5*sigma) goto 205
            summ=summ+temp*exp(-(temp-tt(nday))*(temp-tt(nday))*s22) 
            temp=temp+dtp
          end do
205       continue
        end do
        pdd(i,j)=summ*pddct 
      end do
      end do

c     (*** ablation ***)
      do i=1,nx
      do j=1,ny
c       (* positive degrees required to melt the snow layer *)
        pds=acc(i,j)/csnow 
c       (* maximum amount of super. ice that can be formed *)
        simax=acc(i,j)*csi
c       (* pos. degrees required to melt the superimposed ice *)
        pdsi=simax/cice
        if (pdd(i,j).le.csi*pds) then
          sif=pdd(i,j)*csnow
        else 
          sif=simax 
        endif

c       surface temperature
        ts(i,j)=(tann(i,j)+26.6*sif)
        ts(i,j)=min(0.0,ts(i,j))
        tshelf(i,j)=ts(i,j)

c       mass balance
        if (pdd(i,j).le.csi*pds)  bm(i,j)=acc(i,j) 
        if ((csi*pds.lt.pdd(i,j)).and.(pdd(i,j).le.pds))
     1   bm(i,j)=acc(i,j)+simax-pdd(i,j)*csnow 
        if ((pds.lt.pdd(i,j)).and.(pdd(i,j).le.pds+pdsi))
     1   bm(i,j)=simax-(pdd(i,j)-pds)*cice
        if (pds+pdsi.le.pdd(i,j)) bm(i,j)=(pds+pdsi-pdd(i,j))*cice

c       
c        modif pour l'experience l3a
         if (icouple.eq.2) then
          bm(i,j)=acc(i,j)
          ts(i,j)=tann(i,j)
         endif
          

      end do
      end do

!940   format('%%%% ',a,'   time=',f8.0,' %%%%')
      end