module auscom_ice_mod
!
!==========================================================================

! AusCOM modules:
use auscom_ice_parameters_mod

!MOM4 modules: 
use fms_mod,         only: file_exist
use fms_mod,         only: open_namelist_file, close_file, check_nml_error
use mpp_domains_mod, only: mpp_update_domains, domain2d
use mpp_mod,         only: mpp_broadcast, mpp_pe, mpp_npes
use mpp_mod,         only: FATAL, stdlog, stdout
use mpp_mod,         only: mpp_error 
use mpp_domains_mod, only: mpp_get_compute_domain, &
                           mpp_get_data_domain
use ocean_types_mod, only: ice_ocean_boundary_type, &
                           ocean_public_type
use ocean_types_mod, only: ocean_prog_tracer_type, &
                           ocean_diag_tracer_type, &
                           ocean_Thickness_type, &
                           ocean_grid_type, &
                           ocean_time_type, &
                           ocean_time_steps_type

! --- fix mom5/cice energy leaking (Fabio Dias - Russ Fiedler, 29/6/2017)
use constants_mod, only: rho0r  &  !(m^3/kg)  rho0r == 1.0/rho0
                        ,rho0   &  ! 1.035e3 kg/m^3 (rho_sw)
                        ,hlf      ! 3.34e5  J/kg

use ocean_parameters_mod,     only: rho_cp, cp_ocean ! from ocean_parameters
! - end fix.

#if defined(UNIT_TESTING)
use dump_field, only: dump_field_2d, dump_field_close
#endif

implicit none

real, dimension(:,:), allocatable :: &
          AQICE,     & ! accumulated ice form/melt flux (C*m)
          QICE         ! total column cooling from ice formation (in C*m) 
                       ! (ICEFLUX can be either time accumulated or averaged!)  

integer :: ATIME       ! accumulated time for ice formation calculation (s)  

real, dimension(:,:), allocatable :: &
          AQICE_redist,     & ! accumulated ice form/melt flux (C*m)
          QICE_redist         ! total column cooling from ice formation (in C*m) 
                       ! (ICEFLUX can be either time accumulated or averaged!)  

integer :: ATIME_redist       ! accumulated time for ice formation calculation (s)  

integer :: iisc, iiec, jjsc, jjec
!integer :: iisd, iied, jjsd, jjed, iisc, iiec, jjsc, jjec

real :: salref = 34.7           !psu, NOT as used in POP (i.e., msu)
real :: salice = 4.0            !psu, re-set as in CICE

real :: dt_ocean

contains

!============================================================================
subroutine auscom_ice_init(domain, Time_steps)

implicit none

type(domain2d) :: domain
type(ocean_time_steps_type) :: Time_steps
integer :: ioun, io_status, ierr

! provide for namelist over-ride of default values
  ioun =  open_namelist_file()
  read (ioun,auscom_ice_nml,IOSTAT=io_status)
  write (stdlog(),auscom_ice_nml)
  write (stdout(),*) 'auscom_ice_nml='
  write (stdout(),'(/)')
  write (stdout(),auscom_ice_nml)
  ierr = check_nml_error(io_status,'auscom_ice_nml')
  call close_file (ioun)

call mpp_get_compute_domain(domain,iisc,iiec,jjsc,jjec)

allocate(QICE (iisc:iiec,jjsc:jjec));     QICE(:,:)  = 0
allocate(AQICE(iisc:iiec,jjsc:jjec));     AQICE(:,:) = 0
allocate(QICE_redist (iisc:iiec,jjsc:jjec));     QICE_redist(:,:)  = 0
allocate(AQICE_redist(iisc:iiec,jjsc:jjec));     AQICE_redist(:,:) = 0

ATIME = 0
ATIME_redist = 0


dt_ocean = Time_steps%dtts

end subroutine auscom_ice_init

!============================================================================
subroutine auscom_ice_formation_new(Time,T_prog,Thickness, Frazil)
!
!calculate ice formation (including top layer and 'sub-surface' layers)
!adapted from POP (ice.F, v 1.11, dated 2003/01/28)
!

