source: trunk/SOURCES/ablation_mod.f90 @ 170

!********************************************************************
! Module pour le calcul de l'ablation
! lecture dans le fichier param de la methode de caclul du BM a utiliser
! pdd_type : defini le type de pdd utilise
! pdd_type=0  ! pdd reeh standard
! pdd_type=1  ! pdd fausto
! pdd_type=2  ! pdd Tarasov
module ablation_mod

  use geography, only : nx,ny

  implicit none
  integer :: pdd_type  ! type de calcul du PDD utilise
  logical :: annual    ! T = annuel, F = mensuel

  real :: sigma_ice                           ! variabilite Tday
  real :: Cice                                ! melting factors for ice 
  real :: Csnow                               ! melting factors for snow 
  real :: csi                                 ! proportion of melted water that can refreeze *)

  real, dimension(nx,ny) :: Cice_2D           ! melting factors for ice 2D 
  real, dimension(nx,ny) :: Csnow_2D          ! melting factors for snow 2D
  real, dimension(nx,ny) :: csi_2D            ! proportion of melted water that can refreeze 2D
  REAL, dimension(nx,ny) :: sigma_ice_2D      ! sigma fn altitude 2D

  REAL, PARAMETER :: PI_L= 3.1415926
  REAL, PARAMETER :: DTP=2.0                  ! integrating step for positive degree days (degrees)

  REAL, dimension(nx,ny) :: S22
  INTEGER, PARAMETER :: NYEAR = 12, NYEAR2 = 360 ! number of months/days in 1 year, st. dev. for temp *)
  real, parameter :: PYG=2*PI_L/NYEAR            ! ct for PDD calculation (PYG remplace PY qui existe dans CLIMBER 
  REAL, dimension(nx,ny) :: PDDCT                ! ct for PDD calculation 
  REAL, dimension(nx,ny) :: PDDCT2               ! ct for PDD calculation
  integer :: methode_abl=0                       ! selection methode pdd (0) ou van den Berg 2008 (1)

contains


subroutine init_ablation
  
  use module3d_phy,only:num_param,num_rep_42

  real :: Cice, Csnow, csi
  
! fichiers snapshots
  namelist/ablation/pdd_type,annual,Cice,Csnow,csi,sigma_ice

! formats pour les ecritures dans 42
428 format(A)
  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
  read(num_param,ablation)
  write(num_rep_42,428)'!___________________________________________________________' 
  write(num_rep_42,428) '&ablation                             ! module ablation_mod'
  write(num_rep_42,'(A,I2)')       'pdd_type = ',pdd_type
  write(num_rep_42,'(A,L7)')       'annual = ',annual
  write(num_rep_42,'(A,F12.6)')    'Cice = ',Cice
  write(num_rep_42,'(A,F12.6)')    'Csnow = ',Csnow
  write(num_rep_42,'(A,F12.6)')    'csi = ',csi
  write(num_rep_42,'(A,F12.6)')    'sigma_ice = ',sigma_ice
  write(num_rep_42,*)'/'                      
  write(num_rep_42,428) '! pdd_type : 0 reeh, 1 Fausto, 2 Tarasov'
  write(num_rep_42,428) '! annual : T = annuel, F = mensuel'
  write(num_rep_42,428)'! Cice and Csnow, melting factors for ice and snow '
  write(num_rep_42,428)'! sigma variabilite Tday'
  write(num_rep_42,428)'! csi proportion of melted water that can refreeze '
  write(num_rep_42,*)

! initialisation des variables 2D :
  Cice_2D(:,:)=Cice
  Csnow_2D(:,:)=Csnow
  csi_2D(:,:)=csi
  sigma_ice_2D(:,:)=sigma_ice
  S22=0.5/sigma_ice/sigma_ice
  PDDCT=DTP/sigma_ice/SQRT(2.*PI_L)/NYEAR*365.
  PDDCT2=DTP/sigma_ice/SQRT(2.*PI_L)/NYEAR2*365.
end subroutine init_ablation

subroutine ablation
! calcul de l'ablation avec methode annuelle ou mensuelle (flag dans fichier param)
! remplace l'ancienne subroutine ablation-0.2.f
!     *********************************************************
!     (******** ABLATION ********)
!     *********************************************************


  USE module3d_phy,only:Tjuly,Tann,Tmois,acc,pdd,TS,Tshelf,precip,BM,Abl,S,dice,cl

  IMPLICIT NONE

  real, dimension(nx,ny) :: pds, simax, pdsi, sif
  integer :: i,j,k,mo,nday
  real :: summ
  real :: temp
  real,dimension(365) :: TT             !< air temperature yearly cycle, for PDD
!- SC- Uniquement pour pdd Tarasov
  REAL, DIMENSION(nx,ny) :: pr_ice_eq, snowmelt, cpsurf
  REAL, DIMENSION(nx,ny) :: refr2, refreezed_ice


! pour pdd fausto calcul Cice_2D, csnow_2D, csi_2D,sigma_ice_2D
  IF (pdd_type.eq.1) THEN
      WHERE (Tjuly(:,:).GE.10.)
          Cice_2D(:,:)=7.
          Csnow_2D(:,:)=3.
      ELSEWHERE ((-1.LE.Tjuly(:,:)).AND.(Tjuly(:,:).LT.10.))
          Cice_2D(:,:)=7.+((15.-7.)/((10.+1)**3))*((10.-Tjuly(:,:))**3)
          Csnow_2D(:,:)=3.
      ELSEWHERE
          Cice_2D(:,:)=15.
          Csnow_2D(:,:)=3.
      ENDWHERE
! pour etre en m/an :
      Cice_2D(:,:)=Cice_2D(:,:)/1000.
      Csnow_2D(:,:)=Csnow_2D(:,:)/1000.
!cdc version std  
      sigma_ice_2D(:,:)=1.574+0.0012224*S(:,:)
! calcul de S22, PDDCT et PDDCT2
      S22(:,:)=0.5/SIGMA_ICE_2D(:,:)/SIGMA_ICE_2D(:,:)
      PDDCT(:,:)=DTP/SIGMA_ICE_2D(:,:)/SQRT(2.*PI_L)/NYEAR*365.
      PDDCT2(:,:)=DTP/SIGMA_ICE_2D(:,:)/SQRT(2.*PI_L)/NYEAR2*365.

! calcul du tx de regel
      WHERE (S(:,:).LE.800)
          CSI_2D(:,:)=0.
      ELSEWHERE ((800.LT.S(:,:)).AND.(S(:,:).LT.2000.))
          CSI_2D(:,:)=(S(:,:)-800.)*0.000833
      ELSEWHERE
          CSI_2D(:,:)=1.
      ENDWHERE

  ELSEIF (pdd_type.EQ.2) THEN ! pdd Tarasov

      where (TJULY(:,:).gt.10.)
        Cice_2D(:,:) = 8.3*1e-3
        Csnow_2D(:,:) = 4.3*1e-3
      endwhere
      where(TJULY(:,:).gt.-1..and.TJULY(:,:).lt.10.)
        Cice_2D(:,:) = 1e-3*(8.3+0.0067*(10.-TJULY(:,:))**3)
        Csnow_2D(:,:) = 1e-3*(2.8+0.15*TJULY(:,:))
      endwhere
      where(TJULY(:,:).le.-1.)
        Cice_2D(:,:) = 17.22*1e-3
        Csnow_2D(:,:) = 2.65*1e-3
      endwhere

      sigma_ice_2D(:,:)=sigma_ice
  ENDIF



  IF (annual) THEN
!     (*** positive degree days Tann et Tjuly***)
      DO i=1,nx
        DO j=1,ny
          summ=0.0
          month_reconstr: DO nday=1,nyear   ! reconstitution du cycle annuel
            tt(nday)=tann(i,j)+(tjuly(i,j)-tann(i,j))*COS(pyg*nday) 
            temp=0.0
            k=1
            average_day_reconstr: DO WHILE ((temp.LE.tt(nday)+2.5*sigma_ice_2D(i,j)).AND.k.LE.50) ! pdd d'un jour par mois    
              summ=summ+temp*EXP(-(temp-tt(nday))*(temp-tt(nday))*s22(i,j)) 
              temp=temp+dtp
              k=k+1
            END DO average_day_reconstr
          END DO month_reconstr
          pdd(i,j)=summ*pddct(i,j) 
        END DO
      END DO
  ELSE            ! pdd mensuel
      DO i=1,nx
        DO j=1,ny
          summ=0.0
! On a deja calcule Tmois(i,j,m) : temperature moyenne de chaque mois
          month: DO mo=1,12 ! boucle sur les mois
            temp=0.0          ! variable d'integration
            k=1
            average_day: DO WHILE ((temp.LE.Tmois(i,j,mo)+2.5*sigma_ice_2D(i,j)).AND.k.LE.50) ! pdd d'un jour par mois    
              summ=summ+temp*EXP(-((temp-Tmois(i,j,mo))*(temp-Tmois(i,j,mo)))*s22(i,j))
              temp=temp+dtp     ! pdtp pas d'integration
              k= k+1
            END DO average_day
          END DO month          !  boucle sur le mois
          pdd(i,j)=summ*pddct(i,j)   ! pdd pour toute l'annee
        END DO
      END DO
  ENDIF


!     calcul du Bilan de masse
  IF (methode_abl.EQ.0) THEN
     IF (pdd_type.EQ.1) THEN
        PDS(:,:)=ACC(:,:)/Csnow_2D(:,:)
        SIMAX(:,:)=ACC(:,:)*CSI_2D(:,:)
        PDSI(:,:)=SIMAX(:,:)/Cice_2D(:,:)
! avec regel de 60% puis fonte (2 premiers where) :
!      WHERE (PDD(:,:).LE.CSI_2D*PDS(:,:))
!          BM(:,:)=ACC(:,:)
!          SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)
!      endwhere
!      WHERE ((CSI_2D*PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)))
!          BM(:,:)=ACC(:,:)+SIMAX(:,:)-PDD(:,:)*Csnow_2D
!          SIF(:,:)=SIMAX(:,:)
!      endwhere
        WHERE (PDD(:,:).LE.PDS(:,:))  ! test avec regel de 60% progressif
           BM(:,:)=ACC(:,:)-PDD(:,:)*Csnow_2D*(1-CSI_2D(:,:))
           SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)*(1-CSI_2D(:,:))
        endwhere
        WHERE ((PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)+PDSI(:,:)))
           BM(:,:)=SIMAX(:,:)-(PDD(:,:)-PDS(:,:))*Cice_2D
           SIF(:,:)=SIMAX(:,:)
        endwhere
        WHERE (PDS(:,:)+PDSI(:,:).LE.PDD(:,:))
           BM(:,:)=(PDS(:,:)+PDSI(:,:)-PDD(:,:))*Cice_2D
           SIF(:,:)=SIMAX(:,:)
        endwhere
     ELSEIF (PDD_type.EQ.2) THEN  ! PDD Tarasov
        PDS(:,:)=ACC(:,:)/Csnow_2D(:,:)
        pr_ice_eq(:,:) = amax1(0.,((PRECIP(:,:)/DICE)-ACC(:,:)))  ! precipe liquide (ice equivalent)

        WHERE (PDD(:,:).LE.PDS(:,:))
           snowmelt(:,:) = Csnow_2D(:,:)*PDD(:,:)
        ELSEWHERE
           snowmelt(:,:) = Csnow_2D(:,:)*PDS(:,:) + Cice_2D(:,:)*(PDD(:,:)-PDS(:,:))
        ENDWHERE

        snowmelt(:,:) = amin1(snowmelt(:,:),ACC(:,:))

