Skip to content
Snippets Groups Projects
sbccpl.F90 188 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
      !! ** Purpose :   provide the heat and freshwater fluxes of the ocean-ice system
      !!
      !! ** Method  :   transform the fields received from the atmosphere into
      !!             surface heat and fresh water boundary condition for the
      !!             ice-ocean system. The following fields are provided:
      !!               * total non solar, solar and freshwater fluxes (qns_tot,
      !!             qsr_tot and emp_tot) (total means weighted ice-ocean flux)
      !!             NB: emp_tot include runoffs and calving.
      !!               * fluxes over ice (qns_ice, qsr_ice, emp_ice) where
      !!             emp_ice = sublimation - solid precipitation as liquid
      !!             precipitation are re-routed directly to the ocean and
      !!             calving directly enter the ocean (runoffs are read but included in trasbc.F90)
      !!               * solid precipitation (sprecip), used to add to qns_tot
      !!             the heat lost associated to melting solid precipitation
      !!             over the ocean fraction.
      !!               * heat content of rain, snow and evap can also be provided,
      !!             otherwise heat flux associated with these mass flux are
      !!             guessed (qemp_oce, qemp_ice)
      !!
      !!             - the fluxes have been separated from the stress as
      !!               (a) they are updated at each ice time step compare to
      !!               an update at each coupled time step for the stress, and
      !!               (b) the conservative computation of the fluxes over the
      !!               sea-ice area requires the knowledge of the ice fraction
      !!               after the ice advection and before the ice thermodynamics,
      !!               so that the stress is updated before the ice dynamics
      !!               while the fluxes are updated after it.
      !!
      !! ** Details
      !!             qns_tot = (1-a) * qns_oce + a * qns_ice               => provided
      !!                     + qemp_oce + qemp_ice                         => recalculated and added up to qns
      !!
      !!             qsr_tot = (1-a) * qsr_oce + a * qsr_ice               => provided
      !!
      !!             emp_tot = emp_oce + emp_ice                           => calving is provided and added to emp_tot (and emp_oce).
      !!                                                                      runoff (which includes rivers+icebergs) and iceshelf
      !!                                                                      are provided but not included in emp here. Only runoff will
      !!                                                                      be included in emp in other parts of NEMO code
      !!
      !! ** Note : In case of the ice-atm coupling with conduction fluxes (such as Jules interface for the Met-Office),
      !!              qsr_ice and qns_ice are not provided and they are not supposed to be used in the ice code.
      !!              However, by precaution we also "fake" qns_ice and qsr_ice this way:
      !!              qns_ice = qml_ice + qcn_ice ??
      !!              qsr_ice = qtr_ice_top ??
      !!
      !! ** Action  :   update at each nf_ice time step:
      !!                   qns_tot, qsr_tot  non-solar and solar total heat fluxes
      !!                   qns_ice, qsr_ice  non-solar and solar heat fluxes over the ice
      !!                   emp_tot           total evaporation - precipitation(liquid and solid) (-calving)
      !!                   emp_ice           ice sublimation - solid precipitation over the ice
      !!                   dqns_ice          d(non-solar heat flux)/d(Temperature) over the ice
      !!                   sprecip           solid precipitation over the ocean
      !!----------------------------------------------------------------------
      INTEGER,  INTENT(in)                                ::   kt         ! ocean model time step index (only for a_i_last_couple)
      REAL(wp), INTENT(in)   , DIMENSION(:,:)             ::   picefr     ! ice fraction                [0 to 1]
      !                                                   !!           ! optional arguments, used only in 'mixed oce-ice' case or for Met-Office coupling
      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   palbi      ! all skies ice albedo
      REAL(wp), INTENT(in)   , DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius]
      REAL(wp), INTENT(inout), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] => inout for Met-Office
      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m]
      REAL(wp), INTENT(in)   , DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m]
      !
      INTEGER  ::   ji, jj, jl   ! dummy loop index
      REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw
      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice
      REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice
      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap_ice_total
      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice, zqtr_ice_top, ztsu
      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
      !!----------------------------------------------------------------------
      !
#if defined key_si3 || defined key_cice
      !
      IF( kt == nit000 ) THEN
         ! allocate ice fractions from last coupling time here and not in sbc_cpl_init because of jpl
         IF( .NOT.ALLOCATED(a_i_last_couple) )   ALLOCATE( a_i_last_couple(jpi,jpj,jpl) )
         ! initialize to a_i for the 1st time step
         a_i_last_couple(:,:,:) = a_i(:,:,:)
      ENDIF
      !
      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0)
      ziceld(:,:) = 1._wp - picefr(:,:)
      zcptn (:,:) = rcp * sst_m(:,:)


      !                                                      ! ========================= !
      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
      !                                                      ! ========================= !
      CASE ('coupled')
         IF (ln_scale_ice_flux) THEN
            WHERE( a_i(:,:,:) > 1.e-10_wp )
               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
            ELSEWHERE
               qml_ice(:,:,:) = 0.0_wp
               qcn_ice(:,:,:) = 0.0_wp
            END WHERE
         ELSE
            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
         ENDIF
      END SELECT


Guillaume Samson's avatar
Guillaume Samson committed
      !
      !                                                      ! ========================= !
      !                                                      !    freshwater budget      !   (emp_tot)
      !                                                      ! ========================= !
      !
      !                                                           ! solid Precipitation                                (sprecip)
      !                                                           ! liquid + solid Precipitation                       (tprecip)
      !                                                           ! total Evaporation - total Precipitation            (emp_tot)
      !                                                           ! sublimation - solid precipitation (cell average)   (emp_ice)
      SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
      CASE( 'conservative' )   ! received fields: jpr_rain, jpr_snow, jpr_ievp, jpr_tevp
         zsprecip(:,:) =   frcv(jpr_snow)%z3(:,:,1)                  ! May need to ensure positive here
         ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here
         IF (.not. ln_couple_ocean_evap ) THEN
            zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:)
         END IF
