source: trunk/SOURCES/bilan_eau_mod.f90 @ 192

Last change on this file since 192 was 192, checked in by aquiquet, 6 years ago

Model performance improved with a reduced domain extension for elliptic equation resolution for buttressing estimation

File size: 5.6 KB
Line 
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
19module bilan_eau_mod
20
21
22use module3D_phy
23implicit none
24
25real :: sum_H
26
27
28real,dimension(nx,ny) :: tot_water    !< bilan d'eau
29real,dimension(nx,ny) :: H_beau_old, diff_H_2D, diff_H_water_bilan
30real,dimension(nx,ny) :: Bm_dtt       !< mass balance on ice points accumulated during dtt
31real,dimension(nx,ny) :: Bmelt_dtt    !< basal melting on ice points accumulated during dtt
32real,dimension(nx,ny) :: calv_dtt     !< calving sur dtt (pour calcul bilan d'eau)
33real,dimension(nx,ny) :: archimtab    ! point pose si > 0
34real,dimension(nx,ny) :: grline_dtt   ! grounding line flux during dtt
35
36real :: sum_H_old
37real :: diff_H
38real :: alpha_flot                      ! ro/row
39
40real,dimension(nx,ny) :: bm_dt,bmelt_dt
41
42
43contains
44
45subroutine 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
67end subroutine init_bilan_eau
68
69       
70
71
72subroutine bilan_eau
73
74tot_water(:,:)=0.
75
76bm_dt(:,:)=0.
77bmelt_dt(:,:)=0.
78
79bm_dt(2:nx-1,2:ny-1)=bm(2:nx-1,2:ny-1)*dt
80bmelt_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
96sum_H = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1)
97
98diff_H = diff_H + (sum_H - sum_H_old) ! on calcul la variation de volume a chaque dt
99
100diff_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
102where (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
105endwhere
106
107archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel
108gr_line(:,:)=0
109do 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
119enddo
120
121where (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(:,:)
129endwhere
130
131
132if (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
150endif
151sum_H_old=sum_H
152
153H_beau_old(:,:)=H(:,:)
154
155end subroutine bilan_eau
156
157end module bilan_eau_mod
Note: See TracBrowser for help on using the repository browser.