!
! $Id: ocean_forced_mod.F90 5285 2024-10-28 13:33:29Z abarral $
!
MODULE ocean_forced_mod
!
! This module is used for both the sub-surfaces ocean and sea-ice for the case of a 
! forced ocean,  "ocean=force".
!
  IMPLICIT NONE

CONTAINS
!
!****************************************************************************************
!
  SUBROUTINE ocean_forced_noice( &
       itime, dtime, jour, knon, knindex, &
       p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, &
       temp_air, spechum, &
       AcoefH, AcoefQ, BcoefH, BcoefQ, &
       AcoefU, AcoefV, BcoefU, BcoefV, &
       ps, u1, v1, gustiness, tsurf_in, &
       radsol, snow, agesno, & 
       qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
       tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa &
#ifdef ISO
       ,xtprecip_rain, xtprecip_snow, xtspechum,Roce,rlat, &
       xtsnow,xtevap,h1 &  
#endif            
       )
!
! This subroutine treats the "open ocean", all grid points that are not entierly covered
! by ice.
! The routine receives data from climatologie file limit.nc and does some calculations at the 
! surface. 
!
    USE dimphy
    USE calcul_fluxs_mod
    USE limit_read_mod
    USE mod_grid_phy_lmdz
    USE indice_sol_mod
    USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
    use config_ocean_skin_m, only: activate_ocean_skin
#ifdef ISO
    USE infotrac_phy, ONLY: ntiso,niso
    USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall    
#ifdef ISOVERIF
    USE isotopes_mod, ONLY: iso_eau,ridicule
    !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
    USE isotopes_verif_mod
#endif
#endif
    USE clesphys_mod_h
    USE yomcst_mod_h
    INCLUDE "flux_arp.h"

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: itime, jour, knon
    INTEGER, DIMENSION(klon), INTENT(IN)     :: knindex
    REAL, INTENT(IN)                         :: dtime
    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
    REAL, DIMENSION(klon), INTENT(IN)        :: cdragh, cdragq, cdragm
    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefH, AcoefQ, BcoefH, BcoefQ
    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefU, AcoefV, BcoefU, BcoefV
    REAL, DIMENSION(klon), INTENT(IN)        :: ps
    REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)

#ifdef ISO
    REAL, DIMENSION(ntiso,klon), INTENT(IN)  :: xtprecip_rain, xtprecip_snow
    REAL, DIMENSION(ntiso,klon), INTENT(IN)  :: xtspechum
    REAL, DIMENSION(klon),       INTENT(IN)  :: rlat
#endif

! In/Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(INOUT)     :: radsol
    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno !? put to 0 in ocean
#ifdef ISO     
    REAL, DIMENSION(niso,klon), INTENT(IN)   :: xtsnow
    REAL, DIMENSION(niso,klon), INTENT(INOUT):: Roce
#endif

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
    REAL, intent(out):: sens_prec_liq(:) ! (knon)

#ifdef ISO     
    REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap ! isotopes in evaporation flux
    REAL, DIMENSION(klon),       INTENT(OUT) :: h1 ! just a diagnostic, not useful for the simulation
#endif

! Local variables
!****************************************************************************************
    INTEGER                     :: i, j
    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd
    REAL, DIMENSION(klon)       :: alb_neig, tsurf_lim, zx_sl
    REAL, DIMENSION(klon)       :: u0, v0
    REAL, DIMENSION(klon)       :: u1_lay, v1_lay
    LOGICAL                     :: check=.FALSE.
    REAL, DIMENSION(knon)       :: sens_prec_sol
    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol    
#ifdef ISO   
    REAL, PARAMETER :: t_coup = 273.15      
#endif


!****************************************************************************************
! Start calculation
!****************************************************************************************
    IF (check) WRITE(*,*)' Entering ocean_forced_noice'

#ifdef ISO
#ifdef ISOVERIF
    DO i = 1, knon
      IF (iso_eau > 0) THEN         
        CALL iso_verif_egalite_choix(xtspechum(iso_eau,i), &
     &                  spechum(i),'ocean_forced_mod 111', &
     &                  errmax,errmaxrel)     
        CALL iso_verif_egalite_choix(snow(i), &
     &                  xtsnow(iso_eau,i),'ocean_forced_mod 117', &
     &                  errmax,errmaxrel)
      ENDIF !IF (iso_eau > 0) THEN
    ENDDO !DO i=1,knon
#endif      
#endif 

