source: trunk/NEMO/LIM_SRC/limdmp.F90 @ 508

MODULE limdmp
   !!======================================================================
   !!                       ***  MODULE limdmp   ***
   !!  Ice model : restoring Ice thickness and Fraction leads
   !!======================================================================
   !! History :   2.0  !  04-04 (S. Theetten) Original code
   !!----------------------------------------------------------------------
#if defined key_ice_lim   &&   defined key_tradmp
   !!----------------------------------------------------------------------
   !!   'key_ice_lim'  AND                                LIM sea-ice model
   !!   'key_tradmp'                                                Damping
   !!----------------------------------------------------------------------
   !!----------------------------------------------------------------------
   !!   lim_dmp      : ice model damping
   !!----------------------------------------------------------------------
   USE in_out_manager  ! I/O manager
   USE phycst          ! physical constants
   USE ice
   USE ice_oce
   USE tradmp
   USE dom_oce
   USE oce
   USE daymod          ! calendar
   USE iom
   
   IMPLICIT NONE
   PRIVATE

   PUBLIC   lim_dmp     ! called by ice_step
   
   INTEGER                        ::   nice1, nice2,  &  ! first and second record used
      &                                inumice_dmp       ! logical unit for ice variables (damping)
   REAL(wp), DIMENSION(jpi,jpj)   ::   hicif_dta  ,   &  ! ice thickness at a given time
      &                                frld_dta          ! fraction lead at a given time
   REAL(wp), DIMENSION(jpi,jpj,2) ::   hicif_data ,   &  ! ice thickness data at two consecutive times
      &                                frld_data         ! fraction lead data at two consecutive times

   !! * Substitution
#  include "vectopt_loop_substitute.h90"
   !!----------------------------------------------------------------------
   !!   LIM 2.0 , UCL-LOCEAN-IPSL  (2006)
   !! $Header$
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE lim_dmp(kt)
      !!-------------------------------------------------------------------
      !!                   ***  ROUTINE lim_dmp ***
      !!
      !! ** purpose : ice model damping : restoring ice thickness and 
      !!              fraction leads
      !!
      !! ** method  : the key_tradmp must be used to compute resto(:,:) coef.
      !!---------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt     ! ocean time-step
      !
      INTEGER             ::   ji, jj         ! dummy loop indices
      !!---------------------------------------------------------------------
      !
      CALL dta_lim( kt )

      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1   ! vector opt.
            hicif(ji,jj) = hicif(ji,jj) - rdt_ice * resto(ji,jj,1) * ( hicif(ji,jj) - hicif_dta(ji,jj) )
            frld(ji,jj)  = frld (ji,jj) - rdt_ice * resto(ji,jj,1) * ( frld(ji,jj)  - frld_dta (ji,jj) )  
         END DO
      END DO
      !
   END SUBROUTINE lim_dmp


   SUBROUTINE dta_lim( kt ) 
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE dta_lim  ***
      !!
      !! ** Purpose :   Reads monthly ice thickness and fraction lead  data
      !!
      !! ** Method  :   Read on unit numicedt the interpolated ice variable
      !!      onto the model grid.
      !!      Data begin at january.
      !!      The value is centered at the middle of month.
      !!      In the opa model, kt=1 agree with january 1.
      !!      At each time step, a linear interpolation is applied between
      !!      two monthly values.
      !!      
      !! ** Action  :   define hicif_dta and frld_dta arrays at time-step kt
      !!----------------------------------------------------------------------
      INTEGER, INTENT(in) ::   kt     ! ocean time-step
      !
      INTEGER  ::   imois, iman, i15          ! temporary integers
      REAL(wp) ::   zxy
      !!----------------------------------------------------------------------

      ! 0. Initialization
      ! -----------------
      iman  = INT( raamo )
!!! better but change the results     i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) )
      i15   = nday / 16
      imois = nmonth + i15 - 1
      IF( imois == 0 ) imois = iman
      
      ! 1. First call kt=nit000: Initialization and Open
      ! -----------------------
      IF( kt == nit000 ) THEN
         nice1 = 0
         IF(lwp) WRITE(numout,*)
         IF(lwp) WRITE(numout,*) 'dtalim : Ice thickness and lead fraction  monthly fields'
         IF(lwp) WRITE(numout,*) '~~~~~~'
         IF(lwp) WRITE(numout,*) '             NetCDF FORMAT'
         IF(lwp) WRITE(numout,*)
         ! open file
         CALL iom_open( 'ice_damping_ATL4.nc', inumice_dmp )
      ENDIF


      ! 2. Read monthly file
      ! --------------------
      IF( ( kt == nit000 ) .OR. imois /= nice1 ) THEN
         !
         ! Calendar computation
         nice1 = imois        ! first file record used 
         nice2 = nice1 + 1    ! last  file record used
         nice1 = MOD( nice1, iman )
         nice2 = MOD( nice2, iman )
         IF( nice1 == 0 )   nice1 = iman
         IF( nice2 == 0 )   nice2 = iman
         IF(lwp) WRITE(numout,*) 'first record file used nice1 ', nice1
         IF(lwp) WRITE(numout,*) 'last  record file used nice2 ', nice2
         
         ! Read monthly ice thickness Levitus 
         CALL iom_get( inumice_dmp, jpdom_data, 'iicethic', hicif_data(:,:,1), nice1 ) 
         CALL iom_get( inumice_dmp, jpdom_data, 'iicethic', hicif_data(:,:,2), nice2 ) 
         
         ! Read monthly ice thickness Levitus 
         CALL iom_get( inumice_dmp, jpdom_data, 'ileadfra', frld_data(:,:,1), nice1 ) 
         CALL iom_get( inumice_dmp, jpdom_data, 'ileadfra', frld_data(:,:,2), nice2 ) 
         
         ! The fraction lead read in the file is in fact the 
         ! ice concentration which is 1 - the fraction lead
         frld_data = 1 - frld_data          
         
         IF(lwp) THEN
            WRITE(numout,*)
            WRITE(numout,*) ' Ice thickness month ', nice1,' and ', nice2
            WRITE(numout,*)
            WRITE(numout,*) ' Ice thickness month = ', nice1
            CALL prihre( hicif_data(1,1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
            WRITE(numout,*)
            WRITE(numout,*) ' Fraction lead months ', nice1,' and ', nice2
            WRITE(numout,*)
            WRITE(numout,*) ' Fraction lead month = ', nice1
            CALL prihre( frld_data(1,1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
         ENDIF
         
         ! 2. At every time step compute ice thickness and fraction lead data
         ! ------------------------------------------------------------------
         zxy = FLOAT( nday + 15 - 30 * i15 ) / 30.
         hicif_dta(:,:) = (1.-zxy) * hicif_data(:,:,1) + zxy * hicif_data(:,:,2)
         frld_dta(:,:) = (1.-zxy) * frld_data(:,:,1) + zxy * frld_data(:,:,2)

      ENDIF
      
      IF( kt == nitend )   CALL iom_close( inumice_dmp )
      !
   END SUBROUTINE dta_lim

#else
   !!----------------------------------------------------------------------
   !!   Default option         Empty Module                  No ice damping
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE lim_dmp( kt )        ! Dummy routine
      WRITE(*,*) 'lim_dmp: You should not see this print! error? ', kt
   END SUBROUTINE lim_dmp
#endif

   !!======================================================================
END MODULE limdmp