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 8498 for branches/2017/dev_r8183_ICEMODEL – NEMO

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

changes in style - part2 -

Location:
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
Files:
17 edited

Legend:

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

    r8486 r8498  
    332332 
    333333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total)          [s-1] 
    334    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_thd     !: ice concentration tendency (thermodynamics) [s-1] 
    335    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics)       [s-1] 
    336334 
    337335   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s] 
     
    479477   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dmi_dyn  !: Change in ice mass due to ice dynamics (kg/m2/s) 
    480478   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dms_dyn  !: Change in snow mass due to ice dynamics (kg/m2/s) 
    481    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xmtrp_ice !: X-component of ice mass transport (kg/s) 
    482    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_ymtrp_ice !: Y-component of ice mass transport (kg/s) 
    483    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xmtrp_snw !: X-component of snow mass transport (kg/s) 
    484    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_ymtrp_snw !: Y-component of snow mass transport (kg/s) 
    485    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_xatrp    !: X-component of area transport (m2/s) 
    486    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_yatrp    !: Y-component of area transport (m2/s) 
    487479   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_bo    !: Bottom conduction flux (W/m2) 
    488480   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_su    !: Surface conduction flux (W/m2) 
    489    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_utau_oi  !: X-direction ocean-ice stress 
    490    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_vtau_oi  !: Y-direction ocean-ice stress   
    491    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dssh_dx  !: X-direction sea-surface tilt term (N/m2) 
    492    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dssh_dy  !: X-direction sea-surface tilt term (N/m2) 
    493    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_corstrx  !: X-direction coriolis stress (N/m2) 
    494    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_corstry  !: Y-direction coriolis stress (N/m2) 
    495    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_intstrx  !: X-direction internal stress (N/m2) 
    496    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_intstry  !: Y-direction internal stress (N/m2) 
    497    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_sig1     !: Average normal stress in sea ice    
    498    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_sig2     !: Maximum shear stress in sea ice 
    499481   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_shear    !: Maximum shear of sea-ice velocity field 
    500482 
     
    534516         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     & 
    535517         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,     & 
    536          &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) , rn_amax_2d(jpi,jpj),   & 
     518         &      afx_tot(jpi,jpj) , rn_amax_2d(jpi,jpj),                                         & 
    537519         &      fhtur  (jpi,jpj) , qlead  (jpi,jpj) ,                                           & 
    538520         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  & 
     
    616598      ALLOCATE( t_si (jpi,jpj,jpl)    , tm_si(jpi,jpj)        ,    &  
    617599                diag_dmi_dyn(jpi,jpj) , diag_dms_dyn(jpi,jpj) ,    & 
    618                 diag_xmtrp_ice(jpi,jpj), diag_ymtrp_ice(jpi,jpj),  & 
    619                 diag_xmtrp_snw(jpi,jpj), diag_ymtrp_snw(jpi,jpj),  & 
    620                 diag_xatrp(jpi,jpj)    , diag_yatrp(jpi,jpj)    ,  & 
    621600                diag_fc_bo(jpi,jpj)   , diag_fc_su(jpi,jpj)   ,    & 
    622                 diag_utau_oi(jpi,jpj) , diag_vtau_oi(jpi,jpj) ,    & 
    623                 diag_dssh_dx(jpi,jpj) , diag_dssh_dy(jpi,jpj) ,    & 
    624                 diag_corstrx(jpi,jpj) , diag_corstry(jpi,jpj) ,    & 
    625                 diag_intstrx(jpi,jpj) , diag_intstry(jpi,jpj) ,    & 
    626                 diag_sig1(jpi,jpj)    , diag_sig2(jpi,jpj)    ,    & 
    627601                diag_shear(jpi,jpj)   , STAT = ierr(ii) ) 
    628602 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90

    r8486 r8498  
    2727   REAL(wp), PUBLIC, PARAMETER ::   rn_alb_oce = 0.066   !: ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 
    2828 
    29 !!gm real variable name starting with a "c" NOT DOCTOR !!!! 
    3029   INTEGER  ::   albd_init = 0       ! control flag for initialization   
    31    REAL(wp) , PARAMETER ::   c1        = 0.05    ! snow thickness (only for nn_ice_alb=0) 
    32    REAL(wp) , PARAMETER ::   c2        = 0.10    !  "        " 
    33    REAL(wp) , PARAMETER ::   rcloud    = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0) 
    34    REAL(wp) , PARAMETER ::   r1_c1 = 1. / c1 
    35    REAL(wp) , PARAMETER ::   r1_c2 = 1. / c2 
     30   REAL(wp) , PARAMETER ::   rc1    = 0.05    ! snow thickness (only for nn_ice_alb=0) 
     31   REAL(wp) , PARAMETER ::   rc2    = 0.10    !  "        " 
     32   REAL(wp) , PARAMETER ::   rcloud = 0.06    ! cloud effect on albedo (only-for nn_ice_alb=0) 
     33   REAL(wp) , PARAMETER ::   r1_c1 = 1. / rc1 
     34   REAL(wp) , PARAMETER ::   r1_c2 = 1. / rc2 
    3635  
    3736   !                             !!* namelist namsbc_alb * 
     
    114113      IF ( .NOT. ld_pnd ) THEN   !--- No melt ponds OR radiatively inactive melt ponds 
    115114         ! Bare ice albedo is prescribed, with implicit assumption on pond fraction and depth 
    116          WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = rn_alb_imlt 
    117                                                        ! !!! MV I think we could replace rt0_ice by rt0 and get rid of rt0 
    118          ELSE WHERE                                              ;   zalb(:,:,:) = rn_alb_idry 
    119          END  WHERE 
     115         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0 )   ;   zalb(:,:,:) = rn_alb_imlt 
     116         ELSEWHERE                                           ;   zalb(:,:,:) = rn_alb_idry 
     117         END WHERE 
    120118      ENDIF 
    121119 
     
    127125         ! 
    128126         !                       ! Thickness-dependent bare ice albedo 
    129          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
    130          ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) 
    131          ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              & 
    132             &                                                                 - 0.8608 * ph_ice * ph_ice     & 
    133             &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice 
    134          ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice 
     127         WHERE    ( 1.5  < ph_ice                     )   ;   zalb_it = zalb 
     128         ELSEWHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) 
     129         ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )   ;   zalb_it = 0.2467 + 0.7049 * ph_ice              & 
     130            &                                                                  - 0.8608 * ph_ice * ph_ice     & 
     131            &                                                                  + 0.3812 * ph_ice * ph_ice * ph_ice 
     132         ELSEWHERE                                        ;   zalb_it = 0.1    + 3.6    * ph_ice 
    135133         END WHERE 
    136134         ! 
     
    139137            zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd )  
    140138            !                          ! Snow-free ice albedo is a function of pond fraction 
    141             WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )    
     139            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 )    
    142140               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd  
    143141            END WHERE 
     
    150148               DO ji = 1, jpi 
    151149                  !                    ! Freezing snow 
    152                   ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2 
    153                   zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
     150                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > rc2 
     151                  zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - rc1 ) ) ) 
    154152                  zalb_sf = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  & 
    155153                     &                                   + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * r1_c1  )   & 
     
    157155                     ! 
    158156                  !                    ! Melting snow 
    159                   ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2 
    160                   zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) ) 
     157                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > rc2 
     158                  zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - rc2 ) ) 
    161159                  zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * r1_c2 )   & 
    162160                     &      +         zswitch   *   rn_alb_smlt  
     
    211209            ! 
    212210            !                    ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions 
    213             WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice ) 
     211            WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 ) 
    214212               zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd 
    215213            END WHERE 
     
    256254         !                       !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time) 
    257255         zafrac_pnd = 0._wp 
    258          IF ( ld_pnd ) THEN   
    259             WHERE( ph_snw == 0._wp ) ;  zafrac_pnd = pafrac_pnd ;  END WHERE  ! Snow fully "shades" melt ponds 
     256         IF ( ld_pnd ) THEN 
     257            WHERE( ph_snw == 0._wp )   zafrac_pnd = pafrac_pnd   ! Snow fully "shades" melt ponds 
    260258         ENDIF          
    261259         ! 
    262260         !                       !--- Specific snow fraction (for now, prescribed) 
    263          WHERE     ( ph_snw > 0._wp     ) ;  zafrac_snw = 1. 
    264          ELSE WHERE                       ;  zafrac_snw = 0. 
     261         WHERE( ph_snw > 0._wp )   ;   zafrac_snw = 1. 
     262         ELSEWHERE                 ;   zafrac_snw = 0. 
    265263         END WHERE 
    266264         ! 
     
    273271         z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
    274272         z1_c2 = 1. / 0.05 
    275          WHERE     ( 1.5  < ph_ice                     )  ;  zalb_ice = zalb 
    276          ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_ice = zalb     + ( 0.18 - zalb     ) * z1_c1 *  & 
     273         WHERE    ( 1.5  < ph_ice                     )   ;   zalb_ice = zalb 
     274         ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_ice = zalb       + ( 0.18 - zalb       ) * z1_c1 *  & 
    277275            &                                                                       ( LOG(1.5) - LOG(ph_ice) ) 
    278          ELSE WHERE                                       ;  zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice 
     276         ELSEWHERE                                        ;   zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice 
    279277         END WHERE 
    280278         ! 
     
    282280         z1_c2 = 1. / 0.03 
    283281         ! 
    284          WHERE( pt_ice < rt0_snow ) ; zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 ); 
    285          ELSE WHERE                 ; zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 ); 
     282         WHERE( pt_ice < rt0_snow )   ;   zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 ) 
     283         ELSEWHERE                    ;   zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 ) 
    286284         END WHERE 
    287285         ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icecor.F90

    r8491 r8498  
    2626   USE lib_mpp        ! MPP library 
    2727   USE timing         ! Timing 
     28   USE iom            ! 
    2829 
    2930   IMPLICIT NONE 
     
    5354      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    5455      REAL(wp) ::   zsal, zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b, zzc 
     56      REAL(wp), DIMENSION(jpi,jpj) ::   zafx   ! concentration trends diag 
    5557      !!---------------------------------------------------------------------- 
    5658      IF( nn_timing == 1 )   CALL timing_start('icecor') 
     
    6567      ! 
    6668      ! 
    67       !                             !----------------------------------------------------- 
    68       IF( kn == 2 ) THEN            !  thickness of the smallest category above himin    ! 
     69      IF( kn == 2 ) THEN 
    6970         !                          !----------------------------------------------------- 
    70          ! 
    71          DO jj = 1, jpj  
    72             DO ji = 1, jpi 
    73 !!gm replace this 
    74                rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) )   !0 if no ice and 1 if yes 
    75                ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch 
    76 !!gm by more readable coding (not slower coding since already a IF in the loop): 
    77 !               IF( a_i(ji,jj,1) >= epsi20 )   ht_i(ji,jj,1) = v_i (ji,jj,1) / a_i(ji,jj,1) 
    78 !!gm 
    79                IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN 
    80                   a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin 
    81                ENDIF 
    82             END DO 
    83          END DO 
     71         !                          !  thickness of the smallest category above himin    ! 
     72         !                          !----------------------------------------------------- 
     73         WHERE( a_i(:,:,1) >= epsi20 )   ;   ht_i(:,:,1) = v_i (:,:,1) / a_i(:,:,1) 
     74         ELSEWHERE                       ;   ht_i(:,:,1) = 0._wp 
     75         END WHERE 
     76         WHERE( ht_i(:,:,1) < rn_himin )     a_i (:,:,1) = a_i (:,:,1) * ht_i(:,:,1) / rn_himin 
    8477         ! 
    8578      ENDIF 
    86  
    8779      !                             !----------------------------------------------------- 
    88       at_i(:,:) = a_i(:,:,1)        !  ice concentration should not exceed amax          ! 
    89       DO jl = 2, jpl                !----------------------------------------------------- 
    90          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     80      !                             !  ice concentration should not exceed amax          ! 
     81      !                             !----------------------------------------------------- 
     82      at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
     83      DO jl  = 1, jpl 
     84         WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:) 
    9185      END DO 
    92       ! 
    93 !!gm   Question   it seams to me that we have the following equality (dropping the "(ji,jj)": 
    94 !      ( 1. - ( 1. - rn_amax_2d / at_i ) ) =  ( 1. - ( at_i - rn_amax_2d ) / at_i ) 
    95 !                                          =  ( at_i - ( at_i - rn_amax_2d ) ) / at_i 
    96 !                                          =  ( + rn_amax_2d  ) / at_i 
    97 !                                          =  rn_amax_2d / at_i 
    98 !     No ?  if yes see "!!gm   better" juste below  
    99 !gm 
    100       DO jl  = 1, jpl 
    101          DO jj = 1, jpj 
    102             DO ji = 1, jpi 
    103                IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    104                   a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) ) 
    105 !!gm   better:    a_i(ji,jj,jl) = a_i(ji,jj,jl) * rn_amax_2d(ji,jj) / at_i(ji,jj)            
    106                ENDIF 
    107             END DO 
    108          END DO 
    109       END DO 
    110 !!gm  Other question:  why testing a_i(ji,jj,jl) > 0._wp ?   a_i is >=0, a multiplication by 0 does not change the results.... 
    111 !!gm  so at the end, the loop can be recoded without IF as: 
    112 !      WHERE( at_i(:,:) > rn_amax_2d(:,:) ) 
    113 !         DO jl  = 1, jpl 
    114 !            a_i(:,:,jl) = a_i(:,:,jl) * MAX( rn_amax_2d(:,:), at_i(:,:) ) / at_i(:,:) 
    115 !         END DO 
    116 !      END WHERE 
    117 !!gm  No? 
    11886     
    11987      !                             !----------------------------------------------------- 
    120       IF (  nn_icesal == 2  ) THEN  !  Ice salinity bounds                               ! 
     88      IF (  nn_icesal == 2  ) THEN  !  salinity stays in bounds [Simin,Simax]            ! 
    12189      !                             !----------------------------------------------------- 
    12290         zzc = rhoic * r1_rdtice 
    123          DO jl = 1, jpl                  ! salinity stays in bounds [Simin,Simax] 
     91         DO jl = 1, jpl 
    12492            DO jj = 1, jpj  
    12593               DO ji = 1, jpi 
    126                   IF( v_i(ji,jj,jl) > 0._wp ) THEN   ! clem: useless IF ??? 
    127                      zsal = smv_i(ji,jj,jl) 
    128                      smv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , smv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  ) 
    129                      ! associated salt flux 
    130                      sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * zzc 
    131                   ENDIF 
     94                  zsal = smv_i(ji,jj,jl) 
     95                  smv_i(ji,jj,jl) = MIN(  MAX( rn_simin*v_i(ji,jj,jl) , smv_i(ji,jj,jl) ) , rn_simax*v_i(ji,jj,jl)  ) 
     96                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * zzc   ! associated salt flux 
    13297               END DO 
    13398            END DO 
     
    157122         END DO 
    158123         CALL lbc_lnk_multi( u_ice, 'U', -1., v_ice, 'V', -1. )            ! lateral boundary conditions 
    159          ! 
    160 !!gm  I think masking here is unnecessary, u_ice already masked and we only introduce zeros in the field 
    161          u_ice(:,:) = u_ice(:,:) * umask(:,:,1)                            ! mask velocities 
    162          v_ice(:,:) = v_ice(:,:) * vmask(:,:,1) 
    163124      ENDIF 
    164125 
     
    172133      !                             !----------------------------------------------------- 
    173134      CASE( 1 )                        !--- dyn trend diagnostics 
    174          DO jl  = 1, jpl 
    175             afx_dyn(:,:) = afx_dyn(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 
    176          END DO 
    177135         ! 
    178136!!gm   here I think the number of ice cat is too small to use a SUM instruction... 
     
    188146            END DO 
    189147         END DO 
     148         !                       ! concentration tendency (dynamics) 
     149         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice  
     150         afx_tot(:,:) = zafx(:,:) 
     151         IF( iom_use('afxdyn') )   CALL iom_put( 'afxdyn' , zafx(:,:) ) 
    190152         ! 
    191153      CASE( 2 )                        !--- thermo trend diagnostics & ice aging 
    192154         ! 
    193          DO jl  = 1, jpl 
    194             oa_i(:,:,jl) = oa_i(:,:,jl) +   a_i(:,:,jl)                   * rdt_ice       ! ice natural aging incrementation 
    195             afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice     ! thermo tendancy 
    196          END DO 
    197          afx_tot(:,:) = afx_thd(:,:) + afx_dyn(:,:) 
     155         oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice   ! ice natural aging incrementation 
    198156         ! 
    199157!!gm   here I think the number of ice cat is too small to use a SUM instruction... 
     
    209167            END DO 
    210168         END DO 
     169         !                       ! concentration tendency (total + thermo) 
     170         zafx   (:,:) = SUM( a_i(:,:,:) - a_i_b(:,:,:), dim=3 ) * r1_rdtice 
     171         afx_tot(:,:) = afx_tot(:,:) + zafx(:,:) 
     172         IF( iom_use('afxthd') )   CALL iom_put( 'afxthd' , zafx(:,:) ) 
     173         IF( iom_use('afxtot') )   CALL iom_put( 'afxtot' , afx_tot(:,:) ) 
    211174         ! 
    212175      END SELECT 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90

    r8486 r8498  
    7171      REAL(wp), PARAMETER             :: zconv = 1.e-9 ! convert W to GW and kg to Mt 
    7272      !!---------------------------------------------------------------------- 
    73  
    74 !!gm  Note that glo_sum for a 2D or 3D array use a multiplication by tmask_i(ji,jj) 
    75 !!    so below  the  * tmask(:,:,1) is useless   ===>> I have removed them 
    76 !!    I also move the conversion factor from then glo_sum argument (become a single multiplication 
    77   
     73      ! 
    7874      IF( icount == 0 ) THEN 
    7975         !                          ! salt flux 
     
    271267               !WRITE(numout,*) ' sea-ice stress utau_ice  : ', utau_ice(ji,jj)  
    272268               !WRITE(numout,*) ' sea-ice stress vtau_ice  : ', vtau_ice(ji,jj) 
    273                !WRITE(numout,*) ' oceanic speed u          : ', u_oce(ji,jj) 
    274                !WRITE(numout,*) ' oceanic speed v          : ', v_oce(ji,jj) 
    275269               !WRITE(numout,*) ' sst                      : ', sst_m(ji,jj) 
    276270               !WRITE(numout,*) ' sss                      : ', sss_m(ji,jj) 
     
    556550               WRITE(numout,*) ' utau       : ', utau    (ji,jj)  
    557551               WRITE(numout,*) ' vtau       : ', vtau    (ji,jj) 
    558                WRITE(numout,*) ' oc. vel. u : ', u_oce   (ji,jj) 
    559                WRITE(numout,*) ' oc. vel. v : ', v_oce   (ji,jj) 
    560552            ENDIF 
    561553             
     
    669661      CALL prt_ctl(tab2d_1=utau       , clinfo1= ' utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ') 
    670662      CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ') 
    671       CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ') 
    672663       
    673664   END SUBROUTINE ice_prt3D 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90

    r8486 r8498  
    6565      IF( kt == nit000 .AND. lwp ) THEN 
    6666         WRITE(numout,*) 
    67          WRITE(numout,*)'ice_forcing_tau' 
     67         WRITE(numout,*)'ice_forcing_tau : Surface boundary condition for sea ice (momentum)' 
    6868         WRITE(numout,*)'~~~~~~~~~~~~~~~' 
    6969      ENDIF 
     
    124124      IF( kt == nit000 .AND. lwp ) THEN 
    125125         WRITE(numout,*) 
    126          WRITE(numout,*)'ice_forcing_flx' 
     126         WRITE(numout,*)'ice_forcing_flx : Surface boundary condition for sea ice (flux)' 
    127127         WRITE(numout,*)'~~~~~~~~~~~~~~~' 
    128128      ENDIF 
     
    132132 
    133133      ! albedo depends on cloud fraction because of non-linear spectral effects 
    134       DO jl = 1, jpl 
    135          DO jj = 2, jpjm1 
    136             DO ji = 2, jpim1 
    137134!!gm cldf_ice is a real, DOCTOR naming rule: start with cd means CHARACTER passed in argument ! 
    138                alb_ice(ji,jj,jl) = ( 1. - cldf_ice ) * zalb_cs(ji,jj,jl) + cldf_ice * zalb_os(ji,jj,jl) 
    139             END DO 
    140          END DO 
    141       END DO 
    142       CALL lbc_lnk( alb_ice, 'T', 1. ) 
     135      alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    143136       
    144137      SELECT CASE( ksbc )      !==  fluxes over sea ice  ==! 
     
    198191      INTEGER  ::   jl      ! dummy loop index 
    199192      ! 
    200       REAL(wp), DIMENSION(jpi,jpj) :: zalb_m    ! Mean albedo over all categories 
    201       REAL(wp), DIMENSION(jpi,jpj) :: ztem_m    ! Mean temperature over all categories 
    202       ! 
    203       REAL(wp), DIMENSION(jpi,jpj) :: z_qsr_m   ! Mean solar heat flux over all categories 
    204       REAL(wp), DIMENSION(jpi,jpj) :: z_qns_m   ! Mean non solar heat flux over all categories 
    205       REAL(wp), DIMENSION(jpi,jpj) :: z_evap_m  ! Mean sublimation over all categories 
    206       REAL(wp), DIMENSION(jpi,jpj) :: z_dqn_m   ! Mean d(qns)/dT over all categories 
    207       REAL(wp), DIMENSION(jpi,jpj) :: z_devap_m ! Mean d(evap)/dT over all categories 
     193      REAL(wp), DIMENSION(jpi,jpj) ::   z1_at_i   ! inverse of concentration 
     194      ! 
     195      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qsr_m   ! Mean solar heat flux over all categories 
     196      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_qns_m   ! Mean non solar heat flux over all categories 
     197      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_evap_m  ! Mean sublimation over all categories 
     198      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_dqn_m   ! Mean d(qns)/dT over all categories 
     199      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z_devap_m ! Mean d(evap)/dT over all categories 
     200      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zalb_m    ! Mean albedo over all categories 
     201      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztem_m    ! Mean temperature over all categories 
    208202      !!---------------------------------------------------------------------- 
    209203      ! 
    210204      IF( nn_timing == 1 )  CALL timing_start('ice_flx_dist') 
    211205      ! 
     206      WHERE ( at_i (:,:) > 0._wp )   ; z1_at_i(:,:) = 1._wp / at_i (:,:) 
     207      ELSEWHERE                      ; z1_at_i(:,:) = 0._wp 
     208      END WHERE 
     209       
    212210      SELECT CASE( k_limflx )       !==  averaged on all ice categories  ==! 
    213211      ! 
    214212      CASE( 0 , 1 ) 
    215          z_qns_m  (:,:) = fice_ice_ave( pqns_ice  (:,:,:) ) 
    216          z_qsr_m  (:,:) = fice_ice_ave( pqsr_ice  (:,:,:) ) 
    217          z_dqn_m  (:,:) = fice_ice_ave( pdqn_ice  (:,:,:) ) 
    218          z_evap_m (:,:) = fice_ice_ave( pevap_ice (:,:,:) ) 
    219          z_devap_m(:,:) = fice_ice_ave( pdevap_ice(:,:,:) ) 
    220 !!gm faster coding 
    221 !    REAL(wp), DIMENSION(jpi,jpj) ::   z1_at_i   !  
    222 ! ... 
    223 !      WHERE ( at_i (:,:) > 0._wp )   ; z1_at_i(:,:) = 1._wp / at_i (:,:) 
    224 !      ELSEWHERE                      ; z1_at_i(:,:) = 0._wp 
    225 !      END WHERE 
    226 !      z_qns_m  (:,:) = SUM( a_i(:,:,:) * pqns_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:) 
    227 !      z_qsr_m  (:,:) = SUM( a_i(:,:,:) * pqsr_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:) 
    228 !      z_dqn_m  (:,:) = SUM( a_i(:,:,:) * pdqn_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:) 
    229 !      z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 
    230 !      z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    231 !! and remove the 2 functions: fice_ice_ave and fice_cell_ave 
    232 !!gm 
     213         ! 
     214         ALLOCATE( z_qns_m(jpi,jpj), z_qsr_m(jpi,jpj), z_dqn_m(jpi,jpj), z_evap_m(jpi,jpj), z_devap_m(jpi,jpj) )   
     215         ! 
     216         z_qns_m  (:,:) = SUM( a_i(:,:,:) * pqns_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:) 
     217         z_qsr_m  (:,:) = SUM( a_i(:,:,:) * pqsr_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:) 
     218         z_dqn_m  (:,:) = SUM( a_i(:,:,:) * pdqn_ice  (:,:,:) , dim=3 ) * z1_at_i(:,:) 
     219         z_evap_m (:,:) = SUM( a_i(:,:,:) * pevap_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 
     220         z_devap_m(:,:) = SUM( a_i(:,:,:) * pdevap_ice(:,:,:) , dim=3 ) * z1_at_i(:,:) 
    233221         DO jl = 1, jpl 
    234             pdqn_ice  (:,:,jl) = z_dqn_m  (:,:) 
    235             pdevap_ice(:,:,jl) = z_devap_m(:,:) 
    236222            pqns_ice  (:,:,jl) = z_qns_m (:,:) 
    237223            pqsr_ice  (:,:,jl) = z_qsr_m (:,:) 
     224            pdqn_ice  (:,:,jl) = z_dqn_m  (:,:) 
    238225            pevap_ice (:,:,jl) = z_evap_m(:,:) 
     226            pdevap_ice(:,:,jl) = z_devap_m(:,:) 
    239227         END DO 
    240228         ! 
     229         DEALLOCATE( z_qns_m, z_qsr_m, z_dqn_m, z_evap_m, z_devap_m )   
     230         ! 
    241231      END SELECT 
    242232      ! 
    243233      SELECT CASE( k_limflx )       !==  redistribution on all ice categories  ==! 
     234      ! 
    244235      CASE( 1 , 2 ) 
    245236         ! 
    246          zalb_m(:,:) = fice_ice_ave( palb_ice(:,:,:) ) 
    247          ztem_m(:,:) = fice_ice_ave( ptn_ice (:,:,:) ) 
     237         ALLOCATE( zalb_m(jpi,jpj), ztem_m(jpi,jpj) )   
     238         ! 
     239         zalb_m(:,:) = SUM( a_i(:,:,:) * palb_ice(:,:,:) , dim=3 ) * z1_at_i(:,:) 
     240         ztem_m(:,:) = SUM( a_i(:,:,:) * ptn_ice (:,:,:) , dim=3 ) * z1_at_i(:,:) 
    248241         DO jl = 1, jpl 
    249242            pqns_ice (:,:,jl) = pqns_ice (:,:,jl) + pdqn_ice  (:,:,jl) * ( ptn_ice(:,:,jl) - ztem_m(:,:) ) 
     
    252245         END DO 
    253246         ! 
     247         DEALLOCATE( zalb_m, ztem_m )   
     248         ! 
    254249      END SELECT 
    255250      ! 
     
    257252      ! 
    258253   END SUBROUTINE ice_flx_dist 
    259  
    260 !!gm TO BE REMOVED ====>>>>> 
    261    FUNCTION fice_cell_ave ( ptab ) 
    262       !!-------------------------------------------------------------------------- 
    263       !! * Compute average over categories, for grid cell (ice covered and free ocean) 
    264       !!-------------------------------------------------------------------------- 
    265       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) :: ptab 
    266       REAL(wp), DIMENSION(jpi,jpj)                 :: fice_cell_ave 
    267       INTEGER :: jl 
    268       !!-------------------------------------------------------------------------- 
    269       fice_cell_ave(:,:) = a_i(:,:,1) * ptab (:,:,1) 
    270       DO jl = 2, jpl 
    271          fice_cell_ave(:,:) = fice_cell_ave(:,:) + a_i(:,:,jl) * ptab (:,:,jl) 
    272       END DO 
    273    END FUNCTION fice_cell_ave 
    274  
    275  
    276    FUNCTION fice_ice_ave ( ptab ) 
    277       !!-------------------------------------------------------------------------- 
    278       !! * Compute average over categories, for ice covered part of grid cell 
    279       !!-------------------------------------------------------------------------- 
    280       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in) ::   ptab   ! 
    281       REAL(wp), DIMENSION(jpi,jpj)                 :: fice_ice_ave 
    282       !!-------------------------------------------------------------------------- 
    283       WHERE ( at_i (:,:) > 0.0_wp ) ; fice_ice_ave (:,:) = fice_cell_ave( ptab (:,:,:) ) / at_i (:,:) 
    284       ELSEWHERE                     ; fice_ice_ave (:,:) = 0.0_wp 
    285       END WHERE 
    286    END FUNCTION fice_ice_ave 
    287  
    288 !!gm <<<<<<====  end of TO BE REMOVED  
    289254 
    290255#else 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceist.F90

    r8486 r8498  
    546546         WRITE(numout,*) 
    547547         WRITE(numout,*) 'ice_ist_init : ice parameters inititialisation ' 
    548          WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    549          WRITE(numout,*) '   initialization with ice (T) or not (F)       ln_limini     = ', ln_limini 
    550          WRITE(numout,*) '   ice initialization from a netcdf file      ln_limini_file  = ', ln_limini_file 
    551          WRITE(numout,*) '   threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst 
    552          WRITE(numout,*) '   initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n 
    553          WRITE(numout,*) '   initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s  
    554          WRITE(numout,*) '   initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n 
    555          WRITE(numout,*) '   initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s 
    556          WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n 
    557          WRITE(numout,*) '   initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s 
    558          WRITE(numout,*) '   initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n 
    559          WRITE(numout,*) '   initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s 
    560          WRITE(numout,*) '   initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n 
    561          WRITE(numout,*) '   initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s 
     548         WRITE(numout,*) '~~~~~~~~~~~~' 
     549         WRITE(numout,*) '   Namelist namiceini' 
     550         WRITE(numout,*) '      initialization with ice (T) or not (F)       ln_limini     = ', ln_limini 
     551         WRITE(numout,*) '      ice initialization from a netcdf file      ln_limini_file  = ', ln_limini_file 
     552         WRITE(numout,*) '      threshold water temp. for initial sea-ice    rn_thres_sst  = ', rn_thres_sst 
     553         WRITE(numout,*) '      initial snow thickness in the north          rn_hts_ini_n  = ', rn_hts_ini_n 
     554         WRITE(numout,*) '      initial snow thickness in the south          rn_hts_ini_s  = ', rn_hts_ini_s  
     555         WRITE(numout,*) '      initial ice thickness  in the north          rn_hti_ini_n  = ', rn_hti_ini_n 
     556         WRITE(numout,*) '      initial ice thickness  in the south          rn_hti_ini_s  = ', rn_hti_ini_s 
     557         WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_n  = ', rn_ati_ini_n 
     558         WRITE(numout,*) '      initial ice concentr.  in the north          rn_ati_ini_s  = ', rn_ati_ini_s 
     559         WRITE(numout,*) '      initial  ice salinity  in the north          rn_smi_ini_n  = ', rn_smi_ini_n 
     560         WRITE(numout,*) '      initial  ice salinity  in the south          rn_smi_ini_s  = ', rn_smi_ini_s 
     561         WRITE(numout,*) '      initial  ice/snw temp  in the north          rn_tmi_ini_n  = ', rn_tmi_ini_n 
     562         WRITE(numout,*) '      initial  ice/snw temp  in the south          rn_tmi_ini_s  = ', rn_tmi_ini_s 
    562563      ENDIF 
    563564 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90

    r8486 r8498  
    426426            ENDIF 
    427427            ! 
    428 !!gm use of rswitch not faster as there is already IF in the DO-loop ==>>> use IF which is more readable 
    429 !            rswitch    = MAX( 0._wp , SIGN( 1._wp , v_i_2d(ji,jl1) - epsi10 ) ) 
    430 !            zworkv(ji) = pdvice(ji,jl) / MAX( v_i_2d(ji,jl1), epsi10 ) * rswitch 
    431 !            ! 
    432 !            rswitch    = MAX( 0._wp , SIGN( 1._wp , a_i_2d(ji,jl1) - epsi10 ) ) 
    433 !            zworka(ji) = pdaice(ji,jl) / MAX( a_i_2d(ji,jl1), epsi10 ) * rswitch         
    434  
    435428            IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1) 
    436429            ELSE                                  ;   zworkv(ji) = 0._wp 
     
    439432            ELSE                                  ;   zworka(ji) = 0._wp 
    440433            ENDIF 
    441 !!gm end 
    442434            ! 
    443435            a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas 
  • 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      ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90

    r8486 r8498  
    3131   USE prtctl         ! Print control 
    3232   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     33   USE iom            ! 
    3334#if defined key_agrif 
    3435   USE agrif_lim3_interp 
     
    111112      INTEGER ::   ji, jj       ! dummy loop indices 
    112113      INTEGER ::   jter         ! local integers 
    113       CHARACTER (len=50) ::   charout 
    114114 
    115115      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio 
     
    121121      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! temporary scalars 
    122122 
    123       REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress 
    124123      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity 
    125124      REAL(wp) ::   zintb, zintn                                             ! dummy argument 
    126125      REAL(wp) ::   zfac_x, zfac_y 
     126      REAL(wp) ::   zshear, zdum1, zdum2 
    127127       
    128128      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors 
     
    152152      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
    153153      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity 
     154      !! --- diags 
     155      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zswi, zsig1, zsig2, zsig3 
     156      !! --- SIMIP diags 
     157      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice    
     158      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig2      ! Maximum shear stress in sea ice 
     159      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dx   ! X-direction sea-surface tilt term (N/m2) 
     160      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_dssh_dy   ! X-direction sea-surface tilt term (N/m2) 
     161      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstrx   ! X-direction coriolis stress (N/m2) 
     162      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_corstry   ! Y-direction coriolis stress (N/m2) 
     163      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstrx   ! X-direction internal stress (N/m2) 
     164      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_intstry   ! Y-direction internal stress (N/m2) 
     165      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_utau_oi   ! X-direction ocean-ice stress 
     166      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_vtau_oi   ! Y-direction ocean-ice stress   
     167      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s) 
     168      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s) 
     169      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s) 
     170      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s) 
     171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s) 
     172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)       
    154173      !!------------------------------------------------------------------- 
    155174 
     
    671690 
    672691      !------------------------------------------------------------------------------! 
    673       ! 5) SIMIP diagnostics 
    674       !------------------------------------------------------------------------------! 
    675  
    676 !!gm  encapsulate with a flag (iom_use of the variable or better a flag defined one for all testing if one of the 
    677 !!    diag is output.  NB the diag_...  are should only be ALLOCATED if the flag is true ! 
    678  
    679       DO jj = 2, jpjm1 
    680          DO ji = 2, jpim1 
    681              rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
    682  
    683              ! Stress tensor invariants (normal and shear stress N/m) 
    684              diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress 
    685              diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress 
    686  
    687              ! Stress terms of the momentum equation (N/m2) 
    688              diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term 
    689              diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 
    690  
    691              diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term 
    692              diag_corstry(ji,jj) = zCory(ji,jj) * rswitch 
    693  
    694              diag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term 
    695              diag_intstry(ji,jj) = zfV(ji,jj)   * rswitch 
    696             
    697              diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress 
    698              diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 
    699  
    700              ! 2D ice mass, snow mass, area transport arrays (X, Y) 
    701              zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 
    702              zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
    703  
    704              diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
    705              diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
    706  
    707              diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
    708              diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
    709  
    710              diag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component 
    711              diag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   '' 
    712  
    713          END DO 
    714       END DO 
    715  
    716       CALL lbc_lnk_multi(   diag_sig1   , 'T',  1., diag_sig2   , 'T',  1.,   & 
    717          &                  diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   & 
    718          &                  diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   &  
    719          &                  diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    ) 
    720  
    721       CALL lbc_lnk_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    ) 
    722  
    723       CALL lbc_lnk_multi(   diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1.,   & 
    724          &                  diag_xatrp    , 'U', -1., diag_ymtrp_ice, 'V', -1.,   & 
    725          &                  diag_ymtrp_snw, 'V', -1., diag_yatrp    , 'V', -1.    ) 
    726  
    727       ! 
    728       !------------------------------------------------------------------------------! 
    729       ! 6) Control prints of residual and charge ellipse 
    730       !------------------------------------------------------------------------------! 
    731       ! 
    732       ! print the residual for convergence 
    733       IF(ln_ctl) THEN 
    734          WRITE(charout,FMT="('ice_rhg_evp  : res =',D23.16, ' iter =',I4)") zresm, jter 
    735          CALL prt_ctl_info(charout) 
    736          CALL prt_ctl(tab2d_1=u_ice, clinfo1=' ice_rhg_evp  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :') 
     692      ! 5) diagnostics 
     693      !------------------------------------------------------------------------------! 
     694      ! --- ellipse --- ! 
     695      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN 
    737696         ! 
    738          ! print charge ellipse 
    739          ! This can be desactivated once the user is sure that the stress state 
    740       ! lie on the charge ellipse. See Bouillon et al. (2008) for more details 
    741          IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN 
    742             WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. ' 
    743             CALL prt_ctl_info(charout) 
    744             DO jj = 2, jpjm1 
    745                DO ji = 2, jpim1 
    746                   IF (strength(ji,jj) > 1.0) THEN 
    747                      zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )  
    748                      zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
    749                      WRITE(charout,FMT="('ice_rhg_evp  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 
    750                      CALL prt_ctl_info(charout) 
    751                   ENDIF 
    752                END DO 
    753             END DO 
    754             WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 2000, kt_ice, 0, 0, ' ch. ell. ' 
    755             CALL prt_ctl_info(charout) 
    756          ENDIF 
     697         ALLOCATE( zswi(jpi,jpj) , zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) ) 
     698         ! 
     699         DO jj = 1, jpj 
     700            DO ji = 1, jpi 
     701               zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
     702            END DO 
     703         END DO 
     704          
     705         DO jj = 2, jpjm1 
     706            DO ji = 2, jpim1 
     707               zdum1 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
     708                  &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  & 
     709                  &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 
     710 
     711               zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress   
     712 
     713               zdum2 = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
     714 
     715!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
     716!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
     717!!               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
     718!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
     719               zsig1(ji,jj) = 0.5_wp * zdum2 * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
     720               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                    ! shear stress 
     721               zsig3(ji,jj) = zdum2**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
     722            END DO 
     723         END DO 
     724         CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
     725         ! 
     726         IF( iom_use('isig1') )   CALL iom_put( "isig1" , zsig1 ) 
     727         IF( iom_use('isig2') )   CALL iom_put( "isig2" , zsig2 ) 
     728         IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 ) 
     729         ! 
     730         DEALLOCATE( zswi , zsig1 , zsig2 , zsig3 ) 
     731      ENDIF 
     732       
     733      ! --- SIMIP --- ! 
     734      IF ( iom_use( 'normstr'  ) .OR. iom_use( 'sheastr'  ) .OR. iom_use( 'dssh_dx'  ) .OR. iom_use( 'dssh_dy'  ) .OR. & 
     735         & iom_use( 'corstrx'  ) .OR. iom_use( 'corstry'  ) .OR. iom_use( 'intstrx'  ) .OR. iom_use( 'intstry'  ) .OR. & 
     736         & iom_use( 'utau_oi'  ) .OR. iom_use( 'vtau_oi'  ) .OR. iom_use( 'xmtrpice' ) .OR. iom_use( 'ymtrpice' ) .OR. & 
     737         & iom_use( 'xmtrpsnw' ) .OR. iom_use( 'ymtrpsnw' ) .OR. iom_use( 'xatrp'    ) .OR. iom_use( 'yatrp'    ) ) THEN 
     738 
     739         ALLOCATE( zdiag_sig1     (jpi,jpj) , zdiag_sig2     (jpi,jpj) , zdiag_dssh_dx  (jpi,jpj) , zdiag_dssh_dy  (jpi,jpj) ,  & 
     740            &      zdiag_corstrx  (jpi,jpj) , zdiag_corstry  (jpi,jpj) , zdiag_intstrx  (jpi,jpj) , zdiag_intstry  (jpi,jpj) ,  & 
     741            &      zdiag_utau_oi  (jpi,jpj) , zdiag_vtau_oi  (jpi,jpj) , zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) ,  & 
     742            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp    (jpi,jpj) , zdiag_yatrp    (jpi,jpj) ) 
     743          
     744         DO jj = 2, jpjm1 
     745            DO ji = 2, jpim1 
     746               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
     747                
     748               ! Stress tensor invariants (normal and shear stress N/m) 
     749               zdiag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress 
     750               zdiag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress 
     751                
     752               ! Stress terms of the momentum equation (N/m2) 
     753               zdiag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term 
     754               zdiag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch 
     755                
     756               zdiag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term 
     757               zdiag_corstry(ji,jj) = zCory(ji,jj) * rswitch 
     758                
     759               zdiag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term 
     760               zdiag_intstry(ji,jj) = zfV(ji,jj)   * rswitch 
     761                
     762               zdiag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress 
     763               zdiag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch 
     764                
     765               ! 2D ice mass, snow mass, area transport arrays (X, Y) 
     766               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch 
     767               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch 
     768                
     769               zdiag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component 
     770               zdiag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   '' 
     771                
     772               zdiag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component 
     773               zdiag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   '' 
     774                
     775               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component 
     776               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   '' 
     777                
     778            END DO 
     779         END DO 
     780          
     781         CALL lbc_lnk_multi(   zdiag_sig1   , 'T',  1., zdiag_sig2   , 'T',  1.,   & 
     782            &                  zdiag_dssh_dx, 'U', -1., zdiag_dssh_dy, 'V', -1.,   & 
     783            &                  zdiag_corstrx, 'U', -1., zdiag_corstry, 'V', -1.,   &  
     784            &                  zdiag_intstrx, 'U', -1., zdiag_intstry, 'V', -1.    ) 
     785          
     786         CALL lbc_lnk_multi(   zdiag_utau_oi, 'U', -1., zdiag_vtau_oi, 'V', -1.    ) 
     787          
     788         CALL lbc_lnk_multi(   zdiag_xmtrp_ice, 'U', -1., zdiag_xmtrp_snw, 'U', -1.,   & 
     789            &                  zdiag_xatrp    , 'U', -1., zdiag_ymtrp_ice, 'V', -1.,   & 
     790            &                  zdiag_ymtrp_snw, 'V', -1., zdiag_yatrp    , 'V', -1.    ) 
     791          
     792         IF ( iom_use( 'normstr'  ) )   CALL iom_put( 'normstr'  ,  zdiag_sig1(:,:)      )   ! Normal stress 
     793         IF ( iom_use( 'sheastr'  ) )   CALL iom_put( 'sheastr'  ,  zdiag_sig2(:,:)      )   ! Shear stress 
     794         IF ( iom_use( 'dssh_dx'  ) )   CALL iom_put( 'dssh_dx'  ,  zdiag_dssh_dx(:,:)   )   ! Sea-surface tilt term in force balance (x) 
     795         IF ( iom_use( 'dssh_dy'  ) )   CALL iom_put( 'dssh_dy'  ,  zdiag_dssh_dy(:,:)   )   ! Sea-surface tilt term in force balance (y) 
     796         IF ( iom_use( 'corstrx'  ) )   CALL iom_put( 'corstrx'  ,  zdiag_corstrx(:,:)   )   ! Coriolis force term in force balance (x) 
     797         IF ( iom_use( 'corstry'  ) )   CALL iom_put( 'corstry'  ,  zdiag_corstry(:,:)   )   ! Coriolis force term in force balance (y) 
     798         IF ( iom_use( 'intstrx'  ) )   CALL iom_put( 'intstrx'  ,  zdiag_intstrx(:,:)   )   ! Internal force term in force balance (x) 
     799         IF ( iom_use( 'intstry'  ) )   CALL iom_put( 'intstry'  ,  zdiag_intstry(:,:)   )   ! Internal force term in force balance (y) 
     800         IF ( iom_use( 'utau_oi'  ) )   CALL iom_put( 'utau_oi'  ,  zdiag_utau_oi(:,:)   )   ! Ocean stress term in force balance (x) 
     801         IF ( iom_use( 'vtau_oi'  ) )   CALL iom_put( 'vtau_oi'  ,  zdiag_vtau_oi(:,:)   )   ! Ocean stress term in force balance (y) 
     802         IF ( iom_use( 'xmtrpice' ) )   CALL iom_put( 'xmtrpice' ,  zdiag_xmtrp_ice(:,:) )   ! X-component of sea-ice mass transport (kg/s) 
     803         IF ( iom_use( 'ymtrpice' ) )   CALL iom_put( 'ymtrpice' ,  zdiag_ymtrp_ice(:,:) )   ! Y-component of sea-ice mass transport  
     804         IF ( iom_use( 'xmtrpsnw' ) )   CALL iom_put( 'xmtrpsnw' ,  zdiag_xmtrp_snw(:,:) )   ! X-component of snow mass transport (kg/s) 
     805         IF ( iom_use( 'ymtrpsnw' ) )   CALL iom_put( 'ymtrpsnw' ,  zdiag_ymtrp_snw(:,:) )   ! Y-component of snow mass transport 
     806         IF ( iom_use( 'xatrp'    ) )   CALL iom_put( 'xatrp'    ,  zdiag_xatrp(:,:)     )   ! X-component of ice area transport 
     807         IF ( iom_use( 'yatrp'    ) )   CALL iom_put( 'yatrp'    ,  zdiag_yatrp(:,:)     )   ! Y-component of ice area transport 
     808 
     809         DEALLOCATE( zdiag_sig1      , zdiag_sig2      , zdiag_dssh_dx   , zdiag_dssh_dy   ,  & 
     810            &        zdiag_corstrx   , zdiag_corstry   , zdiag_intstrx   , zdiag_intstry   ,  & 
     811            &        zdiag_utau_oi   , zdiag_vtau_oi   , zdiag_xmtrp_ice , zdiag_ymtrp_ice ,  & 
     812            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp     , zdiag_yatrp     ) 
     813 
    757814      ENDIF 
    758815      ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

    r8486 r8498  
    122122         !                                   !-----------------------! 
    123123         ! 
    124          kt_ice = kt       ! Ice model time step 
     124         kt_ice = kt                              ! -- Ice model time step 
    125125         ! 
    126126# if defined key_agrif 
    127127         IF( .NOT. Agrif_Root() )  lim_nbstep = MOD( lim_nbstep, Agrif_irhot() * Agrif_Parent(nn_fsbc) / nn_fsbc ) + 1 
    128128# endif 
    129  
    130          !                 ! mean surface ocean current at ice velocity point 
    131          u_oce(:,:) = ssu_m(:,:) * umask(:,:,1) 
    132          v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1) 
    133 !!gm the umask, vmask above is useless as ssu_m, ssv_m are build from masked un,vn (used t be here for B-grid sea-ice) 
    134  
    135          !                 ! masked sea surface freezing temperature [Kelvin] (set to rt0 over land) 
    136          CALL eos_fzp( sss_m(:,:) , t_bo(:,:) ) 
     129         u_oce(:,:) = ssu_m(:,:)                  ! -- mean surface ocean current 
     130         v_oce(:,:) = ssv_m(:,:) 
     131         ! 
     132         CALL eos_fzp( sss_m(:,:) , t_bo(:,:) )   ! -- freezing temperature [Kelvin] (set to rt0 over land) 
    137133         t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) + rt0 * ( 1._wp - tmask(:,:,1) ) 
    138134         ! 
    139          CALL ice_bef      ! Store previous ice values 
     135         CALL store_fields                        ! -- Store now ice values 
    140136         ! 
    141137         !------------------------------------------------! 
     
    179175         CALL ice_var_glo2eqv            ! ht_i and ht_s for ice albedo calculation 
    180176         CALL ice_var_agg(1)             ! at_i for coupling  
    181          CALL ice_bef                    ! Store previous ice values 
     177         CALL store_fields               ! Store now ice values 
    182178 
    183179         !------------------------------------------------------! 
     
    223219!!         IF( .NOT. Agrif_Root() )   CALL Agrif_ParentGrid_To_ChildGrid() 
    224220!!# endif 
    225          IF( ln_limdiahsb )           CALL ice_dia( kt )        ! -- Diagnostics and outputs  
    226          ! 
    227                                       CALL ice_wri( 1 )         ! -- Ice outputs  
     221         IF( ln_limdiahsb )         CALL ice_dia( kt )        ! -- Diagnostics and outputs  
     222         ! 
     223                                    CALL ice_wri( 1 )         ! -- Ice outputs  
    228224         ! 
    229225         IF( kt == nit000 .AND. ln_rstart )   &                !!gm  I don't understand the ln_rstart, if needed, add a comment, please 
    230             &                         CALL iom_close( numrir )  ! close input ice restart file 
    231          ! 
    232          IF( lrst_ice )               CALL ice_rst_write( kt )  ! -- Ice restart file  
    233          ! 
    234          IF( ln_limctl )              CALL ice_ctl( kt )        ! alerts in case of model crash 
     226            &                       CALL iom_close( numrir )  ! close input ice restart file 
     227         ! 
     228         IF( lrst_ice )             CALL ice_rst_write( kt )  ! -- Ice restart file  
     229         ! 
     230         IF( ln_limctl )            CALL ice_ctl( kt )        ! alerts in case of model crash 
    235231         ! 
    236232      ENDIF   ! End sea-ice time step only 
     
    241237      ! Update surface ocean stresses (only in ice-dynamic case) otherwise the atm.-ocean stresses are used everywhere 
    242238      !    using before instantaneous surf. currents 
    243       IF( ln_limdyn )                 CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) ) 
     239      IF( ln_limdyn )               CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) ) 
    244240!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    245241      ! 
     
    476472 
    477473 
    478    SUBROUTINE ice_bef 
    479       !!---------------------------------------------------------------------- 
    480       !!                  ***  ROUTINE ice_bef  *** 
     474   SUBROUTINE store_fields 
     475      !!---------------------------------------------------------------------- 
     476      !!                  ***  ROUTINE store_fields  *** 
    481477      !! 
    482478      !! ** purpose :  store ice variables at "before" time step 
     
    525521      CALL lbc_lnk_multi( at_i_b, 'T', 1., u_ice_b , 'U', -1., v_ice_b , 'V', -1. ) 
    526522      ! 
    527    END SUBROUTINE ice_bef 
     523   END SUBROUTINE store_fields 
    528524 
    529525 
     
    570566            ! 
    571567            afx_tot(ji,jj) = 0._wp   ; 
    572             afx_dyn(ji,jj) = 0._wp   ;   afx_thd(ji,jj) = 0._wp 
    573568            ! 
    574569            diag_heat(ji,jj) = 0._wp ;   diag_smvi(ji,jj) = 0._wp 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

    r8486 r8498  
    193193      !     Second step in icethd_dh      :  heat remaining if total melt (zq_rema)  
    194194      !     Third  step in iceupdate.F90  :  heat from ice-ocean mass exchange (zf_mass) + solar 
    195       DO jj = 1, jpj 
    196          DO ji = 1, jpi 
    197             hfx_out(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj)  &  ! Non solar heat flux received by the ocean                
    198                &             - qlead(ji,jj) * r1_rdtice                                    &  ! heat flux taken from the ocean where there is open water ice formation 
    199                &             - at_i(ji,jj) * fhtur(ji,jj)                                  &  ! heat flux taken by turbulence 
    200                &             - at_i(ji,jj) *  fhld(ji,jj)                                     ! heat flux taken during bottom growth/melt  
    201                                                                                               !    (fhld should be 0 while bott growth) 
    202          END DO 
    203       END DO 
    204  
     195      hfx_out(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + qemp_oce(:,:)  &  ! Non solar heat flux received by the ocean                
     196         &           - qlead(:,:) * r1_rdtice                                &  ! heat flux taken from the ocean where there is open water ice formation 
     197         &           - at_i (:,:) * fhtur(:,:)                               &  ! heat flux taken by turbulence 
     198         &           - at_i (:,:) *  fhld(:,:)                                  ! heat flux taken during bottom growth/melt  
     199                                                                                !    (fhld should be 0 while bott growth) 
    205200      !-------------------------------------------------------------------------------------------! 
    206201      ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories 
     
    288283      !!------------------------------------------------------------------- 
    289284      INTEGER  ::   ji, jk   ! dummy loop indices 
    290       REAL(wp) ::   ztmelts, zaaa, zbbb, zccc, zdiscrim  ! local scalar  
     285      REAL(wp) ::   ztmelts, z1_2cp, zbbb, zccc  ! local scalar  
    291286      !!------------------------------------------------------------------- 
    292287      ! Recover ice temperature 
     288      z1_2cp = 1._wp / ( 2._wp * cpic ) 
    293289      DO jk = 1, nlay_i 
    294290         DO ji = 1, nidx 
    295             ztmelts       =  -tmut * s_i_1d(ji,jk) + rt0 
     291            ztmelts       = -tmut * s_i_1d(ji,jk) 
    296292            ! Conversion q(S,T) -> T (second order equation) 
    297             zaaa          =  cpic 
    298             zbbb          =  ( rcp - cpic ) * ( ztmelts - rt0 ) + e_i_1d(ji,jk) * r1_rhoic - lfus 
    299             zccc          =  lfus * ( ztmelts - rt0 ) 
    300             zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 
    301             t_i_1d(ji,jk) =  rt0 - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 
     293            zbbb          = ( rcp - cpic ) * ztmelts + e_i_1d(ji,jk) * r1_rhoic - lfus 
     294            zccc          = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts, 0._wp ) ) 
     295            t_i_1d(ji,jk) = rt0 - ( zbbb + zccc ) * z1_2cp 
    302296             
    303297            ! mask temperature 
    304             rswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )  
    305             t_i_1d(ji,jk) =  rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0 
     298            rswitch       = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) )  
     299            t_i_1d(ji,jk) = rswitch * t_i_1d(ji,jk) + ( 1._wp - rswitch ) * rt0 
    306300         END DO  
    307301      END DO  
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_da.F90

    r8486 r8498  
    3232CONTAINS 
    3333 
    34 !!gm  even comment line of more than 130 character may cause compilation errors 
    35 !!gm         ===>>> reformat the text 
    3634   SUBROUTINE ice_thd_da 
    3735      !!------------------------------------------------------------------- 
     
    4644      !!                   W = m1 * (Tw -Tf)**m2                    --- originally from Josberger 1979 --- 
    4745      !!                      (Tw - Tf) = elevation of water temp above freezing 
    48       !!                      m1 and m2 = (1.6e-6 , 1.36) best fit from field experiment near the coast of Prince Patrick Island (Perovich 1983) => static ice 
    49       !!                      m1 and m2 = (3.0e-6 , 1.36) best fit from MIZEX 84 experiment (Maykut and Perovich 1987) => moving ice 
     46      !!                      m1 and m2 = (1.6e-6 , 1.36) best fit from field experiment near the coast of Prince Patrick Island 
     47      !!                                                                                           (Perovich 1983) => static ice 
     48      !!                      m1 and m2 = (3.0e-6 , 1.36) best fit from MIZEX 84 experiment 
     49      !!                                                                                (Maykut and Perovich 1987) => moving ice 
    5050      !! 
    5151      !!                   P = N * pi * D                           --- from Rothrock and Thorndike 1984 --- 
     
    5959      !!                      Astar = 1 / ( 1 - (Dmin/Dmax)**(1/beta) ) 
    6060      !!                      Dmin = minimum floe diameter (recommended to be 8m +- 20%) 
    61       !!                      Dmax = maximum floe diameter (recommended to be 300m, but it does not impact melting much except for Dmax<100m) 
     61      !!                      Dmax = maximum floe diameter (recommended to be 300m, 
     62      !!                                                    but it does not impact melting much except for Dmax<100m) 
    6263      !!                      beta = 1.0 +-20% (recommended value) 
    6364      !!                           = 0.3 best fit for western Fram Strait and Antarctica 
     
    113114         ! --- Calculate reduction of total sea ice concentration --- ! 
    114115         zdfloe = rn_dmin * ( zastar / ( zastar - at_i_1d(ji) ) )**rn_beta         ! Mean floe caliper diameter [m] 
    115          zperi  = at_i_1d(ji) * rpi / ( zcs * zdfloe )                             ! Mean perimeter of the floe = N*pi*D = (A/cs*D^2)*pi*D [m.m-2] 
     116         ! 
     117         zperi  = at_i_1d(ji) * rpi / ( zcs * zdfloe )                             ! Mean perimeter of the floe [m.m-2] 
     118         !                                                                         !    = N*pi*D = (A/cs*D^2)*pi*D 
    116119         zwlat  = zm1 * ( MAX( 0._wp, sst_1d(ji) - ( t_bo_1d(ji) - rt0 ) ) )**zm2  ! Melt speed rate [m/s] 
    117           
     120         ! 
    118121         zda_tot(ji) = MIN( zwlat * zperi * rdt_ice, at_i_1d(ji) )                 ! sea ice concentration decrease (>0) 
    119122       
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

    r8486 r8498  
    102102      REAL(wp), DIMENSION(jpij) ::   zsnw        ! distribution of snow after wind blowing 
    103103 
    104       REAL(wp) :: zswitch_sal 
     104      REAL(wp) ::   zswitch_sal 
    105105 
    106106      INTEGER  ::   num_iter_max      ! Heat conservation  
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90

    r8486 r8498  
    5353      !!--------------------------------------------------------------------- 
    5454 
    55       !--------------------------------------------------------------------| 
    56       ! 1) salinity constant in time                                       | 
    57       !--------------------------------------------------------------------| 
    58       ! do nothing 
    59  
    60       !----------------------------------------------------------------------| 
    61       !  2) salinity varying in time                                         | 
    62       !----------------------------------------------------------------------| 
    63       IF(  nn_icesal == 2  ) THEN 
    64  
     55      SELECT CASE ( nn_icesal ) 
     56      ! 
     57      !              !---------------------------------------------! 
     58      CASE( 2 )      !  time varying salinity with linear profile  ! 
     59      !              !---------------------------------------------! 
    6560         DO ji = 1, nidx 
    6661 
     
    9691         CALL ice_var_salprof1d 
    9792         ! 
    98       ENDIF  
    99  
    100       !------------------------------------------------------------------------------| 
    101       !  3) vertical profile of salinity, constant in time                           | 
    102       !------------------------------------------------------------------------------| 
    103       IF(  nn_icesal == 3  )   CALL ice_var_salprof1d 
     93      !                 !---------------------------------------------! 
     94      CASE( 3 )         ! constant salinity with a fix profile        ! (Schwarzacher (1959) multiyear salinity profile(mean = 2.30) 
     95      !                 !---------------------------------------------! 
     96         CALL ice_var_salprof1d 
    10497      ! 
     98   END SELECT 
     99   ! 
    105100   END SUBROUTINE ice_thd_sal 
    106101 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90

    r8486 r8498  
    8686      !!   
    8787      !! ** Purpose :   Update the surface ocean boundary condition for heat  
    88       !!              salt and mass over areas where sea-ice is non-zero 
     88      !!                salt and mass over areas where sea-ice is non-zero 
    8989      !!          
    9090      !! ** Action  : - computes the heat and freshwater/salt fluxes 
    91       !!              at the ice-ocean interface. 
     91      !!                at the ice-ocean interface. 
    9292      !!              - Update the ocean sbc 
    9393      !!      
     
    133133         DO ji = 1, jpi 
    134134 
    135             !------------------------------------------! 
    136             !      heat flux at the ocean surface      ! 
    137             !------------------------------------------! 
    138135            ! Solar heat flux reaching the ocean = zqsr (W.m-2)  
    139136            !--------------------------------------------------- 
    140             zqsr = qsr_tot(ji,jj) 
    141             DO jl = 1, jpl 
    142                zqsr = zqsr - a_i_b(ji,jj,jl) * (  qsr_ice(ji,jj,jl) - ftr_ice(ji,jj,jl) )  
    143             END DO 
    144 !!gm  why not like almost everywhere else : 
    145 !!gm        zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * (  qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) ) 
     137            zqsr = qsr_tot(ji,jj) - SUM( a_i_b(ji,jj,:) * ( qsr_ice(ji,jj,:) - ftr_ice(ji,jj,:) ) ) 
    146138 
    147139            ! Total heat flux reaching the ocean = hfx_out (W.m-2)  
     
    204196      END DO 
    205197 
    206       !-----------------------------------------------! 
    207       !   Storing the transmitted variables           ! 
    208       !-----------------------------------------------! 
     198      ! Storing the transmitted variables 
     199      !---------------------------------- 
    209200      fr_i  (:,:)   = at_i(:,:)             ! Sea-ice fraction             
    210201      tn_ice(:,:,:) = t_su(:,:,:)           ! Ice surface temperature                       
    211202 
    212       !------------------------------------------------------------------------! 
    213       !    Snow/ice albedo (only if sent to coupler, useless in forced mode)   ! 
    214       !------------------------------------------------------------------------! 
     203      ! Snow/ice albedo (only if sent to coupler, useless in forced mode) 
     204      !------------------------------------------------------------------ 
    215205      CALL ice_alb( t_su, ht_i, ht_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    216206      ! 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

    r8496 r8498  
    157157      !!------------------------------------------------------------------ 
    158158      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    159       REAL(wp) ::   ze_i, z1_CpR, zdiscrim, zaaa, z1_2aaa             ! local scalars 
    160       REAL(wp) ::   ze_s, zL_Cp , ztmelts , zbbb, zccc                !   -      - 
    161       REAL(wp) ::   zhmax, z1_zhmax, zsm_i, zcpMcpic, zt_i   !   -      - 
    162       REAL(wp) ::   zlay_i, zlay_s   !   -      - 
     159      REAL(wp) ::   ze_i, z1_cp, z1_2cp             ! local scalars 
     160      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      - 
     161      REAL(wp) ::   zhmax, z1_zhmax, zsm_i          !   -      - 
     162      REAL(wp) ::   zlay_i, zlay_s                  !   -      - 
    163163      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i 
    164164      !!------------------------------------------------------------------ 
     
    181181      ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness 
    182182 
    183      zhmax    =          hi_max(jpl) 
     183      zhmax    =          hi_max(jpl) 
    184184      z1_zhmax =  1._wp / hi_max(jpl)                
    185185      WHERE( ht_i(:,:,jpl) > zhmax )               !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area 
     
    205205      !------------------- 
    206206      zlay_i   = REAL( nlay_i , wp )    ! number of layers 
    207       zaaa     = cpic                   ! Conversion q(S,T) -> T (second order equation) 
    208       z1_2aaa  = 1._wp / ( 2._wp * zaaa ) 
    209       zcpMcpic = rcp - cpic 
     207      z1_2cp  = 1._wp / ( 2._wp * cpic ) 
    210208      DO jl = 1, jpl 
    211209         DO jk = 1, nlay_i 
     
    214212                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area  
    215213                     ! 
    216                      ze_i    =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3] 
    217                      ztmelts = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C] 
    218                      ! 
    219                      zbbb     =  zcpMcpic * ztmelts + ze_i * r1_rhoic - lfus 
    220                      zccc     =  lfus * ztmelts 
    221                      zdiscrim =  SQRT(  MAX( zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp)  ) 
    222                      zt_i     = - ( zbbb + zdiscrim ) * z1_2aaa                          ! [C] 
    223                      t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( zt_i , ztmelts )  ) + rt0   ! [K] with bounds: -100 < zt_i < ztmelts 
     214                     ze_i             =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3] 
     215                     ztmelts          = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C] 
     216                     ! Conversion q(S,T) -> T (second order equation) 
     217                     zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus 
     218                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) ) 
     219                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * z1_2cp , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts 
    224220                     ! 
    225221                  ELSE                                !--- no ice 
    226                      t_i(ji,jj,jk,jl) =  rt0 - 100._wp                                   ! impose 173.15 K (i.e. -100 C) 
     222                     t_i(ji,jj,jk,jl) = rt0 - 100._wp                                   ! impose 173.15 K (i.e. -100 C) 
    227223!!clem: I think we should impose rt0 instead 
    228224                  ENDIF 
     
    235231      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.) imposed everywhere) 
    236232      !-------------------- 
    237       z1_CpR = 1._wp / ( cpic * rhosn ) 
    238       zL_Cp  = lfus  /   cpic 
    239233      zlay_s = REAL( nlay_s , wp ) 
     234      z1_cp  = 1._wp / cpic 
    240235      DO jk = 1, nlay_s 
    241236         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area 
    242             t_s(:,:,jk,:) = MAX(  -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp )  ) + rt0        
     237            t_s(:,:,jk,:) = MAX( -100._wp , MIN( z1_cp * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0 
    243238         ELSEWHERE                           !--- no ice 
    244239            t_s(:,:,jk,:) =  rt0 - 100._wp         ! impose 173.15 K (i.e. -100 C) 
    245240         END WHERE 
    246241      END DO 
    247 !!gm perhaps better like this (?) :  (if this compile, yes! one test instead of nlay_s tests) 
    248 !      WHERE( v_s(:,:,:) > epsi20 )        !--- icy area 
    249 !         DO jk = 1, nlay_s 
    250 !            t_s(:,:,jk,:) = MAX(  -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp )  ) + rt0        
    251 !         END DO 
    252 !      ELSEWHERE                           !--- no ice 
    253 !         DO jk = 1, nlay_s 
    254 !            t_s(:,:,jk,:) =  rt0 - 100._wp      ! impose 173.15 K (i.e. -100 C) 
    255 !         END DO 
    256 !      END WHERE 
    257 !!gm end better ? 
    258242 
    259243      ! integrated values  
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90

    r8488 r8498  
    5454      INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices 
    5555      REAL(wp) ::  z2da, z2db, ztmp, zrho1, zrho2, zmiss_val 
    56       REAL(wp) ::  zs12, zshear 
    5756      REAL(wp), DIMENSION(jpi,jpj)     ::  z2d, zswi, zmiss !  2D workspace 
    5857      REAL(wp), DIMENSION(jpi,jpj)     ::  zfb              ! ice freeboard 
    5958      REAL(wp), DIMENSION(jpi,jpj)     ::  zamask, zamask15 ! 15% concentration mask 
    60       REAL(wp), DIMENSION(jpi,jpj)     ::  zsig1, zsig2, zsig3 
    6159      REAL(wp), DIMENSION(jpi,jpj,jpl) ::  zswi2, zmiss2 
    6260      ! 
     
    196194      CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean       
    197195  
    198       CALL iom_put( "afxtot"     , afx_tot              )        ! concentration tendency (total) 
    199       CALL iom_put( "afxdyn"     , afx_dyn              )        ! concentration tendency (dynamics) 
    200       CALL iom_put( "afxthd"     , afx_thd              )        ! concentration tendency (thermo) 
    201  
    202196      CALL iom_put ('hfxthd'     , hfx_thd(:,:)         )   !   
    203197      CALL iom_put ('hfxdyn'     , hfx_dyn(:,:)         )   !   
     
    219213      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  
    220214 
    221 !!gm ====>>>>>  THIS should be moved in icerhg_evp    (generalize this everywhere it is possible and logic...) 
    222       ! specific outputs for EVP rheology 
    223       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN 
    224          zsig1(:,:) = 0._wp; zsig2(:,:) = 0._wp; zsig3(:,:) = 0._wp; 
    225          DO jj = 2, jpjm1 
    226             DO ji = 2, jpim1 
    227                zs12 = ( zswi(ji-1,jj) * stress12_i(ji-1,jj) + zswi(ji  ,jj-1) * stress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point 
    228                   &     zswi(ji  ,jj) * stress12_i(ji  ,jj) + zswi(ji-1,jj-1) * stress12_i(ji-1,jj-1) )  & 
    229                   &   / MAX( 1._wp, zswi(ji-1,jj) + zswi(ji,jj-1) + zswi(ji,jj) + zswi(ji-1,jj-1) ) 
    230  
    231                zshear = SQRT( stress2_i(ji,jj) * stress2_i(ji,jj) + 4._wp * zs12 * zs12 ) ! shear stress   
    232  
    233                z2da = zswi(ji,jj) / MAX( 1._wp, strength(ji,jj) ) 
    234  
    235 !!               zsig1(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002) 
    236 !!               zsig2(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002) 
    237 !!               zsig3(ji,jj) = z2da**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress 
    238 !!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 
    239                zsig1(ji,jj) = 0.5_wp * z2da * ( stress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015 
    240                zsig2(ji,jj) = 0.5_wp * z2da * ( zshear )                    ! shear stress 
    241                zsig3(ji,jj) = z2da**2 * ( ( stress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 
    242             END DO 
    243          END DO 
    244          CALL lbc_lnk_multi( zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. ) 
    245          CALL iom_put( "isig1" , zsig1 ) 
    246          CALL iom_put( "isig2" , zsig2 ) 
    247          CALL iom_put( "isig3" , zsig3 ) 
    248       ENDIF 
    249 !!gm  <<<<<<======= end 
    250  
    251215!!gm ====>>>>>  THIS should be moved where at_ip, vt_ip are computed fro the last time in the time-step  (limmpd) 
    252216      ! MV MP 2016 
     
    339303      IF ( iom_use( "vice_mv"  ) ) CALL iom_put( "vice_mv"     ,   v_ice(:,:)               * zswi(:,:) + zmiss(:,:) )   ! ice velocity v component 
    340304       
    341       IF ( iom_use( "xmtrpice" ) ) CALL iom_put( "xmtrpice"     ,  diag_xmtrp_ice(:,:)      )                            ! X-component of sea-ice mass transport (kg/s) 
    342       IF ( iom_use( "ymtrpice" ) ) CALL iom_put( "ymtrpice"     ,  diag_ymtrp_ice(:,:)      )                            ! Y-component of sea-ice mass transport  
    343  
    344       IF ( iom_use( "xmtrpsnw" ) ) CALL iom_put( "xmtrpsnw"     ,  diag_xmtrp_snw(:,:)      )                            ! X-component of snow mass transport (kg/s) 
    345       IF ( iom_use( "ymtrpsnw" ) ) CALL iom_put( "ymtrpsnw"     ,  diag_ymtrp_snw(:,:)      )                            ! Y-component of snow mass transport 
    346  
    347       IF ( iom_use( "xatrp"    ) ) CALL iom_put( "xatrp"        ,  diag_xatrp(:,:)              )                        ! X-component of ice area transport 
    348       IF ( iom_use( "yatrp"    ) ) CALL iom_put( "yatrp"        ,  diag_yatrp(:,:)              )                        ! Y-component of ice area transport 
    349  
    350305      IF ( iom_use( "utau_ice" ) ) CALL iom_put( "utau_ice"     ,  utau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (x) 
    351306      IF ( iom_use( "vtau_ice" ) ) CALL iom_put( "vtau_ice"     ,  vtau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (y) 
    352307 
    353       IF ( iom_use( "utau_oi"  ) ) CALL iom_put( "utau_oi"     ,   diag_utau_oi(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Ocean stress term in force balance (x) 
    354       IF ( iom_use( "vtau_oi"  ) ) CALL iom_put( "vtau_oi"     ,   diag_vtau_oi(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Ocean stress term in force balance (y) 
    355308 
    356309      IF ( iom_use( "icestr"   ) ) CALL iom_put( "icestr"      ,   strength(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Ice strength 
    357  
    358       IF ( iom_use( "dssh_dx"  ) ) CALL iom_put( "dssh_dx"     ,   diag_dssh_dx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Sea-surface tilt term in force balance (x) 
    359       IF ( iom_use( "dssh_dy"  ) ) CALL iom_put( "dssh_dy"     ,   diag_dssh_dy(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Sea-surface tilt term in force balance (y) 
    360  
    361       IF ( iom_use( "corstrx"  ) ) CALL iom_put( "corstrx"     ,   diag_corstrx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Coriolis force term in force balance (x) 
    362       IF ( iom_use( "corstry"  ) ) CALL iom_put( "corstry"     ,   diag_corstry(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Coriolis force term in force balance (y) 
    363  
    364       IF ( iom_use( "intstrx"  ) ) CALL iom_put( "intstrx"     ,   diag_intstrx(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Internal force term in force balance (x) 
    365       IF ( iom_use( "intstry"  ) ) CALL iom_put( "intstry"     ,   diag_intstry(:,:)        * zswi(:,:) + zmiss(:,:) )   ! Internal force term in force balance (y) 
    366  
    367       IF ( iom_use( "normstr"  ) ) CALL iom_put( "normstr"     ,   diag_sig1(:,:)           * zswi(:,:) + zmiss(:,:) )   ! Normal stress 
    368       IF ( iom_use( "sheastr"  ) ) CALL iom_put( "sheastr"     ,   diag_sig2(:,:)           * zswi(:,:) + zmiss(:,:) )   ! Shear stress 
    369310 
    370311      !-------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.