!****************************************************************************************
! 1)    
! Read sea-surface temperature from file limit.nc
!
!****************************************************************************************
!--sb:
!!jyg    if (knon.eq.1) then ! single-column model
    if (klon_glo.eq.1) then ! single-column model
      ! EV: now surface Tin flux_arp.h
      !CALL read_tsurf1d(knon,tsurf_lim) ! new
       DO i = 1, knon
        tsurf_lim(i) = tg
       ENDDO

    else ! GCM
      CALL limit_read_sst(knon,knindex,tsurf_lim &
#ifdef ISO
     &     ,Roce,rlat &
#endif     
     &     )
    endif ! knon
!sb--

!****************************************************************************************
! 2)
! Flux calculation
!
!****************************************************************************************
! Set some variables for calcul_fluxs
    !cal = 0.
    !beta = 1.
    !dif_grnd = 0.
    
    
    ! EV: use calbeta to calculate beta
    ! Need to initialize qsurf for calbeta but it is not modified by this routine
    qsurf(:)=0.
    CALL calbeta(dtime, is_oce, knon, snow, qsurf, beta, cal, dif_grnd)


    alb_neig(:) = 0.
    agesno(:) = 0.
    lat_prec_liq = 0.; lat_prec_sol = 0.

! Suppose zero surface speed
    u0(:)=0.0
    v0(:)=0.0
    u1_lay(:) = u1(:) - u0(:)
    v1_lay(:) = v1(:) - v0(:)

! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf
    CALL calcul_fluxs(knon, is_oce, dtime, &
         merge(tsurf_in, tsurf_lim, activate_ocean_skin == 2), p1lay, cal, &
         beta, cdragh, cdragq, ps, &
         precip_rain, precip_snow, snow, qsurf,  &
         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
         f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
    if (activate_ocean_skin == 2) tsurf_new = tsurf_lim

    do j = 1, knon
      i = knindex(j)
      sens_prec_liq_o(i,1) = sens_prec_liq(j)
      sens_prec_sol_o(i,1) = sens_prec_sol(j)
      lat_prec_liq_o(i,1) = lat_prec_liq(j)
      lat_prec_sol_o(i,1) = lat_prec_sol(j)
    enddo


! - Flux calculation at first modele level for U and V
    CALL calcul_flux_wind(knon, dtime, &
         u0, v0, u1, v1, gustiness, cdragm, &
         AcoefU, AcoefV, BcoefU, BcoefV, &
         p1lay, temp_air, &
         flux_u1, flux_v1)  

#ifdef ISO     
    CALL calcul_iso_surf_oce_vectall(klon, knon,t_coup, &
     &    ps,tsurf_new,spechum,u1_lay, v1_lay, xtspechum, &
     &    evap, Roce,xtevap,h1 &
#ifdef ISOTRAC
     &    ,knindex &
#endif
     &    )
#endif         

#ifdef ISO
#ifdef ISOVERIF
!          write(*,*) 'ocean_forced_mod 176: sortie de ocean_forced_noice'
    IF (iso_eau > 0) THEN
      DO i = 1, knon               
        CALL iso_verif_egalite_choix(snow(i), &
     &          xtsnow(iso_eau,i),'ocean_forced_mod 180', &
     &          errmax,errmaxrel)
      ENDDO ! DO j=1,knon
    ENDIF !IF (iso_eau > 0) THEN
#endif
#endif   

  END SUBROUTINE ocean_forced_noice
!
!***************************************************************************************
!
  SUBROUTINE ocean_forced_ice( &
       itime, dtime, jour, knon, knindex, &
       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air,spechum, &
       AcoefH, AcoefQ, BcoefH, BcoefQ, &
       AcoefU, AcoefV, BcoefU, BcoefV, &
       ps, u1, v1, gustiness, &
       radsol, snow, qsol, agesno, tsoil, &
       qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
       tsurf_new, dflux_s, dflux_l, rhoa &
#ifdef ISO
       ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &
       xtsnow, xtsol,xtevap,Rland_ice &  
#endif            
       )
!
! This subroutine treats the ocean where there is ice. 
! The routine reads data from climatologie file and does flux calculations at the 
! surface.
!
    USE dimphy
    USE geometry_mod, ONLY: longitude,latitude
    USE calcul_fluxs_mod
    USE surface_data,     ONLY : calice, calsno
    USE limit_read_mod
    USE fonte_neige_mod,  ONLY : fonte_neige
    USE indice_sol_mod
    USE phys_output_var_mod, ONLY : sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o
#ifdef ISO
    USE infotrac_phy, ONLY: niso, ntiso
    USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall 
#ifdef ISOVERIF
    USE isotopes_mod, ONLY: iso_eau,ridicule
    !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix
    USE isotopes_verif_mod
#endif
#endif
    USE clesphys_mod_h
    USE yomcst_mod_h
USE dimsoil_mod_h, ONLY: nsoilmx

!   INCLUDE "indicesol.h"

    INCLUDE "flux_arp.h"

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                  :: itime, jour, knon
    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
    REAL, INTENT(IN)                     :: dtime
    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
    REAL, DIMENSION(klon), INTENT(IN)    :: ps
    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
#ifdef ISO
    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum
    REAL, DIMENSION(niso,klon),  INTENT(IN) :: Roce 
    REAL, DIMENSION(niso,klon),  INTENT(IN) :: Rland_ice
#endif

! In/Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
#ifdef ISO     
    REAL, DIMENSION(niso,klon), INTENT(INOUT)     :: xtsnow
    REAL, DIMENSION(niso,klon), INTENT(IN)        :: xtsol
#endif

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
#ifdef ISO     
    REAL, DIMENSION(ntiso,klon), INTENT(OUT)      :: xtevap
#endif      

! Local variables
!****************************************************************************************
    LOGICAL                     :: check=.FALSE.
    INTEGER                     :: i, j
    REAL                        :: zfra
    REAL, PARAMETER             :: t_grnd=271.35
    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd, capsol
    REAL, DIMENSION(klon)       :: alb_neig, tsurf_tmp
    REAL, DIMENSION(klon)       :: soilcap, soilflux
    REAL, DIMENSION(klon)       :: u0, v0
    REAL, DIMENSION(klon)       :: u1_lay, v1_lay
    REAL, DIMENSION(knon)       :: sens_prec_liq, sens_prec_sol
    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol    

#ifdef ISO
    REAL, PARAMETER :: t_coup = 273.15
    REAL, DIMENSION(klon) :: fq_fonte_diag
    REAL, DIMENSION(klon) :: fqfonte_diag
    REAL, DIMENSION(klon) :: snow_evap_diag 
    REAL, DIMENSION(klon) :: fqcalving_diag 
    REAL, DIMENSION(klon) :: run_off_lic_diag
    REAL :: coeff_rel_diag
    REAL :: max_eau_sol_diag  
    REAL, DIMENSION(klon) :: runoff_diag    
    INTEGER IXT
    REAL, DIMENSION(niso,klon) :: xtsnow_prec, xtsol_prec
    REAL, DIMENSION(klon) :: snow_prec, qsol_prec  
#endif

!****************************************************************************************
! Start calculation
!****************************************************************************************
    IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon 

!****************************************************************************************
! 1) 
! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, flux_u1, flux_v1
!                    dflux_s, dflux_l and qsurf
!****************************************************************************************

    tsurf_tmp(:) = tsurf_in(:)

! calculate the parameters cal, beta, capsol and dif_grnd and then recalculate cal
    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)

    
    IF (soil_model) THEN 
! update tsoil and calculate soilcap and soilflux
       CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, qsol, &
        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), tsoil,soilcap, soilflux)
       cal(1:knon) = RCPD / soilcap(1:knon)
       radsol(1:knon) = radsol(1:knon)  + soilflux(1:knon)
       dif_grnd = 1.0 / tau_gl
    ELSE 
       dif_grnd = 1.0 / tau_gl
       cal = RCPD * calice
       WHERE (snow > 0.0) cal = RCPD * calsno 
    ENDIF

!    beta = 1.0
    lat_prec_liq = 0.; lat_prec_sol = 0.

! Suppose zero surface speed
    u0(:)=0.0
    v0(:)=0.0
    u1_lay(:) = u1(:) - u0(:)
    v1_lay(:) = v1(:) - v0(:)
    CALL calcul_fluxs(knon, is_sic, dtime, &
         tsurf_tmp, p1lay, cal, beta, cdragh, cdragh, ps, &
         precip_rain, precip_snow, snow, qsurf,  &
         radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
         f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
         tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, &
         sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa)
    do j = 1, knon
      i = knindex(j)
      sens_prec_liq_o(i,2) = sens_prec_liq(j)
      sens_prec_sol_o(i,2) = sens_prec_sol(j)
      lat_prec_liq_o(i,2) = lat_prec_liq(j)
      lat_prec_sol_o(i,2) = lat_prec_sol(j)
    enddo

