source: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90 @ 5051

MODULE limcat_1D
   !!======================================================================
   !!                     ***  MODULE  limcat_1D  ***
   !!  Used for LIM3 to convert cell averages of ice thickness, snow thickness
   !!  and ice cover into a prescribed distribution over the cell. 
   !!  (Example of application: BDY forcings when input are cell averaged)  
   !!======================================================================
   !! History :   -   ! Original code from M. Vancoppenolle (?)
   !!                 ! 2011-12 (C. Rousset) rewritten for clarity
   !!----------------------------------------------------------------------
#if defined key_lim3
   !!----------------------------------------------------------------------
   !!   'key_lim3' :                                    LIM3 sea-ice model
   !!----------------------------------------------------------------------
   !!   lim_cat_1D      :  main subroutine
   !!----------------------------------------------------------------------
   !! Modules used
   USE phycst
   USE oce             ! dynamics and tracers variables
   USE dom_oce
   USE sbc_oce         ! Surface boundary condition: ocean fields
   USE par_ice         ! ice parameters
   USE ice             ! ice variables
   USE eosbn2          ! equation of state
   USE in_out_manager
   USE dom_ice
   USE ice
   USE lbclnk
   USE timing          ! Timing
   USE wrk_nemo

   IMPLICIT NONE
   PRIVATE

   !! Accessibility
   PUBLIC lim_cat_1D      

CONTAINS

   SUBROUTINE lim_cat_1D( zhti, zhts, zai, zht_i, zht_s, za_i )
      !! Local variables
      INTEGER  :: ji, jk, jl             ! dummy loop indices
      INTEGER  :: ijpij, i_fill, jl0  
      REAL(wp) :: zarg, zV, zconv, zdh
      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
      INTEGER , POINTER, DIMENSION(:)         ::   itest
 
      IF( nn_timing == 1 )  CALL timing_start('limcat_1D')

      CALL wrk_alloc( 4, itest )
      !--------------------------------------------------------------------
      ! initialisation of variables
      !--------------------------------------------------------------------
      ijpij = SIZE(zhti,1)
      zht_i(1:ijpij,1:jpl) = 0._wp
      zht_s(1:ijpij,1:jpl) = 0._wp
      za_i (1:ijpij,1:jpl) = 0._wp

      !------------------------------------------------------------------------------------
      ! Distribute ice concentration and thickness into the categories
      !------------------------------------------------------------------------------------
      ! Method: we first try to fill the jpl ice categories bounded by thicknesses
      !         hmax(0:jpl) with a gaussian distribution, and check whether the distribution
      !         fulfills volume and area conservation, positivity and ice categories bounds.
      !         In other words, if ice input is too thin, the last category (jpl)
      !         cannot be filled, so we try to fill jpl-1 categories...
      !         And so forth iteratively until the number of categories filled
      !         fulfills ice volume concervation between input and output (itest=4) 
      !--------------------------------------------------------------------------------------

      ! ----------------------------------------
      ! distribution over the jpl ice categories
      ! ----------------------------------------
      DO ji = 1, ijpij
         
         IF( zhti(ji) > 0._wp ) THEN

         ! initialisation of tests
         itest(:)  = 0
         
         i_fill = jpl + 1                                             !====================================
         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories  
            ! iteration                                               !====================================
            i_fill = i_fill - 1
            
            ! initialisation of ice variables for each try
            zht_i(ji,1:jpl) = 0._wp
            za_i (ji,1:jpl) = 0._wp
            
            ! *** case very thin ice: fill only category 1
            IF ( i_fill == 1 ) THEN
               zht_i(ji,1) = zhti(ji)
               za_i (ji,1) = zai (ji)

            ! *** case ice is thicker: fill categories >1
            ELSE

               ! Fill ice thicknesses except the last one (i_fill) by hmean 
               DO jl = 1, i_fill - 1
                  zht_i(ji,jl) = hi_mean(jl)
               END DO
               
               ! find which category (jl0) the input ice thickness falls into
               jl0 = i_fill
               DO jl = 1, i_fill
                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
                     jl0 = jl
           CYCLE
                  ENDIF
               END DO
               
               ! Concentrations in the (i_fill-1) categories 
               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
               DO jl = 1, i_fill - 1
                  IF ( jl == jl0 ) CYCLE
                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
               END DO
               
               ! Concentration in the last (i_fill) category
               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
               
               ! Ice thickness in the last (i_fill) category
               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 
               
            ENDIF ! case ice is thick or thin
            
            !---------------------
            ! Compatibility tests
            !--------------------- 
            ! Test 1: area conservation
            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
            IF ( zconv < epsi06 ) itest(1) = 1
            
            ! Test 2: volume conservation
            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
            IF ( zconv < epsi06 ) itest(2) = 1
            
            ! Test 3: thickness of the last category is in-bounds ?
            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
            
            ! Test 4: positivity of ice concentrations
            itest(4) = 1
            DO jl = 1, i_fill
               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
            END DO            
                                                           !============================
         END DO                                            ! end iteration on categories
                                                           !============================
         ENDIF ! if zhti > 0
      END DO ! i loop

      ! ------------------------------------------------
      ! Adding Snow in each category where za_i is not 0
      ! ------------------------------------------------ 
      DO jl = 1, jpl
         DO ji = 1, ijpij
            IF( za_i(ji,jl) > 0._wp ) THEN
               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
               ! In case snow load is in excess that would lead to transformation from snow to ice
               ! Then, transfer the snow excess into the ice (different from limthd_dh)
               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
               ! recompute ht_i, ht_s avoiding out of bounds values
               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn )
            ENDIF
         ENDDO
      ENDDO

      CALL wrk_dealloc( 4, itest )
      !
      IF( nn_timing == 1 )  CALL timing_stop('limcat_1D')
     
    END SUBROUTINE lim_cat_1D
    
#else
   !!----------------------------------------------------------------------
   !!   Default option :         Empty module          NO LIM sea-ice model
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE lim_cat_1D          ! Empty routine
   END SUBROUTINE lim_cat_1D
#endif

   !!======================================================================
END MODULE limcat_1D