! Calculate Frazil as a diagnstic tracer. This should allow us to balance heat budgets

! Using pointers to array sections. Can use array syntax. Much clearer. POTICE now an ARRAY

implicit none

type(ocean_time_type), intent(in)           :: Time
type(ocean_prog_tracer_type), intent(inout),target :: T_prog(:)
type(ocean_Thickness_type), intent(in),target      :: Thickness
type(ocean_diag_tracer_type), intent(inout),target :: Frazil
!
!
real :: cp_over_lhfusion
       !cp_over_lhfusion = rho_sw*cp_sw/(latent_heat_fusion*rho_fw)
       !   (/deg C)        (J/m^3/deg C)      (J/kg)     (1000kg/m^3)
real :: epsilon = 1.0e-20       !as in POP

integer :: i, k, index_temp, index_salt, taup1
integer :: num_prog_tracers

real, pointer :: PTR_TEMP(:,:),PTR_SALT(:,:),PTR_THICK(:,:), PTR_FRAZIL(:,:)

real, dimension(:,:),allocatable ::  POTICE, TEMP_BEFORE

taup1=Time%taup1

! fix auscom now calculating cp_over_lhfusion:
cp_over_lhfusion = rho_cp/hlf/1000.0

num_prog_tracers = size(T_prog)
do i=1, num_prog_tracers
   if (T_prog(i)%name == 'temp') index_temp = i
   if (T_prog(i)%name == 'salt') index_salt = i
enddo

allocate(POTICE(iisc:iiec,jjsc:jjec),TEMP_BEFORE(iisc:iiec,jjsc:jjec))

!  initialize flux to zero
QICE   = 0.0

!-----------------------------------------------------------------------
! compute frazil ice formation for sub-surface layers. if ice
! forms in lower layers but layers above are warm - the heat is
! used to melt the ice. the ice formation occurs at salinity, Si.
! this volume is replaced with an equal volume at the salinity of
! the layer above. the total ice heat flux is accumulated.
!
! (???) WARNING: unless a monotone advection scheme is in place,
! advective errors could lead to temps that are far below freezing
! in some locations and this scheme will form lots of ice.
! ice formation should be limited to the top layer (kmxice=1)
! if the advection scheme is not monotone.
!-----------------------------------------------------------------------

POTICE = 0.0

do k = kmxice, 1, -1
   PTR_TEMP => T_prog(index_temp)%field(iisc:iiec,jjsc:jjec,k,taup1)
   PTR_SALT => T_prog(index_salt)%field(iisc:iiec,jjsc:jjec,k,taup1)
   PTR_THICK => Thickness%dzt(iisc:iiec,jjsc:jjec,k)
   PTR_FRAZIL => Frazil%field(iisc:iiec,jjsc:jjec,k)

#if defined(UNIT_TESTING)
    call dump_field_2d('ice_formation.input.temp', mpp_pe(), PTR_TEMP)
    call dump_field_2d('ice_formation.input.salt', mpp_pe(), PTR_SALT)
    call dump_field_2d('ice_formation.input.thickness', mpp_pe(), PTR_THICK)
    call dump_field_2d('ice_formation.input.frazil', mpp_pe(), PTR_FRAZIL)
