source: trunk/NEMO/LIM_SRC_2/limdmp_2.F90 @ 1850

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

   PUBLIC   lim_dmp_2     ! called by ice_step_2
   
   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)
   !! $Id$
   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
   !!----------------------------------------------------------------------

CONTAINS

   SUBROUTINE lim_dmp_2(kt)
      !!-------------------------------------------------------------------
      !!                   ***  ROUTINE lim_dmp_2 ***
      !!
      !! ** 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_2( 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_2


   SUBROUTINE dta_lim_2( kt ) 
      !!----------------------------------------------------------------------
      !!                   ***  ROUTINE dta_lim_2  ***
      !!
      !! ** 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

         CALL FLUSH(numout)

      ENDIF
        
      ! 3. 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)

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

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

   !!======================================================================
END MODULE limdmp_2