Guillaume Samson's avatar
Guillaume Samson committed
      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp
         zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1)
         zemp_ice(:,:) = frcv(jpr_semp)%z3(:,:,1) * picefr(:,:)
         zsprecip(:,:) = frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_semp)%z3(:,:,1)
         ztprecip(:,:) = frcv(jpr_semp)%z3(:,:,1) - frcv(jpr_sbpr)%z3(:,:,1) + zsprecip(:,:)
      CASE( 'none'      )       ! Not available as for now: needs additional coding below when computing zevap_oce
      !                         ! since fields received are not defined with none option
         CALL ctl_stop( 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl' )
      CASE default              ! Default
         CALL ctl_stop( 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_emp value in namelist namsbc_cpl' )
Guillaume Samson's avatar
Guillaume Samson committed
      END SELECT

      ! --- evaporation over ice (kg/m2/s) --- !
      IF (ln_scale_ice_flux) THEN ! typically met-office requirements
         IF (sn_rcv_emp%clcat == 'yes') THEN
            WHERE( a_i(:,:,:) > 1.e-10 )  ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
            ELSEWHERE                     ; zevap_ice(:,:,:) = 0._wp
            END WHERE
            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
            END WHERE
         ELSE
            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * SUM( a_i_last_couple, dim=3 ) / picefr(:,:)
            ELSEWHERE                     ; zevap_ice(:,:,1) = 0._wp
            END WHERE
            zevap_ice_total(:,:) = zevap_ice(:,:,1)
            DO jl = 2, jpl
               zevap_ice(:,:,jl) = zevap_ice(:,:,1)
            ENDDO
         ENDIF
      ELSE
         IF (sn_rcv_emp%clcat == 'yes') THEN
            zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl)
            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:)
            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp
            END WHERE
         ELSE
            zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1)
            zevap_ice_total(:,:) = zevap_ice(:,:,1)
            DO jl = 2, jpl
               zevap_ice(:,:,jl) = zevap_ice(:,:,1)
            ENDDO
         ENDIF
      ENDIF

      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN
         ! For conservative case zemp_ice has not been defined yet. Do it now.
         zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:)
      ENDIF

      ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing)
      zsnw(:,:) = 0._wp   ;   CALL ice_var_snwblow( ziceld, zsnw )

      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- !
      zemp_ice(:,:) = zemp_ice(:,:) + zsprecip(:,:) * ( picefr(:,:) - zsnw(:,:) )  ! emp_ice = A * sublimation - zsnw * sprecip

      IF ( ln_couple_ocean_evap ) THEN
         zemp_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_rain)%z3(:,:,1) & !Ocean evap minus rain (as all rain goes straight to ocean in GC5)
                        & - zsprecip(:,:) * ( 1._wp - zsnw(:,:) )        !subtract snow in leads after correction for blowing snow
         zemp_tot(:,:) = zemp_oce(:,:) + zemp_ice(:,:)
         zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1)
      ELSE
         zemp_oce(:,:) = zemp_tot(:,:) - zemp_ice(:,:)                                ! emp_oce = emp_tot - emp_ice

         ! --- evaporation over ocean (used later for qemp) --- !
         zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:)
      END IF
Guillaume Samson's avatar
Guillaume Samson committed

      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0
      ! therefore, sublimation is not redistributed over the ice categories when no subgrid scale fluxes are provided by atm.
      zdevap_ice(:,:) = 0._wp

      ! --- Continental fluxes --- !
      IF( srcv(jpr_rnf)%laction ) THEN   ! 2D runoffs (included in emp later on)
Guillaume Samson's avatar
Guillaume Samson committed
         rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
      ENDIF
      IF( srcv(jpr_rnf_1d)%laction ) THEN ! 1D runoff
         CALL cpl_rnf_1d_to_2d(frcv(jpr_rnf_1d)%z3(:,:,:))
      ENDIF
Guillaume Samson's avatar
Guillaume Samson committed
      IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot and emp_oce)
         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
         zemp_oce(:,:) = zemp_oce(:,:) - frcv(jpr_cal)%z3(:,:,1)
      ENDIF
      IF( srcv(jpr_icb)%laction ) THEN   ! iceberg added to runoffs
         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
         rnf(:,:)    = rnf(:,:) + fwficb(:,:)
      ENDIF
      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf > 0 mean melting)
        fwfisf_oasis(:,:) = frcv(jpr_isf)%z3(:,:,1)
      ENDIF

      IF( ln_mixcpl ) THEN
         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
         emp_oce(:,:) = emp_oce(:,:) * xcplmask(:,:,0) + zemp_oce(:,:) * zmsk(:,:)
         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
         DO jl = 1, jpl
            evap_ice (:,:,jl) = evap_ice (:,:,jl) * xcplmask(:,:,0) + zevap_ice (:,:,jl) * zmsk(:,:)
            devap_ice(:,:,jl) = devap_ice(:,:,jl) * xcplmask(:,:,0) + zdevap_ice(:,:)    * zmsk(:,:)
         END DO
      ELSE
         emp_tot (:,:)   = zemp_tot (:,:)
         emp_ice (:,:)   = zemp_ice (:,:)
         emp_oce (:,:)   = zemp_oce (:,:)
         sprecip (:,:)   = zsprecip (:,:)
         tprecip (:,:)   = ztprecip (:,:)
         evap_ice(:,:,:) = zevap_ice(:,:,:)
         DO jl = 1, jpl
            devap_ice(:,:,jl) = zdevap_ice(:,:)
         END DO
      ENDIF

