source: trunk/SOURCES/ablation_mod.f90 @ 23

Last change on this file since 23 was 17, checked in by dumas, 9 years ago

Correction bug dans ablation_mod.f90 : dans version précédentes, le SMB est calculé par le PDD Tarasov quel que soit la valeur de PDD_type. Nouvelle version de output_greeneem_mod-0.4.f90 avec calcul du volume sur le Groenland (nouveau découpage plus précis sans Islande et sans Elesmer. climat_forcage_mois_mod.f90 : Lecture du climat t2m et precip dans fichiers NetCDF.

File size: 12.2 KB
Line 
1!********************************************************************
2! Module pour le calcul de l'ablation
3! lecture dans le fichier param de la methode de caclul du BM a utiliser
4! pdd_type : defini le type de pdd utilise
5! pdd_type=0  ! pdd reeh standard
6! pdd_type=1  ! pdd fausto
7! pdd_type=2  ! pdd Tarasov
8module ablation_mod
9
10  use geography, only : nx,ny
11
12  implicit none
13  integer :: pdd_type  ! type de calcul du PDD utilise
14  logical :: annual    ! T = annuel, F = mensuel
15
16  real :: sigma_ice                           ! variabilite Tday
17  real :: Cice                                ! melting factors for ice
18  real :: Csnow                               ! melting factors for snow
19  real :: csi                                 ! proportion of melted water that can refreeze *)
20
21  real, dimension(nx,ny) :: Cice_2D           ! melting factors for ice 2D
22  real, dimension(nx,ny) :: Csnow_2D          ! melting factors for snow 2D
23  real, dimension(nx,ny) :: csi_2D            ! proportion of melted water that can refreeze 2D
24  REAL, dimension(nx,ny) :: sigma_ice_2D      ! sigma fn altitude 2D
25
26  REAL, PARAMETER :: PI_L= 3.1415926
27  REAL, PARAMETER :: DTP=2.0                  ! integrating step for positive degree days (degrees)
28
29  REAL, dimension(nx,ny) :: S22
30  INTEGER, PARAMETER :: NYEAR = 12, NYEAR2 = 360 ! number of months/days in 1 year, st. dev. for temp *)
31  real, parameter :: PYG=2*PI_L/NYEAR            ! ct for PDD calculation (PYG remplace PY qui existe dans CLIMBER
32  REAL, dimension(nx,ny) :: PDDCT                ! ct for PDD calculation
33  REAL, dimension(nx,ny) :: PDDCT2               ! ct for PDD calculation
34
35contains
36
37
38subroutine init_ablation
39 
40  use module3d_phy,only:num_param,num_rep_42
41
42  real :: Cice, Csnow, csi
43 
44! fichiers snapshots
45  namelist/ablation/pdd_type,annual,Cice,Csnow,csi,sigma_ice
46
47! formats pour les ecritures dans 42
48428 format(A)
49  rewind(num_param)        ! pour revenir au debut du fichier param_list.dat
50  read(num_param,ablation)
51  write(num_rep_42,428)'!___________________________________________________________' 
52  write(num_rep_42,428) '&ablation                             ! module ablation_mod'
53  write(num_rep_42,'(A,I2)')       'pdd_type = ',pdd_type
54  write(num_rep_42,'(A,L7)')       'annual = ',annual
55  write(num_rep_42,'(A,F12.6)')    'Cice = ',Cice
56  write(num_rep_42,'(A,F12.6)')    'Csnow = ',Csnow
57  write(num_rep_42,'(A,F12.6)')    'csi = ',csi
58  write(num_rep_42,'(A,F12.6)')    'sigma_ice = ',sigma_ice
59  write(num_rep_42,*)'/'                     
60  write(num_rep_42,428) '! pdd_type : 0 reeh, 1 Fausto, 2 Tarasov'
61  write(num_rep_42,428) '! annual : T = annuel, F = mensuel'
62  write(num_rep_42,428)'! Cice and Csnow, melting factors for ice and snow '
63  write(num_rep_42,428)'! sigma variabilite Tday'
64  write(num_rep_42,428)'! csi proportion of melted water that can refreeze '
65  write(num_rep_42,*)
66
67! initialisation des variables 2D :
68  Cice_2D(:,:)=Cice
69  Csnow_2D(:,:)=Csnow
70  csi_2D(:,:)=csi
71  sigma_ice_2D(:,:)=sigma_ice
72  S22=0.5/sigma_ice/sigma_ice
73  PDDCT=DTP/sigma_ice/SQRT(2.*PI_L)/NYEAR*365.
74  PDDCT2=DTP/sigma_ice/SQRT(2.*PI_L)/NYEAR2*365.
75end subroutine init_ablation
76
77subroutine ablation
78! calcul de l'ablation avec methode annuelle ou mensuelle (flag dans fichier param)
79! remplace l'ancienne subroutine ablation-0.2.f
80!     *********************************************************
81!     (******** ABLATION ********)
82!     *********************************************************
83
84
85  USE module3d_phy,only:Tjuly,Tann,Tmois,acc,pdd,TS,Tshelf,precip,BM,S,dice,cl
86
87  IMPLICIT NONE
88
89  real, dimension(nx,ny) :: pds, simax, pdsi, sif
90  integer :: i,j,k,mo,nday
91  real :: summ
92  real :: temp
93  real,dimension(365) :: TT             !< air temperature yearly cycle, for PDD
94!- SC- Uniquement pour pdd Tarasov
95  REAL, DIMENSION(nx,ny) :: pr_ice_eq, snowmelt, cpsurf
96  REAL, DIMENSION(nx,ny) :: refr2, refreezed_ice
97
98
99! pour pdd fausto calcul Cice_2D, csnow_2D, csi_2D,sigma_ice_2D
100  IF (pdd_type.eq.1) THEN
101      WHERE (Tjuly(:,:).GE.10.)
102          Cice_2D(:,:)=7.
103          Csnow_2D(:,:)=3.
104      ELSEWHERE ((-1.LE.Tjuly(:,:)).AND.(Tjuly(:,:).LT.10.))
105          Cice_2D(:,:)=7.+((15.-7.)/((10.+1)**3))*((10.-Tjuly(:,:))**3)
106          Csnow_2D(:,:)=3.
107      ELSEWHERE
108          Cice_2D(:,:)=15.
109          Csnow_2D(:,:)=3.
110      ENDWHERE
111! pour etre en m/an :
112      Cice_2D(:,:)=Cice_2D(:,:)/1000.
113      Csnow_2D(:,:)=Csnow_2D(:,:)/1000.
114!cdc version std 
115      sigma_ice_2D(:,:)=1.574+0.0012224*S(:,:)
116! calcul de S22, PDDCT et PDDCT2
117      S22(:,:)=0.5/SIGMA_ICE_2D(:,:)/SIGMA_ICE_2D(:,:)
118      PDDCT(:,:)=DTP/SIGMA_ICE_2D(:,:)/SQRT(2.*PI_L)/NYEAR*365.
119      PDDCT2(:,:)=DTP/SIGMA_ICE_2D(:,:)/SQRT(2.*PI_L)/NYEAR2*365.
120
121! calcul du tx de regel
122      WHERE (S(:,:).LE.800)
123          CSI_2D(:,:)=0.
124      ELSEWHERE ((800.LT.S(:,:)).AND.(S(:,:).LT.2000.))
125          CSI_2D(:,:)=(S(:,:)-800.)*0.000833
126      ELSEWHERE
127          CSI_2D(:,:)=1.
128      ENDWHERE
129
130   ELSEIF (pdd_type.EQ.2) THEN ! pdd Tarasov
131
132      where (TJULY(:,:).gt.10.)
133         Cice_2D(:,:) = 8.3*1e-3
134         Csnow_2D(:,:) = 4.3*1e-3
135      endwhere
136      where(TJULY(:,:).gt.-1..and.TJULY(:,:).lt.10.)
137         Cice_2D(:,:) = 1e-3*(8.3+0.0067*(10.-TJULY(:,:))**3)
138         Csnow_2D(:,:) = 1e-3*(2.8+0.15*TJULY(:,:))
139      endwhere
140      where(TJULY(:,:).le.-1.)
141         Cice_2D(:,:) = 17.22*1e-3
142         Csnow_2D(:,:) = 2.65*1e-3
143      endwhere
144
145      sigma_ice_2D(:,:)=sigma_ice
146   ENDIF
147
148
149
150  IF (annual) THEN
151!     (*** positive degree days Tann et Tjuly***)
152      DO i=1,nx
153        DO j=1,ny
154          summ=0.0
155          month_reconstr: DO nday=1,nyear   ! reconstitution du cycle annuel
156            tt(nday)=tann(i,j)+(tjuly(i,j)-tann(i,j))*COS(pyg*nday) 
157            temp=0.0
158            k=1
159            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   
160              summ=summ+temp*EXP(-(temp-tt(nday))*(temp-tt(nday))*s22(i,j)) 
161              temp=temp+dtp
162              k=k+1
163            END DO average_day_reconstr
164          END DO month_reconstr
165          pdd(i,j)=summ*pddct(i,j) 
166        END DO
167      END DO
168  ELSE            ! pdd mensuel
169      DO i=1,nx
170        DO j=1,ny
171          summ=0.0
172! On a deja calcule Tmois(i,j,m) : temperature moyenne de chaque mois
173          month: DO mo=1,12 ! boucle sur les mois
174            temp=0.0          ! variable d'integration
175            k=1
176            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   
177              summ=summ+temp*EXP(-((temp-Tmois(i,j,mo))*(temp-Tmois(i,j,mo)))*s22(i,j))
178              temp=temp+dtp     ! pdtp pas d'integration
179              k= k+1
180            END DO average_day
181          END DO month          !  boucle sur le mois
182          pdd(i,j)=summ*pddct(i,j)   ! pdd pour toute l'annee
183        END DO
184      END DO
185  ENDIF
186
187
188!     calcul du Bilan de masse
189  IF (pdd_type.EQ.1) THEN
190      PDS(:,:)=ACC(:,:)/Csnow_2D(:,:)
191      SIMAX(:,:)=ACC(:,:)*CSI_2D(:,:)
192      PDSI(:,:)=SIMAX(:,:)/Cice_2D(:,:)
193! avec regel de 60% puis fonte (2 premiers where) :
194!      WHERE (PDD(:,:).LE.CSI_2D*PDS(:,:))
195!          BM(:,:)=ACC(:,:)
196!          SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)
197!      endwhere
198!      WHERE ((CSI_2D*PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)))
199!          BM(:,:)=ACC(:,:)+SIMAX(:,:)-PDD(:,:)*Csnow_2D
200!          SIF(:,:)=SIMAX(:,:)
201!      endwhere
202      WHERE (PDD(:,:).LE.PDS(:,:))  ! test avec regel de 60% progressif
203          BM(:,:)=ACC(:,:)-PDD(:,:)*Csnow_2D*(1-CSI_2D(:,:))
204          SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)*(1-CSI_2D(:,:))
205      endwhere
206      WHERE ((PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)+PDSI(:,:)))
207          BM(:,:)=SIMAX(:,:)-(PDD(:,:)-PDS(:,:))*Cice_2D
208          SIF(:,:)=SIMAX(:,:)
209      endwhere
210      WHERE (PDS(:,:)+PDSI(:,:).LE.PDD(:,:))
211          BM(:,:)=(PDS(:,:)+PDSI(:,:)-PDD(:,:))*Cice_2D
212          SIF(:,:)=SIMAX(:,:)
213      endwhere
214 ELSEIF (PDD_type.EQ.2) THEN  ! PDD Tarasov
215        PDS(:,:)=ACC(:,:)/Csnow_2D(:,:)
216        pr_ice_eq(:,:) = amax1(0.,((PRECIP(:,:)/DICE)-ACC(:,:)))  ! precipe liquide (ice equivalent)
217
218        WHERE (PDD(:,:).LE.PDS(:,:))
219           snowmelt(:,:) = Csnow_2D(:,:)*PDD(:,:)
220        ELSEWHERE
221           snowmelt(:,:) = Csnow_2D(:,:)*PDS(:,:) + Cice_2D(:,:)*(PDD(:,:)-PDS(:,:))
222        ENDWHERE
223
224        snowmelt(:,:) = amin1(snowmelt(:,:),ACC(:,:))
225
226! Deux formules possibles pour la capacité calorifique (en J/kg.K):
227! Formule 1 correspond à celle figurant dans prop-thermiques_mod.f90
228! Formule 2 : celle de Tarasov
229        cpsurf(:,:) = 2115.3+7.79293*TANN(:,:) ! Formule Catherine
230!      cpsurf(:,:) = 152.5 + 7.122*TANN(:,:) ! Formule Tarasov
231
232!        where(snowmelt(:,:).lt.ACC(:,:))
233           refr2(:,:) = 2.2*(ACC(:,:)-snowmelt(:,:))-(cpsurf(:,:)/CL)*amin1(TANN(:,:),0.)
234           refreezed_ice(:,:) = amin1(pr_ice_eq(:,:)+snowmelt(:,:),refr2(:,:))
235!        elsewhere
236!           refr2(:,:) = -(cpsurf(:,:)/CL)*amin1(TANN(:,:),0.)
237!           refreezed_ice(:,:) = amin1(pr_ice_eq(:,:)+snowmelt(:,:),refr2(:,:))
238!        endwhere
239        SIMAX(:,:)=refreezed_ice(:,:)
240        PDSI(:,:)=SIMAX(:,:)/Cice_2D(:,:)
241        WHERE (PDD(:,:).LE.PDS(:,:))       
242           BM(:,:)=ACC(:,:)-PDD(:,:)*Csnow_2D(:,:) + SIMAX(:,:)
243           SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)*(1-SIMAX(:,:))
244        endwhere
245        WHERE ((PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)+PDSI(:,:)))
246           BM(:,:)=SIMAX(:,:)-(PDD(:,:)-PDS(:,:))*Cice_2D(:,:)
247           SIF(:,:)=SIMAX(:,:)
248        endwhere
249        WHERE (PDS(:,:)+PDSI(:,:).LE.PDD(:,:))
250           BM(:,:)=(PDS(:,:)+PDSI(:,:)-PDD(:,:))*Cice_2D(:,:)
251           SIF(:,:)=SIMAX(:,:)
252        endwhere
253
254 ELSE ! pdd standard reeh
255!       (* Positive degrees required to melt the snow layer *)
256        PDS(:,:)=ACC(:,:)/Csnow_2D(:,:)
257!       (* Maximum amount of super. ice that can be formed *)
258        SIMAX(:,:)=ACC(:,:)*CSI_2D(:,:)
259!       (* Pos. degrees required to melt the superimposed ice *)
260        PDSI(:,:)=SIMAX(:,:)/Cice_2D(:,:)
261! avec regel de 60% puis fonte (2 premiers where) :
262      WHERE (PDD(:,:).LE.CSI_2D*PDS(:,:))
263          BM(:,:)=ACC(:,:)
264          SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)
265      endwhere
266      WHERE ((CSI_2D*PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)))
267          BM(:,:)=ACC(:,:)+SIMAX(:,:)-PDD(:,:)*Csnow_2D
268          SIF(:,:)=SIMAX(:,:)
269      endwhere
270!       WHERE (PDD(:,:).LE.PDS(:,:))                  ! test avec regel de 60% progressif|
271!           BM(:,:)=ACC(:,:)-PDD(:,:)*Csnow_2D*(1-CSI_2D)   ! remplace les 2 premiers where    |
272!           SIF(:,:)=PDD(:,:)*Csnow_2D(:,:)*(1-CSI_2D(:,:)) !----------------------------------|
273!      endwhere
274      WHERE ((PDS(:,:).LT.PDD(:,:)).AND.(PDD(:,:).LE.PDS(:,:)+PDSI(:,:)))
275          BM(:,:)=SIMAX(:,:)-(PDD(:,:)-PDS(:,:))*Cice_2D
276          SIF(:,:)=SIMAX(:,:)
277      endwhere
278      WHERE (PDS(:,:)+PDSI(:,:).LE.PDD(:,:))
279          BM(:,:)=(PDS(:,:)+PDSI(:,:)-PDD(:,:))*Cice_2D
280          SIF(:,:)=SIMAX(:,:)
281      endwhere
282
283  ENDIF
284
285! calcul de la temperature de surface (utilisee dans icetemp) :
286   TS(:,:)=(TANN(:,:)+26.6*SIF(:,:))
287   TS(:,:)=min(0.0,TS(:,:))
288   tshelf(:,:)=TS(:,:)
289
290!!$!     (*** ablation ***)
291!!$  DO i=1,nx
292!!$    DO j=1,ny
293!!$!       (* positive degrees required to melt the snow layer *)
294!!$      pds=acc(i,j)/csnow
295!!$!       (* maximum amount of super. ice that can be formed *)
296!!$      simax=acc(i,j)*csi
297!!$!       (* pos. degrees required to melt the superimposed ice *)
298!!$      pdsi=simax/cice
299!!$      IF (pdd(i,j).LE.csi*pds) THEN
300!!$          sif=pdd(i,j)*csnow
301!!$      ELSE
302!!$          sif=simax
303!!$      ENDIF
304!!$
305!!$!       surface temperature
306!!$      ts(i,j)=(tann(i,j)+26.6*sif)
307!!$      ts(i,j)=MIN(0.0,ts(i,j))
308!!$      tshelf(i,j)=ts(i,j)
309!!$
310!!$!       mass balance
311!!$      IF (pdd(i,j).LE.csi*pds)  bm(i,j)=acc(i,j) ! toute la fonte est recuperee en regel
312!!$      IF ((csi*pds.LT.pdd(i,j)).AND.(pdd(i,j).LE.pds)) & ! fonte de la neige uniquement
313!!$         bm(i,j)=acc(i,j)+simax-pdd(i,j)*csnow
314!!$      IF ((pds.LT.pdd(i,j)).AND.(pdd(i,j).LE.pds+pdsi)) & ! fonte toute neige et une partie glace regel
315!!$         bm(i,j)=simax-(pdd(i,j)-pds)*cice
316!!$      IF (pds+pdsi.LE.pdd(i,j)) bm(i,j)=(pds+pdsi-pdd(i,j))*cice !fonte de tt neige tt regel + glace annees precedentes
317
318         
319!940   format('%%%% ',a,'   time=',f8.0,' %%%%')
320END SUBROUTINE ABLATION
321
322end module ablation_mod
Note: See TracBrowser for help on using the repository browser.