source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limcat_2D.F90 @ 4317

MODULE limcat_2D
   !!======================================================================
   !!                     ***  MODULE  limcat_2D  ***
   !!  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      :  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

   IMPLICIT NONE
   PRIVATE

   !! Accessibility
   PUBLIC lim_cat_2D      

CONTAINS

   SUBROUTINE lim_cat_2D(zhti,zhts,zai,zht_i,zht_s,za_i)
      !! Local variables
      INTEGER  :: ji, jj, jk, jl             ! dummy loop indices
      INTEGER  :: ijpi, ijpj, 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
  
      !--------------------------------------------------------------------
      ! initialisation of variables
      !--------------------------------------------------------------------
      ijpi = SIZE(zhti,1)
      ijpj = SIZE(zhti,2)
      zht_i(1:ijpi,1:ijpj,1:jpl) = 0.d0
      zht_s(1:ijpi,1:ijpj,1:jpl) = 0.d0
      za_i (1:ijpi,1:ijpj,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) 
      !--------------------------------------------------------------------------------------
 
      ! ----------------------------------------
      ! distribution over the jpl ice categories
      ! ----------------------------------------
      DO jj = 1, ijpj
         DO ji = 1, ijpi

            ! snow thickness in each category 
            zht_s(ji,jj,1:jpl) = zhts(ji,jj)

            ! 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,jj,1:jpl) = 0.d0
               za_i (ji,jj,1:jpl) = 0.d0

               ! *** case very thin ice: fill only category 1
               IF ( i_fill == 1 ) THEN
                  zht_i(ji,jj,1) = zhti(ji,jj)
                  za_i (ji,jj,1) = zai (ji,jj)
            
               ! *** 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,jj,jl) = ( hi_max(jl) + hi_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,jj) >= hi_max(jl-1) ) .AND. ( zhti(ji,jj) < hi_max(jl) ) ) THEN
                        jl0 = jl
                        CYCLE
                     ENDIF
                  END DO

                  ! Concentrations in the (i_fill-1) categories 
                  za_i(ji,jj,jl0) = zai(ji,jj) / SQRT(REAL(jpl))
                  DO jl = 1, i_fill - 1
                     IF ( jl == jl0 ) CYCLE
                     zarg           = ( zht_i(ji,jj,jl) - zhti(ji,jj) ) / ( zhti(ji,jj) / 2. )
                     za_i(ji,jj,jl) =   za_i (ji,jj,jl0) * EXP(-zarg**2)
                  END DO 

                  ! Concentration in the last (i_fill) category
                  za_i(ji,jj,i_fill) = zai(ji,jj) - SUM( za_i(ji,jj,1:i_fill-1) )

                  ! Ice thickness in the last (i_fill) category
                  zV = SUM( za_i(ji,jj,1:i_fill-1) * zht_i(ji,jj,1:i_fill-1) )
                  zht_i(ji,jj,i_fill) = ( zhti(ji,jj)*zai(ji,jj) -  zV ) / za_i(ji,jj,i_fill) 
 
               ENDIF ! case ice is thick or thin

               !---------------------
               ! Compatibility tests
               !--------------------- 
               ! Test 1: area conservation
               zconv = ABS( zai(ji,jj) - SUM( za_i(ji,jj,1:jpl) ) )
               IF ( zconv < epsi06 ) ztest_1 = 1
 
               ! Test 2: volume conservation
               zconv = ABS( zhti(ji,jj)*zai(ji,jj) - SUM( za_i(ji,jj,1:jpl)*zht_i(ji,jj,1:jpl) ) )
               IF ( zconv < epsi06 ) ztest_2 = 1
 
               ! Test 3: thickness of the last category is in-bounds ?
               IF ( zht_i(ji,jj,i_fill) >= hi_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,jj,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 .NE. 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
            !ENDIF

         END DO ! i loop
      END DO ! j loop

   END SUBROUTINE lim_cat_2D

#else
   !!----------------------------------------------------------------------
   !!   Default option :         Empty module          NO LIM sea-ice model
   !!----------------------------------------------------------------------
CONTAINS
   SUBROUTINE lim_cat_2D          ! Empty routine
   END SUBROUTINE lim_cat_2D
#endif

   !!======================================================================
END MODULE limcat_2D