source: NEMO/branches/2019/ENHANCE-02_ISF_nemo/src/OCE/ISF/isfpar.F90 @ 11395

MODULE isfpar
   !!======================================================================
   !!                       ***  MODULE  isfpar  ***
   !! ice shelf module :  update ocean boundary condition under ice
   !!                   shelf
   !!======================================================================
   !! History :  3.2  !  2011-02  (C.Harris  ) Original code isf cav
   !!            X.X  !  2006-02  (C. Wang   ) Original code bg03
   !!            3.4  !  2013-03  (P. Mathiot) Merging + parametrization
   !!            4.0  !  2019-09  (P. Mathiot) Restructuration
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   isf_par       : compute ice shelf melt using a prametrisation of ice shelf cavities
   !!----------------------------------------------------------------------
   USE oce            ! ocean dynamics and tracers
   USE isf
   USE isfparmlt
   USE isftbl
   USE isfdiags
   USE isfutils
   USE dom_oce        ! ocean space and time domain
   USE phycst         ! physical constants
   USE eosbn2         ! equation of state
   !
   USE in_out_manager ! I/O manager
   USE iom            ! I/O library
   USE fldread        ! read input field at current time step
   USE lbclnk         !

   IMPLICIT NONE
   PRIVATE

   PUBLIC   isf_par, isf_par_init

   !!----------------------------------------------------------------------
   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
   !! $Id: sbcisf.F90 10536 2019-01-16 19:21:09Z mathiot $
   !! Software governed by the CeCILL license (see ./LICENSE)
   !!----------------------------------------------------------------------
CONTAINS
 
   SUBROUTINE isf_par( kt, ptsc, pqfwf )
      !!---------------------------------------------------------------------
      !!                     ***  ROUTINE sbc_isf_cav  ***
      !!
      !! ** Purpose :   
      !!
      !! ** Method  : 
      !!
      !! ** Action  :   
      !!                
      !!---------------------------------------------------------------------
      !!-------------------------- OUT --------------------------------------
      REAL(wp), DIMENSION(jpi,jpj)     , INTENT(inout) :: pqfwf
      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(inout) :: ptsc
      !!-------------------------- IN  --------------------------------------
      INTEGER, INTENT(in) ::   kt                                           ! ocean time step
      !!---------------------------------------------------------------------
      !!---------------------------------------------------------------------
      REAL(wp), DIMENSION(jpi,jpj) :: zqoce, zqhc, zqlat, zqh
      !!---------------------------------------------------------------------
      !
      ! compute misfkb_par, rhisf_tbl
      rhisf_tbl_par(:,:) = rhisf0_tbl_par(:,:)
      CALL isf_tbl_lvl( ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
      !
      ! compute heat content, latent heat and melt fluxes (2d)
      CALL isfpar_mlt( kt, zqhc, zqoce, pqfwf  )
      !
      ! compute heat and water flux  (change signe directly in the melt subroutine)
      pqfwf(:,:) = pqfwf(:,:) * mskisf_par(:,:)
      zqoce(:,:) = zqoce(:,:) * mskisf_par(:,:)
      zqhc (:,:) = zqhc(:,:)  * mskisf_par(:,:)
      !
      ! compute heat content flux
      zqlat(:,:) = pqfwf(:,:) * rLfusisf    ! 2d latent heat flux (W/m2)
      !
      ! total heat flux
      zqh(:,:) = ( zqhc (:,:) + zqoce(:,:) )
      !
      ! lbclnk on melt
      CALL lbc_lnk_multi( 'isfmlt', zqh, 'T', 1., pqfwf, 'T', 1.)
      !
      ! output fluxes
      CALL isf_diags_flx( misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, 'par', pqfwf, zqoce, zqlat, zqhc)
      !
      ! set temperature content
      ptsc(:,:,jp_tem) = zqh(:,:) * r1_rau0_rcp
      !
   END SUBROUTINE isf_par

   SUBROUTINE isf_par_init
      !!---------------------------------------------------------------------
      !!                  ***  ROUTINE isf_par_init  ***
      !!
      !! ** Purpose : 
      !!
      !! ** Method  :  
      !!----------------------------------------------------------------------
      INTEGER               :: ierr
      REAL(wp), DIMENSION(jpi,jpj) :: ztblmax, ztblmin
      !!----------------------------------------------------------------------
      !
      ! allocation
      CALL isf_alloc_par()
      !
      ! define isf tbl tickness, top and bottom indice
      CALL read_2dcstdta(TRIM(sn_isfpar_zmax%clname), TRIM(sn_isfpar_zmax%clvar), ztblmax)
      CALL read_2dcstdta(TRIM(sn_isfpar_zmin%clname), TRIM(sn_isfpar_zmin%clvar), ztblmin)
      !
      ! mask ice shelf parametrisation location
      ztblmax(:,:) = ztblmax(:,:) * ssmask(:,:)
      ztblmin(:,:) = ztblmin(:,:) * ssmask(:,:)
      !
      ! if param used under an ice shelf overwrite ztblmax by the ice shelf draft
      WHERE ( risfdep > 0._wp )
         ztblmin(:,:) = risfdep(:,:)
      END WHERE
      !
      ! compute ktop
      CALL isftbl_ktop(ztblmin, misfkt_par)
      !
      ! initial tbl thickness
      rhisf0_tbl_par(:,:) = ztblmax(:,:) - ztblmin(:,:)
      !
      ! compute misfkb_par, rhisf_tbl
      rhisf_tbl(:,:) = rhisf0_tbl_par(:,:)
      CALL isf_tbl_lvl( ht_n, e3t_n, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par )
      !
      ! define iceshelf parametrisation mask
      mskisf_par = 0
      WHERE ( ztblmax > 0._wp )
         mskisf_par(:,:) = 1._wp
      END WHERE
      !
      SELECT CASE ( TRIM(cn_isfpar_mlt) )
         !
      CASE ( 'spe' )

         ALLOCATE( sf_isfpar_fwf(1), STAT=ierr )
         ALLOCATE( sf_isfpar_fwf(1)%fnow(jpi,jpj,1), sf_isfpar_fwf(1)%fdta(jpi,jpj,1,2) )
         CALL fld_fill( sf_isfpar_fwf, (/ sn_isfpar_fwf /), cn_dirisf, 'isf_par_init', 'read fresh water flux isf data', 'namisf' )

         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) '      ==>>>   ice melt read from forcing field (cn_isfmlt_par = spe)'

      CASE ( 'bg03' )
         !
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) '      ==>>>   bg03 parametrisation (cn_isfmlt_par = bg03)'
         !
         ! read effective length
         CALL read_2dcstdta(TRIM(sn_isfpar_Leff%clname), TRIM(sn_isfpar_Leff%clvar), risfLeff)
         risfLeff = risfLeff*1000.0_wp           !: convertion in m
         !
      CASE ( 'oasis' )

         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) '      ==>>>    isf melt provided by OASIS (cn_isfmlt_par = oasis)'

      CASE DEFAULT
         CALL ctl_stop( 'sbc_isf_init: wrong value of nn_isf' )
      END SELECT
      !
   END SUBROUTINE isf_par_init

   END MODULE isfpar