Ignore:
Timestamp:
2015-06-19T17:18:00+02:00 (5 years ago)
Author:
davestorkey
Message:

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r5234 r5443  
    1818   USE thd_ice          ! LIM thermodynamics 
    1919   USE ice              ! LIM variables 
    20    USE par_ice          ! LIM parameters 
    2120   USE dom_ice          ! LIM domain 
    22    USE limthd_lac       ! LIM 
    2321   USE limvar           ! LIM 
    24    USE in_out_manager   ! I/O manager 
    2522   USE lbclnk           ! lateral boundary condition - MPP exchanges 
    2623   USE lib_mpp          ! MPP library 
    2724   USE wrk_nemo         ! work arrays 
    2825   USE prtctl           ! Print control 
    29   ! Check budget (Rousset) 
     26 
     27   USE in_out_manager   ! I/O manager 
    3028   USE iom              ! I/O manager 
    3129   USE lib_fortran      ! glob_sum 
     
    4038   PUBLIC   lim_itd_me_icestrength 
    4139   PUBLIC   lim_itd_me_init 
    42    PUBLIC   lim_itd_me_zapsmall 
    43    PUBLIC   lim_itd_me_alloc        ! called by iceini.F90 
     40   PUBLIC   lim_itd_me_alloc        ! called by sbc_lim_init  
    4441 
    4542   !----------------------------------------------------------------------- 
     
    125122      !!  and Elizabeth C. Hunke, LANL are gratefully acknowledged 
    126123      !!--------------------------------------------------------------------! 
    127       INTEGER ::   ji, jj, jk, jl   ! dummy loop index 
    128       INTEGER ::   niter, nitermax = 20   ! local integer  
    129       LOGICAL  ::   asum_error            ! flag for asum .ne. 1 
     124      INTEGER  ::   ji, jj, jk, jl        ! dummy loop index 
     125      INTEGER  ::   niter                 ! local integer  
    130126      INTEGER  ::   iterate_ridging       ! if true, repeat the ridging 
    131       REAL(wp) ::   w1, tmpfac            ! local scalar 
     127      REAL(wp) ::   za, zfac              ! local scalar 
    132128      CHARACTER (len = 15) ::   fieldid 
    133       REAL(wp), POINTER, DIMENSION(:,:) ::   closing_net     ! net rate at which area is removed    (1/s) 
    134                                                              ! (ridging ice area - area of new ridges) / dt 
    135       REAL(wp), POINTER, DIMENSION(:,:) ::   divu_adv        ! divu as implied by transport scheme  (1/s) 
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   opning          ! rate of opening due to divergence/shear 
    137       REAL(wp), POINTER, DIMENSION(:,:) ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
    138       REAL(wp), POINTER, DIMENSION(:,:) ::   msnow_mlt       ! mass of snow added to ocean (kg m-2) 
    139       REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2) 
    140       REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
     129      REAL(wp), POINTER, DIMENSION(:,:)   ::   closing_net     ! net rate at which area is removed    (1/s) 
     130                                                               ! (ridging ice area - area of new ridges) / dt 
     131      REAL(wp), POINTER, DIMENSION(:,:)   ::   divu_adv        ! divu as implied by transport scheme  (1/s) 
     132      REAL(wp), POINTER, DIMENSION(:,:)   ::   opning          ! rate of opening due to divergence/shear 
     133      REAL(wp), POINTER, DIMENSION(:,:)   ::   closing_gross   ! rate at which area removed, not counting area of new ridges 
     134      REAL(wp), POINTER, DIMENSION(:,:)   ::   msnow_mlt       ! mass of snow added to ocean (kg m-2) 
     135      REAL(wp), POINTER, DIMENSION(:,:)   ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2) 
     136      REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
     137      ! 
     138      INTEGER, PARAMETER ::   nitermax = 20     
    141139      ! 
    142140      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
     
    144142      IF( nn_timing == 1 )  CALL timing_start('limitd_me') 
    145143 
    146       CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
     144      CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    147145 
    148146      IF(ln_ctl) THEN 
     
    156154      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    157155 
     156      CALL lim_var_zapsmall 
     157      CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
     158 
    158159      !-----------------------------------------------------------------------------! 
    159160      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    160161      !-----------------------------------------------------------------------------! 
    161       Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0                ! proport const for PE 
     162      Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0             ! proport const for PE 
    162163      ! 
    163164      CALL lim_itd_me_ridgeprep                                    ! prepare ridging 
     
    193194            !  (thick, newly ridged ice). 
    194195 
    195             closing_net(ji,jj) = Cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
     196            closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 
    196197 
    197198            ! 2.2 divu_adv 
     
    237238               ! Reduce the closing rate if more than 100% of the open water  
    238239               ! would be removed.  Reduce the opening rate proportionately. 
    239                IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 
    240                   w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 
    241                   IF ( w1 .GT. ato_i(ji,jj)) THEN 
    242                      tmpfac = ato_i(ji,jj) / w1 
    243                      closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 
    244                      opning(ji,jj) = opning(ji,jj) * tmpfac 
    245                   ENDIF !w1 
    246                ENDIF !at0i and athorn 
    247  
    248             END DO ! ji 
    249          END DO ! jj 
     240               za   = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 
     241               IF( za > epsi20 ) THEN 
     242                  zfac = MIN( 1._wp, ato_i(ji,jj) / za ) 
     243                  closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 
     244                  opning       (ji,jj) = opning       (ji,jj) * zfac 
     245               ENDIF 
     246 
     247            END DO 
     248         END DO 
    250249 
    251250         ! correction to closing rate / opening if excessive ice removal 
     
    253252         ! Reduce the closing rate if more than 100% of any ice category  
    254253         ! would be removed.  Reduce the opening rate proportionately. 
    255  
    256254         DO jl = 1, jpl 
    257255            DO jj = 1, jpj 
    258256               DO ji = 1, jpi 
    259                   IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN 
    260                      w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
    261                      IF ( w1  >  a_i(ji,jj,jl) ) THEN 
    262                         tmpfac = a_i(ji,jj,jl) / w1 
    263                         closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 
    264                         opning       (ji,jj) = opning       (ji,jj) * tmpfac 
    265                      ENDIF 
     257                  za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
     258                  IF( za  >  epsi20 ) THEN 
     259                     zfac = MIN( 1._wp, a_i(ji,jj,jl) / za ) 
     260                     closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 
     261                     opning       (ji,jj) = opning       (ji,jj) * zfac 
    266262                  ENDIF 
    267                END DO !ji 
    268             END DO ! jj 
    269          END DO !jl 
     263               END DO 
     264            END DO 
     265         END DO 
    270266 
    271267         ! 3.3 Redistribute area, volume, and energy. 
     
    276272         ! 3.4 Compute total area of ice plus open water after ridging. 
    277273         !-----------------------------------------------------------------------------! 
    278  
    279          CALL lim_itd_me_asumr 
     274         ! This is in general not equal to one because of divergence during transport 
     275         asum(:,:) = ato_i(:,:) 
     276         DO jl = 1, jpl 
     277            asum(:,:) = asum(:,:) + a_i(:,:,jl) 
     278         END DO 
    280279 
    281280         ! 3.5 Do we keep on iterating ??? 
     
    288287         DO jj = 1, jpj 
    289288            DO ji = 1, jpi 
    290                IF (ABS(asum(ji,jj) - kamax ) .LT. epsi10) THEN 
     289               IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN 
    291290                  closing_net(ji,jj) = 0._wp 
    292291                  opning     (ji,jj) = 0._wp 
     
    324323      ! Convert ridging rate diagnostics to correct units. 
    325324      ! Update fresh water and heat fluxes due to snow melt. 
    326  
    327       asum_error = .false.  
    328  
    329325      DO jj = 1, jpj 
    330326         DO ji = 1, jpi 
    331  
    332             IF(ABS(asum(ji,jj) - kamax) > epsi10 ) asum_error = .true. 
    333327 
    334328            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice 
     
    341335            !-----------------------------------------------------------------------------! 
    342336            wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
    343             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice  ! heat sink for ocean (<0, W.m-2) 
     337            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean (<0, W.m-2) 
    344338 
    345339         END DO 
     
    347341 
    348342      ! Check if there is a ridging error 
    349       DO jj = 1, jpj 
    350          DO ji = 1, jpi 
    351             IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug 
    352                WRITE(numout,*) ' ' 
    353                WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
    354                WRITE(numout,*) ' limitd_me ' 
    355                WRITE(numout,*) ' POINT : ', ji, jj 
    356                WRITE(numout,*) ' jpl, a_i, athorn ' 
    357                WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 
    358                DO jl = 1, jpl 
    359                   WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 
    360                END DO 
    361             ENDIF  ! asum 
    362  
    363          END DO !ji 
    364       END DO !jj 
     343      IF( lwp ) THEN 
     344         DO jj = 1, jpj 
     345            DO ji = 1, jpi 
     346               IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug 
     347                  WRITE(numout,*) ' ' 
     348                  WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
     349                  WRITE(numout,*) ' limitd_me ' 
     350                  WRITE(numout,*) ' POINT : ', ji, jj 
     351                  WRITE(numout,*) ' jpl, a_i, athorn ' 
     352                  WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 
     353                  DO jl = 1, jpl 
     354                     WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 
     355                  END DO 
     356               ENDIF 
     357            END DO 
     358         END DO 
     359      END IF 
    365360 
    366361      ! Conservation check 
     
    371366      ENDIF 
    372367 
     368      CALL lim_var_agg( 1 )  
     369 
    373370      !-----------------------------------------------------------------------------! 
    374       ! 6) Updating state variables and trend terms (done in limupdate) 
     371      ! control prints 
    375372      !-----------------------------------------------------------------------------! 
    376       CALL lim_var_glo2eqv 
    377       CALL lim_itd_me_zapsmall 
    378  
    379  
    380       IF(ln_ctl) THEN     ! Control print 
     373      IF(ln_ctl) THEN  
     374         CALL lim_var_glo2eqv 
     375 
    381376         CALL prt_ctl_info(' ') 
    382377         CALL prt_ctl_info(' - Cell values : ') 
    383378         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    384          CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_me  : cell area :') 
     379         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me  : cell area :') 
    385380         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me  : at_i      :') 
    386381         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me  : vt_i      :') 
     
    436431      !!---------------------------------------------------------------------- 
    437432      INTEGER, INTENT(in) ::   kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979) 
    438  
    439       INTEGER ::   ji,jj, jl   ! dummy loop indices 
    440       INTEGER ::   ksmooth     ! smoothing the resistance to deformation 
    441       INTEGER ::   numts_rm    ! number of time steps for the P smoothing 
    442       REAL(wp) ::   hi, zw1, zp, zdummy, zzc, z1_3   ! local scalars 
     433      INTEGER             ::   ji,jj, jl   ! dummy loop indices 
     434      INTEGER             ::   ksmooth     ! smoothing the resistance to deformation 
     435      INTEGER             ::   numts_rm    ! number of time steps for the P smoothing 
     436      REAL(wp)            ::   zhi, zp, z1_3  ! local scalars 
    443437      REAL(wp), POINTER, DIMENSION(:,:) ::   zworka   ! temporary array used here 
    444438      !!---------------------------------------------------------------------- 
     
    466460                  ! 
    467461                  IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN 
    468                      hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     462                     zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    469463                     !---------------------------- 
    470464                     ! PE loss from deforming ice 
    471465                     !---------------------------- 
    472                      strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * hi * hi 
     466                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi 
    473467 
    474468                     !-------------------------- 
    475469                     ! PE gain from rafting ice 
    476470                     !-------------------------- 
    477                      strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * hi * hi 
     471                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi 
    478472 
    479473                     !---------------------------- 
    480474                     ! PE gain from ridging ice 
    481475                     !---------------------------- 
    482                      strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl)/krdg(ji,jj,jl)     & 
    483                         * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )    
    484 !!gm Optimization:  (a**3-b**3)/(a-b) = a*a+ab+b*b   ==> less costly operations even if a**3 is replaced by a*a*a...                     
    485                   ENDIF            ! aicen > epsi10 
     476                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl)     & 
     477                        * z1_3 * ( hrmax(ji,jj,jl)**2 + hrmin(ji,jj,jl)**2 +  hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )   
     478                        !!(a**3-b**3)/(a-b) = a*a+ab+b*b                       
     479                  ENDIF 
    486480                  ! 
    487                END DO ! ji 
    488             END DO !jj 
    489          END DO !jl 
    490  
    491          zzc = Cf * Cp     ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and Cf accounts for frictional dissipation 
    492          strength(:,:) = zzc * strength(:,:) / aksum(:,:) 
    493  
     481               END DO 
     482            END DO 
     483         END DO 
     484    
     485         strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) 
     486                         ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 
    494487         ksmooth = 1 
    495488 
     
    499492      ELSE                      ! kstrngth ne 1:  Hibler (1979) form 
    500493         ! 
    501          strength(:,:) = Pstar * vt_i(:,:) * EXP( - C_rhg * ( 1._wp - at_i(:,:) )  ) 
     494         strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) )  ) 
    502495         ! 
    503496         ksmooth = 1 
     
    511504      ! CAN BE REMOVED 
    512505      ! 
    513       IF ( brinstren_swi == 1 ) THEN 
     506      IF( ln_icestr_bvf ) THEN 
    514507 
    515508         DO jj = 1, jpj 
    516509            DO ji = 1, jpi 
    517                IF ( bv_i(ji,jj) .GT. 0.0 ) THEN 
    518                   zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 ) 
    519                ELSE 
    520                   zdummy = 0.0 
    521                ENDIF 
    522510               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0))) 
    523             END DO              ! j 
    524          END DO                 ! i 
     511            END DO 
     512         END DO 
    525513 
    526514      ENDIF 
     
    538526         CALL lbc_lnk( strength, 'T', 1. ) 
    539527 
    540          DO jj = 2, jpj - 1 
    541             DO ji = 2, jpi - 1 
    542                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is 
    543                   ! present 
    544                   zworka(ji,jj) = 4.0 * strength(ji,jj)              & 
    545                      &          + strength(ji-1,jj) * tms(ji-1,jj) &   
    546                      &          + strength(ji+1,jj) * tms(ji+1,jj) &   
    547                      &          + strength(ji,jj-1) * tms(ji,jj-1) &   
    548                      &          + strength(ji,jj+1) * tms(ji,jj+1)     
    549  
    550                   zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1) 
    551                   zworka(ji,jj) = zworka(ji,jj) / zw1 
     528         DO jj = 2, jpjm1 
     529            DO ji = 2, jpim1 
     530               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN  
     531                  zworka(ji,jj) = ( 4.0 * strength(ji,jj)              & 
     532                     &                  + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) &   
     533                     &                  + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) & 
     534                     &            ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 
    552535               ELSE 
    553536                  zworka(ji,jj) = 0._wp 
     
    556539         END DO 
    557540 
    558          DO jj = 2, jpj - 1 
    559             DO ji = 2, jpi - 1 
     541         DO jj = 2, jpjm1 
     542            DO ji = 2, jpim1 
    560543               strength(ji,jj) = zworka(ji,jj) 
    561544            END DO 
     
    563546         CALL lbc_lnk( strength, 'T', 1. ) 
    564547 
    565       ENDIF ! ksmooth 
     548      ENDIF 
    566549 
    567550      !-------------------- 
     
    580563         DO jj = 1, jpj - 1 
    581564            DO ji = 1, jpi - 1 
    582                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN       ! ice is present 
     565               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN  
    583566                  numts_rm = 1 ! number of time steps for the running mean 
    584                   IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 
    585                   IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 
     567                  IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 
     568                  IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 
    586569                  zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm 
    587570                  strp2(ji,jj) = strp1(ji,jj) 
     
    612595      !!---------------------------------------------------------------------! 
    613596      INTEGER ::   ji,jj, jl    ! dummy loop indices 
    614       REAL(wp) ::   Gstari, astari, hi, hrmean, zdummy   ! local scalar 
     597      REAL(wp) ::   Gstari, astari, zhi, hrmean, zdummy   ! local scalar 
    615598      REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka    ! temporary array used here 
    616599      REAL(wp), POINTER, DIMENSION(:,:,:) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n 
     
    620603      CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 
    621604 
    622       Gstari     = 1.0/Gstar     
    623       astari     = 1.0/astar     
     605      Gstari     = 1.0/rn_gstar     
     606      astari     = 1.0/rn_astar     
    624607      aksum(:,:)    = 0.0 
    625608      athorn(:,:,:) = 0.0 
     
    632615 
    633616      !     ! Zero out categories with very small areas 
    634       CALL lim_itd_me_zapsmall 
     617      CALL lim_var_zapsmall 
    635618 
    636619      !------------------------------------------------------------------------------! 
     
    639622 
    640623      ! Compute total area of ice plus open water. 
    641       CALL lim_itd_me_asumr 
    642       ! This is in general not equal to one  
    643       ! because of divergence during transport 
     624      ! This is in general not equal to one because of divergence during transport 
     625      asum(:,:) = ato_i(:,:) 
     626      DO jl = 1, jpl 
     627         asum(:,:) = asum(:,:) + a_i(:,:,jl) 
     628      END DO 
    644629 
    645630      ! Compute cumulative thickness distribution function 
     
    649634 
    650635      Gsum(:,:,-1) = 0._wp 
    651  
    652       DO jj = 1, jpj 
    653          DO ji = 1, jpi 
    654             IF( ato_i(ji,jj) > epsi10 ) THEN   ;   Gsum(ji,jj,0) = ato_i(ji,jj) 
    655             ELSE                               ;   Gsum(ji,jj,0) = 0._wp 
    656             ENDIF 
    657          END DO 
    658       END DO 
     636      Gsum(:,:,0 ) = ato_i(:,:) 
    659637 
    660638      ! for each value of h, you have to add ice concentration then 
    661639      DO jl = 1, jpl 
    662          DO jj = 1, jpj  
    663             DO ji = 1, jpi 
    664                IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 
    665                ELSE                                   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 
    666                ENDIF 
    667             END DO 
    668          END DO 
     640         Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 
    669641      END DO 
    670642 
     
    687659      !----------------------------------------------------------------- 
    688660 
    689       IF( partfun_swi == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975) 
     661      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975) 
    690662         DO jl = 0, jpl     
    691663            DO jj = 1, jpj  
    692664               DO ji = 1, jpi 
    693                   IF( Gsum(ji,jj,jl) < Gstar) THEN 
    694                      athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * & 
    695                         (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari) 
    696                   ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN 
    697                      athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) *  & 
    698                         (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari) 
     665                  IF( Gsum(ji,jj,jl) < rn_gstar) THEN 
     666                     athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * & 
     667                        &                        ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari ) 
     668                  ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN 
     669                     athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) *  & 
     670                        &                        ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari ) 
    699671                  ELSE 
    700672                     athorn(ji,jj,jl) = 0.0 
    701673                  ENDIF 
    702                END DO ! ji 
    703             END DO ! jj 
    704          END DO ! jl  
     674               END DO 
     675            END DO 
     676         END DO 
    705677 
    706678      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007) 
    707679         !                         
    708680         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array 
    709  
    710681         DO jl = -1, jpl 
    711682            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 
    712          END DO !jl 
     683         END DO 
    713684         DO jl = 0, jpl 
    714685             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 
    715686         END DO 
    716687         ! 
    717       ENDIF ! partfun_swi 
    718  
    719       IF( raft_swi == 1 ) THEN      ! Ridging and rafting ice participation functions 
     688      ENDIF 
     689 
     690      IF( ln_rafting ) THEN      ! Ridging and rafting ice participation functions 
    720691         ! 
    721692         DO jl = 1, jpl 
    722693            DO jj = 1, jpj  
    723694               DO ji = 1, jpi 
    724                   IF ( athorn(ji,jj,jl) .GT. 0._wp ) THEN 
     695                  IF ( athorn(ji,jj,jl) > 0._wp ) THEN 
    725696!!gm  TANH( -X ) = - TANH( X )  so can be computed only 1 time.... 
    726                      aridge(ji,jj,jl) = ( TANH (  Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 
    727                      araft (ji,jj,jl) = ( TANH ( -Craft * ( ht_i(ji,jj,jl) - hparmeter ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 
     697                     aridge(ji,jj,jl) = ( TANH (  rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 
     698                     araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 
    728699                     IF ( araft(ji,jj,jl) < epsi06 )   araft(ji,jj,jl)  = 0._wp 
    729700                     aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 ) 
    730                   ENDIF ! athorn 
    731                END DO ! ji 
    732             END DO ! jj 
    733          END DO ! jl 
    734  
    735       ELSE  ! raft_swi = 0 
     701                  ENDIF 
     702               END DO 
     703            END DO 
     704         END DO 
     705 
     706      ELSE 
    736707         ! 
    737708         DO jl = 1, jpl 
     
    741712      ENDIF 
    742713 
    743       IF ( raft_swi == 1 ) THEN 
    744  
    745          IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN 
     714      IF( ln_rafting ) THEN 
     715 
     716         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN 
    746717            DO jl = 1, jpl 
    747718               DO jj = 1, jpj 
    748719                  DO ji = 1, jpi 
    749                      IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN 
     720                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN 
    750721                        WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 
    751722                        WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl 
     
    793764            DO ji = 1, jpi 
    794765 
    795                IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 
    796                   hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    797                   hrmean          = MAX(SQRT(Hstar*hi), hi*krdgmin) 
    798                   hrmin(ji,jj,jl) = MIN(2.0*hi, 0.5*(hrmean + hi)) 
     766               IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN 
     767                  zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     768                  hrmean          = MAX(SQRT(rn_hstar*zhi), zhi*krdgmin) 
     769                  hrmin(ji,jj,jl) = MIN(2.0*zhi, 0.5*(hrmean + zhi)) 
    799770                  hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl) 
    800                   hraft(ji,jj,jl) = kraft*hi 
    801                   krdg(ji,jj,jl)  = hrmean / hi 
     771                  hraft(ji,jj,jl) = kraft*zhi 
     772                  krdg(ji,jj,jl)  = hrmean / zhi 
    802773               ELSE 
    803774                  hraft(ji,jj,jl) = 0.0 
     
    807778               ENDIF 
    808779 
    809             END DO ! ji 
    810          END DO ! jj 
    811       END DO ! jl 
     780            END DO 
     781         END DO 
     782      END DO 
    812783 
    813784      ! Normalization factor : aksum, ensures mass conservation 
     
    841812      LOGICAL, PARAMETER ::   l_conservation_check = .true.  ! if true, check conservation (useful for debugging) 
    842813      ! 
    843       LOGICAL ::   neg_ato_i      ! flag for ato_i(i,j) < -puny 
    844       LOGICAL ::   large_afrac    ! flag for afrac > 1 
    845       LOGICAL ::   large_afrft    ! flag for afrac > 1 
    846814      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices 
    847815      INTEGER ::   ij                ! horizontal index, combines i and j loops 
    848816      INTEGER ::   icells            ! number of cells with aicen > puny 
    849       REAL(wp) ::   hL, hR, farea, zdummy, zdummy0, ztmelts    ! left and right limits of integration 
    850       REAL(wp) ::   zsstK            ! SST in Kelvin 
     817      REAL(wp) ::   hL, hR, farea, ztmelts    ! left and right limits of integration 
    851818 
    852819      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices 
     
    864831      REAL(wp), POINTER, DIMENSION(:,:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges 
    865832      REAL(wp), POINTER, DIMENSION(:,:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice 
    866       REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! areal age content of ridged & rifging ice 
    867833      REAL(wp), POINTER, DIMENSION(:,:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2 
    868834 
     
    873839      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg2   ! sal*volume of new ridges 
    874840      REAL(wp), POINTER, DIMENSION(:,:) ::   smsw    ! sal*volume of water trapped into ridges 
     841      REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! ice age of ice ridged 
    875842 
    876843      REAL(wp), POINTER, DIMENSION(:,:) ::   afrft            ! fraction of category area rafted 
     
    878845      REAL(wp), POINTER, DIMENSION(:,:) ::   virft , vsrft    ! ice & snow volume of rafting ice 
    879846      REAL(wp), POINTER, DIMENSION(:,:) ::   esrft , smrft    ! snow energy & salinity of rafting ice 
    880       REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! areal age content of rafted ice & rafting ice 
     847      REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! ice age of ice rafted 
    881848 
    882849      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eirft      ! ice energy of rafting ice 
     
    886853      !!---------------------------------------------------------------------- 
    887854 
    888       CALL wrk_alloc( (jpi+1)*(jpj+1),      indxi, indxj ) 
    889       CALL wrk_alloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final ) 
    890       CALL wrk_alloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 ) 
    891       CALL wrk_alloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw ) 
    892       CALL wrk_alloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    893       CALL wrk_alloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
    894       CALL wrk_alloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw ) 
    895       CALL wrk_alloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 
     855      CALL wrk_alloc( (jpi+1)*(jpj+1),       indxi, indxj ) 
     856      CALL wrk_alloc( jpi, jpj,              vice_init, vice_final, eice_init, eice_final ) 
     857      CALL wrk_alloc( jpi, jpj,              afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 
     858      CALL wrk_alloc( jpi, jpj,              vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 
     859      CALL wrk_alloc( jpi, jpj,              afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
     860      CALL wrk_alloc( jpi, jpj, jpl,         aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
     861      CALL wrk_alloc( jpi, jpj, nlay_i,      eirft, erdg1, erdg2, ersw ) 
     862      CALL wrk_alloc( jpi, jpj, nlay_i, jpl, eicen_init ) 
    896863 
    897864      ! Conservation check 
     
    901868         CALL lim_column_sum        (jpl,    v_i,       vice_init ) 
    902869         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_init ) 
    903          DO ji = mi0(jiindx), mi1(jiindx) 
    904             DO jj = mj0(jjindx), mj1(jjindx) 
     870         DO ji = mi0(iiceprt), mi1(iiceprt) 
     871            DO jj = mj0(jiceprt), mj1(jiceprt) 
    905872               WRITE(numout,*) ' vice_init  : ', vice_init(ji,jj) 
    906873               WRITE(numout,*) ' eice_init  : ', eice_init(ji,jj) 
     
    912879      ! 1) Compute change in open water area due to closing and opening. 
    913880      !------------------------------------------------------------------------------- 
    914  
    915       neg_ato_i = .false. 
    916  
    917881      DO jj = 1, jpj 
    918882         DO ji = 1, jpi 
    919883            ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        & 
    920884               &                        + opning(ji,jj)                          * rdt_ice 
    921             IF( ato_i(ji,jj) < -epsi10 ) THEN 
    922                neg_ato_i = .TRUE. 
    923             ELSEIF( ato_i(ji,jj) < 0._wp ) THEN    ! roundoff error 
     885            IF    ( ato_i(ji,jj) < -epsi10 ) THEN    ! there is a bug 
     886               IF(lwp)   WRITE(numout,*) 'Ridging error: ato_i < 0 -- ato_i : ',ato_i(ji,jj) 
     887            ELSEIF( ato_i(ji,jj) < 0._wp   ) THEN    ! roundoff error 
    924888               ato_i(ji,jj) = 0._wp 
    925889            ENDIF 
    926          END DO !jj 
    927       END DO !ji 
    928  
    929       ! if negative open water area alert it 
    930       IF( neg_ato_i ) THEN       ! there is a bug 
    931          DO jj = 1, jpj  
    932             DO ji = 1, jpi 
    933                IF( ato_i(ji,jj) < -epsi10 ) THEN  
    934                   WRITE(numout,*) ''   
    935                   WRITE(numout,*) 'Ridging error: ato_i < 0' 
    936                   WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) 
    937                ENDIF               ! ato_i < -epsi10 
    938             END DO 
    939          END DO 
    940       ENDIF 
     890         END DO 
     891      END DO 
    941892 
    942893      !----------------------------------------------------------------- 
    943894      ! 2) Save initial state variables 
    944895      !----------------------------------------------------------------- 
    945  
    946       DO jl = 1, jpl 
    947          aicen_init(:,:,jl) = a_i(:,:,jl) 
    948          vicen_init(:,:,jl) = v_i(:,:,jl) 
    949          vsnwn_init(:,:,jl) = v_s(:,:,jl) 
    950          ! 
    951          smv_i_init(:,:,jl) = smv_i(:,:,jl) 
    952          oa_i_init (:,:,jl) = oa_i (:,:,jl) 
    953       END DO !jl 
    954  
    955       esnwn_init(:,:,:) = e_s(:,:,1,:) 
    956  
    957       DO jl = 1, jpl   
    958          DO jk = 1, nlay_i 
    959             eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl) 
    960          END DO 
    961       END DO 
     896      aicen_init(:,:,:)   = a_i  (:,:,:) 
     897      vicen_init(:,:,:)   = v_i  (:,:,:) 
     898      vsnwn_init(:,:,:)   = v_s  (:,:,:) 
     899      smv_i_init(:,:,:)   = smv_i(:,:,:) 
     900      esnwn_init(:,:,:)   = e_s  (:,:,1,:) 
     901      eicen_init(:,:,:,:) = e_i  (:,:,:,:) 
     902      oa_i_init (:,:,:)   = oa_i (:,:,:) 
    962903 
    963904      ! 
     
    982923                  indxi(icells) = ji 
    983924                  indxj(icells) = jj 
    984                ENDIF ! test on a_icen_init  
    985             END DO ! ji 
    986          END DO ! jj 
    987  
    988          large_afrac = .false. 
    989          large_afrft = .false. 
    990  
    991 !CDIR NODEP 
     925               ENDIF 
     926            END DO 
     927         END DO 
     928 
    992929         DO ij = 1, icells 
    993930            ji = indxi(ij) 
     
    1003940            arft2(ji,jj) = arft1(ji,jj) / kraft 
    1004941 
    1005             oirdg1(ji,jj)= aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice 
    1006             oirft1(ji,jj)= araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice 
    1007             oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1) 
    1008             oirft2(ji,jj)= oirft1(ji,jj) / kraft 
    1009  
    1010942            !--------------------------------------------------------------- 
    1011943            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1  
     
    1015947            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 
    1016948 
    1017             IF (afrac(ji,jj) > kamax + epsi10) THEN  !riging 
    1018                large_afrac = .true. 
    1019             ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error 
     949            IF( afrac(ji,jj) > kamax + epsi10 ) THEN  ! there is a bug 
     950               IF(lwp)   WRITE(numout,*) ' ardg > a_i -- ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1) 
     951            ELSEIF( afrac(ji,jj) > kamax ) THEN       ! roundoff error 
    1020952               afrac(ji,jj) = kamax 
    1021953            ENDIF 
    1022             IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting 
    1023                large_afrft = .true. 
    1024             ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error 
     954 
     955            IF( afrft(ji,jj) > kamax + epsi10 ) THEN ! there is a bug 
     956               IF(lwp)   WRITE(numout,*) ' arft > a_i -- arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1)  
     957            ELSEIF( afrft(ji,jj) > kamax) THEN       ! roundoff error 
    1025958               afrft(ji,jj) = kamax 
    1026959            ENDIF 
     
    1031964            !-------------------------------------------------------------------------- 
    1032965            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) 
    1033             vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por ) 
    1034             vsw  (ji,jj) = vrdg1(ji,jj) * ridge_por 
    1035  
    1036             vsrdg(ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 
    1037             esrdg(ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 
    1038             srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 
    1039             srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) !! MV HC 2014 this line seems useless 
     966            vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg ) 
     967            vsw  (ji,jj) = vrdg1(ji,jj) * rn_por_rdg 
     968 
     969            vsrdg (ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 
     970            esrdg (ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 
     971            srdg1 (ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 
     972            oirdg1(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) 
     973            oirdg2(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) / krdg(ji,jj,jl1)  
    1040974 
    1041975            ! rafting volumes, heat contents ... 
    1042             virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 
    1043             vsrft(ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 
    1044             esrft(ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 
    1045             smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)  
     976            virft (ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 
     977            vsrft (ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 
     978            esrft (ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 
     979            smrft (ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj)  
     980            oirft1(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj)  
     981            oirft2(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) / kraft  
    1046982 
    1047983            ! substract everything 
    1048             a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1(ji,jj)  - arft1(ji,jj) 
    1049             v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1(ji,jj)  - virft(ji,jj) 
    1050             v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg(ji,jj)  - vsrft(ji,jj) 
    1051             e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj)  - esrft(ji,jj) 
     984            a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1 (ji,jj) - arft1 (ji,jj) 
     985            v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1 (ji,jj) - virft (ji,jj) 
     986            v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg (ji,jj) - vsrft (ji,jj) 
     987            e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg (ji,jj) - esrft (ji,jj) 
     988            smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1 (ji,jj) - smrft (ji,jj) 
    1052989            oa_i(ji,jj,jl1)  = oa_i(ji,jj,jl1)  - oirdg1(ji,jj) - oirft1(ji,jj) 
    1053             smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj)  - smrft(ji,jj) 
    1054990 
    1055991            !----------------------------------------------------------------- 
    1056992            ! 3.5) Compute properties of new ridges 
    1057993            !----------------------------------------------------------------- 
    1058             !------------- 
     994            !--------- 
    1059995            ! Salinity 
    1060             !------------- 
     996            !--------- 
    1061997            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014 
    1062998            srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge 
    1063999 
    1064             !srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity 
     1000            !srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity 
    10651001             
    10661002            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 
    1067             wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice   ! gurvan: increase in ice volume du to seawater frozen in voids              
     1003            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice   ! increase in ice volume du to seawater frozen in voids              
    10681004 
    10691005            !------------------------------------             
     
    10911027            !           ij looping 1-icells 
    10921028 
    1093             msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg)   &   ! rafting included 
    1094                &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 
    1095  
    1096             ! in 1e-9 Joules (same as e_s) 
    1097             esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included 
    1098                &                                - esrft(ji,jj)*(1.0-fsnowrft)           
     1029            msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-rn_fsnowrdg)   &   ! rafting included 
     1030               &                                + rhosn*vsrft(ji,jj)*(1.0-rn_fsnowrft) 
     1031 
     1032            ! in J/m2 (same as e_s) 
     1033            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-rn_fsnowrdg)         &   !rafting included 
     1034               &                                - esrft(ji,jj)*(1.0-rn_fsnowrft)           
    10991035 
    11001036            !----------------------------------------------------------------- 
     
    11091045            dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 
    11101046 
    1111          END DO                 ! ij 
     1047         END DO 
    11121048 
    11131049         !-------------------------------------------------------------------- 
     
    11161052         !-------------------------------------------------------------------- 
    11171053         DO jk = 1, nlay_i 
    1118 !CDIR NODEP 
    11191054            DO ij = 1, icells 
    11201055               ji = indxi(ij) 
     
    11281063               ! enthalpy of the trapped seawater (J/m2, >0) 
    11291064               ! clem: if sst>0, then ersw <0 (is that possible?) 
    1130                zsstK  = sst_m(ji,jj) + rt0 
    1131                ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i ) 
     1065               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i 
    11321066 
    11331067               ! heat flux to the ocean 
    11341068               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux  
    11351069 
    1136                ! Correct dimensions to avoid big values 
    1137                ersw(ji,jj,jk)   = ersw(ji,jj,jk) / unit_fac 
    1138  
    1139                ! Mutliply by ice volume, and divide by number of layers to get heat content in 1.e9 J 
    1140                ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean  
    1141                !! MV HC 2014 
    1142                ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) 
    1143  
     1070               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 
    11441071               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
    11451072 
    1146             END DO ! ij 
    1147          END DO !jk 
     1073            END DO 
     1074         END DO 
    11481075 
    11491076 
    11501077         IF( con_i ) THEN 
    11511078            DO jk = 1, nlay_i 
    1152 !CDIR NODEP 
    11531079               DO ij = 1, icells 
    11541080                  ji = indxi(ij) 
    11551081                  jj = indxj(ij) 
    11561082                  eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk) 
    1157                END DO ! ij 
    1158             END DO !jk 
    1159          ENDIF 
    1160  
    1161          IF( large_afrac ) THEN   ! there is a bug 
    1162 !CDIR NODEP 
    1163             DO ij = 1, icells 
    1164                ji = indxi(ij) 
    1165                jj = indxj(ij) 
    1166                IF( afrac(ji,jj) > kamax + epsi10 ) THEN  
    1167                   WRITE(numout,*) '' 
    1168                   WRITE(numout,*) ' ardg > a_i' 
    1169                   WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1) 
    1170                ENDIF 
    1171             END DO 
    1172          ENDIF 
    1173          IF( large_afrft ) THEN  ! there is a bug 
    1174 !CDIR NODEP 
    1175             DO ij = 1, icells 
    1176                ji = indxi(ij) 
    1177                jj = indxj(ij) 
    1178                IF( afrft(ji,jj) > kamax + epsi10 ) THEN  
    1179                   WRITE(numout,*) '' 
    1180                   WRITE(numout,*) ' arft > a_i' 
    1181                   WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1) 
    1182                ENDIF 
     1083               END DO 
    11831084            END DO 
    11841085         ENDIF 
     
    11901091         DO jl2  = 1, jpl  
    11911092            ! over categories to which ridged ice is transferred 
    1192 !CDIR NODEP 
    11931093            DO ij = 1, icells 
    11941094               ji = indxi(ij) 
     
    11991099               ! Transfer area, volume, and energy accordingly. 
    12001100 
    1201                IF( hrmin(ji,jj,jl1) >= hi_max(jl2)  .OR.        & 
    1202                    hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 
     1101               IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 
    12031102                  hL = 0._wp 
    12041103                  hR = 0._wp 
     
    12141113               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ardg2 (ji,jj) * farea 
    12151114               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj) 
    1216                v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg 
    1217                e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg 
     1115               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg 
     1116               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg 
    12181117               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + srdg2 (ji,jj) * fvol(ji,jj) 
    12191118               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea 
    12201119 
    1221             END DO ! ij 
     1120            END DO 
    12221121 
    12231122            ! Transfer ice energy to category jl2 by ridging 
    12241123            DO jk = 1, nlay_i 
    1225 !CDIR NODEP 
    12261124               DO ij = 1, icells 
    12271125                  ji = indxi(ij) 
    12281126                  jj = indxj(ij) 
    1229                   e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj)*erdg2(ji,jj,jk) 
     1127                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk) 
    12301128               END DO 
    12311129            END DO 
     
    12351133         DO jl2 = 1, jpl  
    12361134 
    1237 !CDIR NODEP 
    12381135            DO ij = 1, icells 
    12391136               ji = indxi(ij) 
     
    12421139               ! thickness category jl2, transfer area, volume, and energy accordingly. 
    12431140               ! 
    1244                IF( hraft(ji,jj,jl1) <= hi_max(jl2)  .AND.        & 
    1245                    hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
     1141               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN 
    12461142                  a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj) 
    12471143                  v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj) 
    1248                   v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * fsnowrft 
    1249                   e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft 
     1144                  v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * rn_fsnowrft 
     1145                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft 
    12501146                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)     
    1251                   oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)     
    1252                ENDIF ! hraft 
     1147                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj) 
     1148               ENDIF 
    12531149               ! 
    1254             END DO ! ij 
     1150            END DO 
    12551151 
    12561152            ! Transfer rafted ice energy to category jl2  
    12571153            DO jk = 1, nlay_i 
    1258 !CDIR NODEP 
    12591154               DO ij = 1, icells 
    12601155                  ji = indxi(ij) 
    12611156                  jj = indxj(ij) 
    1262                   IF(  hraft(ji,jj,jl1)  <=  hi_max(jl2)   .AND.        & 
    1263                        hraft(ji,jj,jl1)  >   hi_max(jl2-1)  ) THEN 
     1157                  IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1)  ) THEN 
    12641158                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 
    12651159                  ENDIF 
    1266                END DO           ! ij 
    1267             END DO !jk 
    1268  
    1269          END DO ! jl2 
     1160               END DO 
     1161            END DO 
     1162 
     1163         END DO 
    12701164 
    12711165      END DO ! jl1 (deforming categories) 
     
    12811175         CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid)  
    12821176 
    1283          DO ji = mi0(jiindx), mi1(jiindx) 
    1284             DO jj = mj0(jjindx), mj1(jjindx) 
     1177         DO ji = mi0(iiceprt), mi1(iiceprt) 
     1178            DO jj = mj0(jiceprt), mj1(jiceprt) 
    12851179               WRITE(numout,*) ' vice_init  : ', vice_init (ji,jj) 
    12861180               WRITE(numout,*) ' vice_final : ', vice_final(ji,jj) 
     
    12911185      ENDIF 
    12921186      ! 
    1293       CALL wrk_dealloc( (jpi+1)*(jpj+1),      indxi, indxj ) 
    1294       CALL wrk_dealloc( jpi, jpj,             vice_init, vice_final, eice_init, eice_final ) 
    1295       CALL wrk_dealloc( jpi, jpj,             afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 ) 
    1296       CALL wrk_dealloc( jpi, jpj,             vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw ) 
    1297       CALL wrk_dealloc( jpi, jpj,             afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
    1298       CALL wrk_dealloc( jpi, jpj, jpl,        aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
    1299       CALL wrk_dealloc( jpi, jpj, nlay_i+1,      eirft, erdg1, erdg2, ersw ) 
    1300       CALL wrk_dealloc( jpi, jpj, nlay_i+1, jpl, eicen_init ) 
     1187      CALL wrk_dealloc( (jpi+1)*(jpj+1),        indxi, indxj ) 
     1188      CALL wrk_dealloc( jpi, jpj,               vice_init, vice_final, eice_init, eice_final ) 
     1189      CALL wrk_dealloc( jpi, jpj,               afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 
     1190      CALL wrk_dealloc( jpi, jpj,               vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 
     1191      CALL wrk_dealloc( jpi, jpj,               afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 
     1192      CALL wrk_dealloc( jpi, jpj, jpl,          aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 
     1193      CALL wrk_dealloc( jpi, jpj, nlay_i,       eirft, erdg1, erdg2, ersw ) 
     1194      CALL wrk_dealloc( jpi, jpj, nlay_i, jpl, eicen_init ) 
    13011195      ! 
    13021196   END SUBROUTINE lim_itd_me_ridgeshift 
    1303  
    1304  
    1305    SUBROUTINE lim_itd_me_asumr 
    1306       !!----------------------------------------------------------------------------- 
    1307       !!                ***  ROUTINE lim_itd_me_asumr *** 
    1308       !! 
    1309       !! ** Purpose :   finds total fractional area 
    1310       !! 
    1311       !! ** Method  :   Find the total area of ice plus open water in each grid cell. 
    1312       !!              This is similar to the aggregate_area subroutine except that the 
    1313       !!              total area can be greater than 1, so the open water area is  
    1314       !!              included in the sum instead of being computed as a residual.  
    1315       !!----------------------------------------------------------------------------- 
    1316       INTEGER ::   jl   ! dummy loop index 
    1317       !!----------------------------------------------------------------------------- 
    1318       ! 
    1319       asum(:,:) = ato_i(:,:)                    ! open water 
    1320       DO jl = 1, jpl                            ! ice categories 
    1321          asum(:,:) = asum(:,:) + a_i(:,:,jl) 
    1322       END DO 
    1323       ! 
    1324    END SUBROUTINE lim_itd_me_asumr 
    1325  
    13261197 
    13271198   SUBROUTINE lim_itd_me_init 
     
    13391210      !!------------------------------------------------------------------- 
    13401211      INTEGER :: ios                 ! Local integer output status for namelist read 
    1341       NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,              &  
    1342         &                   Gstar, astar, Hstar, raft_swi, hparmeter, Craft, ridge_por, & 
    1343         &                   partfun_swi, brinstren_swi 
     1212      NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft,              &  
     1213        &                   rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, & 
     1214        &                   nn_partfun 
    13441215      !!------------------------------------------------------------------- 
    13451216      ! 
     
    13571228         WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution ' 
    13581229         WRITE(numout,*)' ~~~~~~~~~~~~~~~' 
    1359          WRITE(numout,*)'   Switch choosing the ice redistribution scheme           ridge_scheme_swi', ridge_scheme_swi  
    1360          WRITE(numout,*)'   Fraction of shear energy contributing to ridging        Cs              ', Cs  
    1361          WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging      Cf              ', Cf  
    1362          WRITE(numout,*)'   Fraction of snow volume conserved during ridging        fsnowrdg        ', fsnowrdg  
    1363          WRITE(numout,*)'   Fraction of snow volume conserved during ridging        fsnowrft        ', fsnowrft  
    1364          WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  Gstar           ', Gstar 
    1365          WRITE(numout,*)'   Equivalent to G* for an exponential part function       astar           ', astar 
    1366          WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     Hstar           ', Hstar 
    1367          WRITE(numout,*)'   Rafting of ice sheets or not                            raft_swi        ', raft_swi 
    1368          WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       hparmeter       ', hparmeter 
    1369          WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  Craft           ', Craft   
    1370          WRITE(numout,*)'   Initial porosity of ridges                              ridge_por       ', ridge_por 
    1371          WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    partfun_swi     ', partfun_swi 
    1372          WRITE(numout,*)'   Switch for including brine volume in ice strength comp. brinstren_swi   ', brinstren_swi 
     1230         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs  
     1231         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg  
     1232         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft  
     1233         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar 
     1234         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar 
     1235         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar 
     1236         WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting 
     1237         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft 
     1238         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft   
     1239         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg 
     1240         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun 
    13731241      ENDIF 
    13741242      ! 
    13751243   END SUBROUTINE lim_itd_me_init 
    1376  
    1377  
    1378    SUBROUTINE lim_itd_me_zapsmall 
    1379       !!------------------------------------------------------------------- 
    1380       !!                   ***  ROUTINE lim_itd_me_zapsmall *** 
    1381       !! 
    1382       !! ** Purpose :   Remove too small sea ice areas and correct salt fluxes 
    1383       !! 
    1384       !! history : 
    1385       !! author: William H. Lipscomb, LANL 
    1386       !! Nov 2003:  Modified by Julie Schramm to conserve volume and energy 
    1387       !! Sept 2004: Modified by William Lipscomb; replaced normalize_state with 
    1388       !!            additions to local freshwater, salt, and heat fluxes 
    1389       !!  9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code 
    1390       !!------------------------------------------------------------------- 
    1391       INTEGER ::   ji, jj, jl, jk   ! dummy loop indices 
    1392       INTEGER ::   icells           ! number of cells with ice to zap 
    1393  
    1394       REAL(wp), POINTER, DIMENSION(:,:) ::   zmask   ! 2D workspace 
    1395       REAL(wp)                          ::   zmask_glo, zsal, zvi, zvs, zei, zes 
    1396 !!gm      REAL(wp) ::   xtmp      ! temporary variable 
    1397       !!------------------------------------------------------------------- 
    1398  
    1399       CALL wrk_alloc( jpi, jpj, zmask ) 
    1400  
    1401       ! to be sure that at_i is the sum of a_i(jl) 
    1402       at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    1403  
    1404       DO jl = 1, jpl 
    1405          !----------------------------------------------------------------- 
    1406          ! Count categories to be zapped. 
    1407          !----------------------------------------------------------------- 
    1408          icells = 0 
    1409          zmask(:,:)  = 0._wp 
    1410          DO jj = 1, jpj 
    1411             DO ji = 1, jpi 
    1412                IF( a_i(ji,jj,jl) <= epsi10 .OR. v_i(ji,jj,jl) <= epsi10 .OR. at_i(ji,jj) <= epsi10 ) THEN 
    1413                   zmask(ji,jj) = 1._wp 
    1414                ENDIF 
    1415             END DO 
    1416          END DO 
    1417          !zmask_glo = glob_sum(zmask) 
    1418          !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean ' 
    1419  
    1420          !----------------------------------------------------------------- 
    1421          ! Zap ice energy and use ocean heat to melt ice 
    1422          !----------------------------------------------------------------- 
    1423  
    1424          DO jk = 1, nlay_i 
    1425             DO jj = 1 , jpj 
    1426                DO ji = 1 , jpi 
    1427                   zei  = e_i(ji,jj,jk,jl) 
    1428                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) 
    1429                   t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) + rtt * zmask(ji,jj) 
    1430                   ! update exchanges with ocean 
    1431                   hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    1432                END DO 
    1433             END DO 
    1434          END DO 
    1435  
    1436          DO jj = 1 , jpj 
    1437             DO ji = 1 , jpi 
    1438                 
    1439                zsal = smv_i(ji,jj,jl) 
    1440                zvi  = v_i(ji,jj,jl) 
    1441                zvs  = v_s(ji,jj,jl) 
    1442                zes  = e_s(ji,jj,1,jl) 
    1443                !----------------------------------------------------------------- 
    1444                ! Zap snow energy and use ocean heat to melt snow 
    1445                !----------------------------------------------------------------- 
    1446                !           xtmp = esnon(i,j,n) / dt ! < 0 
    1447                !           fhnet(i,j)      = fhnet(i,j)      + xtmp 
    1448                !           fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp 
    1449                ! xtmp is greater than 0 
    1450                ! fluxes are positive to the ocean 
    1451                ! here the flux has to be negative for the ocean 
    1452                t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 
    1453  
    1454                !----------------------------------------------------------------- 
    1455                ! zap ice and snow volume, add water and salt to ocean 
    1456                !----------------------------------------------------------------- 
    1457                ato_i(ji,jj)    = a_i  (ji,jj,jl) *           zmask(ji,jj)   + ato_i(ji,jj) 
    1458                a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1459                v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1460                v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1461                t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 
    1462                oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1463                smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    1464                e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) ) 
    1465                ! additional condition 
    1466                IF( v_s(ji,jj,jl) <= epsi10 ) THEN 
    1467                   v_s(ji,jj,jl)   = 0._wp 
    1468                   e_s(ji,jj,1,jl) = 0._wp 
    1469                ENDIF 
    1470                ! update exchanges with ocean 
    1471                sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
    1472                wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
    1473                wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
    1474                hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    1475             END DO 
    1476          END DO 
    1477       END DO ! jl  
    1478  
    1479       ! to be sure that at_i is the sum of a_i(jl) 
    1480       at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    1481       ! 
    1482       CALL wrk_dealloc( jpi, jpj, zmask ) 
    1483       ! 
    1484    END SUBROUTINE lim_itd_me_zapsmall 
    14851244 
    14861245#else 
     
    14931252   SUBROUTINE lim_itd_me_icestrength 
    14941253   END SUBROUTINE lim_itd_me_icestrength 
    1495    SUBROUTINE lim_itd_me_sort 
    1496    END SUBROUTINE lim_itd_me_sort 
    14971254   SUBROUTINE lim_itd_me_init 
    14981255   END SUBROUTINE lim_itd_me_init 
    1499    SUBROUTINE lim_itd_me_zapsmall 
    1500    END SUBROUTINE lim_itd_me_zapsmall 
    15011256#endif 
    15021257   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.