Skip to content
Snippets Groups Projects
icevar.F90 64.2 KiB
Newer Older
Guillaume Samson's avatar
Guillaume Samson committed
         END DO
      END DO
      !
      DEALLOCATE( z1_hti )
      !
      ! == temperature and salinity == !
      DO jl = 1, jpl
         pt_i (:,jl) = ptmi (:)
         pt_s (:,jl) = ptms (:)
         pt_su(:,jl) = ptmsu(:)
         ps_i (:,jl) = psmi (:)
      END DO
      !
      ! == ponds == !
      ALLOCATE( zfra(idim) )
      ! keep the same pond fraction atip/ati for each category
      WHERE( pati(:) /= 0._wp )   ;   zfra(:) = patip(:) / pati(:)
      ELSEWHERE                   ;   zfra(:) = 0._wp
      END WHERE
      DO jl = 1, jpl
         pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
      END DO
      ! keep the same v_ip/v_i ratio for each category
      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) )
      ELSEWHERE                                 ;   zfra(:) = 0._wp
      END WHERE
      DO jl = 1, jpl
         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
         ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp
         END WHERE
      END DO
      ! keep the same v_il/v_i ratio for each category
      WHERE( ( phti(:) * pati(:) ) /= 0._wp )   ;   zfra(:) = ( phtil(:) * patip(:) ) / ( phti(:) * pati(:) )
      ELSEWHERE                                 ;   zfra(:) = 0._wp
      END WHERE
      DO jl = 1, jpl
         WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
         ELSEWHERE                       ;   ph_il(:,jl) = 0._wp
         END WHERE
      END DO
      DEALLOCATE( zfra )
      !
   END SUBROUTINE ice_var_itd_1cMc

   SUBROUTINE ice_var_itd_NcMc( phti, phts, pati ,                             ph_i, ph_s, pa_i, &
      &                         ptmi, ptms, ptmsu, psmi, patip, phtip, phtil,  pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il )
      !!-------------------------------------------------------------------
      !!
      !! ** Purpose :  converting N-cat ice to jpl ice categories
      !!
      !!                  ice thickness distribution follows a gaussian law
      !!               around the concentration of the most likely ice thickness
      !!                           (similar as iceistate.F90)
      !!
      !! ** Method:   Iterative procedure
      !!
      !!               1) Fill ice cat that correspond to input thicknesses
      !!                  Find the lowest(jlmin) and highest(jlmax) cat that are filled
      !!
      !!               2) Expand the filling to the cat jlmin-1 and jlmax+1
      !!                   by removing 25% ice area from jlmin and jlmax (resp.)
      !!
      !!               3) Expand the filling to the empty cat between jlmin and jlmax
      !!                   by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)
      !!                      b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)
      !!
      !! ** Arguments : phti: N-cat ice thickness
      !!                phts: N-cat snow depth
      !!                pati: N-cat ice concentration
      !!
      !! ** Output    : jpl-cat
      !!
      !!  (Example of application: BDY forcings when inputs have N-cat /= jpl)
      !!-------------------------------------------------------------------
      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables
      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   ph_i, ph_s, pa_i    ! output ice/snow variables
      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   ptmi, ptms, ptmsu, psmi, patip, phtip, phtil    ! input  ice/snow temp & sal & ponds
      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip, ph_il    ! output ice/snow temp & sal & ponds
      !
      INTEGER , ALLOCATABLE, DIMENSION(:,:) ::   jlfil, jlfil2
      INTEGER , ALLOCATABLE, DIMENSION(:)   ::   jlmax, jlmin
      REAL(wp), ALLOCATABLE, DIMENSION(:)   ::   z1_ai, z1_vi, z1_vs, ztmp, zfra
      !
      REAL(wp), PARAMETER ::   ztrans = 0.25_wp
      INTEGER  ::   ji, jl, jl1, jl2
      INTEGER  ::   idim, icat
      !!-------------------------------------------------------------------
      !
      idim = SIZE( phti, 1 )
      icat = SIZE( phti, 2 )
      !
      ! == thickness and concentration == !
      !                                 ! ---------------------- !
      IF( icat == jpl ) THEN            ! input cat = output cat !
         !                              ! ---------------------- !
         ph_i(:,:) = phti(:,:)
         ph_s(:,:) = phts(:,:)
         pa_i(:,:) = pati(:,:)
         !
         ! == temperature and salinity and ponds == !
         pt_i (:,:) = ptmi (:,:)
         pt_s (:,:) = ptms (:,:)
         pt_su(:,:) = ptmsu(:,:)
         ps_i (:,:) = psmi (:,:)
         pa_ip(:,:) = patip(:,:)
         ph_ip(:,:) = phtip(:,:)
         ph_il(:,:) = phtil(:,:)
         !                              ! ---------------------- !
      ELSEIF( icat == 1 ) THEN          ! input cat = 1          !
         !                              ! ---------------------- !
         CALL  ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), &
            &                    ph_i(:,:), ph_s(:,:), pa_i (:,:), &
            &                    ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), phtil(:,1), &
            &                    pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:), ph_il(:,:)  )
         !                              ! ---------------------- !
      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         !
         !                              ! ---------------------- !
         CALL  ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), &
            &                    ph_i(:,1), ph_s(:,1), pa_i (:,1), &
            &                    ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), phtil(:,:), &
            &                    pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1), ph_il(:,1)  )
         !                              ! ----------------------- !
      ELSE                              ! input cat /= output cat !
         !                              ! ----------------------- !

         ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
         ALLOCATE( jlmin(idim), jlmax(idim) )

         ! --- initialize output fields to 0 --- !
         ph_i(1:idim,1:jpl) = 0._wp
         ph_s(1:idim,1:jpl) = 0._wp
         pa_i(1:idim,1:jpl) = 0._wp
         !
         ! --- fill the categories --- !
         !     find where cat-input = cat-output and fill cat-output fields
         jlmax(:) = 0
         jlmin(:) = 999
         jlfil(:,:) = 0
         DO jl1 = 1, jpl
            DO jl2 = 1, icat
               DO ji = 1, idim
                  IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN
                     ! fill the right category
                     ph_i(ji,jl1) = phti(ji,jl2)
                     ph_s(ji,jl1) = phts(ji,jl2)
                     pa_i(ji,jl1) = pati(ji,jl2)
                     ! record categories that are filled
                     jlmax(ji) = MAX( jlmax(ji), jl1 )
                     jlmin(ji) = MIN( jlmin(ji), jl1 )
                     jlfil(ji,jl1) = jl1
                  ENDIF
               END DO
            END DO
         END DO
         !
         ! --- fill the gaps between categories --- !
         !     transfer from categories filled at the previous step to the empty ones in between
         DO ji = 1, idim
            jl1 = jlmin(ji)
            jl2 = jlmax(ji)
            IF( jl1 > 1 ) THEN
               ! fill the lower cat (jl1-1)
               pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1)
               ph_i(ji,jl1-1) = hi_mean(jl1-1)
               ! remove from cat jl1
               pa_i(ji,jl1  ) = ( 1._wp - ztrans ) * pa_i(ji,jl1)
            ENDIF
            IF( jl2 < jpl ) THEN
               ! fill the upper cat (jl2+1)
               pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2)
               ph_i(ji,jl2+1) = hi_mean(jl2+1)
               ! remove from cat jl2
               pa_i(ji,jl2  ) = ( 1._wp - ztrans ) * pa_i(ji,jl2)
            ENDIF
         END DO
         !
         jlfil2(:,:) = jlfil(:,:)
         ! fill categories from low to high
         DO jl = 2, jpl-1
            DO ji = 1, idim
               IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
                  ! fill high
                  pa_i(ji,jl) = ztrans * pa_i(ji,jl-1)
                  ph_i(ji,jl) = hi_mean(jl)
                  jlfil(ji,jl) = jl
                  ! remove low
                  pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1)
               ENDIF
            END DO
         END DO
         !
         ! fill categories from high to low
         DO jl = jpl-1, 2, -1
            DO ji = 1, idim
               IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
                  ! fill low
                  pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1)
                  ph_i(ji,jl) = hi_mean(jl)
                  jlfil2(ji,jl) = jl
                  ! remove high
                  pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1)
               ENDIF
            END DO
         END DO
         !
         DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
         DEALLOCATE( jlmin, jlmax )
         !
         ! == temperature and salinity == !
         !
         ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) )
         !
         WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp )               ;   z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 )
         ELSEWHERE                                               ;   z1_ai(:) = 0._wp
         END WHERE
         WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp )   ;   z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 )
         ELSEWHERE                                               ;   z1_vi(:) = 0._wp
         END WHERE
         WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp )   ;   z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 )
         ELSEWHERE                                               ;   z1_vs(:) = 0._wp
         END WHERE
         !
         ! fill all the categories with the same value
         ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
         DO jl = 1, jpl
            pt_i (:,jl) = ztmp(:)
         END DO
         ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:)
         DO jl = 1, jpl
            pt_s (:,jl) = ztmp(:)
         END DO
         ztmp(:) = SUM( ptmsu(:,:) * pati(:,:)            , dim=2 ) * z1_ai(:)
         DO jl = 1, jpl
            pt_su(:,jl) = ztmp(:)
         END DO
         ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:)
         DO jl = 1, jpl
            ps_i (:,jl) = ztmp(:)
         END DO
         !
         DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp )
         !
         ! == ponds == !
         ALLOCATE( zfra(idim) )
         ! keep the same pond fraction atip/ati for each category
         WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp )   ;   zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 )
         ELSEWHERE                                   ;   zfra(:) = 0._wp
         END WHERE
         DO jl = 1, jpl
            pa_ip(:,jl) = zfra(:) * pa_i(:,jl)
         END DO
         ! keep the same v_ip/v_i ratio for each category
         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
            zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
         ELSEWHERE
            zfra(:) = 0._wp
         END WHERE
         DO jl = 1, jpl
            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
            ELSEWHERE                       ;   ph_ip(:,jl) = 0._wp
            END WHERE
         END DO
         ! keep the same v_il/v_i ratio for each category
         WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp )
            zfra(:) = SUM( phtil(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 )
         ELSEWHERE
            zfra(:) = 0._wp
         END WHERE
         DO jl = 1, jpl
            WHERE( pa_ip(:,jl) /= 0._wp )   ;   ph_il(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl)
            ELSEWHERE                       ;   ph_il(:,jl) = 0._wp
            END WHERE
         END DO
         DEALLOCATE( zfra )
         !
      ENDIF
      !
   END SUBROUTINE ice_var_itd_NcMc

   !!-------------------------------------------------------------------
   !! INTERFACE ice_var_snwfra
   !!
   !! ** Purpose :  fraction of ice covered by snow
   !!
   !! ** Method  :  In absence of proper snow model on top of sea ice,
   !!               we argue that snow does not cover the whole ice because
   !!               of wind blowing...
   !!
   !! ** Arguments : ph_s: snow thickness
   !!
   !! ** Output    : pa_s_fra: fraction of ice covered by snow
   !!
   !!-------------------------------------------------------------------
   SUBROUTINE ice_var_snwfra_3d( ph_s, pa_s_fra )
      REAL(wp), DIMENSION(A2D(0),jpl), INTENT(in   ) ::   ph_s        ! snow thickness
      REAL(wp), DIMENSION(A2D(0),jpl), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
Guillaume Samson's avatar
Guillaume Samson committed
      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover
         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
         ELSEWHERE             ; pa_s_fra = 0._wp
         END WHERE
      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style)
         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style)
         pa_s_fra = ph_s / ( ph_s + 0.02_wp )
      ENDIF
   END SUBROUTINE ice_var_snwfra_3d

   SUBROUTINE ice_var_snwfra_2d( ph_s, pa_s_fra )
      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   ph_s        ! snow thickness
      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover
         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
         ELSEWHERE             ; pa_s_fra = 0._wp
         END WHERE
      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style)
         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style)
         pa_s_fra = ph_s / ( ph_s + 0.02_wp )
      ENDIF
   END SUBROUTINE ice_var_snwfra_2d

   SUBROUTINE ice_var_snwfra_1d( ph_s, pa_s_fra )
      REAL(wp), DIMENSION(:), INTENT(in   ) ::   ph_s        ! snow thickness
      REAL(wp), DIMENSION(:), INTENT(  out) ::   pa_s_fra    ! ice fraction covered by snow
      IF    ( nn_snwfra == 0 ) THEN   ! basic 0 or 1 snow cover
         WHERE( ph_s > 0._wp ) ; pa_s_fra = 1._wp
         ELSEWHERE             ; pa_s_fra = 0._wp
         END WHERE
      ELSEIF( nn_snwfra == 1 ) THEN   ! snow cover depends on hsnow (met-office style)
         pa_s_fra = 1._wp - EXP( -0.2_wp * rhos * ph_s )
      ELSEIF( nn_snwfra == 2 ) THEN   ! snow cover depends on hsnow (cice style)
         pa_s_fra = ph_s / ( ph_s + 0.02_wp )
      ENDIF
   END SUBROUTINE ice_var_snwfra_1d

   !!--------------------------------------------------------------------------
   !! INTERFACE ice_var_snwblow
   !!
   !! ** Purpose :   Compute distribution of precip over the ice
   !!
   !!                Snow accumulation in one thermodynamic time step
   !!                snowfall is partitionned between leads and ice.
   !!                If snow fall was uniform, a fraction (1-at_i) would fall into leads
   !!                but because of the winds, more snow falls on leads than on sea ice
   !!                and a greater fraction (1-at_i)^beta of the total mass of snow
   !!                (beta < 1) falls in leads.
   !!                In reality, beta depends on wind speed,
   !!                and should decrease with increasing wind speed but here, it is
   !!                considered as a constant. an average value is 0.66
   !!--------------------------------------------------------------------------
!!gm  I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
   SUBROUTINE ice_var_snwblow_2d( pin, pout )
      REAL(wp), DIMENSION(A2D(0)), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
      REAL(wp), DIMENSION(A2D(0)), INTENT(inout) :: pout
Guillaume Samson's avatar
Guillaume Samson committed
      pout = ( 1._wp - ( pin )**rn_snwblow )
   END SUBROUTINE ice_var_snwblow_2d

   SUBROUTINE ice_var_snwblow_1d( pin, pout )
      REAL(wp), DIMENSION(:), INTENT(in   ) :: pin
      REAL(wp), DIMENSION(:), INTENT(inout) :: pout
      pout = ( 1._wp - ( pin )**rn_snwblow )
   END SUBROUTINE ice_var_snwblow_1d

#else
   !!----------------------------------------------------------------------
   !!   Default option         Dummy module           NO SI3 sea-ice model
   !!----------------------------------------------------------------------
#endif

   !!======================================================================
END MODULE icevar