!! for CICE ??
!!$      zsnw(:,:) = picefr(:,:)
!!$      ! --- Continental fluxes --- !
!!$      IF( srcv(jpr_rnf)%laction ) THEN   ! runoffs (included in emp later on)
!!$         rnf(:,:) = frcv(jpr_rnf)%z3(:,:,1)
!!$      ENDIF
!!$      IF( srcv(jpr_cal)%laction ) THEN   ! calving (put in emp_tot)
!!$         zemp_tot(:,:) = zemp_tot(:,:) - frcv(jpr_cal)%z3(:,:,1)
!!$      ENDIF
!!$      IF( srcv(jpr_icb)%laction ) THEN   ! iceberg added to runoffs
!!$         fwficb(:,:) = frcv(jpr_icb)%z3(:,:,1)
!!$         rnf(:,:)    = rnf(:,:) + fwficb(:,:)
!!$      ENDIF
!!$      IF( srcv(jpr_isf)%laction ) THEN   ! iceshelf (fwfisf >0 mean melting)
!!$        fwfisf_oasis(:,:) = frcv(jpr_isf)%z3(:,:,1)
!!$      ENDIF
!!$      !
!!$      IF( ln_mixcpl ) THEN
!!$         emp_tot(:,:) = emp_tot(:,:) * xcplmask(:,:,0) + zemp_tot(:,:) * zmsk(:,:)
!!$         emp_ice(:,:) = emp_ice(:,:) * xcplmask(:,:,0) + zemp_ice(:,:) * zmsk(:,:)
!!$         sprecip(:,:) = sprecip(:,:) * xcplmask(:,:,0) + zsprecip(:,:) * zmsk(:,:)
!!$         tprecip(:,:) = tprecip(:,:) * xcplmask(:,:,0) + ztprecip(:,:) * zmsk(:,:)
!!$      ELSE
!!$         emp_tot(:,:) =                                  zemp_tot(:,:)
!!$         emp_ice(:,:) =                                  zemp_ice(:,:)
!!$         sprecip(:,:) =                                  zsprecip(:,:)
!!$         tprecip(:,:) =                                  ztprecip(:,:)
!!$      ENDIF
      !
      ! outputs
      IF( srcv(jpr_cal)%laction )    CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving
      IF( srcv(jpr_icb)%laction )    CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs
      IF( iom_use('snowpre') )       CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow
      IF( iom_use('precip') )        CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation
      IF( iom_use('rain') )          CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation
      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average)
      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average)
      IF( iom_use('rain_ao_cea') )   CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * ziceld(:,:)         )  ! liquid precipitation over ocean (cell average)
      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea' , zevap_ice_total(:,:) * picefr(:,:) * tmask(:,:,1)     )  ! Sublimation over sea-ice (cell average)
      IF( iom_use('evap_ao_cea') )   CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  &
         &                                                         - zevap_ice_total(:,:) * picefr(:,:) ) * tmask(:,:,1) )  ! ice-free oce evap (cell average)
      ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf
!!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff
!!      IF( srcv(jpr_isf)%laction )   CALL iom_put( 'iceshelf_cea', fwfisf(:,:) * tmask(:,:,1)                         )  ! iceshelf
      !
      !                                                      ! ========================= !
      SELECT CASE( TRIM( sn_rcv_iceflx%cldes ) )             !  ice topmelt and botmelt  !
      !                                                      ! ========================= !
      CASE ('coupled')
         IF (ln_scale_ice_flux) THEN
            WHERE( a_i(:,:,:) > 1.e-10_wp )
               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:)
            ELSEWHERE
               qml_ice(:,:,:) = 0.0_wp
               qcn_ice(:,:,:) = 0.0_wp
            END WHERE
         ELSE
            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:)
            qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:)
         ENDIF
      END SELECT
      !
      !                                                      ! ========================= !
      SELECT CASE( TRIM( sn_rcv_qns%cldes ) )                !   non solar heat fluxes   !   (qns)
      !                                                      ! ========================= !
      CASE( 'oce only' )         ! the required field is directly provided
         ! Get the sea ice non solar heat flux from conductive, melting and sublimation fluxes
         IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN
            zqns_ice(:,:,:) = qml_ice(:,:,:) + qcn_ice(:,:,:)
         ELSE
            zqns_ice(:,:,:) = 0._wp
         ENDIF
         ! Calculate the total non solar heat flux. The ocean only non solar heat flux (zqns_oce) will be recalculated after this CASE
         ! statement to be consistent with other coupling methods even though .zqns_oce = frcv(jpr_qnsoce)%z3(:,:,1)
         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1) + SUM( zqns_ice(:,:,:) * a_i(:,:,:), dim=3 )
      CASE( 'conservative' )     ! the required fields are directly provided
         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1)
         IF( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
            zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl)
         ELSE
            DO jl = 1, jpl
               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) ! Set all category values equal
            END DO
         ENDIF
      CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes
         zqns_tot(:,:) =  ziceld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1)
         IF( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN
            DO jl=1,jpl
               zqns_tot(:,:   ) = zqns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)
               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,jl)
            ENDDO
         ELSE
            zqns_tot(:,:) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1)
            DO jl = 1, jpl
               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1)
            END DO
         ENDIF
      CASE( 'mixed oce-ice' )    ! the ice flux is cumputed from the total flux, the SST and ice informations
! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
         zqns_tot(:,:  ) = frcv(jpr_qnsmix)%z3(:,:,1)
         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
            DO jl = 1, jpl
               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:,jl)    &
                  &             + frcv(jpr_dqnsdt)%z3(:,:,jl) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
                  &                                             + pist(:,:,jl) * picefr(:,:) ) )
            END DO
         ELSE
            DO jl = 1, jpl
               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:, 1)    &
                  &             + frcv(jpr_dqnsdt)%z3(:,:, 1) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   &
                  &                                             + pist(:,:,jl) * picefr(:,:) ) )
            END DO
         ENDIF
      END SELECT
      !
      ! --- calving (removed from qns_tot) --- !
      IF( srcv(jpr_cal)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_cal)%z3(:,:,1) * rLfus  ! remove latent heat of calving
                                                                                                     ! we suppose it melts at 0deg, though it should be temp. of surrounding ocean
      ! --- iceberg (removed from qns_tot) --- !
      IF( srcv(jpr_icb)%laction )   zqns_tot(:,:) = zqns_tot(:,:) - frcv(jpr_icb)%z3(:,:,1) * rLfus  ! remove latent heat of iceberg melting

      ! --- non solar flux over ocean --- !
      !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
      zqns_oce = 0._wp



      WHERE( ziceld /= 0._wp )   zqns_oce(:,:) = ( zqns_tot(:,:) - SUM( a_i * zqns_ice, dim=3 ) ) / ziceld(:,:)
Guillaume Samson's avatar
Guillaume Samson committed
      ! Heat content per unit mass of snow (J/kg)
      WHERE( SUM( a_i, dim=3 ) > 1.e-10 )   ;   zcptsnw(:,:) = rcpi * SUM( (tn_ice - rt0) * a_i, dim=3 ) / SUM( a_i, dim=3 )
      ELSEWHERE                             ;   zcptsnw(:,:) = zcptn(:,:)
      ENDWHERE
      ! Heat content per unit mass of rain (J/kg)
      zcptrain(:,:) = rcp * ( SUM( (tn_ice(:,:,:) - rt0) * a_i(:,:,:), dim=3 ) + sst_m(:,:) * ziceld(:,:) )

      ! --- enthalpy of snow precip over ice in J/m3 (to be used in 1D-thermo) --- !
      zqprec_ice(:,:) = rhos * ( zcptsnw(:,:) - rLfus )

      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- !
      DO jl = 1, jpl
         ! zqevap_ice(:,:,jl) = 0._wp ! should be -evap * ( ( Tice - rt0 ) * rcpi ) but atm. does not take it into account

         ! How much enthalpy is stored in sublimating snow and ice. At this stage we don't know if it is snow or ice that is
         ! sublimating so we will use the combined snow and ice layer temperature t1_ice.
         zqevap_ice(:,:,jl) = -zevap_ice(:,:,jl) * ( ( rt0 - t1_ice(:,:,jl) ) * rcpi + rLfus )

