New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 5053 – NEMO

Changeset 5053


Ignore:
Timestamp:
2015-02-03T18:11:02+01:00 (9 years ago)
Author:
clem
Message:

LIM3 cleaning continues. No change in the physics besides the introduction of the monocategory sea ice

Location:
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO
Files:
1 deleted
19 edited

Legend:

Unmodified
Added
Removed

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90

    r4306 r5053  
    5151   REAL(wp), PUBLIC ::   ahi0        !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
    5252   REAL(wp), PUBLIC ::   alphaevp    !: coefficient for the solution of EVP int. stresses 
    53    REAL(wp), PUBLIC ::   hminrhg = 0.001_wp    !: clem : ice volume (a*h in m) below which ice velocity is set to ocean velocity 
    5453 
    5554   REAL(wp), PUBLIC ::   usecc2                !:  = 1.0 / ( ecc * ecc ) 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90

    r4624 r5053  
    228228         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   & 
    229229         &                c_rhg, etamn, creepl, ecc, ahi0,                  & 
    230          &                nevp, telast, alphaevp, hminrhg 
     230         &                nevp, telast, alphaevp 
    231231      !!------------------------------------------------------------------- 
    232232                     
    … …  
    262262         WRITE(numout,*) '       timescale for elastic waves telast = ', telast 
    263263         WRITE(numout,*) '       coefficient for the solution of int. stresses alphaevp = ', alphaevp 
    264          WRITE(numout,*) '       min ice thickness for rheology calculations     hminrhg = ', hminrhg 
    265264      ENDIF 
    266265      ! 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r5051 r5053  
    178178   REAL(wp), PUBLIC ::   ahi0             !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
    179179   REAL(wp), PUBLIC ::   relast           !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp)  
    180    REAL(wp), PUBLIC ::   hminrhg          !: ice volume (a*h, in m) below which ice velocity is set to ocean velocity 
    181180 
    182181   !                                     !!** ice-salinity namelist (namicesal) ** 
    … …  
    300299   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice 
    301300 
    302    ! lateral melting (mono-category) 
    303    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_melt !: net melting (m) 
    304  
    305301   !!-------------------------------------------------------------------------- 
    306302   !! * Ice global state variables 
    … …  
    370366   !!-------------------------------------------------------------------------- 
    371367   !                                                  !!: ** Namelist namicerun read in sbc_lim_init ** 
     368!clem   INTEGER               , PUBLIC ::   jpl             !: number of ice  categories  
     369!clem   INTEGER               , PUBLIC ::   nlay_i          !: number of ice  layers  
     370!clem   INTEGER               , PUBLIC ::   nlay_s          !: number of snow layers  
    372371   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input) 
    373372   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output) 
    … …  
    439438         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,  qlead  (jpi,jpj) ,     & 
    440439         &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) , & 
    441          &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl) , dh_i_melt(jpi,jpj,jpl) , & 
     440         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl),    & 
    442441         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,      & 
    443442         &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,   & 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r5047 r5053  
    3131   USE timing          ! Timing 
    3232   USE limcons        ! conservation tests 
     33   USE limvar 
    3334 
    3435   IMPLICIT NONE 
    … …  
    7576      CALL wrk_alloc( jpi, jpj, zu_io, zv_io ) 
    7677      CALL wrk_alloc( jpj, zswitch, zmsk ) 
     78 
     79      CALL lim_var_agg(1)             ! aggregate ice categories 
    7780 
    7881      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only) 
    … …  
    241244      !!------------------------------------------------------------------- 
    242245      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    243       NAMELIST/namicedyn/ cw, pstar,   & 
    244          &                c_rhg, creepl, ecc, ahi0,     & 
    245          &                nevp, relast, hminrhg 
     246      NAMELIST/namicedyn/ cw, pstar, c_rhg, creepl, ecc, ahi0, nevp, relast 
    246247      !!------------------------------------------------------------------- 
    247248 
    … …  
    267268         WRITE(numout,*) '   number of iterations for subcycling              nevp   = ', nevp 
    268269         WRITE(numout,*) '   ratio of elastic timescale over ice time step    relast = ', relast 
    269          WRITE(numout,*) '   min ice thickness for rheology calculations     hminrhg = ', hminrhg 
    270270      ENDIF 
    271271      ! 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r5051 r5053  
    375375      CALL lim_var_glo2eqv 
    376376      CALL lim_var_zapsmall 
     377      CALL lim_var_agg( 1 )  
    377378 
    378379 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r5051 r5053  
    1313   !!   'key_lim3' :                                   LIM3 sea-ice model 
    1414   !!---------------------------------------------------------------------- 
    15    !!   lim_itd_th       : thermodynamics of ice thickness distribution 
    1615   !!   lim_itd_th_rem   : 
    1716   !!   lim_itd_th_reb   : 
    … …  
    2625   USE ice              ! LIM-3 variables 
    2726   USE par_ice          ! LIM-3 parameters 
    28    USE limthd_lac       ! LIM-3 lateral accretion 
    2927   USE limvar           ! LIM-3 variables 
    3028   USE limcons          ! LIM-3 conservation 
    … …  
    3432   USE wrk_nemo         ! work arrays 
    3533   USE lib_fortran      ! to use key_nosignedzero 
    36    USE timing          ! Timing 
    37    USE limcons        ! conservation tests 
     34   USE limcons          ! conservation tests 
    3835 
    3936   IMPLICIT NONE 
    4037   PRIVATE 
    4138 
    42    PUBLIC   lim_itd_th         ! called by ice_stp 
    4339   PUBLIC   lim_itd_th_rem 
    4440   PUBLIC   lim_itd_th_reb 
    … …  
    5248   !!---------------------------------------------------------------------- 
    5349CONTAINS 
    54  
    55    SUBROUTINE lim_itd_th( kt ) 
    56       !!------------------------------------------------------------------ 
    57       !!                ***  ROUTINE lim_itd_th *** 
    58       !! 
    59       !! ** Purpose :   computes the thermodynamics of ice thickness distribution 
    60       !! 
    61       !! ** Method  : 
    62       !!------------------------------------------------------------------ 
    63       INTEGER, INTENT(in) ::   kt   ! time step index 
    64       ! 
    65       INTEGER ::   ji, jj, jk, jl   ! dummy loop index          
    66       ! 
    67       REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    68       !!------------------------------------------------------------------ 
    69       IF( nn_timing == 1 )  CALL timing_start('limitd_th') 
    70  
    71       ! conservation test 
    72       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    73  
    74       IF( kt == nit000 .AND. lwp ) THEN 
    75          WRITE(numout,*) 
    76          WRITE(numout,*) 'lim_itd_th  : Thermodynamics of the ice thickness distribution' 
    77          WRITE(numout,*) '~~~~~~~~~~~' 
    78       ENDIF 
    79  
    80       !------------------------------------------------------------------------------| 
    81       !  1) Transport of ice between thickness categories.                           | 
    82       !------------------------------------------------------------------------------| 
    83       ! Given thermodynamic growth rates, transport ice between 
    84       ! thickness categories. 
    85       IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt ) 
    86       ! 
    87       CALL lim_var_glo2eqv    ! only for info 
    88       CALL lim_var_agg(1) 
    89  
    90       !------------------------------------------------------------------------------| 
    91       !  3) Add frazil ice growing in leads. 
    92       !------------------------------------------------------------------------------| 
    93       CALL lim_thd_lac 
    94       CALL lim_var_glo2eqv    ! only for info 
    95  
    96       ! MV: Could put lateral melting here, would be better I think ??? 
    97  
    98       
    99       IF(ln_ctl) THEN   ! Control print 
    100          CALL prt_ctl_info(' ') 
    101          CALL prt_ctl_info(' - Cell values : ') 
    102          CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
    103          CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th  : cell area :') 
    104          CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :') 
    105          CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :') 
    106          CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th  : vt_s      :') 
    107          DO jl = 1, jpl 
    108             CALL prt_ctl_info(' ') 
    109             CALL prt_ctl_info(' - Category : ', ivar1=jl) 
    110             CALL prt_ctl_info('   ~~~~~~~~~~') 
    111             CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : a_i      : ') 
    112             CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_i     : ') 
    113             CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_s     : ') 
    114             CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_i      : ') 
    115             CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_s      : ') 
    116             CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : e_s      : ') 
    117             CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_th  : t_su     : ') 
    118             CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : t_snow   : ') 
    119             CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ') 
    120             CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ') 
    121             DO jk = 1, nlay_i 
    122                CALL prt_ctl_info(' ') 
    123                CALL prt_ctl_info(' - Layer : ', ivar1=jk) 
    124                CALL prt_ctl_info('   ~~~~~~~') 
    125                CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ') 
    126                CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ') 
    127             END DO 
    128          END DO 
    129       ENDIF 
    130       ! 
    131       ! conservation test 
    132       IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    133       ! 
    134      IF( nn_timing == 1 )  CALL timing_stop('limitd_th') 
    135    END SUBROUTINE lim_itd_th 
    136    ! 
    13750 
    13851   SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 
    … …  
    15669      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars 
    15770      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      - 
    158       REAL(wp) ::   zx3,             zareamin          !   -      - 
     71      REAL(wp) ::   zx3         
    15972      CHARACTER (len = 15) :: fieldid 
    16073 
    … …  
    191104      CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
    192105 
    193       zareamin = epsi10   !minimum area in thickness categories tolerated by the conceptors of the model 
    194  
    195106      !!---------------------------------------------------------------------------------------------- 
    196107      !! 0) Conservation checkand changes in each ice category 
    … …  
    242153      DO jj = 1, jpj 
    243154         DO ji = 1, jpi 
    244             IF ( at_i(ji,jj) .gt. zareamin ) THEN 
     155            IF ( at_i(ji,jj) > epsi10 ) THEN 
    245156               nbrem         = nbrem + 1 
    246157               nind_i(nbrem) = ji 
    … …  
    1014925   !!---------------------------------------------------------------------- 
    1015926CONTAINS 
    1016    SUBROUTINE lim_itd_th           ! Empty routines 
    1017    END SUBROUTINE lim_itd_th 
    1018    SUBROUTINE lim_itd_th_ini 
    1019    END SUBROUTINE lim_itd_th_ini 
    1020927   SUBROUTINE lim_itd_th_rem 
    1021928   END SUBROUTINE lim_itd_th_rem 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r5051 r5053  
    152152 
    153153      REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter 
     154      REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 
    154155      !!------------------------------------------------------------------- 
    155156 
    … …  
    161162#if  defined key_lim2 && ! defined key_lim2_vp 
    162163# if defined key_agrif 
    163      USE ice_2, vt_s => hsnm 
    164      USE ice_2, vt_i => hicm 
     164      USE ice_2, vt_s => hsnm 
     165      USE ice_2, vt_i => hicm 
    165166# else 
    166      vt_s => hsnm 
    167      vt_i => hicm 
     167      vt_s => hsnm 
     168      vt_i => hicm 
    168169# endif 
    169      at_i(:,:) = 1. - frld(:,:) 
     170      at_i(:,:) = 1. - frld(:,:) 
    170171#endif 
    171172#if defined key_agrif && defined key_lim2  
    172     CALL agrif_rhg_lim2_load      ! First interpolation of coarse values 
     173      CALL agrif_rhg_lim2_load      ! First interpolation of coarse values 
    173174#endif 
    174175      ! 
    … …  
    189190#endif 
    190191 
    191 !CDIR NOVERRCHK 
    192192      DO jj = k_j1 , k_jpj       ! Ice mass and temp variables 
    193 !CDIR NOVERRCHK 
    194193         DO ji = 1 , jpi 
    195194#if defined key_lim3 
    … …  
    206205      ! Ice strength on grid cell corners (zpreshc) 
    207206      ! needed for calculation of shear stress  
    208 !CDIR NOVERRCHK 
    209207      DO jj = k_j1+1, k_jpj-1 
    210 !CDIR NOVERRCHK 
    211208         DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 
    212             zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 
    213                &              tms(ji,jj+1)   * wght(ji+1,jj+1,1,2) + & 
    214                &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + & 
    215                &              tms(ji,jj)     * wght(ji+1,jj+1,1,1) 
    216             zpreshc(ji,jj) = (  zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + & 
    217                &                zpresh(ji,jj+1)   * wght(ji+1,jj+1,1,2) + & 
    218                &                zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + &  
    219                &                zpresh(ji,jj)     * wght(ji+1,jj+1,1,1)   & 
     209            zstms          =  tms(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + tms(ji,jj+1) * wght(ji+1,jj+1,1,2) +   & 
     210               &              tms(ji+1,jj)   * wght(ji+1,jj+1,2,1) + tms(ji,jj)   * wght(ji+1,jj+1,1,1) 
     211            zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   & 
     212               &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     & 
    220213               &             ) / MAX( zstms, zepsi ) 
    221214         END DO 
    222215      END DO 
    223  
    224216      CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 
    225217      ! 
    … …  
    242234 
    243235      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
    244           !                                             
    245           ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
    246           !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 
     236         !                                             
     237         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1} 
     238         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1} 
    247239         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp      
    248           ! 
    249           ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
    250           !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
     240         ! 
     241         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1} 
     242         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1}) 
    251243         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    252           ! 
     244         ! 
    253245         zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:)  ) * r1_rau0 
    254           ! 
     246         ! 
    255247      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
    256248         zpice(:,:) = ssh_m(:,:) 
    … …  
    363355               ! 
    364356               ! 
    365                divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
    366                   &             -e2u(ji-1,jj)*u_ice(ji-1,jj)                  & 
    367                   &             +e1v(ji,jj)*v_ice(ji,jj)                      & 
    368                   &             -e1v(ji,jj-1)*v_ice(ji,jj-1)                  & 
    369                   &             )                                             & 
    370                   &            / area(ji,jj) 
    371  
    372                zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    & 
    373                   &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                & 
    374                   &           )*e2t(ji,jj)*e2t(ji,jj)                      & 
    375                   &          -( v_ice(ji,jj)/e1v(ji,jj)                    & 
    376                   &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                & 
    377                   &           )*e1t(ji,jj)*e1t(ji,jj)                      & 
    378                   &         )                                              & 
    379                   &        / area(ji,jj) 
     357               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     358                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     359                  &            ) / area(ji,jj) 
     360 
     361               zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     362                  &         - ( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     363                  &         ) / area(ji,jj) 
    380364 
    381365               ! 
    382                zds(ji,jj) = ( ( u_ice(ji,jj+1)/e1u(ji,jj+1)                & 
    383                   &            -u_ice(ji,jj)/e1u(ji,jj)                    & 
    384                   &           )*e1f(ji,jj)*e1f(ji,jj)                      & 
    385                   &          +( v_ice(ji+1,jj)/e2v(ji+1,jj)                & 
    386                   &            -v_ice(ji,jj)/e2v(ji,jj)                    & 
    387                   &           )*e2f(ji,jj)*e2f(ji,jj)                      & 
    388                   &         )                                              & 
    389                   &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) & 
    390                   &        * tmi(ji,jj) * tmi(ji,jj+1)                     & 
    391                   &        * tmi(ji+1,jj) * tmi(ji+1,jj+1) 
    392  
    393  
    394                v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   & 
    395                   &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 
    396                   &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)  
    397  
    398                u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   & 
    399                   &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 
    400                   &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
     366               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     367                  &         + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     368                  &         ) / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2._wp - tmf(ji,jj) )   & 
     369                  &         * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1) 
     370 
     371 
     372               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji,jj)   + v_ice(ji,jj-1)   ) * e1t(ji+1,jj)   & 
     373                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji,jj) )   & 
     374                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)  
     375 
     376               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj)   + u_ice(ji-1,jj)   ) * e2t(ji,jj+1)   & 
     377                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj) )   & 
     378                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 
    401379 
    402380            END DO 
    … …  
    404382         CALL lbc_lnk( v_ice1, 'U', -1. )   ;   CALL lbc_lnk( u_ice2, 'V', -1. )      ! lateral boundary cond. 
    405383 
    406 !CDIR NOVERRCHK 
    407384         DO jj = k_j1+1, k_jpj-1 
    408 !CDIR NOVERRCHK 
    409385            DO ji = fs_2, fs_jpim1 
    410386 
    411387               !- Calculate Delta at centre of grid cells 
    412                zdst      = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          & 
    413                   &          - e2u(ji-1, jj) * v_ice1(ji-1,jj)          & 
    414                   &          + e1v(ji, jj  ) * u_ice2(ji,jj  )          & 
    415                   &          - e1v(ji, jj-1) * u_ice2(ji,jj-1)          & 
    416                   &          )                                          & 
    417                   &         / area(ji,jj) 
    418  
    419                delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )   
     388               zdst      = ( e2u(ji, jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   & 
     389                  &        + e1v(ji, jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   & 
     390                  &        ) / area(ji,jj) 
     391 
     392               delta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) * usecc2 )   
    420393               delta_i(ji,jj) = delta + creepl 
    421394               !-Calculate stress tensor components zs1 and zs2  
    … …  
    432405         CALL lbc_lnk( zs2(:,:), 'T', 1. ) 
    433406 
    434 !CDIR NOVERRCHK 
    435407         DO jj = k_j1+1, k_jpj-1 
    436 !CDIR NOVERRCHK 
    437408            DO ji = fs_2, fs_jpim1 
    438409               !- Calculate Delta on corners 
    … …  
    490461         IF (MOD(jter,2).eq.0) THEN  
    491462 
    492 !CDIR NOVERRCHK 
    493463            DO jj = k_j1+1, k_jpj-1 
    494 !CDIR NOVERRCHK 
    495464               DO ji = fs_2, fs_jpim1 
    496465                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 
    … …  
    520489#endif          
    521490 
    522 !CDIR NOVERRCHK 
    523491            DO jj = k_j1+1, k_jpj-1 
    524 !CDIR NOVERRCHK 
    525492               DO ji = fs_2, fs_jpim1 
    526493 
    … …  
    551518 
    552519         ELSE  
    553 !CDIR NOVERRCHK 
    554520            DO jj = k_j1+1, k_jpj-1 
    555 !CDIR NOVERRCHK 
    556521               DO ji = fs_2, fs_jpim1 
    557522                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj) 
    … …  
    581546#endif          
    582547 
    583 !CDIR NOVERRCHK 
    584548            DO jj = k_j1+1, k_jpj-1 
    585 !CDIR NOVERRCHK 
    586549               DO ji = fs_2, fs_jpim1 
    587550                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj) 
    … …  
    628591      ! 4) Prevent ice velocities when the ice is thin 
    629592      !------------------------------------------------------------------------------! 
    630       ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 
    631       ! ocean velocity,  
    632       ! This prevents high velocity when ice is thin 
    633 !CDIR NOVERRCHK 
     593      ! If the ice volume is below zvmin then ice velocity should equal the 
     594      ! ocean velocity. This prevents high velocity when ice is thin 
    634595      DO jj = k_j1+1, k_jpj-1 
    635 !CDIR NOVERRCHK 
    636596         DO ji = fs_2, fs_jpim1 
    637             IF ( vt_i(ji,jj) .LE. hminrhg ) THEN 
     597            IF ( vt_i(ji,jj) <= zvmin ) THEN 
    638598               u_ice(ji,jj) = u_oce(ji,jj) 
    639599               v_ice(ji,jj) = v_oce(ji,jj) 
    … …  
    655615      DO jj = k_j1+1, k_jpj-1  
    656616         DO ji = fs_2, fs_jpim1 
    657             IF ( vt_i(ji,jj) .LE. hminrhg ) THEN 
    658                v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   & 
    659                   &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 
    660                   &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj) 
    661  
    662                u_ice2(ji,jj)  = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj+1)   & 
    663                   &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj)) & 
    664                   &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj) 
     617            IF ( vt_i(ji,jj) <= zvmin ) THEN 
     618               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     & 
     619                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )  & 
     620                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 
     621 
     622               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
     623                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )  & 
     624                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj) 
    665625            ENDIF  
    666626         END DO 
    … …  
    671631 
    672632      ! Recompute delta, shear and div, inputs for mechanical redistribution  
    673 !CDIR NOVERRCHK 
    674633      DO jj = k_j1+1, k_jpj-1 
    675 !CDIR NOVERRCHK 
    676634         DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi 
    677635            !- divu_i(:,:), zdt(:,:): divergence and tension at centre  
    678636            !- zds(:,:): shear on northeast corner of grid cells 
    679             IF ( vt_i(ji,jj) .LE. hminrhg ) THEN 
     637            IF ( vt_i(ji,jj) <= zvmin ) THEN 
    680638 
    681639               divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
    … …  
    722680      ! 5) Store stress tensor and its invariants 
    723681      !------------------------------------------------------------------------------! 
    724       ! 
    725682      ! * Invariants of the stress tensor are required for limitd_me 
    726683      !   (accelerates convergence and improves stability) 
    727684      DO jj = k_j1+1, k_jpj-1 
    728685         DO ji = fs_2, fs_jpim1 
    729             ! begin TECLIM change  
    730686            zdst= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &    
    731687               &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &    
    … …  
    733689               &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj)  
    734690            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 
    735             ! end TECLIM change 
    736691         END DO 
    737692      END DO 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90

    r5051 r5053  
    2727   USE wrk_nemo       ! work arrays 
    2828   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     29   USE limctl 
    2930 
    3031   IMPLICIT NONE 
    … …  
    8687      ENDIF 
    8788      ! 
     89      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
    8890   END SUBROUTINE lim_rst_opn 
    8991 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r5051 r5053  
    4343   USE iom 
    4444   USE domvvl           ! Variable volume 
     45   USE limctl 
    4546 
    4647   IMPLICIT NONE 
    … …  
    225226      ENDIF 
    226227 
     228      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
    227229 
    228230      IF(ln_ctl) THEN 
    … …  
    270272      ! 
    271273      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step) 
    272 !CDIR NOVERRCHK 
    273274         DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point) 
    274 !CDIR NOVERRCHK 
    275275            DO ji = fs_2, fs_jpim1 
    276276               !                                               ! 2*(U_ice-U_oce) at T-point 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5051 r5053  
    3434   USE limthd_sal     ! LIM: thermodynamics, ice salinity 
    3535   USE limthd_ent     ! LIM: thermodynamics, ice enthalpy redistribution 
     36   USE limthd_lac     ! LIM-3 lateral accretion 
     37   USE limitd_th      ! remapping thickness distribution 
    3638   USE limtab         ! LIM: 1D <==> 2D transformation 
    3739   USE limvar         ! LIM: sea-ice variables 
    … …  
    4446   USE timing         ! Timing 
    4547   USE limcons        ! conservation tests 
     48   USE limctl 
    4649 
    4750   IMPLICIT NONE 
    … …  
    8083      !! ** References :  
    8184      !!--------------------------------------------------------------------- 
    82       INTEGER, INTENT(in) ::   kt    ! number of iteration 
     85      INTEGER, INTENT(in) :: kt    ! number of iteration 
    8386      !! 
    8487      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices 
    8588      INTEGER  :: nbpb             ! nb of icy pts for vertical thermo calculations 
    86       INTEGER  :: nbplm            ! nb of icy pts for lateral melting calculations (mono-cat) 
    8789      INTEGER  :: ii, ij           ! temporary dummy loop index 
    88       REAL(wp) :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04) 
    89       REAL(wp) :: zch        = 0.0057_wp    ! heat transfer coefficient 
    90       REAL(wp) :: zareamin  
    9190      REAL(wp) :: zfric_u, zqld, zqfr 
    92       ! 
    9391      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
     92      REAL(wp), PARAMETER :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04) 
     93      REAL(wp), PARAMETER :: zch        = 0.0057_wp    ! heat transfer coefficient 
    9494      ! 
    9595      REAL(wp), POINTER, DIMENSION(:,:) ::  zqsr, zqns 
    … …  
    106106      !------------------------------------------------------------------------! 
    107107      ftr_ice(:,:,:) = 0._wp  ! part of solar radiation transmitted through the ice 
    108  
    109108 
    110109      !-------------------- 
    … …  
    162161      ENDIF 
    163162 
    164 !CDIR NOVERRCHK 
    165163      DO jj = 1, jpj 
    166 !CDIR NOVERRCHK 
    167164         DO ji = 1, jpi 
    168             rswitch          = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 
     165            rswitch  = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 
    169166            ! 
    170167            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
    … …  
    260257         ENDIF 
    261258 
    262          zareamin = epsi10 
    263259         nbpb = 0 
    264260         DO jj = 1, jpj 
    265261            DO ji = 1, jpi 
    266                IF ( a_i(ji,jj,jl) .gt. zareamin ) THEN      
     262               IF ( a_i(ji,jj,jl) > epsi10 ) THEN      
    267263                  nbpb      = nbpb  + 1 
    268264                  npb(nbpb) = (jj - 1) * jpi + ji 
    … …  
    442438              CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj ) 
    443439              CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb), jpi, jpj ) 
    444  
    445 !clem            IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN 
    446 !clem              CALL tab_1d_2d( nbpb, dh_i_melt(:,:,jl) , npb, dh_i_melt_1d(1:nbpb) , jpi, jpj ) 
    447 !clem            ENDIF 
    448440            ! 
    449441              CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
    … …  
    460452 
    461453      !------------------------ 
    462       ! 5.1) Ice heat content               
     454      ! Ice heat content               
    463455      !------------------------ 
    464456      ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 
    … …  
    470462 
    471463      !------------------------ 
    472       ! 5.2) Snow heat content               
     464      ! Snow heat content               
    473465      !------------------------ 
    474466      ! Enthalpies are global variables we have to readjust the units (heat content in Joules) 
    … …  
    478470         END DO 
    479471      END DO 
     472  
     473      !------------------------ 
     474      ! Ice natural aging               
     475      !------------------------ 
     476      oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rdt_ice /rday 
    480477 
    481478      !---------------------------------- 
    482       ! 5.3) Change thickness to volume 
     479      ! Change thickness to volume 
    483480      !---------------------------------- 
    484481      CALL lim_var_eqv2glo 
    … …  
    487484      ! 5.4) Diagnostic thermodynamic growth rates 
    488485      !-------------------------------------------- 
     486      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
     487 
    489488      IF(ln_ctl) THEN            ! Control print 
    490489         CALL prt_ctl_info(' ') 
    … …  
    522521      CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 
    523522 
     523      !------------------------------------------------------------------------------| 
     524      !  1) Transport of ice between thickness categories.                           | 
     525      !------------------------------------------------------------------------------| 
     526      ! Given thermodynamic growth rates, transport ice between 
     527      ! thickness categories. 
     528      IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt ) 
     529      ! 
     530      CALL lim_var_glo2eqv    ! only for info 
     531      CALL lim_var_agg(1) 
     532 
     533      !------------------------------------------------------------------------------| 
     534      !  3) Add frazil ice growing in leads. 
     535      !------------------------------------------------------------------------------| 
     536      CALL lim_thd_lac 
     537      CALL lim_var_glo2eqv    ! only for info 
     538       
     539      IF(ln_ctl) THEN   ! Control print 
     540         CALL prt_ctl_info(' ') 
     541         CALL prt_ctl_info(' - Cell values : ') 
     542         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ') 
     543         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th  : cell area :') 
     544         CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :') 
     545         CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :') 
     546         CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th  : vt_s      :') 
     547         DO jl = 1, jpl 
     548            CALL prt_ctl_info(' ') 
     549            CALL prt_ctl_info(' - Category : ', ivar1=jl) 
     550            CALL prt_ctl_info('   ~~~~~~~~~~') 
     551            CALL prt_ctl(tab2d_1=a_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : a_i      : ') 
     552            CALL prt_ctl(tab2d_1=ht_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_i     : ') 
     553            CALL prt_ctl(tab2d_1=ht_s  (:,:,jl)   , clinfo1= ' lim_itd_th  : ht_s     : ') 
     554            CALL prt_ctl(tab2d_1=v_i   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_i      : ') 
     555            CALL prt_ctl(tab2d_1=v_s   (:,:,jl)   , clinfo1= ' lim_itd_th  : v_s      : ') 
     556            CALL prt_ctl(tab2d_1=e_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : e_s      : ') 
     557            CALL prt_ctl(tab2d_1=t_su  (:,:,jl)   , clinfo1= ' lim_itd_th  : t_su     : ') 
     558            CALL prt_ctl(tab2d_1=t_s   (:,:,1,jl) , clinfo1= ' lim_itd_th  : t_snow   : ') 
     559            CALL prt_ctl(tab2d_1=sm_i  (:,:,jl)   , clinfo1= ' lim_itd_th  : sm_i     : ') 
     560            CALL prt_ctl(tab2d_1=smv_i (:,:,jl)   , clinfo1= ' lim_itd_th  : smv_i    : ') 
     561            DO jk = 1, nlay_i 
     562               CALL prt_ctl_info(' ') 
     563               CALL prt_ctl_info(' - Layer : ', ivar1=jk) 
     564               CALL prt_ctl_info('   ~~~~~~~') 
     565               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : t_i      : ') 
     566               CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th  : e_i      : ') 
     567            END DO 
     568         END DO 
     569      ENDIF 
    524570      ! 
    525571      ! conservation test 
    526572      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    527       ! 
     573 
    528574      IF( nn_timing == 1 )  CALL timing_stop('limthd') 
    529575 
    … …  
    571617      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop 
    572618      INTEGER             ::   ji            ! dummy loop indices 
    573  
    574       WRITE(numout,*) ' Lateral melting ON ' 
     619      REAL(wp)            ::   zhi_bef       ! ice thickness before thermo 
     620      REAL(wp)            ::   zdh_mel       ! net melting 
     621 
    575622      DO ji = kideb, kiut 
    576          IF( ht_i_1d(ji) > epsi10 .AND. dh_i_melt_1d(ji) < 0._wp ) THEN      
    577             a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * dh_i_melt_1d(ji) / ( 2._wp * ht_i_1d(ji) ) )  
    578          END IF 
     623         zhi_bef = ht_i_1d(ji) - ( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
     624         zdh_mel = MIN( 0._wp,     dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 
     625         IF( zdh_mel < 0._wp )   & 
     626            &   a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) )  
    579627      END DO 
    580628      at_i_1d(:) = a_i_1d(:) 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5051 r5053  
    689689      END DO 
    690690       
    691       ! 
    692       !------------------------------------------------------------- 
    693       ! Net melting (input for thin ice melting, mono-category case) 
    694       !------------------------------------------------------------- 
    695       ! 
    696       IF ( ( ( nn_monocat == 1 ) .OR. ( nn_monocat == 4 ) ) .AND. ( jpl == 1 ) ) THEN 
    697          DO ji = kideb, kiut 
    698             dh_i_melt_1d(ji) = MIN( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji), 0._wp ) 
    699          END DO 
    700       ENDIF 
    701  
    702691      CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
    703692      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r5051 r5053  
    2222   USE limadv         ! ice advection 
    2323   USE limhdf         ! ice horizontal diffusion 
    24    USE limvar         ! clem for ice thickness correction 
     24   USE limvar         !  
    2525   ! 
    2626   USE in_out_manager ! I/O manager 
    … …  
    3232   USE timing         ! Timing 
    3333   USE limcons        ! conservation tests 
     34   USE limctl         ! control prints 
    3435 
    3536   IMPLICIT NONE 
    … …  
    455456      ENDIF 
    456457 
     458      CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
     459 
     460      ! ------------------------------------------------- 
     461      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
     462      ! ------------------------------------------------- 
    457463      IF(ln_ctl) THEN   ! Control print 
    458464         CALL prt_ctl_info(' ') 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r5051 r5053  
    5959CONTAINS 
    6060 
    61    SUBROUTINE lim_update1 
     61   SUBROUTINE lim_update1( kt ) 
    6262      !!------------------------------------------------------------------- 
    6363      !!               ***  ROUTINE lim_update1  *** 
    … …  
    6767      !!                 
    6868      !!--------------------------------------------------------------------- 
     69      INTEGER, INTENT(in) :: kt    ! number of iteration 
    6970      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    7071      INTEGER  ::   i_ice_switch 
    … …  
    7677 
    7778      IF( ln_limdyn ) THEN  
     79 
     80      IF( kt == nit000 .AND. lwp ) THEN 
     81         WRITE(numout,*) ' lim_update1 '  
     82         WRITE(numout,*) ' ~~~~~~~~~~~ ' 
     83      ENDIF 
    7884 
    7985      ! conservation test 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r5051 r5053  
    4141   USE timing          ! Timing 
    4242   USE limcons        ! conservation tests 
     43   USE limctl 
    4344 
    4445   IMPLICIT NONE 
    … …  
    5657CONTAINS 
    5758 
    58    SUBROUTINE lim_update2 
     59   SUBROUTINE lim_update2( kt ) 
    5960      !!------------------------------------------------------------------- 
    6061      !!               ***  ROUTINE lim_update2  *** 
    … …  
    6465      !! 
    6566      !!--------------------------------------------------------------------- 
     67      INTEGER, INTENT(in) :: kt    ! number of iteration 
    6668      INTEGER  ::   ji, jj, jk, jl    ! dummy loop indices 
    6769      INTEGER  ::   i_ice_switch 
    … …  
    7173      !!------------------------------------------------------------------- 
    7274      IF( nn_timing == 1 )  CALL timing_start('limupdate2') 
     75 
     76      IF( kt == nit000 .AND. lwp ) THEN 
     77         WRITE(numout,*) ' lim_update2 ' 
     78         WRITE(numout,*) ' ~~~~~~~~~~~ ' 
     79      ENDIF 
    7380 
    7481      ! conservation test 
    … …  
    208215      ! conservation test 
    209216      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     217 
     218      IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
    210219 
    211220      IF(ln_ctl) THEN   ! Control print 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r5051 r5053  
    3030   !!====================================================================== 
    3131   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code 
    32    !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
     32   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
    3333   !!---------------------------------------------------------------------- 
    3434#if defined key_lim3 
    3535   !!---------------------------------------------------------------------- 
    3636   !!   'key_lim3'                                      LIM3 sea-ice model 
    37    !!---------------------------------------------------------------------- 
    38    !!   lim_var_agg       :  
    39    !!   lim_var_glo2eqv   : 
    40    !!   lim_var_eqv2glo   : 
    41    !!   lim_var_salprof   :  
    42    !!   lim_var_salprof1d : 
    43    !!   lim_var_bv        : 
    4437   !!---------------------------------------------------------------------- 
    4538   USE par_oce        ! ocean parameters 
    … …  
    5851   PRIVATE 
    5952 
    60    PUBLIC   lim_var_agg          ! 
    61    PUBLIC   lim_var_glo2eqv      ! 
    62    PUBLIC   lim_var_eqv2glo      ! 
    63    PUBLIC   lim_var_salprof      ! 
    64    PUBLIC   lim_var_icetm        ! 
    65    PUBLIC   lim_var_bv           ! 
    66    PUBLIC   lim_var_salprof1d    ! 
     53   PUBLIC   lim_var_agg           
     54   PUBLIC   lim_var_glo2eqv       
     55   PUBLIC   lim_var_eqv2glo       
     56   PUBLIC   lim_var_salprof       
     57   PUBLIC   lim_var_icetm         
     58   PUBLIC   lim_var_bv            
     59   PUBLIC   lim_var_salprof1d     
    6760   PUBLIC   lim_var_zapsmall 
     61   PUBLIC   lim_var_itd 
    6862 
    6963   !!---------------------------------------------------------------------- 
    70    !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     64   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011) 
    7165   !! $Id$ 
    7266   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    … …  
    130124            DO jj = 1, jpj 
    131125               DO ji = 1, jpi 
    132                   et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                       ! snow heat content 
     126                  et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                           ! snow heat content 
    133127                  rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )  
    134128                  smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch   ! ice salinity 
    … …  
    192186      ! Ice temperatures 
    193187      !------------------- 
    194 !CDIR NOVERRCHK 
    195       DO jl = 1, jpl 
    196 !CDIR NOVERRCHK 
     188      DO jl = 1, jpl 
    197189         DO jk = 1, nlay_i 
    198 !CDIR NOVERRCHK 
    199             DO jj = 1, jpj 
    200 !CDIR NOVERRCHK 
     190            DO jj = 1, jpj 
    201191               DO ji = 1, jpi 
    202192                  !                                                              ! Energy of melting q(S,T) [J.m-3] 
    203193                  rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes 
    204194                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp)  
    205                   zq_i    = zq_i * unit_fac                             !convert units 
    206                   ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt                       ! Ice layer melt temperature 
     195                  zq_i    = zq_i * unit_fac                             ! convert units 
     196                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt              ! Ice layer melt temperature 
    207197                  ! 
    208198                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation) 
    … …  
    289279      !!              alpha-weighted linear combination of s_inf and s_zero 
    290280      !! 
    291       !! ** References : Vancoppenolle et al., 2007 (in preparation) 
     281      !! ** References : Vancoppenolle et al., 2007 
    292282      !!------------------------------------------------------------------ 
    293283      INTEGER  ::   ji, jj, jk, jl   ! dummy loop index 
    … …  
    368358         ! 
    369359         DO jl = 1, jpl 
    370 !CDIR NOVERRCHK 
    371360            DO jk = 1, nlay_i 
    372361               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp) 
    … …  
    451440      ! 
    452441      INTEGER  ::   ji, jk    ! dummy loop indices 
    453       INTEGER  ::   ii, ij  ! local integers 
     442      INTEGER  ::   ii, ij    ! local integers 
    454443      REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars 
    455444      REAL(wp) ::   zalpha, zswi0, zswi01, zswibal, zs_zero              !   -      - 
    … …  
    483472         dummy_fac2 = 1._wp / REAL(nlay_i,wp) 
    484473 
    485 !CDIR NOVERRCHK 
    486474         DO jk = 1, nlay_i 
    487 !CDIR NOVERRCHK 
    488475            DO ji = kideb, kiut 
    489476               ii =  MOD( npb(ji) - 1 , jpi ) + 1 
    … …  
    502489               ! weighting the profile 
    503490               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 
    504             END DO ! ji 
    505          END DO ! jk 
    506  
    507       ENDIF ! num_sal 
     491            END DO  
     492         END DO  
     493 
     494      ENDIF  
    508495 
    509496      !------------------------------------------------------- 
    … …  
    515502         sm_i_1d(:) = 2.30_wp 
    516503         ! 
    517 !CDIR NOVERRCHK 
    518504         DO jk = 1, nlay_i 
    519505            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp) 
    … …  
    603589      ! 
    604590   END SUBROUTINE lim_var_zapsmall 
     591 
     592   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i ) 
     593      !!------------------------------------------------------------------ 
     594      !!                ***  ROUTINE lim_var_itd   *** 
     595      !! 
     596      !! ** Purpose :  converting 1-cat ice to multiple ice categories 
     597      !! 
     598      !!                  ice thickness distribution follows a gaussian law 
     599      !!               around the concentration of the most likely ice thickness 
     600      !!                           (similar as limistate.F90) 
     601      !! 
     602      !! ** Method:   Iterative procedure 
     603      !!                 
     604      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 
     605      !! 
     606      !!               2) Check whether the distribution conserves area and volume, positivity and 
     607      !!                  category boundaries 
     608      !!               
     609      !!               3) If not (input ice is too thin), the last category is empty and 
     610      !!                  the number of categories is reduced (jpl-1) 
     611      !! 
     612      !!               4) Iterate until ok (SUM(itest(:) = 4) 
     613      !! 
     614      !! ** Arguments : zhti: 1-cat ice thickness 
     615      !!                zhts: 1-cat snow depth 
     616      !!                zai : 1-cat ice concentration 
     617      !! 
     618      !! ** Output    : jpl-cat  
     619      !! 
     620      !!  (Example of application: BDY forcings when input are cell averaged)   
     621      !! 
     622      !!------------------------------------------------------------------- 
     623      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code 
     624      !!                    2014    (C. Rousset)        Rewriting 
     625      !!------------------------------------------------------------------- 
     626      !! Local variables 
     627      INTEGER  :: ji, jk, jl             ! dummy loop indices 
     628      INTEGER  :: ijpij, i_fill, jl0   
     629      REAL(wp) :: zarg, zV, zconv, zdh 
     630      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables 
     631      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables 
     632      INTEGER , POINTER, DIMENSION(:)         ::   itest 
     633  
     634      CALL wrk_alloc( 4, itest ) 
     635      !-------------------------------------------------------------------- 
     636      ! initialisation of variables 
     637      !-------------------------------------------------------------------- 
     638      ijpij = SIZE(zhti,1) 
     639      zht_i(1:ijpij,1:jpl) = 0._wp 
     640      zht_s(1:ijpij,1:jpl) = 0._wp 
     641      za_i (1:ijpij,1:jpl) = 0._wp 
     642 
     643      ! ---------------------------------------- 
     644      ! distribution over the jpl ice categories 
     645      ! ---------------------------------------- 
     646      DO ji = 1, ijpij 
     647          
     648         IF( zhti(ji) > 0._wp ) THEN 
     649 
     650         ! initialisation of tests 
     651         itest(:)  = 0 
     652          
     653         i_fill = jpl + 1                                             !==================================== 
     654         DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories   
     655            ! iteration                                               !==================================== 
     656            i_fill = i_fill - 1 
     657             
     658            ! initialisation of ice variables for each try 
     659            zht_i(ji,1:jpl) = 0._wp 
     660            za_i (ji,1:jpl) = 0._wp 
     661             
     662            ! *** case very thin ice: fill only category 1 
     663            IF ( i_fill == 1 ) THEN 
     664               zht_i(ji,1) = zhti(ji) 
     665               za_i (ji,1) = zai (ji) 
     666 
     667            ! *** case ice is thicker: fill categories >1 
     668            ELSE 
     669 
     670               ! Fill ice thicknesses except the last one (i_fill) by hmean  
     671               DO jl = 1, i_fill - 1 
     672                  zht_i(ji,jl) = hi_mean(jl) 
     673               END DO 
     674                
     675               ! find which category (jl0) the input ice thickness falls into 
     676               jl0 = i_fill 
     677               DO jl = 1, i_fill 
     678                  IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 
     679                     jl0 = jl 
     680           CYCLE 
     681                  ENDIF 
     682               END DO 
     683                
     684               ! Concentrations in the (i_fill-1) categories  
     685               za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 
     686               DO jl = 1, i_fill - 1 
     687                  IF ( jl == jl0 ) CYCLE 
     688                  zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 
     689                  za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2) 
     690               END DO 
     691                
     692               ! Concentration in the last (i_fill) category 
     693               za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 
     694                
     695               ! Ice thickness in the last (i_fill) category 
     696               zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 
     697               zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill)  
     698                
     699            ENDIF ! case ice is thick or thin 
     700             
     701            !--------------------- 
     702            ! Compatibility tests 
     703            !---------------------  
     704            ! Test 1: area conservation 
     705            zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 
     706            IF ( zconv < epsi06 ) itest(1) = 1 
     707             
     708            ! Test 2: volume conservation 
     709            zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 
     710            IF ( zconv < epsi06 ) itest(2) = 1 
     711             
     712            ! Test 3: thickness of the last category is in-bounds ? 
     713            IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 
     714             
     715            ! Test 4: positivity of ice concentrations 
     716            itest(4) = 1 
     717            DO jl = 1, i_fill 
     718               IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 
     719            END DO             
     720                                                           !============================ 
     721         END DO                                            ! end iteration on categories 
     722                                                           !============================ 
     723         ENDIF ! if zhti > 0 
     724      END DO ! i loop 
     725 
     726      ! ------------------------------------------------ 
     727      ! Adding Snow in each category where za_i is not 0 
     728      ! ------------------------------------------------  
     729      DO jl = 1, jpl 
     730         DO ji = 1, ijpij 
     731            IF( za_i(ji,jl) > 0._wp ) THEN 
     732               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 
     733               ! In case snow load is in excess that would lead to transformation from snow to ice 
     734               ! Then, transfer the snow excess into the ice (different from limthd_dh) 
     735               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 )  
     736               ! recompute ht_i, ht_s avoiding out of bounds values 
     737               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh ) 
     738               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn ) 
     739            ENDIF 
     740         ENDDO 
     741      ENDDO 
     742 
     743      CALL wrk_dealloc( 4, itest ) 
     744      ! 
     745    END SUBROUTINE lim_var_itd 
     746 
    605747 
    606748#else 
    … …  
    623765   SUBROUTINE lim_var_zapsmall 
    624766   END SUBROUTINE lim_var_zapsmall 
     767   SUBROUTINE lim_var_itd 
     768   END SUBROUTINE lim_var_itd 
    625769#endif 
    626770 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r5048 r5053  
    110110   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott     !: Ice bottom accretion/ablation  [m] 
    111111   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice    !: Snow ice formation             [m of ice] 
    112    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_melt_1d  !: Net melting [m], for mono-cat lateral melting 
    113112   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sm_i_1d       !: Ice bulk salinity [ppt] 
    114113   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new       !: Salinity of new ice at the bottom 
    … …  
    167166         &      ht_s_1d    (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &     
    168167         &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &     
    169          &      dh_snowice(jpij) , dh_i_melt_1d(jpij) , & 
     168         &      dh_snowice(jpij) , & 
    170169         &      sm_i_1d   (jpij) , s_i_new  (jpij) ,    & 
    171170         &      t_s_1d(jpij,nlay_s),                                       & 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90

    r4689 r5053  
    3535   USE par_ice 
    3636   USE ice 
    37    USE limcat_1D          ! redistribute ice input into categories 
     37   USE limvar          ! redistribute ice input into categories 
    3838#endif 
    3939   USE sbcapr 
    … …  
    380380#if defined key_lim3 
    381381               IF( .NOT. ll_bdylim3 .AND. cn_ice_lim(ib_bdy) /= 'none' .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is lim2 type) 
    382                 CALL lim_cat_1D ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
     382                CALL lim_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), & 
    383383                                  & dta_bdy(ib_bdy)%ht_i,     dta_bdy(ib_bdy)%ht_s,     dta_bdy(ib_bdy)%a_i     ) 
    384384               ENDIF 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90

    r4990 r5053  
    6060      !!---------------------------------------------------------------------- 
    6161      INTEGER, INTENT( in ) :: kt     ! Main time step counter 
    62       !! 
    6362      INTEGER               :: ib_bdy ! Loop index 
     63 
    6464      DO ib_bdy=1, nb_bdy 
    6565 
    … …  
    7272            CALL ctl_stop( 'bdy_ice_lim : unrecognised option for open boundaries for ice fields' ) 
    7373         END SELECT 
    74       ENDDO 
     74 
     75      END DO 
    7576 
    7677   END SUBROUTINE bdy_ice_lim 
    … …  
    247248            END DO 
    248249 
    249  
    250          END DO !jb 
     250         END DO 
    251251  
    252          CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )                                         ! lateral boundary conditions 
     252         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy ) 
    253253         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy ) 
    254254         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy ) 

  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r5051 r5053  
    3939   USE limtrp          ! Ice transport 
    4040   USE limthd          ! Ice thermodynamics 
    41    USE limitd_th       ! Thermodynamics on ice thickness distribution  
    4241   USE limitd_me       ! Mechanics on ice thickness distribution 
    4342   USE limsbc          ! sea surface boundary condition 
    … …  
    109108      INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE, =5 COUPLED) 
    110109      !! 
    111       INTEGER  ::   jl      ! dummy loop index 
    112       REAL(wp) ::   zcoef   ! local scalar 
     110      INTEGER  ::   jl                 ! dummy loop index 
    113111      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky 
    114112      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice          ! mean ice albedo (for coupled) 
    … …  
    117115      IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim') 
    118116 
    119       !                                        !----------------------! 
    120       IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  ! 
    121          !                                     !----------------------! 
    122          !                                           !  Bulk Formulae ! 
    123          !                                           !----------------! 
    124          ! 
     117      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only 
     118         !-----------------------!                                            
     119         ! --- Bulk Formulae --- !                                            
     120         !-----------------------! 
    125121         u_oce(:,:) = ssu_m(:,:) * umask(:,:,1)      ! mean surface ocean current at ice velocity point 
    126122         v_oce(:,:) = ssv_m(:,:) * vmask(:,:,1)      ! (C-grid dynamics :  U- & V-points as the ocean) 
    … …  
    133129         CALL albedo_ice( t_su, ht_i, ht_s, zalb_cs, zalb_os )  ! cloud-sky and overcast-sky ice albedos 
    134130 
     131         ! CORE and COUPLED bulk formulations 
    135132         SELECT CASE( kblk ) 
    136          CASE( jp_core , jp_cpl )   ! CORE and COUPLED bulk formulations 
     133         CASE( jp_core , jp_cpl ) 
    137134 
    138135            ! albedo depends on cloud fraction because of non-linear spectral effects 
    … …  
    183180         END SELECT 
    184181          
    185          !                                           !----------------------! 
    186          !                                           ! LIM-3  time-stepping ! 
    187          !                                           !----------------------! 
    188          !  
     182         !------------------------------! 
     183         ! --- LIM-3 main time-step --- ! 
     184         !------------------------------! 
    189185         numit = numit + nn_fsbc                     ! Ice model time step 
    190186         ! 
    … …  
    200196         v_ice_b(:,:)     = v_ice(:,:) 
    201197 
    202                           CALL sbc_lim_flx0          ! set diag of mass, heat and salt fluxes to 0 
     198                          CALL sbc_lim_diag0         ! set diag of mass, heat and salt fluxes to 0 
    203199 
    204200                          CALL lim_rst_opn( kt )     ! Open Ice restart file 
    205201         ! 
    206          IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print 
    207202         ! ---------------------------------------------- 
    208203         ! ice dynamics and transport (except in 1D case) 
    … …  
    211206 
    212207                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics ) 
     208 
    213209                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion ) 
    214                           CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
    215          IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
     210 
    216211         IF( nn_monocat /= 2 )   & 
    217212            &             CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting) 
    218                           CALL lim_var_agg( 1 )  
     213 
    219214#if defined key_bdy 
    220215                          ! bdy ice thermo  
    … …  
    225220         IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermo bdy - ' )   ! control print 
    226221#endif 
    227                           CALL lim_update1 
     222 
     223                          CALL lim_update1( kt ) 
     224 
    228225         ENDIF 
    229 !                         !- Change old values for new values 
    230                           u_ice_b(:,:)     = u_ice(:,:) 
    231                           v_ice_b(:,:)     = v_ice(:,:) 
    232                           a_i_b  (:,:,:)   = a_i  (:,:,:) 
    233                           v_s_b  (:,:,:)   = v_s  (:,:,:) 
    234                           v_i_b  (:,:,:)   = v_i  (:,:,:) 
    235                           e_s_b  (:,:,:,:) = e_s  (:,:,:,:) 
    236                           e_i_b  (:,:,:,:) = e_i  (:,:,:,:) 
    237                           oa_i_b (:,:,:)   = oa_i (:,:,:) 
    238                           smv_i_b(:,:,:)   = smv_i(:,:,:) 
     226 
     227         !- Change old values for new values 
     228         u_ice_b(:,:)     = u_ice(:,:) 
     229         v_ice_b(:,:)     = v_ice(:,:) 
     230         a_i_b  (:,:,:)   = a_i  (:,:,:) 
     231         v_s_b  (:,:,:)   = v_s  (:,:,:) 
     232         v_i_b  (:,:,:)   = v_i  (:,:,:) 
     233         e_s_b  (:,:,:,:) = e_s  (:,:,:,:) 
     234         e_i_b  (:,:,:,:) = e_i  (:,:,:,:) 
     235         oa_i_b (:,:,:)   = oa_i (:,:,:) 
     236         smv_i_b(:,:,:)   = smv_i(:,:,:) 
    239237  
    240238         ! ---------------------------------------------- 
    … …  
    253251                             IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
    254252                          &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
    255                            ! Latent heat flux is forced to 0 in coupled : 
    256                            !  it is included in qns (non-solar heat flux) 
     253                             ! Latent heat flux is forced to 0 in coupled: it is included in qns (non-solar heat flux) 
    257254                             qla_ice  (:,:,:) = 0._wp 
    258255                             dqla_ice (:,:,:) = 0._wp 
    259256                          END SELECT 
    260257                          ! 
    261                           CALL lim_var_bv                 ! bulk brine volume (diag) 
    262258                          CALL lim_thd( kt )              ! Ice thermodynamics  
    263                           zcoef = rdt_ice /rday           !  Ice natural aging 
    264                           oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef 
    265          IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 1, ' - ice thermodyn. - ' )   ! control print 
    266                           CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion 
    267                           CALL lim_update2                ! Global variables update 
    268  
    269          IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 2, ' - Final state - ' )   ! control print 
    270          ! 
    271                           CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes 
    272          ! 
    273          IF( ln_nicep )   CALL lim_prt( kt, jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print 
     259 
     260                          CALL lim_update2( kt )          ! Global variables update 
     261         ! 
     262                          CALL lim_sbc_flx( kt )          ! Update surface ocean mass, heat and salt fluxes 
    274263         ! 
    275264         IF(ln_limdiaout) CALL lim_diahsb                 ! Diagnostics and outputs  
    276265 
    277                           CALL lim_wri( 1  )              ! Ice outputs  
     266                          CALL lim_wri( 1 )               ! Ice outputs  
    278267 
    279268         IF( kt == nit000 .AND. ln_rstart )   & 
    280             &             CALL iom_close( numrir )        ! clem: close input ice restart file 
     269            &             CALL iom_close( numrir )        ! close input ice restart file 
    281270         ! 
    282271         IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file  
    … …  
    315304      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM-3) time stepping' 
    316305      ! 
    317       ! 
    318306      !                                ! Allocate the ice arrays 
    319307      ierr =        ice_alloc        ()      ! ice variables 
    … …  
    343331      CALL lim_thd_sal_init            ! set ice salinity parameters 
    344332      ! 
    345       rdt_ice   = nn_fsbc * rdttra(1)  ! sea-ice timestep 
    346       r1_rdtice = 1._wp / rdt_ice      ! sea-ice timestep inverse 
    347       ! 
    348333      CALL lim_msh                     ! ice mesh initialization 
    349334      ! 
    350       CALL lim_itd_init                 ! ice thickness distribution initialization 
     335      CALL lim_itd_init                ! ice thickness distribution initialization 
    351336      ! 
    352337      CALL lim_itd_me_init             ! ice thickness distribution initialization 
    … …  
    356341         numit = nit000 - 1 
    357342         CALL lim_istate 
    358          CALL lim_var_agg(1) 
    359          CALL lim_var_glo2eqv 
    360343      ELSE                                    ! start from a restart file 
    361344         CALL lim_rst_read 
    362345         numit = nit000 - 1 
    363          CALL lim_var_agg(1) 
    364          CALL lim_var_glo2eqv 
    365346      ENDIF 
     347      CALL lim_var_agg(1) 
     348      CALL lim_var_glo2eqv 
    366349      ! 
    367350      CALL lim_sbc_init                 ! ice surface boundary condition    
    … …  
    390373      !! ** input   :   Namelist namicerun 
    391374      !!------------------------------------------------------------------- 
    392       NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, ln_nicep, ln_limdiahsb, ln_limdiaout 
     375!clem      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, cn_icerst_in, cn_icerst_out,   & 
     376!clem         &                ln_limdyn, amax, ln_nicep, ln_limdiahsb, ln_limdiaout 
     377      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out,   & 
     378         &                ln_limdyn, amax, ln_nicep, ln_limdiahsb, ln_limdiaout 
    393379      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    394380      !!------------------------------------------------------------------- 
    … …  
    408394         WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice' 
    409395         WRITE(numout,*) ' ~~~~~~' 
     396!clem         WRITE(numout,*) '   number of ice  categories                               = ', jpl 
     397!clem         WRITE(numout,*) '   number of ice  layers                                   = ', nlay_i 
     398!clem         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s 
    410399         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn 
    411400         WRITE(numout,*) '   maximum ice concentration                               = ', amax  
    … …  
    423412         IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    424413      ENDIF 
     414      ! 
     415      ! sea-ice timestep and inverse 
     416      rdt_ice   = nn_fsbc * rdttra(1)   
     417      r1_rdtice = 1._wp / rdt_ice  
    425418      ! 
    426419   END SUBROUTINE ice_run 
    … …  
    577570   END SUBROUTINE ice_lim_flx 
    578571 
    579    SUBROUTINE sbc_lim_flx0 
     572   SUBROUTINE sbc_lim_diag0 
    580573      !!---------------------------------------------------------------------- 
    581574      !!                  ***  ROUTINE sbc_lim_flx0  *** 
    … …  
    610603      afx_dyn(:,:) = 0._wp   ;   afx_thd(:,:) = 0._wp 
    611604       
    612    END SUBROUTINE sbc_lim_flx0 
     605   END SUBROUTINE sbc_lim_diag0 
    613606       
    614607   FUNCTION fice_cell_ave ( ptab ) 
    … …  
    623616       
    624617      DO jl = 1, jpl 
    625          fice_cell_ave (:,:) = fice_cell_ave (:,:) & 
    626             &                  + a_i (:,:,jl) * ptab (:,:,jl) 
     618         fice_cell_ave (:,:) = fice_cell_ave (:,:) + a_i (:,:,jl) * ptab (:,:,jl) 
    627619      END DO 
    628620       
Note: See TracChangeset for help on using the changeset viewer.