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 8559 – NEMO

Changeset 8559


Ignore:
Timestamp:
2017-09-22T16:55:24+02:00 (6 years ago)
Author:
clem
Message:

almost useless commits

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/ice1D.F90

    r8554 r8559  
    100100   !                                                     ! to reintegrate longwave flux inside the ice thermodynamics 
    101101   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice 
    102    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hicol_1d      !: Ice collection thickness accumulated in leads 
    103102 
    104103   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_su_1d       !: <==> the 2D  t_su 
     
    197196         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  & 
    198197         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  & 
    199          &      sfx_lam_1d (jpij) , sfx_dyn_1d(jpij), hicol_1d   (jpij)    , STAT=ierr(ii) ) 
     198         &      sfx_lam_1d (jpij) , sfx_dyn_1d(jpij)  , STAT=ierr(ii) ) 
    200199      ! 
    201200      ii = ii + 1 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

    r8534 r8559  
    293293      DO jk = 1, nlay_i 
    294294         DO ji = 1, nidx 
    295             ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K] 
    296              
    297             IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
     295            ztmelts           = - tmut * s_i_1d(ji,jk)          ! Melting point of layer k [C] 
     296             
     297            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting 
    298298 
    299299               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]        
     
    320320                
    321321               zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0] 
    322                zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0] 
     322               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0] 
    323323               zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0 
    324324                
     
    430430                                  + ( 1. - zswitch_sal ) * sm_i_1d(ji)  
    431431               ! New ice growth 
    432                ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) 
     432               ztmelts            = - tmut * s_i_new(ji)                 ! New ice melting point (C) 
    433433 
    434434               zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    435435                
    436                zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
    437                   &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & 
    438                   &               + rcp  * ( ztmelts-rt0 )           
     436               zEi                = cpic * ( zt_i_new - (ztmelts+rt0) ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
     437                  &               - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) )   & 
     438                  &               + rcp  * ztmelts           
    439439 
    440440               zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     
    450450            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0) 
    451451 
    452             ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K) 
     452            ztmelts        = - tmut * s_i_new(ji)                ! New ice melting point (C) 
    453453             
    454454            zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 
    455455             
    456             zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
    457                &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   & 
    458                &               + rcp  * ( ztmelts-rt0 )           
     456            zEi            = cpic * ( zt_i_new - (ztmelts+rt0) ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
     457               &               - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) )   & 
     458               &               + rcp  * ztmelts           
    459459             
    460460            zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     
    490490            IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting  
    491491 
    492                ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K) 
    493  
    494                IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 
     492               ztmelts = - tmut * s_i_1d(ji,jk)  ! Melting point of layer jk (C) 
     493 
     494               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting 
    495495 
    496496                  zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0) 
     
    520520 
    521521                  zEi             = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 
    522                   zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0) 
     522                  zEw             = rcp * ztmelts              ! Specific enthalpy of meltwater (J/kg, <0) 
    523523                  zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0) 
    524524 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90

    r8534 r8559  
    8080      INTEGER  ::   ji,jj,jk,jl      ! dummy loop indices 
    8181      INTEGER  ::   iter     !   -       - 
    82       REAL(wp) ::   ztmelts, zdv, zfrazb, zweight, zde                          ! local scalars 
     82      REAL(wp) ::   ztmelts, zfrazb, zweight, zde                          ! local scalars 
    8383      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      - 
    8484      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      - 
     
    102102      REAL(wp), DIMENSION(jpij) ::   zdv_res     ! residual volume in case of excessive heat budget 
    103103      REAL(wp), DIMENSION(jpij) ::   zda_res     ! residual area in case of excessive heat budget 
    104       REAL(wp), DIMENSION(jpij) ::   zat_i_1d    ! total ice fraction     
    105104      REAL(wp), DIMENSION(jpij) ::   zv_frazb    ! accretion of frazil ice at the ice bottom 
    106105      REAL(wp), DIMENSION(jpij) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector) 
     
    108107      REAL(wp), DIMENSION(jpij,jpl) ::   zv_b      ! old volume of ice in category jl 
    109108      REAL(wp), DIMENSION(jpij,jpl) ::   za_b      ! old area of ice in category jl 
    110       REAL(wp), DIMENSION(jpij,jpl) ::   za_i_1d   ! 1-D version of a_i 
    111       REAL(wp), DIMENSION(jpij,jpl) ::   zv_i_1d   ! 1-D version of v_i 
    112       REAL(wp), DIMENSION(jpij,jpl) ::   zsmv_i_1d ! 1-D version of smv_i 
    113  
    114       REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_1d !: 1-D version of e_i 
     109 
     110      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d !: 1-D version of e_i 
    115111 
    116112      REAL(wp), DIMENSION(jpi,jpj) ::   zvrel     ! relative ice / frazil velocity 
     
    243239      IF ( nidx > 0 ) THEN 
    244240 
    245          CALL tab_2d_1d( nidx, idxice(1:nidx), zat_i_1d  (1:nidx)     , at_i           ) 
     241         CALL tab_2d_1d( nidx, idxice(1:nidx), at_i_1d (1:nidx)      , at_i         ) 
     242         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i  (:,:,:) ) 
     243         CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) ) 
     244         CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) ) 
    246245         DO jl = 1, jpl 
    247             CALL tab_2d_1d( nidx, idxice(1:nidx), za_i_1d  (1:nidx,jl), a_i  (:,:,jl)  ) 
    248             CALL tab_2d_1d( nidx, idxice(1:nidx), zv_i_1d  (1:nidx,jl), v_i  (:,:,jl)  ) 
    249             CALL tab_2d_1d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i(:,:,jl)  ) 
    250246            DO jk = 1, nlay_i 
    251                CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl)   ) 
    252             END DO 
    253          END DO 
    254  
    255          CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d  (1:nidx)     , qlead       ) 
    256          CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d   (1:nidx)     , t_bo        ) 
    257          CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx)     , sfx_opw     ) 
    258          CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx)     , wfx_opw     ) 
    259          CALL tab_2d_1d( nidx, idxice(1:nidx), hicol_1d  (1:nidx)     , hicol       ) 
    260          CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx)     , zvrel       ) 
    261  
    262          CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx)     , hfx_thd     ) 
    263          CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx)     , hfx_opw     ) 
    264          CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx)     , rn_amax_2d  ) 
    265          CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d    (1:nidx)     , sss_m       ) 
     247               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) 
     248            END DO 
     249         END DO 
     250         CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d  (1:nidx) , qlead       ) 
     251         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d   (1:nidx) , t_bo        ) 
     252         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx) , sfx_opw     ) 
     253         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx) , wfx_opw     ) 
     254         CALL tab_2d_1d( nidx, idxice(1:nidx), zh_newice (1:nidx) , hicol       ) 
     255         CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx) , zvrel       ) 
     256 
     257         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx) , hfx_thd     ) 
     258         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx) , hfx_opw     ) 
     259         CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx) , rn_amax_2d  ) 
     260         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d    (1:nidx) , sss_m       ) 
    266261 
    267262         !------------------------------------------------------------------------------| 
     
    269264         !------------------------------------------------------------------------------| 
    270265         DO jl = 1, jpl 
    271             DO jk = 1, nlay_i 
    272                DO ji = 1, nidx 
    273                   IF( zv_i_1d(ji,jl) > 0._wp )   ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) / zv_i_1d(ji,jl) * REAL( nlay_i ) 
    274                END DO 
     266            DO jk = 1, nlay_i                
     267               WHERE( v_i_2d(1:nidx,jl) > 0._wp ) 
     268                  ze_i_2d(1:nidx,jk,jl) = ze_i_2d(1:nidx,jk,jl) / v_i_2d(1:nidx,jl) * REAL( nlay_i ) 
     269               ELSEWHERE 
     270                  ze_i_2d(1:nidx,jk,jl) = 0._wp 
     271               END WHERE 
    275272            END DO 
    276273         END DO 
     
    282279         ! Keep old ice areas and volume in memory 
    283280         !----------------------------------------- 
    284          zv_b(1:nidx,:) = zv_i_1d(1:nidx,:)  
    285          za_b(1:nidx,:) = za_i_1d(1:nidx,:) 
    286  
    287          !---------------------- 
    288          ! Thickness of new ice 
    289          !---------------------- 
    290          zh_newice(1:nidx) = hicol_1d(1:nidx) 
     281         zv_b(1:nidx,:) = v_i_2d(1:nidx,:)  
     282         za_b(1:nidx,:) = a_i_2d(1:nidx,:) 
    291283 
    292284         !---------------------- 
     
    309301         ! We assume that new ice is formed at the seawater freezing point 
    310302         DO ji = 1, nidx 
    311             ztmelts       = - tmut * zs_newice(ji) + rt0                  ! Melting point (K) 
    312             ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                                         & 
    313                &                       + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & 
    314                &                       - rcp  *         ( ztmelts - rt0 ) ) 
     303            ztmelts       = - tmut * zs_newice(ji)                  ! Melting point (C) 
     304            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     & 
     305               &                       + lfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   & 
     306               &                       - rcp  *         ztmelts ) 
    315307         END DO 
    316308 
     
    318310         ! Age of new ice 
    319311         !---------------- 
    320          DO ji = 1, nidx 
    321             zo_newice(ji) = 0._wp 
    322          END DO 
     312         zo_newice(1:nidx) = 0._wp 
    323313 
    324314         !------------------- 
     
    350340         END DO 
    351341          
    352          zv_frazb(:) = 0._wp 
     342         zv_frazb(1:nidx) = 0._wp 
    353343         IF( ln_frazil ) THEN 
    354344            ! A fraction zfrazb of frazil ice is accreted at the ice bottom 
    355345            DO ji = 1, nidx 
    356                rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 
     346               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) ) 
    357347               zfrazb        = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz 
    358348               zv_frazb(ji)  =         zfrazb   * zv_newice(ji) 
     
    378368         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
    379369         DO ji = 1, nidx 
    380             IF ( za_newice(ji) >  ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN 
    381                zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) ) 
     370            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN 
     371               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) ) 
    382372               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji)  
    383373               za_newice(ji) = za_newice(ji) - zda_res  (ji) 
     
    390380 
    391381         ! find which category to fill 
    392          zat_i_1d(:) = 0._wp 
    393382         DO jl = 1, jpl 
    394383            DO ji = 1, nidx 
    395384               IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN 
    396                   za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji) 
    397                   zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji) 
    398                   jcat    (ji)    = jl 
     385                  a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji) 
     386                  v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice(ji) 
     387                  jcat(ji) = jl 
    399388               ENDIF 
    400                zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl) 
    401             END DO 
    402          END DO 
     389            END DO 
     390         END DO 
     391         at_i_1d(1:nidx) = SUM( a_i_2d(1:nidx,:), dim=2 ) 
    403392 
    404393         ! Heat content 
     
    411400            DO ji = 1, nidx 
    412401               jl = jcat(ji) 
    413                rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 
    414                ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    & 
    415                   &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) )  & 
    416                   &        * rswitch / MAX( zv_i_1d(ji,jl), epsi20 ) 
     402               rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) ) 
     403               ze_i_2d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    & 
     404                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) )  & 
     405                  &        * rswitch / MAX( v_i_2d(ji,jl), epsi20 ) 
    417406            END DO 
    418407         END DO 
     
    428417            DO jk = 1, nlay_i 
    429418               DO ji = 1, nidx 
    430                   h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i 
    431                   eh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 
     419                  h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i 
     420                  eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk) 
    432421               END DO 
    433422            END DO 
     
    435424            ! new volumes including lateral/bottom accretion + residual 
    436425            DO ji = 1, nidx 
    437                rswitch        = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 
    438                zv_newfra      = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 
    439                za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl)                
    440                zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 
     426               rswitch        = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) ) 
     427               zv_newfra     = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 ) 
     428               a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)                
     429               v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra 
    441430               ! for remapping 
    442431               h_i_old (ji,nlay_i+1) = zv_newfra 
     
    444433            ENDDO 
    445434            ! --- Ice enthalpy remapping --- ! 
    446             CALL ice_thd_ent( ze_i_1d(1:nidx,:,jl) )  
     435            CALL ice_thd_ent( ze_i_2d(1:nidx,:,jl) )  
    447436         ENDDO 
    448437 
     
    452441         DO jl = 1, jpl 
    453442            DO ji = 1, nidx 
    454                zdv   = zv_i_1d(ji,jl) - zv_b(ji,jl) 
    455                zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 
     443               smv_i_2d(ji,jl) = smv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) ) 
    456444            END DO 
    457445         END DO 
     
    462450         DO jl = 1, jpl 
    463451            DO jk = 1, nlay_i 
    464                DO ji = 1, nidx 
    465                   ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) * zv_i_1d(ji,jl) * r1_nlay_i  
    466                END DO 
     452               ze_i_2d(1:nidx,jk,jl) = ze_i_2d(1:nidx,jk,jl) * v_i_2d(1:nidx,jl) * r1_nlay_i  
    467453            END DO 
    468454         END DO 
     
    470456         ! 7) Change 2D vectors to 1D vectors  
    471457         !------------------------------------------------------------------------------! 
    472          DO jl = 1, jpl 
    473             CALL tab_1d_2d( nidx, idxice(1:nidx), za_i_1d (1:nidx,jl), a_i (:,:,jl) ) 
    474             CALL tab_1d_2d( nidx, idxice(1:nidx), zv_i_1d (1:nidx,jl), v_i (:,:,jl) ) 
    475             CALL tab_1d_2d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i (:,:,jl)   ) 
     458         CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i  (:,:,:) ) 
     459         CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) ) 
     460         CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) ) 
     461          DO jl = 1, jpl 
    476462            DO jk = 1, nlay_i 
    477                CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) 
    478             END DO 
    479          END DO 
    480          CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw  ) 
    481          CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw  ) 
    482          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd  ) 
    483          CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw  ) 
     463               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) 
     464            END DO 
     465         END DO 
     466         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw ) 
     467         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw ) 
     468         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd ) 
     469         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw ) 
    484470         ! 
    485471      ENDIF ! nidx > 0 
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90

    r8534 r8559  
    5959      REAL(wp) ::   iflush, igravdr              ! local scalars 
    6060      REAL(wp) ::   zs_sni, zsm_i_gd, zsm_i_fl, zsm_i_si, zsm_i_bg   ! local scalars 
     61      REAL(wp) ::   z1_time_gd, z1_time_fl 
    6162      !!--------------------------------------------------------------------- 
    6263 
     
    6667      CASE( 2 )       !  time varying salinity with linear profile  ! 
    6768      !               !---------------------------------------------! 
     69         z1_time_gd = 1._wp / rn_time_gd * rdt_ice 
     70         z1_time_fl = 1._wp / rn_time_fl * rdt_ice 
     71         ! 
    6872         DO ji = 1, nidx 
    6973 
     
    7175            !  Update ice salinity from snow-ice and bottom growth 
    7276            !--------------------------------------------------------- 
    73             zs_sni   = sss_1d(ji) * ( rhoic - rhosn ) * r1_rhoic   ! Salinity of snow ice 
    74             rswitch  = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 
    75             zsm_i_si = ( zs_sni      - sm_i_1d(ji) ) *             dh_snowice(ji)  / MAX( ht_i_1d(ji), epsi20 ) * rswitch ! snow-ice     
    76             zsm_i_bg = ( s_i_new(ji) - sm_i_1d(ji) ) * MAX( 0._wp, dh_i_bott(ji) ) / MAX( ht_i_1d(ji), epsi20 ) * rswitch ! bottom growth 
    77  
    78             ! Update salinity (nb: salt flux already included in icethd_dh) 
    79             sm_i_1d(ji) = sm_i_1d(ji) + zsm_i_bg + zsm_i_si 
     77            IF( ht_i_1d(ji) > 0._wp ) THEN 
     78               zs_sni   = sss_1d(ji) * ( rhoic - rhosn ) * r1_rhoic                                 ! Salinity of snow ice 
     79               zsm_i_si = ( zs_sni      - sm_i_1d(ji) ) *             dh_snowice(ji)  / ht_i_1d(ji) ! snow-ice     
     80               zsm_i_bg = ( s_i_new(ji) - sm_i_1d(ji) ) * MAX( 0._wp, dh_i_bott(ji) ) / ht_i_1d(ji) ! bottom growth 
     81               ! Update salinity (nb: salt flux already included in icethd_dh) 
     82               sm_i_1d(ji) = sm_i_1d(ji) + zsm_i_bg + zsm_i_si 
     83            ENDIF 
    8084 
    8185            IF( ld_sal ) THEN 
     
    8589               iflush   = MAX( 0._wp , SIGN( 1._wp , t_su_1d(ji) - rt0         ) )  ! =1 if summer  
    8690               igravdr  = MAX( 0._wp , SIGN( 1._wp , t_bo_1d(ji) - t_su_1d(ji) ) )  ! =1 if t_su < t_bo 
    87                zsm_i_gd = - igravdr * MAX( sm_i_1d(ji) - rn_sal_gd , 0._wp ) / rn_time_gd * rdt_ice  ! gravity drainage  
    88                zsm_i_fl = - iflush  * MAX( sm_i_1d(ji) - rn_sal_fl , 0._wp ) / rn_time_fl * rdt_ice  ! flushing 
     91 
     92               zsm_i_gd = - igravdr * MAX( sm_i_1d(ji) - rn_sal_gd , 0._wp ) * z1_time_gd  ! gravity drainage  
     93               zsm_i_fl = - iflush  * MAX( sm_i_1d(ji) - rn_sal_fl , 0._wp ) * z1_time_fl  ! flushing 
    8994                
    9095               ! Update salinity    
  • branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

    r8534 r8559  
    4242   !!   ice_var_eqv2glo   : transform from VEQV to VGLO 
    4343   !!   ice_var_salprof   : salinity profile in the ice 
    44    !!   ice_var_bv        : brine volume 
    4544   !!   ice_var_salprof1d : salinity profile in the ice 1D 
    4645   !!   ice_var_zapsmall  : remove very small area and volume 
    4746   !!   ice_var_itd       : convert 1-cat to multiple cat 
     47   !!   ice_var_bv        : brine volume 
    4848   !!---------------------------------------------------------------------- 
    4949   USE dom_oce        ! ocean space and time domain 
     
    6464   PUBLIC   ice_var_eqv2glo       
    6565   PUBLIC   ice_var_salprof       
    66    PUBLIC   ice_var_bv            
    6766   PUBLIC   ice_var_salprof1d     
    6867   PUBLIC   ice_var_zapsmall 
    6968   PUBLIC   ice_var_itd 
     69   PUBLIC   ice_var_bv            
    7070 
    7171   !!---------------------------------------------------------------------- 
     
    170170      END WHERE 
    171171      ! 
     172      WHERE( v_i(:,:,:) > epsi20 )   ;   z1_v_i(:,:,:) = 1._wp / v_i(:,:,:) 
     173      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp 
     174      END WHERE 
     175      ! 
    172176      ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness 
    173177 
     
    177181         ht_i  (:,:,jpl) = zhmax 
    178182         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax  
    179          z1_a_i(:,:,jpl) = zhmax / v_i(:,:,jpl)          ! NB: v_i always /=0 as ht_i > hi_max 
     183         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)          ! NB: v_i always /=0 as ht_i > hi_max 
    180184      END WHERE 
    181185 
     
    185189 
    186190      IF( nn_icesal == 2 ) THEN                    !--- salinity (with a minimum value imposed everywhere) 
    187          WHERE( v_i(:,:,:) > epsi20 )   ;   sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) ) 
     191         WHERE( v_i(:,:,:) > epsi20 )   ;   sm_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, smv_i(:,:,:) * z1_v_i(:,:,:) ) ) 
    188192         ELSEWHERE                      ;   sm_i(:,:,:) = rn_simin 
    189193         END WHERE 
     
    202206                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area  
    203207                     ! 
    204                      ze_i             =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3] 
     208                     ze_i             =   e_i(ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3] 
    205209                     ztmelts          = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C] 
    206210                     ! Conversion q(S,T) -> T (second order equation) 
     
    295299         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) ) 
    296300         ! 
    297          DO jk = 1, nlay_i 
    298             s_i(:,:,jk,:)  = sm_i(:,:,:) 
     301         DO jl = 1, jpl 
     302            DO jk = 1, nlay_i 
     303               s_i(:,:,jk,jl)  = sm_i(:,:,jl) 
     304            END DO 
    299305         END DO 
    300306         !                                      ! Slope of the linear profile  
     
    353359   END SUBROUTINE ice_var_salprof 
    354360 
    355  
    356    SUBROUTINE ice_var_bv 
    357       !!------------------------------------------------------------------- 
    358       !!                ***  ROUTINE ice_var_bv *** 
    359       !! 
    360       !! ** Purpose :   computes mean brine volume (%) in sea ice 
    361       !! 
    362       !! ** Method  : e = - 0.054 * S (ppt) / T (C) 
    363       !! 
    364       !! References : Vancoppenolle et al., JGR, 2007 
    365       !!------------------------------------------------------------------- 
    366       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    367       !!------------------------------------------------------------------- 
    368       ! 
    369 !!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done 
    370 !!   instead of setting everything to zero as just below 
    371       bv_i (:,:,:) = 0._wp 
    372       DO jl = 1, jpl 
    373          DO jk = 1, nlay_i 
    374             WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )    
    375                bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 
    376             END WHERE 
    377          END DO 
    378       END DO 
    379       WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 
    380       ELSEWHERE                     ;   bvm_i(:,:) = 0._wp 
    381       END WHERE 
    382       ! 
    383    END SUBROUTINE ice_var_bv 
    384  
    385  
    386361   SUBROUTINE ice_var_salprof1d 
    387362      !!------------------------------------------------------------------- 
     
    393368      INTEGER  ::   ji, jk    ! dummy loop indices 
    394369      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars 
    395       REAL(wp) ::   zalpha, zs, zs0              !   -      - 
    396       ! 
    397       REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s   ! 
     370      REAL(wp) ::   zs, zs0              !   -      - 
     371      ! 
     372      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   ! 
    398373      REAL(wp), PARAMETER :: zsi0 = 3.5_wp 
    399374      REAL(wp), PARAMETER :: zsi1 = 4.5_wp 
     
    405380      CASE( 1 )       !  constant salinity in time and space  ! 
    406381         !            !---------------------------------------! 
    407          s_i_1d(:,:) = rn_icesal 
     382         s_i_1d(1:nidx,:) = rn_icesal 
    408383         ! 
    409384         !            !---------------------------------------------! 
     
    411386         !            !---------------------------------------------! 
    412387         ! 
    413          ALLOCATE( z_slope_s(jpij) ) 
    414          ! 
    415          DO ji = 1, nidx          ! Slope of the linear profile 
    416             rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 
    417             z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) ) 
    418          END DO 
    419  
     388         ALLOCATE( z_slope_s(jpij), zalpha(jpij) ) 
     389         ! 
     390         !                                      ! Slope of the linear profile  
     391         WHERE( ht_i_1d(1:nidx) > epsi20 )   ;   z_slope_s(1:nidx) = 2._wp * sm_i_1d(1:nidx) / ht_i_1d(1:nidx) 
     392         ELSEWHERE                           ;   z_slope_s(1:nidx) = 0._wp 
     393         END WHERE 
     394          
    420395         z1_dS = 1._wp / ( zsi1 - zsi0 ) 
     396         DO ji = 1, nidx 
     397            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp  )  ) 
     398            !                             ! force a constant profile when SSS too low (Baltic Sea) 
     399            IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp 
     400         END DO 
     401         ! 
     402         ! Computation of the profile 
    421403         DO jk = 1, nlay_i 
    422404            DO ji = 1, nidx 
    423                zalpha = MAX(  0._wp , MIN(  ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp  )  ) 
    424                IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) )   zalpha = 0._wp               ! cst profile when SSS too low (Baltic Sea) 
    425                ! 
    426                zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i   ! linear profile with 0 surface value 
    427                zs  = zalpha * zs0 + ( 1._wp - zalpha ) * sm_i_1d(ji)                      ! weighting the profile 
    428                s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )                     ! bounding salinity 
    429             END DO 
    430          END DO 
    431          ! 
    432          DEALLOCATE( z_slope_s ) 
     405               !                          ! linear profile with 0 surface value 
     406               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 
     407               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * sm_i_1d(ji) 
     408               s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) ) 
     409            END DO 
     410         END DO 
     411         ! 
     412         DEALLOCATE( z_slope_s, zalpha ) 
    433413 
    434414         !            !-------------------------------------------! 
     
    436416         !            !-------------------------------------------!                                   (mean = 2.30) 
    437417         ! 
    438          sm_i_1d(:) = 2.30_wp 
     418         sm_i_1d(1:nidx) = 2.30_wp 
    439419         ! 
    440420!!gm cf remark in ice_var_salprof routine, CASE( 3 ) 
     
    459439      !!------------------------------------------------------------------- 
    460440      INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices 
    461       REAL(wp) ::   zsal, zvi, zvs, zei, zes, zvp 
     441      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch 
    462442      !!------------------------------------------------------------------- 
    463443      ! 
     
    467447         ! Zap ice energy and use ocean heat to melt ice 
    468448         !----------------------------------------------------------------- 
     449         WHERE( a_i(:,:,jl) > epsi20 )   ;   ht_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl) 
     450         ELSEWHERE                       ;   ht_i(:,:,jl) = 0._wp 
     451         END WHERE 
     452         ! 
     453         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. ht_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp 
     454         ELSEWHERE                                                                            ;   zswitch(:,:) = 1._wp 
     455         END WHERE 
     456          
    469457         DO jk = 1, nlay_i 
    470458            DO jj = 1 , jpj 
    471459               DO ji = 1 , jpi 
    472 !!gm I think we can do better/faster  :  to be discussed 
    473                   rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
    474                   rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 
    475                   rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  & 
    476                      &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 
    477                   zei              = e_i(ji,jj,jk,jl) 
    478                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch 
    479                   t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch ) 
    480460                  ! update exchanges with ocean 
    481                   hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 
     461                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 
     462                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj) 
     463                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 
    482464               END DO 
    483465            END DO 
     
    486468         DO jj = 1 , jpj 
    487469            DO ji = 1 , jpi 
    488                rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 
    489                rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 
    490                rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  & 
    491                   &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 
    492                zsal = smv_i(ji,jj,  jl) 
    493                zvi  = v_i  (ji,jj,  jl) 
    494                zvs  = v_s  (ji,jj,  jl) 
    495                zes  = e_s  (ji,jj,1,jl) 
    496                IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw )   zvp  = v_ip (ji,jj  ,jl) 
     470               ! update exchanges with ocean 
     471               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * smv_i(ji,jj,jl)   * rhoic * r1_rdtice 
     472               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i  (ji,jj,jl)   * rhoic * r1_rdtice 
     473               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s  (ji,jj,jl)   * rhosn * r1_rdtice 
     474               hfx_res(ji,jj)  = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s  (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0 
    497475               !----------------------------------------------------------------- 
    498476               ! Zap snow energy  
    499477               !----------------------------------------------------------------- 
    500                t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch ) 
    501                e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch 
     478               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 
     479               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj) 
    502480 
    503481               !----------------------------------------------------------------- 
    504482               ! zap ice and snow volume, add water and salt to ocean 
    505483               !----------------------------------------------------------------- 
    506                ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj) 
    507                a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch 
    508                v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch 
    509                v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch 
    510                t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch ) 
    511                oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch 
    512                smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch 
    513  
    514                ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch 
    515                ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch 
    516  
    517                ! MV MP 2016 
    518                IF ( ln_pnd ) THEN  
    519                   a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch 
    520                   v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch 
    521                   IF ( ln_pnd_fw )   wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_ip(ji,jj,jl)  - zvp  ) * rhofw * r1_rdtice 
    522                ENDIF 
    523                ! END MV MP 2016 
    524  
    525                ! update exchanges with ocean 
    526                sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
    527                wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
    528                wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
    529                hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 
    530             END DO 
    531          END DO 
     484               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj) 
     485               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * zswitch(ji,jj) 
     486               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * zswitch(ji,jj) 
     487               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * zswitch(ji,jj) 
     488               t_su (ji,jj,jl) = t_su (ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) ) 
     489               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * zswitch(ji,jj) 
     490               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * zswitch(ji,jj) 
     491 
     492               ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * zswitch(ji,jj) 
     493               ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * zswitch(ji,jj) 
     494 
     495            END DO 
     496         END DO 
     497 
     498         IF( ln_pnd ) THEN  
     499            DO jj = 1 , jpj 
     500               DO ji = 1 , jpi 
     501                  IF( ln_pnd_fw )   & 
     502                     &   wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice 
     503                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
     504                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
     505               END DO 
     506            END DO 
     507         ENDIF 
     508          
    532509      END DO  
    533510 
     
    688665    END SUBROUTINE ice_var_itd 
    689666 
     667 
     668    SUBROUTINE ice_var_bv 
     669      !!------------------------------------------------------------------- 
     670      !!                ***  ROUTINE ice_var_bv *** 
     671      !! 
     672      !! ** Purpose :   computes mean brine volume (%) in sea ice 
     673      !! 
     674      !! ** Method  : e = - 0.054 * S (ppt) / T (C) 
     675      !! 
     676      !! References : Vancoppenolle et al., JGR, 2007 
     677      !!------------------------------------------------------------------- 
     678      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
     679      !!------------------------------------------------------------------- 
     680      ! 
     681!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done 
     682!!   instead of setting everything to zero as just below 
     683      bv_i (:,:,:) = 0._wp 
     684      DO jl = 1, jpl 
     685         DO jk = 1, nlay_i 
     686            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )    
     687               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 
     688            END WHERE 
     689         END DO 
     690      END DO 
     691      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 
     692      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp 
     693      END WHERE 
     694      ! 
     695   END SUBROUTINE ice_var_bv 
     696 
     697 
    690698#else 
    691699   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.