#endif

   ! Save original temperature

   TEMP_BEFORE = PTR_TEMP


   !***
   !*** potice is the potential amount of ice formation
   !*** (potice>0) or melting (potice<0) in layer k
   !***

   POTICE = (-0.054* PTR_SALT - PTR_TEMP) * PTR_THICK

   !***
   !*** if POTICE < 0, use the heat to melt any ice from lower layers
   !*** if POTICE > 0, keep on freezing (QICE < 0)
   !***

   POTICE = max(POTICE, QICE)

   !***
   !*** adjust tracer values based on freeze/melt
   !***

   PTR_TEMP = PTR_TEMP + POTICE / PTR_THICK

   if (iceform_adj_salt) then 
      PTR_SALT = PTR_SALT + (salref - salice) * POTICE * cp_over_lhfusion / PTR_THICK
   endif

   QICE = QICE - POTICE  ! accumulate (vertically) freezing potential


 !-----------------------------------------------------------------------
 !
 ! let any residual heat in the upper layer melt previously formed ice
 !
 !-----------------------------------------------------------------------

  if ( k == 1 ) then
     AQICE = AQICE + QICE       !in degC*m

 !-----------------------------------------------------------------------
 !
 ! recalculate freezing potential based on adjusted T, S.  only interested 
 ! in melt potential now (POTICE < 0) -------- use this melt to offset any 
 ! accumulated freezing (AQICE < 0) and adjust T,S to reflect this melting 
 !
 !-----------------------------------------------------------------------

     POTICE = (-0.054* PTR_SALT - PTR_TEMP) * PTR_THICK

 !surface layer ice form[>0]/melt[<0] potential

     POTICE = max(POTICE, AQICE)
 !tricky......
     AQICE = AQICE - POTICE

 !
 ! adjust T, S again:
 !

    PTR_TEMP = PTR_TEMP + POTICE / PTR_THICK

    if (iceform_adj_salt) then 
       PTR_SALT = PTR_SALT + (salref - salice) * POTICE * cp_over_lhfusion / PTR_THICK
    endif

   

 !-----------------------------------------------------------------------
 !
 ! compute the heat flux for input to the sea ice model. 
 ! note that either:
 !   AQICE<0 and SST=Tfreeze => qflux = -AQICE >0 (net ice made), or  
 !   AQICE=0 and SST>Tfreeze => qflux ~ Tfreeze-SST <0 (melt potential)
 !
 ! REMARK: qflux ( QICE below) is needed only if it's time to communicate
 !         with oasis. however, in order to isolate this subroutine from 
 !         any coupling details, qflux is computed at every ocean time 
 !         step (but only the last calculated one is sent to oasis/ice). 
 !         Note qflux will be converted to W/m^2 later in get_o2i_fields
 !         when coupling actually happens.
 !-----------------------------------------------------------------------

    POTICE = (-0.054* PTR_SALT - PTR_TEMP) * PTR_THICK
    QICE = POTICE - AQICE 

  endif

 !-----------------------------------------------------------------------
 !  Calculate total change to Temperature and therefore heat due to Frazil
 !-----------------------------------------------------------------------

  PTR_FRAZIL = rho_cp * frazil_factor * ( PTR_TEMP - TEMP_BEFORE) * PTR_THICK 

#if defined(UNIT_TESTING)
    call dump_field_2d('ice_formation.output.frazil', mpp_pe(), PTR_FRAZIL)
    call dump_field_2d('ice_formation.output.qice', mpp_pe(), QICE)
#endif

enddo  !loop k = kmxice, 1, -1

ATIME = ATIME + dt_ocean

#if defined(UNIT_TESTING)
    call dump_field_2d('ice_formation.output.atime', mpp_pe(), &
                       reshape((/ float(ATIME) /), (/ 1, 1 /)))
#endif

deallocate(TEMP_BEFORE,POTICE)

end subroutine auscom_ice_formation_new

!============================================================================
subroutine auscom_ice_formation_new_redist(Time,T_prog,Thickness, Frazil)

! Special version for the redistributin of heat FAFMIP experiments.
! Note that this routine needs to be called before the regular version
!It is possoble to change the salt. We save the incoming salt and reset it
!before exiting.
!
!calculate ice formation (including top layer and 'sub-surface' layers)
!adapted from POP (ice.F, v 1.11, dated 2003/01/28)
!

! Calculate Frazil as a diagnstic tracer. This should allow us to balance heat budgets

! Using pointers to array sections. Can use array syntax. Much clearer. POTICE now an ARRAY

implicit none

type(ocean_time_type), intent(in)           :: Time
type(ocean_prog_tracer_type), intent(inout),target :: T_prog(:)
type(ocean_Thickness_type), intent(in),target      :: Thickness
type(ocean_diag_tracer_type), intent(inout),target :: Frazil
!
!
real :: cp_over_lhfusion
       !cp_over_lhfusion = rho_sw*cp_sw/(latent_heat_fusion*rho_fw)
       !   (/deg C)        (J/m^3/deg C)      (J/kg)     (1000kg/m^3)
