MODULE bulk

  use com_cpl
  use ncdf
  use store

  implicit none

  private

  public :: bulk_dqnsdt

  integer, parameter :: jpi=362  !! ORCA1 x-points
  integer, parameter :: jpj=292  !! ORCA1 y=points
  integer, parameter :: jpl=1    !! ice categories
  integer, parameter :: wp=kind(0.d0)

  !--- CORE bulk parameters
  real(wp), parameter ::   rhoa =    1.22        ! air density
  real(wp), parameter ::   cpa  = 1000.5         ! specific heat of air
  real(wp), parameter ::   Lv   =    2.5e6       ! latent heat of vaporization
  real(wp), parameter ::   Ls   =    2.839e6     ! latent heat of sublimation
  real(wp), parameter ::   Stef =    5.67e-8     ! Stefan Boltzmann constant
  real(wp), parameter ::   Cice =    1.63e-3     ! transfer coefficient over ice

CONTAINS

    subroutine bulk_dqnsdt(dqns, dqla, dqlw, dqsb)
      real(kind=8), pointer :: dqns(:), dqla(:), dqlw(:), dqsb(:)

      !--- Local
      real(wp), dimension(jpi,jpj,jpl)   :: pst      ! ice temp [K]
      real(wp), dimension(jpi,jpj)       :: z_wnds_t ! wind speed at T point [m/s]
      real(wp), dimension(jpi, jpj, jpl) :: p_dqns   ! non solar heat sensistivity  (T-point) [W/m2/K]
      real(wp), dimension(jpi, jpj, jpl) :: p_dqla   ! latent    heat sensistivity  (T-point) [W/m2/K]
      real(wp), dimension(jpi, jpj, jpl) :: p_dqlw   ! long wave heat sensitivity over ice
      real(wp), dimension(jpi, jpj, jpl) :: p_dqsb   ! sensible  heat sensitivity over ice

      integer :: nwrds, sz
      real(kind=8), pointer :: wrk1(:), wrk2(:)
      integer(kind=8) :: store_time, cpl_time

      nullify(wrk1,wrk2)

      if ( jpi /= nlon_o+olap_o .or. jpj /= nlat_o ) then
        write(6,*)"bulk_dqnsdt: Dimension size mismatch."
        write(6,*)" jpi,jpj = ",jpi,jpj
        write(6,*)" expecting ",nlon_o+olap_o,nlat_o
        call xit("bulk_dqnsdt",-1)
      endif

      !--- Read ice temperature from storage
      !--- If ice temperature is sent from NEMO it will be one coupling time step
      !--- behind the current time step
      cpl_time = cpl_elapsed_time_secs - cpl_time_step_secs
      store_time = max(0_8, cpl_time)
      nwrds = store_readRecord("O_TepIce", ocn_gid, wrk1, time=store_time)
      if ( nwrds < 0 ) then
        !--- The field was not found in storage
        write(6,*)"prepare_nemo_fields: O_TepIce is missing from storage."
        call xit("bulk_dqnsdt",-1)
      endif
      pst(1:jpi,1:jpj,1) = reshape( wrk1, (/ jpi,jpj /) )