! - Flux calculation at first modele level for U and V
    CALL calcul_flux_wind(knon, dtime, &
         u0, v0, u1, v1, gustiness, cdragm, &
         AcoefU, AcoefV, BcoefU, BcoefV, &
         p1lay, temp_air, &
         flux_u1, flux_v1)  

!****************************************************************************************
! 2)
! Calculations due to snow and runoff
!
!****************************************************************************************
#ifdef ISO
   ! verif
#ifdef ISOVERIF
    DO i = 1, knon
      IF (iso_eau > 0) THEN
        IF (snow(i) > ridicule) THEN
          CALL iso_verif_egalite_choix(xtsnow(iso_eau,i),snow(i), &
   &              'interfsurf 964',errmax,errmaxrel)
        ENDIF !IF ((snow(i) > ridicule)) THEN
      ENDIF !IF (iso_eau > 0) THEN      
    ENDDO !DO i=1,knon  
#endif
   ! end verif

    DO i = 1, knon
      snow_prec(i) = snow(i)
      DO ixt = 1, niso
      xtsnow_prec(ixt,i) = xtsnow(ixt,i)
      ENDDO !DO ixt=1,niso
      ! initialisation:
      fq_fonte_diag(i) = 0.0
      fqfonte_diag(i)  = 0.0
      snow_evap_diag(i)= 0.0
    ENDDO !DO i=1,knon
#endif


    CALL fonte_neige( knon, is_sic, knindex, dtime, &
         tsurf_tmp, precip_rain, precip_snow, &
         snow, qsol, tsurf_new, evap &
#ifdef ISO    
     &  ,fq_fonte_diag,fqfonte_diag,snow_evap_diag,fqcalving_diag   &
     &  ,max_eau_sol_diag,runoff_diag,run_off_lic_diag,coeff_rel_diag   &
#endif
     &   )


#ifdef ISO
! isotopes: tout est externalisé
!#ifdef ISOVERIF
!        write(*,*) 'ocean_forced_mod 377: call calcul_iso_surf_sic_vectall'
!        write(*,*) 'klon,knon=',klon,knon
!#endif
    CALL calcul_iso_surf_sic_vectall(klon,knon, &
     &          evap,snow_evap_diag,Tsurf_new,Roce,snow, &
     &          fq_fonte_diag,fqfonte_diag,dtime,t_coup, &
     &          precip_snow,xtprecip_snow,xtprecip_rain, snow_prec,xtsnow_prec, &
     &          xtspechum,spechum,ps, &
     &          xtevap,xtsnow,fqcalving_diag, &
     &          knindex,is_sic,run_off_lic_diag,coeff_rel_diag,Rland_ice &
     &   )
#ifdef ISOVERIF
        !write(*,*) 'ocean_forced_mod 391: sortie calcul_iso_surf_sic_vectall'
    IF (iso_eau > 0) THEN
      DO i = 1, knon  
        CALL iso_verif_egalite_choix(snow(i), &
     &           xtsnow(iso_eau,i),'ocean_forced_mod 396', &
     &           errmax,errmaxrel)
      ENDDO ! DO j=1,knon
    ENDIF !IF (iso_eau > 0) then
#endif 
!#ifdef ISOVERIF
#endif   
!#ifdef ISO
    
! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
! 
    CALL albsno(klon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:))  

    WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0.

    alb1_new(:) = 0.0
    DO i=1, knon
       zfra = MAX(0.0,MIN(1.0,snow(i)/(snow(i)+10.0)))
       alb1_new(i) = alb_neig(i) * zfra +  0.6 * (1.0-zfra)
    ENDDO

    alb2_new(:) = alb1_new(:)

  END SUBROUTINE ocean_forced_ice

!************************************************************************
! 1D case
!************************************************************************
!  SUBROUTINE read_tsurf1d(knon,sst_out)
!
! This subroutine specifies the surface temperature to be used in 1D simulations
!
!      USE dimphy, ONLY : klon
!
!      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
!      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
!
!       INTEGER :: i
! COMMON defined in lmdz1d.F:
!       real ts_cur
!       common /sst_forcing/ts_cur
!
!       DO i = 1, knon
!        sst_out(i) = ts_cur
!       ENDDO
!
!      END SUBROUTINE read_tsurf1d
!
!
!************************************************************************
END MODULE ocean_forced_mod