real :: epsilon = 1.0e-20       !as in POP

integer :: i, k
integer :: index_temp_redist = -1
integer ::  index_salt = -1
integer :: taup1
integer :: num_prog_tracers

real, pointer :: PTR_TEMP(:,:),PTR_SALT(:,:),PTR_THICK(:,:), PTR_FRAZIL(:,:)

real, dimension(:,:),allocatable ::  POTICE, TEMP_BEFORE, SALT_BEFORE

taup1=Time%taup1

! fix auscom now calculating cp_over_lhfusion:
cp_over_lhfusion = rho_cp/hlf/1000.0

num_prog_tracers = size(T_prog)
do i=1, num_prog_tracers
   if (T_prog(i)%name == 'redist_heat') index_temp_redist = i
   if (T_prog(i)%name == 'salt') index_salt = i
enddo

if ( index_temp_redist == -1 .or. index_salt == -1 ) then
    call mpp_error(FATAL, '==>Error: auscom_ice_mod: Redistributed tracer indices not found')
endif

allocate(POTICE(iisc:iiec,jjsc:jjec),TEMP_BEFORE(iisc:iiec,jjsc:jjec))
allocate(SALT_BEFORE(iisc:iiec,jjsc:jjec))

!  initialize flux to zero
QICE_redist   = 0.0

!-----------------------------------------------------------------------
! compute frazil ice formation for sub-surface layers. if ice
! forms in lower layers but layers above are warm - the heat is
! used to melt the ice. the ice formation occurs at salinity, Si.
! this volume is replaced with an equal volume at the salinity of
! the layer above. the total ice heat flux is accumulated.
!
! (???) WARNING: unless a monotone advection scheme is in place,
! advective errors could lead to temps that are far below freezing
! in some locations and this scheme will form lots of ice.
! ice formation should be limited to the top layer (kmxice=1)
! if the advection scheme is not monotone.
!-----------------------------------------------------------------------

POTICE = 0.0

do k = kmxice, 1, -1
   PTR_TEMP => T_prog(index_temp_redist)%field(iisc:iiec,jjsc:jjec,k,taup1)
   PTR_SALT => T_prog(index_salt)%field(iisc:iiec,jjsc:jjec,k,taup1)
   PTR_THICK => Thickness%dzt(iisc:iiec,jjsc:jjec,k)
   PTR_FRAZIL => Frazil%field(iisc:iiec,jjsc:jjec,k)


   ! Save original temperature

   TEMP_BEFORE = PTR_TEMP
   SALT_BEFORE = PTR_SALT


   !***
   !*** potice is the potential amount of ice formation
   !*** (potice>0) or melting (potice<0) in layer k
   !***

   POTICE = (-0.054* PTR_SALT - PTR_TEMP) * PTR_THICK

   !***
   !*** if POTICE < 0, use the heat to melt any ice from lower layers
   !*** if POTICE > 0, keep on freezing (QICE < 0)
   !***

   POTICE = max(POTICE, QICE_redist)

   !***
   !*** adjust tracer values based on freeze/melt
   !***

   PTR_TEMP = PTR_TEMP + POTICE / PTR_THICK
   if (iceform_adj_salt) then 
      PTR_SALT = PTR_SALT + (salref - salice) * POTICE * cp_over_lhfusion / PTR_THICK
   endif

   QICE_redist = QICE_redist - POTICE  ! accumulate (vertically) freezing potential


 !-----------------------------------------------------------------------
 !
 ! let any residual heat in the upper layer melt previously formed ice
 !
 !-----------------------------------------------------------------------

  if ( k == 1 ) then
     AQICE_redist = AQICE_redist + QICE_redist       !in degC*m

 !-----------------------------------------------------------------------
 !
 ! recalculate freezing potential based on adjusted T, S.  only interested 
 ! in melt potential now (POTICE < 0) -------- use this melt to offset any 
 ! accumulated freezing (AQICE < 0) and adjust T,S to reflect this melting 
 !
 !-----------------------------------------------------------------------

     POTICE = (-0.054* PTR_SALT - PTR_TEMP) * PTR_THICK

 !surface layer ice form[>0]/melt[<0] potential

     POTICE = max(POTICE, AQICE_redist)
 !tricky......
     AQICE_redist = AQICE_redist - POTICE

 !
 ! adjust T, S again:
 !

    PTR_TEMP = PTR_TEMP + POTICE / PTR_THICK

    if (iceform_adj_salt) then 
       PTR_SALT = PTR_SALT + (salref - salice) * POTICE * cp_over_lhfusion / PTR_THICK
    endif

   

 !-----------------------------------------------------------------------
 !
 ! compute the heat flux for input to the sea ice model. 
 ! note that either:
 !   AQICE<0 and SST=Tfreeze => qflux = -AQICE >0 (net ice made), or  
 !   AQICE=0 and SST>Tfreeze => qflux ~ Tfreeze-SST <0 (melt potential)
 !
 ! REMARK: qflux ( QICE below) is needed only if it's time to communicate
 !         with oasis. however, in order to isolate this subroutine from 
 !         any coupling details, qflux is computed at every ocean time 
 !         step (but only the last calculated one is sent to oasis/ice). 
 !         Note qflux will be converted to W/m^2 later in get_o2i_fields
 !         when coupling actually happens.
 !-----------------------------------------------------------------------

    POTICE = (-0.054* PTR_SALT - PTR_TEMP) * PTR_THICK
    QICE_redist = POTICE - AQICE_redist 

  endif

 !-----------------------------------------------------------------------
 !  Calculate total change to Temperature and therefore heat due to Frazil
 !-----------------------------------------------------------------------

  PTR_FRAZIL = rho_cp * frazil_factor * ( PTR_TEMP - TEMP_BEFORE) * PTR_THICK 