Guillaume Samson's avatar
Guillaume Samson committed
      END DO

      ! --- heat flux associated with emp (W/m2) --- !
      zqemp_oce(:,:) = -  zevap_oce(:,:)                                      *   zcptn   (:,:)   &        ! evap
         &             + ( ztprecip(:,:) - zsprecip(:,:) )                    *   zcptrain(:,:)   &        ! liquid precip
         &             +   zsprecip(:,:)                   * ( 1._wp - zsnw ) * ( zcptsnw (:,:) - rLfus )  ! solid precip over ocean + snow melting
      zqemp_ice(:,:) =     zsprecip(:,:)                   * zsnw             * ( zcptsnw (:,:) - rLfus )  ! solid precip over ice (qevap_ice=0 since atm. does not take it into account)
!!    zqemp_ice(:,:) = -   frcv(jpr_ievp)%z3(:,:,1)        * picefr(:,:)      *   zcptsnw (:,:)   &        ! ice evap
!!       &             +   zsprecip(:,:)                   * zsnw             * zqprec_ice(:,:) * r1_rhos  ! solid precip over ice

      ! --- total non solar flux (including evap/precip) --- !
      zqns_tot(:,:) = zqns_tot(:,:) + zqemp_ice(:,:) + zqemp_oce(:,:)

      ! --- in case both coupled/forced are active, we must mix values --- !
      IF( ln_mixcpl ) THEN
         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) + zqns_tot(:,:)* zmsk(:,:)
         qns_oce(:,:) = qns_oce(:,:) * xcplmask(:,:,0) + zqns_oce(:,:)* zmsk(:,:)
Guillaume Samson's avatar
Guillaume Samson committed
         DO jl=1,jpl
            qns_ice  (:,:,jl) = qns_ice  (:,:,jl) * xcplmask(:,:,0) +  zqns_ice  (:,:,jl)* zmsk(:,:)
            qevap_ice(:,:,jl) = qevap_ice(:,:,jl) * xcplmask(:,:,0) +  zqevap_ice(:,:,jl)* zmsk(:,:)
         ENDDO
         qprec_ice(:,:) = qprec_ice(:,:) * xcplmask(:,:,0) + zqprec_ice(:,:)* zmsk(:,:)
         qemp_oce (:,:) =  qemp_oce(:,:) * xcplmask(:,:,0) +  zqemp_oce(:,:)* zmsk(:,:)
         qemp_ice (:,:) =  qemp_ice(:,:) * xcplmask(:,:,0) +  zqemp_ice(:,:)* zmsk(:,:)
      ELSE
         qns_tot  (:,:  ) = zqns_tot  (:,:  )
         qns_oce  (:,:  ) = zqns_oce  (:,:  )
         qns_ice  (:,:,:) = zqns_ice  (:,:,:)
         qevap_ice(:,:,:) = zqevap_ice(:,:,:)
         qprec_ice(:,:  ) = zqprec_ice(:,:  )
         qemp_oce (:,:  ) = zqemp_oce (:,:  )
         qemp_ice (:,:  ) = zqemp_ice (:,:  )
      ENDIF

!! for CICE ??
!!$      ! --- non solar flux over ocean --- !
!!$      zcptsnw (:,:) = zcptn(:,:)
!!$      zcptrain(:,:) = zcptn(:,:)
!!$
!!$      ! clem: this formulation is certainly wrong... but better than it was...
!!$      zqns_tot(:,:) = zqns_tot(:,:)                             &          ! zqns_tot update over free ocean with:
!!$         &          - (  ziceld(:,:) * zsprecip(:,:) * rLfus )  &          ! remove the latent heat flux of solid precip. melting
!!$         &          - (  zemp_tot(:,:)                          &          ! remove the heat content of mass flux (assumed to be at SST)
!!$         &             - zemp_ice(:,:) ) * zcptn(:,:)
!!$
!!$     IF( ln_mixcpl ) THEN
!!$         qns_tot(:,:) = qns(:,:) * ziceld(:,:) + SUM( qns_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
!!$         qns_tot(:,:) = qns_tot(:,:) * xcplmask(:,:,0) +  zqns_tot(:,:)* zmsk(:,:)
!!$         DO jl=1,jpl
!!$            qns_ice(:,:,jl) = qns_ice(:,:,jl) * xcplmask(:,:,0) +  zqns_ice(:,:,jl)* zmsk(:,:)
!!$         ENDDO
!!$      ELSE
!!$         qns_tot(:,:  ) = zqns_tot(:,:  )
!!$         qns_ice(:,:,:) = zqns_ice(:,:,:)
!!$      ENDIF

      ! outputs
      IF ( srcv(jpr_cal)%laction ) CALL iom_put('hflx_cal_cea' , - frcv(jpr_cal)%z3(:,:,1) * rLfus ) ! latent heat from calving
      IF ( srcv(jpr_icb)%laction ) CALL iom_put('hflx_icb_cea' , - frcv(jpr_icb)%z3(:,:,1) * rLfus ) ! latent heat from icebergs melting
      IF (        iom_use('hflx_rain_cea') )    &                                                    ! heat flux from rain (cell average)
         &   CALL iom_put('hflx_rain_cea' , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )

      IF ( ln_couple_ocean_evap ) THEN
         IF ( iom_use(   'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , frcv(jpr_tevp)%z3(:,:,1) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
      ELSE
         IF ( iom_use(   'hflx_evap_cea') ) CALL iom_put('hflx_evap_cea' , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) )  &
           &                         * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average)
      END IF

Guillaume Samson's avatar
Guillaume Samson committed
      IF (        iom_use('hflx_prec_cea') )    &                                                    ! heat flux from all precip (cell avg)
         &   CALL iom_put('hflx_prec_cea' ,    sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  &
         &                                 + ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )
      IF (        iom_use('hflx_snow_cea') )    &                                                    ! heat flux from snow (cell average)
         &   CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )
      IF (        iom_use('hflx_snow_ao_cea') ) &                                                    ! heat flux from snow (over ocean)
         &   CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) * ( 1._wp - zsnw(:,:) ) )
      IF (        iom_use('hflx_snow_ai_cea') ) &                                                    ! heat flux from snow (over ice)
         &   CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) *  zsnw(:,:) )
      IF(         iom_use('hflx_subl_cea') )    &                                                    ! heat flux from sublimation
         &   CALL iom_put('hflx_subl_cea' ,   SUM( qevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) * tmask(:,:,1) )
      ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp.
      !
      !                                                      ! ========================= !
      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
      !                                                      ! ========================= !
      CASE ('coupled')
         IF( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
         ELSE
            ! Set all category values equal for the moment
            DO jl=1,jpl
               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
            ENDDO
         ENDIF
      CASE( 'none' )
         zdqns_ice(:,:,:) = 0._wp
      END SELECT

      IF( ln_mixcpl ) THEN
         DO jl=1,jpl
            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
         ENDDO
      ELSE
         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
      ENDIF
      !                                                      ! ========================= !
      SELECT CASE( TRIM( sn_rcv_qsr%cldes ) )                !      solar heat fluxes    !   (qsr)
      !                                                      ! ========================= !
      CASE( 'oce only' )
         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
         ! For the Met Office the only sea ice solar flux is the transmitted qsr which is added onto zqsr_ice
         ! further down. Therefore start zqsr_ice off at zero.
         zqsr_ice(:,:,:) = 0._wp
      CASE( 'conservative' )
         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
         IF( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
            zqsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl)
         ELSE
            ! Set all category values equal for the moment
            DO jl = 1, jpl
               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
            END DO
         ENDIF
      CASE( 'oce and ice' )
         zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1)
         IF( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
            DO jl = 1, jpl
               zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)
               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,jl)
            END DO
         ELSE
            zqsr_tot(:,:) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1)
            DO jl = 1, jpl
               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1)
            END DO
         ENDIF
      CASE( 'mixed oce-ice' )
         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
