source: branches/iLoveclim/SOURCES/bmelt-ismip6-param_mod.f90 @ 254

!> \file bmelt-ismip6-param.f90
!! bmelt computed from the non-local formulation suggested in ismip6
!<

!> \namespace bmelt-ismip6-param
!! Module for sub-shelf basal melting (grounded or ice shelves)
!! \author aquiquet
!! \date April 2019
!! @note from ismip6 suggested parametrisation
!! Should be chosen in the module_choix
!! @note Used modules
!! @note   - use module3D_phy
!! @note   - use netcdf
!! @note   - use io_netcdf_grisli
!<

module  bmelt_ismip6_param_mod

  use module3D_phy,only: nx,ny,dx,dy,ro,rofresh,row,cl,S,H,sealevel_2d,flot,bmelt,dirnameinp,num_param,num_rep_42
  ! note: the geom. (nx,ny,dx,dy) come from module_geoplace
  ! note: the densities come from param_phy_mod
  use netcdf
  use io_netcdf_grisli

  implicit none

  integer, parameter :: nzoc = 30 !< number of vertical levels in the ocean
  real*8,parameter :: cpw    = 3974.d0 !Specific heat of sea water (J/kg/K)

  real*8, dimension (nx,ny,nzoc) :: thermal_forcing !< thermal forcing            , input
  real*8, dimension (nzoc)       :: zoc             !< depth of oceanic levels    , input
  real*8, dimension (nx,ny)      :: basinNumber     !< IMBIE2 oc. basin identifier, input
  real*8, dimension (nx,ny)      :: deltaT_basin    !< deltaT                     , input
  real*8                         :: gamma0          !< gamma0                     , input
  
  real*8, dimension (nx,ny)      :: mesh_area       !< grid cell area
  real*8, dimension (nx,ny)      :: ice_draft       !< ice draft (S-H)

  real*8 :: cste !< a pre-computed constant
  integer :: Nbasin                                 !< number of oceanic basin
  
