Changeset 9452


Ignore:
Timestamp:
2018-03-30T18:12:48+02:00 (3 years ago)
Author:
clem
Message:

ensure perfect conservation in ridging/rafting routine which failed sometimes because of unrealistic ice thicknesses from advection scheme (>99m). Maybe we should prevent it to happen but it would mean hiding a defect

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90

    r9271 r9452  
    106106      !! 
    107107      !! ** Method  :   Steps : 
    108       !!       1) Thickness categories boundaries, ice / o.w. concentrations 
    109       !!          Ridge preparation 
    110       !!       2) Dynamical inputs (closing rate, divu_adv, opning) 
    111       !!       3) Ridging iteration 
    112       !!       4) Ridging diagnostics 
    113       !!       5) Heat, salt and freshwater fluxes 
    114       !!       6) Compute increments of tate variables and come back to old values 
     108      !!       0) Identify grid cells with ice 
     109      !!       1) Calculate closing rate, divergence and opening 
     110      !!       2) Identify grid cells with ridging 
     111      !!       3) Start ridging iterations 
     112      !!          - prep = ridged and rafted ice + closing_gross 
     113      !!          - shift = move ice from one category to another 
     114      !! 
     115      !! ** Details 
     116      !!    step1: The net rate of closing is due to convergence and shear, based on Flato and Hibler (1995). 
     117      !!           The energy dissipation rate is equal to the net closing rate times the ice strength. 
     118      !! 
     119      !!    step3: The gross closing rate is equal to the first two terms (open 
     120      !!           water closing and thin ice ridging) without the third term 
     121      !!           (thick, newly ridged ice). 
    115122      !! 
    116123      !! References :   Flato, G. M., and W. D. Hibler III, 1995, JGR, 100, 18,611-18,626. 
     
    123130      !! 
    124131      !!     This routine is based on CICE code and authors William H. Lipscomb, 
    125       !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged 
     132      !!     and Elizabeth C. Hunke, LANL are gratefully acknowledged 
    126133      !!------------------------------------------------------------------- 
    127134      INTEGER, INTENT(in) ::   kt     ! number of iteration 
     
    153160      !-------------------------------- 
    154161      npti = 0   ;   nptidx(:) = 0 
     162      ipti = 0   ;   iptidx(:) = 0 
    155163      DO jj = 1, jpj 
    156164         DO ji = 1, jpi 
    157             IF ( SUM(a_i(ji,jj,:)) > 0._wp ) THEN 
     165            IF ( at_i(ji,jj) > 0._wp ) THEN 
    158166               npti           = npti + 1 
    159167               nptidx( npti ) = (jj - 1) * jpi + ji 
     
    161169         END DO 
    162170      END DO 
    163  
     171       
     172      !-------------------------------------------------------- 
     173      ! 1) Dynamical inputs (closing rate, divergence, opening) 
     174      !-------------------------------------------------------- 
    164175      IF( npti > 0 ) THEN 
    165176         
     
    173184 
    174185         DO ji = 1, npti 
    175             !-----------------------------------------------------------------------------! 
    176             ! 2) Dynamical inputs (closing rate, divu_adv, opning) 
    177             !-----------------------------------------------------------------------------! 
     186            ! closing_net = rate at which open water area is removed + ice area removed by ridging  
     187            !                                                        - ice area added in new ridges 
     188            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
    178189            ! 
    179             ! 2.1 closing_net 
    180             !----------------- 
    181             ! Compute the net rate of closing due to convergence  
    182             ! and shear, based on Flato and Hibler (1995). 
    183             !  
    184             ! The energy dissipation rate is equal to the net closing rate 
    185             ! times the ice strength. 
     190            ! divergence given by the advection scheme 
     191            !   (which may not be equal to divu as computed from the velocity field) 
     192            zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 
    186193            ! 
    187             ! NOTE: The NET closing rate is equal to the rate that open water  
    188             !  area is removed, plus the rate at which ice area is removed by  
    189             !  ridging, minus the rate at which area is added in new ridges. 
    190             !  The GROSS closing rate is equal to the first two terms (open 
    191             !  water closing and thin ice ridging) without the third term 
    192             !  (thick, newly ridged ice). 
    193  
    194             closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 
    195  
    196             ! 2.2 divu_adv 
    197             !-------------- 
    198             ! Compute divu_adv, the divergence rate given by the transport/ 
    199             ! advection scheme, which may not be equal to divu as computed  
    200             ! from the velocity field. 
    201             ! 
    202             ! If divu_adv < 0, make sure the closing rate is large enough 
    203             ! to give asum = 1.0 after ridging. 
    204             zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 
    205  
    206             IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) 
    207  
    208             ! 2.3 opning 
    209             !------------ 
    210             ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging. 
     194            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough 
     195            !                                                                                        ! to give asum = 1.0 after ridging 
     196            ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 
    211197            opning(ji) = closing_net(ji) + zdivu_adv(ji) 
    212198         END DO 
    213  
    214199         ! 
    215          !------------------------------------------------ 
    216          ! 2.1) Identify grid cells with nonzero ridging 
    217          !------------------------------------------------ 
    218          CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) 
    219          ipti = 0 ; iptidx(:) = 0 
     200         !------------------------------------ 
     201         ! 2) Identify grid cells with ridging 
     202         !------------------------------------ 
     203         CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 
     204 
    220205         DO ji = 1, npti 
    221             IF( SUM( ABS(apartf(ji,1:jpl)) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 
     206            IF( SUM( apartf(ji,1:jpl) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 
    222207               ipti = ipti + 1 
    223208               iptidx     (ipti)   = nptidx     (ji) 
    224                a_i_2d     (ipti,:) = a_i_2d     (ji,:)  ! adjust to new indices 
    225                v_i_2d     (ipti,:) = v_i_2d     (ji,:)  ! adjust to new indices 
    226                ato_i_1d   (ipti)   = ato_i_1d   (ji)    ! adjust to new indices 
    227                closing_net(ipti)   = closing_net(ji)    ! adjust to new indices 
    228                zdivu_adv  (ipti)   = zdivu_adv  (ji)    ! adjust to new indices 
    229                opning     (ipti)   = opning     (ji)    ! adjust to new indices 
     209               ! adjust to new indices 
     210               a_i_2d     (ipti,:) = a_i_2d     (ji,:) 
     211               v_i_2d     (ipti,:) = v_i_2d     (ji,:) 
     212               ato_i_1d   (ipti)   = ato_i_1d   (ji) 
     213               closing_net(ipti)   = closing_net(ji) 
     214               zdivu_adv  (ipti)   = zdivu_adv  (ji) 
     215               opning     (ipti)   = opning     (ji) 
    230216            ENDIF 
    231217         END DO 
    232          nptidx(:) = iptidx(:) 
    233          npti      = ipti 
    234218 
    235219      ENDIF 
    236220 
    237       !-------------- 
    238       ! Start ridging 
    239       !-------------- 
     221      ! grid cells with ridging 
     222      nptidx(:) = iptidx(:) 
     223      npti      = ipti 
     224 
     225      !----------------- 
     226      ! 3) Start ridging 
     227      !----------------- 
    240228      IF( npti > 0 ) THEN 
    241229          
    242          !----------- 
    243          ! 2D <==> 1D 
    244          !----------- 
    245          ! fields used but not modified 
    246          CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) ) 
    247          CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) ) 
    248          ! the following fields are modified in this routine 
    249          !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 
    250          !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) ) 
    251          !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) ) 
    252          CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 
    253          CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 
    254          CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 
    255          CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    256          CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    257          DO jl = 1, jpl 
    258             DO jk = 1, nlay_s 
    259                CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 
    260             END DO 
    261             DO jk = 1, nlay_i 
    262                CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 
    263             END DO 
    264          END DO 
    265          CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) ) 
    266          CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) ) 
    267          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) ) 
    268          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) ) 
    269          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    270          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
    271           
    272          !-----------------------------------------------------------------------------! 
    273          ! 3) Ridging iteration 
    274          !-----------------------------------------------------------------------------! 
     230         CALL ice_dyn_1d2d( 1 )            ! --- Move to 1D arrays --- ! 
     231 
    275232         iter            = 1 
    276233         iterate_ridging = 1       
    277          DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax ) 
    278  
    279             CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) 
    280  
    281             ! 3.2 Redistribute area, volume, and energy. 
    282             !-----------------------------------------------------------------------------! 
     234         !                                                        !----------------------! 
     235         DO WHILE( iterate_ridging > 0 .AND. iter < jp_itermax )  !  ridging iterations  ! 
     236            !                                                     !----------------------! 
     237            ! Calculate participation function (apartf) 
     238            !       and transfer      function 
     239            !       and closing_gross (+correction on opening) 
     240            CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d, closing_net ) 
     241 
     242            ! Redistribute area, volume, and energy between categories 
    283243            CALL rdgrft_shift 
    284244 
    285             ! 3.4 Do we keep on iterating? 
    286             !-----------------------------------------------------------------------------! 
    287             ! Check whether asum = 1.  If not (because the closing and opening 
    288             ! rates were reduced above), ridge again with new rates. 
     245            ! Do we keep on iterating? 
     246            !------------------------- 
     247            ! Check whether a_i + ato_i = 0 
     248            ! If not, because the closing and opening rates were reduced above, ridge again with new rates 
    289249            iterate_ridging = 0 
    290250            DO ji = 1, npti 
     
    307267         END DO 
    308268 
    309          !----------- 
    310          ! 1D <==> 2D 
    311          !----------- 
    312          CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 
    313          CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) ) 
    314          CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 
    315          CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 
    316          CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 
    317          CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 
    318          CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    319          CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    320          DO jl = 1, jpl 
    321             DO jk = 1, nlay_s 
    322                CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 
    323             END DO 
    324             DO jk = 1, nlay_i 
    325                CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 
    326             END DO 
    327          END DO 
    328          CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) ) 
    329          CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) ) 
    330          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) ) 
    331          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) ) 
    332          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    333          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
    334  
    335       ENDIF ! npti > 0 
     269         CALL ice_dyn_1d2d( 2 )            ! --- Move to 2D arrays --- ! 
     270 
     271      ENDIF 
    336272    
    337273      CALL ice_var_agg( 1 )  
     
    345281 
    346282 
    347    SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i ) 
     283   SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i, pclosing_net ) 
    348284      !!------------------------------------------------------------------- 
    349285      !!                ***  ROUTINE rdgrft_prep *** 
     
    352288      !! 
    353289      !! ** Method  :   Compute the thickness distribution of the ice and open water  
    354       !!              participating in ridging and of the resulting ridges. 
    355       !!------------------------------------------------------------------- 
    356       REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i  
     290      !!                participating in ridging and of the resulting ridges. 
     291      !!------------------------------------------------------------------- 
     292      REAL(wp), DIMENSION(:)  , INTENT(in) ::   pato_i, pclosing_net  
    357293      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pa_i, pv_i  
    358294      !! 
     
    372308      END WHERE 
    373309 
     310      ! 1) Participation function (apartf): a(h) = b(h).g(h) 
    374311      !----------------------------------------------------------------- 
    375       ! 1) Participation function: a(h) = b(h).g(h) (apartf) 
    376       !----------------------------------------------------------------- 
    377       ! Compute the participation function apartf; this is analogous to 
    378       ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 
    379       ! area lost from category n due to ridging/closing 
    380       ! apartf(n)   = total area lost due to ridging/closing 
    381       ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).  
    382       ! 
    383       ! The expressions for apartf are found by integrating b(h)g(h) between 
    384       ! the category boundaries. 
     312      ! Compute the participation function = total area lost due to ridging/closing 
     313      ! This is analogous to 
     314      !   a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 
     315      !   assuming b(h) = (2/Gstar) * (1 - G(h)/Gstar).  
     316      ! 
     317      ! apartf = integrating b(h)g(h) between the category boundaries 
    385318      ! apartf is always >= 0 and SUM(apartf(0:jpl))=1 
    386319      !----------------------------------------------------------------- 
     
    393326      ELSEWHERE                        ;   z1_asum(1:npti) = 0._wp 
    394327      END WHERE 
     328      ! 
    395329      ! Compute cumulative thickness distribution function 
    396330      ! Compute the cumulative thickness distribution function zGsum, 
    397331      ! where zGsum(n) is the fractional area in categories 0 to n. 
    398       ! initial value (in h = 0) equals open water area 
     332      ! initial value (in h = 0) = open water area 
    399333      zGsum(1:npti,-1) = 0._wp 
    400334      zGsum(1:npti,0 ) = pato_i(1:npti) * z1_asum(1:npti) 
     
    465399      ENDIF 
    466400 
    467       !----------------------------------------------------------------- 
    468401      ! 2) Transfer function 
    469402      !----------------------------------------------------------------- 
     
    479412      !  constrained by hrmin <= (zhmean + hi)/2. 
    480413      !  
    481       ! The maximum ridging thickness, hrmax, is determined by 
    482       !  zhmean and hrmin. 
     414      ! The maximum ridging thickness, hrmax, is determined by zhmean and hrmin. 
    483415      ! 
    484416      ! These modifications have the effect of reducing the ice strength 
    485       ! (relative to the Hibler formulation) when very thick ice is 
    486       ! ridging. 
     417      ! (relative to the Hibler formulation) when very thick ice is ridging. 
    487418      ! 
    488419      ! zaksum = net area removed/ total area removed 
     
    490421      !         net area removed = total area removed - area of new ridges 
    491422      !----------------------------------------------------------------- 
    492  
    493423      zfac = 1._wp / hi_hrft 
    494424      zaksum(1:npti) = apartf(1:npti,0) 
    495       ! Transfer function 
    496       DO jl = 1, jpl !all categories have a specific transfer function 
     425      ! 
     426      DO jl = 1, jpl 
    497427         DO ji = 1, npti 
    498428            IF ( apartf(ji,jl) > 0._wp ) THEN 
     
    515445      END DO 
    516446      ! 
    517       ! 3.1 closing_gross 
    518       !-----------------------------------------------------------------------------! 
    519       ! Based on the ITD of ridging and ridged ice, convert the net 
    520       !  closing rate to a gross closing rate.   
     447      ! 3) closing_gross 
     448      !----------------- 
     449      ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.   
    521450      ! NOTE: 0 < aksum <= 1 
    522       WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = closing_net(1:npti) / zaksum(1:npti) 
     451      WHERE( zaksum(1:npti) > 0._wp )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 
    523452      ELSEWHERE                         ;   closing_gross(1:npti) = 0._wp 
    524453      END WHERE 
    525454       
     455      ! correction to closing rate if excessive ice removal 
     456      !---------------------------------------------------- 
     457      ! Reduce the closing rate if more than 100% of any ice category would be removed 
     458      ! Reduce the opening rate in proportion 
     459      DO jl = 1, jpl 
     460         DO ji = 1, npti 
     461            zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 
     462            IF( zfac > pa_i(ji,jl) ) THEN 
     463               closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 
     464            ENDIF 
     465         END DO 
     466      END DO       
     467 
     468      ! 4) correction to opening if excessive open water removal 
     469      !--------------------------------------------------------- 
     470      ! Reduce the closing rate if more than 100% of the open water would be removed 
     471      ! Reduce the opening rate in proportion 
    526472      DO ji = 1, npti   
    527          ! correction to closing rate and opening if closing rate is excessive 
    528          !--------------------------------------------------------------------- 
    529          ! Reduce the closing rate if more than 100% of the open water  
    530          ! would be removed.  Reduce the opening rate proportionately. 
    531          zfac = ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice 
    532          IF    ( zfac < 0._wp .AND. zfac > - pato_i(ji) ) THEN                  ! would lead to negative ato_i 
     473         zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice 
     474         IF( zfac < 0._wp ) THEN           ! would lead to negative ato_i 
    533475            opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice  
    534          ELSEIF( zfac > 0._wp .AND. zfac > ( zasum(ji) - pato_i(ji) ) ) THEN  ! would lead to ato_i > asum 
     476         ELSEIF( zfac > zasum(ji) ) THEN   ! would lead to ato_i > asum 
    535477            opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice  
    536478         ENDIF 
    537479      END DO 
    538        
    539       ! correction to closing rate / opening if excessive ice removal 
    540       !--------------------------------------------------------------- 
    541       ! Reduce the closing rate if more than 100% of any ice category  
    542       ! would be removed.  Reduce the opening rate proportionately. 
    543       DO jl = 1, jpl 
    544          DO ji = 1, npti 
    545             zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 
    546 !!            IF( zfac > pa_i(ji,jl) .AND. zfac > 0._wp ) THEN 
    547             IF( zfac > pa_i(ji,jl) ) THEN 
    548                closing_gross(ji) = closing_gross(ji) * pa_i(ji,jl) / zfac 
    549             ENDIF 
    550          END DO 
    551       END DO       
    552480      ! 
    553481   END SUBROUTINE rdgrft_prep 
     
    561489      !! 
    562490      !! ** Method  :   Remove area, volume, and energy from each ridging category 
    563       !!              and add to thicker ice categories. 
     491      !!                and add to thicker ice categories. 
    564492      !!------------------------------------------------------------------- 
    565493      ! 
     
    581509      REAL(wp), DIMENSION(jpij,nlay_s) ::   esrdg     ! enth*volume of new ridges       
    582510      REAL(wp), DIMENSION(jpij,nlay_i) ::   eirdg     ! enth*volume of new ridges 
    583       !!------------------------------------------------------------------- 
    584  
    585       !------------------------------------------------------------------------------- 
    586       ! 1) Compute change in open water area due to closing and opening. 
    587       !------------------------------------------------------------------------------- 
     511      ! 
     512      INTEGER , DIMENSION(jpij) ::   itest_rdg, itest_rft   ! test for conservation 
     513      !!------------------------------------------------------------------- 
     514 
     515      ! 1) Change in open water area due to closing and opening 
     516      !-------------------------------------------------------- 
    588517      DO ji = 1, npti 
    589518         ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice ) 
    590519      END DO 
    591520       
    592       !----------------------------------------------------------------- 
    593       ! 2) compute categories participating to ridging/rafting (jl1)  
    594       !----------------------------------------------------------------- 
     521      ! 2) compute categories in which ice is removed (jl1)  
     522      !---------------------------------------------------- 
    595523      DO jl1 = 1, jpl 
    596524 
     
    599527         DO ji = 1, npti 
    600528 
    601             !------------------------------------------------ 
    602             ! 2.1) Identify grid cells with nonzero ridging 
    603             !------------------------------------------------ 
    604             IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 
     529            IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging 
    605530 
    606531               z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 
    607532 
    608                !-------------------------------------------------------------------- 
    609                ! 2.2) Compute area of ridging ice (airdg1) and of new ridge (airdg2) 
    610                !-------------------------------------------------------------------- 
     533               ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 
    611534               airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice 
    612535               airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice 
     
    615538               airft2(ji) = airft1 * hi_hrft 
    616539 
    617                !--------------------------------------------------------------- 
    618                ! 2.3) Compute ridging /rafting fractions, make sure afrdg <=1  
    619                !--------------------------------------------------------------- 
     540               ! ridging /rafting fractions 
    620541               afrdg = airdg1 * z1_ai(ji) 
    621542               afrft = airft1 * z1_ai(ji) 
     
    625546               ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji)   ! clem: if sst>0, then ersw <0 (is that possible?) 
    626547 
    627                !--------------------------------------------------------------- 
    628                ! 2.4) Compute ridging ice and new ridges (vi, vs, sm, oi, es, ei) 
    629                !--------------------------------------------------------------- 
     548               ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 
    630549               virdg1     = v_i_2d (ji,jl1)   * afrdg 
    631550               virdg2(ji) = v_i_2d (ji,jl1)   * afrdg * ( 1. + rn_porordg ) 
     
    651570               ENDIF 
    652571 
    653                !----------------------------------------------------------------- 
    654                ! 2.5) Ice-ocean exchanges 
    655                !----------------------------------------------------------------- 
     572               ! Ice-ocean exchanges associated with ice porosity 
    656573               wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids 
    657574               sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_rdtice 
     
    659576 
    660577               ! Put the snow lost by ridging into the ocean 
    661                !------------------------------------------------ 
    662578               !  Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow. 
    663579               wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg )   &   ! fresh water source for ocean 
     
    665581 
    666582               ! Put the melt pond water into the ocean 
    667                !------------------------------------------             
    668583               ! clem: I think the following lines must be commented since there 
    669584               !       is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.) 
     
    674589 
    675590               ! virtual salt flux to keep salinity constant 
    676                !------------------------------------------             
    677591               IF( nn_icesal /= 2 )  THEN 
    678592                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i 
     
    681595               ENDIF 
    682596 
    683  
    684                !------------------------------------------             
    685                ! 2.6) update jl1 (removing ridged/rafted area) 
    686                !------------------------------------------             
     597               ! Remove area, volume of new ridge to each category jl1 
     598               !------------------------------------------------------ 
    687599               a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1    - airft1 
    688600               v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1    - virft(ji) 
     
    705617                  afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji) 
    706618                  afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji) 
    707                   ! Compute ridging ice and new ridges for es 
     619                  ! Compute ridging /rafting ice and new ridges for es 
    708620                  esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg 
    709621                  esrft(ji,jk) = ze_s_2d (ji,jk,jl1) * afrft 
     
    711623                  hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg )   &                 ! heat sink for ocean (<0, W.m-2) 
    712624                     &                                - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice 
    713                   ! Update jl1 
     625                  ! 
     626                  ! Remove energy of new ridge to each category jl1 
     627                  !------------------------------------------------- 
    714628                  ze_s_2d(ji,jk,jl1) = ze_s_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )  
    715629               ENDIF 
     
    727641                  eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i 
    728642                  eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft 
    729                   ! Update jl1 
     643                  ! 
     644                  ! Remove energy of new ridge to each category jl1 
     645                  !------------------------------------------------- 
    730646                  ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft )  
    731647               ENDIF 
     
    733649         END DO 
    734650          
    735  
    736          !------------------------------------------------------------------------------- 
    737          ! 3) Add area, volume, and energy of new ridge to each category jl2 
    738          !------------------------------------------------------------------------------- 
     651         ! 3) compute categories in which ice is added (jl2)  
     652         !-------------------------------------------------- 
     653         itest_rdg(1:npti) = 0 
     654         itest_rft(1:npti) = 0 
    739655         DO jl2  = 1, jpl  
    740656            ! 
     
    743659               IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 
    744660 
    745                   ! Compute the fraction of ridged ice area and volume going to thickness category jl2. 
     661                  ! Compute the fraction of ridged ice area and volume going to thickness category jl2 
    746662                  IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN 
    747                      !----------------------------------------------------------------- 
    748                      ! 2.8 Compute quantities used to apportion ice among categories 
    749                      ! in the n2 loop below 
    750                      !----------------------------------------------------------------- 
    751663                     hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) 
    752664                     hR = MIN( hrmax(ji,jl1), hi_max(jl2)   ) 
    753665                     farea    = ( hR      - hL      ) / ( hrmax(ji,jl1)                 - hrmin(ji,jl1)                 ) 
    754666                     fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) ) 
     667                     ! 
     668                     itest_rdg(ji) = 1   ! test for conservation 
    755669                  ELSE 
    756670                     farea    = 0._wp  
     
    759673 
    760674                  ! Compute the fraction of rafted ice area and volume going to thickness category jl2 
    761                   !!gm see above               IF( hraft(ji) <= hi_max(jl2) .AND. hraft(ji) >  hi_max(jl2-1) ) THEN 
    762                   IF( hi_max(jl2-1) < hraft(ji,jl1) .AND. hraft(ji,jl1) <= hi_max(jl2)  ) THEN   ;   zswitch(ji) = 1._wp 
    763                   ELSE                                                                           ;   zswitch(ji) = 0._wp 
     675                  IF( hraft(ji,jl1) <= hi_max(jl2) .AND. hraft(ji,jl1) >  hi_max(jl2-1) ) THEN 
     676                     zswitch(ji) = 1._wp 
     677                     ! 
     678                     itest_rft(ji) = 1   ! test for conservation 
     679                  ELSE 
     680                     zswitch(ji) = 0._wp 
    764681                  ENDIF 
    765682                  ! 
     683                  ! Patch to ensure perfect conservation if ice thickness goes mad 
     684                  ! Sometimes thickness is larger than hi_max(jpl) because of advection scheme (for very small areas) 
     685                  ! Then ice volume is removed from one category but the ridging/rafting scheme 
     686                  ! does not know where to move it, leading to a conservation issue.   
     687                  IF( itest_rdg(ji) == 0 .AND. jl2 == jpl ) THEN   ;   farea = 1._wp   ;   fvol(ji) = 1._wp   ;   ENDIF 
     688                  IF( itest_rft(ji) == 0 .AND. jl2 == jpl )      zswitch(ji) = 1._wp 
     689                  ! 
     690                  ! Add area, volume of new ridge to category jl2 
     691                  !---------------------------------------------- 
    766692                  a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea    + airft2(ji) * zswitch(ji) ) 
    767693                  oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ( oirdg2(ji) * farea    + oirft2(ji) * zswitch(ji) ) 
     
    780706 
    781707            END DO 
    782             ! for snow enthalpy 
     708            ! Add snow energy of new ridge to category jl2 
     709            !--------------------------------------------- 
    783710            DO jk = 1, nlay_s 
    784711               DO ji = 1, npti 
     
    788715               END DO 
    789716            END DO 
    790             ! for ice enthalpy 
     717            ! Add ice energy of new ridge to category jl2 
     718            !-------------------------------------------- 
    791719            DO jk = 1, nlay_i 
    792720               DO ji = 1, npti 
     
    886814   END SUBROUTINE ice_strength 
    887815 
     816    
     817   SUBROUTINE ice_dyn_1d2d( kn ) 
     818      !!----------------------------------------------------------------------- 
     819      !!                   ***  ROUTINE ice_dyn_1d2d ***  
     820      !!                  
     821      !! ** Purpose :   move arrays from 1d to 2d and the reverse 
     822      !!----------------------------------------------------------------------- 
     823      INTEGER, INTENT(in) ::   kn   ! 1= from 2D to 1D   ;   2= from 1D to 2D 
     824      ! 
     825      INTEGER ::   jl, jk   ! dummy loop indices 
     826      !!----------------------------------------------------------------------- 
     827      ! 
     828      SELECT CASE( kn ) 
     829      !                    !---------------------! 
     830      CASE( 1 )            !==  from 2D to 1D  ==! 
     831         !                 !---------------------! 
     832         ! fields used but not modified 
     833         CALL tab_2d_1d( npti, nptidx(1:npti), sss_1d(1:npti), sss_m(:,:) ) 
     834         CALL tab_2d_1d( npti, nptidx(1:npti), sst_1d(1:npti), sst_m(:,:) ) 
     835         ! the following fields are modified in this routine 
     836         !!CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 
     837         !!CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d(1:npti,1:jpl), a_i(:,:,:) ) 
     838         !!CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i  (:,:,:) ) 
     839         CALL tab_3d_2d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 
     840         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 
     841         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 
     842         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
     843         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
     844         DO jl = 1, jpl 
     845            DO jk = 1, nlay_s 
     846               CALL tab_2d_1d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 
     847            END DO 
     848            DO jk = 1, nlay_i 
     849               CALL tab_2d_1d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 
     850            END DO 
     851         END DO 
     852         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) ) 
     853         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) ) 
     854         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) ) 
     855         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) ) 
     856         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
     857         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
     858 
     859         !                 !---------------------! 
     860      CASE( 2 )            !==  from 1D to 2D  ==! 
     861         !                 !---------------------! 
     862         CALL tab_1d_2d( npti, nptidx(1:npti), ato_i_1d(1:npti), ato_i(:,:) ) 
     863         CALL tab_2d_3d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:) ) 
     864         CALL tab_2d_3d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:) ) 
     865         CALL tab_2d_3d( npti, nptidx(1:npti), v_s_2d (1:npti,1:jpl), v_s (:,:,:) ) 
     866         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 
     867         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 
     868         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
     869         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
     870         DO jl = 1, jpl 
     871            DO jk = 1, nlay_s 
     872               CALL tab_1d_2d( npti, nptidx(1:npti), ze_s_2d(1:npti,jk,jl), e_s(:,:,jk,jl) ) 
     873            END DO 
     874            DO jk = 1, nlay_i 
     875               CALL tab_1d_2d( npti, nptidx(1:npti), ze_i_2d(1:npti,jk,jl), e_i(:,:,jk,jl) ) 
     876            END DO 
     877         END DO 
     878         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_dyn_1d    (1:npti), sfx_dyn    (:,:) ) 
     879         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bri_1d    (1:npti), sfx_bri    (:,:) ) 
     880         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_dyn_1d    (1:npti), wfx_dyn    (:,:) ) 
     881         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) ) 
     882         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
     883         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
     884         ! 
     885      END SELECT 
     886      ! 
     887   END SUBROUTINE ice_dyn_1d2d 
     888    
    888889 
    889890   SUBROUTINE ice_dyn_rdgrft_init 
Note: See TracChangeset for help on using the changeset viewer.