Changeset 8355


Ignore:
Timestamp:
2017-07-20T10:20:49+02:00 (3 years ago)
Author:
clem
Message:

simplifications

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

Legend:

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

    r8341 r8355  
    2828   PRIVATE 
    2929 
    30    PUBLIC   lim_column_sum 
    31    PUBLIC   lim_column_sum_energy 
    32    PUBLIC   lim_cons_check 
    3330   PUBLIC   lim_cons_hsm 
    3431   PUBLIC   lim_cons_final 
     
    4037   !!---------------------------------------------------------------------- 
    4138CONTAINS 
    42  
    43    SUBROUTINE lim_column_sum( ksum, pin, pout ) 
    44       !!------------------------------------------------------------------- 
    45       !!               ***  ROUTINE lim_column_sum *** 
    46       !! 
    47       !! ** Purpose : Compute the sum of xin over nsum categories 
    48       !! 
    49       !! ** Method  : Arithmetics 
    50       !! 
    51       !! ** Action  : Gets xin(ji,jj,jl) and computes xout(ji,jj) 
    52       !!--------------------------------------------------------------------- 
    53       INTEGER                   , INTENT(in   ) ::   ksum   ! number of categories/layers 
    54       REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pin    ! input field 
    55       REAL(wp), DIMENSION(:,:)  , INTENT(  out) ::   pout   ! output field 
    56       ! 
    57       INTEGER ::   jl   ! dummy loop indices 
    58       !!--------------------------------------------------------------------- 
    59       ! 
    60       pout(:,:) = pin(:,:,1) 
    61       DO jl = 2, ksum 
    62          pout(:,:) = pout(:,:) + pin(:,:,jl) 
    63       END DO 
    64       ! 
    65    END SUBROUTINE lim_column_sum 
    66  
    67  
    68    SUBROUTINE lim_column_sum_energy( ksum, klay, pin, pout) 
    69       !!------------------------------------------------------------------- 
    70       !!               ***  ROUTINE lim_column_sum_energy *** 
    71       !! 
    72       !! ** Purpose : Compute the sum of xin over nsum categories 
    73       !!              and nlay layers 
    74       !! 
    75       !! ** Method  : Arithmetics 
    76       !!--------------------------------------------------------------------- 
    77       INTEGER                                , INTENT(in   ) ::   ksum   !: number of categories 
    78       INTEGER                                , INTENT(in   ) ::   klay   !: number of vertical layers 
    79       REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl), INTENT(in   ) ::   pin    !: input field 
    80       REAL(wp), DIMENSION(jpi,jpj)           , INTENT(  out) ::   pout   !: output field 
    81       ! 
    82       INTEGER ::   jk, jl   ! dummy loop indices 
    83       !!--------------------------------------------------------------------- 
    84       ! 
    85       pout(:,:) = 0._wp 
    86       DO jl = 1, ksum 
    87          DO jk = 2, klay  
    88             pout(:,:) = pout(:,:) + pin(:,:,jk,jl) 
    89          END DO 
    90       END DO 
    91       ! 
    92    END SUBROUTINE lim_column_sum_energy 
    93  
    94  
    95    SUBROUTINE lim_cons_check( px1, px2, pmax_err, cd_fieldid ) 
    96       !!------------------------------------------------------------------- 
    97       !!               ***  ROUTINE lim_cons_check *** 
    98       !! 
    99       !! ** Purpose : Test the conservation of a certain variable 
    100       !!              For each physical grid cell, check that initial  
    101       !!              and final values 
    102       !!              of a conserved field are equal to within a small value. 
    103       !! 
    104       !! ** Method  : 
    105       !!--------------------------------------------------------------------- 
    106       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   px1          !: initial field 
    107       REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   px2          !: final field 
    108       REAL(wp)                , INTENT(in   ) ::   pmax_err     !: max allowed error 
    109       CHARACTER(len=15)       , INTENT(in   ) ::   cd_fieldid   !: field identifyer 
    110       ! 
    111       INTEGER  ::   ji, jj          ! dummy loop indices 
    112       INTEGER  ::   inb_error       ! number of g.c where there is a cons. error 
    113       LOGICAL  ::   llconserv_err   ! = .true. if conservation check failed 
    114       REAL(wp) ::   zmean_error     ! mean error on error points 
    115       !!--------------------------------------------------------------------- 
    116       ! 
    117       IF(lwp) WRITE(numout,*) ' lim_cons_check ' 
    118       IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    119  
    120       llconserv_err = .FALSE. 
    121       inb_error     = 0 
    122       zmean_error   = 0._wp 
    123       IF( MAXVAL( px2(:,:) - px1(:,:) ) > pmax_err )   llconserv_err = .TRUE. 
    124  
    125       IF( llconserv_err ) THEN 
    126          DO jj = 1, jpj  
    127             DO ji = 1, jpi 
    128                IF( ABS( px2(ji,jj) - px1(ji,jj) ) > pmax_err ) THEN 
    129                   inb_error   = inb_error + 1 
    130                   zmean_error = zmean_error + ABS( px2(ji,jj) - px1(ji,jj) ) 
    131                   ! 
    132                   IF(lwp) THEN 
    133                      WRITE (numout,*) ' ALERTE 99 ' 
    134                      WRITE (numout,*) ' Conservation error: ', cd_fieldid 
    135                      WRITE (numout,*) ' Point             : ', ji, jj  
    136                      WRITE (numout,*) ' lat, lon          : ', gphit(ji,jj), glamt(ji,jj) 
    137                      WRITE (numout,*) ' Initial value     : ', px1(ji,jj) 
    138                      WRITE (numout,*) ' Final value       : ', px2(ji,jj) 
    139                      WRITE (numout,*) ' Difference        : ', px2(ji,jj) - px1(ji,jj) 
    140                   ENDIF 
    141                ENDIF 
    142             END DO 
    143          END DO 
    144          ! 
    145       ENDIF 
    146       IF(lk_mpp)   CALL mpp_sum( inb_error   ) 
    147       IF(lk_mpp)   CALL mpp_sum( zmean_error ) 
    148       ! 
    149       IF( inb_error > 0 .AND. lwp ) THEN 
    150          zmean_error = zmean_error / REAL( inb_error, wp ) 
    151          WRITE(numout,*) ' Conservation check for : ', cd_fieldid 
    152          WRITE(numout,*) ' Number of error points : ', inb_error 
    153          WRITE(numout,*) ' Mean error on these pts: ', zmean_error 
    154       ENDIF 
    155       ! 
    156    END SUBROUTINE lim_cons_check 
    157  
    15839 
    15940   SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r8233 r8355  
    4444CONTAINS 
    4545 
    46    SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 
     46   SUBROUTINE lim_itd_th_rem( kt ) 
    4747      !!------------------------------------------------------------------ 
    4848      !!                ***  ROUTINE lim_itd_th_rem *** 
     
    5555      !! References : W.H. Lipscomb, JGR 2001 
    5656      !!------------------------------------------------------------------ 
    57       INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point 
    58       INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied 
    5957      INTEGER , INTENT (in) ::   kt      ! Ocean time step  
    6058      ! 
     
    6563      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      - 
    6664      REAL(wp) ::   zx3         
    67       CHARACTER (len = 15) :: fieldid 
    6865 
    6966      INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor   ! donor category index 
     
    7572      REAL(wp), POINTER, DIMENSION(:,:,:) ::   hR          ! left boundary for the ITD for each thickness 
    7673      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zht_i_b     ! old ice thickness 
    77       REAL(wp), POINTER, DIMENSION(:,:,:) ::   dummy_es 
    7874      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice          ! local increment of ice area and volume 
    79       REAL(wp), POINTER, DIMENSION(:)     ::   zvetamin, zvetamax      ! maximum values for etas 
    8075      INTEGER , POINTER, DIMENSION(:)     ::   nind_i, nind_j          ! compressed indices for i/j directions 
    8176      INTEGER                             ::   nbrem                   ! number of cells with ice to transfer 
    8277      REAL(wp)                            ::   zslope                  ! used to compute local thermodynamic "speeds" 
    8378      REAL(wp), POINTER, DIMENSION(:,:)   ::   zhb0, zhb1              ! category boundaries for thinnes categories 
    84       REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_i_init, vt_i_final   !  ice volume summed over categories 
    85       REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_s_init, vt_s_final   !  snow volume summed over categories 
    86       REAL(wp), POINTER, DIMENSION(:,:)   ::   et_i_init, et_i_final   !  ice energy summed over categories 
    87       REAL(wp), POINTER, DIMENSION(:,:)   ::   et_s_init, et_s_final   !  snow energy summed over categories 
    8879      INTEGER , POINTER, DIMENSION(:,:)   ::   zremap_flag      ! compute remapping or not ???? 
    8980      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zhbnew           ! new boundaries of ice categories 
     
    9283      CALL wrk_alloc( jpi,jpj, zremap_flag ) 
    9384      CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 
    94       CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
     85      CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b ) 
    9586      CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    9687      CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    97       CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )    
    9888      CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )  
    99       CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
    100  
    101       !!---------------------------------------------------------------------------------------------- 
    102       !! 0) Conservation checkand changes in each ice category 
    103       !!---------------------------------------------------------------------------------------------- 
    104       IF( con_i ) THEN 
    105          CALL lim_column_sum (jpl,   v_i, vt_i_init) 
    106          CALL lim_column_sum (jpl,   v_s, vt_s_init) 
    107          CALL lim_column_sum_energy (jpl, nlay_i,   e_i, et_i_init) 
    108          dummy_es(:,:,:) = e_s(:,:,1,:) 
    109          CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init) 
    110       ENDIF 
     89      CALL wrk_alloc( jpi,jpj, zhb0, zhb1 ) 
    11190 
    11291      !!---------------------------------------------------------------------------------------------- 
     
    11796         WRITE(numout,*) 'lim_itd_th_rem  : Remapping the ice thickness distribution' 
    11897         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    119          WRITE(numout,*) ' klbnd :       ', klbnd 
    120          WRITE(numout,*) ' kubnd :       ', kubnd 
    12198      ENDIF 
    12299 
    123100      zdhice(:,:,:) = 0._wp 
    124       DO jl = klbnd, kubnd 
     101      DO jl = 1, jpl 
    125102         DO jj = 1, jpj 
    126103            DO ji = 1, jpi 
     
    134111      END DO 
    135112 
    136       !----------------------------------------------------------------------------------------------- 
    137       !  2) Compute fractional ice area in each grid cell 
    138       !----------------------------------------------------------------------------------------------- 
    139       at_i(:,:) = 0._wp 
    140       DO jl = klbnd, kubnd 
    141          at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
    142       END DO 
     113!!      !----------------------------------------------------------------------------------------------- 
     114!!      !  2) Compute fractional ice area in each grid cell 
     115!!      !----------------------------------------------------------------------------------------------- 
     116!!      at_i(:,:) = 0._wp 
     117!!      DO jl = 1, jpl 
     118!!         at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 
     119!!      END DO 
    143120 
    144121      !----------------------------------------------------------------------------------------------- 
     
    160137 
    161138      !----------------------------------------------------------------------------------------------- 
    162       !  4) Compute new category boundaries 
    163       !----------------------------------------------------------------------------------------------- 
    164       !- 4.1 Compute category boundaries 
     139      !  4) Compute new category boundaries: 1:jpl-1 
     140      !----------------------------------------------------------------------------------------------- 
    165141      zhbnew(:,:,:) = 0._wp 
    166142 
    167       DO jl = klbnd, kubnd - 1 
    168          DO ji = 1, nbrem 
    169             ii = nind_i(ji) 
    170             ij = nind_j(ji) 
    171             ! 
    172             zhbnew(ii,ij,jl) = hi_max(jl) 
    173             IF    ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 
    174                !interpolate between adjacent category growth rates 
    175                zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 
    176                zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 
    177             ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN 
    178                zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 
    179             ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN 
    180                zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 
    181             ENDIF 
    182          END DO 
    183  
    184          !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 
    185          DO ji = 1, nbrem 
    186             ii = nind_i(ji) 
    187             ij = nind_j(ji) 
    188  
    189             ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible  
    190             ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
    191             IF    ( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN 
    192                zremap_flag(ii,ij) = 0 
    193             ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN 
    194                zremap_flag(ii,ij) = 0 
    195             ENDIF 
    196  
    197             !- 4.3 Check that each zhbnew does not exceed maximal values hi_max   
    198             IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
    199             IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 
    200             ! clem bug: why is not the following instead? 
    201             !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 
    202             !!IF( zhbnew(ii,ij,jl) > hi_max(jl  ) ) zremap_flag(ii,ij) = 0 
    203   
    204          END DO 
    205  
    206       END DO 
    207  
     143      IF( nbrem > 0 ) THEN 
     144         DO jl = 1, jpl - 1 
     145            DO ji = 1, nbrem 
     146               ii = nind_i(ji) 
     147               ij = nind_j(ji) 
     148               ! 
     149               ! --- New boundary categories Hn* = Hn + Fn*dt --- 
     150               !!    Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - ht_i_b) 
     151               IF    ( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN !! a_i_b(jl+1) & a_i_b(jl) /= 0 
     152                  zslope           = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 
     153                  zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 
     154!!               ELSEIF( a_i_b(ii,ij,jl) >  epsi10 .AND. a_i_b(ii,jj,jl+1) <= epsi10 ) THEN !! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt 
     155               ELSEIF( a_i_b(ii,ij,jl) >  epsi10 ) THEN !! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt 
     156                  zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 
     157!!               ELSEIF( a_i_b(ii,ij,jl) <= epsi10 .AND. a_i_b(ii,ij,jl+1) >  epsi10 ) THEN !! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt 
     158               ELSEIF( a_i_b(ii,ij,jl+1) > epsi10 ) THEN !! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt 
     159                  zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 
     160               ELSE                                                                       !! a_i_b(jl+1) & a_i_b(jl) = 0  
     161                  zhbnew(ii,ij,jl) = hi_max(jl) 
     162               ENDIF 
     163             
     164               ! --- 2 conditions for remapping --- ! 
     165               ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi                
     166               !    Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible  
     167               !          in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     168               IF( a_i(ii,ij,jl  ) > epsi10 .AND. ht_i(ii,ij,jl  ) > ( zhbnew(ii,ij,jl) - epsi10 ) )   zremap_flag(ii,ij) = 0 
     169               IF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) )   zremap_flag(ii,ij) = 0 
     170                
     171               ! 2) Hn-1 < Hn* < Hn+1   
     172               IF( zhbnew(ii,ij,jl) < hi_max(jl-1) )   zremap_flag(ii,ij) = 0 
     173               IF( zhbnew(ii,ij,jl) > hi_max(jl+1) )   zremap_flag(ii,ij) = 0 
     174                
     175            END DO 
     176         END DO 
     177      ENDIF 
     178    
    208179      !----------------------------------------------------------------------------------------------- 
    209180      !  5) Identify cells where ITD is to be remapped 
     
    221192 
    222193      !----------------------------------------------------------------------------------------------- 
    223       !  6) Fill arrays with lowermost / uppermost boundaries of 'new' categories 
     194      !  6) Compute new category boundaries: jpl 
    224195      !----------------------------------------------------------------------------------------------- 
    225196      DO jj = 1, jpj 
     
    228199            zhb1(ji,jj) = hi_max(1) 
    229200 
    230             IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 
    231                zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) 
     201            ! boundary for jpl 
     202            IF( a_i(ji,jj,jpl) > epsi10 ) THEN 
     203               zhbnew(ji,jj,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i(ji,jj,jpl) - 2._wp * zhbnew(ji,jj,jpl-1) ) 
    232204            ELSE 
    233 !clem bug               zhbnew(ji,jj,kubnd) = hi_max(kubnd)   
    234                zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway 
    235             ENDIF 
    236  
    237             ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible  
    238             ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
    239             IF    ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) )  THEN 
    240                zremap_flag(ji,jj) = 0 
    241             ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) )  THEN 
    242                zremap_flag(ji,jj) = 0 
    243             ENDIF 
     205!clem bug               zhbnew(ji,jj,jpl) = hi_max(jpl)   
     206               zhbnew(ji,jj,jpl) = hi_max(jpl-1) ! not used anyway 
     207            ENDIF 
     208 
     209            ! 1 condition for remapping for the 1st category 
     210            ! H0+epsi < h1(t) < H1-epsi  
     211            !    h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible  
     212            !    in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 
     213            IF( zht_i_b(ji,jj,1) < ( zhb0(ji,jj) + epsi10 ) )   zremap_flag(ji,jj) = 0 
     214            IF( zht_i_b(ji,jj,1) > ( zhb1(ji,jj) - epsi10 ) )   zremap_flag(ji,jj) = 0 
    244215 
    245216         END DO 
     
    250221      !----------------------------------------------------------------------------------------------- 
    251222      !- 7.1 g(h) for category 1 at start of time step 
    252       CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   & 
    253          &                  hR(:,:,klbnd), zremap_flag ) 
    254  
    255       !- 7.2 Area lost due to melting of thin ice (first category,  klbnd) 
     223      CALL lim_itd_fitline( 1, zhb0, zhb1, zht_i_b(:,:,1), g0(:,:,1), g1(:,:,1), hL(:,:,1),   & 
     224         &                  hR(:,:,1), zremap_flag ) 
     225 
     226      !- 7.2 Area lost due to melting of thin ice (first category,  1) 
    256227      DO ji = 1, nbrem 
    257228         ii = nind_i(ji)  
    258229         ij = nind_j(ji)  
    259230 
    260          IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 
    261  
    262             zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 
     231         IF( a_i(ii,ij,1) > epsi10 ) THEN 
     232 
     233            zdh0 = zdhice(ii,ij,1) !decrease of ice thickness in the lower category 
    263234 
    264235            IF( zdh0 < 0.0 ) THEN      !remove area from category 1 
    265                zdh0 = MIN( -zdh0, hi_max(klbnd) ) 
     236               zdh0 = MIN( -zdh0, hi_max(1) ) 
    266237               !Integrate g(1) from 0 to dh0 to estimate area melted 
    267                zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 
     238               zetamax = MIN( zdh0, hR(ii,ij,1) ) - hL(ii,ij,1) 
    268239 
    269240               IF( zetamax > 0.0 ) THEN 
    270241                  zx1    = zetamax 
    271242                  zx2    = 0.5 * zetamax * zetamax  
    272                   zda0   = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1                        ! ice area removed 
    273                   zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i                 
     243                  zda0   = g1(ii,ij,1) * zx2 + g0(ii,ij,1) * zx1                        ! ice area removed 
     244                  zdamax = a_i(ii,ij,1) * (1.0 - ht_i(ii,ij,1) / zht_i_b(ii,ij,1) ) ! Constrain new thickness <= ht_i                 
    274245                  zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting  
    275246                                                                                                !     of thin ice (zdamax > 0) 
    276247                  ! Remove area, conserving volume 
    277                   ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 
    278                   a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0 
    279                   v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? 
     248                  ht_i(ii,ij,1) = ht_i(ii,ij,1) * a_i(ii,ij,1) / ( a_i(ii,ij,1) - zda0 ) 
     249                  a_i(ii,ij,1)  = a_i(ii,ij,1) - zda0 
     250                  v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) ! clem-useless ? 
    280251               ENDIF 
    281252 
    282253            ELSE ! if ice accretion zdh0 > 0 
    283254               ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 
    284                zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) )  
     255               zhbnew(ii,ij,0) = MIN( zdh0, hi_max(1) )  
    285256            ENDIF 
    286257 
     
    290261 
    291262      !- 7.3 g(h) for each thickness category   
    292       DO jl = klbnd, kubnd 
     263      DO jl = 1, jpl 
    293264         CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  & 
    294265            &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 
     
    296267 
    297268      !----------------------------------------------------------------------------------------------- 
    298       !  8) Compute area and volume to be shifted across each boundary 
    299       !----------------------------------------------------------------------------------------------- 
    300  
    301       DO jl = klbnd, kubnd - 1 
     269      !  8) Compute area and volume to be shifted across each boundary (Eq. 18) 
     270      !----------------------------------------------------------------------------------------------- 
     271 
     272      DO jl = 1, jpl - 1 
    302273         DO jj = 1, jpj 
    303274            DO ji = 1, jpi 
     
    312283            ij = nind_j(ji) 
    313284 
    314             IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 
    315                ! left and right integration limits in eta space 
    316                zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) 
    317                zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 
     285            ! left and right integration limits in eta space 
     286            IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 
     287               zetamin = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)       ! hi_max(jl) - hL  
     288               zetamax = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) ! hR - hL 
    318289               zdonor(ii,ij,jl) = jl 
    319  
    320             ELSE                                    ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 
    321                ! left and right integration limits in eta space 
    322                zvetamin(ji) = 0.0 
    323                zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1) 
     290            ELSE                                    ! Hn* <= Hn => transfer from jl+1 to jl 
     291               zetamin = 0.0 
     292               zetamax = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)   ! hi_max(jl) - hL 
    324293               zdonor(ii,ij,jl) = jl + 1 
    325  
    326             ENDIF 
    327  
    328             zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 
    329             zetamin = zvetamin(ji) 
     294            ENDIF 
     295            zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin 
    330296 
    331297            zx1  = zetamax - zetamin 
     
    346312      !! 9) Shift ice between categories 
    347313      !!---------------------------------------------------------------------------------------------- 
    348       CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice ) 
     314      CALL lim_itd_shiftice ( zdonor, zdaice, zdvice ) 
    349315 
    350316      !!---------------------------------------------------------------------------------------------- 
     
    366332      END DO 
    367333 
    368       !!---------------------------------------------------------------------------------------------- 
    369       !! 11) Conservation check 
    370       !!---------------------------------------------------------------------------------------------- 
    371       IF ( con_i ) THEN 
    372          CALL lim_column_sum (jpl,   v_i, vt_i_final) 
    373          fieldid = ' v_i : limitd_th ' 
    374          CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)  
    375  
    376          CALL lim_column_sum_energy (jpl, nlay_i,  e_i, et_i_final) 
    377          fieldid = ' e_i : limitd_th ' 
    378          CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid)  
    379  
    380          CALL lim_column_sum (jpl,   v_s, vt_s_final) 
    381          fieldid = ' v_s : limitd_th ' 
    382          CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)  
    383  
    384          dummy_es(:,:,:) = e_s(:,:,1,:) 
    385          CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final) 
    386          fieldid = ' e_s : limitd_th ' 
    387          CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)  
    388       ENDIF 
    389  
    390334      CALL wrk_dealloc( jpi,jpj, zremap_flag ) 
    391335      CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 
    392       CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 
     336      CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b ) 
    393337      CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )    
    394338      CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )    
    395       CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )    
    396339      CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )  
    397       CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
     340      CALL wrk_dealloc( jpi,jpj, zhb0, zhb1 ) 
    398341 
    399342   END SUBROUTINE lim_itd_th_rem 
     
    404347      !!                ***  ROUTINE lim_itd_fitline *** 
    405348      !! 
    406       !! ** Purpose :   fit g(h) with a line using area, volume constraints 
    407       !! 
    408       !! ** Method  :   Fit g(h) with a line, satisfying area and volume constraints. 
    409       !!              To reduce roundoff errors caused by large values of g0 and g1, 
    410       !!              we actually compute g(eta), where eta = h - hL, and hL is the 
    411       !!              left boundary. 
     349      !! ** Purpose :   build g(h) satisfying area and volume constraints (Eq. 6 and 9) 
     350      !! 
     351      !! ** Method  :   g(h) is linear and written as: g(eta) = g1(eta) + g0 
     352      !!                with eta = h - HL 
     353      !!                 
    412354      !!------------------------------------------------------------------ 
    413355      INTEGER                     , INTENT(in   ) ::   num_cat      ! category index 
     
    436378               hR(ji,jj) = HbR(ji,jj) 
    437379 
    438                ! Change hL or hR if hice falls outside central third of range 
    439                zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 
    440                zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 
    441  
    442                IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) 
    443                ELSEIF( hice(ji,jj) > zh23 ) THEN   ;   hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj) 
     380               ! Change hL or hR if hice falls outside central third of range, 
     381               ! so that hice is in the central third of the range [HL HR] 
     382               zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) +       hR(ji,jj) ) 
     383               zh23 = 1.0 / 3.0 * (       hL(ji,jj) + 2.0 * hR(ji,jj) ) 
     384 
     385               IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) ! move HR to the left 
     386               ELSEIF( hice(ji,jj) > zh23 ) THEN   ;   hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj) ! move HL to the right 
    444387               ENDIF 
    445388 
     
    448391               zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 
    449392               zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 
    450                g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) 
    451                g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 
     393               g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )      ! Eq. 14 
     394               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14 
    452395               ! 
    453396            ELSE  ! remap_flag = .false. or a_i < epsi10  
     
    464407 
    465408 
    466    SUBROUTINE lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
     409   SUBROUTINE lim_itd_shiftice( zdonor, zdaice, zdvice ) 
    467410      !!------------------------------------------------------------------ 
    468411      !!                ***  ROUTINE lim_itd_shiftice *** 
     
    473416      !! ** Method  : 
    474417      !!------------------------------------------------------------------ 
    475       INTEGER                           , INTENT(in   ) ::   klbnd    ! Start thickness category index point 
    476       INTEGER                           , INTENT(in   ) ::   kubnd    ! End point on which the  the computation is applied 
    477418      INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in   ) ::   zdonor   ! donor category index 
    478419      REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) ::   zdaice   ! ice area transferred across boundary 
     
    508449      !---------------------------------------------------------------------------------------------- 
    509450 
    510       DO jl = klbnd, kubnd 
     451      DO jl = 1, jpl 
    511452         zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 
    512453      END DO 
     
    516457      !------------------------------------------------------------------------------- 
    517458 
    518       DO jl = klbnd, kubnd - 1 
     459      DO jl = 1, jpl - 1 
    519460         nbrem = 0 
    520461         DO jj = 1, jpj 
     
    539480            ENDIF 
    540481 
    541             !-------------- 
    542482            ! Ice areas 
    543             !-------------- 
    544483            a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 
    545484            a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 
    546485 
    547             !-------------- 
    548486            ! Ice volumes 
    549             !-------------- 
    550487            v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl)  
    551488            v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 
    552489 
    553             !-------------- 
    554490            ! Snow volumes 
    555             !-------------- 
    556491            zdvsnow        = v_s(ii,ij,jl1) * zworka(ii,ij) 
    557492            v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 
    558493            v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow  
    559494 
    560             !-------------------- 
    561495            ! Snow heat content   
    562             !-------------------- 
    563496            zdesnow            = e_s(ii,ij,1,jl1) * zworka(ii,ij) 
    564497            e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow 
    565498            e_s(ii,ij,1,jl2)   = e_s(ii,ij,1,jl2) + zdesnow 
    566499 
    567             !-------------- 
    568500            ! Ice age  
    569             !-------------- 
    570501            zdo_aice           = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 
    571502            oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice 
    572503            oa_i(ii,ij,jl2)    = oa_i(ii,ij,jl2) + zdo_aice 
    573504 
    574             !-------------- 
    575505            ! Ice salinity 
    576             !-------------- 
    577506            zdsm_vice          = smv_i(ii,ij,jl1) * zworka(ii,ij) 
    578507            smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice 
    579508            smv_i(ii,ij,jl2)   = smv_i(ii,ij,jl2) + zdsm_vice 
    580509 
    581             !--------------------- 
    582510            ! Surface temperature 
    583             !--------------------- 
    584511            zdaTsf             = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 
    585512            zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf 
     
    588515            ! MV MP 2016  
    589516            IF ( nn_pnd_scheme > 0 ) THEN 
    590             !--------------------- 
    591517            ! Pond fraction 
    592             !--------------------- 
    593518               zdapnd             = a_ip(ii,ij,jl1) * zdaice(ii,ij,jl) 
    594519               a_ip(ii,ij,jl1)    = a_ip(ii,ij,jl1) - zdapnd 
    595520               a_ip(ii,ij,jl2)    = a_ip(ii,ij,jl2) + zdapnd 
    596521 
    597             !--------------------- 
    598522            ! Pond volume 
    599             !--------------------- 
    600523               zdvpnd             = v_ip(ii,ij,jl1) * zdaice(ii,ij,jl) 
    601524               v_ip(ii,ij,jl1)    = v_ip(ii,ij,jl1) - zdvpnd 
     
    607530         END DO 
    608531 
    609          !------------------ 
    610532         ! Ice heat content 
    611          !------------------ 
    612  
    613533         DO jk = 1, nlay_i 
    614534            DO ji = 1, nbrem 
     
    635555      !----------------------------------------------------------------- 
    636556 
    637       DO jl = klbnd, kubnd 
     557      DO jl = 1, jpl 
    638558         DO jj = 1, jpj 
    639559            DO ji = 1, jpi  
     
    656576    
    657577 
    658    SUBROUTINE lim_itd_th_reb( klbnd, kubnd ) 
     578   SUBROUTINE lim_itd_th_reb 
    659579      !!------------------------------------------------------------------ 
    660580      !!                ***  ROUTINE lim_itd_th_reb *** 
     
    664584      !! ** Method  : 
    665585      !!------------------------------------------------------------------ 
    666       INTEGER , INTENT (in) ::   klbnd   ! Start thickness category index point 
    667       INTEGER , INTENT (in) ::   kubnd   ! End point on which the  the computation is applied 
    668       ! 
    669586      INTEGER ::   ji,jj, jl   ! dummy loop indices 
    670587      INTEGER ::   zshiftflag          ! = .true. if ice must be shifted 
    671       CHARACTER (len = 15) :: fieldid 
    672588 
    673589      INTEGER , POINTER, DIMENSION(:,:,:) ::   zdonor           ! donor category index 
    674590      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zdaice, zdvice   ! ice area and volume transferred 
    675  
    676       REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final   ! ice volume summed over categories 
    677       REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   ! snow volume summed over categories 
    678591      !!------------------------------------------------------------------ 
    679592       
    680593      CALL wrk_alloc( jpi,jpj,jpl, zdonor )   ! interger 
    681594      CALL wrk_alloc( jpi,jpj,jpl, zdaice, zdvice ) 
    682       CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) 
    683595      !      
    684       IF( con_i ) THEN                 ! conservation check 
    685          CALL lim_column_sum (jpl,   v_i, vt_i_init) 
    686          CALL lim_column_sum (jpl,   v_s, vt_s_init) 
    687       ENDIF 
    688  
    689       ! 
    690596      !------------------------------------------------------------------------------ 
    691597      ! 1) Compute ice thickness. 
    692598      !------------------------------------------------------------------------------ 
    693       DO jl = klbnd, kubnd 
     599      DO jl = 1, jpl 
    694600         DO jj = 1, jpj 
    695601            DO ji = 1, jpi  
     
    707613      ! Initialize shift arrays 
    708614      !------------------------- 
    709       DO jl = klbnd, kubnd 
     615      DO jl = 1, jpl 
    710616         zdonor(:,:,jl) = 0 
    711617         zdaice(:,:,jl) = 0._wp 
     
    717623      !------------------------- 
    718624 
    719       DO jl = klbnd, kubnd - 1  ! loop over category boundaries 
     625      DO jl = 1, jpl - 1  ! loop over category boundaries 
    720626 
    721627         !--------------------------------------- 
     
    742648 
    743649         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
    744             CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
     650            CALL lim_itd_shiftice( zdonor, zdaice, zdvice ) 
    745651            ! Reset shift parameters 
    746652            zdonor(:,:,jl) = 0 
     
    755661      !---------------------------- 
    756662 
    757       DO jl = kubnd - 1, 1, -1       ! loop over category boundaries 
     663      DO jl = jpl - 1, 1, -1       ! loop over category boundaries 
    758664 
    759665         !----------------------------------------- 
     
    777683          
    778684         IF( zshiftflag == 1 ) THEN            ! Shift ice between categories 
    779             CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 
     685            CALL lim_itd_shiftice( zdonor, zdaice, zdvice ) 
    780686            ! Reset shift parameters 
    781687            zdonor(:,:,jl) = 0 
     
    785691 
    786692      END DO 
    787  
    788       !------------------------------------------------------------------------------ 
    789       ! 3) Conservation check 
    790       !------------------------------------------------------------------------------ 
    791  
    792       IF( con_i ) THEN 
    793          CALL lim_column_sum (jpl,   v_i, vt_i_final) 
    794          fieldid = ' v_i : limitd_reb ' 
    795          CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)  
    796  
    797          CALL lim_column_sum (jpl,   v_s, vt_s_final) 
    798          fieldid = ' v_s : limitd_reb ' 
    799          CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)  
    800       ENDIF 
    801693      ! 
    802694      CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 
    803695      CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 
    804       CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) 
    805696 
    806697   END SUBROUTINE lim_itd_th_reb 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r8342 r8355  
    319319      IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    320320 
    321       IF( jpl > 1 )      CALL lim_itd_th_rem( 1, jpl, kt ) 
     321      IF( jpl > 1 )      CALL lim_itd_th_rem( kt ) 
    322322 
    323323      IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r8316 r8355  
    104104      ! Rebin categories with thickness out of bounds 
    105105      !---------------------------------------------------- 
    106       IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 
     106      IF ( jpl > 1 ) CALL lim_itd_th_reb 
    107107 
    108108      !----------------- 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r8316 r8355  
    118118      ! Rebin categories with thickness out of bounds 
    119119      !---------------------------------------------------- 
    120       IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl ) 
     120      IF ( jpl > 1 ) CALL lim_itd_th_reb 
    121121 
    122122      !----------------- 
Note: See TracChangeset for help on using the changeset viewer.