Changeset 921


Ignore:
Timestamp:
2008-05-13T10:28:52+02:00 (13 years ago)
Author:
rblod
Message:

Correct indentation and print for debug in LIM3, see ticket #134, step I

Location:
trunk/NEMO
Files:
31 edited

Legend:

Unmodified
Added
Removed

  • trunk/NEMO/LIM_SRC_3/dom_ice.F90

    r834 r921  
    2424   INTEGER, PUBLIC ::   &  !: 
    2525      njeq , njeqm1        !: j-index of the equator if it is inside the domain 
    26       !                    !  (otherwise = jpj+10 (SH) or -10 (SH) ) 
     26   !                    !  (otherwise = jpj+10 (SH) or -10 (SH) ) 
    2727 
    2828   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 

  • trunk/NEMO/LIM_SRC_3/ice.F90

    r904 r921  
    217217      s_i_min  =  0.1 ,  &  !: minimum ice salinity (ppt) 
    218218      s_i_0    =  3.5 ,  &  !: 1st sal. value for the computation of sal .prof. 
    219                             !: (ppt) 
     219                                !: (ppt) 
    220220      s_i_1    =  4.5 ,  &  !: 2nd sal. value for the computation of sal .prof. 
    221                             !: (ppt) 
     221                                !: (ppt) 
    222222      sal_G    = 5.00 ,  &  !: restoring salinity for gravity drainage 
    223                             !: (ppt) 
     223                                !: (ppt) 
    224224      sal_F    = 2.50 ,  &  !: restoring salinity for flushing 
    225                             !: (ppt) 
     225                                !: (ppt) 
    226226      time_G   = 1.728e+06,&!: restoring time constant for gravity drainage 
    227                             !: (= 20 days, in s) 
     227                                !: (= 20 days, in s) 
    228228      time_F   = 8.640e+05,&!: restoring time constant for gravity drainage  
    229                             !: (= 10 days, in s) 
     229                                !: (= 10 days, in s) 
    230230      bulk_sal = 4.0        !: bulk salinity (ppt) in case of constant salinity 
    231231 
    232232   INTEGER , PUBLIC ::   & !!: ** ice-salinity namelist (namicesal) ** 
    233233      num_sal  = 1    ,  &  !: salinity configuration used in the model 
    234                             !: 1 - s constant in space and time 
    235                             !: 2 - prognostic salinity (s(z,t)) 
    236                             !: 3 - salinity profile, constant in time 
    237                             !: 4 - salinity variations affect only ice 
    238                             !      thermodynamics 
     234                                !: 1 - s constant in space and time 
     235                                !: 2 - prognostic salinity (s(z,t)) 
     236                                !: 3 - salinity profile, constant in time 
     237                                !: 4 - salinity variations affect only ice 
     238                                !      thermodynamics 
    239239      sal_prof = 1    ,  &  !: salinity profile or not  
    240240      thcon_i_swi = 1       !: thermal conductivity of Untersteiner (1964) (1) or 
    241                             !: Pringle et al (2007) (2) 
     241   !: Pringle et al (2007) (2) 
    242242 
    243243   REAL(wp), PUBLIC ::   & !!: ** ice-mechanical redistribution namelist (namiceitdme) 
    … …  
    249249      astar = 0.05    ,  & !!: equivalent of G* for an exponential participation function 
    250250      Hstar = 100.0   ,  & !!: thickness that determines the maximal thickness of ridged 
    251                            !!: ice 
     251                                !!: ice 
    252252      hparmeter = 0.75,  & !!: threshold thickness (m) for rafting / ridging  
    253253      Craft = 5.0     ,  & !!: coefficient for smoothness of the hyperbolic tangent in rafting 
    … …  
    256256      betas    = 1.0      , & !:: coef. for partitioning of snowfall between leads and sea ice 
    257257      kappa_i  = 1.0      , & !!: coefficient for the extinction of radiation 
    258                               !!: Grenfell et al. (2006) (m-1) 
     258                                !!: Grenfell et al. (2006) (m-1) 
    259259      nconv_i_thd = 50    , & !!: maximal number of iterations for heat diffusion 
    260260      maxer_i_thd = 1.0e-4    !!: maximal tolerated error (C) for heat diffusion 
    … …  
    264264      raftswi          = 1, & !!: rafting of ice or not                         
    265265      partfun_swi      = 1, & !!: participation function Thorndike et al. JGR75 (0)  
    266                               !!: or Lipscomb et al. JGR07 (1)  
     266                                !!: or Lipscomb et al. JGR07 (1)  
    267267      transfun_swi     = 0, & !!: transfer function of Hibler, MWR80 (0)  
    268                               !!: or Lipscomb et al., 2007 (1) 
     268                                !!: or Lipscomb et al., 2007 (1) 
    269269      brinstren_swi    = 0    !!: use brine volume to diminish ice strength 
    270270 
    … …  
    301301      t_bo   ,   &  !: Sea-Ice bottom temperature (Kelvin)       
    302302      hicifp ,   &  !: Ice production/melting 
    303       !obsolete... can be removed 
     303                                !obsolete... can be removed 
    304304      frld   ,   &  !: Leads fraction = 1-a/totalarea REFERS TO LEAD FRACTION everywhere 
    305                     !: except in the OUTPUTS!!!! 
     305                                !: except in the OUTPUTS!!!! 
    306306      pfrld  ,   &  !: Leads fraction at previous time   
    307307      phicif ,   &  !: Old ice thickness 
    … …  
    328328      fheat_res, &  !: Residual heat flux due to correction of ice thickness 
    329329      fhmec         !: Heat flux due to snow loss during compression 
    330        
     330 
    331331   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::       &  !: 
    332332      albege ,   &  !: Albedo of the snow or ice (only for outputs) 
    … …  
    334334      tauc          !: Cloud optical depth 
    335335 
    336 ! temporary arrays for dummy version of the code 
     336   ! temporary arrays for dummy version of the code 
    337337   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: 
    338338      dh_i_surf2D, dh_i_bott2D, fstbif, fsup2D, focea2D, q_s 
    … …  
    354354      sm_i   ,   &  !: Sea-Ice Bulk salinity (ppt) 
    355355      smv_i  ,   &  !: Sea-Ice Bulk salinity times volume per area (ppt.m) 
    356                     !: this is an extensive variable that has to be transported 
     356                                !: this is an extensive variable that has to be transported 
    357357      o_i    ,   &  !: Sea-Ice Age (days) 
    358358      ov_i   ,   &  !: Sea-Ice Age times volume per area (days.m) 
    … …  
    401401   !!-------------------------------------------------------------------------- 
    402402   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)     ::   &  !: 
    403          sxopw, syopw, sxxopw, syyopw, sxyopw          !: open water in sea ice 
     403      sxopw, syopw, sxxopw, syyopw, sxyopw          !: open water in sea ice 
    404404 
    405405   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   &  !: 
    406          sxice, syice, sxxice, syyice, sxyice,      &  !: ice thickness moments for advection 
    407          sxsn,  sysn,  sxxsn,  syysn,  sxysn,       &  !: snow thickness 
    408          sxa,   sya,   sxxa,   syya,   sxya,        &  !: lead fraction 
    409          sxc0,  syc0,  sxxc0,  syyc0,  sxyc0,       &  !: snow thermal content 
    410          sxsal, sysal, sxxsal, syysal, sxysal,      &  !: ice salinity 
    411          sxage, syage, sxxage, syyage, sxyage          !: ice age 
     406      sxice, syice, sxxice, syyice, sxyice,      &  !: ice thickness moments for advection 
     407      sxsn,  sysn,  sxxsn,  syysn,  sxysn,       &  !: snow thickness 
     408      sxa,   sya,   sxxa,   syya,   sxya,        &  !: lead fraction 
     409      sxc0,  syc0,  sxxc0,  syyc0,  sxyc0,       &  !: snow thermal content 
     410      sxsal, sysal, sxxsal, syysal, sxysal,      &  !: ice salinity 
     411      sxage, syage, sxxage, syyage, sxyage          !: ice age 
    412412 
    413413   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jkmax,jpl) ::   &  !: 
    414          sxe ,  sye ,  sxxe ,  syye ,  sxye            !: ice layers heat content 
     414      sxe ,  sye ,  sxxe ,  syye ,  sxye            !: ice layers heat content 
    415415 
    416416   !!-------------------------------------------------------------------------- 
    … …  
    446446   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jkmax,jpl) ::   &  !: 
    447447      d_e_i_thd, d_e_i_trp 
    448     
     448 
    449449   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::       &  !: ice velocity  
    450450      d_u_ice_dyn, d_v_ice_dyn 
    … …  
    459459   INTEGER, PUBLIC, DIMENSION(jpm,2)              ::   &  !: 
    460460      ice_cat_bounds !: Matrix containing the integer upper and  
    461                      !: lower boundaries of ice thickness categories 
     461   !: lower boundaries of ice thickness categories 
    462462 
    463463   ! REMOVE 
    … …  
    474474   REAL(wp), PUBLIC, DIMENSION(0:jpl,jpm)         ::   &  !: 
    475475      hi_max_typ     !: Boundary of ice thickness categories  
    476                      !:in thickness space (same but specific for each ice type) 
     476   !:in thickness space (same but specific for each ice type) 
     477 
     478   !!-------------------------------------------------------------------------- 
     479   !! * Ice Run 
     480   !!-------------------------------------------------------------------------- 
     481   !! Namelist namicerun read in iceini 
     482   LOGICAL , PUBLIC  ::     & !!! ** init namelist (namicerun) ** 
     483      ln_limdyn   = .TRUE., & !: flag for ice dynamics (T) or not (F) 
     484      ln_nicep    = .TRUE.    !: flag for sea-ice points output (T) or not (F) 
     485   REAL(wp), PUBLIC  ::   &  !: 
     486      hsndif = 0.e0    ,  &  !: computation of temp. in snow (0) or not (9999) 
     487      hicdif = 0.e0    ,  &  !: computation of temp. in ice (0) or not (9999) 
     488      cai    = 1.40e-3 ,  &  !: atmospheric drag over sea ice 
     489      cao    = 1.00e-3       !: atmospheric drag over ocean 
     490   REAL(wp), PUBLIC, DIMENSION(2)  ::  &  !: 
     491      acrit  = (/ 1.e-06 , 1.e-06 /)    !: minimum fraction for leads in 
     492   !                                   !  north and south hemisphere 
    477493 
    478494   !!-------------------------------------------------------------------------- 

  • trunk/NEMO/LIM_SRC_3/iceini.F90

    r888 r921  
    3232 
    3333   !! * Share Module variables 
    34    LOGICAL , PUBLIC  ::     & !!! ** init namelist (namicerun) ** 
    35       ln_limdyn   = .TRUE., & !: flag for ice dynamics (T) or not (F) 
    36       ln_nicep    = .TRUE.    !: flag for sea-ice points output (T) or not (F) 
    3734   INTEGER , PUBLIC  ::   &  !: 
    3835      nstart ,            &  !: iteration number of the begining of the run  
    … …  
    4138      numit                  !: iteration number 
    4239   REAL(wp), PUBLIC  ::   &  !: 
    43       hsndif = 0.e0    ,  &  !: computation of temp. in snow (0) or not (9999) 
    44       hicdif = 0.e0    ,  &  !: computation of temp. in ice (0) or not (9999) 
    45       tpstot           ,  &  !: time of the run in seconds 
    46       cai    = 1.40e-3 ,  &  !: atmospheric drag over sea ice 
    47       cao    = 1.00e-3       !: atmospheric drag over ocean 
    48    REAL(wp), PUBLIC, DIMENSION(2)  ::  &  !: 
    49       acrit  = (/ 1.e-06 , 1.e-06 /)    !: minimum fraction for leads in  
    50       !                                   !  north and south hemisphere 
     40      tpstot                 !: time of the run in seconds 
    5141   !!---------------------------------------------------------------------- 
    5242   !!   LIM 3.0,  UCL-ASTR-LOCEAN-IPSL (2008)  
    … …  
    7262 
    7363      CALL ice_run                    !  read in namelist some run parameters 
    74                   
     64 
    7565      ! Louvain la Neuve Ice model 
    7666      IF( nacc == 1 ) THEN 
    77           dtsd2   = nn_fsbc * rdtmin * 0.5 
    78           rdt_ice = nn_fsbc * rdtmin 
     67         dtsd2   = nn_fsbc * rdtmin * 0.5 
     68         rdt_ice = nn_fsbc * rdtmin 
    7969      ELSE 
    80           dtsd2   = nn_fsbc * rdt * 0.5 
    81           rdt_ice = nn_fsbc * rdt 
     70         dtsd2   = nn_fsbc * rdt * 0.5 
     71         rdt_ice = nn_fsbc * rdt 
    8272      ENDIF 
    8373 
    8474      CALL lim_msh                    ! ice mesh initialization 
    85       
     75 
    8676      CALL lim_itd_ini                ! initialize the ice thickness 
    87                                       ! distribution 
     77      ! distribution 
    8878      ! Initial sea-ice state 
    8979      IF( .NOT.ln_rstart ) THEN 
    … …  
    9282         CALL lim_istate              ! start from rest: sea-ice deduced from sst 
    9383         CALL lim_var_agg(1)          ! aggregate category variables in 
    94                                       ! bulk variables 
     84         ! bulk variables 
    9585         CALL lim_var_glo2eqv         ! convert global variables in equivalent 
    96                                       ! variables 
     86         ! variables 
    9787      ELSE 
    9888         CALL lim_rst_read            ! start from a restart file 
    … …  
    10898      alb_ice(:,:,:) = albege(:,:)      ! sea-ice albedo 
    10999# endif 
    110        
     100 
    111101      nstart = numit  + nn_fsbc       
    112102      nitrun = nitend - nit000 + 1  
    … …  
    138128      REWIND ( numnam_ice ) 
    139129      READ   ( numnam_ice , namicerun ) 
     130      ln_nicep = ln_nicep .AND. lwp 
    140131      IF(lwp) THEN 
    141132         WRITE(numout,*) 
    … …  
    150141         WRITE(numout,*) '   Several ice points in the ice or not in ocean.output = ', ln_nicep 
    151142      ENDIF 
    152        
     143 
    153144   END SUBROUTINE ice_run 
    154145 
    155146   SUBROUTINE lim_itd_ini 
    156         !!------------------------------------------------------------------ 
    157         !!                ***  ROUTINE lim_itd_ini *** 
    158         !! ** Purpose : 
    159         !!            Initializes the ice thickness distribution 
    160         !! ** Method  : 
    161         !!            Very simple. Currently there are no ice types in the 
    162         !!            model... 
    163         !! 
    164         !! ** Arguments : 
    165         !!           kideb , kiut : Starting and ending points on which the 
    166         !!                         the computation is applied 
    167         !! 
    168         !! ** Inputs / Ouputs : (global commons) 
    169         !! 
    170         !! ** External : 
    171         !! 
    172         !! ** References : 
    173         !! 
    174         !! ** History : 
    175         !!           (12-2005) Martin Vancoppenolle 
    176         !! 
    177         !!------------------------------------------------------------------ 
    178         !! * Arguments 
    179  
    180        !! * Local variables 
    181        INTEGER ::   jl,       &   ! ice category dummy loop index 
    182                     jm            ! ice types    dummy loop index 
    183  
    184        REAL(wp)  ::           &  ! constant values 
    185           zeps      =  1.0e-10,   & ! 
    186           zc1                 ,   & ! 
    187           zc2                 ,   & ! 
    188           zc3                 ,   & ! 
    189           zx1 
    190  
    191        WRITE(numout,*) 'lim_itd_ini : Initialization of ice thickness distribution ' 
    192        WRITE(numout,*) '~~~~~~~~~~~~' 
    193  
    194 !!-- End of declarations 
    195 !!------------------------------------------------------------------------------ 
    196  
    197 !------------------------------------------------------------------------------! 
    198 ! 1) Ice thickness distribution parameters initialization     
    199 !------------------------------------------------------------------------------! 
     147      !!------------------------------------------------------------------ 
     148      !!                ***  ROUTINE lim_itd_ini *** 
     149      !! ** Purpose : 
     150      !!            Initializes the ice thickness distribution 
     151      !! ** Method  : 
     152      !!            Very simple. Currently there are no ice types in the 
     153      !!            model... 
     154      !! 
     155      !! ** Arguments : 
     156      !!           kideb , kiut : Starting and ending points on which the 
     157      !!                         the computation is applied 
     158      !! 
     159      !! ** Inputs / Ouputs : (global commons) 
     160      !! 
     161      !! ** External : 
     162      !! 
     163      !! ** References : 
     164      !! 
     165      !! ** History : 
     166      !!           (12-2005) Martin Vancoppenolle 
     167      !! 
     168      !!------------------------------------------------------------------ 
     169      !! * Arguments 
     170 
     171      !! * Local variables 
     172      INTEGER ::   jl,       &   ! ice category dummy loop index 
     173         jm            ! ice types    dummy loop index 
     174 
     175      REAL(wp)  ::           &  ! constant values 
     176         zeps      =  1.0e-10,   & ! 
     177         zc1                 ,   & ! 
     178         zc2                 ,   & ! 
     179         zc3                 ,   & ! 
     180         zx1 
     181 
     182      WRITE(numout,*) 'lim_itd_ini : Initialization of ice thickness distribution ' 
     183      WRITE(numout,*) '~~~~~~~~~~~~' 
     184 
     185      !!-- End of declarations 
     186      !!------------------------------------------------------------------------------ 
     187 
     188      !------------------------------------------------------------------------------! 
     189      ! 1) Ice thickness distribution parameters initialization     
     190      !------------------------------------------------------------------------------! 
    200191 
    201192      !- Types boundaries (integer) 
    … …  
    266257      tn_ice(:,:,:) = t_su(:,:,:) 
    267258 
    268     END SUBROUTINE lim_itd_ini 
     259   END SUBROUTINE lim_itd_ini 
    269260 
    270261#else 

  • trunk/NEMO/LIM_SRC_3/limadv.F90

    r888 r921  
    6666         pdf ,       &  ! ??? 
    6767         pcrh           ! = 1. : lim_adv_x is called before lim_adv_y 
    68          !              ! = 0. : lim_adv_x is called after  lim_adv_y 
     68      !              ! = 0. : lim_adv_x is called after  lim_adv_y 
    6969      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  ::  & 
    7070         put            ! i-direction ice velocity at ocean U-point (m/s) 
    … …  
    114114      !  Calculate fluxes and moments between boxes i<-->i+1               
    115115      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0  
    116 !i bug   DO ji = 1, jpim1  
    117 !i    DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0  
     116         !i bug   DO ji = 1, jpim1  
     117         !i    DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0  
    118118         DO ji = 1, jpi 
    119119            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, put(ji,jj) ) ) 
    … …  
    142142 
    143143      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0. 
    144 !i    DO jj = 1, fs_jpjm1                   !  Flux from i+1 to i when u LT 0. 
     144         !i    DO jj = 1, fs_jpjm1                   !  Flux from i+1 to i when u LT 0. 
    145145         DO ji = 1, fs_jpim1 
    146146            zalf          = MAX( rzero, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj)  
    … …  
    228228      CALL lbc_lnk( psxy, 'T', 1. ) 
    229229 
    230      IF(ln_ctl) THEN 
     230      IF(ln_ctl) THEN 
    231231         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ') 
    232232         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ') 
    233233         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ') 
    234234         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_x: psxy :') 
    235      ENDIF 
     235      ENDIF 
    236236 
    237237   END SUBROUTINE lim_adv_x 
    … …  
    260260         pdf,        &  ! ??? 
    261261         pcrh           ! = 1. : lim_adv_x is called before lim_adv_y 
    262          !              ! = 0. : lim_adv_x is called after  lim_adv_y 
     262      !              ! = 0. : lim_adv_x is called after  lim_adv_y 
    263263      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: & 
    264264         pvt            ! j-direction ice velocity at ocean V-point (m/s) 
    … …  
    285285      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection. 
    286286 
    287        DO jj = 1, jpj 
    288           DO ji = 1, jpi 
    289              zslpmax = MAX( rzero, ps0(ji,jj) ) 
    290              zs1max  = 1.5 * zslpmax 
    291              zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) ) 
    292              zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
    293                 &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  ) 
    294              zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
    295              ps0 (ji,jj) = zslpmax   
    296              psx (ji,jj)  = psx (ji,jj) * zin0 
    297              psxx(ji,jj)  = psxx(ji,jj) * zin0 
    298              psy (ji,jj) = zs1new * zin0 
    299              psyy(ji,jj) = zs2new * zin0 
    300              psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0 
    301           END DO 
    302        END DO 
    303  
    304        !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
    305        psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 
    306  
    307        !  Calculate fluxes and moments between boxes j<-->j+1               
    308 !!bug  DO jj = 2, jpjm1 
    309        DO jj = 1, jpj 
    310           DO ji = 1, jpi 
    311 !!bug     DO ji = 1, jpim1 
    312              !  Flux from j to j+1 WHEN v GT 0    
    313              zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) ) 
    314              zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj) 
    315              zalfq        =  zalf * zalf 
    316              zalf1        =  1.0 - zalf 
    317              zalf1q       =  zalf1 * zalf1 
    318              zfm (ji,jj)  =  zalf  * psm(ji,jj) 
    319              zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) )  
    320              zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) ) 
    321              zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj) 
    322              zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) ) 
    323              zfxy(ji,jj)  =  zalfq * psxy(ji,jj) 
    324              zfxx(ji,jj)  =  zalf  * psxx(ji,jj) 
    325  
    326              !  Readjust moments remaining in the box. 
    327              psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj) 
    328              ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj) 
    329              psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) ) 
    330              psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj) 
    331              psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj) 
    332              psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj) 
    333              psxy(ji,jj)  =  zalf1q * psxy(ji,jj) 
    334           END DO 
    335        END DO 
    336  
    337        DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
    338           DO ji = 1, jpi 
    339 !i     DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
    340 !i        DO ji = 2, jpim1 
    341              zalf          = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1)  
    342              zalg  (ji,jj) = zalf 
    343              zalfq         = zalf * zalf 
    344              zalf1         = 1.0 - zalf 
    345              zalg1 (ji,jj) = zalf1 
    346              zalf1q        = zalf1 * zalf1 
    347              zalg1q(ji,jj) = zalf1q 
    348              zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji,jj+1) 
    349              zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 
    350              zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 
    351              zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * zalfq * psyy(ji,jj+1) 
    352              zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 
    353              zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1) 
    354              zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * psxx(ji,jj+1) 
    355           END DO 
    356        END DO 
    357   
    358        !  Readjust moments remaining in the box.  
    359        DO jj = 2, jpj 
    360           DO ji = 1, jpi 
    361              zbt  =         zbet(ji,jj-1) 
    362              zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    363              psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 
    364              ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 
    365              psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) ) 
    366              psyy(ji,jj) = zalg1 (ji,jj-1)  * zalg1q(ji,jj-1) * psyy(ji,jj) 
    367              psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 
    368              psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) 
    369              psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj) 
    370           END DO 
    371        END DO 
    372  
    373        !   Put the temporary moments into appropriate neighboring boxes.     
    374        DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0. 
    375           DO ji = 1, jpi 
    376              zbt  =         zbet(ji,jj-1) 
    377              zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
    378              psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj)  
    379              zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj)  
    380              zalf1       = 1.0 - zalf 
    381              ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1) 
    382              ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj) 
    383  
    384              psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   & 
    385                 &        + zbt1 * psy(ji,jj)   
    386  
    387              psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             & 
    388                 &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   &  
    389                 &        + zbt1 * psyy(ji,jj) 
    390  
    391              psxy(ji,jj) = zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)              & 
    392                                   + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) )   & 
    393                          + zbt1 * psxy(ji,jj) 
    394              psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 
    395              psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) 
    396           END DO 
    397        END DO 
    398  
    399        DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0. 
    400           DO ji = 1, jpi 
    401              zbt  =         zbet(ji,jj) 
    402              zbt1 = ( 1.0 - zbet(ji,jj) ) 
    403              psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 
    404              zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj) 
    405              zalf1       = 1.0 - zalf 
    406              ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 
    407              ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
    408              psy(ji,jj)  = zbt  * psy(ji,jj)  & 
    409                 &        + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 
    410              psyy(ji,jj) = zbt  * psyy(ji,jj)  & 
    411                 &        + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   & 
    412                 &                 + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) 
    413              psxy(ji,jj) = zbt  * psxy(ji,jj)   & 
    414                 &        + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)   & 
    415                 &                 + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) ) 
    416              psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 
    417              psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 
    418           END DO 
    419        END DO 
     287      DO jj = 1, jpj 
     288         DO ji = 1, jpi 
     289            zslpmax = MAX( rzero, ps0(ji,jj) ) 
     290            zs1max  = 1.5 * zslpmax 
     291            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) ) 
     292            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   & 
     293               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  ) 
     294            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask 
     295            ps0 (ji,jj) = zslpmax   
     296            psx (ji,jj)  = psx (ji,jj) * zin0 
     297            psxx(ji,jj)  = psxx(ji,jj) * zin0 
     298            psy (ji,jj) = zs1new * zin0 
     299            psyy(ji,jj) = zs2new * zin0 
     300            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * zin0 
     301         END DO 
     302      END DO 
     303 
     304      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 
     305      psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 ) 
     306 
     307      !  Calculate fluxes and moments between boxes j<-->j+1               
     308      !!bug  DO jj = 2, jpjm1 
     309      DO jj = 1, jpj 
     310         DO ji = 1, jpi 
     311            !!bug     DO ji = 1, jpim1 
     312            !  Flux from j to j+1 WHEN v GT 0    
     313            zbet(ji,jj)  =  MAX( rzero, SIGN( rone, pvt(ji,jj) ) ) 
     314            zalf         =  MAX( rzero, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj) 
     315            zalfq        =  zalf * zalf 
     316            zalf1        =  1.0 - zalf 
     317            zalf1q       =  zalf1 * zalf1 
     318            zfm (ji,jj)  =  zalf  * psm(ji,jj) 
     319            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) )  
     320            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) ) 
     321            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj) 
     322            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) ) 
     323            zfxy(ji,jj)  =  zalfq * psxy(ji,jj) 
     324            zfxx(ji,jj)  =  zalf  * psxx(ji,jj) 
     325 
     326            !  Readjust moments remaining in the box. 
     327            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj) 
     328            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj) 
     329            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) ) 
     330            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj) 
     331            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj) 
     332            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj) 
     333            psxy(ji,jj)  =  zalf1q * psxy(ji,jj) 
     334         END DO 
     335      END DO 
     336 
     337      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
     338         DO ji = 1, jpi 
     339            !i     DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0. 
     340            !i        DO ji = 2, jpim1 
     341            zalf          = ( MAX(rzero, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1)  
     342            zalg  (ji,jj) = zalf 
     343            zalfq         = zalf * zalf 
     344            zalf1         = 1.0 - zalf 
     345            zalg1 (ji,jj) = zalf1 
     346            zalf1q        = zalf1 * zalf1 
     347            zalg1q(ji,jj) = zalf1q 
     348            zfm   (ji,jj) = zfm (ji,jj) + zalf  * psm(ji,jj+1) 
     349            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0(ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) ) 
     350            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy(ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) ) 
     351            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  * zalfq * psyy(ji,jj+1) 
     352            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx(ji,jj+1) - zalf1 * psxy(ji,jj+1) ) 
     353            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1) 
     354            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  * psxx(ji,jj+1) 
     355         END DO 
     356      END DO 
     357 
     358      !  Readjust moments remaining in the box.  
     359      DO jj = 2, jpj 
     360         DO ji = 1, jpi 
     361            zbt  =         zbet(ji,jj-1) 
     362            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
     363            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) ) 
     364            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) ) 
     365            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) ) 
     366            psyy(ji,jj) = zalg1 (ji,jj-1)  * zalg1q(ji,jj-1) * psyy(ji,jj) 
     367            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) ) 
     368            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) ) 
     369            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj) 
     370         END DO 
     371      END DO 
     372 
     373      !   Put the temporary moments into appropriate neighboring boxes.     
     374      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0. 
     375         DO ji = 1, jpi 
     376            zbt  =         zbet(ji,jj-1) 
     377            zbt1 = ( 1.0 - zbet(ji,jj-1) ) 
     378            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj)  
     379            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj)  
     380            zalf1       = 1.0 - zalf 
     381            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1) 
     382            ps0(ji,jj)  = zbt * (ps0(ji,jj) + zf0(ji,jj-1)) + zbt1 * ps0(ji,jj) 
     383 
     384            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   & 
     385               &        + zbt1 * psy(ji,jj)   
     386 
     387            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             & 
     388               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   &  
     389               &        + zbt1 * psyy(ji,jj) 
     390 
     391            psxy(ji,jj) = zbt  * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)              & 
     392               + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) ) )   & 
     393               + zbt1 * psxy(ji,jj) 
     394            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj) 
     395            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj) 
     396         END DO 
     397      END DO 
     398 
     399      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0. 
     400         DO ji = 1, jpi 
     401            zbt  =         zbet(ji,jj) 
     402            zbt1 = ( 1.0 - zbet(ji,jj) ) 
     403            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) ) 
     404            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj) 
     405            zalf1       = 1.0 - zalf 
     406            ztemp       = -zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj) 
     407            ps0(ji,jj)  = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) ) 
     408            psy(ji,jj)  = zbt  * psy(ji,jj)  & 
     409               &        + zbt1 * ( zalf*zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp ) 
     410            psyy(ji,jj) = zbt  * psyy(ji,jj)  & 
     411               &        + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   & 
     412               &                 + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) ) + ( zalf1 - zalf ) * ztemp ) ) 
     413            psxy(ji,jj) = zbt  * psxy(ji,jj)   & 
     414               &        + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)   & 
     415               &                 + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) ) ) 
     416            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) ) 
     417            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) ) 
     418         END DO 
     419      END DO 
    420420 
    421421      !-- Lateral boundary conditions 
    … …  
    428428      CALL lbc_lnk( psxy, 'T', 1. ) 
    429429 
    430      IF(ln_ctl) THEN 
     430      IF(ln_ctl) THEN 
    431431         CALL prt_ctl(tab2d_1=psm  , clinfo1=' lim_adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ') 
    432432         CALL prt_ctl(tab2d_1=psx  , clinfo1=' lim_adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ') 
    433433         CALL prt_ctl(tab2d_1=psy  , clinfo1=' lim_adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ') 
    434434         CALL prt_ctl(tab2d_1=psxy , clinfo1=' lim_adv_y: psxy :') 
    435      ENDIF 
     435      ENDIF 
    436436 
    437437   END SUBROUTINE lim_adv_y 

  • trunk/NEMO/LIM_SRC_3/limcons.F90

    r834 r921  
    4242CONTAINS 
    4343 
    44 !=============================================================================== 
     44   !=============================================================================== 
    4545 
    4646   SUBROUTINE lim_column_sum(nsum,xin,xout) 
    47 !     !!------------------------------------------------------------------- 
    48 !     !!               ***  ROUTINE lim_column_sum *** 
    49 !     !! 
    50 !     !! ** Purpose : Compute the sum of xin over nsum categories 
    51 !     !! 
    52 !     !! ** Method  : Arithmetics 
    53 !     !! 
    54 !     !! ** Action  : Gets xin(ji,jj,jl) and computes xout(ji,jj) 
    55 !     !! 
    56 !     !! History : 
    57 !     !!   author: William H. Lipscomb, LANL 
    58 !     !!   2.1  !  04-06  (M. Vancoppenolle)   Energy Conservation  
    59 !     !!--------------------------------------------------------------------- 
    60 !     !! * Local variables 
     47      !     !!------------------------------------------------------------------- 
     48      !     !!               ***  ROUTINE lim_column_sum *** 
     49      !     !! 
     50      !     !! ** Purpose : Compute the sum of xin over nsum categories 
     51      !     !! 
     52      !     !! ** Method  : Arithmetics 
     53      !     !! 
     54      !     !! ** Action  : Gets xin(ji,jj,jl) and computes xout(ji,jj) 
     55      !     !! 
     56      !     !! History : 
     57      !     !!   author: William H. Lipscomb, LANL 
     58      !     !!   2.1  !  04-06  (M. Vancoppenolle)   Energy Conservation  
     59      !     !!--------------------------------------------------------------------- 
     60      !     !! * Local variables 
    6161      INTEGER, INTENT(in) ::     & 
    62            nsum                  ! number of categories/layers 
     62         nsum                  ! number of categories/layers 
    6363 
    6464      REAL (wp), DIMENSION(jpi, jpj, jpl), INTENT(IN) ::   & 
    65            xin                   ! input field 
     65         xin                   ! input field 
    6666 
    6767      REAL (wp), DIMENSION(jpi, jpj), INTENT(OUT) ::  & 
    68            xout                  ! output field 
     68         xout                  ! output field 
    6969      INTEGER ::                 & 
    70            ji, jj, jl         ! horizontal indices 
    71  
    72 !     !!--------------------------------------------------------------------- 
    73 !     WRITE(numout,*) ' lim_column_sum ' 
    74 !     WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
     70         ji, jj, jl         ! horizontal indices 
     71 
     72      !     !!--------------------------------------------------------------------- 
     73      !     WRITE(numout,*) ' lim_column_sum ' 
     74      !     WRITE(numout,*) ' ~~~~~~~~~~~~~~ ' 
    7575 
    7676      xout(:,:) = 0.00 
    … …  
    8686   END SUBROUTINE lim_column_sum 
    8787 
    88 !=============================================================================== 
     88   !=============================================================================== 
    8989 
    9090   SUBROUTINE lim_column_sum_energy(nsum,nlay,xin,xout) 
    … …  
    106106      !! * Local variables 
    107107      INTEGER, INTENT(in) ::  & 
    108            nsum,              &  !: number of categories 
    109            nlay                  !: number of vertical layers 
     108         nsum,              &  !: number of categories 
     109         nlay                  !: number of vertical layers 
    110110 
    111111      REAL (wp), DIMENSION(jpi, jpj, jkmax, jpl), INTENT(IN) :: & 
    112            xin                   !: input field 
     112         xin                   !: input field 
    113113 
    114114      REAL (wp), DIMENSION(jpi, jpj), INTENT(OUT) ::  & 
    115            xout                  !: output field 
     115         xout                  !: output field 
    116116 
    117117      INTEGER ::              & 
    118            ji, jj,            &  !: horizontal indices 
    119            jk, jl                !: layer and category  indices 
    120       !!--------------------------------------------------------------------- 
    121  
    122 !     WRITE(numout,*) ' lim_column_sum_energy ' 
    123 !     WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~ ' 
     118         ji, jj,            &  !: horizontal indices 
     119         jk, jl                !: layer and category  indices 
     120      !!--------------------------------------------------------------------- 
     121 
     122      !     WRITE(numout,*) ' lim_column_sum_energy ' 
     123      !     WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~ ' 
    124124 
    125125      xout(:,:) = 0.00 
    … …  
    137137   END SUBROUTINE lim_column_sum_energy 
    138138 
    139 !=============================================================================== 
    140     
     139   !=============================================================================== 
     140 
    141141   SUBROUTINE lim_cons_check(x1, x2, max_err, fieldid) 
    142142      !!------------------------------------------------------------------- 
    … …  
    206206                  WRITE (numout,*) ' Point         : ', ji, jj  
    207207                  WRITE (numout,*) ' lat, lon      : ', gphit(ji,jj), &  
    208                                                         glamt(ji,jj) 
     208                     glamt(ji,jj) 
    209209                  WRITE (numout,*) ' Initial value : ', x1(ji,jj) 
    210210                  WRITE (numout,*) ' Final value   : ', x2(ji,jj) 

  • trunk/NEMO/LIM_SRC_3/limdia.F90

    r895 r921  
    9595      !!------------------------------------------------------------------- 
    9696      !! * Local variables 
    97        INTEGER  ::   jv,ji,jj,jl ! dummy loop indices 
    98        REAL(wp), DIMENSION(jpinfmx) ::  &  
    99           vinfor           ! temporary working space  
    100        REAL(wp) ::    & 
    101           zshift_date   , & ! date from the minimum ice extent 
    102           zday, zday_min, & ! current day, day of minimum extent 
    103           zafy, zamy,     & ! temporary area of fy and my ice 
    104           zindb 
    105        !!------------------------------------------------------------------- 
    106  
    107        ! 0) date from the minimum of ice extent 
    108        !--------------------------------------- 
    109        zday_min = 273.0        ! zday_min = date of minimum extent, here September 30th 
    110        zday = FLOAT(numit-nit000) * rdt_ice / ( 86400.0 * FLOAT(nn_fsbc) ) 
    111        IF (zday.GT.zday_min) THEN  
    112           zshift_date  =  zday - zday_min 
    113        ELSE 
    114           zshift_date  =  zday - (365.0 - zday_min) 
    115        ENDIF 
    116  
    117        IF( numit == nstart )   CALL lim_dia_init   ! initialisation of ice_evolu file       
    118  
    119        ! temporal diagnostics  
    120        vinfor(1) = REAL(numit) 
    121        vinfor(2) = nyear 
    122   
    123        ! put everything to zero 
    124        DO jv = nbvt + 1, nvinfo 
    125           vinfor(jv) = 0.0 
    126        END DO 
    127  
    128        !!------------------------------------------------------------------- 
    129        !! 1) Northern hemisphere 
    130        !!------------------------------------------------------------------- 
    131        !! 1.1) Diagnostics independent on age 
    132        !!------------------------------------ 
    133        DO jj = njeq, jpjm1 
    134           DO ji = fs_2, fs_jpim1   ! vector opt. 
    135              IF( tms(ji,jj) == 1 ) THEN 
    136                 vinfor(3)  = vinfor(3)  + at_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice area 
    137                 IF (at_i(ji,jj).GT.0.15) vinfor(5) = vinfor(5) + aire(ji,jj) / 1.0e12 !ice extent 
    138                 vinfor(7)  = vinfor(7)  + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume 
    139                 vinfor(9)  = vinfor(9)  + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume 
    140                 vinfor(15) = vinfor(15) + ot_i(ji,jj) *vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age 
    141                 vinfor(29) = vinfor(29) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity 
    142                 ! the computation of this diagnostic is not reliable 
    143                 vinfor(31) = vinfor(31) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + &  
    144                                                         v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12  
    145                 vinfor(53) = vinfor(53) + emps(ji,jj)*aire(ji,jj) / 1.0e12 !salt flux 
    146                 vinfor(55) = vinfor(55) + fsbri(ji,jj)*aire(ji,jj) / 1.0e12 !brine drainage flux 
    147                 vinfor(57) = vinfor(57) + fseqv(ji,jj)*aire(ji,jj) / 1.0e12 !equivalent salt flux 
    148                 vinfor(59) = vinfor(59) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SST 
    149                 vinfor(61) = vinfor(61) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SSS 
    150                 vinfor(65) = vinfor(65) + et_s(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12  ! snow temperature 
    151                 vinfor(67) = vinfor(67) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12       ! ice heat content 
    152                 vinfor(69) = vinfor(69) + v_i(ji,jj,1)*aire(ji,jj) / 1.0e12 !ice volume 
    153                 vinfor(71) = vinfor(71) + v_i(ji,jj,2)*aire(ji,jj) / 1.0e12 !ice volume 
    154                 vinfor(73) = vinfor(73) + v_i(ji,jj,3)*aire(ji,jj) / 1.0e12 !ice volume 
    155                 vinfor(75) = vinfor(75) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume 
    156                 vinfor(77) = vinfor(77) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume 
    157                 vinfor(79) = 0.0 
    158                 vinfor(81) = vinfor(81) + emp(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux 
    159              ENDIF 
    160           END DO 
    161        END DO 
    162  
    163        DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
    164           DO jj = njeq, jpjm1 
    165              DO ji = fs_2, fs_jpim1   ! vector opt. 
    166                 IF( tms(ji,jj) == 1 ) THEN 
    167                    vinfor(11) = vinfor(11) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !undef def ice volume 
    168                 ENDIF 
    169              END DO 
    170           END DO 
    171        END DO 
    172  
    173        vinfor(13) = 0.0 
    174  
    175        vinfor(15) = vinfor(15) / MAX(vinfor(7),epsi06) ! these have to be divided by total ice volume to have the 
    176        vinfor(29) = vinfor(29) / MAX(vinfor(7),epsi06) ! right value 
    177        vinfor(31) = SQRT( vinfor(31) / MAX( vinfor(7) , epsi06 ) ) 
    178        vinfor(67) = vinfor(67) / MAX(vinfor(7),epsi06) 
    179  
    180        vinfor(53) = vinfor(53) / MAX(vinfor(5),epsi06) ! these have to be divided by total ice extent to have the 
    181        vinfor(55) = vinfor(55) / MAX(vinfor(5),epsi06) ! right value  
    182        vinfor(57) = vinfor(57) / MAX(vinfor(5),epsi06) !  
    183        vinfor(79) = vinfor(79) / MAX(vinfor(5),epsi06) ! 
    184  
    185        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(3))) ! 
    186        vinfor(59) = zindb*vinfor(59) / MAX(vinfor(3),epsi06) ! divide by ice area 
    187        vinfor(61) = zindb*vinfor(61) / MAX(vinfor(3),epsi06) ! 
    188  
    189        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(9))) ! 
    190        vinfor(65) = zindb*vinfor(65) / MAX(vinfor(9),epsi06) ! divide it by snow volume 
    191  
    192  
    193        DO jl = 1, jpl 
    194           DO jj = njeq, jpjm1 
    195              DO ji = fs_2, fs_jpim1   ! vector opt. 
    196                 IF( tms(ji,jj) == 1 ) THEN 
    197                    vinfor(33) = vinfor(33) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume 
    198                    vinfor(35) = vinfor(35) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume 
    199                 ENDIF 
    200              END DO 
    201           END DO 
    202        END DO 
    203  
    204        DO jj = njeq, jpjm1 
    205           DO ji = fs_2, fs_jpim1   ! vector opt. 
    206                 IF( tms(ji,jj) == 1 ) THEN 
    207                    vinfor(37) = vinfor(37) + diag_sni_gr(ji,jj)*aire(ji,jj) / 1.0e12 !th growth rates 
    208                    vinfor(39) = vinfor(39) + diag_lat_gr(ji,jj)*aire(ji,jj) / 1.0e12  
    209                    vinfor(41) = vinfor(41) + diag_bot_gr(ji,jj)*aire(ji,jj) / 1.0e12 
    210                    vinfor(43) = vinfor(43) + diag_dyn_gr(ji,jj)*aire(ji,jj) / 1.0e12  
    211                    vinfor(45) = vinfor(45) + dv_dt_thd(ji,jj,5)*aire(ji,jj) / 1.0e12 
    212                    vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) / 1.0e12 / rdt_ice ! volume acc in OW 
    213                 ENDIF 
    214           END DO 
    215        END DO 
    216  
    217        DO jl = 1, jpl 
    218           DO jj = njeq, jpjm1 
    219              DO ji = fs_2, fs_jpim1   ! vector opt. 
    220                 IF( tms(ji,jj) == 1 ) THEN 
    221                    vinfor(63) = vinfor(63) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 
    222                 ENDIF 
    223              END DO 
    224           END DO 
    225        END DO 
    226        vinfor(63) = vinfor(63) / MAX(vinfor(3),epsi06) ! these have to be divided by total ice area 
    227  
    228        !! 1.2) Diagnostics dependent on age 
    229        !!------------------------------------ 
    230        DO jj = njeq, jpjm1 
    231           DO ji = fs_2, fs_jpim1   ! vector opt. 
    232              IF( tms(ji,jj) == 1 ) THEN 
    233                 zafy = 0.0 
    234                 zamy = 0.0 
    235                 DO jl = 1, jpl 
    236                    IF ((o_i(ji,jj,jl) - zshift_date).LT.0.0) THEN 
    237                       vinfor(17) = vinfor(17) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice area 
    238                       vinfor(25) = vinfor(25) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice volume 
    239                       vinfor(49) = vinfor(49) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity 
    240                       zafy = zafy + a_i(ji,jj,jl) 
    241                    ENDIF 
    242                    IF ((o_i(ji,jj,jl) - zshift_date).GT.0.0) THEN 
    243                       vinfor(19) = vinfor(19) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12    ! MY ice area 
    244                       vinfor(27) = vinfor(27) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! MY ice volume 
    245                       vinfor(51) = vinfor(51) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !MY ice salinity 
    246                       zamy = zamy + a_i(ji,jj,jl) 
    247                    ENDIF 
    248                 END DO 
    249                 IF ((at_i(ji,jj).GT.0.15).AND.(zafy.GT.zamy)) THEN 
    250                    vinfor(21) = vinfor(21) + aire(ji,jj) / 1.0e12 ! Seasonal ice extent 
    251                 ENDIF 
    252                 IF ((at_i(ji,jj).GT.0.15).AND.(zafy.LE.zamy)) THEN 
    253                    vinfor(23) = vinfor(23) + aire(ji,jj) / 1.0e12 ! Perennial ice extent 
    254                 ENDIF 
    255              ENDIF 
    256           END DO 
    257        END DO 
    258        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(25))) !=0 if no multiyear ice 1 if yes 
    259        vinfor(49) = zindb*vinfor(49) / MAX(vinfor(25),epsi06) 
    260        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(27))) !=0 if no multiyear ice 1 if yes 
    261        vinfor(51) = zindb*vinfor(51) / MAX(vinfor(27),epsi06) 
    262         
    263        !! Fram Strait Export 
    264        !! 83 = area export 
    265        !! 84 = volume export 
    266        !! Fram strait in ORCA2 = 5 points 
    267        !! export = -v_ice*e1t*ddtb*at_i or -v_ice*e1t*ddtb*at_i*h_i 
    268        jj = 136 ! C grid 
    269        vinfor(83) = 0.0 
    270        vinfor(84) = 0.0 
    271        DO ji = 134, 138 
    272           vinfor(83) = vinfor(83) - v_ice(ji,jj) * &  
    273                                       e1t(ji,jj)*at_i(ji,jj)*rdt_ice / 1.0e12 
    274           vinfor(84) = vinfor(84) - v_ice(ji,jj) * &  
    275                                       e1t(ji,jj)*vt_i(ji,jj)*rdt_ice / 1.0e12 
    276        END DO 
    277  
    278        !!------------------------------------------------------------------- 
    279        !! 2) Southern hemisphere 
    280        !!------------------------------------------------------------------- 
    281        !! 2.1) Diagnostics independent on age 
    282        !!------------------------------------ 
    283        DO jj = 2, njeqm1 
    284           DO ji = fs_2, fs_jpim1   ! vector opt. 
    285              IF( tms(ji,jj) == 1 ) THEN 
    286                 vinfor(4)  = vinfor(4)  + at_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice area 
    287                 IF (at_i(ji,jj).GT.0.15) vinfor(6) = vinfor(6) + aire(ji,jj) / 1.0e12 !ice extent 
    288                 vinfor(8)  = vinfor(8)  + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume 
    289                 vinfor(10) = vinfor(10) + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume 
    290                 vinfor(16) = vinfor(16) + ot_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age 
    291                 vinfor(30) = vinfor(30) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity 
    292                 ! this diagnostic is not well computed (weighted by vol instead 
    293                 ! of area) 
    294                 vinfor(32) = vinfor(32) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + &  
    295                                                         v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel 
    296                 vinfor(54) = vinfor(54) + at_i(ji,jj)*emps(ji,jj)*aire(ji,jj) / 1.0e12 ! Total salt flux 
    297                 vinfor(56) = vinfor(56) + at_i(ji,jj)*fsbri(ji,jj)*aire(ji,jj) / 1.0e12 ! Brine drainage salt flux 
    298                 vinfor(58) = vinfor(58) + at_i(ji,jj)*fseqv(ji,jj)*aire(ji,jj) / 1.0e12 ! Equivalent salt flux 
    299                 vinfor(60) = vinfor(60) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SST 
    300                 vinfor(62) = vinfor(62) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SSS 
    301                 vinfor(66) = vinfor(66) + et_s(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! snow temperature 
    302                 vinfor(68) = vinfor(68) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! ice enthalpy 
    303                 vinfor(70) = vinfor(70) + v_i(ji,jj,1)*aire(ji,jj) / 1.0e12 !ice volume 
    304                 vinfor(72) = vinfor(72) + v_i(ji,jj,2)*aire(ji,jj) / 1.0e12 !ice volume 
    305                 vinfor(74) = vinfor(74) + v_i(ji,jj,3)*aire(ji,jj) / 1.0e12 !ice volume 
    306                 vinfor(76) = vinfor(76) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume 
    307                 vinfor(78) = vinfor(78) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume 
    308                 vinfor(80) = 0.0 
    309                 vinfor(82) = vinfor(82) + emp(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux 
    310              ENDIF 
    311           END DO 
    312        END DO 
    313  
    314        DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
    315           DO jj = 2, njeqm1 
    316              DO ji = fs_2, fs_jpim1   ! vector opt. 
    317                 vinfor(12) = vinfor(12) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !undef def ice volume 
    318              END DO 
    319           END DO 
    320        END DO 
    321  
    322        vinfor(14) = 0.0 
    323  
    324        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(8)))  
    325        vinfor(16) = zindb * vinfor(16) / MAX(vinfor(8),epsi06) ! these have to be divided by ice vol 
    326        vinfor(30) = zindb * vinfor(30) / MAX(vinfor(8),epsi06) !  
    327        vinfor(32) = zindb * SQRT( vinfor(32) / MAX( vinfor(8) , epsi06 ) ) 
    328        vinfor(68) = zindb * vinfor(68) / MAX(vinfor(8),epsi06) !  
    329  
    330        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(6)))  
    331        vinfor(54) = zindb * vinfor(54) / MAX(vinfor(6),epsi06) ! these have to be divided by ice extt 
    332        vinfor(56) = zindb * vinfor(56) / MAX(vinfor(6),epsi06) !  
    333        vinfor(58) = zindb * vinfor(58) / MAX(vinfor(6),epsi06) !  
    334        vinfor(80) = zindb * vinfor(80) / MAX(vinfor(6),epsi06) ! 
    335 !      vinfor(84) = vinfor(84) / vinfor(6) ! 
    336   
    337        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4))) ! 
    338        vinfor(60) = zindb*vinfor(60) / ( MAX(vinfor(4), epsi06) ) ! divide by ice area 
    339        vinfor(62) = zindb*vinfor(62) / ( MAX(vinfor(4), epsi06) ) ! 
    340  
    341        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(10))) ! 
    342        vinfor(66) = zindb*vinfor(66) / MAX(vinfor(10),epsi06) ! divide it by snow volume 
    343  
    344        DO jl = 1, jpl 
    345           DO jj = 2, njeqm1 
    346              DO ji = fs_2, fs_jpim1   ! vector opt. 
    347                 IF( tms(ji,jj) == 1 ) THEN 
    348                    vinfor(34) = vinfor(34) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume 
    349                    vinfor(36) = vinfor(36) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume 
    350                 ENDIF 
    351              END DO 
    352           END DO 
    353        END DO 
    354  
    355        DO jj = 2, njeqm1 
    356           DO ji = fs_2, fs_jpim1   ! vector opt. 
    357                 IF( tms(ji,jj) == 1 ) THEN 
    358                    vinfor(38) = vinfor(38) + diag_sni_gr(ji,jj)*aire(ji,jj) / 1.0e12 !th growth rates 
    359                    vinfor(40) = vinfor(40) + diag_lat_gr(ji,jj)*aire(ji,jj) / 1.0e12  
    360                    vinfor(42) = vinfor(42) + diag_bot_gr(ji,jj)*aire(ji,jj) / 1.0e12 
    361                    vinfor(44) = vinfor(44) + diag_dyn_gr(ji,jj)*aire(ji,jj) / 1.0e12  
    362                    vinfor(46) = vinfor(46) + dv_dt_thd(ji,jj,5)*aire(ji,jj) / 1.0e12 
    363                    vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) / 1.0e12 / rdt_ice ! volume acc in OW 
    364                 ENDIF 
    365           END DO 
    366        END DO 
    367  
    368  
    369        DO jl = 1, jpl 
    370           DO jj = 2, njeqm1 
    371              DO ji = fs_2, fs_jpim1   ! vector opt. 
    372                 IF( tms(ji,jj) == 1 ) THEN 
    373                    vinfor(64) = vinfor(64) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 
    374                 ENDIF 
    375              END DO 
    376           END DO 
    377        END DO 
    378        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4))) ! 
    379        vinfor(64) = zindb * vinfor(64) / MAX(vinfor(4),epsi06) ! divide by ice extt 
    380        !! 2.2) Diagnostics dependent on age 
    381        !!------------------------------------ 
    382        DO jj = 2, njeqm1 
    383           DO ji = fs_2, fs_jpim1   ! vector opt. 
    384              IF( tms(ji,jj) == 1 ) THEN 
    385                 zafy = 0.0 
    386                 zamy = 0.0 
    387                 DO jl = 1, jpl 
    388                    IF ((o_i(ji,jj,jl) - zshift_date).LT.0.0) THEN 
    389                       vinfor(18) = vinfor(18) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice area 
    390                       vinfor(26) = vinfor(26) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice volume 
    391                       zafy = zafy + a_i(ji,jj,jl) 
    392                       vinfor(50) = vinfor(50) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity 
    393                    ENDIF 
    394                    IF ((o_i(ji,jj,jl) - zshift_date).GT.0.0) THEN 
    395                       vinfor(20) = vinfor(20) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12    ! MY ice area 
    396                       vinfor(28) = vinfor(28) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 
    397                       vinfor(52) = vinfor(52) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity 
    398                       zamy = zamy + a_i(ji,jj,jl) 
    399                    ENDIF 
    400                 END DO ! jl 
    401                 IF ((at_i(ji,jj).GT.0.15).AND.(zafy.GT.zamy)) THEN 
    402                    vinfor(22) = vinfor(22) + aire(ji,jj) / 1.0e12 ! Seasonal ice extent 
    403                 ENDIF 
    404                 IF ((at_i(ji,jj).GT.0.15).AND.(zafy.LE.zamy)) THEN 
    405                    vinfor(24) = vinfor(24) + aire(ji,jj) / 1.0e12 ! Perennial ice extent 
    406                 ENDIF 
    407              ENDIF ! tms 
    408           END DO ! jj 
    409        END DO ! ji 
    410        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(26))) !=0 if no multiyear ice 1 if yes 
    411        vinfor(50) = zindb*vinfor(50) / MAX(vinfor(26),epsi06) 
    412        zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(28))) !=0 if no multiyear ice 1 if yes 
    413        vinfor(52) = zindb*vinfor(52) / MAX(vinfor(28),epsi06) 
    414  
    415        !  Accumulation before averaging  
    416        DO jv = 1, nvinfo 
    417           vinfom(jv) = vinfom(jv) + vinfor(jv) 
    418        END DO 
    419        naveg = naveg + 1   
    420      
    421        ! oututs on file ice_evolu     
    422 !MV      IF( MOD( numit , ninfo ) == 0 ) THEN 
    423           WRITE(numevo_ice,fmtw) ( titvar(jv), vinfom(jv)/naveg, jv = 1, nvinfo ) 
    424           naveg = 0 
    425           DO jv = 1, nvinfo 
    426              vinfom(jv)=0.0 
    427           END DO 
    428 !MV      ENDIF 
    429    
    430     END SUBROUTINE lim_dia 
    431   
    432     SUBROUTINE lim_dia_init 
    433        !!------------------------------------------------------------------- 
    434        !!                  ***  ROUTINE lim_dia_init  *** 
    435        !!              
    436        !! ** Purpose : Preparation of the file ice_evolu for the output of 
    437        !!      the temporal evolution of key variables 
    438        !! 
    439        !! ** input   : Namelist namicedia 
    440        !! 
    441        !! history : 
    442        !!  8.5  ! 03-08 (C. Ethe) original code 
    443        !!  9.0  ! 08-03 (M. Vancoppenolle) LIM3 
    444        !!------------------------------------------------------------------- 
    445        NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy 
    446  
    447        INTEGER  ::   jv   ,     &  ! dummy loop indice 
    448           &          ntot ,     & 
    449           &          ndeb ,     & 
    450           &          irecl 
    451  
    452        REAL(wp) ::   zxx0, zxx1    ! temporary scalars 
    453  
    454        CHARACTER(len=jpchinf) ::   titinf 
    455        CHARACTER(len=50)      ::   clname 
    456        !!------------------------------------------------------------------- 
    457  
    458  
    459        ! Read Namelist namicedia 
    460        REWIND ( numnam_ice ) 
    461        READ   ( numnam_ice  , namicedia ) 
    462        IF(lwp) THEN 
    463           WRITE(numout,*) 
    464           WRITE(numout,*) 'lim_dia_init : ice parameters for ice diagnostics ' 
    465           WRITE(numout,*) '~~~~~~~~~~~~' 
    466           WRITE(numout,*) '   format of the output values                                 fmtinf = ', fmtinf 
    467           WRITE(numout,*) '   number of variables written in one line                     nfrinf = ', nfrinf  
    468           WRITE(numout,*) '   Instantaneous values of ice evolution or averaging          ntmoy  = ', ntmoy 
    469           WRITE(numout,*) '   frequency of ouputs on file ice_evolu in case of averaging  ninfo  = ', ninfo 
    470        ENDIF 
    471  
    472        ! masked grid cell area 
    473        aire(:,:) = area(:,:) * tms(:,:) 
    474  
    475        ! Titles of ice key variables : 
    476        titvar(1) = 'NoIt'  ! iteration number 
    477        titvar(2) = 'T yr'  ! time step in years 
    478        nbvt = 2            ! number of time variables 
    479  
    480        titvar(3) = 'AI_N'  ! sea ice area in the northern Hemisp.(10^12 km2) 
    481        titvar(4) = 'AI_S'  ! sea ice area in the southern Hemisp.(10^12 km2) 
    482        titvar(5) = 'EI_N'  ! sea ice extent (15%) in the northern Hemisp.(10^12 km2) 
    483        titvar(6) = 'EI_S'  ! sea ice extent (15%) in the southern Hemisp.(10^12 km2) 
    484        titvar(7) = 'VI_N'  ! sea ice volume in the northern Hemisp.(10^3 km3) 
    485        titvar(8) = 'VI_S'  ! sea ice volume in the southern Hemisp.(10^3 km3) 
    486        titvar(9) = 'VS_N'  ! snow volume over sea ice in the northern Hemisp.(10^3 km3) 
    487        titvar(10)= 'VS_S'  ! snow volume over sea ice in the northern Hemisp.(10^3 km3) 
    488        titvar(11)= 'VuIN'  ! undeformed sea ice volume in the northern Hemisp.(10^3 km3) 
    489        titvar(12)= 'VuIS'  ! undeformed sea ice volume in the southern Hemisp.(10^3 km3) 
    490        titvar(13)= 'VdIN'  ! deformed sea ice volume in the northern Hemisp.(10^3 km3) 
    491        titvar(14)= 'VdIS'  ! deformed sea ice volume in the southern Hemisp.(10^3 km3) 
    492        titvar(15)= 'OI_N'  ! sea ice mean age in the northern Hemisp.(years) 
    493        titvar(16)= 'OI_S'  ! sea ice mean age in the southern Hemisp.(years) 
    494        titvar(17)= 'AFYN'  ! total FY ice area northern Hemisp.(10^12 km2) 
    495        titvar(18)= 'AFYS'  ! total FY ice area southern Hemisp.(10^12 km2) 
    496        titvar(19)= 'AMYN'  ! total MY ice area northern Hemisp.(10^12 km2) 
    497        titvar(20)= 'AMYS'  ! total MY ice area southern Hemisp.(10^12 km2) 
    498        titvar(21)= 'EFYN'  ! total FY ice extent northern Hemisp.(10^12 km2) (with more 50% FY ice) 
    499        titvar(22)= 'EFYS'  ! total FY ice extent southern Hemisp.(10^12 km2) (with more 50% FY ice) 
    500        titvar(23)= 'EMYN'  ! total MY ice extent northern Hemisp.(10^12 km2) (with more 50% MY ice) 
    501        titvar(24)= 'EMYS'  ! total MY ice extent southern Hemisp.(10^12 km2) (with more 50% MY ice) 
    502        titvar(25)= 'VFYN'  ! total undeformed FY ice volume northern Hemisp.(10^3 km3)  
    503        titvar(26)= 'VFYS'  ! total undeformed FY ice volume southern Hemisp.(10^3 km3) 
    504        titvar(27)= 'VMYN'  ! total undeformed MY ice volume northern Hemisp.(10^3 km3)  
    505        titvar(28)= 'VMYS'  ! total undeformed MY ice volume southern Hemisp.(10^3 km3)  
    506        titvar(29)= 'IS_N'  ! sea ice mean salinity in the northern hemisphere (ppt)   
    507        titvar(30)= 'IS_S'  ! sea ice mean salinity in the southern hemisphere (ppt)   
    508        titvar(31)= 'IVeN'  ! sea ice mean velocity in the northern hemisphere (m/s)  
    509        titvar(32)= 'IVeS'  ! sea ice mean velocity in the southern hemisphere (m/s)  
    510        titvar(33)= 'DVDN'  ! variation of sea ice volume due to dynamics in the northern hemisphere 
    511        titvar(34)= 'DVDS'  ! variation of sea ice volume due to dynamics in the southern hemisphere 
    512        titvar(35)= 'DVTN'  ! variation of sea ice volume due to thermo in the   northern hemisphere 
    513        titvar(36)= 'DVTS'  ! variation of sea ice volume due to thermo in the   southern hemisphere 
    514        titvar(37)= 'TG1N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 1   
    515        titvar(38)= 'TG1S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 1   
    516        titvar(39)= 'TG2N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 2   
    517        titvar(40)= 'TG2S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 2   
    518        titvar(41)= 'TG3N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 3   
    519        titvar(42)= 'TG3S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 3   
    520        titvar(43)= 'TG4N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 4   
    521        titvar(44)= 'TG4S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 4   
    522        titvar(45)= 'TG5N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 5   
    523        titvar(46)= 'TG5S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 5   
    524        titvar(47)= 'LA_N'  ! lateral accretion growth rate, northern hemisphere 
    525        titvar(48)= 'LA_S'  ! lateral accretion growth rate, southern hemisphere  
    526        titvar(49)= 'SF_N'  ! Salinity FY, NH  
    527        titvar(50)= 'SF_S'  ! Salinity FY, SH  
    528        titvar(51)= 'SF_N'  ! Salinity MY, NH  
    529        titvar(52)= 'SF_S'  ! Salinity MY, SH  
    530        titvar(53)= 'Fs_N'  ! Total salt flux NH 
    531        titvar(54)= 'Fs_S'  ! Total salt flux SH 
    532        titvar(55)= 'FsbN'  ! Salt - brine drainage flux NH 
    533        titvar(56)= 'FsbS'  ! Salt - brine drainage flux SH 
    534        titvar(57)= 'FseN'  ! Salt - Equivalent salt flux NH 
    535        titvar(58)= 'FseS'  ! Salt - Equivalent salt flux SH 
    536        titvar(59)= 'SSTN'  ! SST, NH 
    537        titvar(60)= 'SSTS'  ! SST, SH 
    538        titvar(61)= 'SSSN'  ! SSS, NH 
    539        titvar(62)= 'SSSS'  ! SSS, SH 
    540        titvar(63)= 'TsuN'  ! Tsu, NH 
    541        titvar(64)= 'TsuS'  ! Tsu, SH 
    542        titvar(65)= 'TsnN'  ! Tsn, NH 
    543        titvar(66)= 'TsnS'  ! Tsn, SH 
    544        titvar(67)= 'ei_N'  ! ei, NH 
    545        titvar(68)= 'ei_S'  ! ei, SH 
    546        titvar(69)= 'vi1N'  ! vi1, NH 
    547        titvar(70)= 'vi1S'  ! vi1, SH 
    548        titvar(71)= 'vi2N'  ! vi2, NH 
    549        titvar(72)= 'vi2S'  ! vi2, SH 
    550        titvar(73)= 'vi3N'  ! vi3, NH 
    551        titvar(74)= 'vi3S'  ! vi3, SH 
    552        titvar(75)= 'vi4N'  ! vi4, NH 
    553        titvar(76)= 'vi4S'  ! vi4, SH 
    554        titvar(77)= 'vi5N'  ! vi5, NH 
    555        titvar(78)= 'vi5S'  ! vi5, SH 
    556        titvar(79)= 'vi6N'  ! vi6, NH 
    557        titvar(80)= 'vi6S'  ! vi6, SH 
    558        titvar(81)= 'fmaN'  ! mass flux in the ocean, NH 
    559        titvar(82)= 'fmaS'  ! mass flux in the ocean, SH 
    560        titvar(83)= 'AFSE'  ! Fram Strait Area export 
    561        titvar(84)= 'VFSE'  ! Fram Strait Volume export 
    562        nvinfo = 84 
    563  
    564        ! Definition et Ecriture de l'entete : nombre d'enregistrements  
    565        ndeb   = ( nstart - 1 ) / ninfo 
    566        IF( nstart == 1 ) ndeb = -1 
    567  
    568        nferme = ( nstart - 1 + nitrun) / ninfo 
    569        ntot   = nferme - ndeb 
    570        ndeb   = ninfo * ( 1 + ndeb ) 
    571        nferme = ninfo * nferme 
    572  
    573        ! definition of formats  
    574        WRITE( fmtw  , '(A,I3,A2,I1,A)' )  '(', nfrinf, '(A', jpchsep, ','//fmtinf//'))' 
    575        WRITE( fmtr  , '(A,I3,A,I1,A)'  )  '(', nfrinf, '(', jpchsep, 'X,'//fmtinf//'))' 
    576        WRITE( fmtitr, '(A,I3,A,I1,A)'  )  '(', nvinfo, 'A', jpchinf, ')' 
    577  
    578        ! opening  "ice_evolu" file 
    579        clname = 'ice.evolu' 
    580        irecl = ( jpchinf + 1 ) * nvinfo  
    581        CALL ctlopn( numevo_ice, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',    & 
    582           &         irecl, numout, lwp, 1 ) 
    583  
    584        !- ecriture de 2 lignes d''entete : 
    585        WRITE(numevo_ice,1000) fmtr, fmtw, fmtitr, nvinfo, ntot, 0, nfrinf 
    586        zxx0 = 0.001 * REAL(ninfo) 
    587        zxx1 = 0.001 * REAL(ndeb) 
    588        WRITE(numevo_ice,1111) REAL(jpchinf), 0., zxx1, zxx0, 0., 0., 0 
    589  
    590        !- ecriture de 2 lignes de titre : 
    591        WRITE(numevo_ice,'(A,I8,A,I8,A,I5)')                                      & 
    592           'Evolution chronologique - Experience '//cexper   & 
    593           //'   de', ndeb, ' a', nferme, ' pas', ninfo 
    594        WRITE(numevo_ice,fmtitr) ( titvar(jv), jv = 1, nvinfo ) 
    595  
    596  
    597        !--preparation de "titvar" pour l''ecriture parmi les valeurs numeriques : 
    598        DO  jv = 2 , nvinfo 
    599           titinf     = titvar(jv)(:jpchinf) 
    600           titvar(jv) = '  '//titinf 
    601        END DO 
    602  
    603        !--Initialisation of the arrays for the accumulation 
    604        DO  jv = 1, nvinfo 
    605           vinfom(jv) = 0. 
    606        END DO 
    607        naveg = 0 
    608  
    609 1000   FORMAT( 3(A20),4(1x,I6) ) 
    610 1111   FORMAT( 3(F7.1,1X,F7.3,1X),I3,A )   
    611  
    612     END SUBROUTINE lim_dia_init 
     97      INTEGER  ::   jv,ji,jj,jl ! dummy loop indices 
     98      REAL(wp), DIMENSION(jpinfmx) ::  &  
     99         vinfor           ! temporary working space  
     100      REAL(wp) ::    & 
     101         zshift_date   , & ! date from the minimum ice extent 
     102         zday, zday_min, & ! current day, day of minimum extent 
     103         zafy, zamy,     & ! temporary area of fy and my ice 
     104         zindb 
     105      !!------------------------------------------------------------------- 
     106 
     107      ! 0) date from the minimum of ice extent 
     108      !--------------------------------------- 
     109      zday_min = 273.0        ! zday_min = date of minimum extent, here September 30th 
     110      zday = FLOAT(numit-nit000) * rdt_ice / ( 86400.0 * FLOAT(nn_fsbc) ) 
     111      IF (zday.GT.zday_min) THEN  
     112         zshift_date  =  zday - zday_min 
     113      ELSE 
     114         zshift_date  =  zday - (365.0 - zday_min) 
     115      ENDIF 
     116 
     117      IF( numit == nstart )   CALL lim_dia_init   ! initialisation of ice_evolu file       
     118 
     119      ! temporal diagnostics  
     120      vinfor(1) = REAL(numit) 
     121      vinfor(2) = nyear 
     122 
     123      ! put everything to zero 
     124      DO jv = nbvt + 1, nvinfo 
     125         vinfor(jv) = 0.0 
     126      END DO 
     127 
     128      !!------------------------------------------------------------------- 
     129      !! 1) Northern hemisphere 
     130      !!------------------------------------------------------------------- 
     131      !! 1.1) Diagnostics independent on age 
     132      !!------------------------------------ 
     133      DO jj = njeq, jpjm1 
     134         DO ji = fs_2, fs_jpim1   ! vector opt. 
     135            IF( tms(ji,jj) == 1 ) THEN 
     136               vinfor(3)  = vinfor(3)  + at_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice area 
     137               IF (at_i(ji,jj).GT.0.15) vinfor(5) = vinfor(5) + aire(ji,jj) / 1.0e12 !ice extent 
     138               vinfor(7)  = vinfor(7)  + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume 
     139               vinfor(9)  = vinfor(9)  + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume 
     140               vinfor(15) = vinfor(15) + ot_i(ji,jj) *vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age 
     141               vinfor(29) = vinfor(29) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity 
     142               ! the computation of this diagnostic is not reliable 
     143               vinfor(31) = vinfor(31) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + &  
     144                  v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12  
     145               vinfor(53) = vinfor(53) + emps(ji,jj)*aire(ji,jj) / 1.0e12 !salt flux 
     146               vinfor(55) = vinfor(55) + fsbri(ji,jj)*aire(ji,jj) / 1.0e12 !brine drainage flux 
     147               vinfor(57) = vinfor(57) + fseqv(ji,jj)*aire(ji,jj) / 1.0e12 !equivalent salt flux 
     148               vinfor(59) = vinfor(59) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SST 
     149               vinfor(61) = vinfor(61) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SSS 
     150               vinfor(65) = vinfor(65) + et_s(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12  ! snow temperature 
     151               vinfor(67) = vinfor(67) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12       ! ice heat content 
     152               vinfor(69) = vinfor(69) + v_i(ji,jj,1)*aire(ji,jj) / 1.0e12 !ice volume 
     153               vinfor(71) = vinfor(71) + v_i(ji,jj,2)*aire(ji,jj) / 1.0e12 !ice volume 
     154               vinfor(73) = vinfor(73) + v_i(ji,jj,3)*aire(ji,jj) / 1.0e12 !ice volume 
     155               vinfor(75) = vinfor(75) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume 
     156               vinfor(77) = vinfor(77) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume 
     157               vinfor(79) = 0.0 
     158               vinfor(81) = vinfor(81) + emp(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux 
     159            ENDIF 
     160         END DO 
     161      END DO 
     162 
     163      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
     164         DO jj = njeq, jpjm1 
     165            DO ji = fs_2, fs_jpim1   ! vector opt. 
     166               IF( tms(ji,jj) == 1 ) THEN 
     167                  vinfor(11) = vinfor(11) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !undef def ice volume 
     168               ENDIF 
     169            END DO 
     170         END DO 
     171      END DO 
     172 
     173      vinfor(13) = 0.0 
     174 
     175      vinfor(15) = vinfor(15) / MAX(vinfor(7),epsi06) ! these have to be divided by total ice volume to have the 
     176      vinfor(29) = vinfor(29) / MAX(vinfor(7),epsi06) ! right value 
     177      vinfor(31) = SQRT( vinfor(31) / MAX( vinfor(7) , epsi06 ) ) 
     178      vinfor(67) = vinfor(67) / MAX(vinfor(7),epsi06) 
     179 
     180      vinfor(53) = vinfor(53) / MAX(vinfor(5),epsi06) ! these have to be divided by total ice extent to have the 
     181      vinfor(55) = vinfor(55) / MAX(vinfor(5),epsi06) ! right value  
     182      vinfor(57) = vinfor(57) / MAX(vinfor(5),epsi06) !  
     183      vinfor(79) = vinfor(79) / MAX(vinfor(5),epsi06) ! 
     184 
     185      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(3))) ! 
     186      vinfor(59) = zindb*vinfor(59) / MAX(vinfor(3),epsi06) ! divide by ice area 
     187      vinfor(61) = zindb*vinfor(61) / MAX(vinfor(3),epsi06) ! 
     188 
     189      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(9))) ! 
     190      vinfor(65) = zindb*vinfor(65) / MAX(vinfor(9),epsi06) ! divide it by snow volume 
     191 
     192 
     193      DO jl = 1, jpl 
     194         DO jj = njeq, jpjm1 
     195            DO ji = fs_2, fs_jpim1   ! vector opt. 
     196               IF( tms(ji,jj) == 1 ) THEN 
     197                  vinfor(33) = vinfor(33) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume 
     198                  vinfor(35) = vinfor(35) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume 
     199               ENDIF 
     200            END DO 
     201         END DO 
     202      END DO 
     203 
     204      DO jj = njeq, jpjm1 
     205         DO ji = fs_2, fs_jpim1   ! vector opt. 
     206            IF( tms(ji,jj) == 1 ) THEN 
     207               vinfor(37) = vinfor(37) + diag_sni_gr(ji,jj)*aire(ji,jj) / 1.0e12 !th growth rates 
     208               vinfor(39) = vinfor(39) + diag_lat_gr(ji,jj)*aire(ji,jj) / 1.0e12  
     209               vinfor(41) = vinfor(41) + diag_bot_gr(ji,jj)*aire(ji,jj) / 1.0e12 
     210               vinfor(43) = vinfor(43) + diag_dyn_gr(ji,jj)*aire(ji,jj) / 1.0e12  
     211               vinfor(45) = vinfor(45) + dv_dt_thd(ji,jj,5)*aire(ji,jj) / 1.0e12 
     212               vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) / 1.0e12 / rdt_ice ! volume acc in OW 
     213            ENDIF 
     214         END DO 
     215      END DO 
     216 
     217      DO jl = 1, jpl 
     218         DO jj = njeq, jpjm1 
     219            DO ji = fs_2, fs_jpim1   ! vector opt. 
     220               IF( tms(ji,jj) == 1 ) THEN 
     221                  vinfor(63) = vinfor(63) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 
     222               ENDIF 
     223            END DO 
     224         END DO 
     225      END DO 
     226      vinfor(63) = vinfor(63) / MAX(vinfor(3),epsi06) ! these have to be divided by total ice area 
     227 
     228      !! 1.2) Diagnostics dependent on age 
     229      !!------------------------------------ 
     230      DO jj = njeq, jpjm1 
     231         DO ji = fs_2, fs_jpim1   ! vector opt. 
     232            IF( tms(ji,jj) == 1 ) THEN 
     233               zafy = 0.0 
     234               zamy = 0.0 
     235               DO jl = 1, jpl 
     236                  IF ((o_i(ji,jj,jl) - zshift_date).LT.0.0) THEN 
     237                     vinfor(17) = vinfor(17) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice area 
     238                     vinfor(25) = vinfor(25) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice volume 
     239                     vinfor(49) = vinfor(49) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity 
     240                     zafy = zafy + a_i(ji,jj,jl) 
     241                  ENDIF 
     242                  IF ((o_i(ji,jj,jl) - zshift_date).GT.0.0) THEN 
     243                     vinfor(19) = vinfor(19) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12    ! MY ice area 
     244                     vinfor(27) = vinfor(27) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! MY ice volume 
     245                     vinfor(51) = vinfor(51) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !MY ice salinity 
     246                     zamy = zamy + a_i(ji,jj,jl) 
     247                  ENDIF 
     248               END DO 
     249               IF ((at_i(ji,jj).GT.0.15).AND.(zafy.GT.zamy)) THEN 
     250                  vinfor(21) = vinfor(21) + aire(ji,jj) / 1.0e12 ! Seasonal ice extent 
     251               ENDIF 
     252               IF ((at_i(ji,jj).GT.0.15).AND.(zafy.LE.zamy)) THEN 
     253                  vinfor(23) = vinfor(23) + aire(ji,jj) / 1.0e12 ! Perennial ice extent 
     254               ENDIF 
     255            ENDIF 
     256         END DO 
     257      END DO 
     258      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(25))) !=0 if no multiyear ice 1 if yes 
     259      vinfor(49) = zindb*vinfor(49) / MAX(vinfor(25),epsi06) 
     260      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(27))) !=0 if no multiyear ice 1 if yes 
     261      vinfor(51) = zindb*vinfor(51) / MAX(vinfor(27),epsi06) 
     262 
     263      !! Fram Strait Export 
     264      !! 83 = area export 
     265      !! 84 = volume export 
     266      !! Fram strait in ORCA2 = 5 points 
     267      !! export = -v_ice*e1t*ddtb*at_i or -v_ice*e1t*ddtb*at_i*h_i 
     268      jj = 136 ! C grid 
     269      vinfor(83) = 0.0 
     270      vinfor(84) = 0.0 
     271      DO ji = 134, 138 
     272         vinfor(83) = vinfor(83) - v_ice(ji,jj) * &  
     273            e1t(ji,jj)*at_i(ji,jj)*rdt_ice / 1.0e12 
     274         vinfor(84) = vinfor(84) - v_ice(ji,jj) * &  
     275            e1t(ji,jj)*vt_i(ji,jj)*rdt_ice / 1.0e12 
     276      END DO 
     277 
     278      !!------------------------------------------------------------------- 
     279      !! 2) Southern hemisphere 
     280      !!------------------------------------------------------------------- 
     281      !! 2.1) Diagnostics independent on age 
     282      !!------------------------------------ 
     283      DO jj = 2, njeqm1 
     284         DO ji = fs_2, fs_jpim1   ! vector opt. 
     285            IF( tms(ji,jj) == 1 ) THEN 
     286               vinfor(4)  = vinfor(4)  + at_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice area 
     287               IF (at_i(ji,jj).GT.0.15) vinfor(6) = vinfor(6) + aire(ji,jj) / 1.0e12 !ice extent 
     288               vinfor(8)  = vinfor(8)  + vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !ice volume 
     289               vinfor(10) = vinfor(10) + vt_s(ji,jj)*aire(ji,jj) / 1.0e12 !snow volume 
     290               vinfor(16) = vinfor(16) + ot_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean age 
     291               vinfor(30) = vinfor(30) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) / 1.0e12 !mean salinity 
     292               ! this diagnostic is not well computed (weighted by vol instead 
     293               ! of area) 
     294               vinfor(32) = vinfor(32) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + &  
     295                  v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel 
     296               vinfor(54) = vinfor(54) + at_i(ji,jj)*emps(ji,jj)*aire(ji,jj) / 1.0e12 ! Total salt flux 
     297               vinfor(56) = vinfor(56) + at_i(ji,jj)*fsbri(ji,jj)*aire(ji,jj) / 1.0e12 ! Brine drainage salt flux 
     298               vinfor(58) = vinfor(58) + at_i(ji,jj)*fseqv(ji,jj)*aire(ji,jj) / 1.0e12 ! Equivalent salt flux 
     299               vinfor(60) = vinfor(60) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SST 
     300               vinfor(62) = vinfor(62) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) / 1.0e12  !SSS 
     301               vinfor(66) = vinfor(66) + et_s(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! snow temperature 
     302               vinfor(68) = vinfor(68) + et_i(ji,jj)/1.0e9*aire(ji,jj) / 1.0e12 ! ice enthalpy 
     303               vinfor(70) = vinfor(70) + v_i(ji,jj,1)*aire(ji,jj) / 1.0e12 !ice volume 
     304               vinfor(72) = vinfor(72) + v_i(ji,jj,2)*aire(ji,jj) / 1.0e12 !ice volume 
     305               vinfor(74) = vinfor(74) + v_i(ji,jj,3)*aire(ji,jj) / 1.0e12 !ice volume 
     306               vinfor(76) = vinfor(76) + v_i(ji,jj,4)*aire(ji,jj) / 1.0e12 !ice volume 
     307               vinfor(78) = vinfor(78) + v_i(ji,jj,5)*aire(ji,jj) / 1.0e12 !ice volume 
     308               vinfor(80) = 0.0 
     309               vinfor(82) = vinfor(82) + emp(ji,jj)*aire(ji,jj) / 1.0e12 ! mass flux 
     310            ENDIF 
     311         END DO 
     312      END DO 
     313 
     314      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
     315         DO jj = 2, njeqm1 
     316            DO ji = fs_2, fs_jpim1   ! vector opt. 
     317               vinfor(12) = vinfor(12) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !undef def ice volume 
     318            END DO 
     319         END DO 
     320      END DO 
     321 
     322      vinfor(14) = 0.0 
     323 
     324      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(8)))  
     325      vinfor(16) = zindb * vinfor(16) / MAX(vinfor(8),epsi06) ! these have to be divided by ice vol 
     326      vinfor(30) = zindb * vinfor(30) / MAX(vinfor(8),epsi06) !  
     327      vinfor(32) = zindb * SQRT( vinfor(32) / MAX( vinfor(8) , epsi06 ) ) 
     328      vinfor(68) = zindb * vinfor(68) / MAX(vinfor(8),epsi06) !  
     329 
     330      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(6)))  
     331      vinfor(54) = zindb * vinfor(54) / MAX(vinfor(6),epsi06) ! these have to be divided by ice extt 
     332      vinfor(56) = zindb * vinfor(56) / MAX(vinfor(6),epsi06) !  
     333      vinfor(58) = zindb * vinfor(58) / MAX(vinfor(6),epsi06) !  
     334      vinfor(80) = zindb * vinfor(80) / MAX(vinfor(6),epsi06) ! 
     335      !      vinfor(84) = vinfor(84) / vinfor(6) ! 
     336 
     337      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4))) ! 
     338      vinfor(60) = zindb*vinfor(60) / ( MAX(vinfor(4), epsi06) ) ! divide by ice area 
     339      vinfor(62) = zindb*vinfor(62) / ( MAX(vinfor(4), epsi06) ) ! 
     340 
     341      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(10))) ! 
     342      vinfor(66) = zindb*vinfor(66) / MAX(vinfor(10),epsi06) ! divide it by snow volume 
     343 
     344      DO jl = 1, jpl 
     345         DO jj = 2, njeqm1 
     346            DO ji = fs_2, fs_jpim1   ! vector opt. 
     347               IF( tms(ji,jj) == 1 ) THEN 
     348                  vinfor(34) = vinfor(34) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume 
     349                  vinfor(36) = vinfor(36) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) / 1.0e12 !ice volume 
     350               ENDIF 
     351            END DO 
     352         END DO 
     353      END DO 
     354 
     355      DO jj = 2, njeqm1 
     356         DO ji = fs_2, fs_jpim1   ! vector opt. 
     357            IF( tms(ji,jj) == 1 ) THEN 
     358               vinfor(38) = vinfor(38) + diag_sni_gr(ji,jj)*aire(ji,jj) / 1.0e12 !th growth rates 
     359               vinfor(40) = vinfor(40) + diag_lat_gr(ji,jj)*aire(ji,jj) / 1.0e12  
     360               vinfor(42) = vinfor(42) + diag_bot_gr(ji,jj)*aire(ji,jj) / 1.0e12 
     361               vinfor(44) = vinfor(44) + diag_dyn_gr(ji,jj)*aire(ji,jj) / 1.0e12  
     362               vinfor(46) = vinfor(46) + dv_dt_thd(ji,jj,5)*aire(ji,jj) / 1.0e12 
     363               vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) / 1.0e12 / rdt_ice ! volume acc in OW 
     364            ENDIF 
     365         END DO 
     366      END DO 
     367 
     368 
     369      DO jl = 1, jpl 
     370         DO jj = 2, njeqm1 
     371            DO ji = fs_2, fs_jpim1   ! vector opt. 
     372               IF( tms(ji,jj) == 1 ) THEN 
     373                  vinfor(64) = vinfor(64) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 
     374               ENDIF 
     375            END DO 
     376         END DO 
     377      END DO 
     378      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4))) ! 
     379      vinfor(64) = zindb * vinfor(64) / MAX(vinfor(4),epsi06) ! divide by ice extt 
     380      !! 2.2) Diagnostics dependent on age 
     381      !!------------------------------------ 
     382      DO jj = 2, njeqm1 
     383         DO ji = fs_2, fs_jpim1   ! vector opt. 
     384            IF( tms(ji,jj) == 1 ) THEN 
     385               zafy = 0.0 
     386               zamy = 0.0 
     387               DO jl = 1, jpl 
     388                  IF ((o_i(ji,jj,jl) - zshift_date).LT.0.0) THEN 
     389                     vinfor(18) = vinfor(18) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice area 
     390                     vinfor(26) = vinfor(26) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 ! FY ice volume 
     391                     zafy = zafy + a_i(ji,jj,jl) 
     392                     vinfor(50) = vinfor(50) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity 
     393                  ENDIF 
     394                  IF ((o_i(ji,jj,jl) - zshift_date).GT.0.0) THEN 
     395                     vinfor(20) = vinfor(20) + a_i(ji,jj,jl)*aire(ji,jj) / 1.0e12    ! MY ice area 
     396                     vinfor(28) = vinfor(28) + v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 
     397                     vinfor(52) = vinfor(52) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) / 1.0e12 !FY ice salinity 
     398                     zamy = zamy + a_i(ji,jj,jl) 
     399                  ENDIF 
     400               END DO ! jl 
     401               IF ((at_i(ji,jj).GT.0.15).AND.(zafy.GT.zamy)) THEN 
     402                  vinfor(22) = vinfor(22) + aire(ji,jj) / 1.0e12 ! Seasonal ice extent 
     403               ENDIF 
     404               IF ((at_i(ji,jj).GT.0.15).AND.(zafy.LE.zamy)) THEN 
     405                  vinfor(24) = vinfor(24) + aire(ji,jj) / 1.0e12 ! Perennial ice extent 
     406               ENDIF 
     407            ENDIF ! tms 
     408         END DO ! jj 
     409      END DO ! ji 
     410      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(26))) !=0 if no multiyear ice 1 if yes 
     411      vinfor(50) = zindb*vinfor(50) / MAX(vinfor(26),epsi06) 
     412      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(28))) !=0 if no multiyear ice 1 if yes 
     413      vinfor(52) = zindb*vinfor(52) / MAX(vinfor(28),epsi06) 
     414 
     415      !  Accumulation before averaging  
     416      DO jv = 1, nvinfo 
     417         vinfom(jv) = vinfom(jv) + vinfor(jv) 
     418      END DO 
     419      naveg = naveg + 1   
     420 
     421      ! oututs on file ice_evolu     
     422      !MV      IF( MOD( numit , ninfo ) == 0 ) THEN 
     423      WRITE(numevo_ice,fmtw) ( titvar(jv), vinfom(jv)/naveg, jv = 1, nvinfo ) 
     424      naveg = 0 
     425      DO jv = 1, nvinfo 
     426         vinfom(jv)=0.0 
     427      END DO 
     428      !MV      ENDIF 
     429 
     430   END SUBROUTINE lim_dia 
     431 
     432   SUBROUTINE lim_dia_init 
     433      !!------------------------------------------------------------------- 
     434      !!                  ***  ROUTINE lim_dia_init  *** 
     435      !!              
     436      !! ** Purpose : Preparation of the file ice_evolu for the output of 
     437      !!      the temporal evolution of key variables 
     438      !! 
     439      !! ** input   : Namelist namicedia 
     440      !! 
     441      !! history : 
     442      !!  8.5  ! 03-08 (C. Ethe) original code 
     443      !!  9.0  ! 08-03 (M. Vancoppenolle) LIM3 
     444      !!------------------------------------------------------------------- 
     445      NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy 
     446 
     447      INTEGER  ::   jv   ,     &  ! dummy loop indice 
     448         &          ntot ,     & 
     449         &          ndeb ,     & 
     450         &          irecl 
     451 
     452      REAL(wp) ::   zxx0, zxx1    ! temporary scalars 
     453 
     454      CHARACTER(len=jpchinf) ::   titinf 
     455      CHARACTER(len=50)      ::   clname 
     456      !!------------------------------------------------------------------- 
     457 
     458 
     459      ! Read Namelist namicedia 
     460      REWIND ( numnam_ice ) 
     461      READ   ( numnam_ice  , namicedia ) 
     462      IF(lwp) THEN 
     463         WRITE(numout,*) 
     464         WRITE(numout,*) 'lim_dia_init : ice parameters for ice diagnostics ' 
     465         WRITE(numout,*) '~~~~~~~~~~~~' 
     466         WRITE(numout,*) '   format of the output values                                 fmtinf = ', fmtinf 
     467         WRITE(numout,*) '   number of variables written in one line                     nfrinf = ', nfrinf  
     468         WRITE(numout,*) '   Instantaneous values of ice evolution or averaging          ntmoy  = ', ntmoy 
     469         WRITE(numout,*) '   frequency of ouputs on file ice_evolu in case of averaging  ninfo  = ', ninfo 
     470      ENDIF 
     471 
     472      ! masked grid cell area 
     473      aire(:,:) = area(:,:) * tms(:,:) 
     474 
     475      ! Titles of ice key variables : 
     476      titvar(1) = 'NoIt'  ! iteration number 
     477      titvar(2) = 'T yr'  ! time step in years 
     478      nbvt = 2            ! number of time variables 
     479 
     480      titvar(3) = 'AI_N'  ! sea ice area in the northern Hemisp.(10^12 km2) 
     481      titvar(4) = 'AI_S'  ! sea ice area in the southern Hemisp.(10^12 km2) 
     482      titvar(5) = 'EI_N'  ! sea ice extent (15%) in the northern Hemisp.(10^12 km2) 
     483      titvar(6) = 'EI_S'  ! sea ice extent (15%) in the southern Hemisp.(10^12 km2) 
     484      titvar(7) = 'VI_N'  ! sea ice volume in the northern Hemisp.(10^3 km3) 
     485      titvar(8) = 'VI_S'  ! sea ice volume in the southern Hemisp.(10^3 km3) 
     486      titvar(9) = 'VS_N'  ! snow volume over sea ice in the northern Hemisp.(10^3 km3) 
     487      titvar(10)= 'VS_S'  ! snow volume over sea ice in the northern Hemisp.(10^3 km3) 
     488      titvar(11)= 'VuIN'  ! undeformed sea ice volume in the northern Hemisp.(10^3 km3) 
     489      titvar(12)= 'VuIS'  ! undeformed sea ice volume in the southern Hemisp.(10^3 km3) 
     490      titvar(13)= 'VdIN'  ! deformed sea ice volume in the northern Hemisp.(10^3 km3) 
     491      titvar(14)= 'VdIS'  ! deformed sea ice volume in the southern Hemisp.(10^3 km3) 
     492      titvar(15)= 'OI_N'  ! sea ice mean age in the northern Hemisp.(years) 
     493      titvar(16)= 'OI_S'  ! sea ice mean age in the southern Hemisp.(years) 
     494      titvar(17)= 'AFYN'  ! total FY ice area northern Hemisp.(10^12 km2) 
     495      titvar(18)= 'AFYS'  ! total FY ice area southern Hemisp.(10^12 km2) 
     496      titvar(19)= 'AMYN'  ! total MY ice area northern Hemisp.(10^12 km2) 
     497      titvar(20)= 'AMYS'  ! total MY ice area southern Hemisp.(10^12 km2) 
     498      titvar(21)= 'EFYN'  ! total FY ice extent northern Hemisp.(10^12 km2) (with more 50% FY ice) 
     499      titvar(22)= 'EFYS'  ! total FY ice extent southern Hemisp.(10^12 km2) (with more 50% FY ice) 
     500      titvar(23)= 'EMYN'  ! total MY ice extent northern Hemisp.(10^12 km2) (with more 50% MY ice) 
     501      titvar(24)= 'EMYS'  ! total MY ice extent southern Hemisp.(10^12 km2) (with more 50% MY ice) 
     502      titvar(25)= 'VFYN'  ! total undeformed FY ice volume northern Hemisp.(10^3 km3)  
     503      titvar(26)= 'VFYS'  ! total undeformed FY ice volume southern Hemisp.(10^3 km3) 
     504      titvar(27)= 'VMYN'  ! total undeformed MY ice volume northern Hemisp.(10^3 km3)  
     505      titvar(28)= 'VMYS'  ! total undeformed MY ice volume southern Hemisp.(10^3 km3)  
     506      titvar(29)= 'IS_N'  ! sea ice mean salinity in the northern hemisphere (ppt)   
     507      titvar(30)= 'IS_S'  ! sea ice mean salinity in the southern hemisphere (ppt)   
     508      titvar(31)= 'IVeN'  ! sea ice mean velocity in the northern hemisphere (m/s)  
     509      titvar(32)= 'IVeS'  ! sea ice mean velocity in the southern hemisphere (m/s)  
     510      titvar(33)= 'DVDN'  ! variation of sea ice volume due to dynamics in the northern hemisphere 
     511      titvar(34)= 'DVDS'  ! variation of sea ice volume due to dynamics in the southern hemisphere 
     512      titvar(35)= 'DVTN'  ! variation of sea ice volume due to thermo in the   northern hemisphere 
     513      titvar(36)= 'DVTS'  ! variation of sea ice volume due to thermo in the   southern hemisphere 
     514      titvar(37)= 'TG1N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 1   
     515      titvar(38)= 'TG1S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 1   
     516      titvar(39)= 'TG2N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 2   
     517      titvar(40)= 'TG2S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 2   
     518      titvar(41)= 'TG3N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 3   
     519      titvar(42)= 'TG3S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 3   
     520      titvar(43)= 'TG4N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 4   
     521      titvar(44)= 'TG4S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 4   
     522      titvar(45)= 'TG5N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 5   
     523      titvar(46)= 'TG5S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 5   
     524      titvar(47)= 'LA_N'  ! lateral accretion growth rate, northern hemisphere 
     525      titvar(48)= 'LA_S'  ! lateral accretion growth rate, southern hemisphere  
     526      titvar(49)= 'SF_N'  ! Salinity FY, NH  
     527      titvar(50)= 'SF_S'  ! Salinity FY, SH  
     528      titvar(51)= 'SF_N'  ! Salinity MY, NH  
     529      titvar(52)= 'SF_S'  ! Salinity MY, SH  
     530      titvar(53)= 'Fs_N'  ! Total salt flux NH 
     531      titvar(54)= 'Fs_S'  ! Total salt flux SH 
     532      titvar(55)= 'FsbN'  ! Salt - brine drainage flux NH 
     533      titvar(56)= 'FsbS'  ! Salt - brine drainage flux SH 
     534      titvar(57)= 'FseN'  ! Salt - Equivalent salt flux NH 
     535      titvar(58)= 'FseS'  ! Salt - Equivalent salt flux SH 
     536      titvar(59)= 'SSTN'  ! SST, NH 
     537      titvar(60)= 'SSTS'  ! SST, SH 
     538      titvar(61)= 'SSSN'  ! SSS, NH 
     539      titvar(62)= 'SSSS'  ! SSS, SH 
     540      titvar(63)= 'TsuN'  ! Tsu, NH 
     541      titvar(64)= 'TsuS'  ! Tsu, SH 
     542      titvar(65)= 'TsnN'  ! Tsn, NH 
     543      titvar(66)= 'TsnS'  ! Tsn, SH 
     544      titvar(67)= 'ei_N'  ! ei, NH 
     545      titvar(68)= 'ei_S'  ! ei, SH 
     546      titvar(69)= 'vi1N'  ! vi1, NH 
     547      titvar(70)= 'vi1S'  ! vi1, SH 
     548      titvar(71)= 'vi2N'  ! vi2, NH 
     549      titvar(72)= 'vi2S'  ! vi2, SH 
     550      titvar(73)= 'vi3N'  ! vi3, NH 
     551      titvar(74)= 'vi3S'  ! vi3, SH 
     552      titvar(75)= 'vi4N'  ! vi4, NH 
     553      titvar(76)= 'vi4S'  ! vi4, SH 
     554      titvar(77)= 'vi5N'  ! vi5, NH 
     555      titvar(78)= 'vi5S'  ! vi5, SH 
     556      titvar(79)= 'vi6N'  ! vi6, NH 
     557      titvar(80)= 'vi6S'  ! vi6, SH 
     558      titvar(81)= 'fmaN'  ! mass flux in the ocean, NH 
     559      titvar(82)= 'fmaS'  ! mass flux in the ocean, SH 
     560      titvar(83)= 'AFSE'  ! Fram Strait Area export 
     561      titvar(84)= 'VFSE'  ! Fram Strait Volume export 
     562      nvinfo = 84 
     563 
     564      ! Definition et Ecriture de l'entete : nombre d'enregistrements  
     565      ndeb   = ( nstart - 1 ) / ninfo 
     566      IF( nstart == 1 ) ndeb = -1 
     567 
     568      nferme = ( nstart - 1 + nitrun) / ninfo 
     569      ntot   = nferme - ndeb 
     570      ndeb   = ninfo * ( 1 + ndeb ) 
     571      nferme = ninfo * nferme 
     572 
     573      ! definition of formats  
     574      WRITE( fmtw  , '(A,I3,A2,I1,A)' )  '(', nfrinf, '(A', jpchsep, ','//fmtinf//'))' 
     575      WRITE( fmtr  , '(A,I3,A,I1,A)'  )  '(', nfrinf, '(', jpchsep, 'X,'//fmtinf//'))' 
     576      WRITE( fmtitr, '(A,I3,A,I1,A)'  )  '(', nvinfo, 'A', jpchinf, ')' 
     577 
     578      ! opening  "ice_evolu" file 
     579      clname = 'ice.evolu' 
     580      irecl = ( jpchinf + 1 ) * nvinfo  
     581      CALL ctlopn( numevo_ice, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',    & 
     582         &         irecl, numout, lwp, 1 ) 
     583 
     584      !- ecriture de 2 lignes d''entete : 
     585      WRITE(numevo_ice,1000) fmtr, fmtw, fmtitr, nvinfo, ntot, 0, nfrinf 
     586      zxx0 = 0.001 * REAL(ninfo) 
     587      zxx1 = 0.001 * REAL(ndeb) 
     588      WRITE(numevo_ice,1111) REAL(jpchinf), 0., zxx1, zxx0, 0., 0., 0 
     589 
     590      !- ecriture de 2 lignes de titre : 
     591      WRITE(numevo_ice,'(A,I8,A,I8,A,I5)')                                      & 
     592         'Evolution chronologique - Experience '//cexper   & 
     593         //'   de', ndeb, ' a', nferme, ' pas', ninfo 
     594      WRITE(numevo_ice,fmtitr) ( titvar(jv), jv = 1, nvinfo ) 
     595 
     596 
     597      !--preparation de "titvar" pour l''ecriture parmi les valeurs numeriques : 
     598      DO  jv = 2 , nvinfo 
     599         titinf     = titvar(jv)(:jpchinf) 
     600         titvar(jv) = '  '//titinf 
     601      END DO 
     602 
     603      !--Initialisation of the arrays for the accumulation 
     604      DO  jv = 1, nvinfo 
     605         vinfom(jv) = 0. 
     606      END DO 
     607      naveg = 0 
     608 
     6091000  FORMAT( 3(A20),4(1x,I6) ) 
     6101111  FORMAT( 3(F7.1,1X,F7.3,1X),I3,A )   
     611 
     612   END SUBROUTINE lim_dia_init 
    613613 
    614614#else 

  • trunk/NEMO/LIM_SRC_3/limdyn.F90

    r913 r921  
    4848CONTAINS 
    4949 
    50    SUBROUTINE lim_dyn 
     50   SUBROUTINE lim_dyn( kt ) 
    5151      !!------------------------------------------------------------------- 
    5252      !!               ***  ROUTINE lim_dyn  *** 
    … …  
    6666      !!                   LIM3, EVP, C-grid 
    6767      !!------------------------------------------------------------------------------------ 
     68      INTEGER, INTENT(in) ::   kt     ! number of iteration 
    6869      !! * Local variables 
    6970      INTEGER  ::   ji, jj, jl, ja    ! dummy loop indices 
    … …  
    7576      !!--------------------------------------------------------------------- 
    7677 
    77       WRITE(numout,*) ' lim_dyn : Ice dynamics ' 
    78       WRITE(numout,*) ' ~~~~~~~ ' 
     78      IF( kt == nit000 .AND. lwp ) THEN 
     79         WRITE(numout,*) ' lim_dyn : Ice dynamics ' 
     80         WRITE(numout,*) ' ~~~~~~~ ' 
     81      ENDIF 
    7982 
    8083      IF( numit == nstart  )   CALL lim_dyn_init   ! Initialization (first time-step only) 
    81        
     84 
    8285      IF ( ln_limdyn ) THEN 
    8386 
    … …  
    219222   END SUBROUTINE lim_dyn 
    220223 
    221     SUBROUTINE lim_dyn_init 
     224   SUBROUTINE lim_dyn_init 
    222225      !!------------------------------------------------------------------- 
    223226      !!                  ***  ROUTINE lim_dyn_init  *** 

  • trunk/NEMO/LIM_SRC_3/limhdf.F90

    r888 r921  
    8484      ! Arrays initialization 
    8585      ptab0 (:, : ) = ptab(:,:) 
    86 !bug  zflu (:,jpj) = 0.e0 
    87 !bug  zflv (:,jpj) = 0.e0 
     86      !bug  zflu (:,jpj) = 0.e0 
     87      !bug  zflv (:,jpj) = 0.e0 
    8888      zdiv0(:, 1 ) = 0.e0 
    8989      zdiv0(:,jpj) = 0.e0 

  • trunk/NEMO/LIM_SRC_3/limistate.F90

    r888 r921  
    8686         zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs  
    8787      !-------------------------------------------------------------------- 
    88   
     88 
    8989      !-------------------------------------------------------------------- 
    9090      ! 1) Preliminary things  
    … …  
    113113      zs0 = 34.e0 
    114114      ztf = ABS ( rt0 - 0.0575       * zs0                               & 
    115                &                    + 1.710523e-03 * zs0 * SQRT( zs0 )   & 
    116                &                    - 2.154996e-04 * zs0 *zs0          ) 
     115         &                    + 1.710523e-03 * zs0 * SQRT( zs0 )   & 
     116         &                    - 2.154996e-04 * zs0 *zs0          ) 
    117117 
    118118      ! constants for heat contents 
    … …  
    179179      ! ------------- 
    180180!!! 
    181 ! retour a LIMA_MEC 
    182 !     ! second ice type 
    183 !     zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
    184 !     hi_max(ice_cat_bounds(2,1)-1) = 0.0 
    185  
    186 !     ! here to change !!!! 
    187 !     jm = 2 
    188 !     DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    189 !        zhin (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    190 !        zhin (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + & 
    191 !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0 
    192 !        zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0) 
    193 !        zhis (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    194 !        zhis (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + & 
    195 !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0 
    196 !        zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0) 
    197 !     END DO ! jl 
    198 !     zgfactorn(2) = aginn_d / zgfactorn(2) 
    199 !     zgfactors(2) = agins_d / zgfactors(2) 
    200 !     hi_max(ice_cat_bounds(2,1)-1) = zdummy 
    201 ! END retour a LIMA_MEC 
     181      ! retour a LIMA_MEC 
     182      !     ! second ice type 
     183      !     zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
     184      !     hi_max(ice_cat_bounds(2,1)-1) = 0.0 
     185 
     186      !     ! here to change !!!! 
     187      !     jm = 2 
     188      !     DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
     189      !        zhin (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
     190      !        zhin (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + & 
     191      !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0 
     192      !        zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0) 
     193      !        zhis (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
     194      !        zhis (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + & 
     195      !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0 
     196      !        zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0) 
     197      !     END DO ! jl 
     198      !     zgfactorn(2) = aginn_d / zgfactorn(2) 
     199      !     zgfactors(2) = agins_d / zgfactors(2) 
     200      !     hi_max(ice_cat_bounds(2,1)-1) = zdummy 
     201      ! END retour a LIMA_MEC 
    202202!!! 
    203203      DO jj = 1, jpj 
    … …  
    228228                     zhin(1)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    229229                     a_i(ji,jj,jl)    = zidto * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* &  
    230                                             (zhin(1)-hginn_u)/2.0) , epsi06) 
     230                        (zhin(1)-hginn_u)/2.0) , epsi06) 
    231231                     ! new line 
    232232                     a_i(ji,jj,jl)    = zidto * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) ) 
    … …  
    239239 
    240240!!! 
    241 ! retour a LIMA_MEC 
    242 !              !ridged ice 
    243 !              zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
    244 !              hi_max(ice_cat_bounds(2,1)-1) = 0.0 
    245 !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories 
    246 !                 zhin(2)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    247 !                 a_i(ji,jj,jl)    = zidto * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* & 
    248 !                                    (zhin(2)-hginn_d)/2.0) , epsi06) 
    249 !                 ht_i(ji,jj,jl)   = zidto * zhin(2)  
    250 !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    251 !              END DO 
    252 !              hi_max(ice_cat_bounds(2,1)-1) = zdummy 
    253  
    254 !              !rafted ice 
    255 !              jl = 6 
    256 !              a_i(ji,jj,jl)       = 0.0 
    257 !              ht_i(ji,jj,jl)      = 0.0 
    258 !              v_i(ji,jj,jl)       = 0.0 
    259 ! END retour a LIMA_MEC 
     241               ! retour a LIMA_MEC 
     242               !              !ridged ice 
     243               !              zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
     244               !              hi_max(ice_cat_bounds(2,1)-1) = 0.0 
     245               !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories 
     246               !                 zhin(2)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
     247               !                 a_i(ji,jj,jl)    = zidto * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* & 
     248               !                                    (zhin(2)-hginn_d)/2.0) , epsi06) 
     249               !                 ht_i(ji,jj,jl)   = zidto * zhin(2)  
     250               !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
     251               !              END DO 
     252               !              hi_max(ice_cat_bounds(2,1)-1) = zdummy 
     253 
     254               !              !rafted ice 
     255               !              jl = 6 
     256               !              a_i(ji,jj,jl)       = 0.0 
     257               !              ht_i(ji,jj,jl)      = 0.0 
     258               !              v_i(ji,jj,jl)       = 0.0 
     259               ! END retour a LIMA_MEC 
    260260!!! 
    261261 
    … …  
    279279                  o_i(ji,jj,jl)    = zidto * 1.0   + ( 1.0 - zidto ) 
    280280                  oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl) 
    281                       
     281 
    282282                  !------------------------------ 
    283283                  ! Sea ice surface temperature 
    … …  
    298298                     ! Multiply by volume, so that heat content in 10^9 Joules 
    299299                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 
    300                                         v_s(ji,jj,jl)  / nlay_s 
     300                        v_s(ji,jj,jl)  / nlay_s 
    301301                  END DO !jk 
    302302 
    … …  
    309309                     s_i(ji,jj,jk,jl) = zidto * sinn + ( 1.0 - zidto ) * 0.1 
    310310                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
    311   
    312                   ! heat content per unit volume 
     311 
     312                     ! heat content per unit volume 
    313313                     e_i(ji,jj,jk,jl) = zidto * rhoic * & 
    314                      (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    315                      +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 
    316                      - rcp      * ( ztmelts - rtt ) & 
    317                      ) 
    318  
    319                   ! Correct dimensions to avoid big values 
     314                        (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
     315                        +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 
     316                        - rcp      * ( ztmelts - rtt ) & 
     317                        ) 
     318 
     319                     ! Correct dimensions to avoid big values 
    320320                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    321321 
    322                   ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
     322                     ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    323323                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &  
    324                                   area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 
    325                                   nlay_i 
     324                        area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 
     325                        nlay_i 
    326326                  END DO ! jk 
    327327 
    … …  
    330330            ELSE ! on fcor  
    331331 
    332             !--- Southern hemisphere 
    333             !---------------------------------------------------------------- 
     332               !--- Southern hemisphere 
     333               !---------------------------------------------------------------- 
    334334 
    335335               !----------------------- 
    … …  
    346346 
    347347               ELSE ! several categories 
    348                 
    349                !level ice 
     348 
     349                  !level ice 
    350350                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories 
    351351 
    352352                     zhis(1)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    353353                     a_i(ji,jj,jl) = zidto * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * &  
    354                                         (zhis(1)-hgins_u)/2.0) , epsi06 ) 
     354                        (zhis(1)-hgins_u)/2.0) , epsi06 ) 
    355355                     ! new line square distribution volume conserving 
    356356                     a_i(ji,jj,jl)    = zidto * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) ) 
    357357                     ht_i(ji,jj,jl)   = zidto * zhis(1)  
    358358                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    359                    
     359 
    360360                  END DO ! jl 
    361361 
    … …  
    363363 
    364364!!! 
    365 ! retour a LIMA_MEC 
    366 !              !ridged ice 
    367 !              zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
    368 !              hi_max(ice_cat_bounds(2,1)-1) = 0.0 
    369 !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories 
    370 !                 zhis(2)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    371 !                 a_i(ji,jj,jl) = zidto*MAX( zgfactors(2) * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 ) 
    372 !                 ht_i(ji,jj,jl)   = zidto * zhis(2)  
    373 !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    374 !              END DO 
    375 !              hi_max(ice_cat_bounds(2,1)-1) = zdummy 
    376  
    377 !              !rafted ice 
    378 !              jl = 6 
    379 !              a_i(ji,jj,jl)       = 0.0 
    380 !              ht_i(ji,jj,jl)      = 0.0 
    381 !              v_i(ji,jj,jl)       = 0.0 
    382 ! END retour a LIMA_MEC 
     365               ! retour a LIMA_MEC 
     366               !              !ridged ice 
     367               !              zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
     368               !              hi_max(ice_cat_bounds(2,1)-1) = 0.0 
     369               !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories 
     370               !                 zhis(2)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
     371               !                 a_i(ji,jj,jl) = zidto*MAX( zgfactors(2) * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 ) 
     372               !                 ht_i(ji,jj,jl)   = zidto * zhis(2)  
     373               !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
     374               !              END DO 
     375               !              hi_max(ice_cat_bounds(2,1)-1) = zdummy 
     376 
     377               !              !rafted ice 
     378               !              jl = 6 
     379               !              a_i(ji,jj,jl)       = 0.0 
     380               !              ht_i(ji,jj,jl)      = 0.0 
     381               !              v_i(ji,jj,jl)       = 0.0 
     382               ! END retour a LIMA_MEC 
    383383!!! 
    384384 
    … …  
    424424                     ! Multiply by volume, so that heat content in 10^9 Joules 
    425425                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 
    426                                         v_s(ji,jj,jl)  / nlay_s 
     426                        v_s(ji,jj,jl)  / nlay_s 
    427427                  END DO 
    428428 
    … …  
    435435                     s_i(ji,jj,jk,jl) = zidto * sins + ( 1.0 - zidto ) * 0.1 
    436436                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
    437   
    438                   ! heat content per unit volume 
     437 
     438                     ! heat content per unit volume 
    439439                     e_i(ji,jj,jk,jl) = zidto * rhoic * & 
    440                      (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    441                      +   lfus  * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 
    442                      - rcp      * ( ztmelts - rtt ) & 
    443                      ) 
    444  
    445                   ! Correct dimensions to avoid big values 
     440                        (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
     441                        +   lfus  * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 
     442                        - rcp      * ( ztmelts - rtt ) & 
     443                        ) 
     444 
     445                     ! Correct dimensions to avoid big values 
    446446                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    447447 
    448                   ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
     448                     ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    449449                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &  
    450                                      area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 
    451                                      nlay_i 
     450                        area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 
     451                        nlay_i 
    452452                  END DO !jk 
    453453 
    … …  
    549549      !!----------------------------------------------------------------------------- 
    550550      NAMELIST/namiceini/ ttest, hninn, hginn_u, aginn_u, hginn_d, aginn_d, hnins, & 
    551                           hgins_u, agins_u, hgins_d, agins_d, sinn, sins 
     551         hgins_u, agins_u, hgins_d, agins_d, sinn, sins 
    552552      !!----------------------------------------------------------------------------- 
    553553 
    … …  
    576576         WRITE(numout,*) '   initial  ice salinity       in the south     sins       = ', sins 
    577577      ENDIF 
    578              
     578 
    579579   END SUBROUTINE lim_istate_init 
    580580 

  • trunk/NEMO/LIM_SRC_3/limitd_me.F90

    r903 r921  
    3232   USE prtctl           ! Print control 
    3333   USE lib_mpp 
    34   
     34 
    3535   IMPLICIT NONE 
    3636   PRIVATE 
    … …  
    5353      zone   = 1.e0 
    5454 
    55 !----------------------------------------------------------------------- 
    56 ! Variables shared among ridging subroutines 
    57 !----------------------------------------------------------------------- 
    58       REAL(wp), DIMENSION (jpi,jpj) ::    & 
    59          asum         , & ! sum of total ice and open water area 
    60          aksum            ! ratio of area removed to area ridged 
    61  
    62       REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: &      
    63          athorn           ! participation function; fraction of ridging/ 
    64                           !  closing associated w/ category n 
    65  
    66       REAL(wp), DIMENSION(jpi,jpj,jpl) ::  & 
    67          hrmin      , &   ! minimum ridge thickness 
    68          hrmax      , &   ! maximum ridge thickness 
    69          hraft      , &   ! thickness of rafted ice 
    70          krdg       , &   ! mean ridge thickness/thickness of ridging ice  
    71          aridge     , &   ! participating ice ridging 
    72          araft            ! participating ice rafting 
    73  
    74       REAL(wp), PARAMETER :: & 
    75          krdgmin = 1.1, &    ! min ridge thickness multiplier 
    76          kraft   = 2.0       ! rafting multipliyer 
    77  
    78       REAL(wp) :: &                                
    79          Cp  
    80 ! 
    81 !----------------------------------------------------------------------- 
    82 ! Ridging diagnostic arrays for history files 
    83 !----------------------------------------------------------------------- 
    84 ! 
    85       REAL (wp), DIMENSION(jpi,jpj) :: & 
    86          dardg1dt     , & ! rate of fractional area loss by ridging ice (1/s) 
    87          dardg2dt     , & ! rate of fractional area gain by new ridges (1/s) 
    88          dvirdgdt     , & ! rate of ice volume ridged (m/s) 
    89          opening          ! rate of opening due to divergence/shear (1/s) 
    90                                         
     55   !----------------------------------------------------------------------- 
     56   ! Variables shared among ridging subroutines 
     57   !----------------------------------------------------------------------- 
     58   REAL(wp), DIMENSION (jpi,jpj) ::    & 
     59      asum         , & ! sum of total ice and open water area 
     60      aksum            ! ratio of area removed to area ridged 
     61 
     62   REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: &      
     63      athorn           ! participation function; fraction of ridging/ 
     64   !  closing associated w/ category n 
     65 
     66   REAL(wp), DIMENSION(jpi,jpj,jpl) ::  & 
     67      hrmin      , &   ! minimum ridge thickness 
     68      hrmax      , &   ! maximum ridge thickness 
     69      hraft      , &   ! thickness of rafted ice 
     70      krdg       , &   ! mean ridge thickness/thickness of ridging ice  
     71      aridge     , &   ! participating ice ridging 
     72      araft            ! participating ice rafting 
     73 
     74   REAL(wp), PARAMETER :: & 
     75      krdgmin = 1.1, &    ! min ridge thickness multiplier 
     76      kraft   = 2.0       ! rafting multipliyer 
     77 
     78   REAL(wp) :: &                                
     79      Cp  
     80   ! 
     81   !----------------------------------------------------------------------- 
     82   ! Ridging diagnostic arrays for history files 
     83   !----------------------------------------------------------------------- 
     84   ! 
     85   REAL (wp), DIMENSION(jpi,jpj) :: & 
     86      dardg1dt     , & ! rate of fractional area loss by ridging ice (1/s) 
     87      dardg2dt     , & ! rate of fractional area gain by new ridges (1/s) 
     88      dvirdgdt     , & ! rate of ice volume ridged (m/s) 
     89      opening          ! rate of opening due to divergence/shear (1/s) 
     90 
    9191 
    9292   !!---------------------------------------------------------------------- 
    … …  
    9797CONTAINS 
    9898 
    99 !!-----------------------------------------------------------------------------! 
    100 !!-----------------------------------------------------------------------------! 
     99   !!-----------------------------------------------------------------------------! 
     100   !!-----------------------------------------------------------------------------! 
    101101 
    102102   SUBROUTINE lim_itd_me ! (subroutine 1/6) 
    103         !!---------------------------------------------------------------------! 
    104         !!                ***  ROUTINE lim_itd_me *** 
    105         !! ** Purpose : 
    106         !!        This routine computes the mechanical redistribution 
    107         !!                      of ice thickness 
    108         !! 
    109         !! ** Method  : a very simple method :-) 
    110         !! 
    111         !! ** Arguments : 
    112         !!           kideb , kiut : Starting and ending points on which the  
    113         !!                         the computation is applied 
    114         !! 
    115         !! ** Inputs / Ouputs : (global commons) 
    116         !! 
    117         !! ** External :  
    118         !! 
    119         !! ** Steps : 
    120         !!  1) Thickness categories boundaries, ice / o.w. concentrations 
    121         !!     Ridge preparation 
    122         !!  2) Dynamical inputs (closing rate, divu_adv, opning) 
    123         !!  3) Ridging iteration 
    124         !!  4) Ridging diagnostics 
    125         !!  5) Heat, salt and freshwater fluxes 
    126         !!  6) Compute increments of tate variables and come back to old values 
    127         !! 
    128         !! ** References : There are a lot of references and can be difficult /  
    129         !!                 boring to read 
    130         !! 
    131         !! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength 
    132         !!  in modeling the thickness distribution of Arctic sea ice, 
    133         !!  J. Geophys. Res., 100, 18,611-18,626. 
    134         !! 
    135         !! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice 
    136         !!  cover, Mon. Wea. Rev., 108, 1943-1973, 1980. 
    137         !! 
    138         !! Rothrock, D. A., 1975: The energetics of the plastic deformation of 
    139         !!  pack ice by ridging, J. Geophys. Res., 80, 4514-4519. 
    140         !! 
    141         !! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony,  
    142         !!  1975: The thickness distribution of sea ice, J. Geophys. Res.,  
    143         !!  80, 4501-4513.  
    144         !! 
    145         !! Bitz et al., JGR 2001 
    146         !! 
    147         !! Amundrud and Melling, JGR 2005 
    148         !! 
    149         !! Babko et al., JGR 2002  
    150         !! 
    151         !! ** History : 
    152         !!           This routine is based on CICE code 
    153         !!           and authors William H. Lipscomb, 
    154         !!           and Elizabeth C. Hunke, LANL 
    155         !!           are gratefully acknowledged 
    156         !! 
    157         !!           (02-2006) Martin Vancoppenolle, UCL-ASTR  
    158         !! 
    159         !!--------------------------------------------------------------------! 
    160         !! * Arguments 
    161  
    162         !! * Local variables 
    163         INTEGER ::   ji,       &   ! spatial dummy loop index 
    164                      jj,       &   ! spatial dummy loop index 
    165                      jk,       &   ! vertical layering dummy loop index 
    166                      jl,       &   ! ice category dummy loop index 
    167                      niter,    &   ! iteration counter 
    168                      nitermax = 20 ! max number of ridging iterations  
    169  
    170         REAL(wp)  ::             &  ! constant values 
    171            zeps      =  1.0e-10, & 
    172            epsi10    =  1.0e-10, & 
    173            epsi06    =  1.0e-6 
    174  
    175         REAL(wp), DIMENSION(jpi,jpj) :: & 
    176            closing_net,        &  ! net rate at which area is removed    (1/s) 
    177                                   ! (ridging ice area - area of new ridges) / dt 
    178            divu_adv   ,        &  ! divu as implied by transport scheme  (1/s) 
    179            opning     ,        &  ! rate of opening due to divergence/shear 
    180            closing_gross,      &  ! rate at which area removed, not counting 
    181                                   ! area of new ridges 
    182            msnow_mlt  ,        &  ! mass of snow added to ocean (kg m-2) 
    183            esnow_mlt              ! energy needed to melt snow in ocean (J m-2) 
    184  
    185         REAL(wp) ::            & 
    186            w1,                 &  ! temporary variable 
    187            tmpfac,             &  ! factor by which opening/closing rates are cut 
    188            dti                    ! 1 / dt 
    189  
    190         LOGICAL   ::           & 
    191            asum_error              ! flag for asum .ne. 1 
    192  
    193         INTEGER :: iterate_ridging ! if true, repeat the ridging 
    194  
    195         REAL(wp) ::  &          
    196            big = 1.0e8 
    197  
    198         REAL (wp), DIMENSION(jpi,jpj) :: &  !  
    199            vt_i_init, vt_i_final       !  ice volume summed over categories 
    200  
    201         CHARACTER (len = 15) :: fieldid 
    202  
    203 !!-- End of declarations 
    204 !-----------------------------------------------------------------------------! 
     103      !!---------------------------------------------------------------------! 
     104      !!                ***  ROUTINE lim_itd_me *** 
     105      !! ** Purpose : 
     106      !!        This routine computes the mechanical redistribution 
     107      !!                      of ice thickness 
     108      !! 
     109      !! ** Method  : a very simple method :-) 
     110      !! 
     111      !! ** Arguments : 
     112      !!           kideb , kiut : Starting and ending points on which the  
     113      !!                         the computation is applied 
     114      !! 
     115      !! ** Inputs / Ouputs : (global commons) 
     116      !! 
     117      !! ** External :  
     118      !! 
     119      !! ** Steps : 
     120      !!  1) Thickness categories boundaries, ice / o.w. concentrations 
     121      !!     Ridge preparation 
     122      !!  2) Dynamical inputs (closing rate, divu_adv, opning) 
     123      !!  3) Ridging iteration 
     124      !!  4) Ridging diagnostics 
     125      !!  5) Heat, salt and freshwater fluxes 
     126      !!  6) Compute increments of tate variables and come back to old values 
     127      !! 
     128      !! ** References : There are a lot of references and can be difficult /  
     129      !!                 boring to read 
     130      !! 
     131      !! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength 
     132      !!  in modeling the thickness distribution of Arctic sea ice, 
     133      !!  J. Geophys. Res., 100, 18,611-18,626. 
     134      !! 
     135      !! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice 
     136      !!  cover, Mon. Wea. Rev., 108, 1943-1973, 1980. 
     137      !! 
     138      !! Rothrock, D. A., 1975: The energetics of the plastic deformation of 
     139      !!  pack ice by ridging, J. Geophys. Res., 80, 4514-4519. 
     140      !! 
     141      !! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony,  
     142      !!  1975: The thickness distribution of sea ice, J. Geophys. Res.,  
     143      !!  80, 4501-4513.  
     144      !! 
     145      !! Bitz et al., JGR 2001 
     146      !! 
     147      !! Amundrud and Melling, JGR 2005 
     148      !! 
     149      !! Babko et al., JGR 2002  
     150      !! 
     151      !! ** History : 
     152      !!           This routine is based on CICE code 
     153      !!           and authors William H. Lipscomb, 
     154      !!           and Elizabeth C. Hunke, LANL 
     155      !!           are gratefully acknowledged 
     156      !! 
     157      !!           (02-2006) Martin Vancoppenolle, UCL-ASTR  
     158      !! 
     159      !!--------------------------------------------------------------------! 
     160      !! * Arguments 
     161 
     162      !! * Local variables 
     163      INTEGER ::   ji,       &   ! spatial dummy loop index 
     164         jj,       &   ! spatial dummy loop index 
     165         jk,       &   ! vertical layering dummy loop index 
     166         jl,       &   ! ice category dummy loop index 
     167         niter,    &   ! iteration counter 
     168         nitermax = 20 ! max number of ridging iterations  
     169 
     170      REAL(wp)  ::             &  ! constant values 
     171         zeps      =  1.0e-10, & 
     172         epsi10    =  1.0e-10, & 
     173         epsi06    =  1.0e-6 
     174 
     175      REAL(wp), DIMENSION(jpi,jpj) :: & 
     176         closing_net,        &  ! net rate at which area is removed    (1/s) 
     177                                ! (ridging ice area - area of new ridges) / dt 
     178         divu_adv   ,        &  ! divu as implied by transport scheme  (1/s) 
     179         opning     ,        &  ! rate of opening due to divergence/shear 
     180         closing_gross,      &  ! rate at which area removed, not counting 
     181                                ! area of new ridges 
     182         msnow_mlt  ,        &  ! mass of snow added to ocean (kg m-2) 
     183         esnow_mlt              ! energy needed to melt snow in ocean (J m-2) 
     184 
     185      REAL(wp) ::            & 
     186         w1,                 &  ! temporary variable 
     187         tmpfac,             &  ! factor by which opening/closing rates are cut 
     188         dti                    ! 1 / dt 
     189 
     190      LOGICAL   ::           & 
     191         asum_error              ! flag for asum .ne. 1 
     192 
     193      INTEGER :: iterate_ridging ! if true, repeat the ridging 
     194 
     195      REAL(wp) ::  &          
     196         big = 1.0e8 
     197 
     198      REAL (wp), DIMENSION(jpi,jpj) :: &  !  
     199         vt_i_init, vt_i_final       !  ice volume summed over categories 
     200 
     201      CHARACTER (len = 15) :: fieldid 
     202 
     203      !!-- End of declarations 
     204      !-----------------------------------------------------------------------------! 
    205205 
    206206      IF( numit == nstart  ) CALL lim_itd_me_init ! Initialization (first time-step only) 
    … …  
    211211      ENDIF 
    212212 
    213 !-----------------------------------------------------------------------------! 
    214 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
    215 !-----------------------------------------------------------------------------! 
    216 ! Set hi_max(ncat) to a big value to ensure that all ridged ice  
    217 ! is thinner than hi_max(ncat). 
     213      !-----------------------------------------------------------------------------! 
     214      ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 
     215      !-----------------------------------------------------------------------------! 
     216      ! Set hi_max(ncat) to a big value to ensure that all ridged ice  
     217      ! is thinner than hi_max(ncat). 
    218218 
    219219      hi_max(jpl) = 999.99 
    … …  
    225225      IF ( con_i) CALL lim_column_sum (jpl,   v_i, vt_i_init) 
    226226 
    227 ! Initialize arrays. 
     227      ! Initialize arrays. 
    228228      DO jj = 1, jpj 
    229229         DO ji = 1, jpi 
    230230 
    231          msnow_mlt(ji,jj) = 0.0 
    232          esnow_mlt(ji,jj) = 0.0 
    233          dardg1dt(ji,jj)  = 0.0 
    234          dardg2dt(ji,jj)  = 0.0 
    235          dvirdgdt(ji,jj)  = 0.0 
    236          opening (ji,jj)  = 0.0 
    237  
    238 !-----------------------------------------------------------------------------! 
    239 ! 2) Dynamical inputs (closing rate, divu_adv, opning) 
    240 !-----------------------------------------------------------------------------! 
    241 ! 
    242 ! 2.1 closing_net 
    243 !----------------- 
    244 ! Compute the net rate of closing due to convergence  
    245 ! and shear, based on Flato and Hibler (1995). 
    246 !  
    247 ! The energy dissipation rate is equal to the net closing rate 
    248 ! times the ice strength. 
    249 ! 
    250 ! NOTE: The NET closing rate is equal to the rate that open water  
    251 !  area is removed, plus the rate at which ice area is removed by  
    252 !  ridging, minus the rate at which area is added in new ridges. 
    253 !  The GROSS closing rate is equal to the first two terms (open 
    254 !  water closing and thin ice ridging) without the third term 
    255 !  (thick, newly ridged ice). 
    256  
    257          closing_net(ji,jj) = & 
    258               Cs*0.5*(Delta_i(ji,jj)-ABS(divu_i(ji,jj))) - MIN(divu_i(ji,jj),0.0) 
    259  
    260 ! 2.2 divu_adv 
    261 !-------------- 
    262 ! Compute divu_adv, the divergence rate given by the transport/ 
    263 ! advection scheme, which may not be equal to divu as computed  
    264 ! from the velocity field. 
    265 ! 
    266 ! If divu_adv < 0, make sure the closing rate is large enough 
    267 ! to give asum = 1.0 after ridging. 
    268  
    269          divu_adv(ji,jj) = (1.0-asum(ji,jj)) / rdt_ice  ! asum found in ridgeprep 
    270  
    271          IF (divu_adv(ji,jj) .LT. 0.0) & 
    272               closing_net(ji,jj) = max(closing_net(ji,jj), -divu_adv(ji,jj)) 
    273  
    274 ! 2.3 opning 
    275 !------------ 
    276 ! Compute the (non-negative) opening rate that will give  
    277 ! asum = 1.0 after ridging. 
    278          opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj) 
     231            msnow_mlt(ji,jj) = 0.0 
     232            esnow_mlt(ji,jj) = 0.0 
     233            dardg1dt(ji,jj)  = 0.0 
     234            dardg2dt(ji,jj)  = 0.0 
     235            dvirdgdt(ji,jj)  = 0.0 
     236            opening (ji,jj)  = 0.0 
     237 
     238            !-----------------------------------------------------------------------------! 
     239            ! 2) Dynamical inputs (closing rate, divu_adv, opning) 
     240            !-----------------------------------------------------------------------------! 
     241            ! 
     242            ! 2.1 closing_net 
     243            !----------------- 
     244            ! Compute the net rate of closing due to convergence  
     245            ! and shear, based on Flato and Hibler (1995). 
     246            !  
     247            ! The energy dissipation rate is equal to the net closing rate 
     248            ! times the ice strength. 
     249            ! 
     250            ! NOTE: The NET closing rate is equal to the rate that open water  
     251            !  area is removed, plus the rate at which ice area is removed by  
     252            !  ridging, minus the rate at which area is added in new ridges. 
     253            !  The GROSS closing rate is equal to the first two terms (open 
     254            !  water closing and thin ice ridging) without the third term 
     255            !  (thick, newly ridged ice). 
     256 
     257            closing_net(ji,jj) = & 
     258               Cs*0.5*(Delta_i(ji,jj)-ABS(divu_i(ji,jj))) - MIN(divu_i(ji,jj),0.0) 
     259 
     260            ! 2.2 divu_adv 
     261            !-------------- 
     262            ! Compute divu_adv, the divergence rate given by the transport/ 
     263            ! advection scheme, which may not be equal to divu as computed  
     264            ! from the velocity field. 
     265            ! 
     266            ! If divu_adv < 0, make sure the closing rate is large enough 
     267            ! to give asum = 1.0 after ridging. 
     268 
     269            divu_adv(ji,jj) = (1.0-asum(ji,jj)) / rdt_ice  ! asum found in ridgeprep 
     270 
     271            IF (divu_adv(ji,jj) .LT. 0.0) & 
     272               closing_net(ji,jj) = max(closing_net(ji,jj), -divu_adv(ji,jj)) 
     273 
     274            ! 2.3 opning 
     275            !------------ 
     276            ! Compute the (non-negative) opening rate that will give  
     277            ! asum = 1.0 after ridging. 
     278            opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj) 
    279279 
    280280         END DO 
    281281      END DO 
    282282 
    283 !-----------------------------------------------------------------------------! 
    284 ! 3) Ridging iteration 
    285 !-----------------------------------------------------------------------------! 
     283      !-----------------------------------------------------------------------------! 
     284      ! 3) Ridging iteration 
     285      !-----------------------------------------------------------------------------! 
    286286      niter = 1                 ! iteration counter 
    287287      iterate_ridging = 1 
    … …  
    290290      DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax ) 
    291291 
     292         DO jj = 1, jpj 
     293            DO ji = 1, jpi 
     294 
     295               ! 3.2 closing_gross 
     296               !-----------------------------------------------------------------------------! 
     297               ! Based on the ITD of ridging and ridged ice, convert the net 
     298               !  closing rate to a gross closing rate.   
     299               ! NOTE: 0 < aksum <= 1 
     300               closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj) 
     301 
     302               ! correction to closing rate and opening if closing rate is excessive 
     303               !--------------------------------------------------------------------- 
     304               ! Reduce the closing rate if more than 100% of the open water  
     305               ! would be removed.  Reduce the opening rate proportionately. 
     306               IF ( ato_i(ji,jj) .GT. epsi11 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 
     307                  w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 
     308                  IF ( w1 .GT. ato_i(ji,jj)) THEN 
     309                     tmpfac = ato_i(ji,jj) / w1 
     310                     closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 
     311                     opning(ji,jj) = opning(ji,jj) * tmpfac 
     312                  ENDIF !w1 
     313               ENDIF !at0i and athorn 
     314 
     315            END DO ! ji 
     316         END DO ! jj 
     317 
     318         ! correction to closing rate / opening if excessive ice removal 
     319         !--------------------------------------------------------------- 
     320         ! Reduce the closing rate if more than 100% of any ice category  
     321         ! would be removed.  Reduce the opening rate proportionately. 
     322 
     323         DO jl = 1, jpl 
     324            DO jj = 1, jpj 
     325               DO ji = 1, jpi 
     326                  IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 
     327                     w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
     328                     IF ( w1 .GT. a_i(ji,jj,jl) ) THEN 
     329                        tmpfac = a_i(ji,jj,jl) / w1 
     330                        closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 
     331                        opning(ji,jj) = opning(ji,jj) * tmpfac 
     332                     ENDIF 
     333                  ENDIF 
     334               END DO !ji 
     335            END DO ! jj 
     336         END DO !jl 
     337 
     338         ! 3.3 Redistribute area, volume, and energy. 
     339         !-----------------------------------------------------------------------------! 
     340 
     341         CALL lim_itd_me_ridgeshift (opning,    closing_gross, & 
     342            msnow_mlt, esnow_mlt) 
     343 
     344         ! 3.4 Compute total area of ice plus open water after ridging. 
     345         !-----------------------------------------------------------------------------! 
     346 
     347         CALL lim_itd_me_asumr 
     348 
     349         ! 3.5 Do we keep on iterating ??? 
     350         !-----------------------------------------------------------------------------! 
     351         ! Check whether asum = 1.  If not (because the closing and opening 
     352         ! rates were reduced above), ridge again with new rates. 
     353 
     354         iterate_ridging = 0 
     355 
     356         DO jj = 1, jpj 
     357            DO ji = 1, jpi 
     358               IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN 
     359                  closing_net(ji,jj) = 0.0  
     360                  opning(ji,jj)      = 0.0 
     361               ELSE 
     362                  iterate_ridging    = 1 
     363                  divu_adv(ji,jj)    = (1.0 - asum(ji,jj)) / rdt_ice 
     364                  closing_net(ji,jj) = MAX(0.0, -divu_adv(ji,jj)) 
     365                  opning(ji,jj)      = MAX(0.0, divu_adv(ji,jj)) 
     366               ENDIF 
     367            END DO 
     368         END DO 
     369 
     370         IF( lk_mpp ) CALL mpp_max(iterate_ridging) 
     371 
     372         ! Repeat if necessary. 
     373         ! NOTE: If strength smoothing is turned on, the ridging must be 
     374         ! iterated globally because of the boundary update in the  
     375         ! smoothing. 
     376 
     377         niter = niter + 1 
     378 
     379         IF (iterate_ridging == 1) THEN 
     380            IF (niter .GT. nitermax) THEN 
     381               WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 
     382               WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging 
     383            ENDIF 
     384            CALL lim_itd_me_ridgeprep 
     385         ENDIF 
     386 
     387      END DO !! on the do while over iter 
     388 
     389      !-----------------------------------------------------------------------------! 
     390      ! 4) Ridging diagnostics 
     391      !-----------------------------------------------------------------------------! 
     392      ! Convert ridging rate diagnostics to correct units. 
     393      ! Update fresh water and heat fluxes due to snow melt. 
     394 
     395      dti = 1.0/rdt_ice 
     396 
     397      asum_error = .false.  
     398 
    292399      DO jj = 1, jpj 
    293400         DO ji = 1, jpi 
    294401 
    295 ! 3.2 closing_gross 
    296 !-----------------------------------------------------------------------------! 
    297 ! Based on the ITD of ridging and ridged ice, convert the net 
    298 !  closing rate to a gross closing rate.   
    299 ! NOTE: 0 < aksum <= 1 
    300             closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj) 
    301  
    302 ! correction to closing rate and opening if closing rate is excessive 
    303 !--------------------------------------------------------------------- 
    304 ! Reduce the closing rate if more than 100% of the open water  
    305 ! would be removed.  Reduce the opening rate proportionately. 
    306             IF ( ato_i(ji,jj) .GT. epsi11 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 
    307                w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 
    308                IF ( w1 .GT. ato_i(ji,jj)) THEN 
    309                   tmpfac = ato_i(ji,jj) / w1 
    310                   closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 
    311                   opning(ji,jj) = opning(ji,jj) * tmpfac 
    312                ENDIF !w1 
    313             ENDIF !at0i and athorn 
    314  
    315          END DO ! ji 
    316       END DO ! jj 
    317  
    318 ! correction to closing rate / opening if excessive ice removal 
    319 !--------------------------------------------------------------- 
    320 ! Reduce the closing rate if more than 100% of any ice category  
    321 ! would be removed.  Reduce the opening rate proportionately. 
    322  
    323       DO jl = 1, jpl 
    324          DO jj = 1, jpj 
    325             DO ji = 1, jpi 
    326                IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 
    327                   w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 
    328                   IF ( w1 .GT. a_i(ji,jj,jl) ) THEN 
    329                      tmpfac = a_i(ji,jj,jl) / w1 
    330                      closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 
    331                      opning(ji,jj) = opning(ji,jj) * tmpfac 
    332                   ENDIF 
    333                ENDIF 
    334             END DO !ji 
    335          END DO ! jj 
    336       END DO !jl 
    337  
    338 ! 3.3 Redistribute area, volume, and energy. 
    339 !-----------------------------------------------------------------------------! 
    340  
    341       CALL lim_itd_me_ridgeshift (opning,    closing_gross, & 
    342                         msnow_mlt, esnow_mlt) 
    343  
    344 ! 3.4 Compute total area of ice plus open water after ridging. 
    345 !-----------------------------------------------------------------------------! 
    346  
    347       CALL lim_itd_me_asumr 
    348  
    349 ! 3.5 Do we keep on iterating ??? 
    350 !-----------------------------------------------------------------------------! 
    351 ! Check whether asum = 1.  If not (because the closing and opening 
    352 ! rates were reduced above), ridge again with new rates. 
    353  
    354       iterate_ridging = 0 
    355  
    356       DO jj = 1, jpj 
    357          DO ji = 1, jpi 
    358             IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN 
    359                closing_net(ji,jj) = 0.0  
    360                opning(ji,jj)      = 0.0 
    361             ELSE 
    362                iterate_ridging    = 1 
    363                divu_adv(ji,jj)    = (1.0 - asum(ji,jj)) / rdt_ice 
    364                closing_net(ji,jj) = MAX(0.0, -divu_adv(ji,jj)) 
    365                opning(ji,jj)      = MAX(0.0, divu_adv(ji,jj)) 
    366             ENDIF 
    367          END DO 
    368       END DO 
    369  
    370       IF( lk_mpp ) CALL mpp_max(iterate_ridging) 
    371  
    372 ! Repeat if necessary. 
    373 ! NOTE: If strength smoothing is turned on, the ridging must be 
    374 ! iterated globally because of the boundary update in the  
    375 ! smoothing. 
    376  
    377       niter = niter + 1 
    378  
    379       IF (iterate_ridging == 1) THEN 
    380          IF (niter .GT. nitermax) THEN 
    381             WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 
    382             WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging 
    383          ENDIF 
    384          CALL lim_itd_me_ridgeprep 
    385       ENDIF 
    386  
    387       END DO !! on the do while over iter 
    388  
    389 !-----------------------------------------------------------------------------! 
    390 ! 4) Ridging diagnostics 
    391 !-----------------------------------------------------------------------------! 
    392 ! Convert ridging rate diagnostics to correct units. 
    393 ! Update fresh water and heat fluxes due to snow melt. 
    394  
    395       dti = 1.0/rdt_ice 
    396  
    397       asum_error = .false.  
    398  
    399       DO jj = 1, jpj 
    400          DO ji = 1, jpi 
    401  
    402          IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true. 
    403  
    404          dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti 
    405          dardg2dt(ji,jj) = dardg2dt(ji,jj) * dti 
    406          dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * dti 
    407          opening (ji,jj) = opening (ji,jj) * dti 
    408  
    409 !-----------------------------------------------------------------------------! 
    410 ! 5) Heat, salt and freshwater fluxes 
    411 !-----------------------------------------------------------------------------! 
    412          ! fresh water source for ocean 
    413          fmmec(ji,jj)      = fmmec(ji,jj)      + msnow_mlt(ji,jj)*dti   
    414        
    415          ! heat sink for ocean 
    416          fhmec(ji,jj)      = fhmec(ji,jj)      + esnow_mlt(ji,jj)*dti 
     402            IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true. 
     403 
     404            dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti 
     405            dardg2dt(ji,jj) = dardg2dt(ji,jj) * dti 
     406            dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * dti 
     407            opening (ji,jj) = opening (ji,jj) * dti 
     408 
     409            !-----------------------------------------------------------------------------! 
     410            ! 5) Heat, salt and freshwater fluxes 
     411            !-----------------------------------------------------------------------------! 
     412            ! fresh water source for ocean 
     413            fmmec(ji,jj)      = fmmec(ji,jj)      + msnow_mlt(ji,jj)*dti   
     414 
     415            ! heat sink for ocean 
     416            fhmec(ji,jj)      = fhmec(ji,jj)      + esnow_mlt(ji,jj)*dti 
    417417 
    418418         END DO 
    … …  
    444444      ENDIF 
    445445 
    446 !-----------------------------------------------------------------------------! 
    447 ! 6) Updating state variables and trend terms 
    448 !-----------------------------------------------------------------------------! 
     446      !-----------------------------------------------------------------------------! 
     447      ! 6) Updating state variables and trend terms 
     448      !-----------------------------------------------------------------------------! 
    449449 
    450450      CALL lim_var_glo2eqv 
    … …  
    465465      d_smv_i_trp(:,:,:)   = 0.0 
    466466      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 
    467       d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
     467         d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
    468468 
    469469      IF(ln_ctl) THEN     ! Control print 
    … …  
    513513      oa_i(:,:,:)   = old_oa_i(:,:,:) 
    514514      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &  
    515       smv_i(:,:,:)  = old_smv_i(:,:,:) 
     515         smv_i(:,:,:)  = old_smv_i(:,:,:) 
    516516 
    517517      !----------------------------------------------------! 
    … …  
    528528               DO ji = 1, jpi 
    529529                  IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. & 
    530                        ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN 
    531                       old_e_i(ji,jj,jk,jl)   = d_e_i_trp(ji,jj,jk,jl) 
    532                       d_e_i_trp(ji,jj,jk,jl) = 0.0 
     530                     ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN 
     531                     old_e_i(ji,jj,jk,jl)   = d_e_i_trp(ji,jj,jk,jl) 
     532                     d_e_i_trp(ji,jj,jk,jl) = 0.0 
    533533                  ENDIF 
    534534               END DO 
    … …  
    541541            DO ji = 1, jpi 
    542542               IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. & 
    543                     ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN 
    544                    old_v_i(ji,jj,jl)     = d_v_i_trp(ji,jj,jl) 
    545                    d_v_i_trp(ji,jj,jl)   = 0.0 
    546                    old_a_i(ji,jj,jl)     = d_a_i_trp(ji,jj,jl) 
    547                    d_a_i_trp(ji,jj,jl)   = 0.0 
    548                    old_v_s(ji,jj,jl)     = d_v_s_trp(ji,jj,jl) 
    549                    d_v_s_trp(ji,jj,jl)   = 0.0 
    550                    old_e_s(ji,jj,1,jl)   = d_e_s_trp(ji,jj,1,jl) 
    551                    d_e_s_trp(ji,jj,1,jl) = 0.0 
    552                    old_oa_i(ji,jj,jl)    = d_oa_i_trp(ji,jj,jl) 
    553                    d_oa_i_trp(ji,jj,jl)  = 0.0 
    554                    IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &  
    555                    old_smv_i(ji,jj,jl)   = d_smv_i_trp(ji,jj,jl) 
    556                    d_smv_i_trp(ji,jj,jl) = 0.0 
     543                  ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN 
     544                  old_v_i(ji,jj,jl)     = d_v_i_trp(ji,jj,jl) 
     545                  d_v_i_trp(ji,jj,jl)   = 0.0 
     546                  old_a_i(ji,jj,jl)     = d_a_i_trp(ji,jj,jl) 
     547                  d_a_i_trp(ji,jj,jl)   = 0.0 
     548                  old_v_s(ji,jj,jl)     = d_v_s_trp(ji,jj,jl) 
     549                  d_v_s_trp(ji,jj,jl)   = 0.0 
     550                  old_e_s(ji,jj,1,jl)   = d_e_s_trp(ji,jj,1,jl) 
     551                  d_e_s_trp(ji,jj,1,jl) = 0.0 
     552                  old_oa_i(ji,jj,jl)    = d_oa_i_trp(ji,jj,jl) 
     553                  d_oa_i_trp(ji,jj,jl)  = 0.0 
     554                  IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &  
     555                     old_smv_i(ji,jj,jl)   = d_smv_i_trp(ji,jj,jl) 
     556                  d_smv_i_trp(ji,jj,jl) = 0.0 
    557557               ENDIF 
    558558            END DO 
    559559         END DO 
    560560      END DO 
    561           
     561 
    562562   END SUBROUTINE lim_itd_me 
    563563 
    564 !=============================================================================== 
     564   !=============================================================================== 
    565565 
    566566   SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6) 
    567567 
    568         !!---------------------------------------------------------------------- 
    569         !!                ***  ROUTINE lim_itd_me_icestrength *** 
    570         !! ** Purpose : 
    571         !!        This routine computes ice strength used in dynamics routines 
    572         !!                      of ice thickness 
    573         !! 
    574         !! ** Method  : 
    575         !!       Compute the strength of the ice pack, defined as the energy (J m-2)  
    576         !! dissipated per unit area removed from the ice pack under compression, 
    577         !! and assumed proportional to the change in potential energy caused 
    578         !! by ridging. Note that only Hibler's formulation is stable and that 
    579         !! ice strength has to be smoothed 
    580         !! 
    581         !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using) 
    582         !! 
    583         !! ** External :  
    584         !! 
    585         !! ** References : 
    586         !!                 
    587         !!---------------------------------------------------------------------- 
    588         !! * Arguments 
    589   
     568      !!---------------------------------------------------------------------- 
     569      !!                ***  ROUTINE lim_itd_me_icestrength *** 
     570      !! ** Purpose : 
     571      !!        This routine computes ice strength used in dynamics routines 
     572      !!                      of ice thickness 
     573      !! 
     574      !! ** Method  : 
     575      !!       Compute the strength of the ice pack, defined as the energy (J m-2)  
     576      !! dissipated per unit area removed from the ice pack under compression, 
     577      !! and assumed proportional to the change in potential energy caused 
     578      !! by ridging. Note that only Hibler's formulation is stable and that 
     579      !! ice strength has to be smoothed 
     580      !! 
     581      !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using) 
     582      !! 
     583      !! ** External :  
     584      !! 
     585      !! ** References : 
     586      !!                 
     587      !!---------------------------------------------------------------------- 
     588      !! * Arguments 
     589 
    590590      INTEGER, INTENT(in) :: & 
    591591         kstrngth    ! = 1 for Rothrock formulation, 0 for Hibler (1979) 
    … …  
    606606         zworka              !: temporary array used here 
    607607 
    608 !------------------------------------------------------------------------------! 
    609 ! 1) Initialize 
    610 !------------------------------------------------------------------------------! 
     608      !------------------------------------------------------------------------------! 
     609      ! 1) Initialize 
     610      !------------------------------------------------------------------------------! 
    611611      strength(:,:) = 0.0 
    612612 
    613 !------------------------------------------------------------------------------! 
    614 ! 2) Compute thickness distribution of ridging and ridged ice 
    615 !------------------------------------------------------------------------------! 
     613      !------------------------------------------------------------------------------! 
     614      ! 2) Compute thickness distribution of ridging and ridged ice 
     615      !------------------------------------------------------------------------------! 
    616616      CALL lim_itd_me_ridgeprep 
    617617 
    618 !------------------------------------------------------------------------------! 
    619 ! 3) Rothrock(1975)'s method 
    620 !------------------------------------------------------------------------------! 
     618      !------------------------------------------------------------------------------! 
     619      ! 3) Rothrock(1975)'s method 
     620      !------------------------------------------------------------------------------! 
    621621      IF (kstrngth == 1) then 
    622622 
    … …  
    626626 
    627627                  IF(     ( a_i(ji,jj,jl)    .GT. epsi11 )                     & 
    628                     .AND. ( athorn(ji,jj,jl) .GT. 0.0    ) ) THEN 
     628                     .AND. ( athorn(ji,jj,jl) .GT. 0.0    ) ) THEN 
    629629                     hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
    630630                     !---------------------------- 
    … …  
    632632                     !---------------------------- 
    633633                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) *    & 
    634                      hi * hi 
     634                        hi * hi 
    635635 
    636636                     !-------------------------- 
    … …  
    638638                     !-------------------------- 
    639639                     strength(ji,jj) = strength(ji,jj) + 2.0 * araft(ji,jj,jl) & 
    640                      * hi * hi 
     640                        * hi * hi 
    641641 
    642642                     !---------------------------- 
    … …  
    644644                     !---------------------------- 
    645645                     strength(ji,jj) = strength(ji,jj)                         & 
    646                      + aridge(ji,jj,jl)/krdg(ji,jj,jl)                         & 
    647                      * 1.0/3.0 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3)     & 
    648                      / (hrmax(ji,jj,jl)-hrmin(ji,jj,jl))                       
     646                        + aridge(ji,jj,jl)/krdg(ji,jj,jl)                         & 
     647                        * 1.0/3.0 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3)     & 
     648                        / (hrmax(ji,jj,jl)-hrmin(ji,jj,jl))                       
    649649                  ENDIF            ! aicen > epsi11 
    650650 
    … …  
    656656            DO ji = 1, jpi 
    657657               strength(ji,jj) = Cf * Cp * strength(ji,jj) / aksum(ji,jj)  
    658                           ! Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) 
    659                           ! Cf accounts for frictional dissipation 
    660                 
     658               ! Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) 
     659               ! Cf accounts for frictional dissipation 
     660 
    661661            END DO              ! j 
    662662         END DO                 ! i 
    … …  
    664664         ksmooth = 1 
    665665 
    666 !------------------------------------------------------------------------------! 
    667 ! 4) Hibler (1979)' method 
    668 !------------------------------------------------------------------------------! 
     666         !------------------------------------------------------------------------------! 
     667         ! 4) Hibler (1979)' method 
     668         !------------------------------------------------------------------------------! 
    669669      ELSE                      ! kstrngth ne 1:  Hibler (1979) form 
    670670 
    … …  
    679679      ENDIF                     ! kstrngth 
    680680 
    681 ! 
    682 !------------------------------------------------------------------------------! 
    683 ! 5) Impact of brine volume 
    684 !------------------------------------------------------------------------------! 
    685 ! CAN BE REMOVED 
    686 ! 
     681      ! 
     682      !------------------------------------------------------------------------------! 
     683      ! 5) Impact of brine volume 
     684      !------------------------------------------------------------------------------! 
     685      ! CAN BE REMOVED 
     686      ! 
    687687      IF ( brinstren_swi .EQ. 1 ) THEN 
    688688 
    … …  
    700700      ENDIF 
    701701 
    702 ! 
    703 !------------------------------------------------------------------------------! 
    704 ! 6) Smoothing ice strength 
    705 !------------------------------------------------------------------------------! 
    706 ! 
     702      ! 
     703      !------------------------------------------------------------------------------! 
     704      ! 6) Smoothing ice strength 
     705      !------------------------------------------------------------------------------! 
     706      ! 
    707707      !------------------- 
    708708      ! Spatial smoothing 
    … …  
    715715            DO ji = 2, jpi - 1 
    716716               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is 
    717                                                                      ! present 
     717                  ! present 
    718718                  zworka(ji,jj) = 4.0 * strength(ji,jj)              & 
    719                                   + strength(ji-1,jj) * tms(ji-1,jj) &   
    720                                   + strength(ji+1,jj) * tms(ji+1,jj) &   
    721                                   + strength(ji,jj-1) * tms(ji,jj-1) &   
    722                                   + strength(ji,jj+1) * tms(ji,jj+1)     
     719                     + strength(ji-1,jj) * tms(ji-1,jj) &   
     720                     + strength(ji+1,jj) * tms(ji+1,jj) &   
     721                     + strength(ji,jj-1) * tms(ji,jj-1) &   
     722                     + strength(ji,jj+1) * tms(ji,jj+1)     
    723723 
    724724                  zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj)            & 
    725                             + tms(ji,jj-1) + tms(ji,jj+1) 
     725                     + tms(ji,jj-1) + tms(ji,jj+1) 
    726726                  zworka(ji,jj) = zworka(ji,jj) / zw1 
    727727               ELSE 
    … …  
    749749 
    750750      IF ( ksmooth .EQ. 2 ) THEN 
    751                   
    752           
     751 
     752 
    753753         CALL lbc_lnk( strength, 'T', 1. ) 
    754              
     754 
    755755         DO jj = 1, jpj - 1 
    756756            DO ji = 1, jpi - 1 
    757757               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is 
    758                                                                      ! present 
     758                  ! present 
    759759                  numts_rm = 1 ! number of time steps for the running mean 
    760760                  IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 
    761761                  IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 
    762762                  zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) /   & 
    763                        numts_rm 
     763                     numts_rm 
    764764                  strp2(ji,jj) = strp1(ji,jj) 
    765765                  strp1(ji,jj) = strength(ji,jj) 
    … …  
    771771 
    772772      ENDIF ! ksmooth 
    773        
     773 
    774774      ! Boundary conditions 
    775775      CALL lbc_lnk( strength, 'T', 1. ) 
    776776 
    777       END SUBROUTINE lim_itd_me_icestrength 
    778  
    779 !=============================================================================== 
    780  
    781       SUBROUTINE lim_itd_me_ridgeprep !(subroutine 3/6) 
    782  
    783         !!---------------------------------------------------------------------! 
    784         !!                ***  ROUTINE lim_itd_me_ridgeprep *** 
    785         !! ** Purpose : 
    786         !!         preparation for ridging and strength calculations 
    787         !! 
    788         !! ** Method  : 
    789         !! Compute the thickness distribution of the ice and open water  
    790         !! participating in ridging and of the resulting ridges. 
    791         !! 
    792         !! ** Arguments : 
    793         !! 
    794         !! ** External :  
    795         !! 
    796         !!---------------------------------------------------------------------! 
    797         !! * Arguments 
    798   
     777   END SUBROUTINE lim_itd_me_icestrength 
     778 
     779   !=============================================================================== 
     780 
     781   SUBROUTINE lim_itd_me_ridgeprep !(subroutine 3/6) 
     782 
     783      !!---------------------------------------------------------------------! 
     784      !!                ***  ROUTINE lim_itd_me_ridgeprep *** 
     785      !! ** Purpose : 
     786      !!         preparation for ridging and strength calculations 
     787      !! 
     788      !! ** Method  : 
     789      !! Compute the thickness distribution of the ice and open water  
     790      !! participating in ridging and of the resulting ridges. 
     791      !! 
     792      !! ** Arguments : 
     793      !! 
     794      !! ** External :  
     795      !! 
     796      !!---------------------------------------------------------------------! 
     797      !! * Arguments 
     798 
    799799      INTEGER :: & 
    800800         ji,jj,  &          ! horizontal indices 
    … …  
    820820         epsi06 = 1.0e-6 
    821821 
    822 !------------------------------------------------------------------------------! 
     822      !------------------------------------------------------------------------------! 
    823823 
    824824      Gstari     = 1.0/Gstar     
    … …  
    833833      krdg (:,:,:)  = 1.0 
    834834 
    835 !     ! Zero out categories with very small areas 
     835      !     ! Zero out categories with very small areas 
    836836      CALL lim_itd_me_zapsmall 
    837837 
    838 !------------------------------------------------------------------------------! 
    839 ! 1) Participation function  
    840 !------------------------------------------------------------------------------! 
     838      !------------------------------------------------------------------------------! 
     839      ! 1) Participation function  
     840      !------------------------------------------------------------------------------! 
    841841 
    842842      ! Compute total area of ice plus open water. 
    … …  
    886886            DO ji = 1, jpi 
    887887               Gsum(ji,jj,jl) = Gsum(ji,jj,jl) * zworka(ji,jj) 
    888             END DO  
     888            END DO 
    889889         END DO 
    890890      END DO 
    891891 
    892 ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) 
    893 !-------------------------------------------------------------------------------------------------- 
    894 ! Compute the participation function athorn; this is analogous to 
    895 ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 
    896 ! area lost from category n due to ridging/closing 
    897 ! athorn(n)   = total area lost due to ridging/closing 
    898 ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).  
    899 ! 
    900 ! The expressions for athorn are found by integrating b(h)g(h) between 
    901 ! the category boundaries. 
    902 !----------------------------------------------------------------- 
     892      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) 
     893      !-------------------------------------------------------------------------------------------------- 
     894      ! Compute the participation function athorn; this is analogous to 
     895      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 
     896      ! area lost from category n due to ridging/closing 
     897      ! athorn(n)   = total area lost due to ridging/closing 
     898      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).  
     899      ! 
     900      ! The expressions for athorn are found by integrating b(h)g(h) between 
     901      ! the category boundaries. 
     902      !----------------------------------------------------------------- 
    903903 
    904904      krdg_index = 1 
    … …  
    906906      IF ( krdg_index .EQ. 0 ) THEN 
    907907 
    908       !--- Linear formulation (Thorndike et al., 1975) 
    909       DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates 
    910          DO jj = 1, jpj  
    911             DO ji = 1, jpi 
    912                IF (Gsum(ji,jj,jl) < Gstar) THEN 
    913                   athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * & 
    914                        (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari) 
    915                ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN 
    916                   athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) *  & 
    917                        (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari) 
    918                ELSE 
    919                   athorn(ji,jj,jl) = 0.0 
    920                ENDIF 
    921             END DO ! ji 
    922          END DO ! jj 
    923       END DO ! jl  
     908         !--- Linear formulation (Thorndike et al., 1975) 
     909         DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates 
     910            DO jj = 1, jpj  
     911               DO ji = 1, jpi 
     912                  IF (Gsum(ji,jj,jl) < Gstar) THEN 
     913                     athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * & 
     914                        (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari) 
     915                  ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN 
     916                     athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) *  & 
     917                        (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari) 
     918                  ELSE 
     919                     athorn(ji,jj,jl) = 0.0 
     920                  ENDIF 
     921               END DO ! ji 
     922            END DO ! jj 
     923         END DO ! jl  
    924924 
    925925      ELSE ! krdg_index = 1 
    926        
    927       !--- Exponential, more stable formulation (Lipscomb et al, 2007) 
    928       ! precompute exponential terms using Gsum as a work array 
    929       zdummy = 1.0 / (1.0-EXP(-astari)) 
    930  
    931       DO jl = -1, jpl 
    932          DO jj = 1, jpj 
    933             DO ji = 1, jpi 
    934                Gsum(ji,jj,jl) = EXP(-Gsum(ji,jj,jl)*astari)*zdummy 
    935             END DO !ji 
    936          END DO !jj 
    937       END DO !jl 
    938  
    939       ! compute athorn 
    940       DO jl = 0, ice_cat_bounds(1,2) 
    941          DO jj = 1, jpj 
    942             DO ji = 1, jpi 
    943                athorn(ji,jj,jl) = Gsum(ji,jj,jl-1) - Gsum(ji,jj,jl) 
    944             END DO !ji 
    945          END DO ! jj 
    946       END DO !jl 
     926 
     927         !--- Exponential, more stable formulation (Lipscomb et al, 2007) 
     928         ! precompute exponential terms using Gsum as a work array 
     929         zdummy = 1.0 / (1.0-EXP(-astari)) 
     930 
     931         DO jl = -1, jpl 
     932            DO jj = 1, jpj 
     933               DO ji = 1, jpi 
     934                  Gsum(ji,jj,jl) = EXP(-Gsum(ji,jj,jl)*astari)*zdummy 
     935               END DO !ji 
     936            END DO !jj 
     937         END DO !jl 
     938 
     939         ! compute athorn 
     940         DO jl = 0, ice_cat_bounds(1,2) 
     941            DO jj = 1, jpj 
     942               DO ji = 1, jpi 
     943                  athorn(ji,jj,jl) = Gsum(ji,jj,jl-1) - Gsum(ji,jj,jl) 
     944               END DO !ji 
     945            END DO ! jj 
     946         END DO !jl 
    947947 
    948948      ENDIF ! krdg_index 
    … …  
    956956                  IF ( athorn(ji,jj,jl) .GT. 0.0 ) THEN 
    957957                     aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - & 
    958                                           hparmeter ) ) + 1.0 ) / 2.0 * &  
    959                                           athorn(ji,jj,jl) 
     958                        hparmeter ) ) + 1.0 ) / 2.0 * &  
     959                        athorn(ji,jj,jl) 
    960960                     araft (ji,jj,jl) = ( TANH ( - Craft * ( ht_i(ji,jj,jl) - & 
    961                                           hparmeter ) ) + 1.0 ) / 2.0 * & 
    962                                           athorn(ji,jj,jl) 
     961                        hparmeter ) ) + 1.0 ) / 2.0 * & 
     962                        athorn(ji,jj,jl) 
    963963                     IF ( araft(ji,jj,jl) .LT. epsi06 ) araft(ji,jj,jl)  = 0.0 
    964964                     aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0) 
    … …  
    982982      IF ( raftswi .EQ. 1 ) THEN 
    983983 
    984       IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi11 ) THEN 
    985          DO jl = 1, jpl 
    986             DO jj = 1, jpj 
    987                DO ji = 1, jpi 
    988                   IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. & 
    989                   epsi11 ) THEN 
    990                      WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 
    991                      WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl 
    992                      WRITE(numout,*) ' lat, lon   : ', gphit(ji,jj), glamt(ji,jj) 
    993                      WRITE(numout,*) ' aridge     : ', aridge(ji,jj,1:jpl) 
    994                      WRITE(numout,*) ' araft      : ', araft(ji,jj,1:jpl) 
    995                      WRITE(numout,*) ' athorn     : ', athorn(ji,jj,1:jpl) 
    996                   ENDIF 
     984         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi11 ) THEN 
     985            DO jl = 1, jpl 
     986               DO jj = 1, jpj 
     987                  DO ji = 1, jpi 
     988                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. & 
     989                        epsi11 ) THEN 
     990                        WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 
     991                        WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl 
     992                        WRITE(numout,*) ' lat, lon   : ', gphit(ji,jj), glamt(ji,jj) 
     993                        WRITE(numout,*) ' aridge     : ', aridge(ji,jj,1:jpl) 
     994                        WRITE(numout,*) ' araft      : ', araft(ji,jj,1:jpl) 
     995                        WRITE(numout,*) ' athorn     : ', athorn(ji,jj,1:jpl) 
     996                     ENDIF 
     997                  END DO 
    997998               END DO 
    998999            END DO 
    999          END DO 
     1000         ENDIF 
     1001 
    10001002      ENDIF 
    10011003 
    1002       ENDIF 
    1003  
    1004 !----------------------------------------------------------------- 
    1005 ! 2) Transfer function 
    1006 !----------------------------------------------------------------- 
    1007 ! Compute max and min ridged ice thickness for each ridging category. 
    1008 ! Assume ridged ice is uniformly distributed between hrmin and hrmax. 
    1009 !  
    1010 ! This parameterization is a modified version of Hibler (1980). 
    1011 ! The mean ridging thickness, hrmean, is proportional to hi^(0.5) 
    1012 !  and for very thick ridging ice must be >= krdgmin*hi 
    1013 ! 
    1014 ! The minimum ridging thickness, hrmin, is equal to 2*hi  
    1015 !  (i.e., rafting) and for very thick ridging ice is 
    1016 !  constrained by hrmin <= (hrmean + hi)/2. 
    1017 !  
    1018 ! The maximum ridging thickness, hrmax, is determined by 
    1019 !  hrmean and hrmin. 
    1020 ! 
    1021 ! These modifications have the effect of reducing the ice strength 
    1022 ! (relative to the Hibler formulation) when very thick ice is 
    1023 ! ridging. 
    1024 ! 
    1025 ! aksum = net area removed/ total area removed 
    1026 ! where total area removed = area of ice that ridges 
    1027 !         net area removed = total area removed - area of new ridges 
    1028 !----------------------------------------------------------------- 
     1004      !----------------------------------------------------------------- 
     1005      ! 2) Transfer function 
     1006      !----------------------------------------------------------------- 
     1007      ! Compute max and min ridged ice thickness for each ridging category. 
     1008      ! Assume ridged ice is uniformly distributed between hrmin and hrmax. 
     1009      !  
     1010      ! This parameterization is a modified version of Hibler (1980). 
     1011      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5) 
     1012      !  and for very thick ridging ice must be >= krdgmin*hi 
     1013      ! 
     1014      ! The minimum ridging thickness, hrmin, is equal to 2*hi  
     1015      !  (i.e., rafting) and for very thick ridging ice is 
     1016      !  constrained by hrmin <= (hrmean + hi)/2. 
     1017      !  
     1018      ! The maximum ridging thickness, hrmax, is determined by 
     1019      !  hrmean and hrmin. 
     1020      ! 
     1021      ! These modifications have the effect of reducing the ice strength 
     1022      ! (relative to the Hibler formulation) when very thick ice is 
     1023      ! ridging. 
     1024      ! 
     1025      ! aksum = net area removed/ total area removed 
     1026      ! where total area removed = area of ice that ridges 
     1027      !         net area removed = total area removed - area of new ridges 
     1028      !----------------------------------------------------------------- 
    10291029 
    10301030      ! Transfer function 
    … …  
    10621062            DO ji = 1, jpi 
    10631063               aksum(ji,jj)    = aksum(ji,jj)                          & 
    1064                        + aridge(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl))    & 
    1065                        + araft (ji,jj,jl) * (1.0 - 1.0/kraft) 
     1064                  + aridge(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl))    & 
     1065                  + araft (ji,jj,jl) * (1.0 - 1.0/kraft) 
    10661066            END DO 
    10671067         END DO 
    10681068      END DO 
    10691069 
    1070       END SUBROUTINE lim_itd_me_ridgeprep 
    1071  
    1072 !=============================================================================== 
    1073  
    1074       SUBROUTINE lim_itd_me_ridgeshift(opning,    closing_gross,       & 
    1075                               msnow_mlt, esnow_mlt) ! (subroutine 4/6) 
    1076  
    1077         !!----------------------------------------------------------------------------- 
    1078         !!                ***  ROUTINE lim_itd_me_icestrength *** 
    1079         !! ** Purpose : 
    1080         !!        This routine shift ridging ice among thickness categories 
    1081         !!                      of ice thickness 
    1082         !! 
    1083         !! ** Method  : 
    1084         !! Remove area, volume, and energy from each ridging category 
    1085         !! and add to thicker ice categories. 
    1086         !! 
    1087         !! ** Arguments : 
    1088         !! 
    1089         !! ** Inputs / Ouputs :  
    1090         !! 
    1091         !! ** External :  
    1092         !! 
     1070   END SUBROUTINE lim_itd_me_ridgeprep 
     1071 
     1072   !=============================================================================== 
     1073 
     1074   SUBROUTINE lim_itd_me_ridgeshift(opning,    closing_gross,       & 
     1075      msnow_mlt, esnow_mlt) ! (subroutine 4/6) 
     1076 
     1077      !!----------------------------------------------------------------------------- 
     1078      !!                ***  ROUTINE lim_itd_me_icestrength *** 
     1079      !! ** Purpose : 
     1080      !!        This routine shift ridging ice among thickness categories 
     1081      !!                      of ice thickness 
     1082      !! 
     1083      !! ** Method  : 
     1084      !! Remove area, volume, and energy from each ridging category 
     1085      !! and add to thicker ice categories. 
     1086      !! 
     1087      !! ** Arguments : 
     1088      !! 
     1089      !! ** Inputs / Ouputs :  
     1090      !! 
     1091      !! ** External :  
     1092      !! 
    10931093 
    10941094      REAL (wp), DIMENSION(jpi,jpj), INTENT(IN)   :: & 
    10951095         opning,         & ! rate of opening due to divergence/shear 
    10961096         closing_gross     ! rate at which area removed, not counting 
    1097                            ! area of new ridges 
     1097      ! area of new ridges 
    10981098 
    10991099      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: & 
    … …  
    11761176      LOGICAL, PARAMETER :: & 
    11771177         l_conservation_check = .true.  ! if true, check conservation  
    1178                                         ! (useful for debugging) 
     1178      ! (useful for debugging) 
    11791179      LOGICAL ::         & 
    11801180         neg_ato_i     , &  ! flag for ato_i(i,j) < -puny 
    … …  
    11871187         zindb              ! switch for the presence of ridge poros or not 
    11881188 
    1189    !---------------------------------------------------------------------------- 
     1189      !---------------------------------------------------------------------------- 
    11901190 
    11911191      ! Conservation check 
    … …  
    12021202      epsi10 = 1.0d-10 
    12031203 
    1204 !------------------------------------------------------------------------------- 
    1205 ! 1) Compute change in open water area due to closing and opening. 
    1206 !------------------------------------------------------------------------------- 
     1204      !------------------------------------------------------------------------------- 
     1205      ! 1) Compute change in open water area due to closing and opening. 
     1206      !------------------------------------------------------------------------------- 
    12071207 
    12081208      neg_ato_i = .false. 
    … …  
    12111211         DO ji = 1, jpi 
    12121212            ato_i(ji,jj) = ato_i(ji,jj)                                   & 
    1213                     - athorn(ji,jj,0)*closing_gross(ji,jj)*rdt_ice        & 
    1214                     + opning(ji,jj)*rdt_ice 
     1213               - athorn(ji,jj,0)*closing_gross(ji,jj)*rdt_ice        & 
     1214               + opning(ji,jj)*rdt_ice 
    12151215            IF (ato_i(ji,jj) .LT. -epsi11) THEN 
    12161216               neg_ato_i = .true. 
    … …  
    12341234      ENDIF                     ! neg_ato_i 
    12351235 
    1236 !----------------------------------------------------------------- 
    1237 ! 2) Save initial state variables 
    1238 !----------------------------------------------------------------- 
     1236      !----------------------------------------------------------------- 
     1237      ! 2) Save initial state variables 
     1238      !----------------------------------------------------------------- 
    12391239 
    12401240      DO jl = 1, jpl 
    … …  
    12521252 
    12531253      esnon_init(:,:,:) = e_s(:,:,1,:) 
    1254              
     1254 
    12551255      DO jl = 1, jpl   
    12561256         DO jk = 1, nlay_i 
    … …  
    12631263      END DO !jl 
    12641264 
    1265 ! 
    1266 !----------------------------------------------------------------- 
    1267 ! 3) Pump everything from ice which is being ridged / rafted 
    1268 !----------------------------------------------------------------- 
    1269 ! Compute the area, volume, and energy of ice ridging in each 
    1270 ! category, along with the area of the resulting ridge. 
     1265      ! 
     1266      !----------------------------------------------------------------- 
     1267      ! 3) Pump everything from ice which is being ridged / rafted 
     1268      !----------------------------------------------------------------- 
     1269      ! Compute the area, volume, and energy of ice ridging in each 
     1270      ! category, along with the area of the resulting ridge. 
    12711271 
    12721272      DO jl1 = 1, jpl !jl1 describes the ridging category 
    12731273 
    1274       !------------------------------------------------ 
    1275       ! 3.1) Identify grid cells with nonzero ridging 
    1276       !------------------------------------------------ 
     1274         !------------------------------------------------ 
     1275         ! 3.1) Identify grid cells with nonzero ridging 
     1276         !------------------------------------------------ 
    12771277 
    12781278         icells = 0 
    … …  
    12801280            DO ji = 1, jpi 
    12811281               IF (aicen_init(ji,jj,jl1) .GT. epsi11 .AND. athorn(ji,jj,jl1) .GT. 0.0       & 
    1282                     .AND. closing_gross(ji,jj) > 0.0) THEN 
     1282                  .AND. closing_gross(ji,jj) > 0.0) THEN 
    12831283                  icells = icells + 1 
    12841284                  indxi(icells) = ji 
    … …  
    12961296            jj = indxj(ij) 
    12971297 
    1298       !-------------------------------------------------------------------- 
    1299       ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2) 
    1300       !-------------------------------------------------------------------- 
     1298            !-------------------------------------------------------------------- 
     1299            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2) 
     1300            !-------------------------------------------------------------------- 
    13011301 
    13021302            ardg1(ji,jj) = aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice 
    … …  
    13101310            oirft2(ji,jj)= oirft1(ji,jj) / kraft 
    13111311 
    1312       !--------------------------------------------------------------- 
    1313       ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1  
    1314       !--------------------------------------------------------------- 
     1312            !--------------------------------------------------------------- 
     1313            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1  
     1314            !--------------------------------------------------------------- 
    13151315 
    13161316            afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging 
    … …  
    13281328            ENDIF 
    13291329 
    1330       !-------------------------------------------------------------------------- 
    1331       ! 3.4) Subtract area, volume, and energy from ridging  
    1332       !     / rafting category n1. 
    1333       !-------------------------------------------------------------------------- 
     1330            !-------------------------------------------------------------------------- 
     1331            ! 3.4) Subtract area, volume, and energy from ridging  
     1332            !     / rafting category n1. 
     1333            !-------------------------------------------------------------------------- 
    13341334            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) /             & 
    1335                            ( 1.0 + ridge_por ) 
     1335               ( 1.0 + ridge_por ) 
    13361336            vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por ) 
    13371337            vsw  (ji,jj) = vrdg1(ji,jj) * ridge_por 
    … …  
    13401340            esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj) 
    13411341            srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) /            & 
    1342                             ( 1. + ridge_por ) 
     1342               ( 1. + ridge_por ) 
    13431343            srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 
    13441344 
    … …  
    13571357            smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj)  - smrft(ji,jj) 
    13581358 
    1359       !----------------------------------------------------------------- 
    1360       ! 3.5) Compute properties of new ridges 
    1361       !----------------------------------------------------------------- 
     1359            !----------------------------------------------------------------- 
     1360            ! 3.5) Compute properties of new ridges 
     1361            !----------------------------------------------------------------- 
    13621362            !------------- 
    13631363            ! Salinity 
    … …  
    13731373            ! salt flux due to ridge creation 
    13741374            fsalt_rpo(ji,jj)  = fsalt_rpo(ji,jj) + &  
    1375             MAX ( zdummy - srdg2(ji,jj) , 0.0 )    & 
    1376             * rhoic / rdt_ice 
     1375               MAX ( zdummy - srdg2(ji,jj) , 0.0 )    & 
     1376               * rhoic / rdt_ice 
    13771377 
    13781378            ! sal times volume for new ridges 
    13791379            srdg2(ji,jj)      = sm_newridge * vrdg2(ji,jj)  
    13801380 
    1381       !------------------------------------             
    1382       ! 3.6 Increment ridging diagnostics 
    1383       !------------------------------------             
    1384  
    1385 !        jl1 looping 1-jpl 
    1386 !           ij looping 1-icells 
     1381            !------------------------------------             
     1382            ! 3.6 Increment ridging diagnostics 
     1383            !------------------------------------             
     1384 
     1385            !        jl1 looping 1-jpl 
     1386            !           ij looping 1-icells 
    13871387 
    13881388            dardg1dt(ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 
    … …  
    13931393            IF (con_i) vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj) 
    13941394 
    1395       !------------------------------------------             
    1396       ! 3.7 Put the snow somewhere in the ocean 
    1397       !------------------------------------------             
    1398  
    1399 !  Place part of the snow lost by ridging into the ocean.  
    1400 !  Note that esnow_mlt < 0; the ocean must cool to melt snow. 
    1401 !  If the ocean temp = Tf already, new ice must grow. 
    1402 !  During the next time step, thermo_rates will determine whether 
    1403 !  the ocean cools or new ice grows. 
    1404 !        jl1 looping 1-jpl 
    1405 !           ij looping 1-icells 
    1406                 
     1395            !------------------------------------------             
     1396            ! 3.7 Put the snow somewhere in the ocean 
     1397            !------------------------------------------             
     1398 
     1399            !  Place part of the snow lost by ridging into the ocean.  
     1400            !  Note that esnow_mlt < 0; the ocean must cool to melt snow. 
     1401            !  If the ocean temp = Tf already, new ice must grow. 
     1402            !  During the next time step, thermo_rates will determine whether 
     1403            !  the ocean cools or new ice grows. 
     1404            !        jl1 looping 1-jpl 
     1405            !           ij looping 1-icells 
     1406 
    14071407            msnow_mlt(ji,jj) = msnow_mlt(ji,jj)                  & 
    1408                            + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg)   & 
    1409                            !rafting included 
    1410                            + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 
     1408               + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg)   & 
     1409                                !rafting included 
     1410               + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 
    14111411 
    14121412            esnow_mlt(ji,jj) = esnow_mlt(ji,jj)                  & 
    1413                            + esrdg(ji,jj)*(1.0-fsnowrdg)         & 
    1414                            !rafting included 
    1415                            + esrft(ji,jj)*(1.0-fsnowrft)           
    1416  
    1417       !----------------------------------------------------------------- 
    1418       ! 3.8 Compute quantities used to apportion ice among categories 
    1419       ! in the n2 loop below 
    1420       !----------------------------------------------------------------- 
    1421  
    1422 !        jl1 looping 1-jpl 
    1423 !           ij looping 1-icells 
     1413               + esrdg(ji,jj)*(1.0-fsnowrdg)         & 
     1414                                !rafting included 
     1415               + esrft(ji,jj)*(1.0-fsnowrft)           
     1416 
     1417            !----------------------------------------------------------------- 
     1418            ! 3.8 Compute quantities used to apportion ice among categories 
     1419            ! in the n2 loop below 
     1420            !----------------------------------------------------------------- 
     1421 
     1422            !        jl1 looping 1-jpl 
     1423            !           ij looping 1-icells 
    14241424 
    14251425            dhr(ji,jj)  = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) 
    14261426            dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1)    & 
    1427                       - hrmin(ji,jj,jl1)   * hrmin(ji,jj,jl1) 
     1427               - hrmin(ji,jj,jl1)   * hrmin(ji,jj,jl1) 
    14281428 
    14291429 
    14301430         END DO                 ! ij 
    14311431 
    1432       !-------------------------------------------------------------------- 
    1433       ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and 
    1434       !      compute ridged ice enthalpy  
    1435       !-------------------------------------------------------------------- 
     1432         !-------------------------------------------------------------------- 
     1433         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and 
     1434         !      compute ridged ice enthalpy  
     1435         !-------------------------------------------------------------------- 
    14361436         DO jk = 1, nlay_i 
    14371437!CDIR NODEP 
    14381438            DO ij = 1, icells 
    1439             ji = indxi(ij) 
    1440             jj = indxj(ij) 
    1441             ! heat content of ridged ice 
    1442             erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / &  
    1443                                    ( 1.0 + ridge_por )  
    1444             eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 
    1445             e_i(ji,jj,jk,jl1)    = e_i(ji,jj,jk,jl1)             & 
    1446                                         - erdg1(ji,jj,jk)        & 
    1447                                         - eirft(ji,jj,jk) 
    1448             ! sea water heat content 
    1449             ztmelts          = - tmut * sss_m(ji,jj) + rtt 
    1450             ! heat content per unit volume 
    1451             zdummy0          = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj) 
    1452  
    1453             ! corrected sea water salinity 
    1454             zindb  = MAX( 0.0, SIGN( 1.0, vsw(ji,jj) - zeps ) ) 
    1455             zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / & 
    1456                      MAX( ridge_por * vsw(ji,jj), zeps ) 
    1457  
    1458             ztmelts          = - tmut * zdummy + rtt 
    1459             ersw(ji,jj,jk)   = - rcp * ( ztmelts - rtt ) * vsw(ji,jj) 
    1460  
    1461             ! heat flux 
    1462             fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / & 
    1463                                      rdt_ice 
    1464  
    1465             ! Correct dimensions to avoid big values 
    1466             ersw(ji,jj,jk)   = ersw(ji,jj,jk) / 1.0d+09 
    1467  
    1468             ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    1469             ersw(ji,jj,jk)   = ersw(ji,jj,jk) * & 
    1470                                area(ji,jj) * vsw(ji,jj) / & 
    1471                                nlay_i 
    1472  
    1473             erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
     1439               ji = indxi(ij) 
     1440               jj = indxj(ij) 
     1441               ! heat content of ridged ice 
     1442               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / &  
     1443                  ( 1.0 + ridge_por )  
     1444               eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 
     1445               e_i(ji,jj,jk,jl1)    = e_i(ji,jj,jk,jl1)             & 
     1446                  - erdg1(ji,jj,jk)        & 
     1447                  - eirft(ji,jj,jk) 
     1448               ! sea water heat content 
     1449               ztmelts          = - tmut * sss_m(ji,jj) + rtt 
     1450               ! heat content per unit volume 
     1451               zdummy0          = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj) 
     1452 
     1453               ! corrected sea water salinity 
     1454               zindb  = MAX( 0.0, SIGN( 1.0, vsw(ji,jj) - zeps ) ) 
     1455               zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / & 
     1456                  MAX( ridge_por * vsw(ji,jj), zeps ) 
     1457 
     1458               ztmelts          = - tmut * zdummy + rtt 
     1459               ersw(ji,jj,jk)   = - rcp * ( ztmelts - rtt ) * vsw(ji,jj) 
     1460 
     1461               ! heat flux 
     1462               fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / & 
     1463                  rdt_ice 
     1464 
     1465               ! Correct dimensions to avoid big values 
     1466               ersw(ji,jj,jk)   = ersw(ji,jj,jk) / 1.0d+09 
     1467 
     1468               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
     1469               ersw(ji,jj,jk)   = ersw(ji,jj,jk) * & 
     1470                  area(ji,jj) * vsw(ji,jj) / & 
     1471                  nlay_i 
     1472 
     1473               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
    14741474            END DO ! ij 
    14751475         END DO !jk 
    … …  
    14831483                  jj = indxj(ij) 
    14841484                  eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - & 
    1485                   erdg1(ji,jj,jk) 
     1485                     erdg1(ji,jj,jk) 
    14861486               END DO ! ij 
    14871487            END DO !jk 
    … …  
    14971497                  WRITE(numout,*) ' ardg > a_i' 
    14981498                  WRITE(numout,*) ' ardg, aicen_init : ', & 
    1499                        ardg1(ji,jj), aicen_init(ji,jj,jl1) 
     1499                     ardg1(ji,jj), aicen_init(ji,jj,jl1) 
    15001500               ENDIF            ! afrac > 1 + puny 
    15011501            ENDDO               ! if 
    … …  
    15101510                  WRITE(numout,*) ' arft > a_i' 
    15111511                  WRITE(numout,*) ' arft, aicen_init : ', & 
    1512                        arft1(ji,jj), aicen_init(ji,jj,jl1) 
     1512                     arft1(ji,jj), aicen_init(ji,jj,jl1) 
    15131513               ENDIF            ! afrft > 1 + puny 
    15141514            ENDDO               ! if 
    15151515         ENDIF                  ! large_afrft 
    15161516 
    1517 !------------------------------------------------------------------------------- 
    1518 ! 4) Add area, volume, and energy of new ridge to each category jl2 
    1519 !------------------------------------------------------------------------------- 
    1520 !        jl1 looping 1-jpl 
     1517         !------------------------------------------------------------------------------- 
     1518         ! 4) Add area, volume, and energy of new ridge to each category jl2 
     1519         !------------------------------------------------------------------------------- 
     1520         !        jl1 looping 1-jpl 
    15211521         DO jl2  = ice_cat_bounds(1,1), ice_cat_bounds(1,2)  
    1522          ! over categories to which ridged ice is transferred 
     1522            ! over categories to which ridged ice is transferred 
    15231523!CDIR NODEP 
    15241524            DO ij = 1, icells 
    … …  
    15311531 
    15321532               IF (hrmin(ji,jj,jl1) .GE. hi_max(jl2) .OR.        & 
    1533                    hrmax(ji,jj,jl1) .LE. hi_max(jl2-1)) THEN 
     1533                  hrmax(ji,jj,jl1) .LE. hi_max(jl2-1)) THEN 
    15341534                  hL = 0.0 
    15351535                  hR = 0.0 
    … …  
    15461546               v_i(ji,jj,jl2)    = v_i(ji,jj,jl2) + fvol(ji,jj) * vrdg2(ji,jj) 
    15471547               v_s(ji,jj,jl2)    = v_s(ji,jj,jl2)                             & 
    1548                                  + fvol(ji,jj) * vsrdg(ji,jj) * fsnowrdg 
     1548                  + fvol(ji,jj) * vsrdg(ji,jj) * fsnowrdg 
    15491549               e_s(ji,jj,1,jl2)  = e_s(ji,jj,1,jl2)                           & 
    1550                                  + fvol(ji,jj) * esrdg(ji,jj) * fsnowrdg 
     1550                  + fvol(ji,jj) * esrdg(ji,jj) * fsnowrdg 
    15511551               smv_i(ji,jj,jl2)  = smv_i(ji,jj,jl2) + fvol(ji,jj) * srdg2(ji,jj) 
    15521552               oa_i(ji,jj,jl2)   = oa_i(ji,jj,jl2)  + farea * oirdg2(ji,jj) 
    … …  
    15611561                  jj = indxj(ij) 
    15621562                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2)          & 
    1563                                     + fvol(ji,jj)*erdg2(ji,jj,jk) 
     1563                     + fvol(ji,jj)*erdg2(ji,jj,jk) 
    15641564               END DO           ! ij 
    15651565            END DO !jk 
    … …  
    15761576               ! Compute the fraction of rafted ice area and volume going to  
    15771577               ! thickness category jl2, transfer area, volume, and energy accordingly. 
    1578              
     1578 
    15791579               IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND.        & 
    1580                    hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 
     1580                  hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 
    15811581                  a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + arft2(ji,jj) 
    15821582                  v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + virft(ji,jj) 
    15831583                  v_s(ji,jj,jl2) = v_s(ji,jj,jl2)                   & 
    1584                                  + vsrft(ji,jj)*fsnowrft 
     1584                     + vsrft(ji,jj)*fsnowrft 
    15851585                  e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2)                   & 
    1586                                  + esrft(ji,jj)*fsnowrft 
     1586                     + esrft(ji,jj)*fsnowrft 
    15871587                  smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2)                 & 
    1588                                  + smrft(ji,jj)     
     1588                     + smrft(ji,jj)     
    15891589                  oa_i(ji,jj,jl2)  = oa_i(ji,jj,jl2)                  & 
    1590                                    + oirft2(ji,jj)     
     1590                     + oirft2(ji,jj)     
    15911591               ENDIF ! hraft 
    15921592 
    … …  
    16001600                  jj = indxj(ij) 
    16011601                  IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND.        & 
    1602                       hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 
     1602                     hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 
    16031603                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2)             & 
    1604                                        + eirft(ji,jj,jk) 
     1604                        + eirft(ji,jj,jk) 
    16051605                  ENDIF 
    16061606               END DO           ! ij 
    … …  
    16281628   END SUBROUTINE lim_itd_me_ridgeshift 
    16291629 
    1630 !============================================================================== 
     1630   !============================================================================== 
    16311631 
    16321632   SUBROUTINE lim_itd_me_asumr !(subroutine 5/6) 
    16331633 
    1634         !!----------------------------------------------------------------------------- 
    1635         !!                ***  ROUTINE lim_itd_me_asumr *** 
    1636         !! ** Purpose : 
    1637         !!        This routine finds total fractional area 
    1638         !! 
    1639         !! ** Method  : 
    1640         !! Find the total area of ice plus open water in each grid cell. 
    1641         !! 
    1642         !! This is similar to the aggregate_area subroutine except that the 
    1643         !! total area can be greater than 1, so the open water area is  
    1644         !! included in the sum instead of being computed as a residual.  
    1645         !! 
    1646         !! ** Arguments : 
     1634      !!----------------------------------------------------------------------------- 
     1635      !!                ***  ROUTINE lim_itd_me_asumr *** 
     1636      !! ** Purpose : 
     1637      !!        This routine finds total fractional area 
     1638      !! 
     1639      !! ** Method  : 
     1640      !! Find the total area of ice plus open water in each grid cell. 
     1641      !! 
     1642      !! This is similar to the aggregate_area subroutine except that the 
     1643      !! total area can be greater than 1, so the open water area is  
     1644      !! included in the sum instead of being computed as a residual.  
     1645      !! 
     1646      !! ** Arguments : 
    16471647 
    16481648      INTEGER :: ji, jj, jl 
    … …  
    16721672   END SUBROUTINE lim_itd_me_asumr 
    16731673 
    1674 !============================================================================== 
     1674   !============================================================================== 
    16751675 
    16761676   SUBROUTINE lim_itd_me_init ! (subroutine 6/6) 
    … …  
    16911691      !!------------------------------------------------------------------- 
    16921692      NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,&  
    1693                             Gstar, astar,                                & 
    1694                             Hstar, raftswi, hparmeter, Craft, ridge_por, & 
    1695                             sal_max_ridge,  partfun_swi, transfun_swi,   & 
    1696                             brinstren_swi 
     1693         Gstar, astar,                                & 
     1694         Hstar, raftswi, hparmeter, Craft, ridge_por, & 
     1695         sal_max_ridge,  partfun_swi, transfun_swi,   & 
     1696         brinstren_swi 
    16971697      !!------------------------------------------------------------------- 
    16981698 
    … …  
    17251725   END SUBROUTINE lim_itd_me_init 
    17261726 
    1727 !============================================================================== 
     1727   !============================================================================== 
    17281728 
    17291729   SUBROUTINE lim_itd_me_zapsmall 
    … …  
    17431743 
    17441744      INTEGER ::   & 
    1745            ji,jj,  & ! horizontal indices 
    1746            jl,     & ! ice category index 
    1747            jk,     & ! ice layer index 
    1748 !           ij,     &   ! combined i/j horizontal index 
    1749            icells      ! number of cells with ice to zap 
    1750  
    1751 !      INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 
    1752 !           indxi,  & ! compressed indices for i/j directions 
    1753 !           indxj 
     1745         ji,jj,  & ! horizontal indices 
     1746         jl,     & ! ice category index 
     1747         jk,     & ! ice layer index 
     1748         !           ij,     &   ! combined i/j horizontal index 
     1749         icells      ! number of cells with ice to zap 
     1750 
     1751      !      INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 
     1752      !           indxi,  & ! compressed indices for i/j directions 
     1753      !           indxj 
    17541754 
    17551755      INTEGER, DIMENSION(jpi,jpj) :: zmask 
    … …  
    17571757 
    17581758      REAL(wp) :: & 
    1759            xtmp      ! temporary variable 
     1759         xtmp      ! temporary variable 
    17601760 
    17611761      DO jl = 1, jpl 
    17621762 
    1763       !----------------------------------------------------------------- 
    1764       ! Count categories to be zapped. 
    1765       ! Abort model in case of negative area. 
    1766       !----------------------------------------------------------------- 
    1767          IF( MINVAL(a_i(:,:,jl)) .LT. -epsi11 ) THEN 
     1763         !----------------------------------------------------------------- 
     1764         ! Count categories to be zapped. 
     1765         ! Abort model in case of negative area. 
     1766         !----------------------------------------------------------------- 
     1767         IF( MINVAL(a_i(:,:,jl)) .LT. -epsi11 .AND. ln_nicep ) THEN 
    17681768            DO jj = 1, jpj 
    17691769               DO ji = 1, jpi 
    … …  
    17741774                  ENDIF 
    17751775               END DO 
    1776             END DO  
     1776            END DO 
    17771777         ENDIF 
    1778   
    1779        icells = 0 
    1780        zmask = 0.e0 
    1781        DO jj = 1, jpj 
    1782          DO ji = 1, jpi 
    1783             IF ( ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0.0)       & 
    1784                                          .OR.                                         & 
    1785                      ( a_i(ji,jj,jl) .GT. 0.0     .AND. a_i(ji,jj,jl) .LE. 1.0e-11 )  & 
    1786                                          .OR.                                         & 
    1787                                          !new line 
    1788                      ( v_i(ji,jj,jl) .EQ. 0.0     .AND. a_i(ji,jj,jl) .GT. 0.0    )   & 
    1789                                          .OR.                                         & 
    1790                      ( v_i(ji,jj,jl) .GT. 0.0     .AND. v_i(ji,jj,jl) .LT. 1.e-12 ) ) THEN 
    1791                 zmask(ji,jj) = 1 
    1792             ENDIF 
     1778 
     1779         icells = 0 
     1780         zmask = 0.e0 
     1781         DO jj = 1, jpj 
     1782            DO ji = 1, jpi 
     1783               IF ( ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0.0)       & 
     1784                  .OR.                                         & 
     1785                  ( a_i(ji,jj,jl) .GT. 0.0     .AND. a_i(ji,jj,jl) .LE. 1.0e-11 )  & 
     1786                  .OR.                                         & 
     1787                                !new line 
     1788                  ( v_i(ji,jj,jl) .EQ. 0.0     .AND. a_i(ji,jj,jl) .GT. 0.0    )   & 
     1789                  .OR.                                         & 
     1790                  ( v_i(ji,jj,jl) .GT. 0.0     .AND. v_i(ji,jj,jl) .LT. 1.e-12 ) ) THEN 
     1791                  zmask(ji,jj) = 1 
     1792               ENDIF 
     1793            END DO 
    17931794         END DO 
    1794          END DO 
    1795          WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean ' 
    1796  
    1797       !----------------------------------------------------------------- 
    1798       ! Zap ice energy and use ocean heat to melt ice 
    1799       !----------------------------------------------------------------- 
     1795         IF( ln_nicep ) WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean ' 
     1796 
     1797         !----------------------------------------------------------------- 
     1798         ! Zap ice energy and use ocean heat to melt ice 
     1799         !----------------------------------------------------------------- 
    18001800 
    18011801         DO jk = 1, nlay_i 
    … …  
    18031803               DO ji = 1 , jpi 
    18041804 
    1805                xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice 
    1806                xtmp = xtmp * unit_fac 
    1807 !              fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
    1808                e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1 - zmask(ji,jj) ) 
     1805                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice 
     1806                  xtmp = xtmp * unit_fac 
     1807                  !              fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
     1808                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1 - zmask(ji,jj) ) 
    18091809               END DO           ! ji 
    18101810            END DO           ! jj 
    … …  
    18141814            DO ji = 1 , jpi 
    18151815 
    1816       !----------------------------------------------------------------- 
    1817       ! Zap snow energy and use ocean heat to melt snow 
    1818       !----------------------------------------------------------------- 
    1819  
    1820 !           xtmp = esnon(i,j,n) / dt ! < 0 
    1821 !           fhnet(i,j)      = fhnet(i,j)      + xtmp 
    1822 !           fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp 
    1823             ! xtmp is greater than 0 
    1824             ! fluxes are positive to the ocean 
    1825             ! here the flux has to be negative for the ocean 
    1826             xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice 
    1827 !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
    1828  
    1829             xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice !RB   ??????? 
    1830  
    1831             t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) ) 
    1832  
    1833       !----------------------------------------------------------------- 
    1834       ! zap ice and snow volume, add water and salt to ocean 
    1835       !----------------------------------------------------------------- 
    1836  
    1837 !           xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt 
    1838 !           fresh(i,j)      = fresh(i,j)      + xtmp 
    1839 !           fresh_hist(i,j) = fresh_hist(i,j) + xtmp 
    1840  
    1841 !           fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj)                  ) * &  
    1842 !                               rhosn * v_s(ji,jj,jl) / rdt_ice 
    1843  
    1844 !           fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * &  
    1845 !                               rhoic * v_i(ji,jj,jl) / rdt_ice 
    1846  
    1847 !           emps(i,j)      = emps(i,j)      + xtmp 
    1848 !           fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp 
    1849  
    1850             ato_i(ji,jj)   = a_i(ji,jj,jl)  * zmask(ji,jj) + ato_i(ji,jj) 
    1851             a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1852             v_i(ji,jj,jl)  = v_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1853             v_s(ji,jj,jl)  = v_s(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1854             t_su(ji,jj,jl) = t_su(ji,jj,jl) * (1 -zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 
    1855             oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    1856             smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
     1816               !----------------------------------------------------------------- 
     1817               ! Zap snow energy and use ocean heat to melt snow 
     1818               !----------------------------------------------------------------- 
     1819 
     1820               !           xtmp = esnon(i,j,n) / dt ! < 0 
     1821               !           fhnet(i,j)      = fhnet(i,j)      + xtmp 
     1822               !           fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp 
     1823               ! xtmp is greater than 0 
     1824               ! fluxes are positive to the ocean 
     1825               ! here the flux has to be negative for the ocean 
     1826               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice 
     1827               !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 
     1828 
     1829               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice !RB   ??????? 
     1830 
     1831               t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) ) 
     1832 
     1833               !----------------------------------------------------------------- 
     1834               ! zap ice and snow volume, add water and salt to ocean 
     1835               !----------------------------------------------------------------- 
     1836 
     1837               !           xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt 
     1838               !           fresh(i,j)      = fresh(i,j)      + xtmp 
     1839               !           fresh_hist(i,j) = fresh_hist(i,j) + xtmp 
     1840 
     1841               !           fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj)                  ) * &  
     1842               !                               rhosn * v_s(ji,jj,jl) / rdt_ice 
     1843 
     1844               !           fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * &  
     1845               !                               rhoic * v_i(ji,jj,jl) / rdt_ice 
     1846 
     1847               !           emps(i,j)      = emps(i,j)      + xtmp 
     1848               !           fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp 
     1849 
     1850               ato_i(ji,jj)   = a_i(ji,jj,jl)  * zmask(ji,jj) + ato_i(ji,jj) 
     1851               a_i(ji,jj,jl)  = a_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
     1852               v_i(ji,jj,jl)  = v_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
     1853               v_s(ji,jj,jl)  = v_s(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
     1854               t_su(ji,jj,jl) = t_su(ji,jj,jl) * (1 -zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 
     1855               oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
     1856               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 
    18571857 
    18581858            END DO                 ! ji 

  • trunk/NEMO/LIM_SRC_3/limitd_th.F90

    r869 r921  
    2828   USE prtctl           ! Print control 
    2929   USE lib_mpp  
    30   
     30 
    3131   IMPLICIT NONE 
    3232   PRIVATE 
    … …  
    5151   !!---------------------------------------------------------------------- 
    5252 
    53 !!---------------------------------------------------------------------------------------------- 
    54 !!---------------------------------------------------------------------------------------------- 
     53   !!---------------------------------------------------------------------------------------------- 
     54   !!---------------------------------------------------------------------------------------------- 
    5555 
    5656CONTAINS 
    5757 
    58    SUBROUTINE lim_itd_th 
    59         !!------------------------------------------------------------------ 
    60         !!                ***  ROUTINE lim_itd_th *** 
    61         !! ** Purpose : 
    62         !!        This routine computes the thermodynamics of ice thickness 
    63         !!         distribution 
    64         !! ** Method  : 
    65         !! 
    66         !! ** Arguments : 
    67         !!           kideb , kiut : Starting and ending points on which the  
    68         !!                         the computation is applied 
    69         !! 
    70         !! ** Inputs / Ouputs : (global commons) 
    71         !! 
    72         !! ** External :  
    73         !! 
    74         !! ** References : 
    75         !! 
    76         !! ** History : 
    77         !!           (12-2005) Martin Vancoppenolle  
    78         !! 
    79         !!------------------------------------------------------------------ 
    80         !! * Arguments 
    81  
    82        !! * Local variables 
    83        INTEGER ::   jl, ja,   &   ! ice category, layers 
    84                     jm,       &   ! ice types    dummy loop index 
    85                     jbnd1,    & 
    86                     jbnd2 
    87  
    88        REAL(wp)  ::           &  ! constant values 
    89           zeps      =  1.0e-10, & 
    90           epsi10    =  1.0e-10 
    91  
    92 !!-- End of declarations 
    93 !!---------------------------------------------------------------------------------------------- 
    94  
    95        IF (lwp) THEN 
    96           WRITE(numout,*) 
    97           WRITE(numout,*) 'lim_itd_th  : Thermodynamics of the ice thickness distribution' 
    98           WRITE(numout,*) '~~~~~~~~~~~' 
    99        ENDIF 
    100  
    101 !------------------------------------------------------------------------------| 
    102 !  1) Transport of ice between thickness categories.                           | 
    103 !------------------------------------------------------------------------------| 
     58   SUBROUTINE lim_itd_th( kt ) 
     59      !!------------------------------------------------------------------ 
     60      !!                ***  ROUTINE lim_itd_th *** 
     61      !! ** Purpose : 
     62      !!        This routine computes the thermodynamics of ice thickness 
     63      !!         distribution 
     64      !! ** Method  : 
     65      !! 
     66      !! ** Arguments : 
     67      !!           kideb , kiut : Starting and ending points on which the  
     68      !!                         the computation is applied 
     69      !! 
     70      !! ** Inputs / Ouputs : (global commons) 
     71      !! 
     72      !! ** External :  
     73      !! 
     74      !! ** References : 
     75      !! 
     76      !! ** History : 
     77      !!           (12-2005) Martin Vancoppenolle  
     78      !! 
     79      !!------------------------------------------------------------------ 
     80      !! * Arguments 
     81      INTEGER, INTENT(in) :: kt 
     82      !! * Local variables 
     83      INTEGER ::   jl, ja,   &   ! ice category, layers 
     84         jm,       &   ! ice types    dummy loop index 
     85         jbnd1,    & 
     86         jbnd2 
     87 
     88      REAL(wp)  ::           &  ! constant values 
     89         zeps      =  1.0e-10, & 
     90         epsi10    =  1.0e-10 
     91 
     92      IF( kt == nit000 .AND. lwp ) THEN 
     93         WRITE(numout,*) 
     94         WRITE(numout,*) 'lim_itd_th  : Thermodynamics of the ice thickness distribution' 
     95         WRITE(numout,*) '~~~~~~~~~~~' 
     96      ENDIF 
     97 
     98      !------------------------------------------------------------------------------| 
     99      !  1) Transport of ice between thickness categories.                           | 
     100      !------------------------------------------------------------------------------| 
    104101      ! Given thermodynamic growth rates, transport ice between 
    105102      ! thickness categories. 
    … …  
    107104         jbnd1 = ice_cat_bounds(jm,1) 
    108105         jbnd2 = ice_cat_bounds(jm,2) 
    109          IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_rem(jbnd1, jbnd2, jm) 
     106         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt ) 
    110107      END DO 
    111108 
    … …  
    113110      CALL lim_var_agg(1) 
    114111 
    115 !------------------------------------------------------------------------------| 
    116 !  3) Add frazil ice growing in leads. 
    117 !------------------------------------------------------------------------------| 
     112      !------------------------------------------------------------------------------| 
     113      !  3) Add frazil ice growing in leads. 
     114      !------------------------------------------------------------------------------| 
    118115 
    119116      CALL lim_thd_lac 
    120117      CALL lim_var_glo2eqv ! only for info 
    121118 
    122 !---------------------------------------------------------------------------------------- 
    123 !  4) Computation of trend terms and get back to old values       
    124 !---------------------------------------------------------------------------------------- 
     119      !---------------------------------------------------------------------------------------- 
     120      !  4) Computation of trend terms and get back to old values       
     121      !---------------------------------------------------------------------------------------- 
    125122 
    126123      !- Trend terms 
    … …  
    133130      d_smv_i_thd(:,:,:) = 0.0 
    134131      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 
    135       d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     132         d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
    136133 
    137134      IF(ln_ctl) THEN   ! Control print 
    … …  
    166163         END DO 
    167164      ENDIF 
    168        
     165 
    169166      !- Recover Old values 
    170167      a_i(:,:,:)         = old_a_i (:,:,:) 
    … …  
    175172 
    176173      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 
    177       smv_i(:,:,:)       = old_smv_i (:,:,:) 
    178  
    179  
    180       END SUBROUTINE lim_itd_th 
    181  
    182 !!---------------------------------------------------------------------------------------------- 
    183 !!---------------------------------------------------------------------------------------------- 
    184  
    185     SUBROUTINE lim_itd_th_rem(klbnd,kubnd,ntyp) 
    186         !!------------------------------------------------------------------ 
    187         !!                ***  ROUTINE lim_itd_th_rem *** 
    188         !! ** Purpose : 
    189         !!        This routine computes the redistribution of ice thickness 
    190         !!        after thermodynamic growth of ice thickness 
    191         !! 
    192         !! ** Method  : Linear remapping  
    193         !! 
    194         !! ** Arguments : 
    195         !!           klbnd, kubnd : Starting and ending category index on which the  
    196         !!                         the computation is applied 
    197         !! 
    198         !! ** Inputs / Ouputs : (global commons) 
    199         !! 
    200         !! ** External :  
    201         !! 
    202         !! ** References : W.H. Lipscomb, JGR 2001 
    203         !! 
    204         !! ** History : 
    205         !!           largely inspired from CICE (c) W. H. Lipscomb and E.C. Hunke 
    206         !!  
    207         !!           (01-2006) Martin Vancoppenolle, UCL-ASTR, translation from 
    208         !!                     CICE 
    209         !!           (06-2006) Adaptation to include salt, age and types 
    210         !!           (04-2007) Mass conservation checked 
    211         !!------------------------------------------------------------------ 
    212         !! * Arguments 
    213  
    214        INTEGER , INTENT (IN) ::  & 
    215           klbnd ,  &  ! Start thickness category index point 
    216           kubnd ,  &  ! End point on which the  the computation is applied 
    217           ntyp        ! Number of the type used 
    218  
    219        !! * Local variables 
    220        INTEGER ::   ji,       &   ! spatial dummy loop index 
    221                     jj,       &   ! spatial dummy loop index 
    222                     jl,       &   ! ice category dummy loop index 
    223                     zji, zjj, &   ! dummy indices used when changing coordinates 
    224                     nd            ! used for thickness categories 
    225  
    226        INTEGER , DIMENSION(jpi,jpj,jpl-1) :: &  
    227                     zdonor        ! donor category index 
    228    
    229        REAL(wp)  ::           &   ! constant values 
    230           zeps      =  1.0e-10 
    231  
    232        REAL(wp)  ::           &  ! constant values for ice enthalpy 
    233           zindb     ,         & 
    234           zareamin  ,         &  ! minimum tolerated area in a thickness category 
    235           zwk1, zwk2,         &  ! all the following are dummy arguments 
    236           zx1, zx2, zx3,      &  ! 
    237           zetamin   ,         &  ! minimum value of eta 
    238           zetamax   ,         &  ! maximum value of eta 
    239           zdh0      ,         &  !  
    240           zda0      ,         &  ! 
    241           zdamax    ,         &  ! 
    242           zhimin 
    243  
    244        REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 
    245           zdhice           ,  &  ! ice thickness increment 
    246           g0               ,  &  ! coefficients for fitting the line of the ITD 
    247           g1               ,  &  ! coefficients for fitting the line of the ITD 
    248           hL               ,  &  ! left boundary for the ITD for each thickness 
    249           hR               ,  &  ! left boundary for the ITD for each thickness 
    250           zht_i_o          ,  &  ! old ice thickness 
    251           dummy_es 
    252  
    253        REAL(wp), DIMENSION(jpi,jpj,jpl-1) :: & 
    254           zdaice           ,  &  ! local increment of ice area  
    255           zdvice                 ! local increment of ice volume 
    256  
    257        REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: & 
    258           zhbnew                 ! new boundaries of ice categories 
    259  
    260        REAL(wp), DIMENSION(jpi,jpj) :: & 
    261           zhb0, zhb1             ! category boundaries for thinnes categories 
    262  
    263        REAL, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 
    264           zvetamin, zvetamax     ! maximum values for etas 
    265   
    266        INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 
    267           nind_i      ,  &  ! compressed indices for i/j directions 
    268           nind_j 
    269  
    270        INTEGER :: & 
    271           nbrem             ! number of cells with ice to transfer 
    272  
    273        LOGICAL, DIMENSION(jpi,jpj) ::   &  !: 
    274           zremap_flag             ! compute remapping or not ???? 
    275         
    276        REAL(wp)  ::           &  ! constant values for ice enthalpy 
    277           zslope                 ! used to compute local thermodynamic "speeds" 
    278  
    279        REAL (wp), DIMENSION(jpi,jpj) :: &  !  
    280                vt_i_init, vt_i_final,   &  !  ice volume summed over categories 
    281                vt_s_init, vt_s_final,   &  !  snow volume summed over categories 
    282                et_i_init, et_i_final,   &  !  ice energy summed over categories 
    283                et_s_init, et_s_final       !  snow energy summed over categories 
    284  
    285        CHARACTER (len = 15) :: fieldid 
    286         
    287 !!-- End of declarations 
    288 !!---------------------------------------------------------------------------------------------- 
     174         smv_i(:,:,:)       = old_smv_i (:,:,:) 
     175 
     176   END SUBROUTINE lim_itd_th 
     177   ! 
     178 
     179   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp, kt ) 
     180      !!------------------------------------------------------------------ 
     181      !!                ***  ROUTINE lim_itd_th_rem *** 
     182      !! ** Purpose : 
     183      !!        This routine computes the redistribution of ice thickness 
     184      !!        after thermodynamic growth of ice thickness 
     185      !! 
     186      !! ** Method  : Linear remapping  
     187      !! 
     188      !! ** Arguments : 
     189      !!           klbnd, kubnd : Starting and ending category index on which the  
     190      !!                         the computation is applied 
     191      !! 
     192      !! ** Inputs / Ouputs : (global commons) 
     193      !! 
     194      !! ** External :  
     195      !! 
     196      !! ** References : W.H. Lipscomb, JGR 2001 
     197      !! 
     198      !! ** History : 
     199      !!           largely inspired from CICE (c) W. H. Lipscomb and E.C. Hunke 
     200      !!  
     201      !!           (01-2006) Martin Vancoppenolle, UCL-ASTR, translation from 
     202      !!                     CICE 
     203      !!           (06-2006) Adaptation to include salt, age and types 
     204      !!           (04-2007) Mass conservation checked 
     205      !!------------------------------------------------------------------ 
     206      !! * Arguments 
     207 
     208      INTEGER , INTENT (IN) ::  & 
     209         klbnd ,  &  ! Start thickness category index point 
     210         kubnd ,  &  ! End point on which the  the computation is applied 
     211         ntyp  ,  &  ! Number of the type used 
     212         kt          ! Ocean time step  
     213 
     214      !! * Local variables 
     215      INTEGER ::   ji,       &   ! spatial dummy loop index 
     216         jj,       &   ! spatial dummy loop index 
     217         jl,       &   ! ice category dummy loop index 
     218         zji, zjj, &   ! dummy indices used when changing coordinates 
     219         nd            ! used for thickness categories 
     220 
     221      INTEGER , DIMENSION(jpi,jpj,jpl-1) :: &  
     222         zdonor        ! donor category index 
     223 
     224      REAL(wp)  ::           &   ! constant values 
     225         zeps      =  1.0e-10 
     226 
     227      REAL(wp)  ::           &  ! constant values for ice enthalpy 
     228         zindb     ,         & 
     229         zareamin  ,         &  ! minimum tolerated area in a thickness category 
     230         zwk1, zwk2,         &  ! all the following are dummy arguments 
     231         zx1, zx2, zx3,      &  ! 
     232         zetamin   ,         &  ! minimum value of eta 
     233         zetamax   ,         &  ! maximum value of eta 
     234         zdh0      ,         &  !  
     235         zda0      ,         &  ! 
     236         zdamax    ,         &  ! 
     237         zhimin 
     238 
     239      REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 
     240         zdhice           ,  &  ! ice thickness increment 
     241         g0               ,  &  ! coefficients for fitting the line of the ITD 
     242         g1               ,  &  ! coefficients for fitting the line of the ITD 
     243         hL               ,  &  ! left boundary for the ITD for each thickness 
     244         hR               ,  &  ! left boundary for the ITD for each thickness 
     245         zht_i_o          ,  &  ! old ice thickness 
     246         dummy_es 
     247 
     248      REAL(wp), DIMENSION(jpi,jpj,jpl-1) :: & 
     249         zdaice           ,  &  ! local increment of ice area  
     250         zdvice                 ! local increment of ice volume 
     251 
     252      REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: & 
     253         zhbnew                 ! new boundaries of ice categories 
     254 
     255      REAL(wp), DIMENSION(jpi,jpj) :: & 
     256         zhb0, zhb1             ! category boundaries for thinnes categories 
     257 
     258      REAL, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 
     259         zvetamin, zvetamax     ! maximum values for etas 
     260 
     261      INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 
     262         nind_i      ,  &  ! compressed indices for i/j directions 
     263         nind_j 
     264 
     265      INTEGER :: & 
     266         nbrem             ! number of cells with ice to transfer 
     267 
     268      LOGICAL, DIMENSION(jpi,jpj) ::   &  !: 
     269         zremap_flag             ! compute remapping or not ???? 
     270 
     271      REAL(wp)  ::           &  ! constant values for ice enthalpy 
     272         zslope                 ! used to compute local thermodynamic "speeds" 
     273 
     274      REAL (wp), DIMENSION(jpi,jpj) :: &  !  
     275         vt_i_init, vt_i_final,   &  !  ice volume summed over categories 
     276         vt_s_init, vt_s_final,   &  !  snow volume summed over categories 
     277         et_i_init, et_i_final,   &  !  ice energy summed over categories 
     278         et_s_init, et_s_final       !  snow energy summed over categories 
     279 
     280      CHARACTER (len = 15) :: fieldid 
     281 
     282      !!-- End of declarations 
     283      !!---------------------------------------------------------------------------------------------- 
    289284      zhimin = 0.1      !minimum ice thickness tolerated by the model 
    290285      zareamin = zeps   !minimum area in thickness categories tolerated by the conceptors of the model 
    291286 
    292 !!---------------------------------------------------------------------------------------------- 
    293 !! 0) Conservation checkand changes in each ice category 
    294 !!---------------------------------------------------------------------------------------------- 
     287      !!---------------------------------------------------------------------------------------------- 
     288      !! 0) Conservation checkand changes in each ice category 
     289      !!---------------------------------------------------------------------------------------------- 
    295290      IF ( con_i ) THEN 
    296291         CALL lim_column_sum (jpl,   v_i, vt_i_init) 
    … …  
    300295         CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init) 
    301296      ENDIF 
    302   
    303 !!---------------------------------------------------------------------------------------------- 
    304 !! 1) Compute thickness and changes in each ice category 
    305 !!---------------------------------------------------------------------------------------------- 
    306        IF (lwp) THEN 
    307        WRITE(numout,*) 
    308        WRITE(numout,*) 'lim_itd_th_rem  : Remapping the ice thickness distribution' 
    309        WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    310        WRITE(numout,*) ' klbnd :       ', klbnd 
    311        WRITE(numout,*) ' kubnd :       ', kubnd 
    312        WRITE(numout,*) ' ntyp  :       ', ntyp  
    313        ENDIF 
    314  
    315        zdhice(:,:,:) = 0.0 
    316        DO jl = klbnd, kubnd 
    317           DO jj = 1, jpj 
    318              DO ji = 1, jpi 
    319                 zindb             = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)))     !0 if no ice and 1 if yes 
    320                 ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),zeps) * zindb 
    321                 zindb             = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl))) !0 if no ice and 1 if yes 
    322                 zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),zeps) * zindb 
    323                 IF (a_i(ji,jj,jl).gt.1e-6) THEN 
    324                    zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)  
    325                 ENDIF 
    326              END DO 
    327           END DO 
    328        END DO 
    329  
    330 !----------------------------------------------------------------------------------------------- 
    331 !  2) Compute fractional ice area in each grid cell 
    332 !----------------------------------------------------------------------------------------------- 
     297 
     298      !!---------------------------------------------------------------------------------------------- 
     299      !! 1) Compute thickness and changes in each ice category 
     300      !!---------------------------------------------------------------------------------------------- 
     301      IF (kt == nit000 .AND. lwp) THEN 
     302         WRITE(numout,*) 
     303         WRITE(numout,*) 'lim_itd_th_rem  : Remapping the ice thickness distribution' 
     304         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
     305         WRITE(numout,*) ' klbnd :       ', klbnd 
     306         WRITE(numout,*) ' kubnd :       ', kubnd 
     307         WRITE(numout,*) ' ntyp  :       ', ntyp  
     308      ENDIF 
     309 
     310      zdhice(:,:,:) = 0.0 
     311      DO jl = klbnd, kubnd 
     312         DO jj = 1, jpj 
     313            DO ji = 1, jpi 
     314               zindb             = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)))     !0 if no ice and 1 if yes 
     315               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),zeps) * zindb 
     316               zindb             = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl))) !0 if no ice and 1 if yes 
     317               zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),zeps) * zindb 
     318               IF (a_i(ji,jj,jl).gt.1e-6) THEN 
     319                  zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)  
     320               ENDIF 
     321            END DO 
     322         END DO 
     323      END DO 
     324 
     325      !----------------------------------------------------------------------------------------------- 
     326      !  2) Compute fractional ice area in each grid cell 
     327      !----------------------------------------------------------------------------------------------- 
    333328      at_i(:,:) = 0.0 
    334329      DO jl = klbnd, kubnd 
    … …  
    340335      END DO 
    341336 
    342 !----------------------------------------------------------------------------------------------- 
    343 !  3) Identify grid cells with ice 
    344 !----------------------------------------------------------------------------------------------- 
     337      !----------------------------------------------------------------------------------------------- 
     338      !  3) Identify grid cells with ice 
     339      !----------------------------------------------------------------------------------------------- 
    345340      nbrem = 0 
    346341      DO jj = 1, jpj 
    … …  
    357352      END DO !jj 
    358353 
    359 !----------------------------------------------------------------------------------------------- 
    360 !  4) Compute new category boundaries 
    361 !----------------------------------------------------------------------------------------------- 
     354      !----------------------------------------------------------------------------------------------- 
     355      !  4) Compute new category boundaries 
     356      !----------------------------------------------------------------------------------------------- 
    362357      !- 4.1 Compute category boundaries 
    363358      ! Tricky trick see limitd_me.F90 
    … …  
    374369            ! 
    375370            IF ( ( zht_i_o(zji,zjj,jl)  .GT.zeps ) .AND. &  
    376                  ( zht_i_o(zji,zjj,jl+1).GT.zeps ) ) THEN 
     371               ( zht_i_o(zji,zjj,jl+1).GT.zeps ) ) THEN 
    377372               !interpolate between adjacent category growth rates 
    378373               zslope = ( zdhice(zji,zjj,jl+1)     - zdhice(zji,zjj,jl) ) / & 
    379                         ( zht_i_o   (zji,zjj,jl+1) - zht_i_o   (zji,zjj,jl) ) 
     374                  ( zht_i_o   (zji,zjj,jl+1) - zht_i_o   (zji,zjj,jl) ) 
    380375               zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) + & 
    381                                     zslope * ( hi_max(jl) - zht_i_o(zji,zjj,jl) ) 
     376                  zslope * ( hi_max(jl) - zht_i_o(zji,zjj,jl) ) 
    382377            ELSEIF (zht_i_o(zji,zjj,jl).gt.zeps) THEN 
    383378               zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) 
    … …  
    391386         ! jl 
    392387 
    393       !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 
     388         !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 
    394389         DO ji = 1, nbrem 
    395390            ! jl, ji 
    … …  
    398393            ! jl, ji 
    399394            IF ( ( a_i(zji,zjj,jl) .GT.zeps) .AND. &  
    400                  ( ht_i(zji,zjj,jl).GE. zhbnew(zji,zjj,jl) ) & 
     395               ( ht_i(zji,zjj,jl).GE. zhbnew(zji,zjj,jl) ) & 
    401396               ) THEN 
    402397               zremap_flag(zji,zjj) = .false. 
    403398            ELSEIF ( ( a_i(zji,zjj,jl+1) .GT. zeps ) .AND. & 
    404                      ( ht_i(zji,zjj,jl+1).LE. zhbnew(zji,zjj,jl) ) & 
    405                    ) THEN 
     399               ( ht_i(zji,zjj,jl+1).LE. zhbnew(zji,zjj,jl) ) & 
     400               ) THEN 
    406401               zremap_flag(zji,zjj) = .false. 
    407402            ENDIF 
    408403 
    409       !- 4.3 Check that each zhbnew does not exceed maximal values hi_max   
     404            !- 4.3 Check that each zhbnew does not exceed maximal values hi_max   
    410405            ! jl, ji 
    411406            IF (zhbnew(zji,zjj,jl).gt.hi_max(jl+1)) THEN 
    … …  
    420415         ! ji 
    421416      END DO !jl 
    422              
    423 !----------------------------------------------------------------------------------------------- 
    424 !  5) Identify cells where ITD is to be remapped 
    425 !----------------------------------------------------------------------------------------------- 
    426      nbrem = 0 
    427      DO jj = 1, jpj 
    428         DO ji = 1, jpi 
    429            IF ( zremap_flag(ji,jj) ) THEN 
    430               nbrem         = nbrem + 1 
    431               nind_i(nbrem) = ji 
    432               nind_j(nbrem) = jj 
    433            ENDIF 
    434         END DO !ji 
    435      END DO !jj 
    436  
    437 !----------------------------------------------------------------------------------------------- 
    438 !  6) Fill arrays with lowermost / uppermost boundaries of 'new' categories 
    439 !----------------------------------------------------------------------------------------------- 
    440      DO jj = 1, jpj 
    441         DO ji = 1, jpi 
    442            zhb0(ji,jj) = hi_max_typ(0,ntyp) ! 0eme 
    443            zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er 
    444  
    445            zhbnew(ji,jj,klbnd-1) = 0.0 
    446             
    447            IF ( a_i(ji,jj,kubnd) .GT. zeps ) THEN 
    448               zhbnew(ji,jj,kubnd) = 3.0*ht_i(ji,jj,kubnd) - 2.0*zhbnew(ji,jj,kubnd-1) 
    449            ELSE 
    450               zhbnew(ji,jj,kubnd) = hi_max(kubnd) 
    451            ENDIF 
    452  
    453            IF ( zhbnew(ji,jj,kubnd) .LT. hi_max(kubnd-1) ) & 
    454               zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 
    455  
    456         END DO !jj 
    457      END DO !jj 
    458  
    459 !----------------------------------------------------------------------------------------------- 
    460 !  7) Compute g(h)  
    461 !----------------------------------------------------------------------------------------------- 
    462      !- 7.1 g(h) for category 1 at start of time step 
    463      CALL lim_itd_fitline(klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd), & 
    464                           g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 
    465                           hR(:,:,klbnd), zremap_flag) 
    466  
    467      !- 7.2 Area lost due to melting of thin ice (first category,  klbnd) 
    468      DO ji = 1, nbrem 
    469         zji = nind_i(ji)  
    470         zjj = nind_j(ji)  
    471        
    472         !ji 
    473         IF (a_i(zji,zjj,klbnd) .gt. zeps) THEN 
    474            zdh0 = zdhice(zji,zjj,klbnd) !decrease of ice thickness in the lower category 
    475            ! ji, a_i > zeps 
    476            IF (zdh0 .lt. 0.0) THEN !remove area from category 1 
    477               ! ji, a_i > zeps; zdh0 < 0 
    478               zdh0 = MIN(-zdh0,hi_max(klbnd)) 
    479          
    480               !Integrate g(1) from 0 to dh0 to estimate area melted 
    481               zetamax = MIN(zdh0,hR(zji,zjj,klbnd)) - hL(zji,zjj,klbnd) 
    482               IF (zetamax.gt.0.0) THEN 
    483                  zx1  = zetamax 
    484                  zx2  = 0.5 * zetamax*zetamax  
    485                  zda0 = g1(zji,zjj,klbnd) * zx2 + g0(zji,zjj,klbnd) * zx1 !ice area removed 
    486               ! Constrain new thickness <= ht_i 
    487                  zdamax = a_i(zji,zjj,klbnd) * &  
    488                           (1.0 - ht_i(zji,zjj,klbnd)/zht_i_o(zji,zjj,klbnd)) ! zdamax > 0 
    489               !ice area lost due to melting of thin ice 
    490                  zda0   = MIN(zda0, zdamax) 
    491  
    492               ! Remove area, conserving volume 
    493                  ht_i(zji,zjj,klbnd) = ht_i(zji,zjj,klbnd) &  
    494                                * a_i(zji,zjj,klbnd) / ( a_i(zji,zjj,klbnd) - zda0 ) 
    495                  a_i(zji,zjj,klbnd)  = a_i(zji,zjj,klbnd) - zda0 
    496                  v_i(zji,zjj,klbnd)  = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd) 
    497               ENDIF     ! zetamax > 0 
    498            ! ji, a_i > zeps 
    499  
    500            ELSE ! if ice accretion 
    501               ! ji, a_i > zeps; zdh0 > 0 
    502               IF ( ntyp .EQ. 1 ) zhbnew(zji,zjj,klbnd-1) = MIN(zdh0,hi_max(klbnd))  
    503               ! zhbnew was 0, and is shifted to the right to account for thin ice 
    504               ! growth in openwater (F0 = f1) 
    505               IF ( ntyp .NE. 1 ) zhbnew(zji,zjj,0) = 0  
    506               ! in other types there is 
    507               ! no open water growth (F0 = 0) 
    508            ENDIF ! zdh0  
    509  
    510            ! a_i > zeps 
    511         ENDIF ! a_i > zeps 
    512  
    513      END DO ! ji 
    514  
    515      !- 7.3 g(h) for each thickness category   
    516      DO jl = klbnd, kubnd 
    517         CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 
    518                              g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl),     & 
    519                              zremap_flag) 
    520      END DO 
    521  
    522 !----------------------------------------------------------------------------------------------- 
    523 !  8) Compute area and volume to be shifted across each boundary 
    524 !----------------------------------------------------------------------------------------------- 
    525  
    526      DO jl = klbnd, kubnd - 1 
    527         DO jj = 1, jpj 
    528            DO ji = 1, jpi 
    529               zdonor(ji,jj,jl) = 0 
    530               zdaice(ji,jj,jl) = 0.0 
    531               zdvice(ji,jj,jl) = 0.0 
    532            END DO 
    533         END DO 
    534  
    535         DO ji = 1, nbrem 
    536            zji = nind_i(ji) 
    537            zjj = nind_j(ji) 
    538             
    539            IF (zhbnew(zji,zjj,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 
    540  
    541               ! left and right integration limits in eta space 
    542               zvetamin(ji) = MAX(hi_max(jl), hL(zji,zjj,jl)) - hL(zji,zjj,jl) 
    543               zvetamax(ji) = MIN(zhbnew(zji,zjj,jl), hR(zji,zjj,jl)) - hL(zji,zjj,jl) 
    544               zdonor(zji,zjj,jl) = jl 
    545  
    546            ELSE  ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 
    547  
    548               ! left and right integration limits in eta space 
    549               zvetamin(ji) = 0.0 
    550               zvetamax(ji) = MIN(hi_max(jl), hR(zji,zjj,jl+1)) - hL(zji,zjj,jl+1) 
    551               zdonor(zji,zjj,jl) = jl + 1 
    552  
    553            ENDIF  ! zhbnew(jl) > hi_max(jl) 
    554  
    555            zetamax = MAX(zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin 
    556            zetamin = zvetamin(ji) 
    557  
    558            zx1  = zetamax - zetamin 
    559            zwk1 = zetamin*zetamin 
    560            zwk2 = zetamax*zetamax 
    561            zx2  = 0.5 * (zwk2 - zwk1) 
    562            zwk1 = zwk1 * zetamin 
    563            zwk2 = zwk2 * zetamax 
    564            zx3  = 1.0/3.0 * (zwk2 - zwk1) 
    565            nd   = zdonor(zji,zjj,jl) 
    566            zdaice(zji,zjj,jl) = g1(zji,zjj,nd)*zx2 + g0(zji,zjj,nd)*zx1 
    567            zdvice(zji,zjj,jl) = g1(zji,zjj,nd)*zx3 + g0(zji,zjj,nd)*zx2 + & 
    568                               zdaice(zji,zjj,jl)*hL(zji,zjj,nd) 
    569  
    570         END DO ! ji 
    571      END DO ! jl klbnd -> kubnd - 1 
    572  
    573 !!---------------------------------------------------------------------------------------------- 
    574 !! 9) Shift ice between categories 
    575 !!---------------------------------------------------------------------------------------------- 
    576      CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice ) 
    577  
    578 !!---------------------------------------------------------------------------------------------- 
    579 !! 10) Make sure ht_i >= minimum ice thickness hi_min 
    580 !!---------------------------------------------------------------------------------------------- 
    581  
    582     DO ji = 1, nbrem 
    583         zji = nind_i(ji) 
    584         zjj = nind_j(ji) 
    585         IF ( ( zhimin .GT. 0.0 ) .AND. &  
    586              ( ( a_i(zji,zjj,1) .GT. zeps ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) & 
    587            ) THEN 
    588            a_i(zji,zjj,1)  = a_i(zji,zjj,1) * ht_i(zji,z