! Deux formules possibles pour la capacité calorifique (en J/kg.K):
! Formule 1 correspond à celle figurant dans prop-thermiques_mod.f90
! Formule 2 : celle de Tarasov
        cpsurf(:,:) = 2115.3+7.79293*TANN(:,:) ! Formule Catherine
!      cpsurf(:,:) = 152.5 + 7.122*TANN(:,:) ! Formule Tarasov

!        where(snowmelt(:,:).lt.ACC(:,:))
        refr2(:,:) = 2.2*(ACC(:,:)-snowmelt(:,:))-(cpsurf(:,:)/CL)*amin1(TANN(:,:),0.)
        refreezed_ice(:,:) = amin1(pr_ice_eq(:,:)+snowmelt(:,:),refr2(:,:))
!        elsewhere
!           refr2(:,:) = -(cpsurf(:,:)/CL)*amin1(TANN(:,:),0.)
!           refreezed_ice(:,:) = amin1(pr_ice_eq(:,:)+snowmelt(:,:),refr2(:,:))
!        endwhere
        SIMAX(:,:)=refreezed_ice(:,:)
        PDSI(:,:)=SIMAX(:,:)/Cice_2D(:,:)
        WHERE (PDD(:,:).LE.PDS(:,:))        
           BM(:,:)=ACC(:,:)-PDD(:,:)*Csnow_2D(:,:) + SIMAX(:,:)
           SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)*(1-SIMAX(:,:))
        endwhere
        WHERE ((PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)+PDSI(:,:)))
           BM(:,:)=SIMAX(:,:)-(PDD(:,:)-PDS(:,:))*Cice_2D(:,:)
           SIF(:,:)=SIMAX(:,:)
        endwhere
        WHERE (PDS(:,:)+PDSI(:,:).LE.PDD(:,:))
           BM(:,:)=(PDS(:,:)+PDSI(:,:)-PDD(:,:))*Cice_2D(:,:)
           SIF(:,:)=SIMAX(:,:)
        endwhere

     ELSE ! pdd standard reeh
