!******************************************************************** ! 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,TS,BM,Abl,S use param_phy_mod,only:dice,cl IMPLICIT NONE real, dimension(nx,ny) :: precip, pdd, 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, totmelt, 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(:,:)=max(ACC(:,:),0.)/Csnow_2D(:,:) SIMAX(:,:)=max(ACC(:,:),0.)*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(:,:)=max(ACC(:,:),0.)/Csnow_2D(:,:) pr_ice_eq(:,:) = amax1(0.,((PRECIP(:,:)/DICE)-max(ACC(:,:),0.))) ! precipe liquide (ice equivalent) WHERE (PDD(:,:).LE.PDS(:,:)) totmelt(:,:) = Csnow_2D(:,:)*PDD(:,:) ELSEWHERE totmelt(:,:) = Csnow_2D(:,:)*PDS(:,:) + Cice_2D(:,:)*(PDD(:,:)-PDS(:,:)) ENDWHERE snowmelt(:,:) = amin1(totmelt(:,:),max(ACC(:,:),0.)) ! 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*(max(ACC(:,:),0.)-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 bm(:,:) = ACC(:,:)+refreezed_ice(:,:)-totmelt(:,:) SIF(:,:)=refreezed_ice(:,:) ELSE ! pdd standard reeh ! (* Positive degrees required to melt the snow layer *) PDS(:,:)=max(ACC(:,:),0.)/Csnow_2D(:,:) ! (* Maximum amount of super. ice that can be formed *) SIMAX(:,:)=max(ACC(:,:),0.)*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(:,:)=max(ACC(:,:),0.)*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(:,:)) Abl(:,:)=BM(:,:)-Acc(:,:) END SUBROUTINE ABLATION end module ablation_mod