Changeset 5048


Ignore:
Timestamp:
2015-02-02T11:28:50+01:00 (7 years ago)
Author:
vancop
Message:

new itd boundaries, namelist changes, mono-category and comments

Location:
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM
Files:
13 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/CONFIG/SHARED/field_def.xml

    r4996 r5048  
    297297         <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   /> 
    298298         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   /> 
     299 
     300         <field id="afxtot"       long_name="area tendency (total)"                                        unit="s-1"     /> 
     301         <field id="afxdyn"       long_name="area tendency (dynamics)"                                     unit="s-1"     /> 
     302         <field id="afxthd"       long_name="area tendency (thermo)"                                       unit="s-1"     /> 
    299303 
    300304         <field id="hfxsum"   long_name="heat fluxes causing surface ice melt"            unit="W/m2"  /> 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r5046 r5048  
    22!! NEMO/LIM3 :  1 - dynamics/advection/thermo          (namicerun) 
    33!! namelists    2 - ice intialisation                  (namiceini) 
    4 !!              3 - ice dynamic                        (namicedyn) 
     4!!              3 - ice dynamics                       (namicedyn) 
    55!!              4 - ice advection                      (namicetrp) 
    6 !!              5 - thermodynamic                      (namicethd) 
     6!!              5 - thermodynamics                     (namicethd) 
    77!!              6 - ice salinity                       (namicesal) 
    88!!              7 - mechanical redistribution of ice   (namiceitdme) 
     
    1717   cn_icerst_out = "restart_ice"   !  suffix of ice restart name (output) 
    1818   ln_limdyn     = .true.          !  ice dynamics (T) or thermodynamics only (F) 
    19    amax          = 0.999           !  maximum ice concentration 
    20    ln_nicep      = .false.         !  Ice points output for debug (yes or no) 
    21    ln_limdiahsb  = .false.          !  check the heat and salt budgets (T) or not (F) 
     19   amax          = 0.999           !  maximum tolerated ice concentration 
     20   ln_nicep      = .false.         !  ice points output for debug (yes or no) 
     21   ln_limdiahsb  = .false.         !  check the heat and salt budgets (T) or not (F) 
    2222   ln_limdiaout  = .true.          !  output the heat and salt budgets (T) or not (F) 
    2323/ 
    2424!----------------------------------------------------------------------- 
    25 &namiceini     !   ice initialisation 
     25&namiceini     !   ice initialization 
    2626!----------------------------------------------------------------------- 
    2727   ln_limini   = .false.   !  activate ice initialization (T) or not (F) 
    2828   thres_sst   =  0.0      !  threshold water temperature for initial sea ice 
    29    hts_ini_n   =  0.3      !  initial snow thickness in the north 
    30    hts_ini_s   =  0.3      !        "            "          south 
    31    hti_ini_n   =  3.0      !  initial ice thickness in the north 
    32    hti_ini_s   =  1.0      !        "            "         south 
    33    ati_ini_n   =  0.9      !  initial ice concentration in the north 
    34    ati_ini_s   =  0.9      !        "            "             south 
    35    smi_ini_n   =  6.3      !  initial ice salinity in the north 
    36    smi_ini_s   =  6.3      !        "            "    south 
    37    tmi_ini_n   =  270.     !  initial ice/snw temp in the north 
    38    tmi_ini_s   =  270.     !  initial ice/snw temp in the south 
     29   hts_ini_n   =  0.3      !  initial snow thickness (North) 
     30   hts_ini_s   =  0.3      !        "            "   (South) 
     31   hti_ini_n   =  3.0      !  initial ice thickness  (North) 
     32   hti_ini_s   =  1.0      !        "            "   (South) 
     33   ati_ini_n   =  0.9      !  initial ice concentration (North) 
     34   ati_ini_s   =  0.9      !        "            "      (South) 
     35   smi_ini_n   =  6.3      !  initial ice salinity   (North) 
     36   smi_ini_s   =  6.3      !        "            "   (South) 
     37   tmi_ini_n   =  270.     !  initial ice/snw temp   (North) 
     38   tmi_ini_s   =  270.     !  initial ice/snw temp   (South) 
    3939/ 
    4040!----------------------------------------------------------------------- 
    41 &namicedyn     !   ice dynamic 
     41&namicedyn     !   ice dynamics 
    4242!----------------------------------------------------------------------- 
    43    cw          =   5.0e-03 !  drag coefficient for oceanic stress 
    44    pstar       =   2.0e+04 !  1st bulk-rheology parameter 
     43   cw          =   5.0e-03 !  ice-ocean drag coefficient 
     44   pstar       =   2.0e+04 !  1st bulk-rheology parameter (N/m) 
    4545   c_rhg       =  20.0     !  2nd bulk-rhelogy parameter 
    46    creepl      =   1.0e-12 !  creep limit 
     46   creepl      =   1.0e-12 !  creep limit (s-1) 
    4747   ecc         =   2.0     !  eccentricity of the elliptical yield curve 
    48    ahi0        = 350.e0    !  horizontal eddy diffusivity coefficient for sea-ice [m2/s] ! depends on resolution 
    49    nevp        = 120       !  number of iterations for subcycling in EVP 
    50    relast      = 0.333     !  ratio of elastic timescale over ice time step (1/3 if nevp=120 ; 1/9 if nevp=300) 
    51    hminrhg     =   0.001   !  ice volume (a*h in m) below which ice velocity equal ocean velocity 
     48   ahi0        = 350.e0    !  horizontal eddy diffusivity coefficient for sea-ice (m2/s), depends on resolution 
     49   nevp        = 120       !  number of EVP subcycles 
     50   relast      =   0.333   !  ratio of elastic timescale to ice time step (1/3 if nevp=120, 1/9 if nevp=300) 
     51   hminrhg     =   0.001   !  ice volume (a*h in m) below which ice velocity equals ocean velocity 
    5252/ 
    5353!----------------------------------------------------------------------- 
    54 &namicethd     !   ice thermodynamic 
     54&namicethd     !   ice thermodynamics 
    5555!----------------------------------------------------------------------- 
    56    hiccrit     = 0.1       !  ice thickness for lateral accretion  
    57                            !         caution 1.0, 1.0 best value to be used!!! (gilles G.)  ???? 
    58    fraz_swi    = 0         !  use of frazil ice collection thickness in function of wind (1.0) or not (0.0) 
     56   hiccrit     = 0.1       !  thickness for new ice formation in open water 
     57   fraz_swi    = 0         !  use frazil ice collection thickness as a function of wind (1.0) or not (0.0) 
    5958   maxfrazb    = 0.0       !  maximum portion of frazil ice collecting at the ice bottom 
    60    vfrazb      = 0.4166667 !  thresold drift speed for frazil ice collecting at the ice bottom 
     59   vfrazb      = 0.417    !  thresold drift speed for frazil ice collecting at the ice bottom 
    6160   Cfrazb      = 5.0       !  squeezing coefficient for frazil ice collecting at the ice bottom 
    6261   hiclim      = 0.10      !  minimum thickness of the first thickness category, must be smaller than hiccrit 
    6362   betas       = 0.6       !  exponent in lead-ice fractionation of snow precipitation 0.66 
    6463                           !        betas = 1 -> equipartition, betas < 1 -> more on leads 
    65    kappa_i     = 1.0       !  extinction radiation parameter in sea ice (1.0) 
     64   kappa_i     = 1.0       !  extinction radiation parameter in sea ice (m-1) 
    6665   nconv_i_thd = 50        !  maximal number of iterations for heat diffusion computation 
    67    maxer_i_thd = 0.0001    !  maximal error in temperature for heat diffusion computation 
    68    thcon_i_swi = 1         !  switch for computation of thermal conductivity in the ice 
    69                            !        (0) Untersteiner (1964), (1) Pringle et al. (2007) 
     66   maxer_i_thd = 0.0001    !  temperature error tolerance after heat diffusion  
     67   thcon_i_swi = 1         !  sea ice thermal conductivity, Untersteiner_64 (0) or Pringle_et_al_07 (1) 
     68   nn_monocat  = 0         !  virtual ITD mono-category parameterizations (1) or not (0); 1 ok for jpl=1 only 
     69                           !  2=replace ridging by piling, 3=activate G(he) only, 4=activate lateral melting only) --- those are temporary test options 
    7070/ 
    7171!----------------------------------------------------------------------- 
    7272&namicesal     !   ice salinity 
    7373!----------------------------------------------------------------------- 
    74    num_sal     =  2        !  salinity option: 1 -> S = bulk_sal 
    75                            !                   2 -> S = S(z,t) with a simple parameterization 
    76                            !                   3 -> S = S(z) profile of Scwharzacher [1959] 
    77                            !                   4 -> S = S(h) Cox and Weeks [1974] 
     74   num_sal     =  2        !  S=bulk_sal (1), S(z,t) - (2), S(z) Schwarzacher_59 (3)  
    7875   bulk_sal    =  4.0      !  if 1 is used, it represents the ice salinity 
    79    sal_G       =  5.00     !  restoring salinity for GD 
    80    time_G      =  1.728e+6 !  restoring time for GD 
    81    sal_F       =  2.00     !  restoring salinity for flushing 
    82    time_F      =  8.640e+5 !  restoring time for flushing 
    83    s_i_max     = 20.0      !  Maximum salinity  
    84    s_i_min     =  0.1      !  Minimum tolerated ice salinity 
     76   sal_G       =  5.00     !  restoring salinity, gravity drainage (g/kg) 
     77   time_G      =  1.728e+6 !  restoring time scale, gravity drainage  (s) 
     78   sal_F       =  2.00     !  restoring salinity, flushing (g/kg) 
     79   time_F      =  8.640e+5 !  restoring time, flushing (s) 
     80   s_i_max     = 20.0      !  maximum tolerated ice salinity (g/kg) 
     81   s_i_min     =  0.1      !  minimum tolerated ice salinity (g/kg) 
    8582/ 
    8683!----------------------------------------------------------------------- 
    8784&namiceitdme   !   parameters for mechanical redistribution of ice  
    8885!----------------------------------------------------------------------- 
    89    ridge_scheme_swi = 0      !  which ridging scheme using (1=Rothrock,else=Hibler79) 
    90    Cs               =   0.50 !  shearing energy contribution to ridging 
    91    Cf               =  17.0  !  ratio of ridging work to PE change in ridging 
    92    fsnowrdg         =   0.5  !  snow fraction that survives in ridging 
    93    fsnowrft         =   0.5  !  snow fraction that survives in rafting 
    94    Gstar            =   0.15 !  fractional area of thin ice being ridged 
    95    astar            =   0.05 !  equivalent of gstar (0.05 for TH75 and 0.03 for weaker ice) 
    96    Hstar            = 100.0  !  parameter determining the maximum thickness of ridged ice 
    97    raft_swi         =   1    !  rafting or not 
    98    hparmeter        =   0.75 !  threshold thickness for rafting or not 
     86   ridge_scheme_swi =   0    !  ice strength parameteriztaion, Hibler_79 (0), Rothrock_75 (1) 
     87   Cs               =   0.50 !  ratio of shearing energy contributing to ridging 
     88   Cf               =  17.0  !  ratio of ridging work to potential energy change in ridging 
     89   fsnowrdg         =   0.5  !  snow volume fraction that survives in ridging 
     90   fsnowrft         =   0.5  !  snow volume fraction that survives in rafting 
     91   Gstar            =   0.15 !  fractional area of thin ice being ridged (partfun_swi = 0) 
     92   astar            =   0.05 !  measures fractional area of thin ice being ridged (partfun_swi = 1) 
     93   Hstar            = 100.0  !  determines the maximum thickness of ridged ice (m) 
     94   raft_swi         =   1    !  rafting activated (1) or not (0) 
     95   hparmeter        =   0.75 !  threshold thickness for rafting (m) 
    9996   Craft            =   5.0  !  coefficient used in the rafting function 
    100    ridge_por        =   0.3  !  initial porosity of the ridged ice (typically 0.30) 
    101    partfun_swi      =   1    !  participation function linear, TH75 (0) or exponential Letal07 (1) 
    102    brinstren_swi    =   0    !  (1) use brine volume to diminish ice strength 
     97   ridge_por        =   0.3  !  porosity of newly ridged ice 
     98   partfun_swi      =   1    !  participation function: linear, Thorndike et al 75 (0) or exponential, Lipscomb 07 (1) 
     99   brinstren_swi    =   0    !  ice strength is a function of brine volume fraction (1) or not (0) 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r5047 r5048  
    165165   REAL(wp), PUBLIC ::   r1_rdtice        !: = 1. / rdt_ice 
    166166 
    167    !                                     !!** ice-dynamic namelist (namicedyn) ** 
     167   ! Clem je sais pas ou les mettre 
     168   INTEGER , PUBLIC ::   nn_hibnd         !: thickness category boundaries: tanh function (1), or h^(-alpha) function (2) 
     169   REAL(wp), PUBLIC ::   rn_hibnd         !: mean thickness of the domain (used to compute category boundaries, nn_hibnd = 2 only) 
     170 
     171   !                                     !!** ice-dynamics namelist (namicedyn) ** 
    168172   INTEGER , PUBLIC ::   nevp             !: number of iterations for subcycling 
    169173   REAL(wp), PUBLIC ::   cw               !: drag coefficient for oceanic stress 
     
    191195   !                                         ! 3 - salinity profile, constant in time 
    192196   INTEGER , PUBLIC ::   thcon_i_swi         !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007) 
     197   INTEGER , PUBLIC ::   nn_monocat          !: virtual ITD mono-category parameterizations (1) or not (0) 
    193198 
    194199   !                                     !!** ice-mechanical redistribution namelist (namiceitdme) 
     
    260265   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res    !: residual component of wfx_ice [kg/m2] 
    261266 
     267   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total) [s-1] 
     268   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_thd     !: ice concentration tendency (thermodynamics) [s-1] 
     269   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics) [s-1] 
     270 
    262271   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s] 
    263272   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bom     !: salt flux due to ice growth/melt                      [PSU/m2/s] 
     
    290299 
    291300   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice 
     301 
     302   ! lateral melting (mono-category) 
     303   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_melt !: net melting (m) 
    292304 
    293305   ! temporary arrays for dummy version of the code 
     
    436448         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     & 
    437449         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,  qlead  (jpi,jpj) ,     & 
    438          &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl) ,      & 
     450         &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) , & 
     451         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl) , dh_i_melt(jpi,jpj,jpl) , & 
    439452         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,      & 
    440453         &      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/iceini.F90

    r5047 r5048  
    172172      !!              limistate (only) and is changed to 99 m in ice_init 
    173173      !!------------------------------------------------------------------ 
    174       INTEGER  ::   jl                   ! dummy loop index 
    175       REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars 
     174      INTEGER  ::   jl                        ! dummy loop index 
     175      REAL(wp) ::   zc1, zc2, zc3, zx1        ! local scalars 
     176      REAL(wp) ::   zhmax, znum, zden, zalpha ! 
    176177      !!------------------------------------------------------------------ 
    177178 
     
    189190      !- Thickness categories boundaries  
    190191      !---------------------------------- 
     192 
     193      !  Clem - je sais pas encore dans quelle namelist les mettre, ca depend des chgts liés à bdy 
     194      nn_hibnd  = 2                !  thickness category boundaries: tanh function (1) h^(-alpha) (2) 
     195      rn_himean = 2.5              !  mean thickness of the domain (used to compute category boundaries, nn_hibnd = 2 only) 
     196 
    191197      hi_max(:) = 0._wp 
    192198 
    193       zc1 =  3._wp / REAL( jpl, wp ) 
    194       zc2 = 10._wp * zc1 
    195       zc3 =  3._wp 
    196  
    197       DO jl = 1, jpl 
    198          zx1 = REAL( jl-1, wp ) / REAL( jpl, wp ) 
    199          hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) ) 
    200       END DO 
     199      SELECT CASE ( nn_hbnd  )        
     200                                   !---------------------- 
     201         CASE (1)                  ! tanh function (CICE) 
     202                                   !---------------------- 
     203 
     204         zc1 =  3._wp / REAL( jpl, wp ) 
     205         zc2 = 10._wp * zc1 
     206         zc3 =  3._wp 
     207 
     208         DO jl = 1, jpl 
     209            zx1 = REAL( jl-1, wp ) / REAL( jpl, wp ) 
     210            hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) ) 
     211         END DO 
     212 
     213                                   !---------------------- 
     214         CASE (2)                  ! h^(-alpha) function 
     215                                   !---------------------- 
     216 
     217         zalpha = 0.05             ! exponent of the transform function 
     218 
     219         zhmax  = 3.*rn_himean 
     220 
     221         DO jl = 1, jpl  
     222            znum = jpl * ( zhmax+1 )**zalpha 
     223            zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl 
     224            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 
     225         END DO 
     226 
     227      END SELECT 
     228 
     229      ! mv- would be nice if we could put the instruction hi_max(jpl) = 99 here 
    201230 
    202231      IF(lwp) WRITE(numout,*) ' Thickness category boundaries ' 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r4990 r5048  
    9393      CALL lim_thd_lac 
    9494      CALL lim_var_glo2eqv    ! only for info 
     95 
     96      ! MV: Could put lateral melting here, would be better I think ??? 
     97 
    9598      
    9699      IF(ln_ctl) THEN   ! Control print 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r5047 r5048  
    8383      !! 
    8484      INTEGER  :: ji, jj, jk, jl   ! dummy loop indices 
    85       INTEGER  :: nbpb             ! nb of icy pts for thermo. cal. 
     85      INTEGER  :: nbpb             ! nb of icy pts for vertical thermo calculations 
     86      INTEGER  :: nbplm            ! nb of icy pts for lateral melting calculations (mono-cat) 
    8687      INTEGER  :: ii, ij           ! temporary dummy loop index 
    8788      REAL(wp) :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04) 
     
    434435              CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj ) 
    435436              CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb)   , jpi, jpj ) 
     437 
     438            IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN 
     439              CALL tab_1d_2d( nbpb, dh_i_melt(:,:,jl) , npb, dh_i_melt_1d(1:nbpb) , jpi, jpj ) 
     440            ENDIF 
    436441            ! 
    437442              CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj) 
     
    471476      !---------------------------------- 
    472477      CALL lim_var_eqv2glo 
     478 
     479      !---------------------------------- 
     480      ! 5.X) Lateral melting 
     481      !---------------------------------- 
     482      !!! declare dh_i_melt (ok), dh_i_melt_1d (ok), nbplm (ok), nplm (ok), zda(ok) 
     483 
     484      IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN 
     485 
     486         WRITE(numout,*) ' Lateral melting ON ' 
     487 
     488         ! select points where lateral melting occurs 
     489         jl = 1 
     490 
     491         nbplm = 0 
     492         DO jj = 1, jpj 
     493            DO ji = 1, jpi 
     494               IF ( ( dh_i_melt(ji,jj,jl)                  .LT.-epsi10 ) .AND.    & 
     495     &              ( ht_i(ji,jj,jl) - dh_i_melt(ji,jj,jl) .GT. epsi10 ) .AND.    &  
     496     &              ( ht_i(ji,jj,jl)                       .GT. epsi10 ) ) THEN      
     497                  nbplm   = nbplm  + 1 
     498                  nplm(nbplm) = (jj - 1) * jpi + ji 
     499               ENDIF 
     500            END DO 
     501         END DO 
     502 
     503         IF( nbplm > 0 ) THEN  ! If there is no net melting, do nothing 
     504 
     505            ! Move to 1D arrays 
     506            !------------------------- 
     507 
     508            CALL tab_2d_1d( nbplm, a_i_1d      (1:nbplm), a_i(:,:,jl)       , jpi, jpj, nplm(1:nbplm) ) 
     509            CALL tab_2d_1d( nbplm, ht_i_1d     (1:nbplm), ht_i(:,:,jl)      , jpi, jpj, nplm(1:nbplm) ) 
     510            CALL tab_2d_1d( nbplm, dh_i_melt_1d(1:nbplm), dh_i_melt(:,:,jl) , jpi, jpj, nplm(1:nbplm) ) 
     511 
     512            ! Compute lateral melting (dA = A/2h dh ) 
     513            DO ji = 1, nbplm 
     514               zda        = a_i_1d(ji) * dh_i_melt_1d(ji) / ( 2._wp * ht_i_1d(ji) ) 
     515               a_i_1d(ji) = a_i_1d(ji) + zda 
     516            END DO 
     517 
     518            ! Move back to 2D arrays 
     519            !------------------------- 
     520            CALL tab_1d_2d( nbplm, a_i (:,:,jl)  , nplm, a_i_1d     (1:nbplm)   , jpi, jpj ) 
     521            at_i(:,:) = a_i(:,:,jl) 
     522 
     523         ENDIF 
     524 
     525      ENDIF 
    473526 
    474527      !-------------------------------------------- 
     
    563616      !!------------------------------------------------------------------- 
    564617      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    565       NAMELIST/namicethd/ hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb,   & 
     618      NAMELIST/namicethd/ hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb,    & 
    566619         &                hiclim, parsub, betas,                          &  
    567          &                kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi 
     620         &                kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi, & 
     621         &                nn_monocat 
    568622      !!------------------------------------------------------------------- 
    569623      ! 
     
    582636902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp ) 
    583637      IF(lwm) WRITE ( numoni, namicethd ) 
     638      ! 
     639      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN  
     640         nn_monocat = 0 
     641         WRITE(numout, *) ' nn_monocat must be 0 in multi-category case ' 
     642      ENDIF 
    584643 
    585644      IF( lk_cpl .AND. parsub /= 0.0 )   CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' ) 
     
    597656         WRITE(numout,*)'      switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub   
    598657         WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          betas        = ', betas 
    599          WRITE(numout,*)'      extinction radiation parameter in sea ice (1.0)         kappa_i      = ', kappa_i 
     658         WRITE(numout,*)'      extinction radiation parameter in sea ice               kappa_i      = ', kappa_i 
    600659         WRITE(numout,*)'      maximal n. of iter. for heat diffusion computation      nconv_i_thd  = ', nconv_i_thd 
    601660         WRITE(numout,*)'      maximal err. on T for heat diffusion computation        maxer_i_thd  = ', maxer_i_thd 
    602661         WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     thcon_i_swi  = ', thcon_i_swi 
    603662         WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i 
     663         WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat 
    604664      ENDIF 
    605665      ! 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r5047 r5048  
    688688         END DO 
    689689      END DO 
     690       
     691      ! 
     692      !------------------------------------------------------------- 
     693      ! Net melting (input for thin ice melting, mono-category case) 
     694      !------------------------------------------------------------- 
     695      ! 
     696      IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 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 
    690701 
    691702      CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r5047 r5048  
    100100      INTEGER ::   nconv       ! number of iterations in iterative procedure 
    101101      INTEGER ::   minnumeqmin, maxnumeqmax 
     102       
    102103      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation 
    103104      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation 
    104105      INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow 
     106       
    105107      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    106108      REAL(wp) ::   zg1       =  2._wp        ! 
     
    113115      REAL(wp) ::   zerritmax   ! current maximal error on temperature  
    114116      REAL(wp) ::   zhsu 
    115       REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! ice melting point 
    116       REAL(wp), POINTER, DIMENSION(:) ::   ztsub       ! old surface temperature (before the iterative procedure ) 
    117       REAL(wp), POINTER, DIMENSION(:) ::   ztsubit     ! surface temperature at previous iteration 
    118       REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness 
    119       REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
    120       REAL(wp), POINTER, DIMENSION(:) ::   zfsw        ! solar radiation absorbed at the surface 
    121       REAL(wp), POINTER, DIMENSION(:) ::   zf          ! surface flux function 
    122       REAL(wp), POINTER, DIMENSION(:) ::   dzf         ! derivative of the surface flux function 
    123       REAL(wp), POINTER, DIMENSION(:) ::   zerrit      ! current error on temperature 
    124       REAL(wp), POINTER, DIMENSION(:) ::   zdifcase    ! case of the equation resolution (1->4) 
    125       REAL(wp), POINTER, DIMENSION(:) ::   zftrice     ! solar radiation transmitted through the ice 
    126       REAL(wp), POINTER, DIMENSION(:) ::   zihic 
    127       REAL(wp), POINTER, DIMENSION(:,:) ::   ztcond_i    ! Ice thermal conductivity 
    128       REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_i    ! Radiation transmitted through the ice 
    129       REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_i    ! Radiation absorbed in the ice 
    130       REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_i    ! Kappa factor in the ice 
    131       REAL(wp), POINTER, DIMENSION(:,:) ::   ztib        ! Old temperature in the ice 
    132       REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_i      ! Eta factor in the ice 
    133       REAL(wp), POINTER, DIMENSION(:,:) ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
    134       REAL(wp), POINTER, DIMENSION(:,:) ::   zspeche_i   ! Ice specific heat 
    135       REAL(wp), POINTER, DIMENSION(:,:) ::   z_i         ! Vertical cotes of the layers in the ice 
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_s    ! Radiation transmited through the snow 
    137       REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_s    ! Radiation absorbed in the snow 
    138       REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_s    ! Kappa factor in the snow 
    139       REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s      ! Eta factor in the snow 
    140       REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
    141       REAL(wp), POINTER, DIMENSION(:,:) ::   ztsb        ! Temporary temperature in the snow 
    142       REAL(wp), POINTER, DIMENSION(:,:) ::   z_s         ! Vertical cotes of the layers in the snow 
    143       REAL(wp), POINTER, DIMENSION(:,:) ::   zswiterm    ! Independent term 
    144       REAL(wp), POINTER, DIMENSION(:,:) ::   zswitbis    ! temporary independent term 
    145       REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis 
    146       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid     ! tridiagonal system terms 
     117       
     118      REAL(wp), POINTER, DIMENSION(:)     ::   ztfs        ! ice melting point 
     119      REAL(wp), POINTER, DIMENSION(:)     ::   ztsub       ! old surface temperature (before the iterative procedure ) 
     120      REAL(wp), POINTER, DIMENSION(:)     ::   ztsubit     ! surface temperature at previous iteration 
     121      REAL(wp), POINTER, DIMENSION(:)     ::   zh_i        ! ice layer thickness 
     122      REAL(wp), POINTER, DIMENSION(:)     ::   zh_s        ! snow layer thickness 
     123      REAL(wp), POINTER, DIMENSION(:)     ::   zfsw        ! solar radiation absorbed at the surface 
     124      REAL(wp), POINTER, DIMENSION(:)     ::   zf          ! surface flux function 
     125      REAL(wp), POINTER, DIMENSION(:)     ::   dzf         ! derivative of the surface flux function 
     126      REAL(wp), POINTER, DIMENSION(:)     ::   zerrit      ! current error on temperature 
     127      REAL(wp), POINTER, DIMENSION(:)     ::   zdifcase    ! case of the equation resolution (1->4) 
     128      REAL(wp), POINTER, DIMENSION(:)     ::   zftrice     ! solar radiation transmitted through the ice 
     129      REAL(wp), POINTER, DIMENSION(:)     ::   zihic 
     130       
     131      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztcond_i    ! Ice thermal conductivity 
     132      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_i    ! Radiation transmitted through the ice 
     133      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_i    ! Radiation absorbed in the ice 
     134      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_i    ! Kappa factor in the ice 
     135      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztib        ! Old temperature in the ice 
     136      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_i      ! Eta factor in the ice 
     137      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
     138      REAL(wp), POINTER, DIMENSION(:,:)   ::   zspeche_i   ! Ice specific heat 
     139      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_i         ! Vertical cotes of the layers in the ice 
     140      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_s    ! Radiation transmited through the snow 
     141      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_s    ! Radiation absorbed in the snow 
     142      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_s    ! Kappa factor in the snow 
     143      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_s      ! Eta factor in the snow 
     144      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence 
     145      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztsb        ! Temporary temperature in the snow 
     146      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_s         ! Vertical cotes of the layers in the snow 
     147      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindterm    ! 'Ind'ependent term 
     148      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindtbis    ! Temporary 'ind'ependent term 
     149      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdiagbis    ! Temporary 'dia'gonal term 
     150      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid       ! Tridiagonal system terms 
     151       
    147152      ! diag errors on heat 
    148       REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 
     153      REAL(wp), POINTER, DIMENSION(:)     :: zdq, zq_ini, zhfx_err 
     154       
     155      ! Mono-category 
     156      REAL(wp)                            :: zepsilon      ! determines thres. above which computation of G(h) is done 
     157      REAL(wp)                            :: zratio_s      ! dummy factor 
     158      REAL(wp)                            :: zratio_i      ! dummy factor 
     159      REAL(wp)                            :: zh_thres      ! thickness thres. for G(h) computation 
     160      REAL(wp)                            :: zhe           ! dummy factor 
     161      REAL(wp)                            :: zswitch       ! dummy switch 
     162      REAL(wp)                            :: zkimean       ! mean sea ice thermal conductivity 
     163      REAL(wp)                            :: zfac          ! dummy factor 
     164      REAL(wp)                            :: zihe          ! dummy factor 
     165      REAL(wp)                            :: zheshth       ! dummy factor 
     166       
     167      REAL(wp), POINTER, DIMENSION(:)     :: zghe          ! G(he), th. conduct enhancement factor, mono-cat 
     168       
    149169      !!------------------------------------------------------------------      
    150170      !  
    151171      CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 
    152172      CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    153       CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic ) 
     173      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 
    154174      CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 
    155175      CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 
    156       CALL wrk_alloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis  ) 
     176      CALL wrk_alloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis  ) 
    157177      CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 
    158178 
     
    200220      ! 
    201221      !------------------------------------------------------------------------------| 
    202       ! 2) Radiations                                                                | 
     222      ! 2) Radiation                                              | 
    203223      !------------------------------------------------------------------------------| 
    204224      ! 
     
    213233      ! zftrice = io.qsr_ice       is below the surface  
    214234      ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
     235      ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 
    215236      zhsu = 0.1_wp ! threshold for the computation of i0 
    216          !fr1_i0_1d = i0 for a thin ice surface 
    217          !fr1_i0_2d = i0 for a thick ice surface 
    218237      DO ji = kideb , kiut 
    219238         ! switches 
     
    339358            END DO 
    340359         ENDIF 
    341          ! 
    342          !------------------------------------------------------------------------------| 
    343          !  5) kappa factors                                                            | 
    344          !------------------------------------------------------------------------------| 
    345          ! 
     360          
     361         ! 
     362         !------------------------------------------------------------------------------| 
     363         !  6) G(he) - enhancement of thermal conductivity in mono-category case        | 
     364         !------------------------------------------------------------------------------| 
     365         ! 
     366         ! Computation of effective thermal conductivity G(h) 
     367         ! Used in mono-category case only to simulate an ITD implicitly 
     368         ! Fichefet and Morales Maqueda, JGR 1997 
     369 
     370         zghe(:) = 0._wp 
     371 
     372         SELECT CASE ( nn_monocat ) 
     373 
     374         CASE (0,2,4) 
     375 
     376            zghe(kideb:kiut) = 1._wp 
     377 
     378         CASE (1,3) ! LIM3 
     379 
     380            zepsilon = 0.1 
     381            zh_thres = EXP( 1._wp ) * zepsilon / 2. 
     382 
     383            DO ji = kideb, kiut 
     384    
     385               ! Mean sea ice thermal conductivity 
     386               zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL(nlay_i+1,wp) 
     387 
     388               ! Effective thickness he (zhe) 
     389               zfac     = 1._wp / ( rcdsn + zkimean ) 
     390               zratio_s = rcdsn   * zfac 
     391               zratio_i = zkimean * zfac 
     392               zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 
     393 
     394               ! G(he) 
     395               zswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if > 
     396               zghe(ji) = ( 1.0 - zswitch ) + zswitch * ( 0.5 + 0.5 * LOG( 2.*zhe / zepsilon ) ) 
     397    
     398               ! Impose G(he) < 2. 
     399               zghe(ji) = MIN( zghe(ji), 2.0 ) 
     400 
     401            END DO 
     402 
     403         END SELECT 
     404 
     405         ! 
     406         !------------------------------------------------------------------------------| 
     407         !  7) kappa factors                                                            | 
     408         !------------------------------------------------------------------------------| 
     409         ! 
     410         !--- Snow 
    346411         DO ji = kideb, kiut 
    347  
    348             !-- Snow kappa factors 
    349             zkappa_s(ji,0)         = rcdsn / MAX(epsi10,zh_s(ji)) 
    350             zkappa_s(ji,nlay_s)    = rcdsn / MAX(epsi10,zh_s(ji)) 
     412            zfac                  =  1. / MAX( epsi10 , zh_s(ji) ) 
     413            zkappa_s(ji,0)        = zghe(ji) * rcdsn * zfac 
     414            zkappa_s(ji,nlay_s)   = zghe(ji) * rcdsn * zfac 
    351415         END DO 
    352416 
    353417         DO jk = 1, nlay_s-1 
    354418            DO ji = kideb , kiut 
    355                zkappa_s(ji,jk)  = 2.0 * rcdsn / & 
    356                   MAX(epsi10,2.0*zh_s(ji)) 
    357             END DO 
    358          END DO 
    359  
     419               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0*zh_s(ji) ) 
     420            END DO 
     421         END DO 
     422 
     423         !--- Ice 
    360424         DO jk = 1, nlay_i-1 
    361425            DO ji = kideb , kiut 
    362                !-- Ice kappa factors 
    363                zkappa_i(ji,jk)  = 2.0*ztcond_i(ji,jk)/ & 
    364                   MAX(epsi10,2.0*zh_i(ji))  
    365             END DO 
    366          END DO 
    367  
    368          DO ji = kideb , kiut 
    369             zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji)) 
    370             zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji)) 
    371             !-- Interface 
    372             zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, & 
    373                (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 
    374             zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 
    375                + zkappa_i(ji,0)*REAL( 1 - isnow(ji) ) 
    376          END DO 
    377          ! 
    378          !------------------------------------------------------------------------------| 
    379          ! 6) Sea ice specific heat, eta factors                                        | 
     426               zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0*zh_i(ji) ) 
     427            END DO 
     428         END DO 
     429 
     430         !--- Snow-ice interface 
     431         DO ji = kideb , kiut 
     432            zfac                  = 1./ MAX( epsi10 , zh_i(ji) ) 
     433            zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac 
     434            zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 
     435            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / &  
     436           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) ) 
     437            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji), wp ) + zkappa_i(ji,0)*REAL( 1 - isnow(ji), wp ) 
     438         END DO 
     439 
     440         ! 
     441         !------------------------------------------------------------------------------| 
     442         ! 8) Sea ice specific heat, eta factors                                        | 
    380443         !------------------------------------------------------------------------------| 
    381444         ! 
     
    383446            DO ji = kideb , kiut 
    384447               ztitemp(ji,jk)   = t_i_1d(ji,jk) 
    385                zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ & 
    386                   MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10) 
    387                zeta_i(ji,jk)    = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 
    388                   epsi10) 
     448               zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ MAX( (t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt) , epsi10 ) 
     449               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic*zspeche_i(ji,jk)*zh_i(ji), epsi10 ) 
    389450            END DO 
    390451         END DO 
     
    396457            END DO 
    397458         END DO 
    398          ! 
    399          !------------------------------------------------------------------------------| 
    400          ! 7) surface flux computation                                                  | 
     459 
     460         ! 
     461         !------------------------------------------------------------------------------| 
     462         ! 9) surface flux computation                                                  | 
    401463         !------------------------------------------------------------------------------| 
    402464         ! 
     
    418480         ! 
    419481         !------------------------------------------------------------------------------| 
    420          ! 8) tridiagonal system terms                                                  | 
     482         ! 10) tridiagonal system terms                                                 | 
    421483         !------------------------------------------------------------------------------| 
    422484         ! 
     
    433495               ztrid(ji,numeq,2) = 0. 
    434496               ztrid(ji,numeq,3) = 0. 
    435                zswiterm(ji,numeq)= 0. 
    436                zswitbis(ji,numeq)= 0. 
     497               zindterm(ji,numeq)= 0. 
     498               zindtbis(ji,numeq)= 0. 
    437499               zdiagbis(ji,numeq)= 0. 
    438500            ENDDO 
     
    446508                  zkappa_i(ji,jk)) 
    447509               ztrid(ji,numeq,3)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk) 
    448                zswiterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* & 
     510               zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* & 
    449511                  zradab_i(ji,jk) 
    450512            END DO 
     
    458520               +  zkappa_i(ji,nlay_i-1) ) 
    459521            ztrid(ji,numeq,3)  =  0.0 
    460             zswiterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 
     522            zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 
    461523               ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 
    462524               *  t_bo_1d(ji) )  
     
    478540                     zkappa_s(ji,jk) ) 
    479541                  ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk) 
    480                   zswiterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* & 
     542                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* & 
    481543                     zradab_s(ji,jk) 
    482544               END DO 
     
    485547               IF ( nlay_i.eq.1 ) THEN 
    486548                  ztrid(ji,nlay_s+2,3)    =  0.0 
    487                   zswiterm(ji,nlay_s+2)   =  zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 
     549                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 
    488550                     t_bo_1d(ji)  
    489551               ENDIF 
     
    502564                  ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 
    503565                  ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 
    504                   zswiterm(ji,1) = dzf(ji)*t_su_1d(ji)   - zf(ji) 
     566                  zindterm(ji,1) = dzf(ji)*t_su_1d(ji)   - zf(ji) 
    505567 
    506568                  !!first layer of snow equation 
     
    508570                  ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 
    509571                  ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1) 
    510                   zswiterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 
     572                  zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 
    511573 
    512574               ELSE  
     
    525587                     zkappa_s(ji,0) * zg1s ) 
    526588                  ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1)  
    527                   zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            & 
     589                  zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            & 
    528590                     ( zradab_s(ji,1) +                         & 
    529591                     zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
     
    549611                  ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1     
    550612                  ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1 
    551                   zswiterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji) 
     613                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji) 
    552614 
    553615                  !!first layer of ice equation 
     
    556618                     + zkappa_i(ji,0) * zg1 ) 
    557619                  ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1)*zkappa_i(ji,1)   
    558                   zswiterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)   
     620                  zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)   
    559621 
    560622                  !!case of only one layer in the ice (surface & ice equations are altered) 
     
    569631                     ztrid(ji,numeqmin(ji)+1,3)  =  0.0 
    570632 
    571                      zswiterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* & 
     633                     zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* & 
    572634                        ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 
    573635                  ENDIF 
     
    589651                     zg1)   
    590652                  ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1) 
    591                   zswiterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 
     653                  zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 
    592654                     zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
    593655 
     
    598660                        zkappa_i(ji,1)) 
    599661                     ztrid(ji,numeqmin(ji),3)  =  0.0 
    600                      zswiterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* & 
     662                     zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* & 
    601663                        (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 
    602664                        + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 
     
    610672         ! 
    611673         !------------------------------------------------------------------------------| 
    612          ! 9) tridiagonal system solving                                                | 
     674         ! 11) tridiagonal system solving                                               | 
    613675         !------------------------------------------------------------------------------| 
    614676         ! 
     
    622684 
    623685         DO ji = kideb , kiut 
    624             zswitbis(ji,numeqmin(ji)) =  zswiterm(ji,numeqmin(ji)) 
     686            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji)) 
    625687            zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2) 
    626688            minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin) 
     
    633695               zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 
    634696                  ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 
    635                zswitbis(ji,numeq)  =  zswiterm(ji,numeq) - ztrid(ji,numeq,1)* & 
    636                   zswitbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 
     697               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1)* & 
     698                  zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 
    637699            END DO 
    638700         END DO 
     
    640702         DO ji = kideb , kiut 
    641703            ! ice temperatures 
    642             t_i_1d(ji,nlay_i)    =  zswitbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 
     704            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 
    643705         END DO 
    644706 
     
    646708            DO ji = kideb , kiut 
    647709               jk    =  numeq - nlay_s - 1 
    648                t_i_1d(ji,jk)  =  (zswitbis(ji,numeq) - ztrid(ji,numeq,3)* & 
     710               t_i_1d(ji,jk)  =  (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 
    649711                  t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 
    650712            END DO 
     
    654716            ! snow temperatures       
    655717            IF (ht_s_1d(ji).GT.0._wp) & 
    656                t_s_1d(ji,nlay_s)     =  (zswitbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 
     718               t_s_1d(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 
    657719               *  t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 
    658720               *        MAX(0.0,SIGN(1.0,ht_s_1d(ji)))  
     
    662724            ztsubit(ji) = t_su_1d(ji) 
    663725            IF( t_su_1d(ji) < ztfs(ji) ) & 
    664                t_su_1d(ji) = ( zswitbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   & 
     726               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   & 
    665727               &          + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
    666728         END DO 
    667729         ! 
    668730         !-------------------------------------------------------------------------- 
    669          !  10) Has the scheme converged ?, end of the iterative procedure         | 
     731         !  12) Has the scheme converged ?, end of the iterative procedure         | 
    670732         !-------------------------------------------------------------------------- 
    671733         ! 
     
    709771      ! 
    710772      !-------------------------------------------------------------------------! 
    711       !   11) Fluxes at the interfaces                                          ! 
     773      !   13) Fluxes at the interfaces                                          ! 
    712774      !-------------------------------------------------------------------------! 
    713775      DO ji = kideb, kiut 
     
    771833      CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 
    772834      CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 
    773       CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic ) 
     835      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 
    774836      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i,   & 
    775837         &              ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 
    776838      CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 
    777       CALL wrk_dealloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis ) 
     839      CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis ) 
    778840      CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 
    779841      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r4990 r5048  
    420420         END DO 
    421421         ! ------------------------------------------------- 
     422          
     423         !-------------------------------------- 
     424         ! Impose a_i < amax in mono-category 
     425         !-------------------------------------- 
     426         ! 
     427         IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2 
     428            DO jj = 1, jpj 
     429               DO ji = 1, jpi 
     430                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), amax ) 
     431               END DO 
     432            END DO 
     433         ENDIF 
    422434 
    423435         ! --- diags --- 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r4990 r5048  
    200200 
    201201      ztmp = rday / rhoic 
    202       CALL iom_put( "vfxres"     , wfx_res * ztmp  )             ! daily prod./melting due to limupdate  
    203       CALL iom_put( "vfxopw"     , wfx_opw * ztmp  )             ! daily lateral thermodynamic ice production 
    204       CALL iom_put( "vfxsni"     , wfx_sni * ztmp  )             ! daily snowice ice production 
    205       CALL iom_put( "vfxbog"     , wfx_bog * ztmp  )             ! daily bottom thermodynamic ice production 
    206       CALL iom_put( "vfxdyn"     , wfx_dyn * ztmp  )             ! daily dynamic ice production (rid/raft) 
    207       CALL iom_put( "vfxsum"     , wfx_sum * ztmp  )             ! surface melt  
    208       CALL iom_put( "vfxbom"     , wfx_bom * ztmp  )             ! bottom melt  
    209       CALL iom_put( "vfxice"     , wfx_ice * ztmp  )             ! total ice growth/melt  
    210       CALL iom_put( "vfxsnw"     , wfx_snw * ztmp  )             ! total snw growth/melt  
    211       CALL iom_put( "vfxsub"     , wfx_sub * ztmp  )             ! sublimation (snow)  
    212       CALL iom_put( "vfxspr"     , wfx_spr * ztmp  )             ! precip (snow)  
     202      CALL iom_put( "vfxres"     , wfx_res * ztmp       )        ! daily prod./melting due to limupdate  
     203      CALL iom_put( "vfxopw"     , wfx_opw * ztmp       )        ! daily lateral thermodynamic ice production 
     204      CALL iom_put( "vfxsni"     , wfx_sni * ztmp       )        ! daily snowice ice production 
     205      CALL iom_put( "vfxbog"     , wfx_bog * ztmp       )        ! daily bottom thermodynamic ice production 
     206      CALL iom_put( "vfxdyn"     , wfx_dyn * ztmp       )        ! daily dynamic ice production (rid/raft) 
     207      CALL iom_put( "vfxsum"     , wfx_sum * ztmp       )        ! surface melt  
     208      CALL iom_put( "vfxbom"     , wfx_bom * ztmp       )        ! bottom melt  
     209      CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt  
     210      CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt  
     211      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow)  
     212      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow) 
     213       
     214      CALL iom_put( "afxtot"     , afx_tot              )        ! concentration tendency (total) 
     215      CALL iom_put( "afxdyn"     , afx_dyn              )        ! concentration tendency (dynamics) 
     216      CALL iom_put( "afxthd"     , afx_thd              )        ! concentration tendency (thermo) 
    213217 
    214218      CALL iom_put ('hfxthd', hfx_thd(:,:) )   !   
     
    216220      CALL iom_put ('hfxres', hfx_res(:,:) )   !   
    217221      CALL iom_put ('hfxout', hfx_out(:,:) )   !   
    218       CALL iom_put ('hfxin' , hfx_in(:,:) )    
     222      CALL iom_put ('hfxin' , hfx_in(:,:) )     
    219223      CALL iom_put ('hfxsnw', hfx_snw(:,:) )   !   
    220224      CALL iom_put ('hfxsub', hfx_sub(:,:) )   !   
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90

    r4873 r5048  
    1212 
    1313   !                                             !!! ice thermodynamics 
    14    INTEGER, PUBLIC, PARAMETER ::   nlay_i   = 5   !: number of ice layers 
     14   INTEGER, PUBLIC, PARAMETER ::   nlay_i   = 2   !: number of ice layers 
    1515   INTEGER, PUBLIC, PARAMETER ::   nlay_s   = 1   !: number of snow layers 
    1616 
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r5047 r5048  
    3535   !: are the variables corresponding to 2d vectors 
    3636 
    37    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npb    !: number of points where computations has to be done 
    38    INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npac   !: correspondance between points (lateral accretion) 
     37   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npb    !: address vector for 1d vertical thermo computations 
     38   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   nplm   !: address vector for mono-category lateral melting 
     39   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npac   !: address vector for new ice formation 
    3940 
    4041   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlead_1d      
     
    109110   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott     !: Ice bottom accretion/ablation  [m] 
    110111   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 
    111113   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sm_i_1d       !: Ice bulk salinity [ppt] 
    112114   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new       !: Salinity of new ice at the bottom 
     
    138140      !!---------------------------------------------------------------------! 
    139141 
    140       ALLOCATE( npb      (jpij) , npac     (jpij),                          & 
     142      ALLOCATE( npb      (jpij) , nplm        (jpij) , npac     (jpij),   & 
    141143         !                                                                  ! 
    142144         &      qlead_1d (jpij) , ftr_ice_1d  (jpij) ,     & 
     
    165167         &      ht_s_1d    (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &     
    166168         &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &     
    167          &      dh_snowice(jpij) , sm_i_1d   (jpij) , s_i_new  (jpij) ,    & 
     169         &      dh_snowice(jpij) , dh_i_melt_1d(jpij) , & 
     170         &      sm_i_1d   (jpij) , s_i_new  (jpij) ,    & 
    168171         &      t_s_1d(jpij,nlay_s),                                       & 
    169172         &      t_i_1d(jpij,nlay_i+1), s_i_1d(jpij,nlay_i+1)                ,     &             
  • branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r4990 r5048  
    191191            CALL sbc_cpl_ice_tau( utau_ice , vtau_ice ) 
    192192 
    193             ! MV -> seb  
    194 !           CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=t_su    ) 
    195  
    196 !           IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   & 
    197 !              &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx ) 
    198 !           ! Latent heat flux is forced to 0 in coupled : 
    199 !           !  it is included in qns (non-solar heat flux) 
    200 !           qla_ice  (:,:,:) = 0._wp 
    201 !           dqla_ice (:,:,:) = 0._wp 
    202             ! END MV -> seb 
    203             ! 
    204193         END SELECT 
    205194          
     
    236225         wfx_spr(:,:) = 0._wp   ;    
    237226 
     227         afx_tot(:,:) = at_i(:,:) ;  afx_dyn(:,:) = 0._wp 
     228         afx_thd(:,:) = 0._wp 
     229 
    238230         hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
    239231         hfx_thd(:,:) = 0._wp   ;    
     
    252244         ! ---------------------------------------------- 
    253245         IF( .NOT. lk_c1d ) THEN 
     246 
     247         IF ( ln_limdyn ) afx_dyn(:,:)     = at_i(:,:) 
     248 
    254249                          CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics ) 
    255250                          CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion ) 
    256251                          CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting 
    257252         IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print 
    258                          CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting) 
     253         IF( nn_monocat /= 2 ) CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting) 
    259254                          CALL lim_var_agg( 1 )  
    260255#if defined key_bdy 
     
    278273                          oa_i_b (:,:,:)   = oa_i (:,:,:) 
    279274                          smv_i_b(:,:,:)   = smv_i(:,:,:) 
     275 
     276         IF ( ln_limdyn ) afx_dyn(:,:)     = ( at_i(:,:) - afx_dyn(:,:) ) * r1_rdtice 
     277                          afx_thd(:,:)     = at_i(:,:) 
    280278  
    281279         ! ---------------------------------------------- 
    282          ! ice thermodynamic 
     280         ! ice thermodynamics 
    283281         ! ---------------------------------------------- 
    284282                          CALL lim_var_glo2eqv            ! equivalent variables 
     
    288286                          phicif(:,:)  = vt_i(:,:) 
    289287 
    290                           ! MV -> seb 
    291288                          SELECT CASE( kblk ) 
    292289                             CASE ( jp_cpl ) 
     
    299296                             dqla_ice (:,:,:) = 0._wp 
    300297                          END SELECT 
    301                           ! END MV -> seb 
    302298                          ! 
    303299                          CALL lim_var_bv                 ! bulk brine volume (diag) 
     
    320316         !                                           ! Diagnostics and outputs  
    321317         IF (ln_limdiaout) CALL lim_diahsb 
     318 
     319                          afx_thd(:,:)     = ( at_i(:,:) - afx_thd(:,:) ) * r1_rdtice 
     320                          afx_tot(:,:)     = ( at_i(:,:) - afx_tot(:,:) ) * r1_rdtice 
    322321 
    323322                          CALL lim_wri( 1  )              ! Ice outputs  
Note: See TracChangeset for help on using the changeset viewer.