Changeset 378

Timestamp:

03/21/23 11:11:44 (13 months ago)

Author:

dumas

Message:

use only in module bilan_eau_mod

File:

• branches/GRISLIv3/SOURCES/bilan_eau_mod.f90 (modified) (1 diff)

branches/GRISLIv3/SOURCES/bilan_eau_mod.f90

@@ -19 +19 @@
 module bilan_eau_mod
 
+  use param_phy_mod, only: ro,row
+  use runparam, only: nt
+  use module3D_phy, only: nx,ny,H,Bsoc,ice,sealevel_2d,gr_line,uxbar,uybar,flot_marais,dtt,&
+       bm,bmelt,dt,dx,dy,isynchro,time,dtmin,ablbord_dtt,water_bilan
 
-use module3D_phy
-implicit none
+  implicit none
 
-real :: sum_H
+  real :: sum_H
 
+  real,dimension(nx,ny) :: tot_water    !< bilan d'eau
+  real,dimension(nx,ny) :: H_beau_old, diff_H_2D, diff_H_water_bilan
+  real,dimension(nx,ny) :: Bm_dtt       !< mass balance on ice points accumulated during dtt
+  real,dimension(nx,ny) :: Bmelt_dtt    !< basal melting on ice points accumulated during dtt
+  real,dimension(nx,ny) :: calv_dtt     !< calving sur dtt (pour calcul bilan d'eau)
+  real,dimension(nx,ny) :: archimtab    ! point pose si > 0
+  real,dimension(nx,ny) :: grline_dtt   ! grounding line flux during dtt
 
-real,dimension(nx,ny) :: tot_water    !< bilan d'eau
-real,dimension(nx,ny) :: H_beau_old, diff_H_2D, diff_H_water_bilan
-real,dimension(nx,ny) :: Bm_dtt       !< mass balance on ice points accumulated during dtt
-real,dimension(nx,ny) :: Bmelt_dtt    !< basal melting on ice points accumulated during dtt
-real,dimension(nx,ny) :: calv_dtt     !< calving sur dtt (pour calcul bilan d'eau)
-real,dimension(nx,ny) :: archimtab    ! point pose si > 0
-real,dimension(nx,ny) :: grline_dtt   ! grounding line flux during dtt
+  real :: sum_H_old
+  real :: diff_H
+  real :: alpha_flot                      ! ro/row
+  real :: dtt_flux   ! 1 an si dtt<1an sinon dtt_loc=dtt 
 
-real :: sum_H_old
-real :: diff_H
-real :: alpha_flot                      ! ro/row
-real :: dtt_flux   ! 1 an si dtt<1an sinon dtt_loc=dtt 
-
-real,dimension(nx,ny) :: bm_dt,bmelt_dt
+  real,dimension(nx,ny) :: bm_dt,bmelt_dt
 
 
 contains
 
-subroutine init_bilan_eau
- ! initialisation des variables
-  diff_H=0.
-  sum_H_old = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1)
-  H_beau_old(:,:)=H(:,:)
-  tot_water(:,:)=0.
-  bm_dt(:,:)=0.
-  bmelt_dt(:,:)=0.
-  alpha_flot=ro/row
-  archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel_2d(:,:)
-  gr_line(:,:)=0
-  do j=1,ny
-    do i=1,nx
-      !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then ! grounded with ice
-      if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice
-        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
-        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
-        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
-        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
-      endif
+  subroutine init_bilan_eau
+
+    integer :: i,j
+    
+    ! initialisation des variables
+    diff_H=0.
+    sum_H_old = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1)
+    H_beau_old(:,:)=H(:,:)
+    tot_water(:,:)=0.
+    bm_dt(:,:)=0.
+    bmelt_dt(:,:)=0.
+    alpha_flot=ro/row
+    archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel_2d(:,:)
+    gr_line(:,:)=0
+    do j=1,ny
+       do i=1,nx
+          !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then ! grounded with ice
+          if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice
+             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
+             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
+             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
+             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
+          endif
+       enddo
     enddo
-  enddo
-  if (dtt.ge.1.) then
-    dtt_flux=dtt
-  else
-    dtt_flux=1.
-  endif
+    if (dtt.ge.1.) then
+       dtt_flux=dtt
+    else
+       dtt_flux=1.
+    endif
 
-end subroutine init_bilan_eau
+  end subroutine init_bilan_eau
 
 
 
 
-subroutine bilan_eau
+  subroutine bilan_eau
 
-tot_water(:,:)=0.
+    integer :: i,j
 
-bm_dt(:,:)=0.
-bmelt_dt(:,:)=0.
+    tot_water(:,:)=0.
 
-bm_dt(2:nx-1,2:ny-1)=bm(2:nx-1,2:ny-1)*dt
-bmelt_dt(2:nx-1,2:ny-1)=bmelt(2:nx-1,2:ny-1)*dt
-!~ where (Bm(:,:).lt.0.)
-!~      bm_dt(:,:)=max(Bm(:,:)*dt,-H(:,:))
-!~ elsewhere
-!~      bm_dt(:,:)=bm(:,:)*dt
-!~ endwhere
+    bm_dt(:,:)=0.
+    bmelt_dt(:,:)=0.
 
