Changeset 11560


Ignore:
Timestamp:
2019-09-18T13:31:36+02:00 (12 months ago)
Author:
clem
Message:

correct the conversion from 1cat to Ncat for ice data at the boundaries when using bdy

Location:
NEMO/trunk/src/ICE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • NEMO/trunk/src/ICE/icedyn_rdgrft.F90

    r11536 r11560  
    278278      ! controls 
    279279      IF( ln_ctl       )   CALL ice_prt3D   ('icedyn_rdgrft')                                                             ! prints 
     280      IF( ln_icectl    )   CALL ice_prt     (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ')                             ! prints 
    280281      IF( ln_icediachk )   CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 
    281282      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icedyn_rdgrft',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) ! conservation 
  • NEMO/trunk/src/ICE/icevar.F90

    r11536 r11560  
    868868      !! ** Purpose :  converting 1-cat ice to jpl ice categories 
    869869      !! 
    870       !!                  ice thickness distribution follows a gaussian law 
    871       !!               around the concentration of the most likely ice thickness 
    872       !!                           (similar as iceistate.F90) 
    873       !! 
    874       !! ** Method:   Iterative procedure 
    875       !!                 
    876       !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 
    877       !! 
    878       !!               2) Check whether the distribution conserves area and volume, positivity and 
    879       !!                  category boundaries 
     870      !! 
     871      !! ** Method:   ice thickness distribution follows a gamma function from Abraham et al. (2015) 
     872      !!              it has the property of conserving total concentration and volume 
    880873      !!               
    881       !!               3) If not (input ice is too thin), the last category is empty and 
    882       !!                  the number of categories is reduced (jpl-1) 
    883       !! 
    884       !!               4) Iterate until ok (SUM(itest(:) = 4) 
    885874      !! 
    886875      !! ** Arguments : phti: 1-cat ice thickness 
     
    890879      !! ** Output    : jpl-cat  
    891880      !! 
    892       !!  (Example of application: BDY forcings when input are cell averaged)   
     881      !!  Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015. 
     882      !!               Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice. 
     883      !!               Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614  
    893884      !!------------------------------------------------------------------- 
    894885      REAL(wp), DIMENSION(:),   INTENT(in)    ::   phti, phts, pati    ! input  ice/snow variables 
     
    897888      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip    ! output ice/snow temp & sal & ponds 
    898889      ! 
    899       INTEGER , DIMENSION(4) ::   itest 
    900       REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra 
     890      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zfra, z1_hti 
    901891      INTEGER  ::   ji, jk, jl 
    902       INTEGER  ::   idim, i_fill, jl0   
    903       REAL(wp) ::   zarg, zV, zconv, zdh, zdv 
    904       !!------------------------------------------------------------------- 
    905       ! 
    906       ! == thickness and concentration == ! 
    907       ! distribution over the jpl ice categories 
    908       !    a gaussian distribution for ice concentration is used 
    909       !    then we check whether the distribution fullfills 
    910       !    volume and area conservation, positivity and ice categories bounds 
     892      INTEGER  ::   idim 
     893      REAL(wp) ::   zv, zdh 
     894      !!------------------------------------------------------------------- 
     895      ! 
    911896      idim = SIZE( phti , 1 ) 
    912897      ! 
     
    915900      pa_i(1:idim,1:jpl) = 0._wp 
    916901      ! 
     902      ALLOCATE( z1_hti(idim) ) 
     903      WHERE( phti(:) /= 0._wp )   ;   z1_hti(:) = 1._wp / phti(:) 
     904      ELSEWHERE                   ;   z1_hti(:) = 0._wp 
     905      END WHERE 
     906      ! 
     907      ! == thickness and concentration == ! 
     908      ! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl) 
     909      DO jl = 1, jpl-1 
     910         DO ji = 1, idim 
     911            ! 
     912            IF( phti(ji) > 0._wp ) THEN 
     913               ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 
     914               pa_i(ji,jl) = pati(ji) * z1_hti(ji) * (  ( phti(ji) + 2.*hi_max(jl-1) ) * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 
     915                  &                                   - ( phti(ji) + 2.*hi_max(jl  ) ) * EXP( -2.*hi_max(jl  )*z1_hti(ji) ) ) 
     916               ! 
     917               ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 
     918               zv = pati(ji) * z1_hti(ji) * (  ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl-1) + 2.*hi_max(jl-1)*hi_max(jl-1) ) & 
     919                  &                            * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 
     920                  &                          - ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl) + 2.*hi_max(jl)*hi_max(jl) ) & 
     921                  &                            * EXP(-2.*hi_max(jl)*z1_hti(ji)) ) 
     922               ! thickness 
     923               IF( pa_i(ji,jl) > epsi06 ) THEN 
     924                  ph_i(ji,jl) = zv / pa_i(ji,jl) 
     925               ELSE 
     926                  ph_i(ji,jl) = 0. 
     927                  pa_i(ji,jl) = 0. 
     928               ENDIF 
     929            ENDIF 
     930            ! 
     931         ENDDO 
     932      ENDDO 
     933      ! 
     934      ! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity 
    917935      DO ji = 1, idim 
    918936         ! 
    919937         IF( phti(ji) > 0._wp ) THEN 
    920             ! 
    921             ! find which category (jl0) the input ice thickness falls into 
    922             jl0 = jpl 
    923             DO jl = 1, jpl 
    924                IF ( ( phti(ji) >= hi_max(jl-1) ) .AND. ( phti(ji) < hi_max(jl) ) ) THEN 
    925                   jl0 = jl 
    926                   CYCLE 
    927                ENDIF 
    928             END DO 
    929             ! 
    930             itest(:) = 0 
    931             i_fill   = jpl + 1                                            !------------------------------------ 
    932             DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )   ! iterative loop on i_fill categories 
    933                !                                                          !------------------------------------ 
    934                i_fill = i_fill - 1 
    935                ! 
    936                ph_i(ji,1:jpl) = 0._wp 
    937                pa_i(ji,1:jpl) = 0._wp 
    938                itest(:)       = 0       
    939                ! 
    940                IF ( i_fill == 1 ) THEN      !-- case very thin ice: fill only category 1 
    941                   ph_i(ji,1) = phti(ji) 
    942                   pa_i(ji,1) = pati (ji) 
    943                ELSE                         !-- case ice is thicker: fill categories >1 
    944                   ! thickness 
    945                   DO jl = 1, i_fill - 1 
    946                      ph_i(ji,jl) = hi_mean(jl) 
    947                   END DO 
    948                   ! 
    949                   ! concentration 
    950                   pa_i(ji,jl0) = pati(ji) / SQRT(REAL(jpl)) 
    951                   DO jl = 1, i_fill - 1 
    952                      IF ( jl /= jl0 ) THEN 
    953                         zarg        = ( ph_i(ji,jl) - phti(ji) ) / ( phti(ji) * 0.5_wp ) 
    954                         pa_i(ji,jl) =   pa_i (ji,jl0) * EXP(-zarg**2) 
    955                      ENDIF 
    956                   END DO 
    957                   ! 
    958                   ! last category 
    959                   pa_i(ji,i_fill) = pati(ji) - SUM( pa_i(ji,1:i_fill-1) ) 
    960                   zV = SUM( pa_i(ji,1:i_fill-1) * ph_i(ji,1:i_fill-1) ) 
    961                   ph_i(ji,i_fill) = ( phti(ji) * pati(ji) - zV ) / MAX( pa_i(ji,i_fill), epsi10 )  
    962                   ! 
    963                   ! correction if concentration of upper cat is greater than lower cat 
    964                   !    (it should be a gaussian around jl0 but sometimes it is not) 
    965                   IF ( jl0 /= jpl ) THEN 
    966                      DO jl = jpl, jl0+1, -1 
    967                         IF ( pa_i(ji,jl) > pa_i(ji,jl-1) ) THEN 
    968                            zdv = ph_i(ji,jl) * pa_i(ji,jl) 
    969                            ph_i(ji,jl    ) = 0._wp 
    970                            pa_i (ji,jl    ) = 0._wp 
    971                            pa_i (ji,1:jl-1) = pa_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * phti(ji), epsi10 ) 
    972                         END IF 
    973                      END DO 
    974                   ENDIF 
    975                   ! 
    976                ENDIF 
    977                ! 
    978                ! Compatibility tests 
    979                zconv = ABS( pati(ji) - SUM( pa_i(ji,1:jpl) ) )  
    980                IF ( zconv < epsi06 )   itest(1) = 1                                        ! Test 1: area conservation 
    981                ! 
    982                zconv = ABS( phti(ji)*pati(ji) - SUM( pa_i(ji,1:jpl)*ph_i(ji,1:jpl) ) ) 
    983                IF ( zconv < epsi06 )   itest(2) = 1                                        ! Test 2: volume conservation 
    984                ! 
    985                IF ( ph_i(ji,i_fill) >= hi_max(i_fill-1) )   itest(3) = 1                   ! Test 3: thickness of the last category is in-bounds ? 
    986                ! 
    987                itest(4) = 1 
    988                DO jl = 1, i_fill 
    989                   IF ( pa_i(ji,jl) < 0._wp ) itest(4) = 0                                  ! Test 4: positivity of ice concentrations 
    990                END DO 
    991                !                                         !---------------------------- 
    992             END DO                                       ! end iteration on categories 
    993             !                                            !---------------------------- 
     938            ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity 
     939            pa_i(ji,jpl) = pati(ji) * z1_hti(ji) * ( phti(ji) + 2.*hi_max(jpl-1) ) * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 
     940 
     941            ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity 
     942            zv = pati(ji) * z1_hti(ji) * ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jpl-1) + 2.*hi_max(jpl-1)*hi_max(jpl-1) ) & 
     943               &                         * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 
     944            ! thickness 
     945            IF( pa_i(ji,jpl) > epsi06 ) THEN 
     946               ph_i(ji,jpl) = zv / pa_i(ji,jpl) 
     947            else 
     948               ph_i(ji,jpl) = 0. 
     949               pa_i(ji,jpl) = 0. 
     950            ENDIF 
    994951         ENDIF 
    995       END DO 
    996  
     952         ! 
     953      ENDDO 
     954      ! 
    997955      ! Add Snow in each category where pa_i is not 0 
    998956      DO jl = 1, jpl 
    999957         DO ji = 1, idim 
    1000958            IF( pa_i(ji,jl) > 0._wp ) THEN 
    1001                ph_s(ji,jl) = ph_i(ji,jl) * ( phts(ji) / phti(ji) ) 
     959               ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji) 
    1002960               ! In case snow load is in excess that would lead to transformation from snow to ice 
    1003961               ! Then, transfer the snow excess into the ice (different from icethd_dh) 
     
    1009967         END DO 
    1010968      END DO 
     969      ! 
     970      DEALLOCATE( z1_hti ) 
    1011971      ! 
    1012972      ! == temperature and salinity == ! 
Note: See TracChangeset for help on using the changeset viewer.