source: trunk/NEMO/OPA_SRC/SBC/sbcice_if.F90 @ 912

MODULE sbcice_if
   !!======================================================================
   !!                       ***  MODULE  sbcice  ***
   !! Surface module :  update surface ocean boundary condition over ice
   !!                   covered area using ice-if model
   !!======================================================================
   !! History :  9.0   !  06-06  (G. Madec)  Original code
   !!----------------------------------------------------------------------

   !!----------------------------------------------------------------------
   !!   sbc_ice_if     : update sbc in ice-covered area
   !!----------------------------------------------------------------------
   USE oce             ! ocean dynamics and tracers
   USE dom_oce         ! ocean space and time domain
   USE phycst          ! physical constants
   USE ocfzpt          ! ocean freezing point
   USE sbc_oce         ! Surface boundary condition: ocean fields
   USE fldread         ! read input field
   USE iom             ! I/O manager library
   USE in_out_manager  ! I/O manager

   IMPLICIT NONE
   PRIVATE

   PUBLIC   sbc_ice_if      ! routine called in sbcmod

   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_ice   ! structure of input ice-cover (file informations, fields read)
   
   !! * Substitutions
#  include "domzgr_substitute.h90"
   !!----------------------------------------------------------------------
   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
   !! $ Id: $
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE sbc_ice_if( kt )
      !!---------------------------------------------------------------------
      !!                     ***  ROUTINE sbc_ice_if  ***
      !!
      !! ** Purpose :   handle surface boundary condition over ice cover area
      !!      when sea-ice model are not used
      !!
      !! ** Method  : - read sea-ice cover climatology
      !!              - blah blah blah, ...
      !!
      !! ** Action  :   qns, qsr:  update heat flux below sea-ice
      !!                emp, emps: update freshwater flux below sea-ice
      !!---------------------------------------------------------------------
      INTEGER, INTENT(in)          ::   kt         ! ocean time step
      !
      CHARACTER(len=100) ::   cn_dir              ! Root directory for location of ice-if files
      TYPE(FLD_N)        ::   sn_ice              ! informations about the fields to be read
      NAMELIST/namsbc_iif/ cn_dir, sn_ice
      !
      INTEGER  ::   ji, jj     ! dummy loop indices
      INTEGER  ::   ierror     ! return error code
      REAL(wp) ::   ztrp, zsice, zt_fzp, zicover_obs, zicover_opa
      REAL(wp) ::   zqri, zqrj, zqrp, zqi
      !!---------------------------------------------------------------------
      !                                         ! ====================== !
      IF( kt == nit000 ) THEN                   !  First call kt=nit000  !
         !                                      ! ====================== !
         ! set file information
         cn_dir = './'        ! directory in which the model is executed
         ! ... default values (NB: frequency positive => hours, negative => months)
         !            !   file    ! frequency !  variable  ! time intep !  clim  ! starting !
         !            !   name    !  (hours)  !   name     !   (T/F)    !  (0/1) !  record  !
         sn_ice = FLD_N('ice_cover',   -12.   ,  'ice_cov'  ,  .TRUE.    ,    1   ,     0    )

         REWIND ( numnam )               ! ... read in namlist namiif
         READ   ( numnam, namsbc_iif )

         ALLOCATE( sf_ice(1), STAT=ierror )
         IF( ierror > 0 ) THEN
            CALL ctl_stop( 'sbc_ice_if: unable to allocate sf_ice structure' )   ;   RETURN
         ENDIF

         ! store namelist information in sf_ice structure
         WRITE(sf_ice(1)%clrootname,'(a,a)' )   TRIM( cn_dir ), TRIM( sn_ice%clname )
         sf_ice(1)%freqh   = sn_ice%freqh
         sf_ice(1)%clvar   = sn_ice%clvar
         sf_ice(1)%ln_tint = sn_ice%ln_tint
         sf_ice(1)%nclim   = sn_ice%nclim
         sf_ice(1)%nstrec  = sn_ice%nstrec

         IF(lwp) THEN      ! control print
            WRITE(numout,*)
            WRITE(numout,*) 'sbc_ice_if : ice-if sea-ice model'
            WRITE(numout,*) '~~~~~~~~~~ '
            WRITE(numout,*) '   ice-cover data in the following file: '
            WRITE(numout,*) '          list of files and frequency (>0: in hours ; <0 in months)'
            WRITE(numout,*) '               root filename: '  , trim( sf_ice(1)%clrootname ),   &
               &                          ' variable name: '  , trim( sf_ice(1)%clvar      )
            WRITE(numout,*) '               frequency: '      ,       sf_ice(1)%freqh       ,   &
               &                          ' time interp: '    ,       sf_ice(1)%ln_tint     ,   &
               &                          ' climatology: '    ,       sf_ice(1)%nclim       ,   &
               &                          ' starting record: ',       sf_ice(1)%nstrec
         ENDIF
         !
      ENDIF

      CALL fld_read( kt, nn_fsbc, sf_ice )           ! Read input fields and provides the
      !                                              ! input fields at the current time-step
      
      IF( MOD( kt-1, nn_fsbc) == 0 ) THEN
         !
         ztrp = -40.             ! restoring terme for temperature (w/m2/k)
         zsice = - 0.04 / 0.8    ! ratio of isohaline compressibility over isotherme compressibility
                                 ! ( d rho / dt ) / ( d rho / ds )      ( s = 34, t = -1.8 )
         ! Flux computation
