1 | !> \file bilan_eau.f90 |
---|
2 | !! Clacul du bilan d'eau sur la calotte |
---|
3 | !< |
---|
4 | |
---|
5 | |
---|
6 | !> MODULE: bilan_eau |
---|
7 | !! Calcule le bilan d'eau sur la grille GRISLI |
---|
8 | !! \author C. Dumas |
---|
9 | !! \date 20/02/2017 |
---|
10 | !! @note Permet de vérifier que le bilan d'eau est ok et qu'il n'y a pas de fuite. |
---|
11 | !> |
---|
12 | |
---|
13 | ! bilan d'eau : |
---|
14 | ! ce qui tombe : Acc |
---|
15 | ! ce qui sort : Bmelt, Abl, calving, ablbord |
---|
16 | ! ce qui est stocké : masse de glace |
---|
17 | ! le total de tous ces termes doit être constant (sur la grille, pas localement) |
---|
18 | |
---|
19 | module bilan_eau_mod |
---|
20 | |
---|
21 | |
---|
22 | use module3D_phy |
---|
23 | implicit none |
---|
24 | |
---|
25 | real :: sum_H |
---|
26 | |
---|
27 | |
---|
28 | real,dimension(nx,ny) :: tot_water !< bilan d'eau |
---|
29 | real,dimension(nx,ny) :: H_beau_old, diff_H_2D, diff_H_water_bilan |
---|
30 | real,dimension(nx,ny) :: Bm_dtt !< mass balance on ice points accumulated during dtt |
---|
31 | real,dimension(nx,ny) :: Bmelt_dtt !< basal melting on ice points accumulated during dtt |
---|
32 | real,dimension(nx,ny) :: calv_dtt !< calving sur dtt (pour calcul bilan d'eau) |
---|
33 | real,dimension(nx,ny) :: archimtab ! point pose si > 0 |
---|
34 | real,dimension(nx,ny) :: grline_dtt ! grounding line flux during dtt |
---|
35 | |
---|
36 | real :: sum_H_old |
---|
37 | real :: diff_H |
---|
38 | real :: alpha_flot ! ro/row |
---|
39 | |
---|
40 | real,dimension(nx,ny) :: bm_dt,bmelt_dt |
---|
41 | |
---|
42 | |
---|
43 | contains |
---|
44 | |
---|
45 | subroutine init_bilan_eau |
---|
46 | ! initialisation des variables |
---|
47 | diff_H=0. |
---|
48 | sum_H_old = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1) |
---|
49 | tot_water(:,:)=0. |
---|
50 | bm_dt(:,:)=0. |
---|
51 | bmelt_dt(:,:)=0. |
---|
52 | alpha_flot=ro/row |
---|
53 | archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel |
---|
54 | gr_line(:,:)=0 |
---|
55 | do j=1,ny |
---|
56 | do i=1,nx |
---|
57 | !afq if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then ! grounded with ice |
---|
58 | if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice |
---|
59 | if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) gr_line(i,j)=1 |
---|
60 | if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) gr_line(i,j)=1 |
---|
61 | if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) gr_line(i,j)=1 |
---|
62 | if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) gr_line(i,j)=1 |
---|
63 | endif |
---|
64 | enddo |
---|
65 | enddo |
---|
66 | |
---|
67 | end subroutine init_bilan_eau |
---|
68 | |
---|
69 | |
---|
70 | |
---|
71 | |
---|
72 | subroutine bilan_eau |
---|
73 | |
---|
74 | tot_water(:,:)=0. |
---|
75 | |
---|
76 | bm_dt(:,:)=0. |
---|
77 | bmelt_dt(:,:)=0. |
---|
78 | |
---|
79 | bm_dt(2:nx-1,2:ny-1)=bm(2:nx-1,2:ny-1)*dt |
---|
80 | bmelt_dt(2:nx-1,2:ny-1)=bmelt(2:nx-1,2:ny-1)*dt |
---|
81 | !~ where (Bm(:,:).lt.0.) |
---|
82 | !~ bm_dt(:,:)=max(Bm(:,:)*dt,-H(:,:)) |
---|
83 | !~ elsewhere |
---|
84 | !~ bm_dt(:,:)=bm(:,:)*dt |
---|
85 | !~ endwhere |
---|
86 | |
---|
87 | !~ where (Bmelt(:,:).lt.0.) |
---|
88 | !~ bmelt_dt(:,:)=max(Bmelt(:,:)*dt,-H(:,:)) |
---|
89 | !~ elsewhere |
---|
90 | !~ bmelt_dt(:,:)=bmelt(:,:)*dt |
---|
91 | !~ endwhere |
---|
92 | |
---|
93 | |
---|
94 | ! on fait la somme sur la grille en excluant les bords : |
---|
95 | |
---|
96 | sum_H = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1) |
---|
97 | |
---|
98 | diff_H = diff_H + (sum_H - sum_H_old) ! on calcul la variation de volume a chaque dt |
---|
99 | |
---|
100 | diff_H_2D(2:nx-1,2:ny-1)=H(2:nx-1,2:ny-1)-H_beau_old(2:nx-1,2:ny-1) |
---|
101 | |
---|
102 | where (ice(2:nx-1,2:ny-1).eq.1) |
---|
103 | Bm_dtt(2:nx-1,2:ny-1) = Bm_dtt(2:nx-1,2:ny-1) + Bm_dt(2:nx-1,2:ny-1) !* dt ! somme Bm sur dt |
---|
104 | bmelt_dtt(2:nx-1,2:ny-1) = bmelt_dtt(2:nx-1,2:ny-1) + bmelt_dt(2:nx-1,2:ny-1) ! * dt ! somme bmelt sur dt |
---|
105 | endwhere |
---|
106 | |
---|
107 | archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel |
---|
108 | gr_line(:,:)=0 |
---|
109 | do j=1,ny |
---|
110 | do i=1,nx |
---|
111 | !afq if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel)) then ! grounded with ice |
---|
112 | if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice |
---|
113 | if (archimtab(i-1,j).LT.0..and.Uxbar(i,j).LT.0..and..not.flot_marais(i-1,j)) gr_line(i,j)=1 |
---|
114 | if (archimtab(i+1,j).LT.0..and.Uxbar(i+1,j).GT.0..and..not.flot_marais(i+1,j)) gr_line(i,j)=1 |
---|
115 | if (archimtab(i,j-1).LT.0..and.Uybar(i,j).LT.0..and..not.flot_marais(i,j-1)) gr_line(i,j)=1 |
---|
116 | if (archimtab(i,j+1).LT.0..and.Uybar(i,j+1).GT.0..and..not.flot_marais(i,j+1)) gr_line(i,j)=1 |
---|
117 | endif |
---|
118 | enddo |
---|
119 | enddo |
---|
120 | |
---|
121 | where (gr_line(:,:).eq.1) |
---|
122 | !~ grline_dtt(:,:)= (((uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))**2+ & |
---|
123 | !~ (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5 & |
---|
124 | !~ *H(:,:) + grline_dtt(:,:) |
---|
125 | grline_dtt(:,:)= - sqrt( & |
---|
126 | ( (uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))*dy/2. )**2 + & |
---|
127 | ( (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))*dx/2. )**2 ) & |
---|
128 | * H(:,:) * dt + grline_dtt(:,:) |
---|
129 | endwhere |
---|
130 | |
---|
131 | |
---|
132 | if (isynchro.eq.1) then |
---|
133 | ! on raisonne en bilan annuel pour simplifier : |
---|
134 | !~ where (ice(2:nx-1,2:ny-1).eq.1) |
---|
135 | !~ tot_water(2:nx-1,2:ny-1) = (Bm_dtt(2:nx-1,2:ny-1) - Bmelt_dtt(2:nx-1,2:ny-1) + calv_dtt(2:nx-1,2:ny-1) - ablbord_dtt(2:nx-1,2:ny-1)) / dtt |
---|
136 | !~ elsewhere |
---|
137 | !~ tot_water(2:nx-1,2:ny-1) = (calv_dtt(2:nx-1,2:ny-1) - ablbord_dtt(2:nx-1,2:ny-1)) / dtt |
---|
138 | !~ endwhere |
---|
139 | !cdc pas besoin du test sur ice ici, il a ete fait avant (et le masque ice varie a chaque dt) |
---|
140 | tot_water(2:nx-1,2:ny-1) = (Bm_dtt(2:nx-1,2:ny-1) - Bmelt_dtt(2:nx-1,2:ny-1) + calv_dtt(2:nx-1,2:ny-1) - ablbord_dtt(2:nx-1,2:ny-1)) / dtt |
---|
141 | |
---|
142 | ! bilan d'eau sur la grille : |
---|
143 | water_bilan=sum(tot_water(:,:)) |
---|
144 | diff_H = diff_H/dtt |
---|
145 | |
---|
146 | !999 format(f0.2,1x,e15.8,1x,i10,8(1x,e15.8)) |
---|
147 | ! write(6,999),time,sum_H,count(ice(:,:)==1),diff_H,water_bilan,sum(calv_dtt(:,:))/dtt,sum(ablbord_dtt(:,:))/dtt,sum(bmelt_dtt(:,:),mask=ice(:,:)==1)/dtt,sum(bm(:,:),mask=ice(:,:)==1),sum(Bm_dtt(:,:))/dtt,sum(bmelt_dtt(:,:))/dtt |
---|
148 | diff_H_water_bilan(2:nx-1,2:ny-1)=tot_water(2:nx-1,2:ny-1)-diff_H_2D(2:nx-1,2:ny-1) |
---|
149 | |
---|
150 | endif |
---|
151 | sum_H_old=sum_H |
---|
152 | |
---|
153 | H_beau_old(:,:)=H(:,:) |
---|
154 | |
---|
155 | end subroutine bilan_eau |
---|
156 | |
---|
157 | end module bilan_eau_mod |
---|