contains


  
  subroutine init_bmelt
    
    ! this routine is used to initialise the sub-shelf basal melting rate


    real*8, dimension(:),     pointer :: tab1d => null()   !< 2d array real pointer, needed for netcdf readings
    real*8, dimension(:,:),   pointer :: tab2d => null()   !< 2d array real pointer, needed for netcdf readings
    real*8, dimension(:,:,:), pointer :: tab3d => null()   !< 3d array real pointer, needed for netcdf readings
    
    character(len=100) :: file_inputs               !< files read
    character(len=100) :: file_TF                   !< files read
    character(len=100) :: file_basinNumbers         !< files read
    character(len=100) :: file_coef                 !< files read

    integer i,j !juste pour les verifs

    namelist/bmelt_ismip6_param/file_TF,file_basinNumbers,file_coef,gamma0
    
    rewind(num_param)        ! loop back at the beginning of the param_list.dat file
    read(num_param,bmelt_ismip6_param)
    write(num_rep_42,'(A)') '!  module bmelt_ismip6_param'
    write(num_rep_42,bmelt_ismip6_param)
    write(num_rep_42,'(A)') '! file_TF: 3d thermal forcing'
    write(num_rep_42,'(A)') '! file_basinNumbers: identifier for the Imbie2 basins'
    write(num_rep_42,'(A)') '! file_coef: deltaT'
    write(num_rep_42,'(A)') '! gamma0: value associated with the deltaT file'

    
    file_inputs=TRIM(DIRNAMEINP)//trim(file_TF)
    call Read_Ncdf_var('z',file_inputs,tab1d)
    zoc(:)  = tab1d(:)
    call Read_Ncdf_var('thermal_forcing',file_inputs,tab3d)
    thermal_forcing(:,:,:)  = tab3d(:,:,:)

    file_inputs=TRIM(DIRNAMEINP)//trim(file_basinNumbers)
    call Read_Ncdf_var('basinNumber',file_inputs,tab2d)
    basinNumber(:,:)  = tab2d(:,:)

    file_inputs=TRIM(DIRNAMEINP)//trim(file_coef)
    call Read_Ncdf_var('deltaT_basin',file_inputs,tab2d)
    deltaT_basin(:,:)  = tab2d(:,:)

    if ( ubound(zoc,dim=1).ne.nzoc) then
       write (*,*) "bmelt-ismip6-param: pb with the number of oceanic layers! abort..."
       STOP
    endif
    
    mesh_area(:,:) = dx*dy ! this could be corrected to account for projection deformation
    
    Nbasin=maxval(basinNumber)+1 ! number of basins
    
    cste = (row*cpw/(ro*cl))**2  ! in K^(-2)

    where (thermal_forcing(:,:,:).lt.0.d0)
       thermal_forcing(:,:,:) = 3.5
    endwhere
    
  end subroutine init_bmelt



  subroutine bmeltshelf

    ! this routine is used to compute the sub-shelf basal melting rate

    real*8, dimension(nx,ny) :: TF_draft
    real*8,allocatable,dimension(:) :: mean_TF, IS_area

    integer :: i,j,k,kinf,ksup,ngr
    real*8  :: bmloc
    
    allocate( mean_TF(Nbasin), IS_area(Nbasin)  )

    mean_TF(:) = 0.d0
    IS_area(:) = 0.d0

    ice_draft(:,:) = S(:,:)-H(:,:)-sealevel_2d(:,:)
    
    do j=1,ny
       do i=1,nx

          if (flot(i,j)) then

             if (H(i,j).gt.200.d0) then !limit on the critical thickness of ice to define the ice shelf mask
                                        ! we should use Hcalv

                ! 1 -  Linear interpolation of the thermal forcing on the ice draft depth :
                ksup=nzoc
                do k=nzoc-1,2,-1
                   if ( zoc(k) .le. ice_draft(i,j) ) ksup = k
                enddo
                kinf = ksup - 1
                if     ( ice_draft(i,j) .gt. zoc(1) ) then
                   TF_draft(i,j) = thermal_forcing(i,j,1)
                elseif ( ice_draft(i,j) .lt. zoc(nzoc) ) then
                   TF_draft(i,j) = thermal_forcing(i,j,nzoc)
                else
                   TF_draft(i,j) = (   (zoc(ksup)-ice_draft(i,j)) * thermal_forcing(i,j,kinf)   &
                        &                   + (ice_draft(i,j)-zoc(kinf)) * thermal_forcing(i,j,ksup) ) / (zoc(ksup)-zoc(kinf))
                endif

                ! 2 -  Mean Thermal forcing in individual basins (NB: fortran norm while basins start at zero) :
                mean_TF( basinNumber(i,j)+1 ) = mean_TF( basinNumber(i,j)+1 ) + mesh_area(i,j) * TF_draft(i,j)
                IS_area( basinNumber(i,j)+1 ) = IS_area( basinNumber(i,j)+1 ) + mesh_area(i,j)

             else
                TF_draft(i,j) = thermal_forcing(i,j,1)
             endif

          else
             
             TF_draft(i,j) = -9999.9d0
             
          endif

       enddo
    enddo

    where ( IS_area(:).gt.0.d0)     ! for all basins that have ice shelves
       mean_TF(:) = mean_TF(:) / IS_area(:)
    elsewhere                       ! we have no floating points over the whole basin
       mean_TF(:) = -9999.d0
    endwhere

    ! 3  - Calculation of melting rate :

    ! melt rate in m/yr (meters of pure water per year) :
    ! [ * rhofw_SI / rhoi_SI to get it in meters of ice per year ]
    do j=1,ny
       do i=1,nx
          if ( TF_draft(i,j) .gt. -5.d0 ) then !floating points
             if ((mean_TF(basinNumber(i,j)+1).gt.-9999.d0).and.(H(i,j).gt.200.)) then !should use Hcoup
                bmelt(i,j) = gamma0 * cste * ( TF_draft(i,j) + deltaT_basin(i,j) )* abs( mean_TF(basinNumber(i,j)+1) + deltaT_basin(i,j) )  
                !if (bmelt(i,j).lt.-0.2) then
                !   write(*,*) "Strong sub-shelf refreezing: ", bmelt(i,j), i,j,flot(i,j),ice_draft(i,j),TF_draft(i,j), mean_TF(basinNumber(i,j)+1),  deltaT_basin(i,j)
                !endif
             else ! in case we have no floating points over a whole basin, we take the local expression
                bmelt(i,j) = gamma0 * cste * max( TF_draft(i,j) + deltaT_basin(i,j) , 0.d0 ) * max( TF_draft(i,j) + deltaT_basin(i,j) , 0.d0 )
             endif
          endif
       enddo
    enddo
    
    do j=1,ny
       do i=1,nx
          if ( TF_draft(i,j) .le. -5.d0 ) then !grounded points
             ngr=0
             bmloc=0.d0
             if (flot(i+1,j)) then
                ngr=ngr+1
                bmloc=bmloc+bmelt(i+1,j)
             endif
             if (flot(i-1,j)) then
                ngr=ngr+1
                bmloc=bmloc+bmelt(i-1,j)
             endif
             if (flot(i,j+1)) then
                ngr=ngr+1
                bmloc=bmloc+bmelt(i,j+1)
             endif
             if (flot(i,j-1)) then
                ngr=ngr+1
                bmloc=bmloc+bmelt(i,j-1)
             endif
             bmelt(i,j)=(ngr/4.)*bmloc+(1-ngr/4.)*bmelt(i,j)
          endif
       enddo
    enddo


  end subroutine bmeltshelf
  

  
end module bmelt_ismip6_param_mod