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