New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 15507 – NEMO

Changeset 15507


Ignore:
Timestamp:
2021-11-15T15:22:24+01:00 (6 months ago)
Author:
emmafiedler
Message:

Working version of Rothrock ice strength

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/NEMO_4.0.4_ice_strength/src/ICE/icedyn_rdgrft.F90

    r15346 r15507  
    4343   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   opning          ! rate of opening due to divergence/shear 
    4444   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
    45    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   zaksum          ! normalisation factor    
     45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   zaksum          ! normalisation factor 
    4646   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:)     ::   zasum           ! sum of ice area plus open water    
    4747   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   apartf          ! participation function; fraction of ridging/closing associated w/ category n 
     
    4949   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmin           ! minimum ridge thickness 
    5050   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrmax           ! maximum ridge thickness 
    51    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hrexp           ! e-folding ridge thickness 
    5251   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hraft           ! thickness of rafted ice 
    5352   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hi_hrdg         ! thickness of ridging ice / mean ridge thickness 
     
    9594      ALLOCATE( closing_net(jpij)  , opning(jpij)      , closing_gross(jpij) ,  & 
    9695         &      apartf(jpij,0:jpl) , hrmin(jpij,jpl)   , hrmax(jpij,jpl)     , zhi(jpij,jpl),  & 
    97          &      zaksum(jpij)       , hraft(jpij,jpl)   , hrexp(jpij,jpl)     , & 
    98          &      hi_hrdg(jpij,jpl)  , araft(jpij,jpl)   , aridge(jpij,jpl)    ,  zasum(jpij), & 
     96         &      zaksum(jpij)       , hraft(jpij,jpl)   , & 
     97         &      hi_hrdg(jpij,jpl)  , araft(jpij,jpl)   , aridge(jpij,jpl)    , zasum(jpij), & 
    9998         &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 
    10099 
     
    146145      INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not 
    147146      REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i 
    148       REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! 1D divu_i as implied by transport scheme 
    149147      ! 
    150148      INTEGER, PARAMETER ::   jp_itermax = 20     
     
    194192            !                                                        - ice area added in new ridges 
    195193            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
    196             ! 
    197  !           zdivu_adv(ji) = (1._wp - zasum(ji)) / rdt_ice  ! Need to do more testing/checking on this 
    198  
    199 !            WRITE(numout,*)'velocity div ',zdivu(ji), 'transport div ',zdivu_adv(ji) 
    200       
     194            !      
    201195            IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough 
    202 !              IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) 
    203  
    204196            !                                                                                ! to give asum = 1.0 after ridging 
    205197            ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 
    206198            opning(ji) = closing_net(ji) + zdivu(ji) 
    207 !            opning(ji) = closing_net(ji) + zdivu_adv(ji) 
    208  
     199            ! 
    209200         END DO 
    210201         ! 
     
    310301      INTEGER  ::   ji, jl                     ! dummy loop indices 
    311302      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zfac   ! local scalar 
    312 !      REAL(wp), DIMENSION(jpij)        ::   zasum, z1_asum           ! sum of a_i+ato_i and inverse  
    313303      REAL(wp), DIMENSION(jpij)        ::   z1_asum                  ! inverse of sum of a_i+ato_i 
    314304      REAL(wp), DIMENSION(jpij,-1:jpl) ::   zGsum                    ! zGsum(n) = sum of areas in categories 0 to n 
     
    324314      ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp 
    325315      END WHERE 
    326       ! Make sure ice thickness is not below minimum in any category (from iceitd.F90; no minimum was causing issues for R75 strength) 
     316      ! Make sure ice thickness is not below the minimum (from iceitd.F90; no minimum was causing issues for R75  
     317      ! strength) 
    327318      WHERE( pa_i(1:npti,:) > epsi10 .AND. zhi(1:npti,:) < rn_himin )  
    328319        pa_i(1:npti,:) = pa_i(1:npti,:) * zhi(1:npti,:) / rn_himin 
    329320        zhi(1:npti,:) = rn_himin 
     321      END WHERE 
     322 
     323      ! Set tiny values of pa_i to zero 
     324      ! If this isn't done, can end up with zhi=0 but apartf>0 
     325      WHERE( pa_i(1:npti,:) <= epsi10 )   ;   pa_i(1:npti,:) = 0._wp 
    330326      END WHERE 
    331327 
     
    344340      ! Compute total area of ice plus open water. 
    345341      ! This is in general not equal to one because of divergence during transport 
    346      ! zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 
    347      ! EF: to match CICE, this should be the following. What is pa_i(:,0)?? Very small. 
    348       zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,1:jpl), dim=2 )    
     342      zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,1:jpl), dim=2 ) 
    349343 
    350344      WHERE( zasum(1:npti) > epsi10 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti) 
     
    359353      zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti) 
    360354      DO jl = 1, jpl 
    361          zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti)  ! sum(1:jl) is ok (and not jpl)  
     355         zGsum(1:npti,jl) = ( pato_i(1:npti) + SUM( pa_i(1:npti,1:jl), dim=2 ) ) * z1_asum(1:npti)  ! sum(1:jl) is correct (and not jpl)  
    362356      END DO 
    363357      ! 
     
    457451               hrmin  (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) ) 
    458452               hrmax  (ji,jl) = 2._wp * zhmean - hrmin(ji,jl) 
    459                hrexp  (ji,jl) = rn_murdg * SQRT( zhi(ji,jl) ) 
    460453               hraft  (ji,jl) = zhi(ji,jl) * zfac 
    461                hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 ) 
     454               hi_hrdg(ji,jl) = zhi(ji,jl) / zhmean 
    462455               ! 
    463456               ! Normalization factor : zaksum, ensures mass conservation 
     
    467460               hrmin  (ji,jl) = 0._wp  
    468461               hrmax  (ji,jl) = 0._wp  
    469                hrexp  (ji,jl) = 0._wp  
    470462               hraft  (ji,jl) = 0._wp  
    471463               hi_hrdg(ji,jl) = 1._wp 
     
    788780      !!              dissipated per unit area removed from the ice pack under compression, 
    789781      !!              and assumed proportional to the change in potential energy caused 
    790       !!              by ridging. [**UPDATE THIS ---> Note that only Hibler's formulation is stable and that 
    791       !!              ice strength has to be smoothed**] 
     782      !!              by ridging. Note that ice strength using Hibler's formulation must be 
     783      !!              smoothed. 
    792784      !!---------------------------------------------------------------------- 
    793785      INTEGER             ::   ji, jj, jl  ! dummy loop indices 
     
    796788      REAL(wp)            ::   zp, z1_3    ! local scalars 
    797789      REAL(wp)            ::   Cp          ! proportional constant for PE 
    798       !REAL(wp)            ::   hi          ! ice thickness of category (m) 
    799790      REAL(wp)            ::   h2rdg       ! mean value of h^2 for new ridge 
    800791      REAL(wp)            ::   dh2rdg      ! change in mean value of h^2 per unit area consumed by ridging  
    801792      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here 
    802793      REAL(wp), DIMENSION(jpi,jpj) ::   zstrp1, zstrp2   ! strength at previous time steps 
    803       REAL(wp), DIMENSION(jpij,jpl) ::   krdg             ! mean ridge thickness / thickness of ridging ice 
    804794      REAL(wp), DIMENSION(jpij)     ::   strength_1d      ! 1-dimensional strength 
    805795      REAL(wp), DIMENSION(jpij,jpl) ::   strength_pe      ! PE component of R75 strength 
     
    809799      IF( ln_str_H79 ) THEN          ! Ice strength => Hibler (1979) method             ! 
    810800      !                              !--------------------------------------------------! 
    811          strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 
    812          ismooth = 1 
     801        strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 
     802        ismooth = 1 
    813803      !                              !--------------------------------------------------! 
    814804      ELSE IF( ln_str_R75 ) THEN     ! Ice strength => Rothrock (1975) method           ! 
    815       !                              !--------------------------------------------------!  
    816       IF ( kt_ice == nit000 ) THEN 
    817  
    818         ! Replace this with input from start file 
    819         strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 
    820         ismooth = 1 
    821  
    822       ELSE 
    823  
    824       !-------------------------------- 
    825       ! 0) Identify grid cells with ice 
    826       !-------------------------------- 
    827       at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    828       vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 
    829       ! 
    830       npti = 0   ;   nptidx(:) = 0 
    831       DO jj = 1, jpj 
    832          DO ji = 1, jpi 
    833             IF ( at_i(ji,jj) > epsi10 ) THEN 
    834                npti           = npti + 1 
    835                nptidx( npti ) = (jj - 1) * jpi + ji 
    836             ENDIF 
    837          END DO 
    838       END DO 
    839  
    840       CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i   ) 
    841       CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i   ) 
    842       CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i ) 
    843  
    844  
    845       !----------------------------------- 
    846       ! Get thickness distribution of ice 
    847       !----------------------------------- 
    848  
    849       CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 
    850        
    851       ! Don't need closing rate here (unless only calculating strength for ridging ice), could split this function up. 
    852  
    853       Cp = 0.5_wp * grav * (rhow-rhoi) * rhoi / rhow  ! put this elsewhere, perhaps ice_prep?   
    854  
    855  
    856       IF ( npti > 0 ) THEN 
    857  
    858       DO jl = 1, jpl              ! categories 
    859         DO ji = 1, npti           ! grid cells 
    860  
    861           IF ( a_i_2d(ji,jl) > epsi10 .and. apartf(ji,jl) > 0._wp ) THEN 
    862  
    863           ! Move this to rdgrft_prep 
    864            krdg(ji,jl) = (hrmin(ji,jl) + hrmax(ji,jl))/2._wp * 1._wp/zhi(ji,jl) 
    865  
    866  
    867          ! Add a logical for distribution    
    868          ! Uniform distribution 
    869            h2rdg = (1._wp/3._wp) * (hrmax(ji,jl)**3._wp - hrmin(ji,jl)**3._wp)  & 
    870                            / (hrmax(ji,jl) - hrmin(ji,jl))         
    871  
    872      
    873          ! Exponential distribution 
    874          !  h2rdg =          hrmin(ji,jl) * hrmin(ji,jl) & 
    875          !         + 2._wp * hrmin(ji,jl) * hrexp(ji,jl) & 
    876          !         + 2._wp * hrexp(ji,jl) * hrexp(ji,jl) 
    877  
    878            dh2rdg = -zhi(ji,jl)*zhi(ji,jl) + h2rdg / krdg(ji,jl) 
    879            
    880            strength_pe(ji,jl) = apartf(ji,jl) * dh2rdg           
    881  
    882           ELSE 
    883  
    884            strength_pe(ji,jl) = 0._wp 
    885  
    886           END IF 
    887  
    888         END DO                    ! grid cells 
    889       END DO                      ! categories 
    890  
    891  
    892       strength_pe_sum(:) = SUM(strength_pe(:,:), dim=2)  ! Sum over all categories 
    893  
    894  
    895       DO ji = 1, npti   
    896          IF ( zaksum(ji) > epsi10 ) THEN   
    897           strength_1d(ji) = rn_Cf * Cp * strength_pe_sum(ji) / zaksum(ji) 
    898         ELSE 
    899           strength_1d(ji) = 0._wp 
    900         END IF 
    901       END DO 
    902           
    903       ! Convert strength_1d to 2d 
    904  
    905       CALL tab_1d_2d( npti, nptidx(1:npti), strength_1d(1:npti), strength(:,:) ) 
    906  
    907       ! Enforce a maximum for strength (- Look into why numbers can be so large) 
    908       WHERE( strength(:,:) > 1.5e5) ; strength(:,:) = 1.5e5 
    909       END WHERE 
    910  
    911       ismooth = 0  ! Buffer size error with lbc_lnk... 
    912  
    913       END IF  ! number of ice points 
    914       END IF  ! timestep 
     805      !                              !--------------------------------------------------! 
     806      ! 
     807        ! Initialise the strength field at each timestep 
     808        strength(:,:) = 0._wp 
     809        ismooth = 0 
     810 
     811        !-------------------------------- 
     812        ! 0) Identify grid cells with ice 
     813        !-------------------------------- 
     814        at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
     815        vt_i(:,:) = SUM( v_i(:,:,:), dim=3 ) 
     816        ato_i(:,:) = 1 - at_i(:,:) 
     817        ! 
     818        npti = 0   ;   nptidx(:) = 0 
     819        DO jj = 1, jpj 
     820           DO ji = 1, jpi 
     821              IF ( at_i(ji,jj) > epsi10 ) THEN 
     822                 npti           = npti + 1 
     823                 nptidx( npti ) = (jj - 1) * jpi + ji 
     824              ENDIF 
     825           END DO 
     826        END DO 
     827 
     828        CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i   ) 
     829        CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i   ) 
     830        CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti)      , ato_i ) 
     831 
     832        !----------------------------------- 
     833        ! Get thickness distribution of ice 
     834        !----------------------------------- 
     835 
     836        Cp = 0.5_wp * grav * (rhow-rhoi) * rhoi / rhow  ! EF: This should go in an initialisation routine but unclear where, won't let me set as a parameter above 
     837 
     838        CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 
     839        ! EF: Note that closing rate is not needed here, could split this function up, and refactor as initial rdgrft_prep is  
     840        ! called here, before going into ridging iterations 
     841 
     842        ! Initialise 
     843        strength_pe_sum(:) = 0._wp 
     844 
     845        IF ( npti > 0 ) THEN 
     846 
     847          DO jl = 1, jpl              ! categories 
     848            DO ji = 1, npti           ! grid cells 
     849              IF ( apartf(ji,jl) > 0._wp ) THEN 
     850 
     851                ! Uniform distribution of ridged ice 
     852                h2rdg = (1._wp/3._wp) * (hrmax(ji,jl)**3._wp - hrmin(ji,jl)**3._wp)  & 
     853                         / (hrmax(ji,jl) - hrmin(ji,jl)) 
     854 
     855                dh2rdg = -zhi(ji,jl)*zhi(ji,jl) + h2rdg * hi_hrdg(ji,jl) 
     856    
     857                strength_pe_sum(ji) = strength_pe_sum(ji) + apartf(ji,jl) * dh2rdg 
     858 
     859              ELSE 
     860 
     861                strength_pe_sum(ji) = strength_pe_sum(ji) + 0._wp 
     862 
     863              END IF 
     864 
     865            END DO                    ! grid cells 
     866          END DO                      ! categories 
     867 
     868 
     869          DO ji = 1, npti 
     870            IF ( zaksum(ji) > epsi10 ) THEN  ! Unclear why is this check should be needed? Occasional x10E-24... 
     871              strength_1d(ji) = rn_Cf * Cp * strength_pe_sum(ji) / zaksum(ji) 
     872            ELSE 
     873              strength_1d(ji) = 0._wp 
     874            END IF 
     875          END DO 
     876 
     877 
     878          ! Enforce a maximum for strength 
     879          WHERE( strength_1d(:) > 2.0e5) ; strength_1d(:) = 2.0e5 
     880          END WHERE 
     881 
     882 
     883          ! Convert strength_1d to 2d 
     884          CALL tab_1d_2d( npti, nptidx(1:npti), strength_1d(1:npti), strength(:,:) ) 
     885 
     886         END IF  ! number of ice points 
    915887 
    916888 
     
    920892         strength(:,:) = 0._wp 
    921893         ismooth = 0 
    922       ENDIF 
     894      ENDIF  ! Ice strength formulation 
    923895      !       
    924  
    925896                                     !--------------------------------------------------! 
    926897      SELECT CASE( ismooth )         ! Smoothing ice strength                           ! 
     
    928899      CASE( 0 )               !--- No smoothing 
    929900 
    930   
    931       !CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 
     901      CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 
    932902 
    933903      CASE( 1 )               !--- Spatial smoothing 
     
    11131083      ENDIF 
    11141084      ! 
     1085      IF ( ( ln_str_H79 .AND. ln_str_R75 ) .OR. ( .NOT.ln_str_H79 .AND. .NOT.ln_str_R75 ) ) THEN 
     1086         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one ice strength formulation (ln_str_H79 or ln_str_R75)' ) 
     1087      ENDIF 
     1088      ! 
    11151089      IF ( ( ln_partf_lin .AND. ln_partf_exp ) .OR. ( .NOT.ln_partf_lin .AND. .NOT.ln_partf_exp ) ) THEN 
    11161090         CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 
Note: See TracChangeset for help on using the changeset viewer.