!       (* Positive degrees required to melt the snow layer *)
        PDS(:,:)=ACC(:,:)/Csnow_2D(:,:)
!       (* Maximum amount of super. ice that can be formed *)
        SIMAX(:,:)=ACC(:,:)*CSI_2D(:,:)
!       (* Pos. degrees required to melt the superimposed ice *)
        PDSI(:,:)=SIMAX(:,:)/Cice_2D(:,:)
! avec regel de 60% puis fonte (2 premiers where) :
        WHERE (PDD(:,:).LE.CSI_2D*PDS(:,:))
           BM(:,:)=ACC(:,:)
           SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)
        endwhere
        WHERE ((CSI_2D*PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)))
           BM(:,:)=ACC(:,:)+SIMAX(:,:)-PDD(:,:)*Csnow_2D
           SIF(:,:)=SIMAX(:,:)
        endwhere
!       WHERE (PDD(:,:).LE.PDS(:,:))                  ! test avec regel de 60% progressif|
!           BM(:,:)=ACC(:,:)-PDD(:,:)*Csnow_2D*(1-CSI_2D)   ! remplace les 2 premiers where    |
!           SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)*(1-CSI_2D(:,:)) !----------------------------------|
!      endwhere
        WHERE ((PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)+PDSI(:,:)))
           BM(:,:)=SIMAX(:,:)-(PDD(:,:)-PDS(:,:))*Cice_2D
           SIF(:,:)=SIMAX(:,:)
        endwhere
        WHERE (PDS(:,:)+PDSI(:,:).LE.PDD(:,:))
           BM(:,:)=(PDS(:,:)+PDSI(:,:)-PDD(:,:))*Cice_2D
           SIF(:,:)=SIMAX(:,:)
        endwhere
     ENDIF
  ELSEIF ( methode_abl.EQ.1 ) THEN ! Insolation Temperature Melt equation, van den Berg, 2008
     SIMAX(:,:)=ACC(:,:)*CSI
     SIF(:,:)=SIMAX(:,:)
  ELSE
     print*,'ablation.f90 : pb methode calcul BM'
     print*,'ATTENTION methode_abl > 1 NON COMPATIBLE AVEC iLOVECLIM' 
     stop
  ENDIF

! calcul de la temperature de surface (utilisee dans icetemp) :
  TS(:,:)=(TANN(:,:)+26.6*SIF(:,:))
  TS(:,:)=min(0.0,TS(:,:))
  tshelf(:,:)=TS(:,:)
  Abl(:,:)=BM(:,:)-Acc(:,:)

END SUBROUTINE ABLATION

end module ablation_mod