where (pst<10.0) pst=10.0

      !--- Read wind speed from storage
      !--- Read O_Wind10 from stroage with the same time as O_TepIce
      cpl_time = cpl_elapsed_time_secs - cpl_time_step_secs
      store_time = max(0_8, cpl_time)
      nwrds = store_readRecord("O_Wind10", atm_gid, wrk1, time=store_time)
      if ( nwrds < 0 ) then
        !--- The field was not found in storage
        write(6,*)"prepare_nemo_fields: O_Wind10 is missing from storage."
        call xit("bulk_dqnsdt",-1)
      endif
      z_wnds_t(1:jpi,1:jpj) = reshape( wrk1, (/ jpi,jpj /) )

      call blk_dqnsdt(pst, z_wnds_t, p_dqns, p_dqla, p_dqlw, p_dqsb)

      sz = size(p_dqns)

      if ( associated(dqns) ) deallocate(dqns)
      allocate( dqns(sz) )
      dqns(1:sz) = reshape( p_dqns, (/ sz /) )

      if ( associated(dqla) ) deallocate(dqla)
      allocate( dqla(sz) )
      dqla(1:sz) = reshape( p_dqla, (/ sz /) )

      if ( associated(dqlw) ) deallocate(dqlw)
      allocate( dqlw(sz) )
      dqlw(1:sz) = reshape( p_dqlw, (/ sz /) )

      if ( associated(dqsb) ) deallocate(dqsb)
      allocate( dqsb(sz) )
      dqsb(1:sz) = reshape( p_dqsb, (/ sz /) )

    end subroutine bulk_dqnsdt

    subroutine read_dqnsdt_orca_from_file(dqnsdt)
      real(kind=8), pointer :: dqnsdt(:)

      !--- Local
      integer(4) :: fid, vid, ncstat, nwrds
      character(128) :: fname
      !--- dQnsdT on the ORCA1 grid
      real(kind=8), pointer :: dqnsdt_tmp(:,:)

      if ( store_GetIndex("dQnsdT", other_gid) > 0 ) then
        !--- Read dQnsdT from storage
        nwrds = store_readRecord("dQnsdT", other_gid, dqnsdt, idx=1)

      else
        !--- Read Dqnsdt_orca from and external file

        !--- This file contains 3 hourly data for dQnsdT
        fname = "/fs/dev/crb/had_data/data/dqns_from_mc_odp_3h_00010701_00010731.nc.001"
        call ncdf_open(fid, fname, "r")

        !--- Read Dqnsdt_orca from the file
        call ncdf_read_var(fid, "dqns", dqnsdt_tmp, nt=1)

        call ncdf_close(fid)

        if ( associated(dqnsdt) ) deallocate(dqnsdt)
        allocate( dqnsdt(size(dqnsdt_tmp)) )
        dqnsdt = reshape(dqnsdt_tmp, (/size(dqnsdt_tmp)/))

        !--- Add this to storage
        call store_AppendRecord("dQnsdT", other_gid, dqnsdt, create=.true.)
      endif

    end subroutine read_dqnsdt_orca_from_file

    SUBROUTINE blk_dqnsdt(pst, z_wnds_t, p_dqns, p_dqla, p_dqlw, p_dqsb)
      !!---------------------------------------------------------------------
      !!                     ***  ROUTINE blk_dqnsdt  ***
      !!
      !! ** Purpose :   provide dQns/dt over sea-ice
      !!
      !! ** Method  :   from CORE bulk formula sbcblk_core.F90
      !!---------------------------------------------------------------------
      !! INPUT !!
      ! ice surface temperature (>0, =rt0 over land) [Kelvin]
      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pst
      ! wind speed at T-point
      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   z_wnds_t

      !! OUTPUT !!
      ! non solar heat sensistivity  (T-point) [W/m2/K]
      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqns
      ! latent    heat sensistivity  (T-point) [W/m2/K]
      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqla
      ! long wave heat sensitivity over ice
      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqlw
      ! sensible  heat sensitivity over ice
      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_dqsb

      REAL(wp) :: zst2, zst3                         ! temporary temps
      REAL(wp) :: zcoef_dqlw, zcoef_dqla, zcoef_dqsb ! coefficients
      INTEGER  :: jl, jj, ji


      zcoef_dqlw   = 4.0 * 0.95 * Stef
      zcoef_dqla   = -Ls * Cice * 11637800. * (-5897.8)
      zcoef_dqsb   = rhoa * cpa * Cice

      DO jl = 1, jpl                       !  Loop over ice categories  !
         !                                 ! ========================== !
         DO jj = 1 , jpj
            DO ji = 1, jpi
               zst2 = pst(ji,jj,jl) * pst(ji,jj,jl)
               zst3 = pst(ji,jj,jl) * zst2
               ! lw sensitivity
               p_dqlw(ji,jj,jl) = zcoef_dqlw * zst3

               ! Latent heat sensitivity for ice (Dqla/Dt)
               p_dqla(ji,jj,jl) = zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) )

               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
               p_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj)

               ! Total non solar heat flux sensitivity for ice
               p_dqns(ji,jj,jl) = - ( p_dqlw(ji,jj,jl) + p_dqsb(ji,jj,jl) + p_dqla(ji,jj,jl) )
            END DO
            !
         END DO
         !
      END DO

    END SUBROUTINE blk_dqnsdt

END MODULE bulk