! ** NEED TO SORT OUT HOW THIS SHOULD WORK IN THE MULTI-CATEGORY CASE - CURRENTLY NOT ALLOWED WHEN INTERFACE INITIALISED **
!       Create solar heat flux over ice using incoming solar heat flux and albedos
!       ( see OASIS3 user guide, 5th edition, p39 )
         IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN
            DO jl = 1, jpl
               zqsr_ice(:,:,jl) = frcv(jpr_qsrmix)%z3(:,:,jl) * ( 1.- palbi(:,:,jl) )   &
                  &            / (  1.- ( alb_oce_mix(:,:   ) * ziceld(:,:)       &
                  &                     + palbi      (:,:,jl) * picefr(:,:) ) )
            END DO
         ELSE
            DO jl = 1, jpl
               zqsr_ice(:,:,jl) = frcv(jpr_qsrmix)%z3(:,:, 1) * ( 1.- palbi(:,:,jl) )   &
                  &            / (  1.- ( alb_oce_mix(:,:   ) * ziceld(:,:)       &
                  &                     + palbi      (:,:,jl) * picefr(:,:) ) )
            END DO
         ENDIF
      CASE( 'none'      )       ! Not available as for now: needs additional coding
      !                         ! since fields received, here zqsr_tot,  are not defined with none option
         CALL ctl_stop('STOP', 'sbccpl/sbc_cpl_ice_flx: some fields are not defined. Change sn_rcv_qsr value in namelist namsbc_cpl')
      END SELECT
      IF( ln_dm2dc .AND. ln_cpl ) THEN   ! modify qsr to include the diurnal cycle
         zqsr_tot(:,:  ) = sbc_dcy( zqsr_tot(:,:  ) )
         DO jl = 1, jpl
            zqsr_ice(:,:,jl) = sbc_dcy( zqsr_ice(:,:,jl) )
         END DO
      ENDIF
      !                                                      ! ========================= !
      !                                                      !      Transmitted Qsr      !   [W/m2]
      !                                                      ! ========================= !
      IF( .NOT.ln_cndflx ) THEN                              !==  No conduction flux as surface forcing  ==!
         !
         IF( nn_qtrice == 0 ) THEN
            ! formulation derived from Grenfell and Maykut (1977), where transmission rate
            !    1) depends on cloudiness
            !       ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81)
            !       !      should be real cloud fraction instead (as in the bulk) but needs to be read from atm.
            !    2) is 0 when there is any snow
            !    3) tends to 1 for thin ice
            ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm
            DO jl = 1, jpl
               WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )       ! linear decrease from hi=0 to 10cm
                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) )
               ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )       ! constant (ztri) when hi>10cm
                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ztri(:,:)
               ELSEWHERE                                                           ! zero when hs>0
                  zqtr_ice_top(:,:,jl) = 0._wp
               END WHERE
            ENDDO
         ELSEIF( nn_qtrice == 1 ) THEN
            ! formulation is derived from the thesis of M. Lebrun (2019).
            !    It represents the best fit using several sets of observations
            !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90)
            zqtr_ice_top(:,:,:) = 0.3_wp * zqsr_ice(:,:,:)
         ENDIF
         !
      ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==!
         !
         SELECT CASE( TRIM( sn_rcv_qtr%cldes ) )
         !
         !      ! ===> here we receive the qtr_ice_top array from the coupler
         CASE ('coupled')

            IF (ln_scale_ice_flux) THEN
               WHERE( a_i(:,:,:) > 0.0_wp ) zqtr_ice_top(:,:,:) = MAX(0._wp, frcv(jpr_qtr)%z3(:,:,:)) * a_i_last_couple(:,:,:) / a_i(:,:,:)
               WHERE( a_i(:,:,:) <= 0.0_wp ) zqtr_ice_top(:,:,:) = 0.0_wp
            ELSE
               zqtr_ice_top(:,:,:) = MAX(0._wp, frcv(jpr_qtr)%z3(:,:,:))
            ENDIF

            ! Add retrieved transmitted solar radiation onto the ice and total solar radiation
            zqsr_ice(:,:,:) = zqsr_ice(:,:,:) + zqtr_ice_top(:,:,:)
            zqsr_tot(:,:) = zqsr_tot(:,:) + SUM( zqtr_ice_top(:,:,:) * a_i(:,:,:), dim=3 )

         !      if we are not getting this data from the coupler then assume zero (fully opaque ice)
         CASE ('none')
            zqtr_ice_top(:,:,:) = 0._wp
         CASE default
            CALL ctl_stop( 'sbc_cpl_ice_flx: Invalid value for sn_rcv_qtr%cldes' )
       END SELECT

         !