! Reset salt to it original value if it has changed since we need to preserve
! its value
  if (iceform_adj_salt) then 
     PTR_SALT = SALT_BEFORE
  endif

enddo  !loop k = kmxice, 1, -1

ATIME_redist = ATIME_redist + dt_ocean

deallocate(TEMP_BEFORE,POTICE, SALT_BEFORE)

end subroutine auscom_ice_formation_new_redist

!=========================================================================
subroutine auscom_ice_heatflux_new(Ocean_sfc)
!
! convert the iceflux m*degC into W/m^2 as required by ice model and set 
! the accumulated AQICE and ATIME back to zero. 
! called once in a coupling interval before calling get_o2i_fields.
! Now also resets FAFMIP values.
! In FAFMIP runs the heat flux is based on the redistributed heat tracer.

implicit none
type (ocean_public_type) :: Ocean_sfc

integer :: i, j

! If ATIME_redist > 0 then we know that the correct heat flux for frazil should
! be based on the redistributed heat in a FAFMIP type 2 heat experiment.

if ( ATIME_redist > 0 ) then

   do j = jjsc, jjec
      do i = iisc, iiec

        Ocean_sfc%frazil(i,j) = QICE_redist(i,j) * frazil_factor * rho_cp / float(ATIME_redist)
  ! here ATIME_redist can be accumulated time (seconds) in a coupling interval 
  ! or just the ocean time step, depending how often the ice_formation_new
  ! is called  (21/07/2008, E. Hunke suggested POP calls ice_foramtion
  !                         once per coupling interval) 

      enddo
   enddo

else

   do j = jjsc, jjec
      do i = iisc, iiec

           Ocean_sfc%frazil(i,j) = QICE(i,j) * frazil_factor * rho_cp / float(ATIME)
  ! here ATIME can be accumulated time (seconds) in a coupling interval 
  ! or just the ocean time step, depending how often the ice_formation_new
  ! is called  (21/07/2008, E. Hunke suggested POP calls ice_foramtion
  !                         once per coupling interval) 

      enddo
   enddo
endif

AQICE = 0.0
ATIME = 0
AQICE_redist = 0.0
ATIME_redist = 0

#if defined(UNIT_TESTING)
    call dump_field_2d('ice_heatflux.output.frazil', mpp_pe(), Ocean_sfc%frazil)
#endif

end subroutine auscom_ice_heatflux_new

!===========================================================================
end module auscom_ice_mod