!CDIR COLLAPSE
         DO jj = 1, jpj
            DO ji = 1, jpi
               ! ... sea surface freezing point temperature [Celcius]
               zt_fzp = (  ( - 0.0575 + 1.710523e-3 * SQRT( sss_m(ji,jj) )   &
            &                         - 2.154996e-4 *       sss_m(ji,jj)   ) * sss_m(ji,jj)  ) * tmask(ji,jj,1)
            
               ! ... indicators : ice cover (obs, ocean model) & hemisphere (=1 north, =-1 south)
               zicover_obs = sf_ice(1)%fnow(ji,jj)                                                ! observed 
               zicover_opa = MAX( 0., SIGN( 1., zt_fzp - sst_m(ji,jj) )  ) * tmask(ji,jj,1)   ! model   

               ! ... avoid over-freezing point temperature
               tn(ji,jj,1) = MAX( tn(ji,jj,1), zt_fzp )

               ! ... solar heat flux : zero below observed ice cover
               qsr(ji,jj) = ( 1. - zicover_obs ) * qsr(ji,jj)

               ! ... non solar heat flux : add a damping term 
               !      - gamma*(t-(tgel-1.))  if observed ice and no opa ice   (zicover_obs=1 zicover_opa=0)
               !      - gamma*min(0,t-tgel)  if observed ice and opa ice      (zicover_obs=1 zicover_opa=1)

               zqri = ztrp * ( tb(ji,jj,1) - ( zt_fzp - 1.) )
               zqrj = ztrp * MIN( 0., tb(ji,jj,1) - zt_fzp )

               zqrp =  ( zicover_obs * ( (1. - zicover_opa ) * zqri    &
                 &                      +      zicover_opa   * zqrj ) ) * tmask(ji,jj,1)

               ! c) net downward heat flux q() = q0 + qrp()
               ! for q0
               ! # qns unchanged              if no climatological ice              (zicover_obs=0)
               ! # qns = zqrp                 if climatological ice and no opa ice  (zicover_obs=1, zicover_opa=0)
               ! # qns = zqrp -2(-4) watt/m2  if climatological ice and opa ice     (zicover_obs=1, zicover_opa=1)
               ! (-2=arctic, -4=antarctic)   
               zqi = -3. + SIGN( 1.e0, ff(ji,jj) )
               qns(ji,jj) = ( ( 1.- zicover_obs ) * qns(ji,jj)   &
                  &          +      zicover_obs   * zicover_opa * zqi ) * tmask(ji,jj,1)   &
                  &       + zqrp
            END DO
         END DO
         !
      ENDIF
      !
   END SUBROUTINE sbc_ice_if

   !!======================================================================
END MODULE sbcice_if