Guillaume Samson's avatar
Guillaume Samson committed
      ENDIF

      IF( ln_mixcpl ) THEN
         qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) + zqsr_tot(:,:) * zmsk(:,:)
         DO jl = 1, jpl
            qsr_ice    (:,:,jl) = qsr_ice    (:,:,jl) * xcplmask(:,:,0) + zqsr_ice    (:,:,jl) * zmsk(:,:)
            qtr_ice_top(:,:,jl) = qtr_ice_top(:,:,jl) * xcplmask(:,:,0) + zqtr_ice_top(:,:,jl) * zmsk(:,:)
         END DO
      ELSE
         qsr_tot    (:,:  ) = zqsr_tot    (:,:  )
         qsr_ice    (:,:,:) = zqsr_ice    (:,:,:)
         qtr_ice_top(:,:,:) = zqtr_ice_top(:,:,:)
      ENDIF
      
      ! --- solar flux over ocean --- !
      !         note: ziceld cannot be = 0 since we limit the ice concentration to amax
      zqsr_oce = 0._wp
      WHERE( ziceld /= 0._wp )  zqsr_oce(:,:) = ( zqsr_tot(:,:) - SUM( a_i * zqsr_ice, dim=3 ) ) / ziceld(:,:)

      IF( ln_mixcpl ) THEN   ;   qsr_oce(:,:) = qsr_oce(:,:) * xcplmask(:,:,0) +  zqsr_oce(:,:)* zmsk(:,:)
      ELSE                   ;   qsr_oce(:,:) = zqsr_oce(:,:)   ;   ENDIF

      IF( ln_mixcpl ) THEN
         qsr_tot(:,:) = qsr(:,:) * ziceld(:,:) + SUM( qsr_ice(:,:,:) * a_i(:,:,:), dim=3 )   ! total flux from blk
         qsr_tot(:,:) = qsr_tot(:,:) * xcplmask(:,:,0) +  zqsr_tot(:,:)* zmsk(:,:)
         DO jl = 1, jpl
            qsr_ice(:,:,jl) = qsr_ice(:,:,jl) * xcplmask(:,:,0) +  zqsr_ice(:,:,jl)* zmsk(:,:)
         END DO
      ELSE
         qsr_tot(:,:  ) = zqsr_tot(:,:  )
         qsr_ice(:,:,:) = zqsr_ice(:,:,:)
      ENDIF

      !                                                      ! ========================= !
      SELECT CASE( TRIM( sn_rcv_dqnsdt%cldes ) )             !          d(qns)/dt        !
      !                                                      ! ========================= !
      CASE ('coupled')
         IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN
            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl)
         ELSE
            ! Set all category values equal for the moment
            DO jl=1,jpl
               zdqns_ice(:,:,jl) = frcv(jpr_dqnsdt)%z3(:,:,1)
            ENDDO
         ENDIF
      CASE( 'none' ) 
         zdqns_ice(:,:,:) = 0._wp
      CASE default
         CALL ctl_stop( 'sbc_cpl_ice_flx: Invalid value for sn_rcv_dqnsdt%cldes' )
      END SELECT
      
      IF( ln_mixcpl ) THEN
         DO jl=1,jpl
            dqns_ice(:,:,jl) = dqns_ice(:,:,jl) * xcplmask(:,:,0) + zdqns_ice(:,:,jl) * zmsk(:,:)
         ENDDO
      ELSE
         dqns_ice(:,:,:) = zdqns_ice(:,:,:)
      ENDIF


Guillaume Samson's avatar
Guillaume Samson committed
      !                                                      ! ================== !
      !                                                      !   ice skin temp.   !
      !                                                      ! ================== !
      ! needed by Met Office
      IF( srcv(jpr_ts_ice)%laction ) THEN
         WHERE    ( frcv(jpr_ts_ice)%z3(:,:,:) > 0.0  )   ;   ztsu(:,:,:) =   0. + rt0
         ELSEWHERE( frcv(jpr_ts_ice)%z3(:,:,:) < -60. )   ;   ztsu(:,:,:) = -60. + rt0
         ELSEWHERE                                        ;   ztsu(:,:,:) = frcv(jpr_ts_ice)%z3(:,:,:) + rt0
         END WHERE
         !
         IF( ln_mixcpl ) THEN
            DO jl=1,jpl
               pist(:,:,jl) = pist(:,:,jl) * xcplmask(:,:,0) + ztsu(:,:,jl) * zmsk(:,:)
            ENDDO
         ELSE
            pist(:,:,:) = ztsu(:,:,:)
         ENDIF
         !
      ENDIF
      !
#endif
      !
   END SUBROUTINE sbc_cpl_ice_flx


   SUBROUTINE sbc_cpl_snd( kt, Kbb, Kmm )
      !!----------------------------------------------------------------------
      !!             ***  ROUTINE sbc_cpl_snd  ***
      !!
      !! ** Purpose :   provide the ocean-ice informations to the atmosphere
      !!
      !! ** Method  :   send to the atmosphere through a call to cpl_snd
      !!              all the needed fields (as defined in sbc_cpl_init)
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt
      INTEGER, INTENT(in) ::   Kbb, Kmm    ! ocean model time level index
      !
      INTEGER ::   ji, jj, jl   ! dummy loop indices      
      INTEGER ::   ikchoix
Guillaume Samson's avatar
Guillaume Samson committed
      INTEGER ::   isec, info   ! local integer
      REAL(wp) ::   zumax, zvmax
      REAL(wp), DIMENSION(jpi,jpj)     ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4
      !!----------------------------------------------------------------------
      !
      isec = ( kt - nit000 ) * NINT( rn_Dt )        ! date of exchanges
      info = OASIS_idle

      zfr_l(:,:) = 1.- fr_i(:,:)
      !                                                      ! ------------------------- !
      !                                                      !    Surface temperature    !   in Kelvin
      !                                                      ! ------------------------- !
      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN

         IF( nn_components == jp_iam_oce ) THEN
            ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)   ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part
         ELSE
            ! we must send the surface potential temperature
            IF( l_useCT )  THEN    
                ztmp1(:,:) =eos_pt_from_ct( CASTSP(ts(:,:,1,jp_tem,Kmm)), CASTSP(ts(:,:,1,jp_sal,Kmm)) )
            ELSE                   
                ztmp1(:,:) = ts(:,:,1,jp_tem,Kmm)
