source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90 @ 4291

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

   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, ztest_1, ztest_2, ztest_3, ztest_4, ztests  
      REAL(wp) :: zarg, zV, zconv
      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
      REAL(wp) :: epsi06 = 1.0e-6
      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
      REAL(wp), DIMENSION(0:jpl)   ::   zhi_max         !:Boundary of ice thickness categories in thickness space
 
      IF( nn_timing == 1 )  CALL timing_start('limcat_1D')
      !--------------------------------------------------------------------
      ! initialisation of variables
      !--------------------------------------------------------------------
      ijpij = SIZE(zhti,1)
      zht_i(1:ijpij,1:jpl) = 0.d0
      zht_s(1:ijpij,1:jpl) = 0.d0
      za_i (1:ijpij,1:jpl) = 0.d0

      !------------------------------------------------------------------------------------
      ! 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 (ztests=4) 
      !--------------------------------------------------------------------------------------

      !- Thickness categories boundaries
      !    hi_max is calculated in iceini.F90 but since limcat_1D.F90 routine 
      !    is called before (in bdydta.F90), one must recalculate it.   
      !    Note clem: there may be a way of doing things cleaner
      !----------------------------------
      zhi_max(:) = 0._wp
      zc1 =  3._wp / REAL( jpl , wp ) ; zc2 = 10._wp * zc1 ; zc3 =  3._wp
      DO jl = 1, jpl
         zx1 = REAL( jl-1 , wp ) / REAL( jpl , wp )
         zhi_max(jl) = zhi_max(jl-1) + zc1 + zc2 * ( 1._wp + TANH( zc3 * ( zx1 - 1._wp ) ) )
      END DO
 
      ! ----------------------------------------
      ! distribution over the jpl ice categories
      ! ----------------------------------------
      DO ji = 1, ijpij
        ! snow thickness in each category 
         zht_s(ji,1:jpl) = zhts(ji)
         
         ! initialisation of tests
         ztest_1 = 0
         ztest_2 = 0
         ztest_3 = 0
         ztest_4 = 0
         ztests  = 0
         
         i_fill = jpl + 1                                    !====================================
         DO WHILE ( ( ztests /= 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.d0
            za_i (ji,1:jpl) = 0.d0
            
            ! *** 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 (hmax-hmin)/2 
               DO jl = 1, i_fill - 1
                  zht_i(ji,jl) = ( zhi_max(jl) + zhi_max(jl-1) ) / 2.
               END DO
               
               ! find which category (jl0) the input ice thickness falls into
               jl0 = i_fill
               DO jl = 1, i_fill
                  IF ( ( zhti(ji) >= zhi_max(jl-1) ) .AND. ( zhti(ji) < zhi_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) / 2. )
                  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 ) ztest_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 ) ztest_2 = 1
            
            ! Test 3: thickness of the last category is in-bounds ?
            IF ( zht_i(ji,i_fill) >= zhi_max(i_fill-1) ) ztest_3 = 1
            
            ! Test 4: positivity of ice concentrations
            ztest_4 = 1
            DO jl = 1, i_fill
               IF ( za_i(ji,jl) < 0.0d0 ) ztest_4 = 0
            END DO
            
            ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
                                                           !============================
         END DO                                            ! end iteration on categories
                                                           !============================
         ! Check if tests have passed (i.e. volume conservation...)
         !IF ( ztests /= 4 ) THEN
         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
         !   WRITE(numout,*) ' !! ALERT categories distribution !!'
         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
         !   WRITE(numout,*) ' *** ztests is not equal to 4 '
         !   WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4
         !   WRITE(numout,*) 'i_fill=',i_fill
         !   WRITE(numout,*) 'zai(ji)=',zai(ji)
         !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:)
         !ENDIF
            
      END DO ! i loop

      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