source: branches/iLoveclim/SOURCES/bmelt_clio_coupl_mod.f90 @ 254

!> \file bmelt_clio_coupl_mod.f90
!! Read the bmelt computed in CLIO (thersf.f) for the 20 vertical layers
!<

!> \namespace  BMELT_CLIO_COUPL_MOD
!! Read the bmelt computed in CLIO and compute grounded basal melt
!! @note Only for coupled version
!! \author Aurelien Quiquet
!! \date October 2016
!! @note Used module
!! @note   - use module3D_phy
!<


MODULE  BMELT_CLIO_COUPL_MOD


  USE module3D_phy,only:nx,ny,i,j,ro,row,h,bsoc,hdot,bmelt,corrbmelt,igrdline,flot,fbm,bmshelfCLIO

implicit none

real, parameter       :: bmgrz_fact = 2. !< bmshelf/bmgrz ratio (fixed)

REAL,dimension(nx,ny) ::  bmgrz     !< fusion basale a la grounding zone
real,dimension(nx,ny) ::  bmshelf   !< fusion basale sous shelf

real,dimension(20)    ::  z_CLIO    !< depth of the CLIO layers (center of layer)
real,dimension(20)    ::  dz_CLIO   !< thickness of the CLIO layers


CONTAINS
!-------------------------------------------------------------------------------

!> SUBROUTINE: init_bmelt
!! This routine does the initialisation of basal melting rates of ice shelves
!<
subroutine init_bmelt

! local variables:
integer :: locdepth ! local depth of the shelf
integer :: noc      ! loop integer on the vertical oceanic layers

! This routine does the initialisation of basal melting rates of ice shelves
!  -  z_CLIO
!  - dz_CLIO
!  - bmshelf and bmgrz
! Called by initial-0.3.f90

!    ecriture dans le fichier parametres
     write(42,*)'fusion basale sous les ice shelves : bmelt from CLIO   '
     write(42,*)'-------------------------------------------------------'

!the depths of the centers of CLIO vertical layers:
     z_CLIO(20)  = 5.00
     z_CLIO(19)  = 15.98
     z_CLIO(18)  = 29.17
     z_CLIO(17)  = 45.20
     z_CLIO(16)  = 64.96
     z_CLIO(15)  = 89.75
     z_CLIO(14)  = 121.52
     z_CLIO(13)  = 163.28
     z_CLIO(12)  = 219.86
     z_CLIO(11)  = 299.26
     z_CLIO(10)  = 415.07
     z_CLIO(9)   = 588.88
     z_CLIO(8)   = 850.19
     z_CLIO(7)   = 1225.11
     z_CLIO(6)   = 1717.90
     z_CLIO(5)   = 2307.36
     z_CLIO(4)   = 2963.25
     z_CLIO(3)   = 3661.11
     z_CLIO(2)   = 4385.22
     z_CLIO(1)   = 5126.18
     
!the thicknesses of the CLIO vertical layers:
     dz_CLIO(20) = 10.00
     dz_CLIO(19) = 11.96
     dz_CLIO(18) = 14.42
     dz_CLIO(17) = 17.64
     dz_CLIO(16) = 21.88
     dz_CLIO(15) = 27.70
     dz_CLIO(14) = 35.84
     dz_CLIO(13) = 47.68
     dz_CLIO(12) = 65.48
     dz_CLIO(11) = 93.38
     dz_CLIO(10) = 138.18
     dz_CLIO(9)  = 209.44
     dz_CLIO(8)  = 313.18
     dz_CLIO(7)  = 436.66
     dz_CLIO(6)  = 548.92
     dz_CLIO(5)  = 630.00
     dz_CLIO(4)  = 681.78
     dz_CLIO(3)  = 713.94
     dz_CLIO(2)  = 734.28
     dz_CLIO(1)  = 747.64

     
     do j=1,ny
        do i=1,nx

! Init: everywhere to surface CLIO melt.
! In CLIO, when we don't have ocean, I assumed a very high bm
!          In this case, in GRISLI, we put constant values (0.2)
           if (bmshelfCLIO(i,j,20).lt.99d0) then! 
              bmshelf(i,j) = real(bmshelfCLIO(i,j,20))
              bmgrz(i,j) = bmshelf(i,j) * bmgrz_fact
           else
              bmshelf(i,j) = 0.2
              bmgrz(i,j) = bmshelf(i,j) * bmgrz_fact
           endif

           if(Bsoc(i,j).lt.-1500) then ! the melt is higher above deep ocean
              bmshelf(i,j) = 20. !bmshelf(i,j) * 10.
              bmgrz(i,j) =   20. !bmgrz(i,j)   * 10.
           endif