Guillaume Samson's avatar
Guillaume Samson committed
            ENDIF
            !
            SELECT CASE( sn_snd_temp%cldes)
            CASE( 'oce only'             )   
                ztmp1(:,:) =   ztmp1(:,:) + rt0
Guillaume Samson's avatar
Guillaume Samson committed
            CASE( 'oce and ice'          )   ;   ztmp1(:,:) =   ztmp1(:,:) + rt0
               SELECT CASE( sn_snd_temp%clcat )
               CASE( 'yes' )
                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl)
               CASE( 'no' )
                  WHERE( SUM( a_i, dim=3 ) /= 0. )
                     ztmp3(:,:,1) = SUM( tn_ice * a_i, dim=3 ) / SUM( a_i, dim=3 )
                  ELSEWHERE
                     ztmp3(:,:,1) = rt0
                  END WHERE
               CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
               END SELECT
            CASE( 'weighted oce and ice' )   ;   ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)
               SELECT CASE( sn_snd_temp%clcat )
               CASE( 'yes' )
                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
               CASE( 'no' )
                  ztmp3(:,:,:) = 0.0
                  DO jl=1,jpl
                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
                  ENDDO
               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
               END SELECT
            CASE( 'oce and weighted ice')    ;   ztmp1(:,:) =   ts(:,:,1,jp_tem,Kmm) + rt0
               SELECT CASE( sn_snd_temp%clcat )
               CASE( 'yes' )
                  ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
               CASE( 'no' )
                  ztmp3(:,:,:) = 0.0
                  DO jl=1,jpl
                     ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
                  ENDDO
               CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%clcat' )
               END SELECT
            CASE( 'mixed oce-ice'        )
               ztmp1(:,:) = ( ztmp1(:,:) + rt0 ) * zfr_l(:,:)
               DO jl=1,jpl
                  ztmp1(:,:) = ztmp1(:,:) + tn_ice(:,:,jl) * a_i(:,:,jl)
               ENDDO
            CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_temp%cldes' )
            END SELECT
         ENDIF
         IF( ssnd(jps_toce)%laction )   CALL cpl_snd( jps_toce, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
         IF( ssnd(jps_tice)%laction )   CALL cpl_snd( jps_tice, isec, ztmp3, info )
         IF( ssnd(jps_tmix)%laction )   CALL cpl_snd( jps_tmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
      ENDIF
      !
      !                                                      ! ------------------------- !
      !                                                      ! 1st layer ice/snow temp.  !
      !                                                      ! ------------------------- !
#if defined key_si3
      ! needed by  Met Office
      IF( ssnd(jps_ttilyr)%laction) THEN
         SELECT CASE( sn_snd_ttilyr%cldes)
         CASE ('weighted ice')
            ztmp3(:,:,1:jpl) = t1_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_ttilyr%cldes' )
         END SELECT
         IF( ssnd(jps_ttilyr)%laction )   CALL cpl_snd( jps_ttilyr, isec, ztmp3, info )
      ENDIF
#endif
      !                                                      ! ------------------------- !
      !                                                      !           Albedo          !
      !                                                      ! ------------------------- !
      IF( ssnd(jps_albice)%laction ) THEN                         ! ice
          SELECT CASE( sn_snd_alb%cldes )
          CASE( 'ice' )
             SELECT CASE( sn_snd_alb%clcat )
             CASE( 'yes' )
                ztmp3(:,:,1:jpl) = alb_ice(:,:,1:jpl)
             CASE( 'no' )
                WHERE( SUM( a_i, dim=3 ) /= 0. )
                   ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 )
                ELSEWHERE
                   ztmp1(:,:) = alb_oce_mix(:,:)
                END WHERE
             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' )
             END SELECT
          CASE( 'weighted ice' )   ;
             SELECT CASE( sn_snd_alb%clcat )
             CASE( 'yes' )
                ztmp3(:,:,1:jpl) =  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
             CASE( 'no' )
                WHERE( fr_i (:,:) > 0. )
                   ztmp1(:,:) = SUM (  alb_ice(:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 )
                ELSEWHERE
                   ztmp1(:,:) = 0.
                END WHERE
             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_ice%clcat' )
             END SELECT
          CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%cldes' )
         END SELECT

         SELECT CASE( sn_snd_alb%clcat )
            CASE( 'yes' )
               CALL cpl_snd( jps_albice, isec, ztmp3, info )      !-> MV this has never been checked in coupled mode
            CASE( 'no'  )
               CALL cpl_snd( jps_albice, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
         END SELECT
      ENDIF

      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean
         ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:)
         DO jl = 1, jpl
            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl)
         END DO
         CALL cpl_snd( jps_albmix, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )
      ENDIF
      !                                                      ! ------------------------- !
      !                                                      !  Ice fraction & Thickness !
      !                                                      ! ------------------------- !
      ! Send ice fraction field to atmosphere
      IF( ssnd(jps_fice)%laction ) THEN
         SELECT CASE( sn_snd_thick%clcat )
         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
         END SELECT
         CALL cpl_snd( jps_fice, isec, ztmp3, info )
      ENDIF

#if defined key_si3 || defined key_cice
      ! If this coupling was successful then save ice fraction for use between coupling points.
      ! This is needed for some calculations where the ice fraction at the last coupling point
      ! is needed.
      IF(  info == OASIS_Sent    .OR. info == OASIS_ToRest .OR. &
         & info == OASIS_SentOut .OR. info == OASIS_ToRestOut ) THEN
         IF ( sn_snd_thick%clcat == 'yes' ) THEN
           a_i_last_couple(:,:,1:jpl) = a_i(:,:,1:jpl)
         ENDIF
      ENDIF
#endif

      IF( ssnd(jps_fice1)%laction ) THEN
         SELECT CASE( sn_snd_thick1%clcat )
         CASE( 'yes' )   ;   ztmp3(:,:,1:jpl) =  a_i(:,:,1:jpl)
         CASE( 'no'  )   ;   ztmp3(:,:,1    ) = fr_i(:,:      )
         CASE default    ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick1%clcat' )
         END SELECT
         CALL cpl_snd( jps_fice1, isec, ztmp3, info )
      ENDIF

      ! Send ice fraction field to OCE (sent by SAS in SAS-OCE coupling)
      IF( ssnd(jps_fice2)%laction ) THEN
         ztmp3(:,:,1) = fr_i(:,:)
         IF( ssnd(jps_fice2)%laction )   CALL cpl_snd( jps_fice2, isec, ztmp3, info )
      ENDIF

      ! Send ice and snow thickness field
      IF( ssnd(jps_hice)%laction .OR. ssnd(jps_hsnw)%laction ) THEN
         SELECT CASE( sn_snd_thick%cldes)
         CASE( 'none'                  )       ! nothing to do
         CASE( 'weighted ice and snow' )
            SELECT CASE( sn_snd_thick%clcat )
            CASE( 'yes' )
               ztmp3(:,:,1:jpl) =  h_i(:,:,1:jpl) * a_i(:,:,1:jpl)
               ztmp4(:,:,1:jpl) =  h_s(:,:,1:jpl) * a_i(:,:,1:jpl)
            CASE( 'no' )
               ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
               DO jl=1,jpl
                  ztmp3(:,:,1) = ztmp3(:,:,1) + h_i(:,:,jl) * a_i(:,:,jl)
                  ztmp4(:,:,1) = ztmp4(:,:,1) + h_s(:,:,jl) * a_i(:,:,jl)
               ENDDO
            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
            END SELECT
         CASE( 'ice and snow'         )
            SELECT CASE( sn_snd_thick%clcat )
            CASE( 'yes' )
               ztmp3(:,:,1:jpl) = h_i(:,:,1:jpl)
               ztmp4(:,:,1:jpl) = h_s(:,:,1:jpl)
            CASE( 'no' )
               WHERE( SUM( a_i, dim=3 ) /= 0. )
                  ztmp3(:,:,1) = SUM( h_i * a_i, dim=3 ) / SUM( a_i, dim=3 )
                  ztmp4(:,:,1) = SUM( h_s * a_i, dim=3 ) / SUM( a_i, dim=3 )
               ELSEWHERE
                 ztmp3(:,:,1) = 0.
                 ztmp4(:,:,1) = 0.
               END WHERE
            CASE default                  ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%clcat' )
            END SELECT
         CASE default                     ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_thick%cldes' )
         END SELECT
         IF( ssnd(jps_hice)%laction )   CALL cpl_snd( jps_hice, isec, ztmp3, info )
         IF( ssnd(jps_hsnw)%laction )   CALL cpl_snd( jps_hsnw, isec, ztmp4, info )
      ENDIF

#if defined key_si3
      !                                                      ! ------------------------- !
      !                                                      !      Ice melt ponds       !
      !                                                      ! ------------------------- !
      ! needed by Met Office: 1) fraction of ponded ice 2) local/actual pond depth
      IF( ssnd(jps_a_p)%laction .OR. ssnd(jps_ht_p)%laction ) THEN
         SELECT CASE( sn_snd_mpnd%cldes)
         CASE( 'ice only' )
            SELECT CASE( sn_snd_mpnd%clcat )
            CASE( 'yes' )
               ztmp3(:,:,1:jpl) =  a_ip_eff(:,:,1:jpl)
               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)
            CASE( 'no' )
               ztmp3(:,:,:) = 0.0
               ztmp4(:,:,:) = 0.0
               DO jl=1,jpl
                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)
                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl)
               ENDDO
            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )
            END SELECT
         CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%cldes' )
         END SELECT
         IF( ssnd(jps_a_p)%laction  )   CALL cpl_snd( jps_a_p , isec, ztmp3, info )
         IF( ssnd(jps_ht_p)%laction )   CALL cpl_snd( jps_ht_p, isec, ztmp4, info )
      ENDIF
      !
      !                                                      ! ------------------------- !
      !                                                      !     Ice conductivity      !
      !                                                      ! ------------------------- !
      ! needed by Met Office
      IF( ssnd(jps_kice)%laction ) THEN
         SELECT CASE( sn_snd_cond%cldes)
         CASE( 'weighted ice' )
            SELECT CASE( sn_snd_cond%clcat )
            CASE( 'yes' )
	       ztmp3(:,:,1:jpl) =  cnd_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
            CASE( 'no' )
               ztmp3(:,:,:) = 0.0
               DO jl=1,jpl
                 ztmp3(:,:,1) = ztmp3(:,:,1) + cnd_ice(:,:,jl) * a_i(:,:,jl)
               ENDDO
            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%clcat' )
            END SELECT
         CASE( 'ice only' )
           ztmp3(:,:,1:jpl) = cnd_ice(:,:,1:jpl)
         CASE default      ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_cond%cldes' )
         END SELECT
         IF( ssnd(jps_kice)%laction )   CALL cpl_snd( jps_kice, isec, ztmp3, info )
      ENDIF