-!~ where (Bmelt(:,:).lt.0.)
-!~      bmelt_dt(:,:)=max(Bmelt(:,:)*dt,-H(:,:))
-!~ elsewhere
-!~      bmelt_dt(:,:)=bmelt(:,:)*dt
-!~ endwhere
+    bm_dt(2:nx-1,2:ny-1)=bm(2:nx-1,2:ny-1)*dt
+    bmelt_dt(2:nx-1,2:ny-1)=bmelt(2:nx-1,2:ny-1)*dt
+    !~ where (Bm(:,:).lt.0.)
+    !~  bm_dt(:,:)=max(Bm(:,:)*dt,-H(:,:))
+    !~ elsewhere
+    !~  bm_dt(:,:)=bm(:,:)*dt
+    !~ endwhere
+
+    !~ where (Bmelt(:,:).lt.0.)
+    !~  bmelt_dt(:,:)=max(Bmelt(:,:)*dt,-H(:,:))
+    !~ elsewhere
+    !~  bmelt_dt(:,:)=bmelt(:,:)*dt
+    !~ endwhere
 
 
-! on fait la somme sur la grille en excluant les bords :
+    ! on fait la somme sur la grille en excluant les bords :
 
-sum_H = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1)
+    sum_H = sum(H(2:nx-1,2:ny-1),mask=ice(2:nx-1,2:ny-1)==1)
 
-diff_H = diff_H + (sum_H - sum_H_old) ! on calcul la variation de volume a chaque dt
+    diff_H = diff_H + (sum_H - sum_H_old) ! on calcul la variation de volume a chaque dt
 
-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)
+    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)
 
-where (ice(2:nx-1,2:ny-1).eq.1)
-        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
-        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
-endwhere
+    where (ice(2:nx-1,2:ny-1).eq.1)
+       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
+       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
+    endwhere
 
-archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel_2d(:,:)
-gr_line(:,:)=0
-do j=1,ny
-  do i=1,nx
-     !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then ! grounded with ice
-     if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice
-      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
-      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
-      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
-      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
+    archimtab(:,:) = Bsoc(:,:)+H(:,:)*alpha_flot - sealevel_2d(:,:)
+    gr_line(:,:)=0
+    do j=1,ny
+       do i=1,nx
+          !afq    if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.).and.(Bsoc(i,j).LE.sealevel_2d(i,j))) then ! grounded with ice
+          if ((H(i,j).gt.0.).and.(archimtab(i,j).GE.0.)) then ! grounded with ice
+             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
+             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
+             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
+             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
+          endif
+       enddo
+    enddo
+
+    where (gr_line(:,:).eq.1)
+       !~   grline_dtt(:,:)= (((uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))**2+ &
+       !~                     (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5 &
+       !~                     *H(:,:) + grline_dtt(:,:) 
+       grline_dtt(:,:)= - sqrt(                                                                        &
+            ( (uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))*dy/2. )**2 +  &
+            ( (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))*dx/2. )**2 )  &
+            * H(:,:) * dt + grline_dtt(:,:)
+    endwhere
+
+    if ((isynchro.eq.1.and.(mod(abs(TIME),1.).lt.dtmin)).or.(nt.eq.1)) then
+       ! on raisonne en bilan annuel pour simplifier :
+       !~       where (ice(2:nx-1,2:ny-1).eq.1)
+       !~               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 
+       !~               elsewhere
+       !~               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
+       !~       endwhere
+       !cdc pas besoin du test sur ice ici, il a ete fait avant (et le masque ice varie a chaque dt)    
+       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_flux
+
+       ! bilan d'eau sur la grille :
+       water_bilan=sum(tot_water(:,:))
+       diff_H = diff_H/dtt_flux
+
+       !999 format(f0.2,1x,e15.8,1x,i10,8(1x,e15.8))
+       !        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
+       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)
+
     endif
-  enddo
-enddo
+    sum_H_old=sum_H
 
-where (gr_line(:,:).eq.1)
-!~   grline_dtt(:,:)= (((uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))**2+ &
-!~                     (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))**2)**0.5)*0.5 &
-!~                     *H(:,:) + grline_dtt(:,:) 
-   grline_dtt(:,:)= - sqrt(                                                                        &
-                           ( (uxbar(:,:)+eoshift(uxbar(:,:),shift=1,boundary=0.,dim=1))*dy/2. )**2 +  &
-                           ( (uybar(:,:)+eoshift(uybar(:,:),shift=1,boundary=0.,dim=2))*dx/2. )**2 )  &
-                           * H(:,:) * dt + grline_dtt(:,:)
-endwhere
+    H_beau_old(:,:)=H(:,:)
 
-if ((isynchro.eq.1.and.(mod(abs(TIME),1.).lt.dtmin)).or.(nt.eq.1)) then
-! on raisonne en bilan annuel pour simplifier :
-!~      where (ice(2:nx-1,2:ny-1).eq.1)
-!~              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 
-!~              elsewhere
-!~              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
-!~      endwhere
-!cdc pas besoin du test sur ice ici, il a ete fait avant (et le masque ice varie a chaque dt)   
-  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_flux
-   
-! bilan d'eau sur la grille :
-  water_bilan=sum(tot_water(:,:))
-  diff_H = diff_H/dtt_flux
-
-!999 format(f0.2,1x,e15.8,1x,i10,8(1x,e15.8))
-!       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
-  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)
-
-endif
-sum_H_old=sum_H
-
-H_beau_old(:,:)=H(:,:)
-
-end subroutine bilan_eau
+  end subroutine bilan_eau
 
 end module bilan_eau_mod