! Now, we look at the depth of the ice shelves to use the right bm
           !if (flot(i,j).and.(H(i,j).gt.1.)) then
           if (flot(i,j)) then
              
              bmelt(i,j) = bmshelf(i,j) ! init
              
              locdepth=ro/row*H(i,j)
              do noc=20,2,-1
                 if ( locdepth .gt. z_CLIO(noc)+dz_CLIO(noc)/2. ) then
                    if (bmshelfCLIO(i,j,noc-1).lt.99d0) then
                       bmelt(i,j) = real(bmshelfCLIO(i,j,noc-1))
                    endif
                 else
                    exit
                 endif
              enddo
              
           endif
            
        enddo
     enddo


      return
      end subroutine init_bmelt

!________________________________________________________________________________

!> SUBROUTINE: bmeltshelf
!! This routine computes the actual basal melting rates
!<
subroutine bmeltshelf

! local variables:
integer :: locdepth ! local depth of the shelf
integer :: noc      ! loop integer on the vertical oceanic layers
integer :: ngr      ! number of floating points, neighbours of a grounded points
real    :: bmsum    ! temporary bm 

do j=1,ny
   do i=1,nx

! Init: everywhere to surface CLIO melt.
! In CLIO, when we don't have ocean, I assumed a very high bm
!          In this case, in GRISLI, we put constant values (0.2)
      if (bmshelfCLIO(i,j,20).lt.99d0) then! 
         bmshelf(i,j) = real(bmshelfCLIO(i,j,20))
         bmgrz(i,j) = bmshelf(i,j) * bmgrz_fact
      else
         bmshelf(i,j) = 0.2
         bmgrz(i,j) = bmshelf(i,j) * bmgrz_fact
      endif

      if(Bsoc(i,j).lt.-1500) then ! the melt is higher above deep ocean
         bmshelf(i,j) = 20. !bmshelf(i,j) * 10.
         bmgrz(i,j) =   20. !bmgrz(i,j)   * 10.
      endif

! Now, we look at the depth of the ice shelves to use the right bm
      if (flot(i,j)) then    ! partie flottante

         bmelt(i,j) = bmshelf(i,j) ! init
         
         locdepth=ro/row*H(i,j)
         do noc=20,2,-1
            if ( locdepth.gt. z_CLIO(noc)+dz_CLIO(noc)/2. ) then
               if (bmshelfCLIO(i,j,noc-1).lt.99d0) then
                  bmelt(i,j) = real(bmshelfCLIO(i,j,noc-1))
               endif
            else
               exit
            endif
         enddo
         
         if (fbm(i,j))  bmelt(i,j)=bmgrz(i,j)


! ATTENTION LE BLOC SUIVANT SERT A FAIRE DES ICE SHELVES STATIONNAIRES
! igrdline est défini dans itnitial1
! afq -- not tested for coupled applications...
         if (igrdline.eq.1) then
            corrbmelt(i,j)=corrbmelt(i,j)+hdot(i,j)*0.8
            bmelt(i,j)=bmelt(i,j)+corrbmelt(i,j)
         endif

      endif ! on flot
   enddo
enddo

do j=2,ny-1
   do i=2,nx-1
      
      if (.not.flot(i,j)) then   ! grounded point, we account for the floating neighbours

         bmsum=0.
         ngr=0
         if (flot(i+1,j)) then
            ngr=ngr+1
            bmsum= bmsum+bmelt(i+1,j)
         endif
         if (flot(i-1,j)) then
            ngr=ngr+1
            bmsum= bmsum+bmelt(i-1,j)
         endif
         if (flot(i,j+1)) then
            ngr=ngr+1
            bmsum= bmsum+bmelt(i,j+1)
         endif
         if (flot(i,j-1)) then
            ngr=ngr+1
            bmsum= bmsum+bmelt(i,j-1)
         endif

! Grounding point basal melting rate is a combined effect of grounded and floating values:

           bmelt(i,j)= ngr/4.*bmgrz(i,j)+(1.-ngr/4.)*bmsum

        endif

     end do
  end do


    return
    end subroutine  bmeltshelf



END MODULE  BMELT_CLIO_COUPL_MOD