#endif

      !                                                      ! ------------------------- !
      !                                                      !  CO2 flux from PISCES     !
      !                                                      ! ------------------------- !
      IF( ssnd(jps_co2)%laction .AND. l_co2cpl )   THEN  
         ztmp1(:,:) = oce_co2(:,:) * 1000.  ! conversion in molC/m2/s 
         CALL cpl_snd( jps_co2, isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ) , info ) 
      ENDIF 
#if defined key_medusa 
      !
      IF (ln_medusa) THEN
      !                                                      ! ---------------------------------------------- !
      !                                                      !  CO2 flux, DMS and chlorophyll from MEDUSA     !
      !                                                      ! ---------------------------------------------- !
         IF ( ssnd(jps_bio_co2)%laction ) THEN
            CALL cpl_snd( jps_bio_co2, isec, RESHAPE( CO2Flux_out_cpl, (/jpi,jpj,1/) ), info )
         ENDIF

         IF ( ssnd(jps_bio_dms)%laction )  THEN
            CALL cpl_snd( jps_bio_dms, isec, RESHAPE( DMS_out_cpl, (/jpi,jpj,1/) ), info )
         ENDIF

         IF ( ssnd(jps_bio_chloro)%laction )  THEN
            CALL cpl_snd( jps_bio_chloro, isec, RESHAPE( chloro_out_cpl, (/jpi,jpj,1/) ), info )
         ENDIF
Guillaume Samson's avatar
Guillaume Samson committed
      ENDIF
      !
Guillaume Samson's avatar
Guillaume Samson committed
      !                                                      ! ------------------------- !
      IF( ssnd(jps_ocx1)%laction ) THEN                      !      Surface current      !
         !                                                   ! ------------------------- !
         !
         !                                                  j+1   j     -----V---F
         ! surface velocity always sent from T point                     !       |
         !  [except for HadGEM3]                                  j      |   T   U