Ignore:
Timestamp:
2017-09-05T19:53:41+02:00 (3 years ago)
Author:
clem
Message:

changes in style - part2 -

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90

    r8486 r8498  
    4848   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging 
    4949   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting 
    50  
     50   ! 
    5151   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier 
    5252   REAL(wp), PARAMETER ::   kraft   = 0.5_wp    ! rafting multipliyer 
    53  
    54 !!gm Cp is  1) not DOCTOR,  
    55 !!          2) misleading name: heat capacity instead of a constant, 
    56 !!          3) recomputed at each time-step, whereas it is stored in the module memory... 
    57 !!             ===>>> compute it one for all inside the IF( kt == nit000 ) (i.e. without the ".AND. lwp") 
    58    REAL(wp)            ::   Cp                  ! ??? !!gm  Not doctor !  
    59     
     53   ! 
     54   REAL(wp)            ::   zdrho               !  
    6055   ! 
    6156   ! 
     
    111106      INTEGER  ::   niter              ! local integer  
    112107      INTEGER  ::   iterate_ridging    ! if =1, repeat the ridging 
    113       REAL(wp) ::   za, zfac, zcs_2    ! local scalar 
    114       CHARACTER (len = 15) ::   fieldid 
     108      REAL(wp) ::   za                 ! local scalar 
    115109      REAL(wp), DIMENSION(jpi,jpj) ::   closing_net     ! net rate at which area is removed    (1/s) 
    116110      !                                                 ! (ridging ice area - area of new ridges) / dt 
     
    125119      IF( nn_timing == 1 )  CALL timing_start('icerdgrft') 
    126120 
    127       IF( kt == nit000 .AND. lwp ) THEN 
    128          WRITE(numout,*) 
    129          WRITE(numout,*)'icerdgrft : ice ridging and rafting' 
    130          WRITE(numout,*)'~~~~~~~~~~' 
    131       ENDIF 
    132 !!gm should be:       
    133 !      IF( kt == nit000 ) THEN 
    134 !         IF(lwp) WRITE(numout,*) 
    135 !         IF(lwp) WRITE(numout,*)'icerdgrft : ???' 
    136 !         IF(lwp) WRITE(numout,*)'~~~~~~~~~~' 
    137 !         ! 
    138 !         Cp    = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE 
    139 !         ! 
    140 !!gm why not adding here zcs_2 computation 
    141 !         ! 
    142 !      ENDIF 
    143 !!gm end 
    144        
     121      IF( kt == nit000 ) THEN 
     122         IF(lwp) WRITE(numout,*) 
     123         IF(lwp) WRITE(numout,*)'icerdgrft : ice ridging and rafting' 
     124         IF(lwp) WRITE(numout,*)'~~~~~~~~~~' 
     125         ! 
     126         zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE 
     127         ! 
     128      ENDIF       
    145129      !                    ! conservation test 
    146130      IF( ln_limdiachk )   CALL ice_cons_hsm(0, 'icerdgrft', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     
    149133      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    150134      !-----------------------------------------------------------------------------! 
    151       Cp    = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0      ! proport const for PE 
    152       zcs_2 = rn_cs * 0.5_wp 
    153135      ! 
    154136      CALL ice_rdgrft_ridgeprep                             ! prepare ridging 
     
    175157            !  (thick, newly ridged ice). 
    176158 
    177             closing_net(ji,jj) = zcs_2 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
     159            closing_net(ji,jj) = rn_cs * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
    178160 
    179161            ! 2.2 divu_adv 
     
    220202               za   = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice 
    221203               IF    ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN                  ! would lead to negative ato_i 
    222                   zfac          = - ato_i(ji,jj) / za 
    223204                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice  
    224205               ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum 
    225                   zfac          = ( asum(ji,jj) - ato_i(ji,jj) ) / za 
    226206                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice  
    227207               ENDIF 
     
    238218                  za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
    239219                  IF( za  >  a_i(ji,jj,jl) ) THEN 
    240                      zfac = a_i(ji,jj,jl) / za 
    241                      closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 
     220                     closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / za 
    242221                  ENDIF 
    243222               END DO 
     
    254233         asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 
    255234 
    256          ! 3.5 Do we keep on iterating ??? 
     235         ! 3.5 Do we keep on iterating? 
    257236         !-----------------------------------------------------------------------------! 
    258237         ! Check whether asum = 1.  If not (because the closing and opening 
     
    316295      !!---------------------------------------------------------------------! 
    317296      INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    318       REAL(wp) ::   Gstari, astari, hrmean, zdummy   ! local scalar     !!gm DOCTOR norme should start with z !!!!! 
    319       REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n 
     297      REAL(wp) ::   z1_gstar, z1_astar, zhmean, zdummy   ! local scalar 
     298      REAL(wp), DIMENSION(jpi,jpj,-1:jpl) ::   zGsum     ! zGsum(n) = sum of areas in categories 0 to n 
    320299      !------------------------------------------------------------------------------! 
    321300 
    322       Gstari        = 1._wp / rn_gstar     
    323       astari        = 1._wp / rn_astar     
    324       aksum(:,:)    = 0._wp 
    325       athorn(:,:,:) = 0._wp 
    326       aridge(:,:,:) = 0._wp 
    327       araft (:,:,:) = 0._wp 
     301      z1_gstar      = 1._wp / rn_gstar     
     302      z1_astar      = 1._wp / rn_astar     
    328303 
    329304      CALL ice_var_zapsmall   ! Zero out categories with very small areas 
    330305 
    331 !      DO jl = 1, jpl          ! Ice thickness needed for rafting 
    332 !         DO jj = 1, jpj 
    333 !            DO ji = 1, jpi 
    334 !!gm               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 
    335 !!gm               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 
    336 !               IF( a_i(ji,jj,jl) >= epsi20 ) THEN   ;   ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 
    337 !               ELSE                                 ;   ht_i(ji,jj,jl) = 0._wp 
    338 !               ENDIF 
    339 !           END DO 
    340 !         END DO 
    341 !      END DO 
    342 !!gm or even better : 
    343 !     !                       ! Ice thickness needed for rafting 
     306      !                       ! Ice thickness needed for rafting 
    344307      WHERE( a_i(:,:,:) >= epsi20 )   ;   ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:) 
    345308      ELSEWHERE                       ;   ht_i(:,:,:) = 0._wp 
    346309      END WHERE 
    347 !!gm end 
    348310 
    349311      !------------------------------------------------------------------------------! 
     
    356318      ! 
    357319      ! Compute cumulative thickness distribution function 
    358       ! Compute the cumulative thickness distribution function Gsum, 
    359       ! where Gsum(n) is the fractional area in categories 0 to n. 
     320      ! Compute the cumulative thickness distribution function zGsum, 
     321      ! where zGsum(n) is the fractional area in categories 0 to n. 
    360322      ! initial value (in h = 0) equals open water area 
    361       Gsum(:,:,-1) = 0._wp 
    362       Gsum(:,:,0 ) = ato_i(:,:) 
    363       ! for each value of h, you have to add ice concentration then 
     323      zGsum(:,:,-1) = 0._wp 
     324      zGsum(:,:,0 ) = ato_i(:,:) / asum(:,:) 
    364325      DO jl = 1, jpl 
    365          Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 
    366       END DO 
    367       ! 
    368       ! Normalize the cumulative distribution to 1 
    369       DO jl = 0, jpl 
    370          Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:) 
     326         zGsum(:,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:) ! sum(1:jl) is ok (and not jpl) 
    371327      END DO 
    372328 
     
    388344            DO jj = 1, jpj  
    389345               DO ji = 1, jpi 
    390                   IF    ( Gsum(ji,jj,jl)   < rn_gstar ) THEN 
    391                      athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * & 
    392                         &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari ) 
    393                   ELSEIF( Gsum(ji,jj,jl-1) < rn_gstar ) THEN 
    394                      athorn(ji,jj,jl) = Gstari * ( rn_gstar       - Gsum(ji,jj,jl-1) ) *  & 
    395                         &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + rn_gstar       ) * Gstari ) 
     346                  IF    ( zGsum(ji,jj,jl)   < rn_gstar ) THEN 
     347                     athorn(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * & 
     348                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar ) 
     349                  ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN 
     350                     athorn(ji,jj,jl) = z1_gstar * ( rn_gstar        - zGsum(ji,jj,jl-1) ) *  & 
     351                        &                          ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar        ) * z1_gstar ) 
    396352                  ELSE 
    397353                     athorn(ji,jj,jl) = 0._wp 
     
    403359      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007) 
    404360         !                         
    405          zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array 
     361         zdummy = 1._wp / ( 1._wp - EXP(-z1_astar) )        ! precompute exponential terms using zGsum as a work array 
    406362         DO jl = -1, jpl 
    407             Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 
     363            zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy 
    408364         END DO 
    409365         DO jl = 0, jpl 
    410              athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 
     366            athorn(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl) 
    411367         END DO 
    412368         ! 
     
    418374            DO jj = 1, jpj  
    419375               DO ji = 1, jpi 
    420                   zdummy           = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) 
    421                   aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl) 
     376                  aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * athorn(ji,jj,jl) 
    422377                  araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl) 
    423378               END DO 
    424379            END DO 
    425380         END DO 
    426       ELSEIF( ln_ridging .AND. .NOT.ln_rafting ) THEN   !- ridging alone 
    427          DO jl = 1, jpl 
    428             aridge(:,:,jl) = athorn(:,:,jl) 
    429          END DO 
    430       ELSEIF( ln_rafting .AND. .NOT.ln_ridging ) THEN   !- rafting alone    
    431          DO jl = 1, jpl 
    432             araft(:,:,jl) = athorn(:,:,jl) 
    433          END DO 
     381      ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN   !- ridging alone 
     382         aridge(:,:,:) = athorn(:,:,1:jpl) 
     383         araft (:,:,:) = 0._wp 
     384      ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN   !- rafting alone    
     385        aridge(:,:,:) = 0._wp 
     386         araft (:,:,:) = athorn(:,:,1:jpl) 
     387      ELSE                                               !- no ridging & no rafting 
     388         aridge(:,:,:) = 0._wp 
     389         araft (:,:,:) = 0._wp          
    434390      ENDIF 
    435391 
     
    441397      !  
    442398      ! This parameterization is a modified version of Hibler (1980). 
    443       ! The mean ridging thickness, hrmean, is proportional to hi^(0.5) 
     399      ! The mean ridging thickness, zhmean, is proportional to hi^(0.5) 
    444400      !  and for very thick ridging ice must be >= krdgmin*hi 
    445401      ! 
    446402      ! The minimum ridging thickness, hrmin, is equal to 2*hi  
    447403      !  (i.e., rafting) and for very thick ridging ice is 
    448       !  constrained by hrmin <= (hrmean + hi)/2. 
     404      !  constrained by hrmin <= (zhmean + hi)/2. 
    449405      !  
    450406      ! The maximum ridging thickness, hrmax, is determined by 
    451       !  hrmean and hrmin. 
     407      !  zhmean and hrmin. 
    452408      ! 
    453409      ! These modifications have the effect of reducing the ice strength 
     
    466422            DO ji = 1, jpi 
    467423               IF ( athorn(ji,jj,jl) > 0._wp ) THEN 
    468                   hrmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin ) 
    469                   hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) ) 
    470                   hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl) 
     424                  zhmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin ) 
     425                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) ) 
     426                  hrmax(ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl) 
    471427                  hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft 
    472                   krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 ) 
     428                  krdg (ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 ) 
    473429                  ! 
    474430                  ! Normalization factor : aksum, ensures mass conservation 
     
    498454      !!---------------------------------------------------------------------- 
    499455      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear 
    500       REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges 
    501       ! 
    502       CHARACTER (len=80) ::   fieldid   ! field identifier     !!gm DOCTOR name wrong 
     456      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area retreats, excluding area of new ridges 
    503457      ! 
    504458      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices 
    505459      INTEGER ::   ij                ! horizontal index, combines i and j loops 
    506460      INTEGER ::   icells            ! number of cells with a_i > puny 
    507       REAL(wp) ::   hL, hR, farea    ! left and right limits of integration 
    508       REAL(wp) ::   zwfx_snw         ! snow mass flux increment 
     461      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration and new area going to jl2 
    509462 
    510463      INTEGER , DIMENSION(jpij) ::   indxi, indxj   ! compressed indices 
    511       REAL(wp), DIMENSION(jpij) ::   zswitch, fvol   ! new ridge volume going to n2 
     464      REAL(wp), DIMENSION(jpij) ::   zswitch, fvol   ! new ridge volume going to jl2 
    512465 
    513466      REAL(wp), DIMENSION(jpij) ::   afrac            ! fraction of category area ridged  
     
    656609            !  During the next time step, thermo_rates will determine whether 
    657610            !  the ocean cools or new ice grows. 
    658             zwfx_snw           = ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   &  
    659                &                 + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean 
    660  
    661             wfx_snw_dyn(ji,jj)  =  wfx_snw_dyn(ji,jj) + zwfx_snw 
     611            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   &  
     612               &                                      + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean 
     613 
     614!!!            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + zwfx_snw 
    662615 
    663616            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         &  
     
    722675 
    723676               ! update jl1 
    724                e_i  (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk) 
     677               e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk) 
    725678 
    726679            END DO 
     
    806759      ! 
    807760      INTEGER             ::   ji,jj, jl   ! dummy loop indices 
    808       INTEGER             ::   ksmooth     ! smoothing the resistance to deformation    !!gm not DOCTOR : start with i !!! 
    809       INTEGER             ::   numts_rm    ! number of time steps for the P smoothing    !!gm not DOCTOR : start with i !!! 
     761      INTEGER             ::   ismooth     ! smoothing the resistance to deformation 
     762      INTEGER             ::   itframe     ! number of time steps for the P smoothing 
    810763      REAL(wp)            ::   zp, z1_3    ! local scalars 
    811764      REAL(wp), DIMENSION(jpi,jpj) ::   zworka           ! temporary array used here 
     
    813766      !!---------------------------------------------------------------------- 
    814767 
    815       !------------------------------------------------------------------------------! 
    816       ! 1) Initialize 
    817       !------------------------------------------------------------------------------! 
    818       strength(:,:) = 0._wp 
    819  
    820       !------------------------------------------------------------------------------! 
    821       ! 2) Compute thickness distribution of ridging and ridged ice 
    822       !------------------------------------------------------------------------------! 
    823       CALL ice_rdgrft_ridgeprep 
    824  
    825       !------------------------------------------------------------------------------! 
    826       ! 3) Rothrock(1975)'s method 
    827       !------------------------------------------------------------------------------! 
    828       IF( kstrngth == 1 ) THEN 
     768      !                              !--------------------------------------------------! 
     769      CALL ice_rdgrft_ridgeprep      ! Thickness distribution of ridging and ridged ice ! 
     770      !                              !--------------------------------------------------! 
     771 
     772      !                              !--------------------------------------------------! 
     773      IF( kstrngth == 1 ) THEN       ! Ice strength => Rothrock (1975) method           ! 
     774      !                              !--------------------------------------------------! 
    829775         z1_3 = 1._wp / 3._wp 
    830776         DO jl = 1, jpl 
    831             DO jj= 1, jpj 
    832                DO ji = 1, jpi 
    833                   ! 
    834                   IF( athorn(ji,jj,jl) > 0._wp ) THEN 
    835                      !---------------------------- 
    836                      ! PE loss from deforming ice 
    837                      !---------------------------- 
    838                      strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl) 
    839  
    840                      !-------------------------- 
    841                      ! PE gain from rafting ice 
    842                      !-------------------------- 
    843                      strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl) 
    844  
    845                      !---------------------------- 
    846                      ! PE gain from ridging ice 
    847                      !---------------------------- 
    848                      strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 *  & 
    849                         &                              ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) +         & 
    850                         &                                hrmin(ji,jj,jl) * hrmin(ji,jj,jl) +         & 
    851                         &                                hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )   
    852                         !!(a**3-b**3)/(a-b) = a*a+ab+b*b                       
    853                   ENDIF 
    854                   ! 
    855                END DO 
    856             END DO 
    857          END DO 
    858     
    859          strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1) 
    860                          ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 
    861          ksmooth = 1 
    862  
    863       !------------------------------------------------------------------------------! 
    864       ! 4) Hibler (1979)' method 
    865       !------------------------------------------------------------------------------! 
    866       ELSE                      ! kstrngth ne 1:  Hibler (1979) form 
    867          ! 
     777            WHERE( athorn(:,:,jl) > 0._wp ) 
     778               strength(:,:) =  -         athorn(:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE loss from deforming ice 
     779                  &             + 2._wp * araft (:,:,jl) * ht_i(:,:,jl) * ht_i(:,:,jl)   &  ! PE gain from rafting ice 
     780                  &             +         aridge(:,:,jl) * krdg(:,:,jl) * z1_3 *   &        ! PE gain from ridging ice 
     781                  &                      ( hrmax(:,:,jl) * hrmax(:,:,jl) +         & 
     782                  &                        hrmin(:,:,jl) * hrmin(:,:,jl) +         & 
     783                  &                        hrmax(:,:,jl) * hrmin(:,:,jl) ) 
     784            ELSEWHERE 
     785               strength(:,:) = 0._wp 
     786            END WHERE 
     787         END DO 
     788         strength(:,:) = rn_pe_rdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1) 
     789                         ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 
     790         ismooth = 1 
     791      !                              !--------------------------------------------------! 
     792      ELSE                           ! Ice strength => Hibler (1979) method             ! 
     793      !                              !--------------------------------------------------! 
    868794         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) * tmask(:,:,1) 
    869795         ! 
    870          ksmooth = 1 
    871          ! 
    872       ENDIF                     ! kstrngth 
    873       ! 
    874       !------------------------------------------------------------------------------! 
    875       ! 5) Impact of brine volume 
    876       !------------------------------------------------------------------------------! 
     796         ismooth = 1 
     797         ! 
     798      ENDIF 
     799      !                              !--------------------------------------------------! 
     800      IF( ln_icestr_bvf ) THEN       ! Impact of brine volume                           ! 
     801      !                              !--------------------------------------------------! 
    877802      ! CAN BE REMOVED 
    878       IF( ln_icestr_bvf ) THEN 
    879803         DO jj = 1, jpj 
    880804            DO ji = 1, jpi 
     
    883807         END DO 
    884808      ENDIF 
    885       ! 
    886       !------------------------------------------------------------------------------! 
    887       ! 6) Smoothing ice strength 
    888       !------------------------------------------------------------------------------! 
    889       SELECT CASE( ksmooth ) 
    890       !                       !------------------- 
    891       CASE( 1 )               ! Spatial smoothing 
    892          !                    !------------------- 
     809      !                              !--------------------------------------------------! 
     810      SELECT CASE( ismooth )         ! Smoothing ice strength                           ! 
     811      !                              !--------------------------------------------------! 
     812      ! 
     813      CASE( 1 )               !--- Spatial smoothing 
    893814         DO jj = 2, jpjm1 
    894815            DO ji = 2, jpim1 
     
    903824            END DO 
    904825         END DO 
    905  
     826          
    906827         DO jj = 2, jpjm1 
    907828            DO ji = 2, jpim1 
     
    911832         CALL lbc_lnk( strength, 'T', 1. ) 
    912833         ! 
    913          !                    !-------------------- 
    914       CASE( 2 )               ! Temporal smoothing 
    915          !                    !-------------------- 
     834      CASE( 2 )               !--- Temporal smoothing 
    916835         IF ( kt_ice == nit000 ) THEN 
    917836            zstrp1(:,:) = 0._wp 
     
    922841            DO ji = 2, jpim1 
    923842               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN  
    924                   numts_rm = 1 ! number of time steps for the running mean 
    925                   IF ( zstrp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 
    926                   IF ( zstrp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 
    927                   zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / numts_rm 
     843                  itframe = 1 ! number of time steps for the running mean 
     844                  IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1 
     845                  IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1 
     846                  zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe 
    928847                  zstrp2  (ji,jj) = zstrp1  (ji,jj) 
    929848                  zstrp1  (ji,jj) = strength(ji,jj) 
     
    932851            END DO 
    933852         END DO 
    934          CALL lbc_lnk( strength, 'T', 1. )      ! Boundary conditions 
     853         CALL lbc_lnk( strength, 'T', 1. ) 
    935854         ! 
    936855      END SELECT 
     
    970889      IF (lwp) THEN                          ! control print 
    971890         WRITE(numout,*) 
    972          WRITE(numout,*)'ice_rdgrft_init : ice parameters for mechanical ice redistribution ' 
    973          WRITE(numout,*)'~~~~~~~~~~~~~~~' 
    974          WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs  
    975          WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun 
    976          WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar 
    977          WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar 
    978          WRITE(numout,*)'   Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging 
    979          WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar 
    980          WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg 
    981          WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg  
    982          WRITE(numout,*)'   Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg  
    983          WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting 
    984          WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft 
    985          WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft   
    986          WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft  
    987          WRITE(numout,*)'   Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft  
     891         WRITE(numout,*) 'ice_rdgrft_init : ice parameters for mechanical ice redistribution ' 
     892         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
     893         WRITE(numout,*) '   Namelist namiceitdme' 
     894         WRITE(numout,*) '      Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs  
     895         WRITE(numout,*) '      Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun 
     896         WRITE(numout,*) '      Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar 
     897         WRITE(numout,*) '      Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar 
     898         WRITE(numout,*) '      Ridging of ice sheets or not                            ln_ridging  = ', ln_ridging 
     899         WRITE(numout,*) '      Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar 
     900         WRITE(numout,*) '      Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg 
     901         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg  
     902         WRITE(numout,*) '      Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg  
     903         WRITE(numout,*) '      Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting 
     904         WRITE(numout,*) '      Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft 
     905         WRITE(numout,*) '      Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft   
     906         WRITE(numout,*) '      Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft  
     907         WRITE(numout,*) '      Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft  
    988908      ENDIF 
    989909      ! 
Note: See TracChangeset for help on using the changeset viewer.